Build a presence/absence mapping from OTU abundance tables. For each OTU, records the set of sample names in which it occurs (abundance > 0).

Dictionary mapping OTU_ID to the set of sample names where present.

Extract all unique sample names from the OTU tables.

Compute the occupancy (number of sites) for each species (OTU). Occupancy Rⱼ is the number of sites in which species j is present.

Compute the binomial coefficient C(n, k) = n! / (k! * (n-k)!). Uses an iterative multiplicative approach to avoid computing large factorials directly. Returns 0 if k < 0 or k > n.

Compute log(C(n, k)) for numerical stability with large values. log(C(n,k)) = Σ_{i=1}^{k} [log(n-k+i) - log(i)]

Log-sum-exp computation for numerically stable summation of potentially very large/small values in log space. log(Σ exp(xᵢ)) = M + log(Σ exp(xᵢ - M)) where M = max(xᵢ)

Compute zeta diversity of a specific order i.

ζᵢ = [ Σⱼ C(Rⱼ, i) ] / C(n, i)

where:

• Rⱼ = occupancy of species j (number of sites where present)
• n = total number of sites
• C(a, b) = binomial coefficient

This is the efficient formula from Hui & McGeoch (2014) that avoids explicitly enumerating all C(n, i) site combinations. The key insight: a species with occupancy Rⱼ contributes to C(Rⱼ, i) of the C(n, i) possible combinations of i sites.

ζᵢ value. Returns 0 if order > number of sites.

Compute zeta diversity using log-space arithmetic for numerical stability with very large numbers of sites or high occupancies.

log(ζᵢ) = log(Σⱼ C(Rⱼ, i)) - log(C(n, i)) = LogSumExp(log(C(Rⱼ, i))) - LogCombination(n, i)

Compute the full zeta diversity series from order 1 to maxOrder. Returns ζ₁, ζ₂, ..., ζₘₐₓₒᵣdₑᵣ.

This is the "zeta decline" - the fundamental pattern that reveals community assembly processes (Hui & McGeoch 2014).

Array where index i-1 contains ζᵢ.

Compute zeta diversity for a specific combination of sites (not averaged). This counts the number of species present in ALL of the specified sites.

Useful for understanding the distribution of zeta values across different site combinations, and for distance-decay analysis.

Number of species present in all specified sites.

Compute zeta ratios: ζᵢ/ζᵢ₋₁ for i = 2, 3, ..., n.

The zeta ratio represents the proportion of species shared among (i-1) sites that are also shared when one more site is added. It indicates the "retention rate" of species as more assemblages are considered.

Ecological interpretation (Hui & McGeoch 2014):

• Increasing ratios suggest structured communities (nestedness)
• Decreasing ratios suggest high turnover
• Constant ratios suggest random assembly

Compute normalized zeta: ζ̃ᵢ = ζᵢ / ζ₁ for all i.

Normalization removes the effect of differences in mean species richness (ζ₁), allowing comparison of zeta decline patterns across communities with different richness levels.

ζ̃₁ = 1.0 always, and ζ̃ᵢ decreases from 1.0 as i increases.

Compute return rates. The return rate at order i is defined as the zeta ratio ζᵢ/ζᵢ₋₁, representing the probability that a species found in (i-1) sites will also be found in the i-th site.

This is equivalent to ComputeZetaRatios but provided separately for conceptual clarity (Hui & McGeoch 2014 discuss return rates in the context of species turnover and metapopulation dynamics).

Compute the commonness index: the area under the normalized zeta decline curve.

Commonness = Σᵢ (ζᵢ / ζ₁) for i = 1 to n

A higher commonness index indicates that more species are widespread (common across many sites), while a lower value indicates that species are rare (restricted to few sites).

The minimum value is 1.0 (only ζ₁ contributes, all higher orders = 0), and the maximum value is n (all species in all sites, ζᵢ = ζ₁ for all i).

Compute Whittaker's beta diversity from zeta diversity.

β = (ζ₁ / ζ₂) - 1

This is the multiplicative beta diversity for pairs of sites, measuring how many times the species composition changes between two typical sites. Higher values indicate greater turnover.

Compute the rarity index based on the exponential decay rate of the zeta decline. A higher decay rate (steeper decline) indicates more rare species (species restricted to few sites).

RarityIndex = b (from exponential fit ζᵢ = a·exp(-b·(i-1)))

Perform ordinary least squares (OLS) linear regression. Fits y = intercept + slope * x. Returns (slope, intercept, r2, sse). Uses only basic Math functions.

Fit zeta decline with a power law model:

ζᵢ = a · i^(-b)

Linearized form: log(ζᵢ) = log(a) - b · log(i) Regression: y = log(ζ), x = log(i)

Ecological meaning (Hui & McGeoch 2014): Power law decline suggests stochastic/neutral community assembly, consistent with log-normal or log-series species abundance distributions. The exponent b relates to the diversity parameter (Fisher's alpha) and species turnover rate.

Fit zeta decline with an exponential model:

ζᵢ = a · exp(-b · (i-1))

Linearized form: log(ζᵢ) = log(a) - b · (i-1) Regression: y = log(ζ), x = (i-1)

Ecological meaning (Hui & McGeoch 2014): Exponential decline suggests niche-based community assembly, consistent with geometric series species abundance distributions. This pattern indicates environmental filtering or competitive hierarchies structuring the community.

Fit zeta decline with a linear model:

ζᵢ = a - b · (i-1)

Regression: y = ζ, x = (i-1)

A linear decline (in arithmetic space) is less commonly observed but can occur in communities with uniform species loss rates.

Select the best-fitting model based on AICc (corrected Akaike Information Criterion). AICc is preferred over AIC for small sample sizes, which is common in zeta diversity analysis where the number of data points equals the number of orders computed.

Returns the name of the best model and its AICc value.

Select the best-fitting model based on R² (coefficient of determination). Higher R² is better.

Perform a permutation test using a frequency-based null model.

The null model randomizes the presence/absence matrix while preserving each species' occupancy (number of sites where present). This tests whether the observed zeta diversity pattern deviates from what would be expected under random spatial assignment of species with the same frequencies.

Method (Hui & McGeoch 2014):

1. For each species, randomly reassign its occurrences to sites (preserving the number of sites it occupies).
2. Recompute the zeta series.
3. Repeat for numPermutations iterations.
4. Compute mean, standard deviation, p-values, and standardized effect sizes (SES) from the null distribution.

Randomly select 'count' items from a list without replacement using Fisher-Yates partial shuffle.

Generate a random number from a Binomial(n, p) distribution using the simple sum-of-Bernoullis method.

Perform complete zeta diversity analysis on OTU table data.

This function computes:

1. The full zeta series (ζ₁ to ζₘₐₓₒᵣdₑᵣ)
2. Zeta ratios (ζᵢ/ζᵢ₋₁)
3. Normalized zeta (ζᵢ/ζ₁)
4. Three model fits (power law, exponential, linear) with AICc/BIC
5. Best model selection
6. Commonness index
7. Beta diversity
8. Optional null model significance testing

Complete ZetaAnalysisResult with all metrics.

Format a ZetaAnalysisResult as a human-readable string for console output or logging.