Activity: Multivariate Statistics (MANOVA)

MANOVA

Hands-on activity: principal component analysis on the Darlingtonia (cobra lily) morphology dataset — standardization, correlations, PCA, loadings, scores, and distances between site centroids.

Author

Bill Perry

Worksheet: Multivariate Statistics

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 same Darlingtonia (cobra lily) morphology dataset from lecture — 89 plants measured on 10 morphological and biomass variables across 4 sites.

▶ Run this in your Script:

library(car)
library(emmeans)
library(psych)
library(vegan)
library(ggfortify)
library(corrplot)
library(FactoMineR)
library(factoextra)
library(ggrepel)
library(plotly)
library(broom)
library(patchwork)
library(janitor)
library(tidyverse)

theme_set(theme_minimal(base_size = 9))

Part 1 · Load and Explore the Data

▶ Run this:

darlingtonia <- read_csv("darlingtonia.csv")

head(darlingtonia)

darlingtonia %>%
  count(site, name = "n_plants")

▶ Run this — visualize the raw data by site (notice the different scales!):

darl_long <- darlingtonia %>%
  pivot_longer(cols = -site,
               names_to = "variable",
               values_to = "value")

ggplot(darl_long, aes(x = site, y = value, fill = site)) +
  geom_boxplot() +
  facet_wrap(~variable, scales = "free_y", ncol = 2) +
  theme_minimal(base_size = 7) +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Raw Data by Site",
       subtitle = "Notice the different scales!")

✏️ Your turn: Which variable(s) have the largest values (dominating the raw scale), and which have the smallest? ________________________

Part 2 · Standardization

▶ Run this — check means and SDs by site before standardizing:

darlingtonia %>%
  group_by(site) %>%
  summarize(
    mean_height = mean(height),
    sd_height = sd(height),
    mean_wingmass = mean(wingmass_g),
    sd_wingmass = sd(wingmass_g)
  )

# After standardization, both will have:
# mean = 0, sd = 1

▶ Run this — standardize all numeric variables and visualize:

darl_numeric <- darlingtonia %>%
  dplyr::select(-site)

# Standardize using scale() - centers (mean=0) and scales (sd=1) each variable
darl_scaled <- scale(darl_numeric)

darl_scaled_df <- as.data.frame(darl_scaled)

darl_scaled_df %>%
  select(height, mouth_diam, wingmass_g) %>%
  pivot_longer(cols = c(height, mouth_diam, wingmass_g), names_to = "variable", values_to = "value") %>%
  ggplot(aes(x = variable, y = value, fill = variable)) +
  geom_boxplot() +
  labs(title = "Standardized Data (same scale)") +
  theme(legend.position = "none")

Part 3 · Correlations

▶ Run this:

cor_matrix <- cor(darl_numeric)

corrplot(cor_matrix,
         method = "color",
         type = "upper",
         order = "hclust",
         tl.col = "black",
         tl.cex = 0.7,
         addCoef.col = "black",
         number.cex = 0.6,
         title = "Correlation Matrix\n(Darlingtonia variables)",
         mar = c(0,0,2,0))

round(cor_matrix, 2)

✏️ Your turn: Which pair of variables has the strongest correlation? Does that make biological sense? ________________________

Part 4 · Running PCA

▶ Run this:

pca_result <- prcomp(darl_numeric,
                     center = TRUE,
                     scale = TRUE)

summary(pca_result)

▶ Run this — a scree plot to decide how many components to keep:

fviz_eig(pca_result,
         addlabels = TRUE,
         ylim = c(0, 50))

✏️ Your turn: Based on the scree plot, how many components would you retain, and why? ________________________

Part 5 · Interpreting Loadings

▶ Run this — view the loadings for the first 3 PCs:

loadings <- pca_result$rotation[, 1:3]
round(loadings, 3)

▶ Run this — PC1 loadings, sorted and visualized:

pc1_loads <- tibble(
  Variable = rownames(pca_result$rotation),
  Loading = pca_result$rotation[, 1]) %>%
  arrange(desc(abs(Loading)))
pc1_loads

pc1_loads %>%
  ggplot(aes(x = reorder(Variable, Loading), y = Loading)) +
  geom_col(aes(fill = Loading > 0), show.legend = FALSE) +
  coord_flip() +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(title = "PC1 Loadings",
       subtitle = "Overall Plant Size",
       x = "Variable",
       y = "Loading") +
  scale_fill_manual(values = c("FALSE" = "coral", "TRUE" = "darkblue")) +
  theme_minimal()

▶ Run this — same for PC2:

pc2_loads <- tibble(
  Variable = rownames(pca_result$rotation),
  Loading = pca_result$rotation[, 2]
) %>%
  arrange(desc(abs(Loading)))

pc2_loads

pc2_loads %>%
  ggplot(aes(x = reorder(Variable, Loading), y = Loading)) +
  geom_col(aes(fill = Loading > 0), show.legend = FALSE) +
  coord_flip() +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(title = "PC2 Loadings",
       subtitle = "Shape: Height vs. Wing Size",
       x = "Variable",
       y = "Loading") +
  scale_fill_manual(values = c("FALSE" = "coral", "TRUE" = "darkblue")) +
  theme_minimal()

✏️ Your turn: Based on the loadings, write a one-sentence biological interpretation of PC1 and PC2 (what does a high score on each mean?).

PC1 = ________________________________________________
PC2 = ________________________________________________

Tip

🚀 If you finish early: Repeat this loadings analysis for PC3 (pca_result$rotation[, 3]). What does PC3 seem to represent?

# Write your code here:

Part 6 · Scores and Ordination

▶ Run this — view PC scores for individual plants:

pc_scores <- as_tibble(pca_result$x) %>%
  mutate(site = darlingtonia$site,
         plant_id = 1:n())

pc_scores %>%
  select(plant_id, site, PC1, PC2, PC3) %>%
  head(10)

▶ Run this — a basic scores plot:

ggplot(pc_scores, aes(x = PC1, y = PC2, color = site)) +
  geom_point(size = 3, alpha = 0.7) +
  stat_ellipse(level = 0.68, linetype = 2) +
  labs(title = "PCA Scores Plot",
       x = paste0("PC1: Overall Size"),
       y = paste0("PC2: Shapeb%)"),
       color = "Site") +
  theme_minimal(base_size = 9) +
  theme(legend.position = "right")

▶ Run this — an interactive 3D version using plotly:

pc_data <- as.data.frame(pca_result$x) %>%
  mutate(site = darlingtonia$site)

plot_ly(pc_data,
        x = ~PC1, y = ~PC2, z = ~PC3,
        color = ~site,
        type = "scatter3d",
        mode = "markers") %>%
  layout(scene = list(
    xaxis = list(title = "PC1"),
    yaxis = list(title = "PC2"),
    zaxis = list(title = "PC3")
  ))

✏️ Your turn: Looking at the scores plot, do the sites appear to separate from each other, or do they overlap heavily? ________________________

Part 7 · Distances Between Sites

▶ Run this — calculate site centroids (mean PC scores) and the Euclidean distance between them:

pc_scores_for_dist <- as.data.frame(pca_result$x) %>%
  mutate(site = darlingtonia$site)

site_centroids <- pc_scores_for_dist %>%
  group_by(site) %>%
  summarise(across(starts_with("PC"), mean)) %>%
  ungroup()

site_centroids_matrix <- site_centroids %>%
  column_to_rownames("site")

dist_matrix <- dist(site_centroids_matrix, method = "euclidean")

dist_matrix

fviz_dist(dist_matrix,
          gradient = list(low = "#00AFBB", high = "#FC4E07"))

✏️ Your turn: Which two sites are the most different from each other (largest centroid distance)? Which two are the most similar? ________________________

Tip

🚀 If you finish early: Run a quick one-way ANOVA testing whether PC1 differs by site (aov(PC1 ~ site, data = pc_scores)), just like we did in lecture. Is the site effect significant?

# Write your code here:

Review and checkpoint

At this point you can:

Explain why standardization matters before PCA
Interpret a correlation matrix in the context of PCA
Run a PCA with prcomp() and read a scree plot
Interpret loadings to give principal components biological meaning
Plot and interpret PC scores, including an interactive 3D version
Calculate Euclidean distances between group centroids in PC space

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_multivariate_statistics.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: Multivariate Statistics

How to use this worksheet

Part 1 · Load and Explore the Data

Part 2 · Standardization

Part 3 · Correlations

Part 4 · Running PCA

Part 5 · Interpreting Loadings

Part 6 · Scores and Ordination

Part 7 · Distances Between Sites

Review and checkpoint

Getting unstuck