Multivariate Ordination

Marc Los Huertos

EA030 · Compiled 2026-10-24

EA030 — Advanced Statistics — Module 19

Multivariate Ordination

Prof. Marc Los Huertos  •  Pomona College

← Back to RTricks

Concepts Data PCA NMDS PERMANOVA RDA Choosing Exercises

By the end of this module you will be able to
  1. Explain what ordination does and why standard univariate tests are inadequate for species community data.
  2. Run and interpret a PCA on continuous environmental variables, including reading a biplot.
  3. Run NMDS on species abundance data, assess ordination quality with a stress value and Shepard diagram, and overlay environmental gradients with envfit().
  4. Test whether community composition differs among groups using PERMANOVA (adonis2()).
  5. Run an RDA to constrain ordination by environmental predictors and report constrained variance.
  6. Choose among PCA, NMDS, and RDA based on data type, distribution, and research question.

What this module covers. A typical vegetation survey records 30 species at 20 sites. You cannot run 30 separate ANOVAs — you would inflate the Type I error rate, lose the correlations among species, and miss the whole point, which is to describe community structure as an integrated whole. Ordination compresses multivariate species or environmental data into a low-dimensional space where patterns of similarity and difference become visible. This module covers the three methods you will encounter most often in ecological literature and EA030 field projects.

1. Core Concepts

1.1 What ordination does

Ordination places sites (rows) in a reduced-dimension space such that sites with similar species composition end up close together, and dissimilar sites end up far apart. Species (columns) can be projected into the same space as arrows or points to show which species drive the separation.

The result is always a biplot — a single figure showing both sites and species (or environmental vectors) simultaneously.

1.2 The three methods

The three ordination methods covered in this module.
Method Full name Best input Distance metric Core question
PCA Principal Components Analysis Continuous variables (env. data) Euclidean (linear) Which env. variables co-vary?
NMDS Non-metric Multidimensional Scaling Community abundance/presence-absence Any (Bray-Curtis recommended) How do communities differ?
RDA Redundancy Analysis Species matrix + environmental predictors Euclidean on fitted values How much community variation do env. predictors explain?

1.3 Distance matrices and the double-zero problem

For species data, Euclidean distance is inappropriate. Two sites that both lack a rare species are scored as “similar” on that species — but shared absence carries no ecological information. Bray-Curtis dissimilarity ignores double-zeros and is the standard for community ecology:

\[BC_{ij} = 1 - \frac{2 \sum_k \min(y_{ik},\, y_{jk})}{\sum_k y_{ik} + \sum_k y_{jk}}\]

where \(y_{ik}\) is the abundance of species \(k\) at site \(i\).

⚠️ Do not run PCA on species abundance data PCA uses Euclidean distance, which is dominated by the many shared-absence zeros in sparse community matrices. Use NMDS or RDA for species data. PCA is appropriate for continuous, roughly normally distributed environmental variables with few zeros.

2. Data

We use two classic datasets built into the vegan package:

library(vegan)
library(ggplot2)
library(dplyr)
library(knitr)
data(dune);      data(dune.env)
data(varespec);  data(varechem)

cat("dune:     ", nrow(dune),    "sites x", ncol(dune),    "species\n")
## dune:      20 sites x 30 species
cat("dune.env: ", nrow(dune.env),"sites x", ncol(dune.env),"variables\n")
## dune.env:  20 sites x 5 variables
cat("varespec: ", nrow(varespec),"sites x", ncol(varespec),"species\n")
## varespec:  24 sites x 44 species
cat("varechem: ", nrow(varechem),"sites x", ncol(varechem),"variables\n")
## varechem:  24 sites x 14 variables
First 6 rows of dune.env. A1 = soil horizon depth (cm); Moisture, Manure = ordinal; Management, Use = factor.
A1 Moisture Management Use Manure
2.8 1 SF Haypastu 4
3.5 1 BF Haypastu 2
4.3 2 SF Haypastu 4
4.2 2 SF Haypastu 4
6.3 1 HF Hayfield 2
4.3 1 HF Haypastu 2

3. Principal Components Analysis

3.1 When to use PCA

PCA is appropriate for the continuous environmental variables in varechem — soil nitrogen, phosphorus, pH, etc. All are numeric, and the question is: which variables co-vary, and can we reduce 14 dimensions to a handful of interpretable axes?

3.2 Running PCA

In vegan, rda() with no right-hand-side formula runs an unconstrained PCA. We scale the variables (z-scores) because they are on very different units (mg/kg vs. percent vs. pH units).

# Scale to unit variance before PCA (important when units differ)
pca <- rda(varechem, scale = TRUE)
summary(pca, display = NULL)   # eigenvalues only — suppress long output
## 
## Call:
## rda(X = varechem, scale = TRUE) 
## 
## Partitioning of correlations:
##               Inertia Proportion
## Total              14          1
## Unconstrained      14          1
## 
## Eigenvalues, and their contribution to the correlations 
## 
## Importance of components:
##                          PC1    PC2    PC3     PC4     PC5     PC6     PC7
## Eigenvalue            5.1916 3.1928 1.6855 1.06899 0.81599 0.70581 0.43641
## Proportion Explained  0.3708 0.2281 0.1204 0.07636 0.05829 0.05042 0.03117
## Cumulative Proportion 0.3708 0.5989 0.7193 0.79564 0.85392 0.90434 0.93551
##                           PC8    PC9    PC10    PC11    PC12     PC13     PC14
## Eigenvalue            0.36884 0.1707 0.14953 0.08526 0.06986 0.035057 0.023615
## Proportion Explained  0.02635 0.0122 0.01068 0.00609 0.00499 0.002504 0.001687
## Cumulative Proportion 0.96185 0.9740 0.98473 0.99082 0.99581 0.998313 1.000000
# Variance explained by each axis
eig      <- eigenvals(pca)
pct_var  <- round(100 * eig / sum(eig), 1)
cum_var  <- cumsum(pct_var)
cat("PC1:", pct_var[1], "%\n")
## PC1: 37.1 %
cat("PC2:", pct_var[2], "%\n")
## PC2: 22.8 %
cat("PC3:", pct_var[3], "%\n")
## PC3: 12 %
cat("Cumulative (PC1-PC3):", cum_var[3], "%\n")
## Cumulative (PC1-PC3): 71.9 %
barplot(pct_var[1:8],
        names.arg = paste0("PC", 1:8),
        col    = "#4a7c5e", border = "white",
        ylab   = "Variance explained (%)",
        main   = "Scree Plot: Soil Chemistry PCA",
        las    = 1)
abline(h = 100 / ncol(varechem), lty = 2,
       col = "#c8832a", lwd = 1.5)
legend("topright", legend = "Equal-contribution threshold",
       lty = 2, col = "#c8832a", bty = "n", cex = 0.85)
Bar chart of variance explained by each PC axis, decreasing from left to right with an elbow after the second or third bar.

Scree plot for the soil chemistry PCA. The elbow suggests retaining two or three principal components, which together explain the majority of variance in the 14-variable dataset.

3.3 The biplot

A biplot shows sites as points and variables as arrows. Arrows pointing in the same direction indicate positively correlated variables; opposite directions indicate negative correlation. Arrow length reflects how strongly each variable loads on the displayed axes.

# Site scores and species (variable) loadings
site_sc <- scores(pca, display = "sites",   choices = 1:2)
var_sc  <- scores(pca, display = "species", choices = 1:2)

# Scale arrows for readability
arrow_scale <- 0.8 * max(abs(site_sc)) / max(abs(var_sc))
var_sc_sc   <- var_sc * arrow_scale

# Plot
par(mar = c(4, 4, 3, 1), las = 1)
plot(site_sc,
     xlab = paste0("PC1 (", pct_var[1], "%)"),
     ylab = paste0("PC2 (", pct_var[2], "%)"),
     pch  = 21, bg = "#4a7c5e", col = "#1e3a2f", cex = 1.3,
     main = "PCA Biplot: Boreal Forest Soil Chemistry")
abline(h = 0, v = 0, col = "#ede8df", lty = 1)
arrows(0, 0, var_sc_sc[, 1], var_sc_sc[, 2],
       length = 0.1, col = "#c8832a", lwd = 1.5)
text(var_sc_sc * 1.12, labels = rownames(var_sc_sc),
     cex = 0.72, col = "#7a5200", font = 2)
text(site_sc, labels = rownames(site_sc), cex = 0.6,
     pos = 3, col = "#1e3a2f")
Ordination biplot with 24 site points and 14 arrows representing soil chemistry variables, with PC1 and PC2 axes labelled with percent variance explained.

PCA biplot for boreal forest soil chemistry (varespec sites). Sites are points; environmental variables are arrows. Sites in the direction of an arrow have higher values for that variable.

📌 Reading a biplot Sites in the upper-right quadrant have high values for variables whose arrows point upper-right. The angle between two arrows approximates the correlation between those variables (small angle = high positive correlation; right angle = uncorrelated; obtuse angle = negative correlation). Project a site point perpendicularly onto an arrow to read off its approximate value on that variable.

3.4 Identifying the axes

Check which variables load most strongly on each axis:

loadings <- scores(pca, display = "species", choices = 1:2)
# Top contributors to PC1
cat("Strongest PC1 loadings:\n")
## Strongest PC1 loadings:
print(sort(abs(loadings[, 1]), decreasing = TRUE)[1:5])
##         K        Ca         P        Mg        Zn 
## 0.9368729 0.9280154 0.9180524 0.9087391 0.8926155
cat("\nStrongest PC2 loadings:\n")
## 
## Strongest PC2 loadings:
print(sort(abs(loadings[, 2]), decreasing = TRUE)[1:5])
##        Al        Fe  Humdepth        pH  Baresoil 
## 0.9719998 0.7748245 0.6824708 0.6517311 0.6386581

4. Non-metric Multidimensional Scaling

4.1 Why NMDS for community data

NMDS does not try to explain variance — it finds a low-dimensional configuration of sites that best preserves the rank order of dissimilarities. Two sites that are highly dissimilar in species composition will end up far apart; similar sites end up close together. Because it works on ranks (not raw distances), it tolerates non-normality, zeros, and skewed abundance distributions.

4.2 Stress — the quality metric

The stress value measures how faithfully the 2D plot represents the true dissimilarity structure. Lower is better:

Rule-of-thumb stress interpretation (Clarke 1993).
Stress value Ordination quality
< 0.05 Excellent
0.05 - 0.10 Good
0.10 - 0.20 Acceptable — minor distortion
> 0.20 Poor — interpret with caution

4.3 Running NMDS

set.seed(42)   # NMDS uses random starts — set seed for reproducibility
nmds <- metaMDS(dune,
                distance = "bray",    # Bray-Curtis dissimilarity
                k        = 2,         # 2 dimensions
                trymax   = 100,       # max random starts
                trace    = FALSE)

cat("Stress:", round(nmds$stress, 4), "\n")
## Stress: 0.1183
cat("Convergent solution:", nmds$converged, "\n")
## Convergent solution: 4

4.4 Shepard diagram — stress diagnostic

The Shepard diagram plots observed dissimilarities against ordination distances. A good ordination shows a monotone relationship with little scatter.

stressplot(nmds,
           main = paste0("Shepard Diagram  (stress = ",
                         round(nmds$stress, 3), ")"))
Scatter plot of observed dissimilarities on the x-axis vs ordination distances on the y-axis, with a monotone step function fit line.

Shepard diagram for the dune meadow NMDS. Points represent site pairs; the stepped line is the non-metric fit. Scatter around the fit line contributes to stress.

4.5 NMDS ordination plot

# Site and species scores
site_nm  <- scores(nmds, display = "sites")
spec_nm  <- scores(nmds, display = "species")

# Colour by Management type
mgmt_col <- c("BF" = "#1e3a2f", "HF" = "#4a7c5e",
              "NM" = "#c8832a", "SF" = "#8aad96")
pt_cols  <- mgmt_col[as.character(dune.env$Management)]

par(mar = c(4, 4, 3, 8), xpd = FALSE)
plot(site_nm,
     xlab = "NMDS1", ylab = "NMDS2",
     pch  = 21, bg = pt_cols, col = "#1e3a2f", cex = 1.6,
     main = paste0("Dune Meadow NMDS  (stress = ",
                   round(nmds$stress, 3), ")"))
abline(h = 0, v = 0, col = "#ede8df")

# Add species names (less frequent ones suppressed for readability)
spec_df <- as.data.frame(spec_nm)
spec_df$label <- rownames(spec_df)
top_sp <- spec_df[order(rowSums(spec_df[, 1:2]^2), decreasing = TRUE)[1:15], ]
text(top_sp$NMDS1, top_sp$NMDS2, labels = top_sp$label,
     cex = 0.6, col = "#7a5200", font = 3)

# Legend
par(xpd = TRUE)
legend(par("usr")[2] + 0.05, par("usr")[4],
       legend = names(mgmt_col),
       pch    = 21,
       pt.bg  = mgmt_col, col = "#1e3a2f",
       title  = "Management", bty = "n", cex = 0.85, pt.cex = 1.4)
Scatter plot of 20 sites in NMDS space with points coloured by four management categories and species names overlaid as text.

NMDS ordination of dune meadow sites coloured by Management type. Sites with similar plant communities cluster together. Management type appears to separate along the ordination axes. 

par(xpd = FALSE)

4.6 Fitting environmental vectors

envfit() projects environmental variables onto the NMDS as arrows, fitting them by correlation with the ordination axes. The \(r^2\) value and permutation p-value indicate how strongly each variable correlates with community composition.

set.seed(42)
ef <- envfit(nmds, dune.env[, c("A1", "Moisture", "Manure")],
             permutations = 999, na.rm = TRUE)
print(ef)
## 
## ***VECTORS
## 
##      NMDS1   NMDS2     r2 Pr(>r)  
## A1 0.96474 0.26321 0.3649  0.023 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
## 
## ***FACTORS:
## 
## Centroids:
##             NMDS1   NMDS2
## Moisture1 -0.5101 -0.0403
## Moisture2 -0.3938  0.0139
## Moisture4  0.2765 -0.4033
## Moisture5  0.6561  0.1476
## Manure0    0.2958  0.5790
## Manure1   -0.2482 -0.0215
## Manure2   -0.3079 -0.1866
## Manure3    0.3101 -0.2470
## Manure4   -0.3463 -0.5583
## 
## Goodness of fit:
##              r2 Pr(>r)   
## Moisture 0.5014  0.002 **
## Manure   0.4247  0.022 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
par(mar = c(4, 4, 3, 8), xpd = FALSE)
plot(site_nm,
     xlab = "NMDS1", ylab = "NMDS2",
     pch  = 21, bg = pt_cols, col = "#1e3a2f", cex = 1.6,
     main = "NMDS with Environmental Vectors")
abline(h = 0, v = 0, col = "#ede8df")
plot(ef, col = "#c8832a", lwd = 2, add = TRUE)
par(xpd = TRUE)
legend(par("usr")[2] + 0.05, par("usr")[4],
       legend = names(mgmt_col),
       pch    = 21,
       pt.bg  = mgmt_col, col = "#1e3a2f",
       title  = "Management", bty = "n", cex = 0.85, pt.cex = 1.4)
NMDS plot with site points coloured by management type and three environmental arrows labelled A1, Moisture, and Manure.

NMDS ordination of dune meadow sites with significant environmental vectors overlaid. Moisture and A1 horizon depth emerge as strong gradients. 

par(xpd = FALSE)

5. PERMANOVA

5.1 Testing group differences in community composition

NMDS shows visual separation by Management — but is it statistically significant? PERMANOVA (Permutational MANOVA, adonis2()) tests whether group centroids in multivariate space differ more than expected by chance, using permutation rather than distributional assumptions.

set.seed(42)
perm_result <- adonis2(
  dune ~ Management,
  data         = dune.env,
  method       = "bray",
  permutations = 999
)
print(perm_result)
## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = dune ~ Management, data = dune.env, permutations = 999, method = "bray")
##          Df SumOfSqs      R2      F Pr(>F)   
## Model     3   1.4686 0.34161 2.7672  0.003 **
## Residual 16   2.8304 0.65839                 
## Total    19   4.2990 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Management explains 34.2% of community variation (R2 = 0.342, p = 0.003).

📌 PERMANOVA vs ANOSIM Both test group differences in community composition. PERMANOVA (adonis2()) is preferred because it returns an R² (proportion of variance explained), handles unbalanced designs, and is less sensitive to within-group dispersion differences. ANOSIM (anosim()) only tests rank differences in among- vs within-group dissimilarities.

5.2 Checking dispersion — a PERMANOVA assumption

PERMANOVA assumes similar within-group dispersion across groups (analogous to homogeneity of variance in ANOVA). betadisper() tests this.

dist_dune <- vegdist(dune, method = "bray")
bd        <- betadisper(dist_dune, dune.env$Management)
anova(bd)   # p > 0.05 = groups have similar dispersion (assumption met)
## Analysis of Variance Table
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq F value Pr(>F)
## Groups     3 0.13831 0.046104  1.9506 0.1622
## Residuals 16 0.37816 0.023635

6. Redundancy Analysis

6.1 Constrained ordination

RDA is a constrained ordination: the axes are linear combinations of environmental predictors, not free to explain any community variation. This directly answers: how much of the community variation is explained by the measured environmental variables?

The variance decomposition is: - Constrained inertia — variation explained by the environmental model - Unconstrained inertia — residual variation not explained by predictors

6.2 Running RDA

# Hellinger transformation reduces the influence of rare species
# and makes Euclidean (RDA) distances ecologically more meaningful
dune_hel <- decostand(dune, method = "hellinger")

# RDA: species ~ environmental predictors
rda_fit <- rda(dune_hel ~ A1 + Moisture + Manure,
               data = dune.env)

summary(rda_fit, display = NULL)
## 
## Call:
## rda(formula = dune_hel ~ A1 + Moisture + Manure, data = dune.env) 
## 
## Partitioning of variance:
##               Inertia Proportion
## Total          0.5559      1.000
## Constrained    0.3513      0.632
## Unconstrained  0.2046      0.368
## 
## Eigenvalues, and their contribution to the variance 
## 
## Importance of components:
##                         RDA1    RDA2    RDA3    RDA4    RDA5    RDA6     RDA7
## Eigenvalue            0.1476 0.09674 0.03638 0.02452 0.01935 0.01222 0.008914
## Proportion Explained  0.2655 0.17402 0.06544 0.04410 0.03481 0.02198 0.016035
## Cumulative Proportion 0.2655 0.43951 0.50494 0.54904 0.58385 0.60583 0.621866
##                           RDA8     PC1     PC2     PC3    PC4     PC5     PC6
## Eigenvalue            0.005648 0.05115 0.03971 0.03461 0.0189 0.01792 0.01219
## Proportion Explained  0.010161 0.09202 0.07143 0.06225 0.0340 0.03224 0.02193
## Cumulative Proportion 0.632027 0.72405 0.79548 0.85773 0.8917 0.92397 0.94591
##                           PC7      PC8      PC9     PC10     PC11
## Eigenvalue            0.01067 0.008518 0.006583 0.003232 0.001066
## Proportion Explained  0.01920 0.015323 0.011842 0.005814 0.001917
## Cumulative Proportion 0.96510 0.980427 0.992269 0.998083 1.000000
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                         RDA1    RDA2    RDA3    RDA4    RDA5    RDA6     RDA7
## Eigenvalue            0.1476 0.09674 0.03638 0.02452 0.01935 0.01222 0.008914
## Proportion Explained  0.4201 0.27533 0.10354 0.06978 0.05508 0.03477 0.025371
## Cumulative Proportion 0.4201 0.69539 0.79893 0.86870 0.92378 0.95855 0.983924
##                           RDA8
## Eigenvalue            0.005648
## Proportion Explained  0.016076
## Cumulative Proportion 1.000000
# Proportion of variance constrained by the environmental model
inertia <- summary(rda_fit, display = NULL)$cont$importance
cat("Total inertia:", round(rda_fit$tot.chi, 4), "\n")
## Total inertia: 0.5559
cat("Constrained:  ",
    round(rda_fit$CCA$tot.chi, 4),
    sprintf("(%.1f%%)\n", 100 * rda_fit$CCA$tot.chi / rda_fit$tot.chi))
## Constrained:   0.3513 (63.2%)
cat("Unconstrained:",
    round(rda_fit$CA$tot.chi, 4),
    sprintf("(%.1f%%)\n", 100 * rda_fit$CA$tot.chi / rda_fit$tot.chi))
## Unconstrained: 0.2046 (36.8%)

6.3 Testing significance of the RDA model

# Global permutation test for the full model
set.seed(42)
anova(rda_fit, permutations = 999)
## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = dune_hel ~ A1 + Moisture + Manure, data = dune.env)
##          Df Variance      F Pr(>F)    
## Model     8  0.35135 2.3617  0.001 ***
## Residual 11  0.20456                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Test each RDA axis individually
set.seed(42)
anova(rda_fit, by = "axis", permutations = 999)
## Permutation test for rda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = dune_hel ~ A1 + Moisture + Manure, data = dune.env)
##          Df Variance      F Pr(>F)    
## RDA1      1 0.147585 7.9363  0.001 ***
## RDA2      1 0.096737 5.6750  0.001 ***
## RDA3      1 0.036377 2.3118  0.305    
## RDA4      1 0.024516 1.6779  0.709    
## RDA5      1 0.019351 1.4190  0.793    
## RDA6      1 0.012217 0.9556  0.958    
## RDA7      1 0.008914 0.7408  0.958    
## RDA8      1 0.005648 0.4970  0.958    
## Residual 11 0.204557                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

6.4 RDA biplot

rda_sites <- scores(rda_fit, display = "sites",    choices = 1:2)
rda_spec  <- scores(rda_fit, display = "species",  choices = 1:2)
rda_env   <- scores(rda_fit, display = "bp",       choices = 1:2)

# Axis labels: proportion of total inertia each RDA axis explains
ax_pct <- round(100 * eigenvals(rda_fit)[1:2] / rda_fit$tot.chi, 1)

# Scale species and env arrows
sp_scale  <- 0.6 * max(abs(rda_sites)) / max(abs(rda_spec))
env_scale <- 0.8 * max(abs(rda_sites)) / max(abs(rda_env))

par(mar = c(4, 4, 3, 8), xpd = FALSE)
plot(rda_sites,
     xlab = paste0("RDA1 (", ax_pct[1], "% of total)"),
     ylab = paste0("RDA2 (", ax_pct[2], "% of total)"),
     pch  = 21, bg = pt_cols, col = "#1e3a2f", cex = 1.5,
     main = "RDA: Dune Vegetation ~ A1 + Moisture + Manure")
abline(h = 0, v = 0, col = "#ede8df")

# Species scores
text(rda_spec * sp_scale,
     labels = rownames(rda_spec),
     cex = 0.55, col = "#2d5a3d", font = 3)

# Environmental arrows
arrows(0, 0, rda_env[, 1] * env_scale, rda_env[, 2] * env_scale,
       length = 0.12, col = "#c8832a", lwd = 2)
text(rda_env * env_scale * 1.15, labels = rownames(rda_env),
     cex = 0.8, col = "#7a5200", font = 2)

par(xpd = TRUE)
legend(par("usr")[2] + 0.05, par("usr")[4],
       legend = names(mgmt_col), pch = 21,
       pt.bg  = mgmt_col, col = "#1e3a2f",
       title  = "Management", bty = "n", cex = 0.85, pt.cex = 1.4)
RDA biplot with 20 coloured site points, species name text, and three environmental arrows labelled A1, Moisture, and Manure, with RDA axis percentages on the axes.

RDA biplot of dune meadow vegetation constrained by soil horizon depth (A1), moisture, and manure input. Site points are coloured by management type. Environmental arrows show the direction of each predictor’s effect on community composition. 

par(xpd = FALSE)

💡 Hellinger transformation for RDA Raw abundance data has many zeros and is dominated by common species. The Hellinger transformation (decostand(x, “hellinger”)) divides each value by the site total, then takes the square root. This gives appropriate weight to both common and rare species and makes Euclidean distance (used internally by RDA) ecologically meaningful.

7. Choosing a Method

Decision guide for multivariate ordination.
Scenario Recommended approach
Explore co-variation among continuous env. variables PCA (scale = TRUE if variables are in different units)
Visualise community similarity among sites, no env. data NMDS with Bray-Curtis
Test whether groups differ in community composition PERMANOVA via adonis2()
Relate community composition to env. predictors RDA (after Hellinger transformation)
Data with many zeros and rare species NMDS or RDA with Hellinger — avoid raw PCA
Presence-absence data NMDS with Jaccard distance; or RDA with decostand(x, ‘pa’)

🔑 Key Takeaways

✏️ Practice Exercises
Exercise 1 — PCA on varechem

The varechem dataset has 14 soil chemistry variables at 24 boreal forest sites. (a) Run PCA with rda(varechem, scale = TRUE). (b) How many axes are needed to explain 70% of variance? (c) Make a biplot. Which soil variables are most strongly positively correlated with each other on PC1? Which point in opposite directions (negatively correlated)?

Show solution 
library(vegan)
data(varechem)
pca_vc <- rda(varechem, scale = TRUE)

# Variance explained
eig <- eigenvals(pca_vc)
cumsum(100 * eig / sum(eig))   # find where cumulative exceeds 70%

# Biplot
biplot(pca_vc, scaling = 2,
       main = "varechem PCA biplot",
       col  = c("#1e3a2f", "#c8832a"))

Typically 3-4 axes are needed for 70%. On PC1, N and Al often load positively together (nutrient-poor acidic soils), while Ca, Mg, and pH load negatively (richer, less acidic soils) — consistent with a soil fertility gradient.

Exercise 2 — NMDS on varespec

Run NMDS on the varespec lichen data with Bray-Curtis distance. (a) Report the stress value and assess quality. (b) Fit the varechem variables onto the ordination with envfit() (permutations = 999, set.seed(42)). (c) Which soil variables are most significantly correlated with lichen community composition? Plot the ordination with significant vectors overlaid.

Show solution 
library(vegan)
data(varespec); data(varechem)
set.seed(42)
nmds_v <- metaMDS(varespec, distance = "bray", k = 2,
                  trymax = 100, trace = FALSE)
cat("Stress:", nmds_v$stress, "\n")

set.seed(42)
ef_v <- envfit(nmds_v, varechem, permutations = 999)
print(ef_v)

plot(nmds_v, display = "sites",
     main = "varespec NMDS with soil chemistry vectors")
plot(ef_v, p.max = 0.05, col = "#c8832a", add = TRUE)

Stress is typically around 0.18 — acceptable. Al, Fe, and humus thickness often emerge as the strongest correlates of lichen community composition, reflecting the importance of soil acidity and organic matter in boreal forest ground-layer communities.

Exercise 3 — PERMANOVA on dune data

Test whether dune meadow communities differ by Use type (Hayfield, Haypastu, Pasture) as well as by Management. (a) Run adonis2(dune ~ Use, data = dune.env, method = “bray”). (b) Check dispersion with betadisper(). (c) Run a model with both predictors: dune ~ Management + Use. Which explains more variance?

Show solution 
library(vegan)
data(dune); data(dune.env)
set.seed(42)

# Use alone
adonis2(dune ~ Use, data = dune.env, method = "bray", permutations = 999)

# Dispersion check
dist_d <- vegdist(dune, method = "bray")
anova(betadisper(dist_d, dune.env$Use))

# Combined model
set.seed(42)
adonis2(dune ~ Management + Use, data = dune.env,
        method = "bray", permutations = 999, by = "margin")

Management typically explains more variance than Use. The combined model shows each predictor’s marginal contribution; compare R2 values to rank importance. Check that betadisper p > 0.05 before trusting PERMANOVA.

Exercise 4 — RDA on varespec and varechem

Run an RDA of varespec (Hellinger-transformed) constrained by Al, N, and pH from varechem. (a) What percentage of community variation is explained by these three soil variables? (b) Test the model significance with anova(). (c) Make a biplot and describe the ecological gradient along RDA1.

Show solution 
library(vegan)
data(varespec); data(varechem)

vs_hel  <- decostand(varespec, "hellinger")
rda_vs  <- rda(vs_hel ~ Al + N + pH, data = varechem)

# Variance explained
cat("Constrained:", round(100 * rda_vs$CCA$tot.chi / rda_vs$tot.chi, 1), "%\n")

# Significance
set.seed(42)
anova(rda_vs, permutations = 999)

# Biplot
plot(rda_vs, scaling = 2,
     main = "RDA: Lichen ~ Al + N + pH",
     col  = c("#1e3a2f", "#c8832a"))

Al and N typically explain 20-35% of lichen community variation. RDA1 often separates acidic, aluminium-rich, nutrient-poor sites (with ericaceous and lichen species) from sites with higher pH and N (supporting more productive mosses and grasses).