Lecture: Principal Component Analysis

Ordination by eigen-decomposition

Principal Component Analysis (PCA) in full — what prcomp() is actually solving, eigenvalues and eigenvectors, axis rotation, standardization, scree plots and component-retention rules, loadings, the arbitrary sign of a component, scores and biplots, worked step-by-step on the classic iris flower dataset.

Part 1 · Where We Left Off

Last Week: The Multivariate Toolkit

Three steps, in order

1. Standardize — because units are not biology.
2. Choose a distance — Euclidean for continuous traits, Bray–Curtis for species counts, Jaccard for presence–absence. This is the modelling decision.
3. Ordinate — draw the resemblances on a readable 2-D map.

Where PCA sits

• PCA is the eigen-decomposition branch of that map: it works on a correlation or covariance matrix, implies Euclidean distance, and assumes linear relationships.
• In exchange it gives you ordered axes, % variance explained, and loadings — which NMDS cannot.

Example NMDS ordination plot from a published lake-ecology study (‘Summer’), with lake samples plotted as colored points (by lake) and shapes (by geographic location) and environmental-variable vectors (dissolved P and N, mixed-layer depth, temperature) overlaid as arrows with their R² and p-values, illustrating a real multivariate ordination used to summarize many variables at once.

An ordination of the kind we are building toward.

Today’s Question

Open the box

• Last week we ran prcomp() and read the answer. Today we ask what it actually solved.
• By the end you will be able to say, precisely:
  • what an eigenvalue is and why it equals a variance
  • what an eigenvector is and why its entries are the loadings
  • why there are exactly p components for p variables
  • why component signs are arbitrary, and what to do about it
  • how to decide how many components to keep

Note

Gotelli & Ellison, A Primer of Ecological Statistics, Ch. 12 — The Analysis of Multivariate Data — is the core reference.

The core equation

For object i and component k:

\[z_{ik} = c_{k1}y_{i1} + c_{k2}y_{i2} + \cdots + c_{kp}y_{ip}\]

• \(y_{i1} \ldots y_{ip}\) — the original values for object i
• \(c_{k1} \ldots c_{kp}\) — the loadings: coefficients saying how much each original variable matters to component k
• \(z_{ik}\) — the score: object i’s position on component k
• PCA’s whole job is to find the \(c\)’s.

Today’s Data: Iris

Four measurements, three species

• 150 flowers, 4 measurements (cm), 50 from each of three species.
• Small enough to follow every number by hand, correlated enough for PCA to have something to do.
• The test we will apply: PCA is never told which species a flower belongs to. Does it recover them anyway?

iris_df <- iris %>% clean_names()

iris_data <- iris_df %>%
  dplyr::select(-species)

The Problem PCA Solves

Four dimensions is already too many

• With 4 variables there are 6 pairwise scatterplots. With 10 there are 45. With 27 there are 351.
• No set of 2-D plots shows you how all the variables move together.
• PCA’s offer: a small number of new axes, built from the old ones, that carry most of the information.
• And it will do it without being told about species.

Part 2 · Eigenvalues, Eigenvectors, and Rotation

PCA Is a Rotation

The 2-variable picture

• Start with two correlated variables — crab burrows and crab biomass. The cloud of points runs diagonally.
• The original axes are a bad description: both of them carry some of the same information.
• PCA keeps the points exactly where they are and turns the axes until the first one lies along the direction of greatest spread.
• Nothing is lost in a rotation. The distances between points are identical before and after.

Scatterplot of number of crab burrows versus total crab biomass, with the original centered axes drawn as a horizontal/vertical cross through the data cloud and two diagonal lines labeled PC1 and PC2 showing the rotated principal-component axes, where PC1 runs along the direction of greatest spread in the data.

The Same Data, New Axes

The same crab burrow and biomass data redrawn after axis rotation, with PC1 and PC2 now shown as the new horizontal and vertical axes and the original burrow-number and biomass axes appearing as diagonal reference lines, illustrating that PCA is literally a rotation of the coordinate system around the data’s center.

After the rotation

• PC1 is now horizontal and carries most of the spread; PC2 is vertical and carries what is left.
• The two new axes are orthogonal, which is the geometric way of saying the two components are uncorrelated.
• The old variables are still in there — they are now the diagonal reference lines, and the angle each one makes with PC1 is exactly what its loading records.
• The rotation is chosen to maximize the variance on PC1. That one sentence is the whole of PCA.

Beyond Two Dimensions

Same idea, more axes

• In 3-D the rotation finds a first axis through the long direction of the cloud, a second at right angles to it through the longest remaining direction, and a third fixed by the first two.
• With p variables you get exactly p components — the rotation uses up every dimension.
• You do not do this by eye. The components are extracted from the p × p covariance or correlation matrix of the original variables, and the extraction is an eigen-decomposition.

3D scatterplot of a point cloud within a cube, with three colored arrows (red ‘1st’, blue ‘2nd’, green ‘3rd’) radiating from the cloud’s center marking the first three principal-component axes, extending the rotation idea from 2D to 3D multivariate space.

Eigenvalues and Eigenvectors

Eigenvalue (\(\lambda\))

• How much variance the component captures, in the units of the input matrix.
• Run PCA on a correlation matrix and each original variable contributes exactly 1 unit of variance, so:
  • the eigenvalues sum to p
  • an eigenvalue of 1 means “this component is worth one average variable”
• % variance for component k: \(\lambda_k / \sum_{j=1}^{p}\lambda_j\)

Eigenvector (c)

• The list of coefficients — one per original variable — that builds the component. There are p entries in each eigenvector and p eigenvectors.
• Each eigenvector has length 1: the loadings’ squares sum to 1.
• Eigenvectors are orthogonal to one another: any two of them have a dot product of 0. That is why the components are uncorrelated.
• Each component has one eigenvalue and one eigenvector: how big, and which way.

Correlation Matrix or Covariance Matrix?

The choice is the standardization choice

• Covariance matrix → variables keep their original variances, so a variable measured in big numbers dominates.
  • Use it when all variables are in the same units and the differences in spread are real biology.
• Correlation matrix → every variable is rescaled to variance 1 before decomposition, so every variable gets an equal vote.
  • Use it when the units differ — almost always, in ecology.
• prcomp(x, scale. = TRUE) is “use the correlation matrix”.

Why correlation matters so much

• PCA earns its keep by collapsing redundancy. If the original variables are uncorrelated, there is no redundancy to collapse:
  • every eigenvalue comes out near 1
  • the “components” are just the original variables, reordered
  • nothing has been reduced
• Correlated variables → fewer components carry more variance. Check the correlation matrix before you run anything.

Part 3 · Assumptions and Standardization

What PCA Assumes

Check these first

1. Linear relationships between variables — PCA sees only straight lines.
2. No extreme outliers — one wild object can define a whole component.
3. Variables are correlated — otherwise there is nothing to collapse.
4. Adequate sample size — at least 3× as many objects as variables; n > 50 is comfortable. We have 150 flowers and 4 variables.
5. Complete cases — prcomp() will not handle missing values for you.

Are the iris variables correlated?

• Petal length and petal width correlate at 0.96. Plenty of redundancy to collapse.

Checking Linearity and Outliers

Relationships are close to straight lines — the setosa cluster makes some panels look bent, but within species they are linear.

A handful of mild sepal-width outliers, nothing extreme. Note how much wider petal length spreads than sepal width — that is the standardization problem in one picture.

Standardization

Why it is needed here

## sepal_length  sepal_width petal_length  petal_width 
##        0.828        0.436        1.765        0.762

• Petal length varies about 4× as much as sepal width does. On a covariance matrix it would dominate PC1 for no better reason than that.
• Z-score: \(Z = (Y_i - \bar{Y})/s\), giving mean 0 and SD 1.
• A sepal length of 5.1 cm becomes -0.9 — just under one SD below average.

iris_scaled <- scale(iris_data)

# check it worked
round(colMeans(iris_scaled), 3)
## sepal_length  sepal_width petal_length  petal_width 
##            0            0            0            0
round(apply(iris_scaled, 2, sd), 3)
## sepal_length  sepal_width petal_length  petal_width 
##            1            1            1            1

Note

scale() here, then prcomp(..., scale. = FALSE), is exactly equivalent to handing the raw data to prcomp(..., scale. = TRUE). We do it in two steps so the standardized values are visible.

Part 4 · Running PCA and Reading the Output

Running the Analysis

The call

iris_pca <- prcomp(iris_scaled,
                   center = FALSE,
                   scale. = FALSE)

• center = FALSE, scale. = FALSE because scale() already did both. Doing it twice is harmless but confusing.
• vegan::rda(iris_scaled) gives the identical decomposition under a different name — PCA is PCA.

What is inside the result

Element	What it holds
sdev	\(\sqrt{\lambda}\) — SD of each component
rotation	the eigenvectors: loadings, variables × components
x	the scores: objects × components
center, scale	what was subtracted and divided

• sdev^2 gives the eigenvalues. Everything else in this lecture is read off rotation and x.

The Summary Table

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.00000

• Standard deviation is \(\sqrt{\lambda}\), not the eigenvalue itself — square it to get the variance.
• PC1 explains 73%, PC2 22.9%; together 95.8% of all the variation in four variables.
• The four eigenvalues sum to 4 — one unit per standardized variable, as promised.

Eigenvalues in a Table

tibble(
  Component  = paste0("PC", 1:4),
  Eigenvalue = round(eigenvalues, 3),
  Percent    = round(prop_variance * 100, 2),
  Cumulative = round(cumsum_variance * 100, 2))
## # A tibble: 4 × 4
##   Component Eigenvalue Percent Cumulative
##   <chr>          <dbl>   <dbl>      <dbl>
## 1 PC1            2.92    73.0        73.0
## 2 PC2            0.914   22.8        95.8
## 3 PC3            0.147    3.67       99.5
## 4 PC4            0.021    0.52      100

Reading it

• PC1’s eigenvalue is 2.92: it is worth about 2.9 of the original variables all by itself.
• PC3 and PC4 together carry 4.2% — essentially measurement noise.
• We have reduced 4 variables to 2 numbers and kept 95.8% of the information. That is dimension reduction working.

Part 5 · How Many Components Do You Keep?

The Scree Plot

What it is

• “Scree” is the geological term for the loose rubble at the base of a cliff. The plot usually shows a cliff then rubble.
• Y-axis: % variance explained. X-axis: component number.
• The elbow method: find where the line stops dropping steeply and flattens out. Keep the components before the flattening.
• Here the cliff runs PC1 → PC2 and the rubble starts at PC3. The elbow says keep two.

Three Rules, and What To Do When They Disagree

# Kaiser criterion: eigenvalue > 1
sum(eigenvalues > 1)
## [1] 1

# first PC reaching 80% cumulative
which(cumsum_variance >= 0.80)[1]
## [1] 2

• Kaiser (\(\lambda > 1\)) says keep 1.
• 80% cumulative says keep 2.
• The elbow says keep 2.

Resolving it

• The rules are guidelines, not tests, and they routinely disagree. Kaiser is known to be strict with few variables — with only 4 variables a second component has to be very strong to clear an eigenvalue of 1, and PC2’s is 0.91.
• Two more criteria settle it in practice:
  1. Does the component have an interpretable biological meaning? PC2 does — see Part 6.
  2. Does keeping it change what you can see? It does: PC2 is what separates the species that PC1 alone confuses.
• We keep two. Say which rule you used and why.

Part 6 · Loadings and Interpretation

What Loadings Are

Recipes for the components

• A loading is one coefficient \(c_{kj}\) — how much original variable j goes into component k.
• Values run from −1 to +1, like correlations.
• |loading| > 0.5 → dominates the component. |loading| ≈ 0 → irrelevant to it.
• Signs within a component are what carry the meaning:
  • all the same sign → a size axis, everything moving together
  • mixed signs → a shape axis, a trade-off

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.524

• Each column is one eigenvector. Square a column and the entries sum to 1.

Eigenvector Properties, Checked

The three claims from Part 2

# 1. unit length
round(colSums(iris_pca$rotation^2), 3)
## PC1 PC2 PC3 PC4 
##   1   1   1   1

# 2. orthogonal: PC1 . PC2 = 0
round(sum(iris_pca$rotation[, 1] *
          iris_pca$rotation[, 2]), 10)
## [1] 0

# 3. scores are uncorrelated
round(cor(iris_pca$x[, 1],
          iris_pca$x[, 2]), 10)
## [1] 0

Why this matters practically

• Unit length is why loadings can be compared across variables within a component — they are competing for a fixed budget of 1.
• Orthogonality is why you can drop PC scores straight into an ANOVA or a regression with no multicollinearity. That is one of PCA’s most useful side effects.
• It is also why PC2 can only describe what PC1 left behind: it is constrained to be at right angles to it.

⚠️ The Sign of a Component Is Arbitrary

Expect this, do not fight it

• An eigenvector defines a line, and a line has no preferred end. Multiply a component’s loadings and its scores by −1 and you have an equally valid answer:
  • same eigenvalue, same % variance
  • same distances between points
  • same biology
• Different software, different R versions, or dropping one flower can all flip a sign. Nothing is wrong when that happens.
• The real error is to flip the loadings and not the scores, or to write “higher PC1” without saying what higher means.

# flip PC2 entirely
flip_load <- -iris_pca$rotation[, 2]
flip_score <- -iris_pca$x[, 2]

round(c(
  var_original = var(iris_pca$x[, 2]),
  var_flipped  = var(flip_score)), 4)
## var_original  var_flipped 
##        0.914        0.914

The rule: before you write a sentence about a component, translate the sign into biology. Say “larger flowers”, not “higher PC1”.

PC1: Overall Flower Size

## sepal_length  sepal_width petal_length  petal_width 
##        0.521       -0.269        0.580        0.565

The pattern

• Three variables load together and strongly: sepal length 0.52, petal length 0.58, petal width 0.56.
• Sepal width opposes them at -0.27 — moderate, and in the other direction.

The interpretation

• PC1 is an overall size axis, with a twist. On this run a high PC1 score means a flower with long sepals, long petals and wide petals, but relatively narrow sepals.
• That is not a bookkeeping quirk — it is real iris biology. The big virginica flowers are big in every dimension except that their sepals are proportionally narrower than setosa’s.
• So PC1 is better described as “large, elongated flowers vs. small, broad-sepalled flowers” than as pure size.
• It carries 73% of all the variation.

PC2: Sepal Width

## sepal_length  sepal_width petal_length  petal_width 
##       -0.377       -0.923       -0.024       -0.067

The pattern

• Sepal width dominates completely at -0.92.
• Sepal length contributes moderately (-0.38).
• The two petal measurements are effectively absent (-0.024 and -0.067).

The interpretation

• PC2 is not a general shape contrast — it is very nearly “sepal width, by itself”, plus a little sepal length.
• Because the loading is negative, on this run:
  • low PC2 → wide sepals
  • high PC2 → narrow sepals
• This is PC1’s leftovers: sepal width was the one variable PC1 could not account for, so PC2 is built almost entirely out of it.
• 22.9% of the variation — the sepal-width story is a real and separate axis of iris variation.

The Loading Plot

Reading the arrows

• Direction = which components the variable builds. Length = how well this plane represents it. The circle is the maximum possible.
• Petal length, petal width and sepal length point right — they build positive PC1 together, and their small mutual angles say they are strongly correlated.
• Sepal width points down and slightly left — nearly at right angles to the petals, which is the plot’s way of saying it carries almost independent information.
• That near-90° angle is why sepal width got its own component.

Part 7 · Scores, Biplots, and Species

Scores: Where Each Flower Lands

From loadings to positions

• A score is an object’s coordinate on a component — the result of the recipe, applied to that object’s standardized values.
• iris_pca$x holds them, 150 rows × 4 components.

## # A tibble: 9 × 5
##   species      PC1   PC2   PC3   PC4
##   <fct>      <dbl> <dbl> <dbl> <dbl>
## 1 setosa     -2.26 -0.48  0.13  0.02
## 2 setosa     -2.07  0.67  0.23  0.1 
## 3 setosa     -2.36  0.34 -0.04  0.03
## 4 versicolor  1.1  -0.86  0.68  0.03
## 5 versicolor  0.73 -0.59  0.09  0   
## 6 versicolor  1.24 -0.61  0.55  0.01
## 7 virginica   1.84 -0.87 -1    -0.05
## 8 virginica   1.15  0.7  -0.53 -0.04
## 9 virginica   2.2  -0.56  0.2   0.06

The Biplot

Points and arrows together

• Points = flowers (scores). Arrows = variables (loadings), on the same plane.
• A flower lying out along an arrow is high in that variable.
• Setosa sits far left — low PC1, so short petals and sepals, but (per PC1’s negative sepal-width loading) broad sepals.
• Virginica sits right — long petals and sepals, narrower sepals.
• Versicolor sits between, overlapping virginica slightly.
• PCA was never told about species. It recovered them from morphology alone.

Did PCA Recover the Species?

as_tibble(iris_pca$x) %>%
  mutate(species = iris_df$species) %>%
  group_by(species) %>%
  summarise(across(PC1:PC2,
                   \(z) round(mean(z), 2)))
## # A tibble: 3 × 3
##   species      PC1   PC2
##   <fct>      <dbl> <dbl>
## 1 setosa     -2.22 -0.29
## 2 versicolor  0.49  0.55
## 3 virginica   1.72 -0.26

Reading it, sign resolved

• PC1 orders the species cleanly: setosa → versicolor → virginica, smallest and broadest-sepalled to largest and narrowest-sepalled.
• PC2 does something different: versicolor stands apart from both others. Setosa and virginica have similar PC2 scores despite sitting at opposite ends of PC1.
• That is exactly the job of a second component — it describes what the first one could not. Dropping PC2 on Kaiser’s advice would have thrown away a real signal.

Part 8 · Other Tools, Reporting, and Wrap-Up

The Same PCA, Three Other Packages

FactoMineR

pca_facto <- PCA(iris_data, graph = FALSE, scale.unit = TRUE)

fviz_pca_biplot(pca_facto, col.ind = iris_df$species,
                palette = "jco", addEllipses = TRUE,
                ellipse.level = 0.95, repel = TRUE,
                geom.ind = "point", labelsize = 3, title = NULL) +
  theme_minimal(base_size = 8)

ggfortify

autoplot(iris_pca, data = iris_df, colour = "species",
         loadings = TRUE, loadings.label = TRUE,
         frame = TRUE, frame.type = "norm") +
  theme_minimal(base_size = 8)

vegan

iris_rda <- rda(iris_scaled)

biplot(iris_rda, scaling = 2)

• All four routes decompose the same matrix. The numbers are identical; only the plotting defaults and the sign conventions differ.

Reporting a PCA

Methods

“The four floral measurements were standardized to zero mean and unit variance and subjected to principal component analysis on the correlation matrix (prcomp(), R 4.6.2). Two components were retained on the basis of the scree plot and cumulative variance (> 80%).”

Results

“PC1 and PC2 together explained 95.8% of the variation (PC1 73%, PC2 22.9%). PC1 loaded positively on petal length (0.58), petal width (0.57) and sepal length (0.52), and negatively on sepal width (−0.27)…”

Always include

1. Whether you used the correlation or covariance matrix
2. % variance for every component you interpret
3. The loadings for retained components
4. How you decided how many to keep
5. The biological reading of each axis, with the sign resolved

Never include

• All p components when 2 are meaningful
• Loadings below |0.3|
• Raw eigenvalues instead of percentages
• Every individual score

When PCA Is and Is Not the Tool

Use PCA when

• Variables are continuous and correlated
• Relationships are linear
• Few zeros, and Euclidean distance is sensible
• The goal is dimension reduction or exploration
• You want uncorrelated variables for a downstream model

Use something else when

• Species composition with many zeros → NMDS or PCoA on Bray–Curtis (week 15)
• Categorical variables → CA or MCA
• You already have a dissimilarity matrix → PCoA or NMDS
• You want to classify or predict → discriminant analysis
• You want to test groups → PERMANOVA (week 15)

Summary

What PCA did

• Rotated 4 correlated axes into 4 uncorrelated ones, ordered by variance.
• PC1 (73%) — large elongated flowers vs. small broad-sepalled ones.
• PC2 (22.9%) — sepal width, almost by itself.
• Together 95.8% of the variation in 2 numbers instead of 4.
• And it separated three species it was never told about.

What to carry forward

• Eigenvalue = variance; eigenvector = loadings. On a correlation matrix the eigenvalues sum to p.
• Standardize unless the units are identical and the spreads are biology.
• Component signs are arbitrary — translate to biology before writing.
• Retention rules disagree; interpretability breaks the tie.
• Next week: when the data are counts and Euclidean distance fails — NMDS, cluster analysis, PERMANOVA and ANOSIM.

Readings

• Gotelli & Ellison, A Primer of Ecological Statistics, Ch. 12 — The Analysis of Multivariate Data
• Quinn & Keough, Experimental Design and Data Analysis for Biologists, Ch. 17 — Principal components and correspondence analysis
• Whitlock & Schluter does not cover ordination.

Data source

• Iris: Anderson, E. (1935), The irises of the Gaspé Peninsula; made famous by Fisher, R.A. (1936), Annals of Eugenics 7:179–188. Shipped with base R as iris.