In STA220 you used factor analysis (FA) on the HAD scale — 14 items measuring depression and anxiety. FA found two latent factors: one loading on depression items, one on anxiety. Both factors were hidden dimensions that explained the correlations among items.
PCA asks a related but different question:
| Factor Analysis (STA220) | PCA (STA230) | |
|---|---|---|
| Goal | Find latent factors that cause correlations | Find directions that capture maximum variance |
| Model | \(\mathbf{X} = \boldsymbol{\Lambda}\mathbf{f} + \boldsymbol{\varepsilon}\) (generative) | Geometric rotation of the data cloud |
| Output | Interpretable factors (depression, anxiety) | Principal components ordered by variance |
| Use case | Psychological measurement, construct validity | Compression, visualisation, multicollinearity |
Both methods reduce many correlated variables to a smaller set of dimensions. PCA is simpler and makes no distributional assumptions — it is purely about variance.
Given \(n\) observations on \(p\) variables, PCA finds a sequence of orthogonal directions (principal components) \(\mathbf{z}_1, \mathbf{z}_2, \ldots\) such that:
\[\mathbf{z}_1 = \mathbf{X}\mathbf{v}_1, \quad \text{where } \mathbf{v}_1 = \arg\max_{\|\mathbf{v}\|=1} \text{Var}(\mathbf{X}\mathbf{v})\]
Each subsequent component captures the maximum remaining variance, orthogonal to all previous ones.
In practice: standardise the variables (subtract mean, divide by SD), then compute the eigendecomposition of the covariance (or correlation) matrix:
\[\mathbf{\Sigma} = \mathbf{V} \mathbf{D} \mathbf{V}^\top\]
Proportion of variance explained by component \(k\): \(\lambda_k / \sum_j \lambda_j\)
| Output | What it tells you |
|---|---|
| Scree plot | Eigenvalues in decreasing order — look for an “elbow” |
| Loadings | How much each original variable contributes to each PC |
| Biplot | Observations (points) and variables (arrows) in PC space |
| Cumulative variance | How many PCs needed to explain e.g. 80% of variance |
Always use scale = TRUE in prcomp(). Variables on different scales (age in years, cholesterol in mmol/L) would otherwise be dominated by whichever has the largest variance — purely an artefact of measurement units, not biology.
We run PCA on five correlated health variables: Age, BMI, Waist circumference, total Cholesterol, and systolic blood pressure (SBP).
vars <- c("Age", "BMI", "Waist", "Cholesterol", "SBP")
X <- nhanes[, vars] |> na.omit() |> scale()
pca <- prcomp(X, center = FALSE, scale. = FALSE) # already scaled
summary(pca)Importance of components:
PC1 PC2 PC3 PC4 PC5
Standard deviation 1.4487 1.1570 0.9649 0.7464 0.27300
Proportion of Variance 0.4197 0.2677 0.1862 0.1114 0.01491
Cumulative Proportion 0.4197 0.6875 0.8737 0.9851 1.00000var_exp <- pca$sdev^2 / sum(pca$sdev^2)
cum_var <- cumsum(var_exp)
df_scree <- data.frame(
PC = factor(paste0("PC", 1:5), levels = paste0("PC", 1:5)),
var = var_exp,
cum = cum_var
)
ggplot(df_scree, aes(x = PC, y = var)) +
geom_col(fill = "#3498db", alpha = 0.8) +
geom_line(aes(y = cum, group = 1), colour = "#e74c3c", linewidth = 1) +
geom_point(aes(y = cum), colour = "#e74c3c", size = 3) +
geom_hline(yintercept = 0.8, linetype = "dashed", colour = "grey50") +
scale_y_continuous(labels = scales::percent_format()) +
labs(title = "Scree plot — NHANES health variables",
subtitle = "Bars: variance per PC | Red line: cumulative",
x = NULL, y = "Variance explained") +
theme_minimal()
round(pca$rotation, 3) PC1 PC2 PC3 PC4 PC5
Age -0.316 -0.591 -0.255 -0.685 0.130
BMI -0.597 0.396 0.087 0.032 0.692
Waist -0.636 0.281 0.045 -0.101 -0.710
Cholesterol -0.106 -0.370 0.922 0.025 0.003
SBP -0.358 -0.527 -0.273 0.721 -0.007biplot(pca, scale = 0, cex = c(0.4, 0.9),
col = c("grey70", "#2c3e50"),
main = "Biplot — PC1 vs PC2")
Interpretation: PC1 captures the shared “adiposity + cardiovascular burden” axis (BMI, Waist, SBP load positively). PC2 separates Age from the metabolic variables. This mirrors what FA would find, but the components are labelled by variance, not by psychological constructs.
R is loading in the background — the exercises below will be interactive once the status bar disappears.
pca$sdev^2 gives eigenvalues; divide by sum(...) for proportionspca$rotation[, 1] gives loadings on the first PCvars2 <- c("BMI", "Waist", "Cholesterol", "SBP", "Age")
dat2 <- na.omit(nhanes[, c(vars2, "Sex")])
pca2 <- prcomp(dat2[, vars2], scale. = TRUE)
scores <- as.data.frame(pca2$x[, 1:2])
scores$Sex <- dat2$Sex
library(ggplot2)
ggplot(scores, aes(PC1, PC2, colour = Sex)) +
geom_point(alpha = 0.3, size = 0.8) +
stat_ellipse() +
theme_minimal() +
labs(title = "PC1 vs PC2 by sex")Yes — this is called Principal Component Regression (PCR). It solves multicollinearity (BMI + Waist from the OLS page) because PCs are orthogonal by construction. The trade-off: PCs are harder to interpret than original variables, and there is no guarantee that the PCs that explain the most variance in X are the ones most predictive of Y. LASSO and Ridge often outperform PCR because they let the response guide variable selection.
Cronbach’s alpha measures internal consistency — whether items within a scale measure the same construct reliably (high inter-item correlation). PCA finds the direction of maximum variance regardless of reliability. On a well-designed scale like HAD, the first PC will align closely with the sum score (which Cronbach’s alpha validates). But PCA on items with low inter-item correlation still produces components — they just explain less variance and are harder to interpret as constructs.
Use the correlation matrix (i.e., scale = TRUE) when variables are on different scales — which is almost always. Using the covariance matrix means variables with larger variances dominate PC1. If you have five health variables measured in different units (years, kg/m², mmol/L, mmHg), the one with the largest raw variance will drive the first component — an artefact of measurement, not biology.
James et al. (2021). An Introduction to Statistical Learning. — Chapter 12.2