So far in this course, we have spent a lot of time thinking about which variables should be included in a model and how those variables can be used to explain patterns in our data. In some situations, though, we may have too many variables describing very similar things. If several variables are strongly correlated with each other, we may be carrying around a lot of repeated information.
For example, suppose we are studying students and collect their scores on algebra, calculus, statistics, physics, chemistry, and several other quantitative courses. A student who performs well in one of these classes may also tend to perform well in many of the others. We could keep every individual variable, but perhaps much of the information could be summarized by a smaller number of underlying measurements such as “general quantitative performance.”
This is the basic idea behind Principal Component Analysis, usually abbreviated as PCA. PCA is an unsupervised learning method that takes a collection of quantitative variables and creates a smaller collection of new variables called principal components. These new variables are linear combinations of our original variables and are designed to capture as much of the variation in the original dataset as possible.
In the previous lecture on K-Means Clustering, our goal was to reduce the number of observations into a few meaningful groups. PCA tackles a different problem: we are trying to reduce the number of features needed to describe those observations.
Let us begin with a small example. Suppose we collected two exam scores from a group of students:
exam1 <- c(52, 58, 61, 65, 68, 72, 74, 79, 83, 88, 91, 95)
exam2 <- c(50, 56, 64, 61, 70, 74, 77, 81, 80, 90, 94, 96)
exam_data <- data.frame(exam1, exam2)
plot(exam_data$exam1, exam_data$exam2, pch=19,
xlab="Exam 1", ylab="Exam 2", main="Two Highly Related Variables")
We can also calculate the correlation between the two variables:
cor(exam_data$exam1, exam_data$exam2)[1] 0.9868283The two exams contain different information, but they are clearly related. Students with high Exam 1 scores also tend to have high Exam 2 scores. This means that keeping both variables may be somewhat redundant. Instead of keeping the two original measurements, imagine drawing a line through the data in the direction where the observations vary the most. We could then describe each student by how far they are located along that line. A student toward one end of the line would generally have lower scores, while a student toward the other end would generally have higher scores.
That new direction is essentially what the first principal component is trying to find. This becomes much more useful when we have a larger number of variables. Suppose we have 20 variables but many of them are strongly correlated. PCA may reveal that most of the information in those 20 variables can be represented reasonably well using only 3 or 4 new variables. Instead of working in 20 dimensions, we may be able to work in a much smaller space.
This process is called dimensionality reduction. There are several reasons we may want to reduce the dimensionality of a dataset:
Of course, reducing the number of variables comes with a cost. If we replace 20 variables with only 3 principal components, we will almost certainly lose some information. PCA is therefore trying to find a compromise: use fewer dimensions while retaining as much information as possible.
To understand PCA, let us return to our two-exam example. We will standardize the two variables and then calculate the principal components using the prcomp() function.
exam_scaled <- scale(exam_data)
exam_pca <- prcomp(exam_data, center=TRUE, scale.=TRUE)
exam_pcaStandard deviations (1, .., p=2):
[1] 1.4095490 0.1147682
Rotation (n x k) = (2 x 2):
PC1 PC2
exam1 0.7071068 -0.7071068
exam2 0.7071068 0.7071068The first principal component, which we will call PC1, points in the direction where the data vary the most. The second principal component, PC2, must be perpendicular to PC1 and captures the largest remaining amount of variation.
We can visualize the two directions:
plot(exam_scaled[,1], exam_scaled[,2], pch=19,
xlab="Standardized Exam 1", ylab="Standardized Exam 2",
main="Principal Component Directions")
abline(a=0, b=exam_pca$rotation[2,1] / exam_pca$rotation[1,1], lwd=2)
abline(a=0, b=exam_pca$rotation[2,2] / exam_pca$rotation[1,2],
lty=2, lwd=2)
PC1 follows the direction where the observations are spread out the most. PC2 is perpendicular to PC1 and describes the remaining variation that PC1 did not capture. In this example, PC1 will largely represent a student’s overall exam performance because both Exam 1 and Exam 2 increase together. PC2 will instead represent differences between the two exams. A student with a relatively high Exam 1 score but a relatively low Exam 2 score would move in one direction along PC2, while the opposite type of student would move in the other direction.
The important idea is that PCA does not simply select one of our existing variables. Instead, it creates entirely new variables by combining the original features.
For standardized variables \(z_1,z_2,\ldots,z_p\), the first principal component has the general form
\[PC_1=a_{11}z_1+a_{12}z_2+\cdots+a_{1p}z_p,\]
where the \(a\) values are weights describing how much each original variable contributes to the component. These weights are called loadings.
The second principal component will have a different set of weights:
\[PC_2=a_{21}z_1+a_{22}z_2+\cdots+a_{2p}z_p.\]
We can see the loadings for our exam example using:
exam_pca$rotation PC1 PC2
exam1 0.7071068 -0.7071068
exam2 0.7071068 0.7071068Suppose PC1 gives similar positive weights to Exam 1 and Exam 2. Then a student with high standardized values on both exams will have a high PC1 score. This is why we might interpret PC1 as something like “overall exam performance.”
We can calculate the actual principal component scores for every student:
exam_pca$x PC1 PC2
[1,] -2.2924575 -0.03706484
[2,] -1.6963350 -0.06074494
[3,] -1.1597562 0.16593265
[4,] -1.0962659 -0.18377879
[5,] -0.5119835 0.09060233
[6,] -0.1145685 0.07481559
[7,] 0.1318425 0.11462575
[8,] 0.5809077 0.04718880
[9,] 0.7398050 -0.20711559
[10,] 1.4750913 0.01166863
[11,] 1.8208561 0.04753211
[12,] 2.1228640 -0.06366169This distinction is important:
We can verify the basic calculation for the first student. Since we standardized the variables before performing PCA, we multiply the standardized Exam 1 and Exam 2 values by their PC1 loadings:
exam_scaled[1,] exam1 exam2
-1.594803 -1.647221 exam_pca$rotation[,1] exam1 exam2
0.7071068 0.7071068 sum(exam_scaled[1,] * exam_pca$rotation[,1])[1] -2.292457exam_pca$x[1,1] PC1
-2.292457 The two values should match apart from very small rounding differences. So, although PCA may initially sound complicated, each PC score is really just a weighted combination of the original variables.
PCA chooses the loadings so that PC1 captures the largest possible amount of variation in the data. PC2 then captures as much of the remaining variation as possible while being perpendicular to PC1. PC3 does the same thing after accounting for PC1 and PC2, and the process continues.
Because the principal components are perpendicular to one another, they are also uncorrelated:
cor(exam_pca$x) PC1 PC2
PC1 1.000000e+00 -5.580649e-16
PC2 -5.580649e-16 1.000000e+00This is one of the reasons PCA can be useful when the original variables are highly correlated. We replace a set of correlated features with a new set of uncorrelated components.
There is a more mathematical explanation involving eigenvalues and eigenvectors of a covariance or correlation matrix, but we do not need to derive those calculations in this course. The important thing for us is understanding what PCA is optimizing and how we can interpret the output.
Let us move to a dataset with more than two variables. We will once again use the iris dataset, which contains four quantitative measurements for 150 flowers.
head(iris) Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 5.1 3.5 1.4 0.2 setosa
2 4.9 3.0 1.4 0.2 setosa
3 4.7 3.2 1.3 0.2 setosa
4 4.6 3.1 1.5 0.2 setosa
5 5.0 3.6 1.4 0.2 setosa
6 5.4 3.9 1.7 0.4 setosairis_data <- iris[,1:4]Before fitting PCA, it is useful to look at the correlations among the variables:
cor(iris_data) Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length 1.0000000 -0.1175698 0.8717538 0.8179411
Sepal.Width -0.1175698 1.0000000 -0.4284401 -0.3661259
Petal.Length 0.8717538 -0.4284401 1.0000000 0.9628654
Petal.Width 0.8179411 -0.3661259 0.9628654 1.0000000Notice that several of the measurements are related to one another. In particular, Petal Length and Petal Width are very strongly correlated. This suggests that there may be some redundant information in the original four variables.
We can perform PCA using prcomp():
iris_pca <- prcomp(iris_data, center=TRUE, scale.=TRUE)
iris_pcaStandard deviations (1, .., p=4):
[1] 1.7083611 0.9560494 0.3830886 0.1439265
Rotation (n x k) = (4 x 4):
PC1 PC2 PC3 PC4
Sepal.Length 0.5210659 -0.37741762 0.7195664 0.2612863
Sepal.Width -0.2693474 -0.92329566 -0.2443818 -0.1235096
Petal.Length 0.5804131 -0.02449161 -0.1421264 -0.8014492
Petal.Width 0.5648565 -0.06694199 -0.6342727 0.5235971There are a few important objects we can pull from the PCA results:
iris_pca$rotation PC1 PC2 PC3 PC4
Sepal.Length 0.5210659 -0.37741762 0.7195664 0.2612863
Sepal.Width -0.2693474 -0.92329566 -0.2443818 -0.1235096
Petal.Length 0.5804131 -0.02449161 -0.1421264 -0.8014492
Petal.Width 0.5648565 -0.06694199 -0.6342727 0.5235971head(iris_pca$x) PC1 PC2 PC3 PC4
[1,] -2.257141 -0.4784238 0.12727962 0.024087508
[2,] -2.074013 0.6718827 0.23382552 0.102662845
[3,] -2.356335 0.3407664 -0.04405390 0.028282305
[4,] -2.291707 0.5953999 -0.09098530 -0.065735340
[5,] -2.381863 -0.6446757 -0.01568565 -0.035802870
[6,] -2.068701 -1.4842053 -0.02687825 0.006586116iris_pca$sdev[1] 1.7083611 0.9560494 0.3830886 0.1439265The rotation object contains the loadings, the x object contains the principal component scores for each flower, and sdev gives the standard deviation of each principal component.
Since the first two principal components are now just two quantitative variables, we can plot all 150 flowers in two dimensions:
plot(iris_pca$x[,1], iris_pca$x[,2], pch=19,
xlab="PC1", ylab="PC2",
main="Iris Data Using the First Two Principal Components")
At this point, PCA has not used the flower species at all. It is still an unsupervised learning method. However, because we happen to know the species in this dataset, we can color the points afterward to see whether the PCA representation reveals any natural structure.
plot(iris_pca$x[,1], iris_pca$x[,2], col=as.numeric(iris$Species), pch=19,
xlab="PC1", ylab="PC2", main="Iris PCA Colored by Species")
This is one of the major uses of PCA. It is impossible to directly make a scatter plot containing four axes, but PCA allows us to compress those four measurements into two new variables and visualize the observations in two dimensions.
We should be careful with our interpretation though. PC1 is not one single original feature and PC2 is not another. Each principal component contains information from all of the original variables.
We can inspect the loadings to determine which variables contribute the most:
round(iris_pca$rotation, 3) PC1 PC2 PC3 PC4
Sepal.Length 0.521 -0.377 0.720 0.261
Sepal.Width -0.269 -0.923 -0.244 -0.124
Petal.Length 0.580 -0.024 -0.142 -0.801
Petal.Width 0.565 -0.067 -0.634 0.524For each component, variables with loadings farther from 0 contribute more strongly to that component. Variables with loadings closer to 0 have relatively little influence.
The sign also matters when comparing variables within the same component. If two variables have loadings with the same sign, larger values of those variables push observations in the same direction along that component. If they have opposite signs, they push observations in opposite directions.
However, the overall sign of a principal component is arbitrary. R could multiply every loading and every score for PC1 by \(-1\) and it would still represent exactly the same principal component. Therefore, we should not become too attached to whether a component points “left” or “right” in a plot. What matters is the relative pattern among the loadings.
For instance, if R gives us
\[PC_1=-0.52(\text{Sepal Length})+0.27(\text{Sepal Width})-0.58(\text{Petal Length})-0.56(\text{Petal Width}),\]
we might find it easier to mentally reverse every sign and describe the same component as increasing with Sepal Length, Petal Length, and Petal Width while decreasing with Sepal Width. The statistical information has not changed.
So far we have created four principal components from the four variables in iris. If we keep all four components, we have not actually reduced the dimensionality of the data. We have simply rewritten the same information using a different coordinate system. The real benefit occurs when the first few principal components explain most of the variation in the original dataset. We can calculate the variance of each principal component by squaring the standard deviation:
pc_variance <- iris_pca$sdev^2
pc_variance[1] 2.91849782 0.91403047 0.14675688 0.02071484We can then divide each variance by the total variance to determine the proportion of variance explained by each component:
prop_variance <- pc_variance / sum(pc_variance)
prop_variance[1] 0.729624454 0.228507618 0.036689219 0.005178709The cumulative proportion tells us how much information we retain if we keep the first 1, 2, 3, and so on principal components:
cumsum(prop_variance)[1] 0.7296245 0.9581321 0.9948213 1.0000000R will calculate this information for us automatically using:
summary(iris_pca)Importance of components:
PC1 PC2 PC3 PC4
Standard deviation 1.7084 0.9560 0.38309 0.14393
Proportion of Variance 0.7296 0.2285 0.03669 0.00518
Cumulative Proportion 0.7296 0.9581 0.99482 1.00000The Proportion of Variance row tells us how much variation each individual component explains, while the Cumulative Proportion row tells us how much is explained by all components up to that point.
For example, if PC1 explains 73% of the variation and PC2 explains another 23%, then using only the first two principal components would retain approximately 96% of the total variation while reducing four original variables down to two.
That is a pretty good trade: we have cut the number of dimensions in half while retaining most of the information in the original data.
We can visualize this using a scree plot:
plot(prop_variance, type="b", pch=19,
xlab="Principal Component", ylab="Proportion of Variance Explained",
main="Scree Plot for Iris PCA")
We could also plot the cumulative proportion of variance:
plot(cumsum(prop_variance), type="b", pch=19, ylim=c(0,1),
xlab="Number of Principal Components",
ylab="Cumulative Proportion of Variance",
main="Cumulative Variance Explained")
There is no universal rule saying that we must retain exactly 80%, 90%, or 95% of the variance. Those cutoffs can be useful guidelines, but the appropriate number of components depends on the goal of the analysis.
Much like choosing \(K\) in K-Means, we are balancing simplicity against information. Keeping more principal components preserves more of the original data, but it also reduces the benefit of dimensionality reduction.
Much like KNN and K-Means, PCA can be heavily affected by the scale of the variables. This is because PCA looks for directions with the largest amount of variance. A variable measured on a very large numerical scale will usually have a larger variance and can dominate the first principal components.
Let us return to the USArrests dataset:
head(USArrests) Murder Assault UrbanPop Rape
Alabama 13.2 236 58 21.2
Alaska 10.0 263 48 44.5
Arizona 8.1 294 80 31.0
Arkansas 8.8 190 50 19.5
California 9.0 276 91 40.6
Colorado 7.9 204 78 38.7apply(USArrests, 2, mean) Murder Assault UrbanPop Rape
7.788 170.760 65.540 21.232 apply(USArrests, 2, sd) Murder Assault UrbanPop Rape
4.355510 83.337661 14.474763 9.366385 The variables are clearly measured on different scales. Assault has much larger values and a much larger standard deviation than the other variables. If we run PCA without scaling:
us_pca_unscaled <- prcomp(USArrests, center=TRUE, scale.=FALSE)
round(us_pca_unscaled$rotation, 3) PC1 PC2 PC3 PC4
Murder 0.042 -0.045 0.080 -0.995
Assault 0.995 -0.059 -0.068 0.039
UrbanPop 0.046 0.977 -0.201 -0.058
Rape 0.075 0.201 0.974 0.072summary(us_pca_unscaled)Importance of components:
PC1 PC2 PC3 PC4
Standard deviation 83.7324 14.21240 6.4894 2.48279
Proportion of Variance 0.9655 0.02782 0.0058 0.00085
Cumulative Proportion 0.9655 0.99335 0.9991 1.00000we will find that the high-variance variables have a very large influence on the principal components.
Instead, we can standardize the features:
us_pca <- prcomp(USArrests, center=TRUE, scale.=TRUE)
round(us_pca$rotation, 3) PC1 PC2 PC3 PC4
Murder -0.536 -0.418 0.341 0.649
Assault -0.583 -0.188 0.268 -0.743
UrbanPop -0.278 0.873 0.378 0.134
Rape -0.543 0.167 -0.818 0.089summary(us_pca)Importance of components:
PC1 PC2 PC3 PC4
Standard deviation 1.5749 0.9949 0.59713 0.41645
Proportion of Variance 0.6201 0.2474 0.08914 0.04336
Cumulative Proportion 0.6201 0.8675 0.95664 1.00000Setting scale.=TRUE essentially tells R to place the variables on comparable scales before finding the principal components. In practice, this means PCA is working with standardized variables rather than the original measurements.
As a general rule, if the variables use substantially different units or numerical scales, scaling is usually a good place to start. If the variables are already measured on the same meaningful scale and differences in variance are important to the problem, we may choose not to scale.
PCA can also be used as a preprocessing step before another machine learning method. Suppose we have a dataset containing many strongly correlated features. Rather than passing all of the original features into a model, we could first use PCA to create a smaller number of uncorrelated components and then fit the model using those components.
Let us connect this idea back to K-Means using mtcars. We will use six quantitative variables that describe the cars:
car_data <- mtcars[,c("mpg", "disp", "hp", "drat", "wt", "qsec")]
round(cor(car_data), 2) mpg disp hp drat wt qsec
mpg 1.00 -0.85 -0.78 0.68 -0.87 0.42
disp -0.85 1.00 0.79 -0.71 0.89 -0.43
hp -0.78 0.79 1.00 -0.45 0.66 -0.71
drat 0.68 -0.71 -0.45 1.00 -0.71 0.09
wt -0.87 0.89 0.66 -0.71 1.00 -0.17
qsec 0.42 -0.43 -0.71 0.09 -0.17 1.00Several of these variables are strongly correlated. For example, heavier cars often have larger displacement and lower fuel efficiency. Rather than clustering on all six original variables, one option is to first perform PCA:
car_pca <- prcomp(car_data, center=TRUE, scale.=TRUE)
summary(car_pca)Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6
Standard deviation 2.0463 1.0715 0.57737 0.39289 0.3533 0.22799
Proportion of Variance 0.6979 0.1913 0.05556 0.02573 0.0208 0.00866
Cumulative Proportion 0.6979 0.8892 0.94481 0.97054 0.9913 1.00000We can inspect the amount of variance explained:
car_prop <- car_pca$sdev^2 / sum(car_pca$sdev^2)
cumsum(car_prop)[1] 0.6978994 0.8892514 0.9448109 0.9705376 0.9913370 1.0000000Suppose we decide that the first three principal components retain enough of the variation for our purpose. We can create a reduced dataset:
car_reduced <- car_pca$x[,1:3]
head(car_reduced) PC1 PC2 PC3
Mazda RX4 -0.8425806 -0.873469391 0.2282783
Mazda RX4 Wag -0.8075041 -0.556341552 0.0126678
Datsun 710 -1.6850448 0.040006569 0.1564937
Hornet 4 Drive -0.0964443 1.294377904 0.5702297
Hornet Sportabout 1.2915096 0.006516693 0.5250741
Valiant 0.2187309 2.005957905 0.7258399Then we can run K-Means using those principal components:
set.seed(101)
car_kmeans_pca <- kmeans(car_reduced, centers=3, nstart=25)
car_kmeans_pca$size[1] 12 6 14Since the first two principal components are easy to visualize, we can color the cars by their K-Means cluster:
plot(car_pca$x[,1], car_pca$x[,2], col=car_kmeans_pca$cluster,
pch=19, xlab="PC1", ylab="PC2", main="K-Means Clusters in PCA Space")
text(car_pca$x[,1], car_pca$x[,2], labels=rownames(mtcars), pos=3, cex=.6)
This combination is fairly common. PCA reduces the dimensionality of the dataset, and then a clustering method tries to identify groups in the reduced space.
However, we should not automatically assume that PCA will always improve a later model. PCA is designed to preserve variation, not necessarily the specific information that another model needs. A lower-variance component could still contain something important for prediction or clustering.
There is also an important warning if we use PCA as part of a supervised learning problem. Suppose we split data into a training set and testing set before fitting a classification or regression model. We should fit PCA using only the training data. Otherwise, information from the testing set would influence the transformation and we would introduce data leakage.
Conceptually, the process would be:
# Fit PCA using the training predictors
pca_train <- prcomp(train_x, center=TRUE, scale.=TRUE)
# Transform the training data
train_pc <- pca_train$x[,1:k]
# Apply the SAME transformation to the testing data
test_pc <- predict(pca_train, newdata=test_x)[,1:k]The exact value of k would represent however many principal components we decided to retain.
This is the same general principle we have followed throughout machine learning: anything that learns information from the data should be learned using the training set and then applied to the testing set. PCA also has several limitations that we should keep in mind. One is that the principal components can be difficult to interpret, meaning the components combine all of the original variables rather than having one obvious meaning.