bio-ecological-genomics-community-ecology

What this skill does

## 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

| Method | Estimand | Strength | Fails when |
|--------|----------|----------|------------|
| CCA | Species-environment unimodal | Standard for chi-square data; fits bell-shaped responses | Linear gradients (use RDA); does not handle short gradients well |
| RDA | Species-environment linear | High power with short gradients; Hellinger-friendly | Long unimodal gradients (>3 SD DCA axis 1); needs no missing data |
| db-RDA / capscale | Constrained ordination on any distance | Flexible (Bray-Curtis, Sorensen, weighted Unifrac); dispersion-robust | Less power than RDA for purely linear gradients |
| PERMANOVA (adonis2) | Centroid difference among groups | Non-parametric; handles any dissimilarity | Sensitive to dispersion difference — MUST run PERMDISP alongside |
| PERMDISP (betadisper) | Dispersion difference among groups | Tests the confound that contaminates PERMANOVA | Low power with small N; report alongside PERMANOVA |
| ANOSIM | Group-difference test | Legacy familiarity | Worse than PERMANOVA for the same use case; biased by unbalanced N |
| Mantel test | Correlation between two distance matrices | Conceptually simple | Low power under spatial autocorrelation; biased; replace with dbRDA or GDM |
| Partial Mantel | Correlation controlling for a third matrix | Conceptually simple | INFLATES Type I error under autocorrelation (worse than basic Mantel) |
| HMSC (Helsinki tradition) | Bayesian JSDM with phylogeny + traits + spatial | Rigorous; ecological-theory priors | Slow for S > 500 species |
| sjSDM (Pichler & Hartig) | Latent-variable-free JSDM via MC + elastic net | Orders of magnitude faster; high-S friendly | Less explicit interpretation than HMSC |
| gjam | Cross-data-type JSDM (counts, presence, continuous) | Integrates heterogeneous data | Different output structure than HMSC; not directly comparable |
| IndVal (Dufrene-Legendre) | Species-group association | Combines specificity AND fidelity | Biased by group size unless `func='IndVal.g'` |
| SES_MPD / SES_MNTD | Phylogenetic structure vs null | Tests whether communities are clustered or overdispersed | Interpretation of sign requires trait-conservatism check (Mayfield & Levine 2010) |
| RLQ + fourth-corner | Trait-environment relationship | Tests species-trait response to environment | Default modeltype gives inflated Type I; use modeltype=6 |

## Decision Tree by Scenario

| Scenario | Recommended approach | Why |
|----------|---------------------|-----|
| Unimodal species responses, gradient > 3 SD (DCA axis 1) | CCA | Linear assumption fails for long gradients |
| Linear species responses, gradient <= 3 SD | RDA on Hellinger-transformed data | RDA assumes linearity; Hellinger solves double-zero problem |
| Bray-Curtis or other non-Euclidean distance preferred | db-RDA | Handles any distance metric while preserving constrained-ordination interpretation |
| Test "do these groups differ in community composition" | PERMANOVA (adonis2) AND PERMDISP (betadisper) | Never PERMANOVA alone; dispersion confound is non-negotiable |
| ANOSIM out of habit | Use PERMANOVA + PERMDISP instead | ANOSIM is biased by unbalanced N and dispersion |
| Mantel test for landscape genetics | dbRDA with spatial covariates OR GDM | Mantel has low power; partial Mantel inflates Type I |
| Joint species modeling with traits and phylogeny | HMSC (S < 500); sjSDM (S > 500) | HMSC encodes theory in priors; sjSDM is the only scalable option for high-S |
| Multi-species "association" structure | JSDM residual covariance | Interpret as ANY of (interaction, shared response to unmeasured driver, dispersal, sampling), NOT pure interaction |
| Indicator species with unbalanced group sizes | `multipatt(..., func = 'IndVal.g')` | `IndVal` is biased by group size; group-equalized form corrects |
| Phylogenetic community structure | SES_MPD with EXPLICIT null model + trait-conservatism test | Cite Mayfield & Levine 2010 for the "clustering = filtering" interpretation trap |
| Trait-environment hypothesis | RLQ + fourth-corner with `modeltype=6` | Default `modeltype=2` or `modeltype=4` gives inflated Type I error |
| Bipartite network nestedness/modularity | NODF2 + modularity with `curveball` null randomization | Strona 2014 curveball is exponentially faster and unbiased for binary matrices |

## 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.

```r
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),