Skip to content

fix: load the package, and restore the docs and test registrations main lost #209

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
hyperpolymath merged 1 commit into from
Oct 23, 2026
Merged

fix: load the package, and restore the docs and test registrations main lost #209

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
160 changes: 159 additions & 1 deletion README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -305,6 +305,33 @@ not exclude samples from this analysis.
multiplicative (CZM) replacement; `replacement_delta` sets the replacement
fraction and optional `bayesian_alpha` overrides GBM's concentration. The
offset is not used. See [zero-handling notes](docs/statistics/zero-handling.md).
- `ancom_bc`: ANCOM-BC sampling-fraction bias-corrected log-linear model
(Lin & Peddada 2020, *Nat. Commun.* 11:3514). Estimates within-group sample
offsets and between-group sampling-fraction bias via an E-M Gaussian mixture
on `log(counts + 1)`, propagating bias-estimation variance into two-sided Wald
tests. The effect is a bias-corrected log2 fold change; useful when asymmetric
shifts or load changes would bias standard CLR centering (needs at least 3
tested taxa; assumes fewer than half are differential).
- `aldex2`: Dirichlet(`count + 0.5`) Monte-Carlo CLR model (Fernandes et al.
2013, 2014). Draws `mc_samples` (default 128, seeded by the study's `seed`)
instances per sample and reports standardized CLR effect sizes with
consistent-sign posterior predictive Welch p-values; suited to sparse or
variable-depth cohorts where deterministic zero imputation is undesirable.
- `songbird`: compositional baseline-category multinomial regression with
reference-frame invariant CLR-centred differentials summing to zero across
fitted taxa (Morton et al. 2019, *Nat. Commun.* 10:2719; Zhang et al. 2017).
Conditions on read totals without extra-multinomial overdispersion; taxa
separated in one group are marked `failed` rather than diverging.
- `multinomial`: one multinomial logistic model for the whole composition
(Songbird's model in pure Julia), which needs no size factor and no zero
replacement. The effect is a log2 fold change in relative abundance against
`mn_reference_taxon`; `mn_penalty` adds an L2 ridge penalty, and a taxon with
no reads in one group is refused rather than reported.
- `dirichlet_multinomial`: the same joint model plus a Dirichlet prior on the
composition (La Rosa et al. 2012), which absorbs between-sample
overdispersion; `dm_overdispersion_method` picks `mle` (profiled) or `mom`.
Both multinomial models bound their Newton fit with `max_iter`, `tolerance`
and `max_seconds`, and a fit that does not converge is refused.

ZINB and hurdle models can overfit sparse data; do not select ZINB solely
because counts are sparse. Compare against NB diagnostics. The methods report
Expand All @@ -316,6 +343,14 @@ histories and is not selectable in the standard analysis configuration or UI.
It currently has no detection covariates. Do not use it as a substitute for
the full workflow tracked in [issue #205](https://github.com/hyperpolymath/MetaManifold-WebUI/issues/205).
See [occupancy model notes](docs/occupancy-models.md).

The multinomial models are the count-based, compositionally coherent
alternative to `nb_glm` (which fits each taxon independently and so ignores the
constraint the composition imposes) and to `clr_welch` (which is Gaussian on
log-ratios and ignores the count likelihood). They are fitted in Julia without R
and without a pseudocount; the dense observed information limits a fit to 1000
taxa, above which raise `min_prevalence` or analyse a coarser rank.

Across the differential methods, `min_prevalence` (default 0) filters taxa by the fraction of
samples with nonzero counts. Benjamini–Hochberg adjustment includes only taxa
that produced a p-value. Failed or filtered taxa have no substitute p-value;
Expand All @@ -326,6 +361,119 @@ Implementation and reference tests: `src/analysis/differential.jl`,
`src/server/routes/analysis.jl`, `bench/compositional_differential/benchmark.jl`,
and `test/unit/test_differential.jl`.

#### Why a count model rather than a t-test

Sequencing measures a sample of the library, not an absolute quantity: a
sample of twice the depth has twice every count, and one taxon can only increase
at another's expense. Comparisons of raw or rarefied abundances therefore carry
the compositionality problem (Gloor et al. 2017). The multinomial models fit the
whole composition at once and condition on each sample's total, so each
coefficient is the log ratio of two taxa. That ratio is unchanged when every
count in a sample is multiplied by a constant, which is why no size factor and no
pseudocount is needed. The intercept is the log-odds of the taxon against the
reference at the reference group; the group coefficient is the change in that
log-odds, i.e. a log fold change in relative abundance. Because a log ratio is
the sum of log ratios through any third taxon, every reported effect is coherent
under compositionality and can be re-expressed against any other taxon.

#### Choosing the reference taxon

Every reported effect is anchored on the reference taxon, so the choice matters
for interpretability and for stability, but not for the fit: the likelihood and
the log ratio between two non-reference taxa are the same whichever taxon is the
reference. An empty `mn_reference_taxon` takes the most abundant taxon that
passes the prevalence filter. An explicit one must be a taxon at the analysed
rank with reads in at least half the samples, and reads in both groups; a taxon
seen in fewer samples makes every ratio it anchors unstable, and a taxon absent
from one group has no finite log fold change, so both are refused with the
reason. Taxa with no reads in one group are removed from the composition and
listed as `failed` with their note, rather than fitted: in a joint model a single
diverging coefficient would otherwise pull the whole fit with it.

#### Overdispersion

The multinomial treats every read as independent, so its variance is the sampling
variance and nothing else. Real samples vary more than that: the same community,
sequenced twice, differs by more than resampling reads from one library. The
Dirichlet-multinomial adds that between-sample variation as a concentration `c`,
with `theta = 1/(1+c)` the model's intra-sample correlation; the count variance
becomes `N p (1 - p) (1 + (N - 1) theta)`. When the data are overdispersed the
multinomial's standard errors are too small, and the Dirichlet-multinomial
corrects them; when they are not, `c` runs to its upper bound, the fit collapses
to the multinomial exactly, and the affected rows are marked `boundary` with a
note saying so. The two models are otherwise the same model on the same scale,
so their coefficients are directly comparable.

#### Which model to use

| Data and question | Method |
| --- | --- |
| Counts, two groups, few samples, expected overdispersion | `nb_glm` (per-taxon negative binomial; the default) |
| Counts, many taxa and samples, a joint composition and a single reference | `multinomial` (Songbird's model) |
| As above, with visible between-sample overdispersion | `dirichlet_multinomial` |
| Composition already transformed, or a quick robust rank-based check | `clr_welch` |
| Excess zeros beyond negative-binomial overdispersion, with evidence for them | `zinb` or `hurdle_nb` (two-part; see the caution above) |
| Asymmetric shifts or a community-wide load change | `ancom_bc` (bias-corrected) |
| Sparse or variable-depth cohorts, Monte-Carlo CLR | `aldex2` |
| Many taxa, a joint multinomial fit, but through R and with CLR-centred effects | `songbird` |

`nb_glm` fits a negative-binomial GLM per taxon with a library-size offset: it
tolerates overdispersion through its own dispersion parameter and is the safer
choice with very few samples per group, but each taxon is fitted independently
and the offsets are estimated, not modelled. The multinomial models borrow
strength across taxa by fitting one model, and their effect has a direct
compositional reading, at the cost of a joint fit that a single separated taxon
can compromise (hence the separation check and the optional ridge penalty), and
of a dense information matrix that caps a fit at 1000 taxa. Songbird's L1
sparsity path is not implemented: a lasso estimate has no valid Wald standard
error, so this implements the L2-penalised model and Wald inference instead.

The multinomial models are fitted in Julia and need no R package.
`max_iter`, `tolerance` and `max_seconds` bound one fit; a fit that runs out of
time fails with an explicit error instead of holding the request open. Runtime
and memory at 100, 1000 and 10000 taxa are measured by
`bench/compositional_scaling/scaling.jl`, which writes `results/scaling.csv`.

Implementation and reference tests: `src/analysis/compositional.jl` (the two
models), `src/analysis/differential.jl`, `src/server/routes/analysis.jl`, and
`test/unit/test_differential.jl`. The reference values in the tests come from an
independent implementation and are quoted to 12 digits; the special functions
are cross-checked against mpmath at 40 digits.

### Evidence Mode

The **Evidence** tab turns a run's outputs into a citable epistemic record. A
claim is a sentence a reader could disagree with, a **scope** naming the part of
the data it is about (its *fibre*: the rows whose value at a rank, or whose
category in a set, is one of the named labels), and a **warrant** saying why it
follows. Claims are stored in `projects/{study}/[<group>/]<run>/evidence/claims.yml`
and can be downloaded as that exact YAML.

Three views share one record:

- **Claims** — the epistemic editor (avec fibre). A claim cannot be saved
without a scope and at least one value, and cannot be marked `supported`
without a measurement, citation or inference warrant: a definition alone
supports nothing.
- **Fibre & echo** — the fibre's share of every sample, its largest members, and
the same fibre seen through a second projection: the other classifier at the
same rank, and the parent rank. Agreement is reported with a verdict;
disagreement is a limit on how strongly the claim may be stated, and an echo
that cannot be computed says so rather than vanishing.
- **Residuals** — where the fibre departs from a null you declare, either
`proportional` (the fibre distributed across samples as the library is) or a
stated `claim_share`, with Pearson residuals, a log2 ratio and a χ² statistic
printed without a p-value. The worst departure is shown first.

The tab's Provenance card states whether `attestation.yml` is absent, degraded
or internally inconsistent before any claim is read. Longer notes on the
definitions, the schema and the limits are kept in `docs/evidence-mode.md`
(documentation under `docs/` is kept local by policy — see `.gitignore`).

Implementation: `src/core/evidence.jl` (schema, residuals, verdicts),
`src/server/routes/evidence.jl` (the SQL), `frontend/src/evidence/`
(editor, fibre visualiser, residual explorer) and `test/unit/test_evidence.jl`.

### Figures and report

The Figures tab lays analysis charts out on a page (A4, letter or a journal
Expand Down Expand Up @@ -662,13 +810,23 @@ analysis:
normalisation: hellinger # hellinger | none | rarefy | srs
depth: 0 # rarefy/srs reads per sample; 0 = smallest non-empty sample
differential:
method: nb_glm # nb_glm | clr_welch | ancom_bc | aldex2 | songbird
method: nb_glm # nb_glm | clr_welch | zinb | hurdle_nb | ancom_bc
# | aldex2 | songbird | multinomial | dirichlet_multinomial
offset: tss # count models: tss | rle
dispersion: ml # nb_glm: ml | glmGamPoi
min_prevalence: 0.0 # fraction of samples a taxon must be seen in
replacement_method: bayesian_multiplicative
# clr_welch: bayesian_multiplicative | multiplicative_replacement
replacement_delta: 0.65 # clr_welch: strictly between 0 and 1
bayesian_alpha: ~ # clr_welch: null (GBM estimate) or a number > 0
mc_samples: 128 # aldex2: Dirichlet Monte-Carlo instances in [16, 4096]
mn_reference_taxon: ~ # multinomial models: null picks the most abundant taxon
mn_penalty: 0.0 # multinomial models: L2 ridge penalty, >= 0
dm_overdispersion_method: mle
# dirichlet_multinomial: mle | mom
max_iter: 500 # multinomial models: Newton iteration cap
tolerance: 1.0e-8 # multinomial models: convergence tolerance
max_seconds: 300.0 # multinomial models: wall-clock budget
nmds:
max_stress: 0.2 # warn if NMDS stress exceeds this value
```
Expand Down
Loading
Loading