Community composition - ordination, RDA and variation partitioning, PERMANOVA with dispersion test, beta-diversity partitioning, indicator species

1. Small designs (about 24 sites or fewer).

   • The jackknife intervals for the shared fraction (coverage 0.855) and for mean nestedness (0.885) are too narrow.
   • The unique fraction of a weak gradient is biased upward by about 10%.
   • With fewer than about 30 sites, report the shared fraction and nestedness as descriptive values without interval claims.

2. Strongly heteroscedastic sites (one group 3× noisier; Section 4d).

   • Jackknife coverage drops to 0.84–0.86 for the marginal R2 of the weaker gradient and for the shared fraction.
   • The PERMANOVA R2 is biased in whichever direction the group-size and dispersion imbalance pushes it.
   • When PERMDISP is significant, treat these intervals and the PERMANOVA R2 with caution.
3. The shared fraction at a true value of 0 has jackknife coverage of 0.83, and the shared fraction has no test (Peres-Neto et al. 2006).

4. Jackknife CIs are not tests.

   • The marginal R2adj CI excludes 0 in 8.5% of null datasets.
   • The PERMANOVA R2 CI almost never excludes 0 under weak effects.
   • Use the permutation p-values for hypothesis tests and the CIs for effect sizes.
5. PERMANOVA confounds location and dispersion in unbalanced designs (validated in Section 4c). This is a property of the method, so the module reports PERMDISP alongside it.
6. No detection or effort model. Truths and estimates are for the community as sampled. Unequal effort among sites shows up as composition differences and nestedness.
7. No restricted permutations for nested or blocked designs (plots within stands). Sites are assumed to be exchangeable.
8. Descriptive outputs have no recovery test: NMDS, DCA, CCA, db-RDA, multi-site β, the dispersion ratio and indicator values. They are checked by the tests in Section 4c and by inspection only.

Module validated 2026-10-19. R 4.5.3, vegan 2.7, betapart, indicspecies. All runs used SLURM on njobvu (16 cores). The raw outputs are in results/community/validation/<run>/ (*_summary.csv, *_reps.csv, HTML report) and the logs are in logs/community_*.log.

This module is descriptive and multivariate. Its estimands are population values of statistics, not parameters of a likelihood. Each truth in simulate() is the value the statistic converges to for an infinite set of sites drawn from the same process: the same species pool, the same distribution of environmental gradients, and groups in the design's proportions. It is computed by large Monte Carlo samples that are independent of the analysed dataset (seed + 7,777,777). See simulate.R.

All truths describe the community as sampled (Poisson counts at the given effort, no detection model). An undersampled survey has more apparent turnover than the latent community, and the truth reflects that by construction.

Each replicate draws a new species pool, so emp_sd in the summaries includes variation in the truth between replicates. Bias and coverage are computed per replicate against that replicate's truth.

Uncertainty: delete-one-site jackknife SE with t_{n−1} 95% intervals for every estimand that has a truth.

• Species: 30 species with Gaussian niches on two gradients.

  • Tolerances: 0.2 (elevation) and 0.4 (soil moisture) of the sampled range, each × logN(0, 0.3²).
  • Peak abundance: logN(log 4, 0.8²).
• Groups: 3 habitat groups with species-specific location effects δ ~ N(0, 0.5²) on log abundance.
• Noise: mean-preserving site × species lognormal noise, SD 0.5.
• Design: 60 sites, gradient correlation 0.4, Poisson counts.
• Recovery uses fast_settings: estimands only, with no permutation tests or ordinations. One replicate takes 25–110 s, mostly the truth Monte Carlo.

The contract passed. With n = 4 some coverage flags fail, as expected at that size (results/community/selftest/).

results/community/validation/default/ (re-run after the engine.R changes as default_v2/; the summary is byte-identical, as expected since dropping empty species does not change any estimand)

200/200 replicates ran (there is no MCMC; "converged" is always TRUE). The coverage bound at n = 200 is 0.904.

Failures in hard_small (coverage bound 0.904):

• R2adj_shared: coverage 0.855. The bias is fine (0.030 vs truth 0.034).
• beta_sne_pair_mean: coverage 0.885. Rel. bias −2.7%.
• R2adj_unique_soil_moisture: rel. bias +10.0% (0.048 vs 0.0435). The tolerance was max(10% = 0.00435, 3 MCSE = 0.0045), so this is a borderline failure. Coverage 0.97.

All other estimands passed:

• R2adj_elevation, R2adj_env_total, R2adj_unique_elevation and R2adj_soil_moisture: coverage 0.945–0.965, abs. rel. bias ≤ 1.9%.
• PERMANOVA R2: coverage 0.97.
• beta_sor and beta_sim: coverage 0.97–0.98.
• turnover_share: coverage 0.905.

Interpretation: with about 24 sites the delete-one-site jackknife is too narrow for the shared fraction and for nestedness. See Known limits.

4a. False-positive rate with effects = 0 (null_effects, nsim = 200)

Settings:

• No group effect: group_effect_sd = 0.
• Soil moisture has no effect on any species (flat niches, tol2 = 1000).
• Gradients are uncorrelated: env_cor = 0.

Truth is about 0 for R2adj_soil_moisture, the unique soil-moisture fraction, the shared fraction and PERMANOVA.

Gate: 10 of 11 pass. The shared fraction fails on coverage at a true value of 0, and the marginal soil-moisture jackknife CI has an FPR of 8.5%. For deciding whether a fraction is non-zero, use the permutation test (the p_value column), not the jackknife CI. The permutation tests' FPRs are in 4c.

4b. The gate can fail: deliberately wrong estimator (misspec_unadjusted_small, nsim = 100)

This run reports unadjusted R² (the vegan rda() summary default) as if it estimated the population share, at 24 sites. The gate FAILS 4 of 11:

The upward bias is about p/(n − 1), as theory predicts. The β-diversity estimands, which this misspecification does not touch, still pass.

4c. PERMANOVA vs PERMDISP: location vs dispersion (dev/community/permanova_rates*.R)

Each scenario has 200 simulated datasets, fitted with full settings:

• adonis2 PERMANOVA (999 permutations)
• betadisper PERMDISP (999 permutations)
• jackknife CI for PERMANOVA R2adj
• IndVal.g

The "n" column is the number of reps that ran. About 4% failed, all because a simulated species was absent from every site, which crashed the indicator table. That was fixed afterwards (see Changes). The failures do not depend on the test outcome. The cells are rejection rates at α = 0.05.

The ×3 rows come from dev/community/permanova_rates_disp3.R and results/community/validation/permanova_rates_disp3/. The mean PERMANOVA R2adj was:

• 0.0160 against a truth of 0.0113 (small group more variable): inflated by 42%.
• 0.0056 against a truth of 0.0090 (large group more variable): deflated by 38%.
• 0.0158 against a truth of 0.0157 (balanced): unbiased.

This is the direction Anderson & Walsh (2013) describe: PERMANOVA is liberal when the smaller group is more dispersed and conservative when the larger group is.

Findings:

• Permutation tests hold their level under the null.

  • PERMANOVA rejects 3.7% (2 groups) and 5.2% (3 groups).
  • PERMDISP rejects 4.7% and 6.8%.
  • Indicator species (any species significant after BH correction) rejects 1.0%.

• Dispersion-only differences produce PERMANOVA "differences in composition".

  • Balanced design: the rate stays near nominal (5.2%).
  • Unbalanced, with the small group the more variable one: PERMANOVA rejects 9.3% at ×2 noise. At ×3 noise the R2 estimate is inflated by 42% when the small group is the more variable one, and deflated by 38% in the mirror design.
  • PERMDISP flags the dispersion difference in 40–79% of datasets at ×2 and in 92–100% at ×3. This is why the module always reports PERMDISP next to PERMANOVA, and why interpret() warns when both are significant, naming the unbalanced case (Anderson & Walsh 2013).
• Simulator caveat (honest). In this simulator, "dispersion only" means extra mean-preserving lognormal noise on expected counts. That adds zeros, so the Bray-Curtis population centroid also moves slightly: truth PERMANOVA R2 is 0.002–0.003 at ×2 and 0.009–0.016 at ×3, not 0. With differing noise, the zero-inflation changes the expected Bray-Curtis profile itself. Part of the excess PERMANOVA rejection is therefore a real but tiny location difference. The ×2 balanced scenario (5.2% rejection) shows the shift is too small to detect on its own.
• The jackknife CI for PERMANOVA R2 is very conservative for detection. In the location scenarios it never excluded 0, while the permutation test rejected 13–17%. Its coverage is still 0.95–0.98. Use the PERMANOVA p-value for testing and the jackknife CI as an interval for the effect size.

4d. Gate run: dispersion-only groups in an unbalanced design (dispersion_unbalanced, nsim = 100)

The PERMANOVA negative control, run through core/run.R recover.

Settings:

• No species-specific group effect (group_effect_sd = 0).
• The small group is 3× noisier (dispersion_mult = 3).
• Group sizes are 48 and 12.

The gate FAILS on the PERMANOVA R2, as intended:

The PERMANOVA effect size takes in the dispersion difference. In the matching 200-rep rate study (Section 4c, row "×3, 48/12"), PERMDISP flags it in 91.5% of datasets.

The two coverage failures are a side finding. With strongly heteroscedastic noise among sites, the jackknife intervals for the marginal soil-moisture R2 and for the shared fraction are too narrow (coverage 0.84–0.86). See Known limits.

dev/community/crosscheck.R was run on 5 simulated datasets at the default design. It compares the module's hand-written estimators, which give the point estimates and make the vectorised jackknife possible, against other packages:

Note: in vegan 2.7, varpart(...)$part$indfract orders the rows as [a] = X1|X2, [b] = X2|X1, [c] = shared. An earlier version of the cross-check script assumed the order [a], shared, [c]. That mismatch, not the module, caused the apparent 0.07 differences; it was corrected.

results/community/validation/crosscheck/crosscheck.csv

core/run.R fit community --config modules/community/example/config.yaml

The data are simulated bird point counts: 60 stations, 30 species and 3 stand types. The truth is in example/example_truth.rds.

All 11 intervals contain the truth. Other results:

• PERMDISP p = 0.78, even though young plantations were simulated with 1.4× noise. This is a low-power, honest non-detection, and the interpretation says "did not differ detectably".
• NMDS stress: 0.163.
• DCA axis 1: 4.6 SD.

The HTML report's interpretation sentences were checked against the table. All 11 figures were rendered and inspected:

• Fixed during review: overlapping labels, clipped axis titles, an empty indicator panel at single-column width, and default ggplot colours in layer-only mappings.
• The DCA sentence was changed, because it wrongly called CCA "the better summary". The Hellinger transformation makes RDA appropriate for long gradients (Legendre & Gallagher 2001).

Power is the share of 95% jackknife CIs excluding 0, at default effect sizes.

What this means for design:

• A strong gradient (about 30% of variation) is detected with 15 sites.
• A moderate one (about 10%) needs about 36 sites for 90% power.
• A weak group effect (truth R2 about 0.04) is rarely excluded from 0 by the conservative jackknife CI. Plan group designs with the PERMANOVA permutation test instead (Section 4c).

The lowest coverages (0.85–0.88) belong to R2adj_shared and beta_sne, consistent with Section 3. results/community/validation/power/

• output.R

  • Vectorised fmt_p. The report crashed on pairwise and indicator tables.
  • Fixed figure sizes, palettes and labels.
  • Corrected the DCA/CCA sentence.
  • Methods text now uses variable labels and spells out "Bray-Curtis".

• engine.R

  • Species with no records are dropped at the start of fit(); this changes no estimand.
  • Made the indicator-table group lookup robust. Without these two changes, about 4% of simulated datasets crashed with full settings.
• dev/community/crosscheck.R: corrected the vegan/ade4 fraction indexing.

The verifier wrote its own simulator and truth code in dev/community/verify/verify_indep.R. Only module$fit() is used from the module; the gate is core's summarise_recovery + recovery_gate.

The verifier's simulator differs from the builder's in every structural choice:

• Counts are negative binomial (size 1.5), not Poisson-lognormal.
• Temperature is Normal (12 + 4·z, °C). The niche is Gaussian in temperature only.
• Canopy (%) = 100·Φ(z2), and species respond log-linearly to canopy, not unimodally.
• Gradient correlation is 0.5.
• 40 sites in 4 unbalanced regions (13/9/9/9), 22 species, region effect SD 0.6.

Truth definitions:

• R2 truths: population R2 of the Hellinger matrix on the reported variables, from 50,000 sites.
• PERMANOVA truth: 1 − Σ π_g E[d²|g,g] / E[d²|mixture], from 40,000 pairs per term.
• β truths: mean over 80,000 random pairs of distinct design sites.

Commands: bin/ejob -c 16 -m 32 -- bin/Rscript dev/community/verify/verify_indep.R <scenario> <nsim>. Results are in results/community/verify/<scenario>/summary_gate.csv.

Larger false-positive study (dev/community/verify/null_rates.R, 1000 datasets, 40 sites, canopy and 3-level region independent of the community, 499 permutations):

The jackknife CI excludes 0 in 2.0%, 1.1%, 5.4% and 1.3% of datasets (marginal, unique, shared, PERMANOVA).

Rerun of the builder's default recovery (results/community/verify/default_rerun, nsim 200) reproduces Section 2 line for line, e.g. elevation 0.3196 vs truth 0.3211, coverage 0.945.

Literature and implementation checks:

• Baselga (2010) components agree with betapart::beta.pair source: β_sne = ((max−min)/(2a+b+c))·(a/(a+min)) = β_sor − β_sim.
• Adjusted R2 (Ezekiel) and the varpart fraction algebra follow Peres-Neto et al. (2006).
• PERMANOVA SS_T = (1/N)Σd², SS_W = Σ_g (1/n_g)Σ_g d² (Anderson 2001).
• dca_axis1_length_sd equals vegan's printed DCA axis length (4.6377 on the example).
• dev/community/verify/xcheck_indep.R compares the module with vegan::varpart, adonis2 and betapart::beta.pair on 3 datasets from the verifier's simulator. Maximum absolute difference: 3.3e-16.

Usability:

• fit on the example config: all 11 truths lie inside their 95% CIs.
• Every interpretation sentence was checked against the estimates. Correlational wording is used, and the as-sampled (no detection model) caveat is present.
• The five required figures (ordination, variance_partition, beta_components, dispersion, indicator_species) plus six optional ones were inspected. Plotted values equal fit$estimates, and axes are in original units.

• Messy data (dev/community/verify/messy/) gives clear, actionable messages:

  • misspelled site id: "in the species file but not in the environment file. Check spelling...";
  • NA cells: "Use 0 for 'not recorded'";
  • NA covariate: names the site;
  • wrong column name: lists the available columns;
  • species-code typos (PSFL/PSFLx): flagged.

Changes made by the verifier:

• io.R

  • Clearer messages for wrong site or variable columns (they list the actual columns).
  • Typo detection for near-duplicate species names.
  • Species and level labels from the config figures: labels block now reach the report text (the indicator sentence said "GCKI" instead of "Golden-crowned Kinglet").

• engine.R

  • read_data notes reach the report.
  • fit() no longer resets the caller's global random stream. It still seeds its own permutations with settings$seed, and the estimates are unchanged. The old behaviour made consecutive simulate-then-fit loops reuse random numbers.

• output.R

  • The PERMANOVA sentence says "did not differ detectably" when p ≥ 0.05.
  • The PERMDISP-not-significant sentence notes PERMDISP's limited power.
  • The IndVal A definition is corrected (group-size-equalised).
  • Removed an unsupported conservation recommendation from the turnover sentence.
  • The marginal PERMANOVA sentence now gives the unadjusted R2 without covariates for comparison.
  • RDA/NMDS arrow labels no longer collide with species labels.
• None of the changes alter any estimate, so the earlier recovery runs stand.

Additional limits found:

• PERMDISP is mildly liberal (7.2% ± 0.8% at α = 0.05 with groups of about 13 sites).
• With presence/absence data, the PERMANOVA R2 point estimate runs about 25% high at 40 sites. Its CI is conservative (coverage 0.98).
• The shared fraction runs about 9% low at 40 sites when it is small.

Verdict: VERIFIED_WITH_LIMITS. The default design passes. The documented limits (small designs, heteroscedasticity, shared fraction at 0, PERMDISP slightly liberal) remain.