Multivariate Data Analysis Module

1) Summarizing & Exploring Your Data

Before conducting a multivariate analysis, first examine the structure and quality of the data, summarize sampling effort, and visualize the main abundance patterns. This helps identify missing observations, unequal sampling effort, dominant or uncommon taxa, and other features that may influence later analytical decisions.

(Download the tutorial script, data and look-up files to the same folder before running the tutorial)

This script demonstrates how to:

  • Reshape community data between wide and long formats.
  • Check data quality and sampling effort.
  • Consolidate taxonomic categories.
  • Calculate abundance summaries across different sampling units.
  • Create stacked abundance bar plots and proportional-composition plots.
  • Calculate univariate diversity metrics with the vegan package

Additional multivariate data prep and summarization resources

2) nMDS Ordination Tutorial & Info

Nonmetric multidimensional scaling, or nMDS, creates a two-dimensional map of differences in community composition. Points plotted close together have more similar assemblages, while those farther apart have more different assemblages.

The tutorial also introduces PERMANOVA (at the end) as a companion analysis for testing differences in multivariate composition.

Before class

Watch these two videos before beginning the tutorial:

  1. Non-Metric Multidimensional Scaling (NMDS) overview

    Focus on how nMDS works, what stress represents, and why the orientation of the axes is arbitrary (not plots are made with base R plotting, not what you’ll be using, but good conceptual overview).

  2. nMDS in R from the Riffomonas channel

    Provides examples of plotting nMDS using ggplot2 and vegan.

Riffomonas is also a useful channel for tidyverse and ggplot2 examples from a microbial ecologist. Microbial and molecular datasets have additional characteristics and challenges, so not every choice will transfer directly to other ecological community data.

nMDS Tutorial files

Place the R script and Excel workbook in the same folder before running the tutorial.

Additional resources

(don’t need to review before class)

What the tutorial covers

The guided analysis uses the invertebrate worksheet from the coral reef dataset and follows this sequence:

  1. Inspect community data and sampling effort.
  2. Average replicate transects within each site and year.
  3. Compare square-root and fourth-root transformations of abundance data.
  4. Calculate taxon prevalence and optionally remove uncommon taxa.
  5. Construct the sample by taxon community (“comm”) table required by vegan.
  6. Calculate Bray-Curtis dissimilarities and a reproducible nMDS with vegan::metaMDS().
  7. Create publication-style plots by assigning color and shape to various variables, add point labels.
  8. Relate assemblage patterns to selected taxa, fitted vectors, species scores, and a stacked abundance plot.
  9. Run introductory PERMANOVA models with vegan::adonis2().

Main analytical choices

Sampling unit

Before running an ordination, decide what each point will represent. A point could represent an individual transect or plot, or replicate samples could first be averaged to a higher level, such as a site, site-year, habitat, or treatment.

This decision should match the research question and study design. Retaining individual replicates shows within-group variation, while averaging focuses the analysis on differences among the higher-level sampling units. Averaging replicates also prevents groups with greater sampling effort from receiving more weight simply because they contain more samples. However, it reduces the number of observations and removes information about variation among replicates. Clearly define and report the sampling unit used in the ordination.

Transformation

The transformation should be selected based on the biological question/objective rather than on which option produces the clearest separation in the ordination.

Before running an nMDS, abundance data are typically transformed to reduce the influence of the most abundant taxa. The choice is commonly between a square-root and fourth-root transformation (4th root typically yields similar result as a log transformation). A square-root transformation retains more information about differences in relative abundance and which taxa dominate, while a fourth-root transformation gives uncommon taxa greater influence and may be preferred when the focus is overall assemblage diversity.

Untransformed abundance and presence/absence represent the two extremes. Untransformed data give dominant taxa the greatest influence, while presence/absence removes all information about abundance.

Choice Relative influence of abundant and uncommon taxa Typical use
No transformation Dominant taxa have the greatest influence When absolute or relative abundance differences among dominant taxa are the primary focus
Square root Moderately reduces the influence of dominant taxa When relative abundance and which taxa dominate remain important
Fourth root More strongly reduces the influence of dominant taxa When uncommon taxa and assemblage-wide biodiversity should receive greater influence
log(x + 1) Strongly reduces dominance and is often similar to a fourth-root transformation An alternative when abundance values span a wide range
Presence/absence Removes all abundance information When only taxon occurrence and assemblage membership are relevant

Rare taxa

Rare taxa are sometimes removed before calculating community dissimilarities because taxa occurring in very few samples may reflect chance observations, inconsistent detection, or differences in taxonomic identification. They may also contribute little information about broader assemblage patterns.

Filtering can be based on prevalence, such as the percentage or number of samples in which a taxon occurs, or on minimum total abundance. The appropriate rule depends on the dataset and research objectives. Rare taxa should be retained when they are biologically important or central to the research question. Any filtering rule should be reported, and conclusions should be checked across a small range of defensible thresholds.

Interpreting an nMDS

Interpret the relative distances among points, not the direction or numerical scale of the axes. Points plotted close together represent samples with more similar assemblages, while points farther apart represent more different assemblages. An nMDS configuration can be rotated, reflected, or shifted without changing these relationships.

Colors, shapes, and labels can identify relevant sampling groups or continuous gradients. Centroids and ellipses can summarize group locations and within-group variation, while selected-taxon abundances, fitted vectors, or species scores can help identify taxa associated with assemblage patterns. These additions support interpretation but should not obscure the underlying distances among samples.

Stress

Community composition is multidimensional because each taxon represents a separate dimension. A two-dimensional nMDS compresses these relationships into a flat plot. Stress measures the resulting distortion, with lower values indicating a better representation.

These are general guidelines, not strict thresholds:

  • 0 to 0.10: Good

  • 0.10 to 0.20: Acceptable

  • 0.20 to 0.30: Potentially acceptable for complex ecological data, but interpret cautiously

  • Greater than 0.30: Poor

Values near zero should also be checked because they can result from too few samples, duplicate samples, or limited variation. Interpret stress alongside convergence and the Shepard plot.

PERMANOVA (using adonis2 from vegan)

Permutational multivariate analysis of variance, or PERMANOVA, tests whether community composition differs among groups or along continuous gradients. It uses a community dissimilarity matrix and permutations rather than assuming multivariate normality. In R, PERMANOVA can be implemented with vegan::adonis2().

Models can include one or multiple predictors. For models with multiple predictors, by = “margin” tests each term after accounting for the others. The sampling unit and permutation structure must match the study design, including any blocking, nesting, or repeated sampling. Because PERMANOVA can also be influenced by differences in within-group dispersion, group dispersions should be examined when interpreting the results.

What to report in a Methods section

A Methods description should identify the sampling unit and averaging procedure; taxon consolidation and filtering; transformation (square or 4th root); Bray-Curtis dissimilarity; final stress (typically on the plot); (if applicable) species-vector display criteria; PERMANOVA terms and number of permutations;

Here’s example methods text from my grad students paper Calderon et. al 2024:

To quantify the overall diet composition of Garibaldi, the weight for each diet content category was averaged across specimens within a given location, reef type, and sex. A similarity matrix was constructed with the Bray-Curtis similarity coefficient using square-root transformed weight values. These relationships were visualized with a two-dimensional non-metric multidimensional scaling plot (nMDS) using the metaMDS function in the ‘vegan’ package (Oksanen et al., 2022). Ordination scores were then visualized as vectors on the nMDS for individual stomach content categories that significantly (p < 0.05) correlated with the diet composition patterns in the nMDS, with arrow length indicating the relative strength of each correlation. The effect of reef type, temperature, and sex on the observed variation in diet composition was assessed using a permutational multivariate analysis of variance (PERMANOVA) with the adonis2 function in the ‘vegan’ package (Oksanen et al., 2022). Because our dataset was unbalanced (i.e., artificial reef types were only present at some locations), the by = “margin” argument was included in the function to analyze the marginal effects of each variable in the model after accounting for all other variables.

Practice adapting the nMDS analysis code to a new data set

Use the terrestrial oribatid mite dataset to practice adapting the tutorial workflow to new community and environmental data. The code in the script provides and prepares the mite data, and gives guidance for exploring categorical and continuous environmental predictors.