Activity: Principal Component Analysis

Ordination

Hands-on activity: principal component analysis on the classic iris flower dataset — assumptions, standardization, loadings, scree plots, scores, and a manually-built biplot.

Author

Bill Perry

Worksheet: Principal Component Analysis

How to use this worksheet

Work through each part in order, at your own pace. Type every line of code yourself into a plain R script — do not copy-paste. Blocks marked ▶ Run this are code you should type and execute. Blocks marked ✏️ Your turn ask you to write, modify, or answer something. Boxes marked 🚀 If you finish early are optional bonus material that goes a bit further than what we covered in lecture.

This worksheet uses the classic iris flower dataset — 4 measurements (sepal length/width, petal length/width) on 3 species of Iris.

▶ Run this in your Script:

library(vegan)
library(corrplot)
library(janitor)
library(patchwork)
library(tidyverse)

What is PCA? PCA (Principal Component Analysis) is a technique to:

Reduce the number of variables in your dataset
Find patterns in high-dimensional data
Create new uncorrelated variables (principal components) from correlated original variables
Visualize complex relationships in multivariate data

When to use PCA: you have multiple continuous variables that may be correlated, too many variables to analyze or visualize easily, need to reduce dimensionality while retaining most information, or want to explore patterns in multivariate data.

Key assumptions of PCA:

Linear relationships between variables
No extreme outliers (can distort results)
Variables should be correlated (if not, PCA won’t reduce dimensions effectively)
Adequate sample size (generally n > 50)
Consider standardization when variables have different scales

Important

Critical first step: Always standardize your data when variables are measured on different scales. This prevents variables with larger values from dominating the analysis.

Part 1 · Iris Data Overview

We’ll analyze the famous iris dataset with measurements from three species: Iris setosa, Iris versicolor, and Iris virginica. Each flower has 4 measurements: sepal length, sepal width, petal length, and petal width.

▶ Run this:

iris_df <- read.csv("data/iris.csv") %>%
  clean_names() %>%
  mutate(ind = row_number()) %>%
  mutate(species_ind = paste(species, ind, sep="_"))

head(iris_df)

# Get numeric values only for PCA
iris_data_df <- iris_df %>%
  dplyr::select(sepal_length, sepal_width, petal_length, petal_width)

# Keep species info for later
iris_species_df <- iris_df %>%
  dplyr::select(species, ind, species_ind)

▶ Run this — visualize the raw measurements by species:

iris_long_df <- iris_df %>%
  pivot_longer(
    cols = c(sepal_length, sepal_width, petal_length, petal_width),
    names_to = "variable",
    values_to = "values"
  )

overview_plot <- iris_long_df %>%
  ggplot(aes(species, values, color = species)) +
  geom_boxplot() +
  facet_wrap(~variable, scales = "free") +
  labs(title = "Iris Measurements by Species",
       x = "Species",
       y = "Measurement Value") +
  theme_minimal(base_size = 7)

overview_plot

Part 2 · Check PCA Assumptions

▶ Run this — check correlations between variables:

cor_matrix <- cor(iris_data_df)
cor_matrix

corrplot(cor_matrix, method = "color", type = "upper",
         addCoef.col = "grey55", tl.cex = 0.8, number.cex = 0.8,
         title = "Correlation Matrix of Iris Variables",
         mar = c(0, 0, 2, 0))

▶ Run this — check for outliers:

outlier_plot <- iris_data_df %>%
  pivot_longer(everything(), names_to = "variable", values_to = "value") %>%
  ggplot(aes(x = variable, y = value)) +
  geom_boxplot() +
  labs(title = "Check for Outliers in Iris Variables",
       x = "Variable",
       y = "Value") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

outlier_plot

✏️ Your turn: Are the iris variables sufficiently correlated for PCA to be useful? Do you see any concerning outliers? ________________________

Part 3 · Standardize the Data

Since our variables have different scales (e.g., petal width ranges 0.1-2.5 while sepal length ranges 4-8), we need to standardize.

▶ Run this:

iris_scaled <- scale(iris_data_df)

iris_scaled_df <- as.data.frame(iris_scaled)

# Check standardization worked
colMeans(iris_scaled_df)
apply(iris_scaled_df, 2, sd)

Part 4 · Perform PCA

▶ Run this:

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

summary(iris_pca_model)

Part 5 · Extract and Understand Results

▶ Run this — eigenvalues and variance explained:

eigenvalues <- iris_pca_model$sdev^2
prop_variance <- eigenvalues / sum(eigenvalues)
cumsum_variance <- cumsum(prop_variance)

pca_summary_df <- data.frame(
  Component = paste0("PC", 1:length(eigenvalues)),
  Eigenvalue = eigenvalues,
  Prop_Variance = prop_variance,
  Cumsum_Variance = cumsum_variance
)

pca_summary_df

▶ Run this — component loadings:

loadings_df <- data.frame(
  Variable = rownames(iris_pca_model$rotation),
  PC1 = iris_pca_model$rotation[, 1],
  PC2 = iris_pca_model$rotation[, 2],
  PC3 = iris_pca_model$rotation[, 3],
  PC4 = iris_pca_model$rotation[, 4]
)

loadings_df

✏️ Your turn: Which variable has the largest loading on PC1? On PC2? ________________________

Part 6 · Determine Number of Components

▶ Run this — a scree plot:

scree_data_df <- data.frame(
  Component = factor(1:4),
  Variance = prop_variance * 100
)

scree_plot <- ggplot(scree_data_df, aes(x = Component, y = Variance)) +
  geom_col(fill = "darkblue") +
  geom_line(aes(group = 1), linewidth = 1) +
  geom_point(size = 3) +
  labs(title = "Scree Plot - Variance Explained by Each Component",
       x = "Principal Component",
       y = "% of Variance Explained") +
  theme_minimal()

scree_plot

▶ Run this — apply decision rules:

# Eigenvalue > 1 rule
components_to_keep <- sum(eigenvalues > 1)
components_to_keep

# Components explaining at least 80% variance
components_80_percent <- which(cumsum_variance >= 0.80)[1]
components_80_percent

Part 7 · Create PCA Scores

▶ Run this:

pca_scores_df <- data.frame(
  iris_pca_model$x,
  species = iris_species_df$species
)

head(pca_scores_df)

Part 8 · Visualize Results

▶ Run this — a PCA scores plot:

scores_plot <- ggplot(pca_scores_df, aes(x = PC1, y = PC2, color = species)) +
  geom_point(size = 3, alpha = 0.7) +
  stat_ellipse(level = 0.68, linetype = 2) +
  labs(title = "PCA Scores Plot",
       subtitle = paste0("PC1 explains ", round(prop_variance[1]*100, 1),
                        "% of variance, PC2 explains ", round(prop_variance[2]*100, 1), "%"),
       x = paste0("PC1 (", round(prop_variance[1]*100, 1), "%)"),
       y = paste0("PC2 (", round(prop_variance[2]*100, 1), "%)"),
       color = "Species") +
  theme_minimal() +
  scale_color_manual(values = c("#00AFBB", "#E7B800", "#FC4E07"))

scores_plot

▶ Run this — a loading plot (arrows only):

loading_data_df <- loadings_df %>%
  dplyr::select(Variable, PC1, PC2)

loading_plot <- 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),
            size = 4) +
  xlim(-1, 1) + ylim(-1, 1) +
  labs(title = "PCA Loading Plot",
       x = paste0("PC1 (", round(prop_variance[1]*100, 1), "%)"),
       y = paste0("PC2 (", round(prop_variance[2]*100, 1), "%)")) +
  theme_minimal() +
  geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5)

loading_plot

▶ Run this — combine scores and loadings into a manual biplot:

arrow_scale <- 3

biplot_plot <- ggplot(pca_scores_df, aes(x = PC1, y = PC2)) +
  geom_point(aes(color = species), size = 2, alpha = 0.6) +
  geom_segment(data = loading_data_df,
               aes(x = 0, y = 0,
                   xend = PC1 * arrow_scale,
                   yend = PC2 * arrow_scale),
               arrow = arrow(length = unit(0.3, "cm")),
               color = "black", linewidth = 0.8) +
  geom_text(data = loading_data_df,
            aes(x = PC1 * arrow_scale * 1.1,
                y = PC2 * arrow_scale * 1.1,
                label = Variable),
            size = 3) +
  stat_ellipse(aes(color = species), level = 0.68, linetype = 2) +
  labs(title = "PCA Biplot - Iris Dataset",
       subtitle = "Points = Individual flowers, Arrows = Original variables",
       x = paste0("PC1 (", round(prop_variance[1]*100, 1), "%)"),
       y = paste0("PC2 (", round(prop_variance[2]*100, 1), "%)"),
       color = "Species") +
  theme_minimal() +
  scale_color_manual(values = c("#00AFBB", "#E7B800", "#FC4E07"))

biplot_plot

✏️ Your turn: In the biplot, which species has the largest (most positive) PC1 scores? Which arrow(s) point in roughly the same direction as that species’ cluster? ________________________

Tip

🚀 If you finish early: arrow_scale <- 3 was chosen somewhat arbitrarily to make the arrows visible against the points. Try changing it to 1 and to 6 and re-running the biplot. What changes, and what stays the same?

# Write your code here:

Part 9 · Interpret Results

▶ Run this:

pc1_loadings <- iris_pca_model$rotation[, 1]
pc1_loadings

pc2_loadings <- iris_pca_model$rotation[, 2]
pc2_loadings

total_variance_2pc <- sum(prop_variance[1:2])
total_variance_2pc

✏️ Your turn: Write a one-sentence biological interpretation of PC1 and of PC2, based on the loadings you just examined.

PC1 = ________________________________________________
PC2 = ________________________________________________

Summary checklist for PCA

When conducting PCA, always follow these steps:

Tip

PCA checklist

Explore your data - check distributions and relationships.
Check correlations - PCA works best with correlated variables.
Check for outliers - they can distort results.
Standardize if needed - essential when variables have different scales.
Run PCA - extract components.
Determine number of components - use scree plot and variance explained.
Interpret loadings - understand what each component represents.
Visualize results - create scores plots and biplots.
Validate interpretation - ensure it makes biological sense.

Key points to remember:

PCA finds new variables (components) that are linear combinations of original variables.
Components are uncorrelated and ordered by variance explained.
Standardization is crucial when variables have different units/scales.
First few components usually capture most variation.
Loadings show how original variables contribute to components.
Scores show where observations fall in the new component space.

Important

Key points from PCA analysis

Check assumptions first - especially correlations and outliers.
Standardize when necessary - prevents scale effects from dominating.
Use multiple criteria to decide number of components (scree plot, eigenvalue > 1, variance explained).
Interpret components based on loadings - what do they represent biologically?
Visualize in 2D using first two components if they explain sufficient variance.
PCA is exploratory - use it to understand patterns, not for hypothesis testing.
Document your choices - why you kept certain components, how you interpreted them.

Remember: PCA is a dimension reduction technique - the goal is to simplify complex data while retaining the important patterns!

Review and checkpoint

At this point you can:

Check PCA assumptions (correlations, outliers, sample size)
Standardize variables before running PCA
Run prcomp() and interpret its summary output
Read a scree plot and apply decision rules for component retention
Interpret loadings to give principal components biological meaning
Build a scores plot, a loading plot, and a combined biplot

Note

📤 What to turn in before next class

Upload both of these to the course management system:

Your code — the scripts/ folder (or just 09_pca.R)
This worksheet, with your written answers

Getting unstuck

When code breaks — and it will, that is normal:

Read the error message out loud. R usually names the line and the problem.
Check the usual suspects: did you run all the library() calls at the top? Spelling? A missing ) or %>% at the start of a line?
?function_name opens the built-in help page.
Bring the exact error (copy-paste it) to class or office hours.

💡 Key idea: Every working scientist googles error messages daily. Getting stuck is not failing — it is the job.

Other Formats

Worksheet: Principal Component Analysis

How to use this worksheet

Part 1 · Iris Data Overview

Part 2 · Check PCA Assumptions

Part 3 · Standardize the Data

Part 4 · Perform PCA

Part 5 · Extract and Understand Results

Part 6 · Determine Number of Components

Part 7 · Create PCA Scores

Part 8 · Visualize Results

Part 9 · Interpret Results

Summary checklist for PCA

Review and checkpoint

Getting unstuck