Lecture 3 PCA Live Demo: USArrests

Data and quick exploration

We perform PCA on the USArrests data set, which is part of the base R package. This data set contains statistics, in arrests per 100,000 residents for assault, murder, and rape in each of the 50 US states in 1973. Also given is the percent of the population living in urban areas.

The rows of the data set contain the 50 states, in alphabetical order.

data(USArrests)
X <- USArrests
states <- row.names(X)
states

 [1] "Alabama"        "Alaska"         "Arizona"        "Arkansas"       "California"     "Colorado"       "Connecticut"   
 [8] "Delaware"       "Florida"        "Georgia"        "Hawaii"         "Idaho"          "Illinois"       "Indiana"       
[15] "Iowa"           "Kansas"         "Kentucky"       "Louisiana"      "Maine"          "Maryland"       "Massachusetts" 
[22] "Michigan"       "Minnesota"      "Mississippi"    "Missouri"       "Montana"        "Nebraska"       "Nevada"        
[29] "New Hampshire"  "New Jersey"     "New Mexico"     "New York"       "North Carolina" "North Dakota"   "Ohio"          
[36] "Oklahoma"       "Oregon"         "Pennsylvania"   "Rhode Island"   "South Carolina" "South Dakota"   "Tennessee"     
[43] "Texas"          "Utah"           "Vermont"        "Virginia"       "Washington"     "West Virginia"  "Wisconsin"     
[50] "Wyoming"       

The columns of the data set contain the four variables.

summary(X)

     Murder         Assault       UrbanPop         Rape     
 Min.   : 0.80   Min.   : 45   Min.   :32.0   Min.   : 7.3  
 1st Qu.: 4.08   1st Qu.:109   1st Qu.:54.5   1st Qu.:15.1  
 Median : 7.25   Median :159   Median :66.0   Median :20.1  
 Mean   : 7.79   Mean   :171   Mean   :65.5   Mean   :21.2  
 3rd Qu.:11.25   3rd Qu.:249   3rd Qu.:77.8   3rd Qu.:26.2  
 Max.   :17.40   Max.   :337   Max.   :91.0   Max.   :46.0  

Visual exploration

We first briefly examine the data. Some variables are highly correlated.

pairs(X, pch = 19, col = "gray40", main = "USArrests: pairwise scatterplots")

The variables also have vastly different means and variances. We see that there are on average three times as many rapes as murders, and more than eight times as many assaults as rapes. The UrbanPop variable measures the percentage of the population in each state living in an urban area, which is not a comparable number to the number of rapes in each state per 100,000 individuals.

boxplot(X, col = c("#7FB3D5", "#F1948A", "#82E0AA", "#F7DC6F"),
        main = "Different variable scales")

PCA: no scaling vs scaling

We now perform principal components analysis using the prcomp() function, which is one of several functions in R that perform PCA.

By default, the prcomp() function centers the variables to have mean zero (we here add center = TRUE to explicitly confirm centering). By using the option scale = TRUE, we scale the variables to have standard deviation one. The output from prcomp() contains a number of useful quantities.

pca_cov <- prcomp(X, center = TRUE, scale. = FALSE)
pca_cor <- prcomp(X, center = TRUE, scale. = TRUE)

names(pca_cov)

[1] "sdev"     "rotation" "center"   "scale"    "x"       

The center and scale components correspond to the means and standard deviations of the variables that were used for scaling prior to implementing PCA.

pca_cor$center

  Murder  Assault UrbanPop     Rape 
   7.788  170.760   65.540   21.232 

pca_cor$scale

  Murder  Assault UrbanPop     Rape 
   4.356   83.338   14.475    9.366 

The rotation matrix provides the principal component loadings; each column of pca$rotation contains the corresponding principal component loading vector. Recall that the principal components are only unique up to a sign change.

pca_cor$rotation

             PC1     PC2     PC3      PC4
Murder   -0.5359 -0.4182  0.3412  0.64923
Assault  -0.5832 -0.1880  0.2681 -0.74341
UrbanPop -0.2782  0.8728  0.3780  0.13388
Rape     -0.5434  0.1673 -0.8178  0.08902

If we failed to scale the variables before performing PCA, then most of the principal components that we observed would be driven by the Assault variable, since it has by far the largest mean and variance. Thus, it is important to standardize the variables to have mean zero and standard deviation one before performing PCA.

library(ggplot2)
library(reshape2)
library(patchwork)

# choose number of PCs

# prepare function
make_df <- function(pca, label) {
  loadings <- pca$rotation
  df <- melt(loadings)
  colnames(df) <- c("Variable", "PC", "Loading")
  df$Method <- label
  df
}

df_cov <- make_df(pca_cov, "Covariance PCA")
df_cor <- make_df(pca_cor, "Correlation PCA")

df_all <- rbind(df_cov, df_cor)

# keep ordering
df_all$Variable <- factor(df_all$Variable,
                         levels = rownames(pca_cov$rotation))

# plot
p <- ggplot(df_all, aes(x = PC, y = Variable, fill = Loading)) +
  geom_tile() +
  scale_fill_gradient2(midpoint = 0) +
  facet_wrap(~Method, nrow = 1) +
  theme_minimal() +
  labs(title = "PCA Loadings: Covariance vs Correlation",
       x = "Principal Component",
       y = "Variable")

p

Scree plots (side-by-side)

The prcomp() function also outputs the standard deviation of each principal component. On the USArrests data set, we can access these standard deviations as follows:

pca_cor$sdev

[1] 1.5749 0.9949 0.5971 0.4164

The variance explained by each principal component is obtained by squaring these. The scree plot shows these variances. We see below that without scaling, most of the variability is explained by the first principal component, which is essentially the Aasult variable.

par(mfrow = c(1, 2))
screeplot(pca_cov, type = "l", main = "Scree: covariance PCA")
screeplot(pca_cor, type = "l", main = "Scree: correlation PCA")

Explained variance and cumulative PVE

pve_cov <- (pca_cov$sdev^2) / sum(pca_cov$sdev^2)
pve_cor <- (pca_cor$sdev^2) / sum(pca_cor$sdev^2)

pve_cor

[1] 0.62006 0.24744 0.08914 0.04336

par(mfrow = c(1, 2))

plot(cumsum(pve_cov), type = "b", pch = 19,
     xlab = "Principal Component",
     ylab = "Cumulative PVE",
     main = "Cumulative PVE (Covariance PCA)",
     ylim = c(0, 1))
abline(h = 0.8, col = "red", lty = 2)

abline(h = 0.9, col = "blue", lty = 2)

plot(cumsum(pve_cor), type = "b", pch = 19,
     xlab = "Principal Component",
     ylab = "Cumulative PVE",
     main = "Cumulative PVE (Correlation PCA)",
     ylim = c(0, 1))

abline(h = 0.8, col = "red", lty = 2)
abline(h = 0.9, col = "blue", lty = 2)

par(mfrow = c(1, 1))

After scaling, we see that the first principal component explains \(62.0 \%\) of the variance in the data, the next principal component explains \(24.7 \%\) of the variance, and so forth.

Number of PCs: Kaiser rule and parallel analysis (scaled PCA)

# Kaiser rule: keep eigenvalues > 1 for correlation-based PCA
screeplot(pca_cor, type = "l", main = "Scree: correlation PCA")
abline(h = 1, col = "red", lty = 2)

## Optional: annotate
text(2, 1.1, "Kaiser threshold = 1", col = "red")

# Parallel analysis
set.seed(24620)
B <- 500
n <- nrow(X)
p <- ncol(X)
eig <- pca_cor$sdev^2

rand_eigs <- replicate(B, {
  X_rand <- matrix(rnorm(n * p), nrow = n, ncol = p)
  eigen(cor(X_rand), symmetric = TRUE, only.values = TRUE)$values
})

rand_mean <- rowMeans(rand_eigs)
rand_95 <- apply(rand_eigs, 1, quantile, probs = 0.95)

pa_cutoff <- apply(rand_eigs, 1, quantile, probs = 0.95)

plot(eig, type = "b", pch = 19,
     xlab = "Principal Component",
     ylab = "Eigenvalue",
     main = "Parallel Analysis")
lines(rand_mean, type = "b", pch = 17, lty = 1, col = 2)

lines(rand_95, lty = 2, col = 4)

• Kaiser suggests keeping 1 PC.
• Parallel analysis (95% baseline) also suggests keeping 1 PC.

Loadings interpretation (scaled PCA)

round(pca_cor$rotation[, 1:2], 3)

            PC1    PC2
Murder   -0.536 -0.418
Assault  -0.583 -0.188
UrbanPop -0.278  0.873
Rape     -0.543  0.167

par(mfrow = c(1, 2))
barplot(pca_cor$rotation[, 1], las = 2, main = "Loadings: PC1", col = "steelblue")
barplot(pca_cor$rotation[, 2], las = 2, main = "Loadings: PC2", col = "tomato")

• PC1 represents overall criminal level.
• PC2 contrasts urbanization with violent crime, especially murder.

Score map (states in first two PCs)

The score matrix is stored in pca$x. The \(k\)th column is the \(k\)th principal component score vector.

scores <- pca_cor$x[, 1:2]
plot(scores, type = "n", xlab = "PC1", ylab = "PC2",
     main = "USArrests states projected to (PC1, PC2)")
text(scores, labels = rownames(X), cex = 0.65, col = "gray30")

abline(h = 0, v = 0, lty = 2, col = "gray70")

Score map + loadings

We can have a more informative visualization by adding the loading directions

biplot(pca_cor, scale = 0, cex = 0.7)

• PC1 represents overall criminal level, states on the left have high criminal rates overall.
• PC2 contrasts urbanization with violent crime: states at the top are more urban, while those at the bottom tend to have higher murder rates relative to their overall crime profile.

Orthogonality check

round(cov(pca_cor$x), 3)

     PC1  PC2   PC3   PC4
PC1 2.48 0.00 0.000 0.000
PC2 0.00 0.99 0.000 0.000
PC3 0.00 0.00 0.357 0.000
PC4 0.00 0.00 0.000 0.173

Wrap-up

• Scaling matters here because variable magnitudes are very different.
• Use multiple criteria to choose PCs (scree/PVE + Kaiser + parallel analysis).
• Interpret loadings (variables) together with score map (states).