diff --git a/DevDocMD/theory/evolutionary_distances.md b/DevDocMD/theory/evolutionary_distances.md index 7b1a3506..6b804241 100644 --- a/DevDocMD/theory/evolutionary_distances.md +++ b/DevDocMD/theory/evolutionary_distances.md @@ -2249,9 +2249,85 @@ F84/TN93 and is included below. `snp-` prefix on every CLI value. **`+Γ` rate-heterogeneity modifier, applicable to `snp-jc`, `snp-k2p`, `snp-k81`, `snp-t92`, `snp-f84`, `snp-hky85`, `snp-tn93`** (not `snp-raw`, nothing to correct; not `snp-logdet`, no standard gamma formulation) — same -formula as the base correction, weighted by a shape parameter `α` supplied -by the user (`--gamma-shape `), not estimated by ML. A modifier on -existing values, not a separate enum arm per distance. +formula as the base correction, weighted by a shape parameter `α` either +supplied by the user (`--gamma-shape `) or estimated from the data +(`--gamma-shape auto`/`estimate`, method-of-moments — not ML; see +"Automatic α estimation" below). A modifier on existing values, not a +separate enum arm per distance. + +### Automatic α estimation (`--gamma-shape auto`) + +**Correction (verified against the primary source, 2026-09-11):** Jin & +Nei (1990) itself (*"Limitations of the Evolutionary Parsimony Method of +Phylogenetic Analysis"*, Mol. Biol. Evol. 7(2):82–102 — the paper this +whole `+Γ` correction is cited from, confirmed algebraically to match this +codebase's `corrected_log`/`k2p` exactly against their eq. A4, general, and +A8, the `a = 1` case) contains **no** data-driven α-estimation procedure. +Their own recommendation (p. 98) is a fixed default: *"we suggest that the +gamma distance with a = 1 [eq. A8] be used. However, one may choose a +different gamma distance, estimating a from data. Wilson et al. (1989) +recently used a distance with a = 1/2 for restriction-site data of +mitochondrial DNA in hominoids."* — i.e. Jin & Nei explicitly punt +data-driven estimation to a *different* paper (Wilson et al. 1989), not +read/verified here. The estimator below is therefore **not** "Jin & Nei's +method" under any framing — that attribution (present in an earlier +revision of this section) was wrong, not just under-cited. + +**Implemented** (`PartitionDispersion`, `obikphylo/src/siblings/algorithms/pairwise.rs`) +as an independent method-of-moments estimator, unrelated to any specific +published procedure: pools substitution counts by **partition** rather +than by genome pair, during the same `reduce_pairwise` pass that builds +`PairwiseTally` (no second scan). + +For partition `i`: `n_i` = substitutions pooled over every genome pair, +`L_i` = eligible loci pooled over every genome pair, `R_i = n_i / L_i`. +Modeling among-site rate heterogeneity the same way as the `+Γ` correction +itself (a `Gamma(α, α)`-distributed, mean-1, multiplicative rate shared by +every locus in a partition — the classical Poisson-Gamma/negative-binomial +mixture, the general identity behind gamma-rate-heterogeneity corrections, +also behind Uzzell & Corbin's (1971) original observation that substitution +counts across sites/regions are over-dispersed relative to Poisson): + +\[ +\mathbb{E}[R_i] = \mu \qquad \mathrm{Var}[R_i] = \frac{\mu}{L_i} + \frac{\mu^2}{\alpha} +\] + +Weighting each partition's squared deviation by its own `L_i` cancels the +Poisson term before attributing what's left to `α`: + +\[ +\hat\mu = \frac{\sum_i n_i}{\sum_i L_i} \qquad +V = \frac{\sum_i L_i (R_i-\hat\mu)^2}{\sum_i L_i} \qquad +\hat\alpha = \frac{\hat\mu^2}{V - \hat\mu/\bar L} +\] + +where `\bar L` is the mean partition size. Returns "no estimate" (falls +back to the uncorrected formula, warns) when fewer than 2 partitions have +data, `\hat\mu \le 0`, or `V` doesn't exceed the Poisson floor +`\hat\mu/\bar L` — no detectable over-dispersion, `α` would be unbounded. + +**Caveat, stated explicitly rather than left implicit**: unlike every +closed-form correction in `snp_distance.rs` (each verified line-by-line +against `ape`'s `dist_dna.c`, and now also against Jin & Nei 1990 directly +for the base `+Γ` formula), this estimator is derived from first +principles (the general Poisson-Gamma/negative-binomial identity) with no +primary-source procedure behind it at all — not Jin & Nei's (confirmed +above), and Wilson et al. (1989), the paper they point to instead, hasn't +been read/verified either. Mathematically self-consistent (re-derived and +checked, not guessed), but a from-scratch method, not a literature +implementation. If `--gamma-shape` needs a value with a literature +pedigree rather than an estimated one, Jin & Nei's own stated default, +`α = 1` (`--gamma-shape 1`), is the better-supported choice today. + +Deliberately **not** gated by `--sankoff-ratio-ceiling` the way +`base_pair_tally` is (same precedent as `cardinality_tally` — see its own +doc comment): that filter excludes individual saturated *pairs* from a +composition estimate computed once at the very end, from the complete +`PairwiseTally`; the partition axis needed here only exists transiently, +one partition at a time, while `PairwiseTally` is still being built — long +before any pair's final SNP ratio (and thus its ratio_ceiling eligibility) +is known. `--exclude-genome` isn't applied either, matching +`reduce_pairwise`'s own raw per-pair fold. **Implemented now: `snp-raw`, `snp-jc`, `snp-k2p`, `snp-k81`, `snp-f81`, `snp-t92`, `snp-tn93`, `snp-tv`, all with `+Γ` except `raw`/`tv`** — see diff --git a/UserDocMD/usage/phylo.md b/UserDocMD/usage/phylo.md index f933ddd4..dd2ed570 100644 --- a/UserDocMD/usage/phylo.md +++ b/UserDocMD/usage/phylo.md @@ -17,7 +17,7 @@ obikmer phylo INDEX [OPTIONS] | Option | Default | Description | |---|---|---| | `--distance` | `jaccard` | See the two tables below for the full list of accepted values | -| `--gamma-shape ALPHA` | none | Rate-heterogeneity correction, for `snp-*` values that support it (see below). No effect on the other values; rejected if given together with a value that doesn't support it | +| `--gamma-shape ALPHA\|auto` | none | Rate-heterogeneity correction, for `snp-*` values that support it (see below). Either a fixed $\alpha$ or `auto`/`estimate` to fit it from the data (see "Automatic $\alpha$ estimation" below). No effect on the other values; rejected if given together with a value that doesn't support it | | `--presence-threshold` | `1` | Minimum count for a kmer to be considered present, for `jaccard`/`mash` on a count index | | `--csv` | off | Write the matrix as plain CSV instead of the default relaxed-PHYLIP format | | `--shared-kmers` | off | Also write the shared-kmer count matrix. Only valid with a whole-index metric, not a `snp-*` value | @@ -121,6 +121,28 @@ $$d = Q$$ `--gamma-shape ALPHA` applies to every value above except `snp-raw` and `snp-tv`: each $-\ln(x)$ term in the formulas above is replaced by $\alpha\left(x^{-1/\alpha}-1\right)$ (the same weight, same $x$). +### Automatic $\alpha$ estimation (`--gamma-shape auto`) + +`--gamma-shape auto` (or the equivalent `--gamma-shape estimate`) fits $\alpha$ from the index itself instead of requiring a user-supplied value, using a method-of-moments estimator computed once, from the same sampling pass that builds the pairwise substitution tally — no extra scan of the index. + +The estimator pools substitution counts by **partition** rather than by genome pair: for partition $i$, let $n_i$ be the total number of substitutions observed across every genome pair, and $L_i$ the total number of eligible loci across every genome pair, in that partition. Define the partition's observed substitution rate: + +$$R_i = \frac{n_i}{L_i}$$ + +Under a single shared substitution rate with no among-site heterogeneity, each $R_i$ would vary only by Poisson sampling noise. Rate heterogeneity is modeled, as elsewhere in this correction, by a $\mathrm{Gamma}(\alpha,\alpha)$-distributed multiplicative rate (mean 1) shared by every locus in a partition — the classical Poisson–Gamma (negative-binomial) mixture. Under that model: + +$$\mathbb{E}[R_i] = \mu \qquad \mathrm{Var}[R_i] = \frac{\mu}{L_i} + \frac{\mu^2}{\alpha}$$ + +where $\mu$ is the pooled substitution rate across every partition. Weighting each partition's squared deviation by its own $L_i$ removes the first (Poisson) term before attributing what's left to genuine rate heterogeneity: + +$$\hat\mu = \frac{\sum_i n_i}{\sum_i L_i} \qquad V = \frac{\sum_i L_i\,(R_i-\hat\mu)^2}{\sum_i L_i} \qquad \bar L = \frac{\sum_i L_i}{\text{number of partitions}}$$ + +$$\hat\alpha = \frac{\hat\mu^2}{V - \hat\mu/\bar L}$$ + +If the measured variance $V$ doesn't exceed the Poisson floor $\hat\mu/\bar L$ (no detectable over-dispersion across partitions — the data are consistent with a single shared rate), $\alpha$ is left undefined: the correction is silently disabled for that run rather than applying a fabricated value, and a warning is logged. When an estimate is produced, it's logged at the `info` level before the distance matrix is computed. + +Note: this is a method-of-moments estimator derived from the standard Poisson–Gamma relationship between substitution counts and gamma-distributed rate variation, applied per-partition — it is not part of Jin & Nei's (1990) original publication, which only defines the `+Γ` distance formula itself and, absent an estimate, recommends the fixed default $\alpha = 1$ (`--gamma-shape 1`) rather than proposing a way to estimate it from data. `alpha < 1` indicates strong among-site rate heterogeneity (many near-invariant loci, a few fast ones); `alpha` growing large makes the correction converge to the uncorrected formula. + ### Output Without `-o`, the matrix goes to stdout in relaxed-PHYLIP format (`n` on the first line, then one `labelvalue...` row per genome). With `--csv`, the format is instead a header row `genome,,,...` followed by one `