# Load the data
darlingtonia <- read_csv("data/darlingtonia.csv")
# View the structure
head(darlingtonia)
## # A tibble: 6 × 11
## site height mouth_diam tube_diam keel_diam wing1_length wing2_length
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 TJH 654 38.4 16.6 6.4 85 76
## 2 TJH 413 22.2 17.2 5.9 55 26
## 3 TJH 610 31.2 19.9 6.7 62 60
## 4 TJH 546 34.4 20.8 6.3 84 79
## 5 TJH 665 30.5 20.4 6.6 60 51
## 6 TJH 665 33.6 19.5 6.6 84 66
## # ℹ 4 more variables: wingsprea <dbl>, hoodmass_g <dbl>, tubemass_g <dbl>,
## # wingmass_g <dbl>
darlingtonia %>%
count(site, name = "n_plants")
## # A tibble: 4 × 2
## site n_plants
## <chr> <int>
## 1 DG 25
## 2 HD 12
## 3 LEH 25
## 4 TJH 25Lecture: Multivariate Statistics
Multivariate intro
Introduction to Multivariate Analysis
Why Study Multivariate Data?
- In ecology, often measure many variables on the same organisms:
- Plant morphology: height, width, leaf size, flower size
- Environmental conditions: temperature, pH, nutrients
- Species composition: presence/absence, abundance
- Challenge: How do we make sense of all these variables at once?
- PCA helps by:
- Reducing many correlated variables to a few key dimensions
- Visualizing patterns in high-dimensional data
- Identifying the most important sources of variation
Note
Gotelli & Ellison, A Primer of Ecological Statistics, Ch. 12 — The Analysis of Multivariate Data, is the core reference for this lecture.
Darlingtonia californica (Cobra Lily)
Our Study System: Darlingtonia
The Cobra Lily Dataset
- Analyze morphological measurements from 89 individual Darlingtonia californica plants collected from 4 sites in the Siskiyou Mountains of Oregon.
- Variables measured (10 total):
- Morphological measurements (mm):
- height - Height of pitcher
- mouth_diam - Mouth diameter
- tube_diam - Tube diameter
- keel_diam - Keel diameter
- wing1_length, wing2_length
- Wing lengths (2 measurements)
- wingsprea - Wing spread
- Biomass measurements (g):
hoodmass_g - Hood mass
tubemass_g - Tube mass
wingmass_g - Wing mass
- Morphological measurements (mm):
Darlingtonia pitcher structure showing measured features
- Sites:
- TJH: T.J. Howell’s Fen
- DG: Day’s Gulch
- LEH: L.E. Hunter’s Fen
- HD: High Divide
Standardizing the Data
Step 2: Z-score Transformation
Visualizing Standardized Data
Understanding Correlations
Why Correlations Matter for PCA
- PCA works best when variables are correlated:
- If variables are uncorrelated → PCA provides no benefit
- If variables are highly correlated → PCA can reduce dimensions effectively
- In data, expect correlations because:
- Larger plants have larger measurements overall
- Wing dimensions likely correlate with each other
- Mass measurements correlate with size measurements
- The correlation matrix shows relationships between all pairs of variables.
- Strong correlations (positive or negative) suggest PCA will be useful!
Strong positive correlations (dark blue) suggest variables measure similar aspects of plant size/shape.
What is Principal Component Analysis?
The Concept of PCA
- PCA creates new variables (components) that are:
- Linear combinations of original variables
- Uncorrelated with each other
- Ordered by importance (amount of variance explained)
- The mathematical idea:
- For each plant i and component k:
- \(Z_{ik} = a_{i1}Y_1 + a_{i2}Y_2 + ... + a_{in}Y_n\)
- Where:
- \(Z_{ik}\) = score on principal component k
- \(Y_j\) = standardized value of original variable j
- \(a_{ij}\) = loading (weight) of variable j on component i
- Goal: Find loadings that maximize variance in first component, then second, etc.
Conceptual diagram of PCA: reducing 2 variables to 1 principal axis
- The first principal axis passes through maximum variance.
- Key properties: - As many components as variables - First components explain most variation - Later components explain residual variation
Scree Plot: Choosing Components
How Many Components to Use?
- A scree plot shows variance explained by each component.
- How to read it:
- Look for an “elbow” - sharp change in slope
- Keep components before the elbow
- Discard components in the “scree” (rubble)
- In this data:
- PC1 explains 49.9% of variance
- PC2 explains 17.7%
- PC3 explains 15.3%
- First 3 PCs explain 77.3% of variance
- Decision: Use PC1, PC2, and PC3 for interpretation.
- remaining 7 components explain 17.1% (mostly noise).
Interpreting PC1: Overall Size
- PC1 Interpretation - Looking at PC1 loadings:
- Key patterns:
- Most variables load positively and similarly: - Height: 0.310 - Mouth diameter: 0.400
- Wing lengths: 0.386, 0.372
- Wing spread: 0.256
- Biomass: 0.397, 0.380, 0.230
- Exception: Keel diameter (-0.177) and tube diameter (~0)
- Biological meaning: PC1 represents overall plant size
- most measurements scale together (bigger plants have bigger everything).
## # A tibble: 10 × 2
## Variable Loading
## <chr> <dbl>
## 1 wingmass_g -0.393
## 2 hoodmass_g -0.385
## 3 wing1_length -0.375
## 4 mouth_diam -0.369
## 5 tubemass_g -0.355
## 6 wing2_length -0.342
## 7 height -0.278
## 8 wingsprea -0.239
## 9 keel_diam 0.202
## 10 tube_diam 0.0402- High PC1 → large plants
- Low PC1 → small plants
Interpreting PC2: Shape Contrast
PC2 Interpretation
- Looking at PC2 loadings:
- Key patterns:
- Height: -0.419 (negative, moderate)
- Mouth: -0.250 (negative, moderate)
- Wing measurements: positive (0.275-0.434)
- Wing spread: 0.434 (positive, strong)
- Biomass: mostly negative
- Biological meaning: PC2 represents shape trade-offs
- tall narrow pitchers vs. short wide pitchers with large wings.
## # A tibble: 10 × 2
## Variable Loading
## <chr> <dbl>
## 1 tube_diam -0.602
## 2 keel_diam -0.491
## 3 height -0.449
## 4 tubemass_g -0.320
## 5 wing1_length 0.171
## 6 wingsprea 0.156
## 7 wing2_length 0.131
## 8 mouth_diam -0.110
## 9 wingmass_g 0.0800
## 10 hoodmass_g -0.0575- High PC2 → short plants with large wings
- Low PC2 → tall narrow plants
Interpreting PC3
PC3 Interpretation
- Looking at PC3 loadings:
- Key patterns:
- Dominated by:
- Tube diameter: 0.743 (very strong positive)
- Keel diameter: 0.576 (strong positive)
- Most other variables near zero
- Biological meaning: PC3 represents tube girth
- width/diameter of tube structure independent of height or wing size.
- Plants with high PC3 scores = fat tubes, allowing to trap larger prey?
## # A tibble: 10 × 2
## Variable Loading
## <chr> <dbl>
## 1 wingsprea -0.549
## 2 wing2_length -0.415
## 3 tube_diam -0.390
## 4 keel_diam -0.328
## 5 wing1_length -0.280
## 6 tubemass_g 0.263
## 7 mouth_diam 0.214
## 8 height 0.198
## 9 hoodmass_g 0.185
## 10 wingmass_g -0.00743- High PC3 → fat tubes
- Low PC3 → skinny tubes
Visualizing PCA Results: Scores Plot
Ordination of Sites
- We can plot PC scores to visualize:
- How plants are distributed in PC space
- Whether sites differ in morphology
- Patterns and groupings
- In this plot:
- Each point = one plant
- Position determined by PC1 and PC2 scores
- Colors indicate collection site
- Observations:
- Sites show some separation
- HD (orange) plants have lower PC1 (smaller)
- TJH (pink) plants variable but lower PC2
- Some overlap among sites
- This is called ordination - ordering objects along axes.
The ordination
PCA Biplot: Combining Scores and Loadings
Reading a Biplot
- A biplot shows both:
- Points = plants (PC scores)
- Arrows = original variables (loadings)
- How to interpret arrows:
- Length = how important the variable is
- Direction = how the variable relates to PCs
- Angle between arrows = correlation
- Small angle = positively correlated
- 90° = uncorrelated
- 180° = negatively correlated
- How to interpret points:
- Plants in the direction of an arrow are high in that variable
- Plants opposite an arrow are low in that variable
The plot
Visualizing Site Differences
Box Plots of PC Scores
HD plants are significantly smaller than plants at other sites.
TJH plants have different shapes
- taller and narrower with smaller wings.
![Boxplot of PC2 (shape) scores by site with jittered points overlaid, showing the TJH site's plants scoring differently, consistent with taller, narrower pitchers with smaller wings.]()
Visualization of the distances
# 1. Prepare PC Scores with Site
pc_scores_for_dist <- as.data.frame(pca_result$x) %>%
mutate(site = darlingtonia$site)
# 2. Calculate the mean PC scores for each site (the site "centroids")
site_centroids <- pc_scores_for_dist %>%
group_by(site) %>%
summarise(across(starts_with("PC"), mean)) %>%
ungroup()
# 3. Set site names as row names BEFORE calculating distances
site_centroids_matrix <- site_centroids %>%
column_to_rownames("site") # This keeps site names as row names
# 4. Calculate the Euclidean distance matrix between site centroids
dist_matrix <- dist(site_centroids_matrix, method = "euclidean")
# Display the distance matrix
dist_matrix
## DG HD LEH
## HD 4.346434
## LEH 1.621198 3.659711
## TJH 1.843120 3.021386 1.960547
# Visualize the distance matrix
fviz_dist(dist_matrix,
gradient = list(low = "#00AFBB", high = "#FC4E07"))
Additional PCA Diagnostics
Checking Assumptions and Quality
- Sampling Adequacy
- Kaiser-Meyer-Olkin (KMO) Test:
- Measures sampling adequacy
- Values > 0.6 = acceptable
- Values > 0.8 = good
Our overall KMO = 0.80 (good!)
# KMO test requires psych package
KMO(cor(darl_numeric))
## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = cor(darl_numeric))
## Overall MSA = 0.76
## MSA for each item =
## height mouth_diam tube_diam keel_diam wing1_length wing2_length
## 0.64 0.84 0.47 0.61 0.85 0.80
## wingsprea hoodmass_g tubemass_g wingmass_g
## 0.76 0.78 0.72 0.862. Outlier Detection
Check for influential observations:
- Functional PCA: - For time series or spatial data - Treats curves as observations - Common in physiology
- Probabilistic PCA: - Includes uncertainty estimation - Can handle missing data - Maximum likelihood framework
- Kernel PCA: - Handles non-linear relationships - Uses kernel trick from machine learning - More flexible but harder to interpret
- These are beyond our scope but worth knowing exist!
❌ Using PCA on categorical data
→ ✓ Use appropriate method (CA, MCA)
❌ Not checking for outliers
→ ✓ Plot scores, check for extreme values
❌ Assuming PCs have biological meaning
→ ✓ Interpret carefully, components are mathematical
❌ Using PCA with too few samples
→ ✓ Need n > 3p (preferably n > 50)
❌ Not reporting variance explained
→ ✓ Always report % variance for each PC
❌ Forgetting that PCA is exploratory
→ ✓ Use for patterns, not hypothesis testing
Alternative PCA Visualization: FactoMineR
# Create dataset with row names
darl_pca_data <- darlingtonia %>%
select(height, mouth_diam, tube_diam, keel_diam,
wing1_length, wing2_length, wingsprea,
hoodmass_g, tubemass_g)
# PCA
pca_facto <- PCA(darl_pca_data,
graph = FALSE,
scale.unit = TRUE)
# Create grouping factor with same length
site_groups <- factor(darlingtonia$site)
# Plot with grouping
fviz_pca_biplot(pca_facto,
col.ind = site_groups, # Use col.ind instead
palette = "jco",
addEllipses = TRUE,
ellipse.level = 0.95,
repel = TRUE)
Quick PCA Plot with ggfortify
PCA with vegan: rda() and Distance Metrics
# Prepare data matrix (samples × variables)
darl_matrix <- darlingtonia %>%
select(height, mouth_diam, tube_diam, keel_diam,
wing1_length, wing2_length, wingsprea,
hoodmass_g, tubemass_g) %>%
scale()
# PCA with vegan (uses SVD, similar to prcomp)
pca_vegan <- rda(darl_matrix)
# Quick plot
biplot(pca_vegan, scaling = 2)
# Calculate distances between sites
site_matrix <- darlingtonia %>%
select(height, mouth_diam, wingmass_g) %>%
mutate(site = darlingtonia$site) %>%
group_by(site) %>%
summarise(across(everything(), mean)) %>%
column_to_rownames("site")
# Euclidean distance (many distance metrics available!)
dist_euclidean <- vegdist(site_matrix, method = "euclidean")
dist_bray <- vegdist(site_matrix, method = "bray") # Bray-Curtis
dist_manhattan <- vegdist(site_matrix, method = "manhattan")
# Create a nice color palette
site_factor <- factor(darlingtonia$site)
n_sites <- length(levels(site_factor))
colors <- brewer.pal(min(n_sites, 8), "Set1") # Use Set1 palette
site_colors <- colors[as.numeric(site_factor)]
# Ordination plot with sites
ordiplot(pca_vegan, type = "none")
points(pca_vegan, display = "sites", col = site_colors, pch = 19, cex = 1.5)
ordispider(pca_vegan, groups = darlingtonia$site, col = colors)
ordiellipse(pca_vegan, groups = darlingtonia$site, draw = "polygon",
col = colors, alpha = 0.2)
# Add legend
legend("topright", legend = levels(site_factor),
col = colors, pch = 19, bty = "n")
Testing Group Differences: PERMANOVA
# Analysis
pca_vegan <- rda(darl_matrix)
# darl_matrix is SCALED, so it contains negative values. adonis2() defaults
# to Bray-Curtis, which is undefined for negatives -- R warns "results may be
# meaningless". For standardized continuous data the right choice is Euclidean.
permanova <- adonis2(darl_matrix ~ site, data = darlingtonia,
method = "euclidean")
# Visualization
fviz_pca_biplot(pca_result,
habillage = darlingtonia$site,
addEllipses = TRUE) +
labs(title = paste0("PERMANOVA: p = ",
round(permanova$`Pr(>F)`[1], 3)))
# Scale the data properly
darl_scaled <- darlingtonia %>%
select(height, mouth_diam, tube_diam, keel_diam,
wing1_length, wing2_length, wingsprea,
hoodmass_g, tubemass_g) %>%
scale() %>%
as.data.frame() # Convert back to data frame
# PERMANOVA
permanova_scaled <- adonis2(darl_scaled ~ site,
data = darlingtonia,
method = "euclidean",
permutations = 999)
permanova_scaled
## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = darl_scaled ~ site, data = darlingtonia, permutations = 999, method = "euclidean")
## Df SumOfSqs R2 F Pr(>F)
## Model 3 189.72 0.24512 8.9836 0.001 ***
## Residual 83 584.28 0.75488
## Total 86 774.00 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ordination with Convex Hulls
# Ordination plot with convex hulls (simpler)
ordiplot(pca_vegan, type = "none")
ordihull(pca_vegan, groups = darlingtonia$site,
draw = "polygon", col = 1:n_sites, alpha = 50)
points(pca_vegan, display = "sites", pch = 21,
bg = as.numeric(factor(darlingtonia$site)), cex = 1.5)
legend("topright", legend = levels(factor(darlingtonia$site)),
pt.bg = 1:n_sites, pch = 21, bty = "n")