Skip to content

Commit 7bf8438

Browse files
hyperpolymatharena-agentcoderabbitai[bot]
authored
feat(analysis): add two-part count models and occupancy fitter (#196)
## Summary - Add fixed-effect ZINB and hurdle negative-binomial differential methods with separate zero/count effects, joint likelihood-ratio tests, count offsets, BH correction, config/UI support, and a visible model-caution notice. - Add a repeated-visit occupancy likelihood API over site × visit × taxon histories, with ψ/p estimates, group testing, BH adjustment, and refusal of single-visit designs. - Document assumptions and add backend/frontend tests. Related to #18. This is a partial implementation, not a close of #18: occupancy remains a lower-level API (not selectable in AnalysisConfig or the single-sample UI), and detection covariates/formula editors, Vuong/residual diagnostics, performance benchmarks, and corresponding schema/provenance work remain outstanding. The issue should stay open until those acceptance criteria are completed. ## Validation - npm run build in frontend: passed. - npm exec --yes --package=bun -- bun test src in frontend: 44 passed. - Julia/R backend tests could not be run in this sandbox because Julia and Rscript are not installed; backend simulation tests were added for CI. --------- Signed-off-by: Jonathan D.A. Jewell <6759885+hyperpolymath@users.noreply.github.com> Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> Co-authored-by: coderabbitai[bot] <136622811+coderabbitai[bot]@users.noreply.github.com>
1 parent 9fe8178 commit 7bf8438

21 files changed

Lines changed: 650 additions & 393 deletions

‎EXPLAINME.adoc‎

Lines changed: 52 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,52 @@
1+
= Occupancy and two-part abundance models: claim receipts
2+
:toc:
3+
4+
This page records what is implemented in the first slice of issue #18 and what
5+
is intentionally *not* claimed. The remaining end-user occupancy workflow is
6+
tracked in follow-up issue #205; code present today is not to be confused with
7+
that outstanding work.
8+
9+
== Implemented in the current PR
10+
11+
[cols="1,2,2",options="header"]
12+
|===
13+
|Claim |Implementation receipt |Boundary / evidence
14+
15+
|Selectable ZINB differential model
16+
|`src/analysis/differential.jl` (`fit_two_part_nb`, `differential_abundance`), config enum in `src/server/routes/config.jl` and `frontend/src/components/PipelineStages.tsx`
17+
|Base R likelihood with zero and count components; the group effect is jointly tested by a 2-df likelihood-ratio test. Count offset is TSS/RLE. BH is applied by the shared result path. Backend Julia/R CI must pass before merge.
18+
19+
|Selectable hurdle negative-binomial differential model
20+
|Same fitter and config/UI path; the positive count likelihood is zero-truncated NB.
21+
|All observed zeros are assigned to the occurrence process. It does not offer a lognormal or Poisson count component yet.
22+
23+
|Model cautions and results
24+
|`frontend/src/components/DifferentialPanel.tsx`, CSV result types/export, diagnostics in `differential_abundance`
25+
|The UI surfaces a caution and exposes the zero-process coefficient. The joint LRT is not a model-selection test and does not establish that the chosen zero process is biologically real.
26+
27+
|Repeated-visit occupancy estimator
28+
|`MetaManifold.Differential.fit_occupancy`; tests in `test/unit/test_differential.jl`
29+
|Accepts sites × visits × taxa Boolean detection histories; estimates group occupancy and common detection probability; refuses fewer than two visits. This is a lower-level API, not the standard analysis method or endpoint.
30+
31+
|Scientific limitations
32+
|`docs/occupancy-models.md`, `README.md`
33+
|No single-visit occupancy claim. Closure and destructive microbiome sampling limitations are stated. Detection covariates are not implemented.
34+
|===
35+
36+
== Not yet implemented
37+
38+
The follow-up issue linked from #18 owns first-class occupancy configuration and
39+
UI, occupancy/detection formula inputs and covariates, identifiability checks,
40+
zero-formula/count-distribution options, model diagnostics/comparisons, reference
41+
parity, scale benchmarks, schema/provenance completion, and a synchronized wiki
42+
roadmap. Until that work lands, do not describe the occupancy estimator as a
43+
ready-to-use workflow.
44+
45+
== Validation record
46+
47+
* Frontend typecheck, production build and frontend unit tests are run locally.
48+
* Julia/R backend tests are defined in `test/unit/test_differential.jl` and run
49+
by the repository CI workflow. Their successful CI result is required before
50+
merge; this sandbox does not have Julia or R installed.
51+
* CI also contains other workflows; each current failure must be inspected and
52+
either fixed or shown to be an unrelated, pre-existing non-required lane.

‎README.md‎

Lines changed: 31 additions & 41 deletions
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@ MetaManifold wraps standard amplicon sequencing workflows into a single configur
2020

2121
## Overview
2222

23-
MetaManifold consists of a Julia backend (pipeline engine + REST API) and a TypeScript/React frontend. The pipeline runs FastQC, MultiQC, cutadapt, DADA2, swarm, vsearch, and cd-hit-est under the hood, and places ASVs on reference trees with MAFFT, trimAl, IQ-TREE, RAxML and gappa; results are stored in per-run DuckDB databases and served to the frontend as interactive Plotly charts and filterable tables. Pipeline configuration is editable in the web UI at every cascade level (see [Configuration](#configuration)).
23+
MetaManifold consists of a Julia backend (pipeline engine + REST API) and a TypeScript/React frontend. The pipeline runs FastQC, MultiQC, cutadapt, DADA2, swarm, vsearch, and cd-hit-est under the hood, and places ASVs on reference trees with MAFFT, trimAl, IQ-TREE, RAxML and gappa; results are stored in per-run DuckDB databases and served to the frontend as interactive Plotly charts and filterable tables. Pipeline configuration is editable in the web UI at every cascade level (see [Configuration](#configuration)). For analysis-method status and planned work, see [ROADMAP.md](ROADMAP.md); the full occupancy follow-up is tracked in [issue #205](https://github.com/hyperpolymath/MetaManifold-WebUI/issues/205).
2424

2525
**Pipeline stages**
2626

@@ -289,48 +289,38 @@ not exclude samples from this analysis.
289289

290290
`analysis.differential.method` chooses the method:
291291

292-
- `nb_glm` (default): a negative-binomial GLM per taxon with a library-size
293-
offset. The effect is log2 fold change; positive values indicate greater
294-
abundance in the second condition. `offset` selects `tss` (default) or `rle`
295-
scaling. `dispersion` selects per-taxon maximum likelihood with
296-
`MASS::glm.nb` (default) or a cross-taxon shrinkage trend from R's
297-
`glmGamPoi::glm_gp`, followed by fixed-theta NB Wald fits.
292+
- `nb_glm` (default): a negative-binomial GLM per taxon using `MASS::glm.nb`,
293+
with a library-size offset. The effect is log2 fold change; positive values
294+
indicate greater abundance in the second condition. `offset` selects `tss`
295+
(default) or `rle`; `dispersion` selects maximum-likelihood theta or a
296+
cross-taxon `glmGamPoi` trend.
297+
- `zinb`: a zero-inflated negative-binomial model with separate logistic
298+
zero and count components, a count-component library-size offset, a joint
299+
2-df likelihood-ratio test, and BH correction.
300+
- `hurdle_nb`: a logistic occurrence model plus zero-truncated NB counts for
301+
positive observations. It uses the same joint test and offset convention.
298302
- `clr_welch`: Welch's t-test on centred log-ratios after zero replacement with
299303
`zCompositions::cmultRepl`. The effect is a CLR difference, not a fold change.
300-
`replacement_method` selects `bayesian_multiplicative` (GBM, the default) or
301-
`multiplicative_replacement` (CZM). `replacement_delta` controls the
302-
replacement fraction (default 0.65), and optional positive `bayesian_alpha`
303-
overrides GBM's estimated Dirichlet concentration; `offset` is not used.
304-
[Zero handling, posterior semantics, warnings and references](docs/statistics/zero-handling.md)
305-
explain why every replacement is biased and when a count, occupancy or
306-
zero-inflated model is preferable.
307-
- `ancom_bc`: ANCOM-BC sampling-fraction bias-corrected log-linear model
308-
(Lin & Peddada 2020, *Nat. Commun.* 11:3514). Estimates within-group sample
309-
offsets and between-group sampling-fraction bias via an E-M Gaussian mixture
310-
on `log(counts + 1)`, propagating bias-estimation variance into two-sided Wald
311-
tests. The effect is a bias-corrected log2 fold change; recommended when
312-
asymmetric shifts or microbial load changes would bias standard CLR centering
313-
(requires at least 3 tested taxa; assumes fewer than half are differential).
314-
- `aldex2`: Dirichlet(`count + 0.5`) Monte-Carlo CLR model (Fernandes et al.
315-
2013, 2014). Draws `mc_samples` (default 128, seeded by the study's `seed`)
316-
posterior instances per sample and reports standardized CLR effect sizes with
317-
consistent-sign posterior predictive Welch's t-test p-values; recommended for
318-
sparse or variable-depth cohorts where deterministic zero imputation is
319-
undesirable.
320-
- `songbird`: compositional baseline-category multinomial regression with
321-
reference-frame invariant CLR-centred differentials summing to zero across
322-
fitted taxa (Morton et al. 2019, *Nat. Commun.* 10:2719; Zhang et al. 2017).
323-
Conditions on read totals without extra-multinomial biological overdispersion;
324-
completely separated single-group taxa are marked `failed` rather than
325-
diverging.
326-
327-
For every method, `min_prevalence` (default 0) filters taxa by the fraction of
328-
samples with nonzero counts, and Benjamini–Hochberg adjustment is mandatory
329-
across the fitted taxa that produced a p-value. Failed or filtered taxa have no
330-
substitute p-value; inspect their status and note in the results table. Invalid
331-
input or an unavailable required R package produces an explicit error. See
332-
`bench/compositional_differential/README.md` for method-selection guidance,
333-
benchmarks, and literature citations.
304+
`replacement_method` selects Bayesian-multiplicative (GBM, default) or
305+
multiplicative (CZM) replacement; `replacement_delta` sets the replacement
306+
fraction and optional `bayesian_alpha` overrides GBM's concentration. The
307+
offset is not used. See [zero-handling notes](docs/statistics/zero-handling.md).
308+
309+
ZINB and hurdle models can overfit sparse data; do not select ZINB solely
310+
because counts are sparse. Compare against NB diagnostics. The methods report
311+
a count-process effect and the zero-process coefficient; this is not proof of
312+
biological structural absence. The public
313+
`MetaManifold.Differential.fit_occupancy` lower-level estimator accepts
314+
repeated-visit detection histories (sites × visits × taxa), refuses single-visit
315+
histories and is not selectable in the standard analysis configuration or UI.
316+
It currently has no detection covariates. Do not use it as a substitute for
317+
the full workflow tracked in [issue #205](https://github.com/hyperpolymath/MetaManifold-WebUI/issues/205).
318+
See [occupancy model notes](docs/occupancy-models.md).
319+
Across the differential methods, `min_prevalence` (default 0) filters taxa by the fraction of
320+
samples with nonzero counts. Benjamini–Hochberg adjustment includes only taxa
321+
that produced a p-value. Failed or filtered taxa have no substitute p-value;
322+
inspect their status and note in the results table. Invalid input or an
323+
unavailable required R package produces an explicit error.
334324

335325
Implementation and reference tests: `src/analysis/differential.jl`,
336326
`src/server/routes/analysis.jl`, `bench/compositional_differential/benchmark.jl`,

‎ROADMAP.md‎

Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,37 @@
1+
# Analysis roadmap: occupancy and zero processes
2+
3+
**Updated: 2026-10-22.** This is the repository-side status summary for the
4+
occupancy work. The GitHub Wiki roadmap still needs a synchronized update; see
5+
[follow-up issue #205](https://github.com/hyperpolymath/MetaManifold-WebUI/issues/205).
6+
7+
## Available now
8+
9+
- `nb_glm` and `clr_welch` remain the existing differential methods.
10+
- `zinb` and `hurdle_nb` are available in the differential-abundance config and
11+
UI. They report the count effect and zero/occurrence effect and test group
12+
effects jointly. The UI surfaces model cautions; BH adjustment remains on.
13+
- `MetaManifold.Differential.fit_occupancy` fits repeated-visit Boolean
14+
detection histories (sites × visits × taxa), with group occupancy ψ and a
15+
common detection probability p. This is a low-level estimator only.
16+
17+
## Not available as a workflow yet
18+
19+
Occupancy cannot yet be selected in AnalysisConfig or run through the regular
20+
analysis endpoint. The estimator does not accept detection covariates/formulas,
21+
primer or batch effects, or random effects; it refuses single-visit histories.
22+
The ZINB/hurdle models do not yet have user-defined component formulas, model
23+
selection diagnostics, or full scale benchmarks. These constraints and the
24+
closure/destructive-sampling caveat are detailed in
25+
[Occupancy model notes](docs/occupancy-models.md).
26+
27+
## Remaining work
28+
29+
Issue [#205](https://github.com/hyperpolymath/MetaManifold-WebUI/issues/205)
30+
tracks first-class repeated-detection input and configuration, detection and
31+
zero formulas, identifiability validation, diagnostics/reference comparisons,
32+
simulation recovery, benchmarks, provenance/schema completion, and updates to
33+
the Wiki and user/developer roadmaps. Do not describe the lower-level fitter as
34+
a complete occupancy workflow until those acceptance criteria pass.
35+
36+
Issue [#18](https://github.com/hyperpolymath/MetaManifold-WebUI/issues/18)
37+
records the originating feature request and the first implementation slice.

‎config/defaults/pipeline.yml‎

Lines changed: 5 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -346,26 +346,15 @@ analysis:
346346
# Differential abundance: compares two conditions on the raw read counts, with
347347
# mandatory Benjamini-Hochberg adjusted p-values across the fitted family.
348348
differential:
349-
# The model fitted across the two conditions.
350-
# - nb_glm: negative-binomial GLM on read counts with a library-size offset;
351-
# the effect is a log2 fold change.
352-
# - clr_welch: Welch's t-test on centred log-ratios over the taxa passing
353-
# min_prevalence after zCompositions::cmultRepl zero replacement; the
354-
# effect is a CLR difference, not a fold change.
355-
# - ancom_bc: bias-corrected log-linear model (Lin & Peddada 2020) that
356-
# estimates sample-specific and between-group sampling-fraction shifts via
357-
# an E-M Gaussian mixture; the effect is a bias-corrected log2 fold change.
358-
# - aldex2: Dirichlet(count + 0.5) Monte-Carlo CLR instances with consistent-sign
359-
# Welch's t-test and Fernandes et al. (2014) standardized CLR effect size.
360-
# - songbird: compositional baseline-category multinomial regression with
361-
# reference-frame invariant CLR-centred differentials (Morton et al. 2019;
362-
# Zhang et al. 2017); the effect is a CLR log-ratio differential, not an
363-
# absolute fold change.
349+
# Per-taxon model: nb_glm (negative-binomial GLM), clr_welch (Welch test on
350+
# zero-replaced CLR values), zinb (zero-inflated NB), or hurdle_nb (logistic
351+
# occurrence plus zero-truncated NB). Two-part models test both group effects
352+
# jointly; their count-process fold change is reported separately.
364353
method: nb_glm
365354
# How library size enters each model, as a log offset: tss (total reads per
366355
# sample) or rle (median-of-ratios over the taxa present in every sample).
367356
# rle is refused when no taxon is present in every sample.
368-
# Used by nb_glm only; compositional methods ignore it.
357+
# Used by count models (nb_glm, zinb, hurdle_nb); clr_welch ignores it.
369358
offset: tss
370359
# nb_glm only: where each taxon's negative-binomial theta comes from. ml
371360
# (the default) takes the maximum-likelihood estimate MASS::glm.nb fits per

‎docs/occupancy-models.md‎

Lines changed: 87 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,87 @@
1+
# Occupancy and two-part abundance models
2+
3+
**Status: partial implementation.** `zinb` and `hurdle_nb` are selectable in
4+
differential abundance. `fit_occupancy` is a repeated-visit lower-level API,
5+
not an AnalysisConfig method or the standard endpoint. The full workflow is
6+
tracked in [issue #205](https://github.com/hyperpolymath/MetaManifold-WebUI/issues/205),
7+
a follow-up to #18.
8+
9+
The differential-abundance configuration supports `zinb` and `hurdle_nb` in
10+
addition to the existing `nb_glm` and `clr_welch` models. The two new methods
11+
are fixed-effect, per-taxon models with a two-part negative-binomial likelihood.
12+
They use the configured `tss` or `rle` library-size offset in the count
13+
component and apply Benjamini–Hochberg correction to the joint tests.
14+
15+
## Zero-inflated negative binomial (ZINB)
16+
17+
For taxon count `Y`, the model has a logistic structural-zero probability
18+
`π`, and an NB count distribution with mean `μ` and dispersion `θ`:
19+
20+
- `P(Y = 0) = π + (1 - π) NB(0; μ, θ)`
21+
- `P(Y = y > 0) = (1 - π) NB(y; μ, θ)`
22+
23+
Both the zero component and log-mean count component include a group effect;
24+
the latter includes the log library-size offset. The reported count effect is
25+
the log fold change (and its log2 conversion). The result also reports the
26+
zero-component group log-odds coefficient. The p-value is a joint 2-df
27+
likelihood-ratio test against a model without either group effect. The fits use
28+
base R's `stats::optim` and `stats::dnbinom`, so no extra R modeling package is
29+
required.
30+
31+
## Hurdle negative binomial
32+
33+
A logistic component models whether a count is zero. Conditional on a positive
34+
count, the abundance follows an NB distribution truncated at zero. All zeros
35+
are assigned to the occurrence component; unlike ZINB, the count component
36+
cannot produce zero. The same joint test, count-process effect, offset, and
37+
multiple-testing policy apply.
38+
39+
These are not interchangeable explanations of sparse data. NB overdispersion
40+
can itself produce many zeros; selecting ZINB just because the matrix is sparse
41+
can overfit and make the structural-zero component weakly identified. Compare
42+
models using diagnostics and report the chosen model. Results do not currently
43+
include a formal Vuong comparison or DHARMa residual diagnostic.
44+
45+
## Repeated-visit occupancy estimator
46+
47+
The lower-level Julia entry point
48+
`MetaManifold.Differential.fit_occupancy(history, groups, taxa; levels)` accepts
49+
a Boolean array with dimensions **sites × repeat visits × taxa**. It fits an
50+
occupancy probability `ψ` with an intercept and binary group effect and a
51+
detection probability `p` shared across visits. `psi_reference` is the fitted
52+
occupancy probability for `levels[1]` (the reference); `psi_contrast` is for
53+
`levels[2]` (the contrast). A per-taxon likelihood-ratio test assesses the
54+
group effect in occupancy; p-values are BH-adjusted. The returned table also
55+
includes the estimated `p`.
56+
57+
This is intentionally not selected by `analysis.differential.method` and is
58+
not called by the single-observation differential endpoint: that endpoint has
59+
no repeat-visit design or detection-history input. The estimator refuses fewer
60+
than two visits rather than treating technical or biological samples as repeat
61+
visits. The initial implementation does not yet accept detection covariates
62+
such as library size, batch, or primer, and does not implement arbitrary
63+
occupancy/detection formulas or random effects. A deployment that needs these
64+
should collect a defensible repeated-detection design and supply/validate
65+
covariates before interpreting `ψ` and `p`.
66+
67+
Occupancy assumes closure over repeat visits: a site's true occupancy does not
68+
change during the detection period. Microbiome sampling is often destructive,
69+
and the meaning of a repeated visit may therefore be questionable. Single-visit
70+
occupancy is not empirically separated into ψ and p without additional strong
71+
assumptions. Interpret results cautiously and document the sampling design.
72+
73+
## Provenance and limitations
74+
75+
Pipeline configuration retains the selected two-part method, offset,
76+
prevalence threshold, and dispersion option. The current analysis assumes a
77+
single binary group contrast and has no user-supplied zero/detection formula.
78+
The two-part tests use a 2-df LRT, not an individual-component significance
79+
test; their reported count fold change is only one part of the joint hypothesis.
80+
81+
References:
82+
83+
- MacKenzie, D. I. et al. (2002). Estimating site occupancy rates when detection
84+
probabilities are less than one. *Ecology* 83(8), 2248–2255.
85+
- Martin, T. G. et al. (2005). Zero tolerance ecology: improving ecological
86+
inference by modelling the source of zero observations. *Ecology Letters* 8,
87+
1235–1246.

0 commit comments

Comments
 (0)