Activity: NMDS & PERMANOVA

NMDS & PERMANOVA

Hands-on activity: NMDS ordination, PERMANOVA, ANOSIM, and environmental fitting, using the iris dataset treated as a community-style abundance matrix.

Author

Bill Perry

Worksheet: NMDS & PERMANOVA

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 iris dataset, treating the four flower measurements as if they were abundances of different “species” at different “sites” — a handy stand-in for real community data.

▶ Run this in your Script:

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

theme_set(theme_minimal(base_size = 9))

What is NMDS? NMDS (Non-metric Multidimensional Scaling) is an ordination technique that:

Visualizes dissimilarity between objects in reduced dimensions
Preserves rank order of distances, not exact distances
Works well with non-linear ecological relationships
Makes few assumptions about data structure

When to use NMDS: community data (species abundance or presence/absence matrices), non-linear relationships (when PCA assumptions are violated), or complex ecological gradients (multiple environmental factors affecting communities).

Key concepts of NMDS:

Dissimilarity matrices instead of covariance
Stress values measure goodness of fit (< 0.2 is acceptable)
Iterative algorithm to find optimal configuration
No eigenvalues - axes have no inherent meaning
Rank-based - preserves order, not exact distances

Important

Critical first step: Always check your stress value! Stress < 0.1 is excellent, 0.1-0.2 is good, > 0.2 is poor representation.

Part 1 · Data Preparation and Exploration

We’ll use the iris dataset for this analysis, treating it as if the measurements represent abundances of different “species” at different “sites”.

▶ Run this:

iris_df <- read_csv("data/iris.csv") %>%
  clean_names()

head(iris_df)

▶ Run this — separate the grouping variable from the numeric “community” matrix:

iris_species_df <- iris_df %>%
  dplyr::select(species)

iris_numeric_df <- iris_df %>%
  dplyr::select(-species)

str(iris_numeric_df)

▶ Run this — a quick pairs plot to see relationships:

iris_pairs_plot <- iris_df %>%
  dplyr::select(-species) %>%
  pairs(main = "Iris Measurements Relationships")

Part 2 · Running NMDS

Step 1: Calculate a distance matrix

# For iris data, we'll use Euclidean distance since these are measurements
# (Bray-Curtis is more common for real abundance/count data)
iris_dist <- dist(iris_numeric_df, method = "euclidean")

iris_dist[1:5]

Step 2: Run NMDS

set.seed(123)
iris_nmds_model <- metaMDS(iris_numeric_df,
                           distance = "euclidean",
                           k = 2,
                           trymax = 100)

iris_nmds_model

✏️ Your turn: Fill in from your own output — Stress value: ______. Is that excellent (<0.1), good (0.1–0.2), or poor (>0.2)? ________________________

Step 3: Extract NMDS scores

nmds_scores_df <- as.data.frame(iris_nmds_model$points) %>%
  rename(nmds1 = MDS1, nmds2 = MDS2) %>%
  bind_cols(iris_species_df)

head(nmds_scores_df)

Step 4: Create an NMDS plot

nmds_basic_plot <- ggplot(nmds_scores_df, aes(x = nmds1, y = nmds2, color = species)) +
  geom_point(size = 3, alpha = 0.8) +
  labs(title = "NMDS Ordination of Iris Data",
       subtitle = paste("Stress =", round(iris_nmds_model$stress, 3)),
       x = "NMDS1",
       y = "NMDS2",
       color = "Species") +
  theme_minimal()

nmds_basic_plot

Step 5: Add confidence ellipses

nmds_ellipse_plot <- ggplot(nmds_scores_df, aes(x = nmds1, y = nmds2, color = species)) +
  geom_point(size = 3, alpha = 0.8) +
  stat_ellipse(level = 0.95, linewidth = 1) +
  labs(title = "NMDS Ordination with 95% Confidence Ellipses",
       subtitle = paste("Stress =", round(iris_nmds_model$stress, 3)),
       x = "NMDS1",
       y = "NMDS2",
       color = "Species") +
  theme_minimal()

nmds_ellipse_plot

Step 6: Stress plot (Shepard diagram)

stressplot(iris_nmds_model, main = "Shepard Diagram: Ordination vs Original Distances")

✏️ Your turn: Do the three species separate cleanly on the NMDS plot, or do any overlap? ________________________

Part 3 · PERMANOVA Analysis

PERMANOVA (Permutational Multivariate Analysis of Variance) tests whether groups have different multivariate centroids using permutation tests.

Step 1: Run PERMANOVA

set.seed(456)
iris_permanova_model <- adonis2(iris_numeric_df ~ species,
                                data = iris_df,
                                method = "euclidean",
                                permutations = 999)

iris_permanova_model

✏️ Your turn: Fill in from your own output — F-statistic: ______, R² (variance explained): ______, p-value: ______. Do species differ significantly in multivariate space?

Step 2: Check homogeneity of dispersions

Before interpreting PERMANOVA, we need to check if groups have similar multivariate spread.

iris_dist_full <- dist(iris_numeric_df)
dispersion_model <- betadisper(iris_dist_full, iris_df$species)

dispersion_test <- anova(dispersion_model)
dispersion_test

Step 3: Visualize dispersions

plot(dispersion_model, main = "Multivariate Dispersion by Species")

Step 4: Pairwise PERMANOVA

pairwise_permanova <- function(data_matrix, groups, distance_method = "euclidean") {
  group_levels <- unique(groups)
  comparisons <- combn(group_levels, 2)

  results_df <- data.frame(
    group1 = character(),
    group2 = character(),
    f_statistic = numeric(),
    r_squared = numeric(),
    p_value = numeric()
  )

  for(i in 1:ncol(comparisons)) {
    group1 <- comparisons[1, i]
    group2 <- comparisons[2, i]
    subset_indices <- which(groups %in% c(group1, group2))

    subset_data <- data_matrix[subset_indices, ]
    subset_groups <- groups[subset_indices]

    temp_result <- adonis2(subset_data ~ subset_groups,
                          method = distance_method,
                          permutations = 999)

    results_df <- rbind(results_df, data.frame(
      group1 = group1,
      group2 = group2,
      f_statistic = temp_result$F[1],
      r_squared = temp_result$R2[1],
      p_value = temp_result$"Pr(>F)"[1]
    ))
  }

  results_df$p_adjusted <- p.adjust(results_df$p_value, method = "bonferroni")

  return(results_df)
}

pairwise_results_df <- pairwise_permanova(iris_numeric_df, iris_df$species)
pairwise_results_df

✏️ Your turn: Which pair of species is the most different (largest F-statistic)? Which pair is the least different? ________________________

Tip

🚀 If you finish early: Does the dispersion test in Step 2 suggest the three species have similar multivariate spread? If not, how would that affect how you interpret the PERMANOVA result from Step 1?

# Write your notes/code here:

Part 4 · ANOSIM Analysis

ANOSIM (Analysis of Similarities) tests whether there is a significant difference between groups using rank dissimilarities.

Step 1: Run ANOSIM

set.seed(789)
iris_anosim_model <- anosim(iris_dist_full, iris_df$species, permutations = 999)

iris_anosim_model

✏️ Your turn: Fill in from your own output — R statistic: ______, p-value: ______. R close to 1 indicates strong separation between groups; how strong is the separation here?

Step 2: Plot ANOSIM results

plot(iris_anosim_model, main = "ANOSIM Results: Distribution of Permuted R Statistics")

Part 5 · Environmental Fitting (Optional)

If we had environmental variables, we could fit them to the ordination. For demonstration, let’s use petal_length as a stand-in “environmental” variable.

▶ Run this:

env_data_df <- data.frame(petal_length = iris_df$petal_length)

env_fit_model <- envfit(iris_nmds_model, env_data_df, permutations = 999)
env_fit_model

▶ Run this — visualize the environmental vector on the ordination:

env_coords_df <- as.data.frame(env_fit_model$vectors$arrows * 2)  # Scale for visibility
env_coords_df$variable <- rownames(env_coords_df)

nmds_env_plot <- ggplot(nmds_scores_df, aes(x = nmds1, y = nmds2, color = species)) +
  geom_point(size = 3, alpha = 0.8) +
  stat_ellipse(level = 0.95, linewidth = 1, alpha = 0.3) +
  geom_segment(data = env_coords_df,
               aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
               arrow = arrow(length = unit(0.3, "cm")),
               color = "black", linewidth = 1) +
  geom_text(data = env_coords_df,
            aes(x = NMDS1 * 1.1, y = NMDS2 * 1.1, label = variable),
            color = "black", size = 4) +
  labs(title = "NMDS with Environmental Vector",
       subtitle = paste("Stress =", round(iris_nmds_model$stress, 3)),
       x = "NMDS1",
       y = "NMDS2",
       color = "Species") +
  theme_minimal()

nmds_env_plot

Tip

🚀 If you finish early: Repeat the environmental fitting with sepal_width instead of petal_length. Does it point in a similar direction, or a very different one?

# Write your code here:

Summary checklist for NMDS and PERMANOVA

Tip

Analysis checklist

Prepare your data - ensure numeric matrix format.
Choose appropriate distance measure - Bray-Curtis for abundance data, Euclidean for measurement data.
Run NMDS with sufficient iterations.
Check stress value - must be < 0.2.
Create ordination plots with groups identified.
Test homogeneity of dispersions before PERMANOVA.
Run PERMANOVA to test group differences.
Consider ANOSIM as complementary test.
Fit environmental variables if available.

Key points to remember:

NMDS preserves rank order of distances, not exact values.
Stress < 0.2 is acceptable, < 0.1 is excellent.
PERMANOVA tests centroids, ANOSIM tests overlap.
Check dispersion homogeneity - violated assumption affects interpretation.
Multiple comparisons require p-value adjustment.
Axes have no inherent meaning in NMDS (unlike PCA).
Use appropriate distance measures for your data type.

Important

Key takeaways from NMDS/PERMANOVA analysis

NMDS is flexible - works with any distance measure and makes few assumptions.
Stress indicates fit quality - always report and check this value.
PERMANOVA is powerful but assumes homogeneous dispersions.
ANOSIM is complementary - provides different perspective on group separation.
Visualization is crucial - always plot your ordination results.
Environmental fitting helps interpret ecological patterns.
Permutation tests avoid distributional assumptions.

Remember: NMDS is iterative and may find different solutions - always set a seed for reproducibility!

Review and checkpoint

At this point you can:

Explain when NMDS is preferred over PCA
Run metaMDS() and interpret the stress value
Build an NMDS ordination plot with confidence ellipses
Run and interpret PERMANOVA, including checking dispersion homogeneity with betadisper()
Run and interpret ANOSIM as a complementary test
Fit an environmental (or other continuous) vector onto an ordination with envfit()

Note

📤 What to turn in before next class

Upload both of these to the course management system:

Your code — the scripts/ folder (or just 10_nmds_permanova.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: NMDS & PERMANOVA

How to use this worksheet

Part 1 · Data Preparation and Exploration

Part 2 · Running NMDS

Step 1: Calculate a distance matrix

Step 2: Run NMDS

Step 3: Extract NMDS scores

Step 4: Create an NMDS plot

Step 5: Add confidence ellipses

Step 6: Stress plot (Shepard diagram)

Part 3 · PERMANOVA Analysis

Step 1: Run PERMANOVA

Step 2: Check homogeneity of dispersions

Step 3: Visualize dispersions

Step 4: Pairwise PERMANOVA

Part 4 · ANOSIM Analysis

Step 1: Run ANOSIM

Step 2: Plot ANOSIM results

Part 5 · Environmental Fitting (Optional)

Summary checklist for NMDS and PERMANOVA

Review and checkpoint

Getting unstuck