Activity: Principal Component Analysis

Ordination

Hands-on activity: principal component analysis on the classic iris flower dataset — assumptions, standardization, eigenvalues and eigenvectors, component retention rules, loadings and the arbitrary sign of a component, scores, and a biplot built by hand.

Author

Bill Perry

In-class Activity 19: Principal Component Analysis

Today’s Objectives

Check the four things PCA assumes before you run it
Standardize, and say why that means “use the correlation matrix”
Run prcomp() and read every piece of what it returns
Verify by hand that eigenvectors have unit length and are orthogonal
Apply three component-retention rules, and resolve them when they disagree
Interpret loadings as biology, with the arbitrary sign resolved
Build a biplot from scratch, and understand why the arrows need scaling
Say whether PCA recovered a grouping it was never told about

📂 Files for today

📥 Script skeleton (.R)

Save the skeleton as scripts/19_pca.R.
No data file today. iris ships with base R — Anderson’s 1935 measurements of 150 flowers, made famous by Fisher in 1936.

What you hand in

All of today’s code is right here in this activity. Type it into your script as we go and run it chunk by chunk. The skeleton has the plumbing and nothing else. Every line of statistics you type yourself.
Section 7 has no code here either. By then you will have built the pieces of a biplot separately, so you assemble it on your own.
The Q boxes are things to think about, not inline homework. They are not graded.
The write-up at the bottom — sections 9, 10 and 11 — is graded. Those answers go into the script as # comments.
Hand in scripts/19_pca.R. Run it top to bottom (Ctrl/Cmd + Shift + Enter) first — it must run with no errors.

🔮 Predict before you run

PCA will not be told which species each flower belongs to. Before you plot anything, commit out loud: will the three species come apart on PC1 and PC2, or will they smear together? All three? Two of three?

The invisible(NULL) lines in the skeleton

Every section ends with invisible(NULL). It is a bumper — press Ctrl/Cmd + Enter on a section’s last real line and the cursor stops there instead of jumping past the next heading. It does nothing and prints nothing.

Keep the plots plain today. The one exception is the final figure in section 11, which is the graded deliverable.

Where this fits

Last week you ran a PCA and read the answer. Today you open the box.

	What you did	What you do today
Week 13	ran `prcomp()`, read the output	—
Today	—	eigenvalues, eigenvectors, retention rules, biplots from scratch
Week 15	—	what to do when Euclidean distance fails

Setup

▶ Run this in your Script:

# ── Setup (1) ───────────────────────────────────────────────────────────
library(vegan)
library(corrplot)
library(GGally)
library(factoextra)
library(janitor)
library(patchwork)
library(tidyverse)

theme_set(theme_minimal(base_size = 11))

Block A · Assumptions and standardization

After lecture Parts 1–3.

1 · The iris data

▶ Run this:

# ── Part 1.1 · Load and look ────────────────────────────────────────────
iris_df <- iris %>% clean_names()

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

glimpse(iris_df)

iris_df %>% count(species)

▶ Run this — the four measurements, by species:

# ── Part 1.2 · The raw variables ────────────────────────────────────────
iris_df %>%
  pivot_longer(-species, names_to = "variable", values_to = "value") %>%
  ggplot(aes(species, value, fill = species)) +
  geom_boxplot() +
  facet_wrap(~variable, scales = "free_y")

✏️ Q1

Which measurement separates the species most cleanly on its own? Which separates them least? If one variable already did the whole job, what would you need PCA for?

2 · Check the assumptions

▶ Run this — are the variables correlated?

# ── Part 2.1 · Correlations ─────────────────────────────────────────────
cor_matrix <- cor(iris_data)

round(cor_matrix, 3)

corrplot(cor_matrix, method = "color", type = "upper",
         addCoef.col = "grey25")

▶ Run this — are the relationships linear, and are there outliers?

# ── Part 2.2 · Linearity and outliers ───────────────────────────────────
ggpairs(iris_data,
        upper = list(continuous = "cor"),
        lower = list(continuous = "smooth"))

iris_data %>%
  pivot_longer(everything(), names_to = "variable", values_to = "value") %>%
  ggplot(aes(variable, value)) +
  geom_boxplot()

✏️ Q2

Petal length and petal width correlate at about 0.96. That is PCA’s best friend — explain why in one sentence, in terms of what PCA is trying to do.

✏️ Q3

Several of the ggpairs panels look slightly bent rather than straight. Look at the colours in the panel for petal length vs. sepal length. Is the bend evidence of a genuinely non-linear relationship, or is something else producing it?

3 · Standardize

▶ Run this:

# ── Part 3.1 · How different are the scales? ────────────────────────────
round(apply(iris_data, 2, sd), 3)

▶ Run this:

# ── Part 3.2 · Standardize ──────────────────────────────────────────────
iris_scaled <- scale(iris_data)

round(colMeans(iris_scaled), 3)
round(apply(iris_scaled, 2, sd), 3)

✏️ Q4

Petal length varies about four times as much as sepal width does. If you skipped standardization — ran PCA on the covariance matrix instead — which variable would dominate PC1, and would that be a biological finding or a bookkeeping accident?

Block B · Running PCA and choosing components

After lecture Parts 4–5.

4 · Run it, and look inside

▶ Run this:

# ── Part 4.1 · Perform PCA ──────────────────────────────────────────────
iris_pca_model <- prcomp(iris_scaled, center = FALSE, scale. = FALSE)

summary(iris_pca_model)

names(iris_pca_model)

▶ Run this — eigenvalues and variance explained:

# ── Part 4.2 · Eigenvalues ──────────────────────────────────────────────
eigenvalues     <- iris_pca_model$sdev^2
prop_variance   <- eigenvalues / sum(eigenvalues)
cumsum_variance <- cumsum(prop_variance)

tibble(Component  = paste0("PC", 1:4),
       Eigenvalue = round(eigenvalues, 3),
       Percent    = round(prop_variance * 100, 2),
       Cumulative = round(cumsum_variance * 100, 2))

# the eigenvalues should sum to the number of variables
sum(eigenvalues)

✏️ Q5

The four eigenvalues sum to exactly 4. Why must that be true when PCA runs on a correlation matrix? What would they sum to for the ten-variable cobra lily data from last week?

✏️ Q6

summary() reports a standard deviation for each component, not an eigenvalue. What is the relationship between the two? Check it in R.

▶ Run this — verify the eigenvector properties from lecture:

# ── Part 4.3 · Eigenvector properties ───────────────────────────────────
# 1. each eigenvector has length 1
round(colSums(iris_pca_model$rotation^2), 3)

# 2. PC1 and PC2 are orthogonal (dot product = 0)
round(sum(iris_pca_model$rotation[, 1] * iris_pca_model$rotation[, 2]), 10)

# 3. and so the scores are uncorrelated
round(cor(iris_pca_model$x[, 1], iris_pca_model$x[, 2]), 10)

✏️ Q7

Property 3 is the reason you are allowed to drop PC scores straight into a multiple regression or an ANOVA. What problem does it spare you, and why would the four original iris variables have caused that problem?

5 · How many components?

▶ Run this — a scree plot:

# ── Part 5.1 · Scree plot ───────────────────────────────────────────────
fviz_eig(iris_pca_model, addlabels = TRUE, ylim = c(0, 85))

▶ Run this — three rules:

# ── Part 5.2 · Retention rules ──────────────────────────────────────────
# Kaiser: eigenvalue > 1
sum(eigenvalues > 1)

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

✏️ Q8

Kaiser says keep one component; the 80% rule and the elbow both say keep two. These are guidelines, not tests. Give two reasons to keep PC2 anyway — one about interpretability, one about what you would be able to see.

Block C · Loadings, scores, and biplots

After lecture Parts 6–8.

6 · Loadings, and the sign trap

▶ Run this:

# ── Part 6.1 · Loadings ─────────────────────────────────────────────────
round(iris_pca_model$rotation, 3)

⚠️ Before you write a word about PC1

Look at the signs on your screen. An eigenvector defines a line, and a line has no preferred end — your loadings may come back with the opposite signs to your neighbour’s, or to the lecture slides, and both are correct.

What is not correct is writing “virginica scores higher on PC1” and stopping. Translate the sign into biology first.

✏️ Q9

On PC1, three variables load one way and sepal width loads the other. Write one sentence saying what a high PC1 score means for a flower, in plain English, with the sign on your screen resolved.

✏️ Q10

Look at PC2’s loadings. One variable dominates it almost completely and two are nearly zero. Which? Given that PC2 is built out of whatever PC1 could not explain, why is that the variable you would have predicted?

▶ Run this — a loading plot, arrows only:

# ── Part 6.2 · Loading plot ─────────────────────────────────────────────
loading_data_df <- tibble(
  Variable = rownames(iris_pca_model$rotation),
  PC1 = iris_pca_model$rotation[, 1],
  PC2 = iris_pca_model$rotation[, 2])

ggplot(loading_data_df, aes(x = 0, y = 0)) +
  geom_segment(aes(xend = PC1, yend = PC2),
               arrow = arrow(length = unit(0.3, "cm")),
               color = "red", linewidth = 1) +
  geom_text(aes(x = PC1 * 1.1, y = PC2 * 1.1, label = Variable)) +
  xlim(-1.2, 1.2) + ylim(-1.2, 1.2) +
  geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5)

✏️ Q11

Two arrows at a small angle mean the variables are correlated; two near 90° mean they are nearly independent. Which arrow sits closest to 90° from the petal arrows, and does the correlation matrix from section 2 agree?

7 · Your turn: build a biplot

There is no code for this section. You now have the two halves — the points and the arrows. Put them on one plot.

Work out and write the code for each of these, in your script under # ── Part 7.x:

Build pca_scores_df from iris_pca_model$x with the species column attached, and plot PC1 against PC2 coloured by species, with stat_ellipse(level = 0.68).
Add the loading arrows from section 6 on top of it with geom_segment(), multiplying both PC1 and PC2 by a constant arrow_scale <- 3 so the arrows are visible against the points.
Label the arrows, and give both axes labels that name the component and its percentage of variance.

Tip

If the arrows vanish into the middle of the point cloud, your arrow_scale is too small; if they shoot off the plot, it is too large.

✏️ Q12

Try arrow_scale <- 1 and then arrow_scale <- 6. What changes about the figure, and what stays exactly the same? Is the scaling constant part of the analysis, or part of the drawing?

✏️ Q13

Which species sits furthest along PC1 in the positive direction? Which arrows point that way? Using your answer to Q9, describe that species’ flowers without using the words “PC1” or “score”.

8 · Did PCA find the species?

▶ Run this:

# ── Part 8.1 · Mean scores by species ───────────────────────────────────
as_tibble(iris_pca_model$x) %>%
  mutate(species = iris_df$species) %>%
  group_by(species) %>%
  summarise(across(PC1:PC2, \(z) round(mean(z), 2)))

✏️ Q14

PC1 puts the three species in order. PC2 does something different — one species stands apart from the other two, and it is not the same one. Which species, and what does that tell you about keeping PC2 after Kaiser told you to drop it?

✏️ Q15

PCA was never given the species column. It recovered the groups from morphology alone. Does that make PCA a classification method? Explain what it would take to turn this into an actual test of whether the species differ.

Review and checkpoint

At this point you can:

Check PCA’s assumptions — correlation, linearity, outliers, sample size
Standardize, and say why that means using the correlation matrix
Run prcomp() and name every element it returns
Verify unit-length and orthogonality of eigenvectors in R
Apply three retention rules and resolve them when they disagree
Interpret loadings as biology with the sign resolved
Build a scores plot, a loading plot, and a biplot from the two
Say what an ordination can and cannot tell you about groups

Write-up — this part is graded

Type these answers into scripts/19_pca.R as # comments, in the boxes at the bottom of the skeleton.

9 · The components

Write three to four sentences. Say how many components you retained and why, what percentage of the variance they carry, and what each one means biologically — with the sign on your screen resolved. No sentence may contain the phrase “higher PC1” without saying what higher means.

10 · The sign problem

Write two to three sentences. A classmate runs the same code and gets PC1 loadings with the opposite signs to yours. Explain to them why neither of you is wrong, what would have to be true for one of you to be wrong, and what the two of you must do before comparing your written interpretations.

11 · The figure

Produce one publication-quality biplot of the iris PCA.

This one gets the full treatment — axis labels naming the component and its percentage, labelled arrows, a legend with a real title, a caption-worthy title, and a clean theme. Save it with ggsave() as a PDF.

Note

📤 What to turn in before next class

scripts/19_pca.R, running top to bottom with no errors
Your biplot from section 11, as a PDF

Other Formats

In-class Activity 19: Principal Component Analysis

Today’s Objectives

Where this fits

Setup

Block A · Assumptions and standardization

1 · The iris data

2 · Check the assumptions

3 · Standardize

Block B · Running PCA and choosing components

4 · Run it, and look inside

5 · How many components?

Block C · Loadings, scores, and biplots

6 · Loadings, and the sign trap

7 · Your turn: build a biplot

8 · Did PCA find the species?

Review and checkpoint

Write-up — this part is graded

9 · The components

10 · The sign problem

11 · The figure