feat(analysis): add two-part count models and occupancy fitter

EXPLAINME.adoc

@@ -0,0 +1,52 @@
+= Occupancy and two-part abundance models: claim receipts
+:toc:
+
+This page records what is implemented in the first slice of issue #18 and what
+is intentionally *not* claimed. The remaining end-user occupancy workflow is
+tracked in follow-up issue #205; code present today is not to be confused with
+that outstanding work.
+
+== Implemented in the current PR
+
+[cols="1,2,2",options="header"]
+|===
+|Claim |Implementation receipt |Boundary / evidence
+
+|Selectable ZINB differential model
+|`src/analysis/differential.jl` (`fit_two_part_nb`, `differential_abundance`), config enum in `src/server/routes/config.jl` and `frontend/src/components/PipelineStages.tsx`
+|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.
+
+|Selectable hurdle negative-binomial differential model
+|Same fitter and config/UI path; the positive count likelihood is zero-truncated NB.
+|All observed zeros are assigned to the occurrence process. It does not offer a lognormal or Poisson count component yet.
+
+|Model cautions and results
+|`frontend/src/components/DifferentialPanel.tsx`, CSV result types/export, diagnostics in `differential_abundance`
+|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.
+
+|Repeated-visit occupancy estimator
+|`MetaManifold.Differential.fit_occupancy`; tests in `test/unit/test_differential.jl`
+|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.
+
+|Scientific limitations
+|`docs/occupancy-models.md`, `README.md`
+|No single-visit occupancy claim. Closure and destructive microbiome sampling limitations are stated. Detection covariates are not implemented.
+|===
+
+== Not yet implemented
+
+The follow-up issue linked from #18 owns first-class occupancy configuration and
+UI, occupancy/detection formula inputs and covariates, identifiability checks,
+zero-formula/count-distribution options, model diagnostics/comparisons, reference
+parity, scale benchmarks, schema/provenance completion, and a synchronized wiki
+roadmap. Until that work lands, do not describe the occupancy estimator as a
+ready-to-use workflow.
+
+== Validation record
+
+* Frontend typecheck, production build and frontend unit tests are run locally.
+* Julia/R backend tests are defined in `test/unit/test_differential.jl` and run
+  by the repository CI workflow. Their successful CI result is required before
+  merge; this sandbox does not have Julia or R installed.
+* CI also contains other workflows; each current failure must be inspected and
+  either fixed or shown to be an unrelated, pre-existing non-required lane.

README.md

@@ -20,7 +20,7 @@ MetaManifold wraps standard amplicon sequencing workflows into a single configur
 
 ## Overview
 
-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)).
+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).
 
 **Pipeline stages**
 
@@ -289,48 +289,38 @@ not exclude samples from this analysis.
 
 `analysis.differential.method` chooses the method:
 
-- `nb_glm` (default): a negative-binomial GLM per taxon with a library-size
-  offset. The effect is log2 fold change; positive values indicate greater
-  abundance in the second condition. `offset` selects `tss` (default) or `rle`
-  scaling. `dispersion` selects per-taxon maximum likelihood with
-  `MASS::glm.nb` (default) or a cross-taxon shrinkage trend from R's
-  `glmGamPoi::glm_gp`, followed by fixed-theta NB Wald fits.
+- `nb_glm` (default): a negative-binomial GLM per taxon using `MASS::glm.nb`,
+  with a library-size offset. The effect is log2 fold change; positive values
+  indicate greater abundance in the second condition. `offset` selects `tss`
+  (default) or `rle`; `dispersion` selects maximum-likelihood theta or a
+  cross-taxon `glmGamPoi` trend.
+- `zinb`: a zero-inflated negative-binomial model with separate logistic
+  zero and count components, a count-component library-size offset, a joint
+  2-df likelihood-ratio test, and BH correction.
+- `hurdle_nb`: a logistic occurrence model plus zero-truncated NB counts for
+  positive observations. It uses the same joint test and offset convention.
 - `clr_welch`: Welch's t-test on centred log-ratios after zero replacement with
   `zCompositions::cmultRepl`. The effect is a CLR difference, not a fold change.
-  `replacement_method` selects `bayesian_multiplicative` (GBM, the default) or
-  `multiplicative_replacement` (CZM). `replacement_delta` controls the
-  replacement fraction (default 0.65), and optional positive `bayesian_alpha`
-  overrides GBM's estimated Dirichlet concentration; `offset` is not used.
-  [Zero handling, posterior semantics, warnings and references](docs/statistics/zero-handling.md)
-  explain why every replacement is biased and when a count, occupancy or
-  zero-inflated model is preferable.
-- `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; recommended when
-  asymmetric shifts or microbial load changes would bias standard CLR centering
-  (requires 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`)
-  posterior instances per sample and reports standardized CLR effect sizes with
-  consistent-sign posterior predictive Welch's t-test p-values; recommended for
-  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 biological overdispersion;
-  completely separated single-group taxa are marked `failed` rather than
-  diverging.
-
-For every method, `min_prevalence` (default 0) filters taxa by the fraction of
-samples with nonzero counts, and Benjamini–Hochberg adjustment is mandatory
-across the fitted taxa that produced a p-value. Failed or filtered taxa have no
-substitute p-value; inspect their status and note in the results table. Invalid
-input or an unavailable required R package produces an explicit error. See
-`bench/compositional_differential/README.md` for method-selection guidance,
-benchmarks, and literature citations.
+  `replacement_method` selects Bayesian-multiplicative (GBM, default) or
+  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).
+
+ZINB and hurdle models can overfit sparse data; do not select ZINB solely
+because counts are sparse. Compare against NB diagnostics. The methods report
+a count-process effect and the zero-process coefficient; this is not proof of
+biological structural absence. The public
+`MetaManifold.Differential.fit_occupancy` lower-level estimator accepts
+repeated-visit detection histories (sites × visits × taxa), refuses single-visit
+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).
+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;
+inspect their status and note in the results table. Invalid input or an
+unavailable required R package produces an explicit error.
 
 Implementation and reference tests: `src/analysis/differential.jl`,
 `src/server/routes/analysis.jl`, `bench/compositional_differential/benchmark.jl`,

ROADMAP.md

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

config/defaults/pipeline.yml

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

docs/occupancy-models.md

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