Bio Ecological Genomics Community Ecology

Version Compatibility

Reference examples tested with: vegan 2.6+, indicspecies 1.7+, Hmsc 3.0+, sjSDM 1.0+, ade4 1.7+, picante 1.8+, ggplot2 3.5+

Before using code patterns, verify installed versions match. If versions differ:

• R: packageVersion('<pkg>') then ?function_name to verify parameters

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Community Ecology

"Test how environmental gradients structure my species communities" -> Constrained ordination (CCA / RDA / db-RDA) with explicit dispersion testing alongside PERMANOVA, indicator-species analysis with group-size correction, Joint Species Distribution Models for residual covariance and prediction, and trait-environment testing with statistically corrected permutation schemes.

• R: vegan::rda(), vegan::dbrda(), vegan::adonis2() for ordination and PERMANOVA
• R: vegan::betadisper() for the mandatory PERMDISP companion to PERMANOVA
• R: Hmsc or sjSDM for joint species distribution modeling
• R: indicspecies::multipatt(..., func = 'IndVal.g') for indicator species

The Single Most Important Modern Insight -- Always run PERMDISP alongside PERMANOVA

Anderson & Walsh 2013 Ecol Monogr 83(4):557-574 established that PERMANOVA's pseudo-F is sensitive to dispersion heterogeneity — a significant PERMANOVA can reflect centroid difference, dispersion difference, or both. Running PERMANOVA without PERMDISP is the single most common methodological failure in community-ecology papers. If betadisper is significant alongside a significant PERMANOVA, the location-difference conclusion is not supported by the data alone; the two could be entirely a dispersion artifact.

A second cornerstone insight: residual species-species covariance in Joint Species Distribution Models is NOT a clean estimator of biotic interaction (Pollock 2014, reaffirmed by Zurell 2018 and Poggiato 2021). It reflects unmeasured covariates, dispersal limitation, sampling artifacts, AND any genuine biotic interactions, in unknown proportions. Skills that interpret residual covariance as interaction are over-interpreting.

Algorithmic Taxonomy

Decision Tree by Scenario

CCA vs RDA — Gradient-Length Decision

Goal: Choose between unimodal (CCA) and linear (RDA) constrained ordination based on the dominant gradient length in the species data.

Approach: Run a Detrended Correspondence Analysis (DCA) on the raw community matrix; use the axis-1 SD length as the gradient-length metric. > 3 SD suggests unimodal CCA; < 3 SD suggests linear RDA on Hellinger-transformed data; 2-3 SD is a gray zone where either is defensible.

library(vegan)

# Step 1: Check gradient length
dca <- decorana(species_matrix)
dca  # axis 1 length in SD units

# Step 2a: Long gradient -> CCA
cca_result <- cca(species_matrix ~ temperature + precipitation + pH + elevation,
                  data = env_data)
anova(cca_result, by = 'margin', permutations = 999)

# Step 2b: Short gradient -> RDA with Hellinger transformation
# Hellinger (Legendre & Gallagher 2001) is MANDATORY before RDA on community data
species_hell <- decostand(species_matrix, method = 'hellinger')
rda_result <- rda(species_hell ~ temperature + precipitation + pH + elevation,
                  data = env_data)
RsquareAdj(rda_result)$adj.r.squared
anova(rda_result, by = 'margin', permutations = 999)

# Forward selection with adjusted R-squared criterion (Peres-Neto 2006 Ecology 87:2614)
rda_null <- rda(species_hell ~ 1, data = env_data)
rda_full <- rda(species_hell ~ ., data = env_data)
rda_sel <- ordiR2step(rda_null, scope = formula(rda_full),
                      direction = 'forward', permutations = 999)

# VIF check: > 10 indicates problematic multicollinearity
vif.cca(rda_sel)

PERMANOVA with the Mandatory PERMDISP Companion

Goal: Test whether community composition differs across groups while detecting the dispersion-heterogeneity confound.

Approach: Run adonis2 (modern PERMANOVA per Anderson 2001 Austral Ecol) on Hellinger or Bray-Curtis distances; THEN run betadisper (PERMDISP per Anderson 2006 Biometrics) to test whether group dispersions are unequal. Report BOTH results. If betadisper is significant, the adonis2 conclusion of centroid difference is not supported — the apparent "group difference" may be entirely a dispersion artifact.

library(vegan)

# Distance matrix
bray_dist <- vegdist(species_matrix, method = 'bray')

# PERMANOVA via adonis2 (modern API; adonis() is deprecated)
# by='margin' for unbalanced designs (sequential SS gives wrong result)
permanova <- adonis2(bray_dist ~ habitat + soil_pH, data = env_data,
                     by = 'margin', permutations = 999)
permanova

# MANDATORY companion: PERMDISP via betadisper
# Tests homogeneity of multivariate dispersions across groups
disp <- betadisper(bray_dist, env_data$habitat)
disp_test <- permutest(disp, permutations = 999)
disp_test

# Interpretation rule:
# PERMANOVA p < 0.05 AND betadisper p > 0.05 -> location difference is real
# PERMANOVA p < 0.05 AND betadisper p < 0.05 -> CONFOUNDED, cannot conclude location difference
# PERMANOVA p > 0.05 -> no group difference detected

# Visualize dispersions
plot(disp)
boxplot(disp)

# For pairwise group comparisons, do NOT use Bonferroni-corrected pairwise PERMANOVA
# (permutation tests with overlapping sets do not give nominal FWER from Bonferroni);
# use pairwiseAdonis::pairwise.adonis2 with FDR correction instead
# install.packages('pairwiseAdonis')

Joint Species Distribution Models — HMSC vs sjSDM

Goal: Model species occurrences jointly to capture environmental responses, traits, phylogeny, and residual covariance among species.

Approach: For S < 500 species with rich theory and traits/phylogeny: use HMSC (Ovaskainen 2017 Ecol Lett 20:561-576; current R package Tikhonov 2020 Methods Ecol Evol 11:442-447) for Bayesian inference with explicit ecological priors. For S > 500 species (modern metabarcoding datasets): use sjSDM (Pichler & Hartig 2021 Methods Ecol Evol 12:2159-2173) which is orders of magnitude faster via Monte Carlo approximation of the joint likelihood with elastic-net regularization. The Wilkinson 2019 Methods Ecol Evol 10:198-211 benchmark is essential reading before picking among HMSC, sjSDM, gjam, and BayesComm. Do NOT interpret residual species-species covariance as biotic interaction (see Zurell 2018 Ecography 41:1812-1819; Poggiato 2021 Trends Ecol Evol 36:391-401).

library(Hmsc)

# HMSC for moderate-S Bayesian JSDM
# X: site x environment data
# Y: site x species presence-absence or abundance
# distr: 'probit' for presence-absence; 'lognormal poisson' for counts
m <- Hmsc(Y = species_matrix, XData = env_data, XFormula = ~ temperature + soil_pH,
         distr = 'probit')

# Sample posterior (production runs need thin >= 100, transient >= 1000, samples >= 1000)
m <- sampleMcmc(m, thin = 10, samples = 1000, transient = 5000, nChains = 4)

# Variance partitioning into fixed (environment), random (latent factors), traits, phylogeny
VP <- computeVariancePartitioning(m)
plotVariancePartitioning(m, VP = VP)

# Predict to new environment
pred <- predict(m, XData = new_env_data, expected = TRUE)

# For S > 500: switch to sjSDM
# library(sjSDM)
# m_sj <- sjSDM(Y = species_matrix, env = linear(env_data, ~ temperature + soil_pH),
#               family = binomial('probit'))
# summary(m_sj); plot(m_sj)

Indicator Species with Group-Size Correction

Goal: Identify species statistically associated with site groups using an indicator value that combines specificity and fidelity, corrected for unequal group sizes.

Approach: Run indicspecies::multipatt with func = 'IndVal.g' (the group-size-equalized form per De Caceres & Legendre 2009 Ecology 90:3566-3574). The original IndVal is biased toward larger groups; the .g form corrects this. For continuous-vs-categorical or rank-based associations, use func = 'r.g' (point-biserial correlation, group-equalized).

library(indicspecies)

# IndVal.g: group-size-equalized indicator value
# multipatt tests species-group associations with permutation
# duleg=TRUE: tests only single-group associations (not combinations)
# duleg=FALSE: tests species against all group combinations (more powerful but more tests)
mp <- multipatt(species_matrix, site_groups,
                func = 'IndVal.g',       # Group-equalized; NOT 'IndVal' (biased)
                duleg = TRUE,
                control = how(nperm = 999))

summary(mp, alpha = 0.05)

# Extract significant indicators sorted by p-value
sig <- mp$sign[!is.na(mp$sign$p.value) & mp$sign$p.value < 0.05, ]
sig[order(sig$p.value), ]

Mantel Replacement — dbRDA with Spatial Covariates

Goal: Test whether community/genetic distance correlates with environmental distance while controlling for spatial autocorrelation, replacing the low-power and bias-prone Mantel framework.

Approach: Use db-RDA with spatial predictors (PCNM eigenvectors or raw coordinates) as Condition() (Legendre & Fortin 2010 Mol Ecol Resour 10:831-844). Partial Mantel inflates Type I error under autocorrelation — do not use for landscape data.

library(vegan)
library(adespatial)

# Build spatial predictors: PCNM (principal coordinates of neighbor matrices)
geo_dist <- dist(coords[, c('longitude', 'latitude')])
pcnm <- pcnm(geo_dist)
pcnm_vars <- pcnm$vectors  # significant axes

# db-RDA controlling for spatial structure
dbrda_result <- dbrda(bray_dist ~ temperature + precipitation +
                       Condition(as.matrix(pcnm_vars)),
                      data = env_data, add = 'lingoes')

# Test marginal significance of environment AFTER conditioning out space
anova(dbrda_result, by = 'margin', permutations = 999)

# DO NOT use partial Mantel for this question; cite Legendre & Fortin 2010

Per-Method Failure Modes

PERMANOVA significant, betadisper also significant -> conclusion not supported

Trigger: Reporting PERMANOVA p < 0.05 as evidence of "community composition differs across groups" without running betadisper, OR running betadisper and finding it significant but ignoring the result.

Mechanism: PERMANOVA's pseudo-F responds to BOTH centroid shifts AND dispersion heterogeneity (Anderson & Walsh 2013). When groups differ in dispersion but not centroid, PERMANOVA can return p < 0.05 from the dispersion difference alone.

Symptom: Reviewer asks "was dispersion checked?"; PCoA visualization shows overlapping centroids but one group is more dispersed; betadisper p < 0.05.

Fix: Report both PERMANOVA and betadisper. If betadisper is significant, the conclusion must be reframed: "Groups differ in either centroid or dispersion, with dispersion heterogeneity present." Consider db-RDA which is more dispersion-robust.

Mantel test reports low p but reflects spatial autocorrelation

Trigger: Using Mantel or partial Mantel to test "is genetic distance correlated with environmental distance" in a spatially-structured landscape.

Mechanism: Mantel statistics are biased downward by spatial autocorrelation in either distance matrix (Legendre & Fortin 2010). Partial Mantel further inflates Type I error rates under autocorrelation (see Guillot & Rousset 2013 Methods Ecol Evol 4:336-344 for the formal demonstration).

Symptom: Mantel p < 0.001 but no detectable signal when re-tested with dbRDA conditioning on geographic distance.

Fix: Use db-RDA with PCNM spatial eigenvectors as Condition(), OR GDM (Generalized Dissimilarity Modeling).