Activity: Multivariate Statistics
Distance and ordination
In-class Activity 18: Multivariate Statistics
Today’s Objectives
- Say why multivariate data cannot be handled one variable at a time
- Standardize a set of variables, and say when you should not
- Build a dissimilarity matrix and explain why it is the hub of everything downstream
- Choose between Euclidean, Manhattan, Bray–Curtis and Jaccard, and defend the choice
- Run a PCA, read a scree plot, and interpret loadings as biology
- Translate an arbitrary component sign into a correct sentence
- Run an NMDS on count data and read its stress
- Test group differences with PERMANOVA, and check its dispersion assumption
📥 Script skeleton (.R) 📥 Darlingtonia data 📥 Doubs fish counts 📥 Doubs environment
- Save the skeleton as
scripts/18_multivariate_statistics.R, and all three.csvfiles in yourdata/folder. - You will meet the Doubs data again in week 15, where NMDS gets the full treatment.
- 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 — libraries and the three data reads — and nothing else. Every line of statistics you type yourself.
- Section 7 has no code here either. By then you will have run a whole ordination workflow once, so you start the second one on your own.
- The Q boxes are things to think about, not inline homework. Read them, talk them over, answer them in your head. 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/18_multivariate_statistics.R. Run it top to bottom (Ctrl/Cmd + Shift + Enter) first — it must run with no errors.
Today has two decisions in it. First: which distance metric does each dataset need? Second: before you look at any output, will the four Darlingtonia sites separate, or overlap? Commit to an answer out loud before each one.
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 your cursor jumps to the next line that has code on it — without the bumper, that is the next section, way past the heading. With it, the cursor stops right there and you keep your place. It does nothing and prints nothing.
Keep the plots plain today — no labels, colors, or themes beyond what the code here gives you. The one exception is the final figure in section 11, which is the graded deliverable.
Where this fits
Every model from week 6 to week 12 had one response variable. Today there is no single response: a plant is its whole shape, and a river site is its whole fish community.
| What it is | Example today | |
|---|---|---|
| Weeks 6–12 | one response, variance partitioned | lifespan ~ treatment + thorax |
| Today | many responses at once, summarized | 10 traits → 3 components |
| Weeks 14–15 | the machinery, in full | PCA; NMDS and cluster analysis |
Setup
▶ Run this in your Script:
# ── Setup (1) ───────────────────────────────────────────────────────────
library(vegan)
library(psych)
library(emmeans)
library(corrplot)
library(factoextra)
library(ggrepel)
library(tidyverse)
# Permutation tests resample at random. Fixing the seed makes your
# numbers reproducible -- yours will match your neighbour's.
set.seed(42)
theme_set(theme_minimal(base_size = 11))Block A · Data, standardization, and distance
After lecture Parts 1–3.
1 · Load and look
▶ Run this:
# ── Part 1.1 · Load the two datasets ────────────────────────────────────
darlingtonia <- read_csv("data/darlingtonia.csv")
doubs_spe <- read_csv("data/DoubsSpe.csv") %>% rename(site = 1)
# das is distance from the river's source, so it orders the
# sites upstream -> downstream. We cut it into three reaches.
doubs_env <- read_csv("data/DoubsEnv.csv") %>%
rename(site = 1) %>%
mutate(reach = factor(
case_when(das < 100 ~ "Upstream",
das < 250 ~ "Midstream",
TRUE ~ "Downstream"),
levels = c("Upstream", "Midstream", "Downstream")))
darlingtonia %>% count(site)▶ Run this — the raw variables, on their own scales:
# ── Part 1.2 · The scale problem, in a picture ──────────────────────────
darlingtonia %>%
pivot_longer(-site, names_to = "variable", values_to = "value") %>%
ggplot(aes(site, value, fill = site)) +
geom_boxplot() +
facet_wrap(~variable, scales = "free_y")Which variable has the largest numbers? Which has the smallest? If you computed a distance on these raw values, which variable would decide the answer almost single-handedly?
2 · Standardize
▶ Run this:
# ── Part 2.1 · How bad is the scale problem? ────────────────────────────
darlingtonia %>%
summarize(sd_height = sd(height),
sd_wingmass = sd(wingmass_g),
ratio = sd(height) / sd(wingmass_g))▶ Run this — z-score every variable:
# ── Part 2.2 · Standardize ──────────────────────────────────────────────
darl_numeric <- darlingtonia %>% dplyr::select(-site)
darl_scaled <- scale(darl_numeric)
round(colMeans(darl_scaled), 3)
round(apply(darl_scaled, 2, sd), 3)The means came back as 0 and the SDs as 1 — by construction, not by luck. So what information did standardizing destroy, and why are we willing to lose it here?
3 · Distance, and the choice of metric
▶ Run this — the distance between two plants, two ways:
# ── Part 3.1 · Euclidean distance between two objects ───────────────────
sqrt(sum((darl_scaled[1, ] - darl_scaled[2, ])^2))
as.matrix(dist(darl_scaled))[1, 2]▶ Run this — the corner of the full dissimilarity matrix:
# ── Part 3.2 · The dissimilarity matrix ─────────────────────────────────
d_darl <- dist(darl_scaled)
round(as.matrix(d_darl)[1:5, 1:5], 2)The matrix is symmetric with zeros down the diagonal. Why must both of those be true? What would it mean if a diagonal entry were not zero?
▶ Run this — now the fish counts, which are a different animal entirely:
# ── Part 3.3 · Why community data is different ──────────────────────────
doubs_mat <- as.matrix(doubs_spe %>% dplyr::select(-site))
# what fraction of the table is zeros?
mean(doubs_mat == 0)
# sites 1 and 29 -- how many species does each hold?
rowSums(doubs_mat > 0)[c(1, 29)]▶ Run this — the same pair of sites, under four metrics:
# ── Part 3.4 · Four metrics, one pair of sites ──────────────────────────
e <- as.matrix(vegdist(doubs_mat, method = "euclidean"))
m <- as.matrix(vegdist(doubs_mat, method = "manhattan"))
b <- as.matrix(vegdist(doubs_mat, method = "bray"))
j <- as.matrix(vegdist(doubs_mat > 0, method = "jaccard"))
c(euclidean = e[1, 29],
manhattan = m[1, 29],
bray = b[1, 29],
jaccard = j[1, 29])Sites 1 and 29 share no species at all. Bray–Curtis and Jaccard both return exactly 1.0. What does 1.0 mean on those scales, and why can Euclidean distance never give you that information?
Site 1 holds one species; site 29 holds 26. Most of the 27 species are absent from both. Euclidean distance treats every shared absence as a point of similarity. Give one ecological reason a species might be missing from a site that has nothing to do with the site resembling any other site.
Block B · A first PCA
After lecture Parts 4–5.
4 · Correlations, then PCA
▶ Run this — does PCA have anything to collapse?
# ── Part 4.1 · Correlations ─────────────────────────────────────────────
cor_matrix <- cor(darl_numeric)
corrplot(cor_matrix, method = "color", type = "upper",
order = "hclust", addCoef.col = "black",
tl.cex = 0.7, number.cex = 0.6)Find the strongest correlation in the matrix. Does it make biological sense? Now find a pair that is close to zero — what does that pair tell you about the number of components you are likely to need?
▶ Run this:
# ── Part 4.2 · Run the PCA ──────────────────────────────────────────────
pca_result <- prcomp(darl_numeric, center = TRUE, scale. = TRUE)
summary(pca_result)▶ Run this — a scree plot:
# ── Part 4.3 · Scree plot ───────────────────────────────────────────────
fviz_eig(pca_result, addlabels = TRUE, ylim = c(0, 60))Where is the elbow? How many components would you keep, and what is your cumulative percentage at that point? Would the Kaiser rule (eigenvalue > 1) give you the same answer? Check it:
sum(pca_result$sdev^2 > 1)5 · Loadings, and the sign trap
▶ Run this:
# ── Part 5.1 · Loadings for the first three components ──────────────────
round(pca_result$rotation[, 1:3], 3)▶ Run this — PC1, as a picture:
# ── Part 5.2 · PC1 loadings ─────────────────────────────────────────────
tibble(Variable = rownames(pca_result$rotation),
Loading = pca_result$rotation[, 1]) %>%
ggplot(aes(reorder(Variable, Loading), Loading)) +
geom_col(aes(fill = Loading > 0), show.legend = FALSE) +
coord_flip() +
geom_hline(yintercept = 0, linetype = "dashed")Look at the sign of the loadings on your screen. If they are mostly negative, then a high PC1 score means a small plant. If they are mostly positive, a high score means a large one.
Either answer is correct — an eigenvector defines a line, and a line has no preferred end. What is not correct is writing “site X scores higher on PC1” and leaving it there. Translate it into biology first.
Write one sentence saying what a high PC1 score means for a cobra lily, in plain biological English, with the sign on your screen resolved. No component numbers in the sentence.
▶ Run this — the other two components:
# ── Part 5.3 · PC2 and PC3 loadings ─────────────────────────────────────
round(pca_result$rotation[, 2], 3)
round(pca_result$rotation[, 3], 3)PC2 and PC3 are each dominated by a small group of variables, not by all ten. Name the variables that dominate each, and give each component a two-word biological name.
6 · Scores, ordination, and a biplot
▶ Run this:
# ── Part 6.1 · Scores ───────────────────────────────────────────────────
pc_scores <- as_tibble(pca_result$x) %>%
mutate(site = darlingtonia$site)
pc_scores %>% dplyr::select(site, PC1, PC2, PC3) %>% head(10)▶ Run this — the ordination:
# ── Part 6.2 · Scores plot ──────────────────────────────────────────────
ggplot(pc_scores, aes(PC1, PC2, color = site)) +
geom_point(size = 3, alpha = 0.7) +
stat_ellipse(level = 0.68, linetype = 2)▶ Run this — points and arrows together:
# ── Part 6.3 · Biplot ───────────────────────────────────────────────────
fviz_pca_biplot(pca_result,
geom.ind = "point",
col.ind = darlingtonia$site,
addEllipses = TRUE,
label = "var",
repel = TRUE)Two arrows that lie almost on top of each other mean the two variables are strongly correlated; two at 90° mean they are nearly independent. Find one pair of each in your biplot. Does the correlation matrix from section 4 agree with you?
Which site’s ellipse sits furthest from the others along PC1? Using your answer to Q8, say what that means about the plants at that site — in biology, not in scores.
Block C · NMDS and testing groups
After lecture Parts 6–7.
7 · Your turn: an NMDS on the fish data
There is no code for this section. You have now run an entire ordination workflow once. Do it again on the Doubs fish, where the right metric is not Euclidean.
Work out and write the code for each of these, in your script under # ── Part 7.x:
- Run an NMDS on
doubs_matusing Bray–Curtis dissimilarity in 2 dimensions. (metaMDS()— look at?metaMDSfor thedistanceandkarguments.) - Pull out and print the stress.
- Extract the site scores with
scores(nmds, display = "sites"), join ondoubs_env$reach, and plot NMDS1 against NMDS2 coloured by reach.
If metaMDS() prints pages of iteration output, add trace = 0.
What is your stress value? Using the rules of thumb from lecture (< 0.05 excellent, < 0.1 good, < 0.2 usable, > 0.3 meaningless), how much do you trust this 2-D map?
Your NMDS plot has no percentages on its axes and the PCA plot did. Why not? What does NMDS report instead of variance explained, and what exactly is it measuring?
8 · Testing group differences
▶ Run this — route 1, ANOVA on component scores:
# ── Part 8.1 · ANOVA on PC scores ───────────────────────────────────────
anova_pc1 <- aov(PC1 ~ site, data = pc_scores)
summary(anova_pc1)
multcomp::cld(emmeans(anova_pc1, ~ site), Letters = letters)▶ Run this — route 2, PERMANOVA on the distance matrix:
# ── Part 8.2 · PERMANOVA ────────────────────────────────────────────────
adonis2(darl_scaled ~ site, data = darlingtonia, method = "euclidean")▶ Run this — and its one real assumption:
# ── Part 8.3 · Dispersion ───────────────────────────────────────────────
disp <- betadisper(vegdist(darl_scaled, method = "euclidean"),
darlingtonia$site)
anova(disp)
round(tapply(disp$distances, darlingtonia$site, mean), 2)PERMANOVA tested all ten variables at once; the ANOVA tested one component at a time. Name one thing each can tell you that the other cannot.
Is the dispersion test significant? Which site is the most variable around its own centroid? Given that result, how should you word the PERMANOVA conclusion so that it is honest?
Review and checkpoint
At this point you can:
Write-up — this part is graded
Type these answers into scripts/18_multivariate_statistics.R as # comments, in the boxes at the bottom of the skeleton.
9 · The metric decision
Write three to four sentences. You used Euclidean distance on the cobra lilies and Bray–Curtis on the fish. Justify each choice by naming the property of the data that drove it. Then say what would have gone wrong had you swapped them.
10 · The site result
Write four to five sentences reporting the Darlingtonia result as you would in a paper:
- how many components you retained and what percentage of variance they carry
- what each retained component means biologically, with the sign resolved
- the PERMANOVA result, with F, R² and p
- the dispersion check, and what it does to your confidence
11 · The figure
Produce one publication-quality figure: the Darlingtonia PCA ordination, coloured by site, with ellipses.
This one gets the full treatment — axis labels that name the component and its percentage, a legend with a real title, a caption-worthy title, and a clean theme. Save it with ggsave() as a PDF.
📤 What to turn in before next class
scripts/18_multivariate_statistics.R, running top to bottom with no errors- Your figure from section 11, as a PDF