Version Compatibility

Reference examples tested with: iNEXT 3.0+, iNEXT.3D 1.0+, vegan 2.6+, betapart 1.6+, 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.

Biodiversity Metrics

"Calculate species diversity for my ecological samples" -> Compute Hill-number diversity (numbers-equivalent of richness, Shannon, Simpson) with coverage-based rarefaction/extrapolation, choose the richness estimator appropriate to the singleton/doubleton signature of the data, and partition beta diversity into turnover and nestedness components with a documented partition framework.

• R: iNEXT::iNEXT() for coverage-based rarefaction/extrapolation
• R: betapart::beta.multi() for Baselga turnover/nestedness partition
• R: picante::ses.mpd() for phylogenetic-community SES with explicit null model

The Single Most Important Modern Insight -- Standardize by COVERAGE not by sample size

Chao & Jost 2012 Ecology 93(12):2533-2547 established that comparing diversity across assemblages by rarefying to a common sample size systematically biases comparisons whenever assemblages differ in underlying diversity: a 100-read rarefaction of a 50-species community is essentially saturated (coverage approximately 99%) while the same 100 reads from a 500-species community covers only approximately 40% of the underlying diversity. The two rarefied diversities are not measuring the same thing. Coverage-based rarefaction with iNEXT is the postdoc-grade default; sample-size rarefaction is now considered a methodological anti-pattern for cross-site comparison.

A second insight pairs with this: raw Shannon and Simpson indices are entropies, NOT diversities. Jost 2006 Oikos 113(2):363-375 showed that only their numbers-equivalents (exp(H), 1/D) are comparable as effective species counts and have the intuitive "doubling property" (merging two equally diverse equally abundant assemblages doubles the diversity). Reporting raw Shannon = 3.2 vs 3.0 hides whether the difference is large or trivial; reporting 1D = 24.5 vs 20.1 species-equivalents makes the 22% gap visible.

Algorithmic Taxonomy

Method	Estimand	Strength	Fails when
Hill numbers q=0,1,2 (iNEXT)	Effective species count at order q	Unifies richness, Shannon, Simpson; doubling property	None; report all three q values
Chao1 = S + f1^2/(2*f2)	Lower bound on asymptotic richness	Non-parametric; works at low coverage	Singletons dominated by PCR error (amplicon data); f2 near zero
Chao1bc = S + f1(f1-1)/(2(f2+1))	Bias-corrected Chao1	Stable when f2 = 0	Same singleton-bias issue
ACE	Asymptotic richness using all rare classes (f1...f10)	Uses more rare-class information than Chao1	Choice of "rare" cutoff (default 10) is arbitrary
Jackknife1 = S + f1*(n-1)/n	Asymptotic richness via resampling	Robust when f2 = 0	Sensitive to singleton count alone
Coverage-based rarefaction (iNEXT)	Diversity at standardized completeness	Correct comparison across sites with unequal effort	Extrapolation beyond 2x reference size is unreliable
Sample-size rarefaction	Diversity at common n	Legacy familiar	Systematically biased when communities differ in true diversity
Faith's PD	Sum of branch lengths on spanning tree	Captures evolutionary distinctness	Sensitive to richness; report alongside SES_PD
Rao's Q	Pairwise functional/phylogenetic distance * abundance	Unifies taxonomic and functional/phylogenetic diversity	Trait dimensionality artifacts (see Maire 2015)
FRic / FEve / FDiv	Functional richness/evenness/divergence	Multidimensional trait coverage	FRic inflates with collinear traits; optimize axis count
Baselga beta partition	Sorensen = turnover (Simpson) + nestedness	Decomposes beta into ecological processes	Not unique; Podani partition gives different components
Podani / Carvalho partition	Alternative turnover + richness-difference	Conceptually distinct from Baselga	Same data, different conclusion possible; report both
Hellinger transform + PCA/RDA	Solves double-zero problem	Standard for community ordination	None; required preprocessing

Decision Tree by Scenario

Scenario	Recommended approach	Why
Cross-site diversity comparison with unequal sampling effort	Coverage-based rarefaction in iNEXT to common C (typically 0.95)	Sample-size rarefaction is biased when assemblages differ in true diversity
Reporting diversity for publication	Hill numbers q=0,1,2 (numbers-equivalents)	Comparable, additive, doubling property; raw entropies are not
Amplicon/eDNA data with many singletons	Skip Chao1 OR use after careful denoising; report Good's coverage alongside	Singletons in amplicon data are dominated by PCR error, not undersampling
Singleton-heavy real data (f1 >> f2)	Chao1 with wide CI; cross-check with ACE	Chao1 variance grows as f1^4; ACE uses more rare-class information
f2 = 0 (no doubletons)	Chao1bc or jackknife1	Original Chao1 is undefined; bias-corrected form is the standard
Diversity at much larger sample size than observed	Stop extrapolation at 2x reference size (the doubling rule)	iNEXT will silently extrapolate further; variance grows superlinearly beyond
Beta diversity for two assemblages	Sorensen or Jaccard with Baselga partition; document choice	Bray-Curtis is not a true metric and breaks some downstream methods
Beta diversity reported as turnover/nestedness	Run BOTH Baselga AND Podani partitions and report both	The partition is not mathematically unique; one alone hides ambiguity
Before PCA/RDA on community data	Hellinger transformation first	Solves double-zero problem; raw Euclidean PCA on counts is malpractice
Phylogenetic community structure (NRI/NTI)	SES_MPD with explicit null.model justification	Default null may not match question; document choice
Functional diversity (FRic)	Optimize PCoA axis count via Maire 2015 mAD; report chosen k	FRic biased by trait-axis count; defaults (k=2-3) are too few
Amplicon/compositional data	Either presence/absence metrics (Chao2, Sorensen) OR explicit CLR-based diversity per Gloor 2017 Front Microbiol 8:2224	Hill-numbers assume absolute abundances; amplicon counts are compositional

Hill Numbers and Why Raw Entropies Are Not Diversities

Goal: Report diversity values that are interpretable as effective species counts and additive under standard partitions.

Approach: Compute Hill numbers qD at q = 0, 1, 2 (Jost 2006 Oikos 113:363-375; iNEXT software paper Hsieh, Ma, Chao 2016 Methods Ecol Evol 7:1451-1456). q controls sensitivity to rare vs common species in a continuous parametric family: q=0 weights all species equally (richness); q=1 is the geometric-mean weighting (Shannon-equivalent exp(H)); q=2 weights down rare species rapidly (Simpson-equivalent 1/D). qD has the doubling property — merging two equally diverse equally abundant assemblages exactly doubles qD. Raw Shannon and Simpson values do NOT.

library(iNEXT)
library(vegan)

abundance_data <- list(
    site_A = c(100, 45, 23, 12, 8, 5, 3, 2, 1, 1),
    site_B = c(80, 60, 40, 30, 20, 15, 10, 8, 5, 3, 2, 1),
    site_C = c(200, 10, 5, 2, 1, 1, 1)
)

# Hill numbers q=0,1,2 with coverage-based rarefaction/extrapolation
# nboot=200 is the publication-quality floor; default nboot=50 is too few for CIs
result <- iNEXT(abundance_data, q = c(0, 1, 2), datatype = 'abundance', nboot = 200)

# Report effective species counts at standardized coverage (postdoc-grade default)
# coverage=0.95: 95% of individuals in the underlying community belong to detected species
est <- estimateD(abundance_data, datatype = 'abundance', base = 'coverage', level = 0.95)
est

# For vegan equivalence: numbers-equivalent of Shannon
shannon_eff <- exp(diversity(community_matrix, index = 'shannon'))
# Inverse Simpson is already in numbers-equivalent units
invsimp <- diversity(community_matrix, index = 'invsimpson')

Asymptotic Richness — Chao1 is a Lower Bound, NOT a Point Estimate

Goal: Estimate the lower bound on true richness when sampling is incomplete and report it correctly.

Approach: Compute Chao1 from the singleton/doubleton ratio (Chao 1984), check that singletons are real biology rather than PCR/sequencing error, and report as a LOWER BOUND with Good's coverage to indicate sampling adequacy. Switch to ACE or jackknife1 when f2 is near zero or singletons are unreliable.

The original Chao1 derivation (Chao 1984) under a Gamma-Poisson mixture model produces a NON-PARAMETRIC LOWER BOUND on richness, not a point estimate. The bound is tight only under specific homogeneity assumptions. Reporting "Chao1 estimated richness = 320" without "at least" is widespread in published literature but statistically wrong.

library(iNEXT)

# AsyEst returns Chao1 (q=0), Chao-Shannon (q=1), Chao-Simpson (q=2) with CIs
asymp <- iNEXT(abundance_data, q = c(0, 1, 2), datatype = 'abundance')$AsyEst

# CRITICAL: also report Good's coverage to indicate whether the bound is informative
# Coverage < 0.85 means heavily under-sampled; Chao1 CI will be huge and bound is uninformative
coverage <- estimateD(abundance_data, datatype = 'abundance',
                      base = 'coverage', level = 0.95)
# Inspect estimateD output's Coverage column at the OBSERVED sample size

# Singleton check: if f1 >> f2 and singletons are likely PCR artifacts (amplicon data),
# do NOT report Chao1 — it will measure "how much PCR error" not "how much undiscovered diversity"
f1_check <- sapply(abundance_data, function(x) sum(x == 1))
f2_check <- sapply(abundance_data, function(x) sum(x == 2))
cat('Singletons f1:', f1_check, '\n')
cat('Doubletons f2:', f2_check, '\n')
cat('f1^2/(2*f2) Chao1 bound term:', f1_check^2 / (2 * pmax(f2_check, 0.5)), '\n')

Coverage-Based Rarefaction with the Doubling Rule

Goal: Compare diversity across sites at a common sampling completeness, with extrapolation bounded by statistical reliability.

Approach: Use iNEXT's coverage-based rarefaction-extrapolation interpolating to the minimum coverage across sites; bound extrapolation at 2x the reference sample size (Chao et al. 2014 Ecol Monogr 84:45-67 derive variance bounds that grow superlinearly beyond 2x).

# Type 3 (diversity vs coverage) is the cross-site comparison plot
# nboot=200 minimum for publication-quality CIs (default 50 is too few)
ggiNEXT(result, type = 3) + theme_bw() +
    labs(title = 'Coverage-Based Hill-Number Diversity')

# The doubling rule: endpoint should be at most 2 * max(observed sample size)
# iNEXT default endpoint = 2 * max(sample size); do NOT override above this
# Beyond 2x, the variance estimate grows superlinearly and CIs become unreliable

# To extract numerical results at standardized coverage:
est_95 <- estimateD(abundance_data, datatype = 'abundance',
                    base = 'coverage', level = 0.95)

The Double-Zero Problem and Hellinger Transformation

Goal: Compute community dissimilarity in a way that does not treat shared absences as evidence of similarity.

Approach: Apply the Hellinger transformation (Legendre & Gallagher 2001 Oecologia 129:271-280) before computing Euclidean distance, PCA, or RDA. This is the single most important preprocessing decision for community ordination.

Standard Euclidean distance and Pearson correlation treat "both samples have 0 abundance of species X" as evidence of similarity. In community ecology this is wrong — two deserts both lacking a rainforest species are not thereby similar. Hellinger transformation removes this artifact while preserving total-abundance information.

library(vegan)

# Hellinger transformation: y_ij' = sqrt(y_ij / row_sum_i)
# Euclidean distance on Hellinger-transformed data = Hellinger distance
species_hell <- decostand(community_matrix, method = 'hellinger')

# Now PCA/RDA on the transformed matrix is biologically meaningful
pca_result <- rda(species_hell)

# Alternative: chord transformation (similar effect, slightly different scaling)
species_chord <- decostand(community_matrix, method = 'normalize')

Beta Diversity Partition — Run BOTH Baselga AND Podani

For the broader "multiple meanings of beta diversity" framework, see Anderson et al. 2011 Ecol Lett 14:19-28.

Goal: Decompose total beta diversity into ecologically interpretable components, acknowledging that the partition is not mathematically unique.

Approach: Compute the Baselga partition (Sorensen = turnover + nestedness) with betapart, and the alternative Podani/Carvalho partition (richness-difference framework), and report both. The two frameworks give different ecological interpretations of the same data; presenting only one hides ambiguity.

library(betapart)

pa_matrix <- ifelse(community_matrix > 0, 1, 0)

# --- Baselga partition (Baselga 2010 Glob Ecol Biogeogr 19:134-143) ---
# beta.SIM (turnover/Simpson) + beta.SNE (nestedness) = beta.SOR (total Sorensen)
pair_sor <- beta.pair(pa_matrix, index.family = 'sorensen')
multi_sor <- beta.multi(pa_matrix, index.family = 'sorensen')

# --- Abundance-based: Bray-Curtis balanced + gradient ---
# beta.bray.bal (balanced variation; analogous to turnover)
# beta.bray.gra (abundance gradient; analogous to nestedness)
pair_abund <- beta.pair.abund(community_matrix, index.family = 'bray')

# --- Podani/Carvalho framework (different decomposition of same data) ---
# Use betapart::beta.pair with the .fam family option, or carvalho package
# These produce richness-difference components instead of nestedness
# Reporting only Baselga without acknowledging Podani exists is incomplete

Phylogenetic Diversity — Faith's PD, MPD, MNTD with Null-Model Choice

Goal: Quantify evolutionary distinctness and phylogenetic community structure with an explicit, justified null model.

Approach: Compute Faith's PD on a community matrix and ultrametric phylogeny (Faith 1992 Biol Conserv 61:1-10), then SES_MPD and SES_MNTD (Webb 2002 Annu Rev Ecol Syst 33:475-505) standardizing against a documented null model. The choice of null model is the dominant scientific decision — different nulls answer different ecological questions.

library(picante)

# Faith's PD: sum of branch lengths spanning the focal species
# include.root=TRUE counts root branch; FALSE excludes (matters for within-clade comparisons)
faith_pd <- pd(community_matrix, phylo_tree, include.root = TRUE)

# SES_MPD with EXPLICIT null model (do not accept default silently)
# null.model='taxa.labels': shuffles species across tree, holds sample richness constant
# null.model='richness': randomizes within sample, preserves species occurrence frequencies
# null.model='independentswap': preserves both row and column sums of community matrix
# Each null answers a different question; document the choice in methods
ses_mpd_taxa <- ses.mpd(community_matrix, cophenetic(phylo_tree),
                        null.model = 'taxa.labels', runs = 999, iterations = 1000)
ses_mpd_indep <- ses.mpd(community_matrix, cophenetic(phylo_tree),
                         null.model = 'independentswap', runs = 999, iterations = 1000)

# SES_MPD < 0 (NRI > 0): phylogenetic clustering
# SES_MPD > 0 (NRI < 0): phylogenetic overdispersion
# DO NOT infer "clustering = environmental filtering" from sign alone;
# Mayfield & Levine 2010 Ecol Lett 13:1085-1093 showed competition can cluster
# when traits track phylogeny and similar species coexist via R*-rule dynamics

Functional Diversity with Trait-Axis Dimensionality Optimization

Goal: Compute multidimensional functional diversity without inflating values via correlated trait axes.

Approach: Run Maire et al. 2015 Glob Ecol Biogeogr 24:728-740 trait-space-quality assessment (mAD metric) to choose the number of PCoA axes that minimizes deviation between original trait distances and axis-reconstructed distances. Use that k for FRic, FEve, FDiv. Default k=2-3 typically biases FRic; optimum is usually 4-6 for typical trait datasets.

library(FD)
library(mFD)

# Trait dissimilarity from a trait matrix (Gower for mixed quantitative/categorical)
trait_dist <- gowdis(trait_matrix)

# mFD optimizes the number of trait axes via Maire 2015 mAD criterion
# Smaller mAD = better fidelity between trait distances and PCoA-axis distances
quality_funct_space <- quality.fspaces(trait_dist, maxdim_pcoa = 10,
                                       fdist_scaling = TRUE, fdendro = NULL)
best_k <- which.min(quality_funct_space$quality_fspaces$mad)
cat('Maire-optimal number of trait axes:', best_k, '\n')

# Compute FD with the optimal k
fd_result <- dbFD(trait_matrix, community_matrix, m = best_k, corr = 'cailliez')
fd_result$FRic   # functional richness (convex hull volume)
fd_result$FEve   # functional evenness
fd_result$FDiv   # functional divergence

Per-Method Failure Modes

Singleton-driven Chao1 inflation in amplicon data

Trigger: Computing Chao1 directly on ASV/OTU tables that include singletons of suspected PCR-error origin.

Mechanism: Chao1 assumes singletons are biologically real rare species; under that assumption, more singletons relative to doubletons signals more undiscovered diversity. PCR/sequencing error produces many singletons that look like rare species to Chao1, inflating the bound.

Symptom: Chao1 dramatically exceeds observed richness (e.g., Chao1 = 500 from S_obs = 100) with extremely wide confidence intervals; Good's coverage < 0.85 despite > 10,000 reads per sample.

Fix: Denoise first (DADA2 / UNOISE3 / swarm v2) to remove PCR-error variants, then either skip Chao1 entirely (Callahan 2017 ASV philosophy) or report observed ASV count + Good's coverage instead. Alternative: report Chao2 from incidence (presence across replicates), which is robust to PCR-error singletons.

Rarefaction across mismatched assemblage sizes

Trigger: Sample-size rarefaction to the smallest assemblage when sites differ in true diversity.

Mechanism: A 100-read rarefaction of a 50-species community is nearly saturated (coverage approximately 99%); the same 100 reads from a 500-species community covers approximately 40% of true diversity. The two rarefied values measure different completeness levels.

Symptom: Rarefied richness ordering disagrees with intuitive site-by-site comparison; small-sample-size sites appear artificially equal in rarefied diversity to high-diversity sites.

Fix: Use coverage-based rarefaction via estimateD(..., base = 'coverage', level = 0.95). The 0.95 coverage target is the modern default; 0.99 for high-precision work, 0.85 minimum for accepting a site into the comparison.

Hellinger forgotten before PCA/RDA on community data

Trigger: Running rda(species_matrix) or prcomp(species_matrix) on raw species counts without prior transformation.

Mechanism: The double-zero problem inflates dissimilarity for sample pairs sharing many absent species. Raw-count PCA projects samples primarily by sequencing depth (PC1 = library size proxy), not by community composition.

Symptom: PC1 axis correlates strongly with sample read totals; biological gradients appear only on PC3 or later; ordination is dominated by sites with extreme sample sizes.

Fix: decostand(matrix, method = 'hellinger') before ordination. Always.

FRic inflation from collinear trait axes

Trigger: Computing FRic from 8-10 raw correlated traits without dimensionality optimization.

Mechanism: FRic = convex-hull volume in trait PCoA space. Adding correlated traits increases the apparent dimensionality of the trait space without adding ecological information; convex-hull volume inflates accordingly.

Symptom: FRic increases when adding new correlated traits to the analysis; FRic comparisons across studies that use different trait counts are inconsistent.

Fix: Run mFD::quality.fspaces() and use the mAD-optimal k. Report k in methods.

SES_MPD interpretation by sign alone

Trigger: Reporting "SES_MPD < 0 indicates environmental filtering" without testing alternative explanations.

Mechanism: Mayfield & Levine 2010 Ecol Lett 13(9):1085-1093 showed competition can produce phylogenetic clustering when ecologically similar species (which tend to be closely related) coexist via R*-rule trait differences. The clustering = filtering / overdispersion = competition dichotomy from Webb 2002 is incomplete.

Symptom: Phylogenetic clustering interpretation is challenged at review; trait-similarity vs phylogenetic-similarity correlation has not been examined.

Fix: Report SES_MPD with one explicit null model AND test trait conservatism (Blomberg's K, Pagel's lambda) — if traits track phylogeny strongly, clustering may reflect either filtering OR competition; cite Mayfield & Levine 2010 in interpretation.

Quantitative Thresholds

Threshold	Value	Source / rationale
Coverage target for cross-site comparison	C = 0.95	Postdoc-grade convention; 0.99 for high-precision, 0.85 minimum to accept site
iNEXT extrapolation limit	endpoint <= 2x reference sample size	Chao et al. 2014 doubling rule; variance grows superlinearly beyond
iNEXT bootstrap floor for CIs	nboot = 200	Default 50 is too few for publication-quality CIs
Chao1 reliability	f2 > 0 AND singletons biological	If f2 = 0, switch to Chao1bc or jackknife1
Hellinger before ordination	Always for community data	Legendre & Gallagher 2001; non-negotiable
SES_MPD significance		SES
Functional-diversity axis count k	Maire 2015 mAD-optimal	Default k=2-3 too few; typically 4-6 optimal
Phylogenetic-tree requirement	Ultrametric	If not, ape::chronos() or BEAST-based dating

Common errors

Error	Cause	Solution
Chao1 reports infinity or NaN	f2 = 0 (no doubletons)	Use Chao1bc form or jackknife1
iNEXT extrapolation curve flat then explodes	Extrapolated beyond doubling-rule limit	Set endpoint <= 2 * max(sample sizes)
Hellinger PCA gives flat results	Decostand not applied or applied to wrong axis	decostand(matrix, method = 'hellinger') rows = sites
ses.mpd returns all p-values approximately 0.5	Wrong null model for question (e.g., taxa.labels on a single-richness dataset)	Choose null that varies the quantity tested
FRic dramatically increases with new traits	Correlated trait axes inflating convex hull	Optimize axis count with mFD::quality.fspaces
beta.multi returns NaN for nestedness	Communities entirely disjoint (no shared species)	beta.SNE -> 0 trivially; interpret as pure turnover
Bray-Curtis flagged for triangle-inequality violation	Bray-Curtis is not a true metric	Switch to Sorensen (a metric) for downstream methods requiring metricity

References

• Chao A, Jost L (2012) Coverage-based rarefaction and extrapolation. Ecology 93(12):2533-2547. doi:10.1890/11-1952.1
• Chao A, Gotelli NJ, Hsieh TC, Sander EL, Ma KH, Colwell RK, Ellison AM (2014) Rarefaction and extrapolation with Hill numbers. Ecol Monogr 84(1):45-67. doi:10.1890/13-0133.1
• Hsieh TC, Ma KH, Chao A (2016) iNEXT: an R package for rarefaction and extrapolation of species diversity. Methods Ecol Evol 7(12):1451-1456. doi:10.1111/2041-210X.12613
• Jost L (2006) Entropy and diversity. Oikos 113(2):363-375. doi:10.1111/j.2006.0030-1299.14714.x
• Faith DP (1992) Conservation evaluation and phylogenetic diversity. Biol Conserv 61(1):1-10. doi:10.1016/0006-3207(92)91201-3
• Legendre P, Gallagher ED (2001) Ecologically meaningful transformations for ordination. Oecologia 129(2):271-280. doi:10.1007/s004420100716
• Baselga A (2010) Partitioning the turnover and nestedness components of beta diversity. Glob Ecol Biogeogr 19(1):134-143. doi:10.1111/j.1466-8238.2009.00490.x
• Anderson MJ, Crist TO, Chase JM et al. (2011) Navigating the multiple meanings of beta diversity. Ecol Lett 14(1):19-28. doi:10.1111/j.1461-0248.2010.01552.x
• Webb CO, Ackerly DD, McPeek MA, Donoghue MJ (2002) Phylogenies and community ecology. Annu Rev Ecol Syst 33:475-505. doi:10.1146/annurev.ecolsys.33.010802.150448
• Mayfield MM, Levine JM (2010) Opposing effects of competitive exclusion on phylogenetic community structure. Ecol Lett 13(9):1085-1093. doi:10.1111/j.1461-0248.2010.01509.x
• Maire E, Grenouillet G, Brosse S, Villeger S (2015) How many dimensions are needed to accurately assess functional diversity? Glob Ecol Biogeogr 24(6):728-740. doi:10.1111/geb.12299
• Chao A (1984) Nonparametric estimation of the number of classes in a population. Scand J Stat 11(4):265-270
• Gloor GB, Macklaim JM, Pawlowsky-Glahn V, Egozcue JJ (2017) Microbiome datasets are compositional: and this is not optional. Front Microbiol 8:2224. doi:10.3389/fmicb.2017.02224

Related Skills

• ecological-genomics/edna-metabarcoding - Generate ASV/species tables prior to diversity analysis
• ecological-genomics/community-ecology - Constrained ordination, indicator species, PERMANOVA on transformed data
• microbiome/diversity-analysis - 16S clinical microbiome diversity metrics with compositional considerations
• data-visualization/ggplot2-fundamentals - Customize diversity plots and rarefaction curves
• phylogenetics/tree-io - Ultrametric tree preparation for PD/MPD/MNTD