From 28bd73242dfd385e83c2e332a6f602254b7f36ef Mon Sep 17 00:00:00 2001
From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com>
Date: Thu, 22 Oct 2026 22:17:33 +0000
Subject: [PATCH 1/6] feat(analysis): add two-part count models and occupancy
 fitter

Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com>
---
 README.md                                     |  17 +-
 config/defaults/pipeline.yml                  |  11 +-
 docs/occupancy-models.md                      |  79 ++++++
 frontend/src/api/types.ts                     |  10 +-
 frontend/src/components/DifferentialPanel.tsx |  13 +-
 .../src/components/PipelineStages.test.ts     |   7 +
 frontend/src/components/PipelineStages.tsx    |   6 +-
 .../src/components/differentialTable.test.ts  |   6 +-
 frontend/src/components/differentialTable.ts  |   6 +-
 src/analysis/differential.jl                  | 260 +++++++++++++++++-
 src/server/routes/config.jl                   |   4 +-
 test/unit/test_differential.jl                |  57 ++++
 12 files changed, 447 insertions(+), 29 deletions(-)
 create mode 100644 docs/occupancy-models.md

diff --git a/README.md b/README.md
index 998fcc15..482ee321 100644
--- a/README.md
+++ b/README.md
@@ -261,12 +261,27 @@ not exclude samples from this analysis.
   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.
+- `zinb`: a zero-inflated negative-binomial model with separate logistic
+  zero and count components, a count-component library-size offset, a joint
+  two-degree-of-freedom 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.
+
+ZINB and hurdle models are more flexible than NB and can overfit sparse data;
+prefer them only when diagnostics support excess zeros. They report a count
+process fold change and the zero-process coefficient separately. The public
+`MetaManifold.Differential.fit_occupancy` lower-level estimator is available for
+explicit repeated-visit detection histories (sites × visits × taxa). It refuses
+single-visit histories because occupancy and detection are not separately
+identified from one observation without strong external assumptions. See
+[occupancy model notes](docs/occupancy-models.md).
+
 - `clr_welch`: Welch's t-test on centred log-ratios after multiplicative zero
   replacement with `zCompositions::cmultRepl`. The effect is a CLR difference,
   not a fold change. `replacement_delta` controls the replacement fraction
   (default 0.65); `offset` is not used.
 
-For either method, `min_prevalence` (default 0) filters taxa by the fraction of
+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
diff --git a/config/defaults/pipeline.yml b/config/defaults/pipeline.yml
index 44b2c882..dbac65e9 100644
--- a/config/defaults/pipeline.yml
+++ b/config/defaults/pipeline.yml
@@ -347,16 +347,15 @@ analysis:
   # comparing two conditions on the raw read counts, with Benjamini-Hochberg
   # adjusted p-values.
   differential:
-    # The model fitted to each taxon. nb_glm: a negative-binomial GLM on the 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 replaces their zeros; the effect is a
-    # difference in CLR units, not a 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; clr_welch ignores it, since log-ratios carry no library size.
+    # 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
diff --git a/docs/occupancy-models.md b/docs/occupancy-models.md
new file mode 100644
index 00000000..7e22880c
--- /dev/null
+++ b/docs/occupancy-models.md
@@ -0,0 +1,79 @@
+# Occupancy and two-part abundance models
+
+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. A per-taxon likelihood-ratio
+test assesses the group effect in occupancy; p-values are BH-adjusted. The
+returned table includes `ψ` in each group and 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.
diff --git a/frontend/src/api/types.ts b/frontend/src/api/types.ts
index 4ea20add..f20aaf15 100644
--- a/frontend/src/api/types.ts
+++ b/frontend/src/api/types.ts
@@ -413,10 +413,12 @@ export interface DifferentialRow {
   padj: number | null
   /** nb_glm only. */
   dispersion_theta: number | null
+  /** Group coefficient in the logistic zero/occurrence component (two-part NB only). */
+  zero_process_log_odds?: number | null
   prevalence: number
 }
 
-export type DifferentialMethod = 'nb_glm' | 'clr_welch'
+export type DifferentialMethod = 'nb_glm' | 'clr_welch' | 'zinb' | 'hurdle_nb'
 
 /** Which row field holds the effect size, and what to call it. */
 export interface DifferentialEffect {
@@ -432,12 +434,12 @@ export interface DifferentialResult {
   groups: { reference: string; contrast: string }
   n_samples: Record<string, number>
   effect: DifferentialEffect
-  /** Empty for clr_welch, which takes no library-size offset. */
+  /** Empty for clr_welch; count models report their size factors. */
   size_factors: { sample: string; group: string; factor: number }[]
   config: {
     method: DifferentialMethod
     min_prevalence: number
-    /** nb_glm only. */
+    /** Count-process offset for nb_glm, zinb and hurdle_nb. */
     offset?: 'tss' | 'rle'
     /** nb_glm only: where each taxon's theta comes from. */
     dispersion?: 'ml' | 'glmGamPoi'
@@ -446,6 +448,8 @@ export interface DifferentialResult {
   }
   diagnostics: {
     n_taxa: number; n_tested: number; n_failed: number; n_boundary: number; n_filtered: number
+    /** Safety note for two-part models. */
+    model_warning?: string
     /** clr_welch only. */
     zero_replacement?: {
       method: string; delta: number; zeros_replaced: number; n_taxa_in_composition: number
diff --git a/frontend/src/components/DifferentialPanel.tsx b/frontend/src/components/DifferentialPanel.tsx
index 825fb812..46976592 100644
--- a/frontend/src/components/DifferentialPanel.tsx
+++ b/frontend/src/components/DifferentialPanel.tsx
@@ -138,6 +138,13 @@ export function DifferentialPanel({ study, runs, option, aggregate }: {
 
           <AnalysisChart study={study} figure={result.figure} />
 
+          {result.diagnostics.model_warning && (
+            <div role="note" style={{ border: '1px solid var(--color-danger)', borderLeftWidth: 4,
+              borderRadius: 4, color: 'var(--color-danger)', padding: '8px 12px', margin: '10px 0' }}>
+              <strong>Model caution:</strong> {result.diagnostics.model_warning}
+            </div>
+          )}
+
           <div style={{ display: 'flex', gap: 6, justifyContent: 'flex-end', margin: '8px 0' }}>
             <button className="btn" onClick={() => saveBlob(csvBlob(result), 'differential_abundance.csv')}>
               Download CSV
@@ -152,7 +159,8 @@ export function DifferentialPanel({ study, runs, option, aggregate }: {
                 <tr style={{ borderBottom: '2px solid var(--color-border)' }}>
                   <th style={{ ...head, textAlign: 'left' }}>{result.rank}</th>
                   <th style={head}>{result.effect.label}</th>
-                  {result.config.method === 'nb_glm' && <th style={head}>Estimate (ln)</th>}
+                  {['nb_glm', 'zinb', 'hurdle_nb'].includes(result.config.method) && <th style={head}>Estimate (ln)</th>}
+                  {['zinb', 'hurdle_nb'].includes(result.config.method) && <th style={head}>Zero-process log-odds</th>}
                   <th style={head}>SE</th>
                   <th style={head}>p</th>
                   <th style={head}>padj (BH)</th>
@@ -164,7 +172,8 @@ export function DifferentialPanel({ study, runs, option, aggregate }: {
                   <tr key={r.taxon} style={{ borderBottom: '1px solid var(--color-border)' }}>
                     <td style={{ ...cell, textAlign: 'left' }}>{r.taxon}</td>
                     <td style={cell}>{formatStat(effectValue(result, r))}</td>
-                    {result.config.method === 'nb_glm' && <td style={cell}>{formatStat(r.estimate)}</td>}
+                    {['nb_glm', 'zinb', 'hurdle_nb'].includes(result.config.method) && <td style={cell}>{formatStat(r.estimate)}</td>}
+                    {['zinb', 'hurdle_nb'].includes(result.config.method) && <td style={cell}>{formatStat(r.zero_process_log_odds ?? null)}</td>}
                     <td style={cell}>{formatStat(r.standard_error)}</td>
                     <td style={cell}>{formatStat(r.pvalue)}</td>
                     <td style={cell}>{formatStat(r.padj)}</td>
diff --git a/frontend/src/components/PipelineStages.test.ts b/frontend/src/components/PipelineStages.test.ts
index 37fa1f24..4f1c33a8 100644
--- a/frontend/src/components/PipelineStages.test.ts
+++ b/frontend/src/components/PipelineStages.test.ts
@@ -7,3 +7,10 @@ test('dispersion uses the enum editor with both supported models', () => {
     kind: 'enum', options: ['ml', 'glmGamPoi'],
   })
 })
+
+
+test('differential model editor exposes supported two-part NB methods', () => {
+  expect(CONFIG_TYPES['analysis.differential.method']).toEqual({
+    kind: 'enum', options: ['nb_glm', 'clr_welch', 'zinb', 'hurdle_nb'],
+  })
+})
diff --git a/frontend/src/components/PipelineStages.tsx b/frontend/src/components/PipelineStages.tsx
index f4cc3fbe..ebac68b5 100644
--- a/frontend/src/components/PipelineStages.tsx
+++ b/frontend/src/components/PipelineStages.tsx
@@ -150,8 +150,8 @@ export const CONFIG_DESCRIPTIONS: Record<string, string> = {
   'analysis.alpha.rarefaction_iterations': 'rarefy draws averaged for each alpha diversity value.',
   'analysis.beta.normalisation':  'Preparation before the Bray-Curtis dissimilarity used by NMDS and PERMANOVA: hellinger (square root of relative abundances), none, rarefy or srs.',
   'analysis.beta.depth':          'Reads per sample for rarefy or srs; 0 uses the smallest non-empty sample.',
-  'analysis.differential.method': 'Model fitted to each taxon: nb_glm (negative-binomial GLM on read counts; effect is a log2 fold change) or clr_welch (Welch\'s t-test on centred log-ratios over the taxa passing min_prevalence, after zCompositions::cmultRepl zero replacement; effect is a CLR difference, not a fold change).',
-  'analysis.differential.offset': 'Library-size offset in each negative-binomial model: tss (total reads per sample) or rle (median-of-ratios over taxa present in every sample; refused when there are none). Used by nb_glm only.',
+  'analysis.differential.method': 'Model per taxon: nb_glm (negative-binomial GLM); clr_welch (Welch test on CLR values after zero replacement); zinb (zero-inflated negative binomial); hurdle_nb (logistic zero process plus zero-truncated negative binomial). Two-part NB p-values jointly test group effects in both processes; use only with evidence for excess zeros and interpret cautiously.',
+  'analysis.differential.offset': 'Library-size offset in the count component of nb_glm, zinb and hurdle_nb: tss (total reads per sample) or rle (median-of-ratios over taxa present in every sample; refused when there are none).',
   'analysis.differential.dispersion': 'nb_glm only: ml estimates theta separately with MASS::glm.nb; glmGamPoi uses a cross-taxon dispersion trend in fixed-scale negative-binomial Wald fits.',
   'analysis.differential.min_prevalence': 'Fraction of samples, 0 to 1, in which a taxon must have reads to be tested; taxa below it are listed as filtered.',
   'analysis.differential.replacement_delta': 'clr_welch only: cmultRepl\'s frac, between 0 and 1; an imputed proportion above its taxon\'s smallest observed proportion is lowered to this fraction of it. Replaced zeros are not measurements.',
@@ -251,7 +251,7 @@ export const CONFIG_TYPES: Record<string, ConfigType> = {
   'analysis.alpha.rarefaction_iterations': { kind: 'int' },
   'analysis.beta.normalisation':  { kind: 'enum', options: ['hellinger', 'none', 'rarefy', 'srs'] },
   'analysis.beta.depth':          { kind: 'int' },
-  'analysis.differential.method': { kind: 'enum', options: ['nb_glm', 'clr_welch'] },
+  'analysis.differential.method': { kind: 'enum', options: ['nb_glm', 'clr_welch', 'zinb', 'hurdle_nb'] },
   'analysis.differential.offset': { kind: 'enum', options: ['tss', 'rle'] },
   'analysis.differential.dispersion': { kind: 'enum', options: ['ml', 'glmGamPoi'] },
   'analysis.differential.min_prevalence': { kind: 'float', min: 0, max: 1, step: 0.05 },
diff --git a/frontend/src/components/differentialTable.test.ts b/frontend/src/components/differentialTable.test.ts
index f2143e07..770a9994 100644
--- a/frontend/src/components/differentialTable.test.ts
+++ b/frontend/src/components/differentialTable.test.ts
@@ -38,9 +38,9 @@ describe('differential CSV', () => {
   test('keeps every taxon, full precision, and leaves missing statistics empty', () => {
     const lines = differentialCsv(result).trimEnd().split('\r\n')
     expect(lines[0]).toContain('run_B vs run_A (reference)')
-    expect(lines[3]).toBe('taxon,status,log2_fold_change,estimate,standard_error,statistic,pvalue,padj,dispersion_theta,prevalence,note')
-    expect(lines[4]).toBe('Bacteroides,ok,2,1.3862943611198906,0.25,5.545,2.9e-8,5.8e-8,12.5,1,')
-    expect(lines[5].startsWith('"Odd, ""quoted""",failed,,,,,,,,0.5,')).toBe(true)
+    expect(lines[3]).toBe('taxon,status,log2_fold_change,estimate,standard_error,statistic,pvalue,padj,dispersion_theta,zero_process_log_odds,prevalence,note')
+    expect(lines[4]).toBe('Bacteroides,ok,2,1.3862943611198906,0.25,5.545,2.9e-8,5.8e-8,12.5,,1,')
+    expect(lines[5].startsWith('"Odd, ""quoted""",failed,,,,,,,,,0.5,')).toBe(true)
     expect(lines).toHaveLength(6)
   })
 })
diff --git a/frontend/src/components/differentialTable.ts b/frontend/src/components/differentialTable.ts
index 90a547a8..42e69c59 100644
--- a/frontend/src/components/differentialTable.ts
+++ b/frontend/src/components/differentialTable.ts
@@ -4,7 +4,7 @@ import type { DifferentialResult, DifferentialRow } from '../api/types'
 
 const COLUMNS = [
   'taxon', 'status', 'log2_fold_change', 'estimate', 'standard_error', 'statistic',
-  'pvalue', 'padj', 'dispersion_theta', 'prevalence', 'note',
+  'pvalue', 'padj', 'dispersion_theta', 'zero_process_log_odds', 'prevalence', 'note',
 ] as const
 
 /** Quote one CSV field when it holds a comma, quote or line break (RFC 4180). */
@@ -27,7 +27,7 @@ export function differentialCsv(result: DifferentialResult): string {
     `# ${result.method}`,
     `# ${configLine(result)}`,
   ].map(l => l.replace(/[\r\n]+/g, ' '))
-  const rows = result.rows.map(r => COLUMNS.map(c => csvField(r[c])).join(','))
+  const rows = result.rows.map(r => COLUMNS.map(c => csvField(r[c] ?? null)).join(','))
   return [...header, COLUMNS.join(','), ...rows].join('\r\n') + '\r\n'
 }
 
@@ -48,7 +48,7 @@ export function configLine(result: DifferentialResult): string {
   const c = result.config
   return c.method === 'clr_welch'
     ? `method clr_welch, replacement_delta ${c.replacement_delta}, min_prevalence ${c.min_prevalence}; estimate is a CLR difference, not a fold change`
-    : `method nb_glm, offset ${c.offset}, min_prevalence ${c.min_prevalence}`
+    : `method ${c.method}, offset ${c.offset}, min_prevalence ${c.min_prevalence}`
 }
 
 /** The effect size of one row, from the field the result names in `effect`. */
diff --git a/src/analysis/differential.jl b/src/analysis/differential.jl
index 66cbf471..7d761dcf 100644
--- a/src/analysis/differential.jl
+++ b/src/analysis/differential.jl
@@ -24,10 +24,10 @@ using ..Analysis: _palette_hex, R_WAIT_SECONDS
 export DifferentialConfig, ScalingRefusal, MASSUnavailable, ZCompositionsUnavailable,
        GlmGamPoiUnavailable, OFFSET_METHODS, METHODS, tss_factors, rle_factors,
        size_factors, bh_adjust, validate_counts, replace_zeros, clr_transform,
-       fit_nb, fit_welch, differential_abundance, volcano_chart
+       fit_nb, fit_two_part_nb, fit_occupancy, fit_welch, differential_abundance, volcano_chart
 
 # Methods accepted under `analysis.differential.method`.
-const METHODS = ("nb_glm", "clr_welch")
+const METHODS = ("nb_glm", "clr_welch", "zinb", "hurdle_nb")
 
 # The per-taxon dispersion `nb_glm` uses: maximum likelihood (MASS::glm.nb's
 # own theta) or glmGamPoi's shrinkage estimator, which shares information
@@ -62,10 +62,12 @@ above its taxon's smallest observed proportion is lowered to this fraction of it
 `dispersion` names the per-taxon dispersion `nb_glm` uses (`"ml"` or
 `"glmGamPoi"`): `"ml"` takes the maximum-likelihood theta from `MASS::glm.nb`,
 `"glmGamPoi"` uses `overdispersion_shrinkage_list.dispersion_trend` from
-`glmGamPoi::glm_gp` for
-the taxa passing the prevalence filter, with `theta = 1 / dispersion` fixed
+`glmGamPoi::glm_gp` for the taxa passing the prevalence filter, with `theta = 1 / dispersion` fixed
 in each negative-binomial GLM. Both options use a Wald test of the group
-coefficient with fixed scale 1; `"glmGamPoi"` is refused for `clr_welch`.
+coefficient with fixed scale 1; `"glmGamPoi"` is refused outside `nb_glm`. `zinb` and `hurdle_nb`
+fit a zero-process and negative-binomial count process per taxon, then use a 2-df likelihood-ratio
+test for the group effect in both processes. Their count-component fold change is reported; the zero
+process effect remains in diagnostics.
 The trend is plugged into an NB Wald model; its estimation uncertainty is not
 propagated. This option does not implement glmGamPoi's quasi-likelihood test.
 Method and offset names are stripped and lowercased; dispersion names are
@@ -514,6 +516,212 @@ function fit_nb(counts::AbstractMatrix{<:Integer}, groups::AbstractVector{<:Abst
               theta     = _num.(raw.theta))
 end
 
+# Two-part NB models are implemented with base R's optim/stats likelihoods: no
+# pscl/glmmTMB dependency is introduced. The joint 2-df LRT tests the group term
+# in both components; a count-part Wald estimate alone would miss zero-process
+# changes. The models are deliberately fixed-effect, per-taxon models.
+const _FIT_TWO_PART_R = raw"""
+da_g01 <- as.numeric(da_group == da_levels[2])
+da_n <- ncol(da_counts)
+da_result <- data.frame(status=rep("ok", da_n), note=rep("", da_n),
+  estimate=rep(NA_real_, da_n), se=rep(NA_real_, da_n), statistic=rep(NA_real_, da_n),
+  pvalue=rep(NA_real_, da_n), theta=rep(NA_real_, da_n), zero_estimate=rep(NA_real_, da_n))
+da_tp_logspace_add <- function(a,b) { m <- pmax(a,b); m + log(exp(a-m)+exp(b-m)) }
+da_tp_nll <- function(par, y, g, off, kind, null=FALSE) {
+  if (null) { b0 <- par[1]; b1 <- 0; z0 <- par[2]; z1 <- 0; lt <- par[3] }
+  else { b0 <- par[1]; b1 <- par[2]; z0 <- par[3]; z1 <- par[4]; lt <- par[5] }
+  if (!is.finite(lt) || lt < -15 || lt > 15) return(1e100)
+  theta <- exp(lt)
+  mu <- exp(pmax(-30, pmin(30, b0 + b1*g + off)))
+  q <- plogis(pmax(-30, pmin(30, z0 + z1*g)))
+  nb <- dnbinom(y, size=theta, mu=mu, log=TRUE)
+  nb0 <- dnbinom(0, size=theta, mu=mu)
+  if (kind == "zinb") {
+    ll <- ifelse(y == 0, da_tp_logspace_add(log(q), log1p(-q) + nb), log1p(-q) + nb)
+  } else {
+    ll <- ifelse(y == 0, log(q), log1p(-q) + nb - log1p(-nb0))
+  }
+  if (any(!is.finite(ll))) return(1e100)
+  -sum(ll)
+}
+for (da_j in seq_len(da_n)) {
+  da_tp_y <- da_counts[,da_j]
+  if (length(unique(da_tp_y)) < 2L) { da_result[da_j,"status"] <- "failed"; da_result[da_j,"note"] <- "counts are constant"; next }
+  da_tp_init <- c(log((sum(da_tp_y)+0.5)/sum(exp(da_offset))), 0,
+                  qlogis(min(.95,max(.05,mean(da_tp_y==0)*.8))), 0, 0)
+  da_tp_fit <- tryCatch(optim(da_tp_init, da_tp_nll, y=da_tp_y, g=da_g01,
+    off=da_offset, kind=da_kind, method="BFGS", hessian=TRUE,
+    control=list(maxit=1000, reltol=1e-9)), error=function(e)e)
+  da_tp_null_init <- da_tp_init[c(1,3,5)]
+  da_tp_nullfit <- tryCatch(optim(da_tp_null_init, da_tp_nll, y=da_tp_y, g=da_g01,
+    off=da_offset, kind=da_kind, null=TRUE, method="BFGS",
+    control=list(maxit=1000, reltol=1e-9)), error=function(e)e)
+  if (inherits(da_tp_fit,"error") || inherits(da_tp_nullfit,"error") ||
+      da_tp_fit$convergence != 0 || da_tp_nullfit$convergence != 0 ||
+      !is.finite(da_tp_fit$value) || !is.finite(da_tp_nullfit$value)) {
+    da_result[da_j,"status"] <- "failed"; da_result[da_j,"note"] <- "two-part likelihood did not converge"; next
+  }
+  da_tp_info <- tryCatch(solve(da_tp_fit$hessian), error=function(e)NULL)
+  if (is.null(da_tp_info) || any(!is.finite(da_tp_info)) || any(diag(da_tp_info)<=0)) {
+    da_result[da_j,"status"] <- "failed"; da_result[da_j,"note"] <- "singular or non-positive likelihood Hessian"; next
+  }
+  da_tp_delta <- da_tp_nullfit$value-da_tp_fit$value
+  if (da_tp_delta < -1e-6) {
+    da_result[da_j,"status"] <- "failed"; da_result[da_j,"note"] <- "full model did not attain the nested null-model likelihood"; next
+  }
+  da_tp_lrt <- max(0,2*da_tp_delta)
+  da_result[da_j,"estimate"] <- da_tp_fit$par[2]
+  da_result[da_j,"se"] <- sqrt(da_tp_info[2,2])
+  da_result[da_j,"statistic"] <- da_tp_lrt
+  da_result[da_j,"pvalue"] <- pchisq(da_tp_lrt, df=2, lower.tail=FALSE)
+  da_result[da_j,"theta"] <- exp(da_tp_fit$par[5])
+  da_result[da_j,"zero_estimate"] <- da_tp_fit$par[4]
+}
+"""
+
+"""
+    fit_two_part_nb(counts, groups, offset; levels, kind) -> DataFrame
+
+Fit per-taxon zero-inflated (`kind="zinb"`) or hurdle (`kind="hurdle_nb"`)
+negative-binomial likelihoods. `offset` is the log size factor for the NB
+count component. Both logistic components and the count mean include the
+binary group effect. The returned p-value is the joint 2-df likelihood-ratio
+test against a model without either group effect; `estimate` is the count
+log-fold change and `zero_estimate` is the group log-odds coefficient for the
+zero/occurrence component. Requires at least two samples in each group.
+"""
+function fit_two_part_nb(counts::AbstractMatrix{<:Integer}, groups::AbstractVector{<:AbstractString},
+                         offset::AbstractVector{<:Real}; levels::NTuple{2,String}, kind::AbstractString)
+    kind in ("zinb", "hurdle_nb") || throw(ArgumentError("kind must be zinb or hurdle_nb"))
+    size(counts,1) == length(groups) == length(offset) || throw(ArgumentError("counts, groups and offset must align"))
+    all(>=(0), counts) || throw(ArgumentError("counts must be non-negative"))
+    all(isfinite, offset) || throw(ArgumentError("count offsets must be finite"))
+    all(g -> g in levels, groups) || throw(ArgumentError("groups must be exactly $(levels)"))
+    all(g -> count(==(g), groups) >= 2, levels) || throw(ArgumentError("both groups need at least two samples for a two-part model"))
+    raw = with_r_lock(; timeout=R_WAIT_SECONDS[]) do
+        RCall.globalEnv[:da_counts] = Matrix{Int}(counts)
+        RCall.globalEnv[:da_group] = String.(groups)
+        RCall.globalEnv[:da_levels] = collect(levels)
+        RCall.globalEnv[:da_offset] = Float64.(offset)
+        RCall.globalEnv[:da_kind] = kind == "zinb" ? "zinb" : "hurdle"
+        try
+            RCall.reval(_FIT_TWO_PART_R)
+            DataFrame(RCall.rcopy(RCall.reval("da_result")))
+        finally
+            RCall.reval("rm(list=intersect(ls(), c('da_counts','da_group','da_levels','da_offset','da_kind','da_g01','da_n','da_result','da_tp_logspace_add','da_tp_nll','da_tp_y','da_tp_init','da_tp_fit','da_tp_null_init','da_tp_nullfit','da_tp_info','da_tp_delta','da_tp_lrt'))); invisible(gc())")
+        end
+    end
+    DataFrame(status=String.(raw.status), note=String.(raw.note), estimate=_num.(raw.estimate),
+              se=_num.(raw.se), statistic=_num.(raw.statistic), pvalue=_num.(raw.pvalue),
+              theta=_num.(raw.theta), zero_estimate=_num.(raw.zero_estimate))
+end
+
+# Single-visit microbiome tables cannot distinguish true absence from missed
+# detection without external assumptions. This estimator therefore accepts an
+# explicit site × visit × taxon detection-history array and refuses fewer than
+# two visits. Detection is constant across visits in v1; ψ has an intercept and
+# group effect. This is intentionally a separate API from differential_abundance
+# until the UI can collect a repeated-visit design instead of inventing one.
+const _FIT_OCCUPANCY_R = raw"""
+da_occ_logspace_add <- function(a,b) { m <- pmax(a,b); m + log(exp(a-m)+exp(b-m)) }
+da_n <- dim(da_history)[3]
+da_occ_result <- data.frame(status=rep("ok",da_n), note=rep("",da_n),
+  estimate=rep(NA_real_,da_n), se=rep(NA_real_,da_n), statistic=rep(NA_real_,da_n),
+  pvalue=rep(NA_real_,da_n), psi_reference=rep(NA_real_,da_n),
+  psi_contrast=rep(NA_real_,da_n), detection_probability=rep(NA_real_,da_n))
+da_occ_nll <- function(par, h, g, null=FALSE) {
+  if (null) { b0 <- par[1]; b1 <- 0; a0 <- par[2] }
+  else { b0 <- par[1]; b1 <- par[2]; a0 <- par[3] }
+  psi <- plogis(pmax(-30,pmin(30,b0+b1*g)))
+  p <- plogis(pmax(-30,pmin(30,a0)))
+  k <- ncol(h); detected <- rowSums(h)>0
+  lp_detect <- rowSums(ifelse(h==1,log(p),log1p(-p)))
+  lp_none <- da_occ_logspace_add(log1p(-psi), log(psi)+k*log1p(-p))
+  ll <- ifelse(detected,log(psi)+lp_detect,lp_none)
+  if(any(!is.finite(ll))) return(1e100)
+  -sum(ll)
+}
+for (da_j in seq_len(da_n)) {
+  da_occ_h <- da_history[,,da_j]; storage.mode(da_occ_h) <- "integer"
+  if (all(rowSums(da_occ_h)==0)) {
+    da_occ_result[da_j,"status"] <- "failed"; da_occ_result[da_j,"note"] <- "no detections in this taxon's history"; next
+  }
+  da_occ_init <- c(qlogis(min(.95,max(.05,mean(rowSums(da_occ_h)>0)*1.5))), 0,
+                   qlogis(min(.95,max(.05,mean(da_occ_h)))))
+  da_occ_fit <- tryCatch(optim(da_occ_init,da_occ_nll,h=da_occ_h,g=da_g01,
+    method="BFGS",hessian=TRUE,control=list(maxit=1000,reltol=1e-10)),error=function(e)e)
+  da_occ_nullfit <- tryCatch(optim(da_occ_init[c(1,3)],da_occ_nll,h=da_occ_h,g=da_g01,
+    null=TRUE,method="BFGS",control=list(maxit=1000,reltol=1e-10)),error=function(e)e)
+  if (inherits(da_occ_fit,"error") || inherits(da_occ_nullfit,"error") ||
+      da_occ_fit$convergence!=0 || da_occ_nullfit$convergence!=0) {
+    da_occ_result[da_j,"status"] <- "failed"; da_occ_result[da_j,"note"] <- "occupancy likelihood did not converge"; next
+  }
+  da_occ_vc <- tryCatch(solve(da_occ_fit$hessian),error=function(e)NULL)
+  if (is.null(da_occ_vc) || any(!is.finite(da_occ_vc)) || da_occ_vc[2,2]<=0) {
+    da_occ_result[da_j,"status"] <- "failed"; da_occ_result[da_j,"note"] <- "singular occupancy likelihood Hessian"; next
+  }
+  da_occ_delta <- da_occ_nullfit$value-da_occ_fit$value
+  if (da_occ_delta < -1e-6) {
+    da_occ_result[da_j,"status"] <- "failed"; da_occ_result[da_j,"note"] <- "full occupancy model did not attain the nested null likelihood"; next
+  }
+  da_occ_lrt <- max(0,2*da_occ_delta)
+  da_occ_result[da_j,"estimate"] <- da_occ_fit$par[2]
+  da_occ_result[da_j,"se"] <- sqrt(da_occ_vc[2,2])
+  da_occ_result[da_j,"statistic"] <- da_occ_lrt
+  da_occ_result[da_j,"pvalue"] <- pchisq(da_occ_lrt,df=1,lower.tail=FALSE)
+  da_occ_result[da_j,"psi_reference"] <- plogis(da_occ_fit$par[1])
+  da_occ_result[da_j,"psi_contrast"] <- plogis(da_occ_fit$par[1]+da_occ_fit$par[2])
+  da_occ_result[da_j,"detection_probability"] <- plogis(da_occ_fit$par[3])
+}
+"""
+
+"""
+    fit_occupancy(history, groups, taxa; levels) -> DataFrame
+
+Single-species occupancy likelihood with imperfect detection and repeated visits.
+`history` is a sites × visits × taxa Bool array. At least two visits are
+required: single-visit histories do not separate occupancy from detection
+without strong external assumptions. Detection probability is constant across
+visits in this initial implementation; occupancy has an intercept and a
+binary reference/contrast effect. Results include ψ estimates, p, the occupancy
+log-odds contrast and an LRT p-value. This is a lower-level API; it is not
+silently applied to ordinary one-observation-per-sample count tables.
+"""
+function fit_occupancy(history::AbstractArray{Bool,3}, groups::AbstractVector{<:AbstractString},
+                       taxa::AbstractVector{<:AbstractString}; levels::NTuple{2,String})
+    sites, visits, n_taxa = size(history)
+    sites > 0 && n_taxa > 0 || throw(ArgumentError("occupancy needs at least one site and one taxon"))
+    visits >= 2 || throw(ArgumentError("occupancy needs at least two repeat visits per site; single-visit detection and occupancy are not separately identifiable"))
+    length(groups) == sites || throw(ArgumentError("occupancy history has $sites sites but $(length(groups)) group labels"))
+    length(taxa) == n_taxa || throw(ArgumentError("occupancy history has $n_taxa taxa but $(length(taxa)) taxon names"))
+    all(g -> g in levels, groups) || throw(ArgumentError("occupancy groups must be exactly $(levels)"))
+    all(g -> count(==(g), groups) >= 2, levels) || throw(ArgumentError("each occupancy group needs at least two sites"))
+    raw = with_r_lock(; timeout=R_WAIT_SECONDS[]) do
+        RCall.globalEnv[:da_history] = Int.(history)
+        RCall.globalEnv[:da_group] = String.(groups)
+        RCall.globalEnv[:da_levels] = collect(levels)
+        RCall.globalEnv[:da_g01] = Float64.(groups .== levels[2])
+        try
+            RCall.reval(_FIT_OCCUPANCY_R)
+            DataFrame(RCall.rcopy(RCall.reval("da_occ_result")))
+        finally
+            RCall.reval("rm(list=intersect(ls(),c('da_history','da_group','da_levels','da_g01','da_n','da_occ_result','da_occ_logspace_add','da_occ_nll','da_occ_h','da_occ_init','da_occ_fit','da_occ_nullfit','da_occ_vc','da_occ_delta','da_occ_lrt'))); invisible(gc())")
+        end
+    end
+    p = [_num(x) for x in raw.pvalue]
+    family = findall(x -> !isnothing(x) && raw.status[x] == "ok", eachindex(p))
+    adjusted = Union{Nothing,Float64}[nothing for _ in eachindex(p)]
+    if !isempty(family)
+        values = bh_adjust(Float64[p[i] for i in family])
+        for (k,i) in enumerate(family); adjusted[i] = values[k]; end
+    end
+    DataFrame(status=String.(raw.status), note=String.(raw.note),
+      occupancy_log_odds=_num.(raw.estimate), standard_error=_num.(raw.se),
+      statistic=_num.(raw.statistic), pvalue=p, padj=adjusted,
+      psi_reference=_num.(raw.psi_reference), psi_contrast=_num.(raw.psi_contrast),
+      detection_probability=_num.(raw.detection_probability), taxon=String.(taxa))
+end
+
 """
     clr_transform(x) -> Matrix{Float64}
 
@@ -761,6 +969,15 @@ change, relative to the geometric mean of the tested taxa. Replaced zeros are
 not measurements; `diagnostics.zero_replacement` reports their share of each
 sample. Each group needs at least 2 samples.
 
+`zinb` and `hurdle_nb` use a count-process negative binomial with a library-size
+offset and a group-specific logistic zero/occurrence process. Their reported
+count effect is the NB coefficient; the p-value is a joint 2-df likelihood-ratio
+test over both group coefficients. ZINB treats zeros as a mixture of structural
+and sampling zeros. Hurdle NB assigns all zeros to its logistic component and
+uses a zero-truncated NB for positive counts. These models can overfit; compare
+against NB diagnostics and do not interpret the zero component as biological
+absence without validation.
+
 For both methods, taxa present in fewer than `config.min_prevalence` of the
 samples are reported as `filtered` and not fitted, and the Benjamini-Hochberg
 family is the fitted taxa that produced a p-value; failed and filtered taxa
@@ -834,7 +1051,7 @@ function differential_abundance(counts::AbstractMatrix{<:Real},
         config_echo = Dict{String,Any}("method" => config.method, "offset" => config.offset,
                                        "min_prevalence" => config.min_prevalence,
                                        "dispersion" => config.dispersion)
-    else
+    elseif config.method == "clr_welch"
         (n_ref >= 2 && n_con >= 2) || throw(ArgumentError(
             "Welch's t-test needs at least 2 samples in each group to estimate its variance; " *
             "'$reference' has $n_ref and '$contrast' has $n_con"))
@@ -875,6 +1092,37 @@ function differential_abundance(counts::AbstractMatrix{<:Real},
             "median_imputed_fraction" => _median(sorted_mass),
             "max_imputed_fraction" => last(sorted_mass),
             "warnings" => replaced.warnings)
+    else
+        factors = size_factors(counts, samples, config.offset)
+        tested = findall(>=(config.min_prevalence), prevalence)
+        fits = isempty(tested) ? nothing :
+            fit_two_part_nb(Int.(counts[:, tested]), groups, log.(factors); levels,
+                            kind=config.method)
+        if !isnothing(fits)
+            for (k, j) in enumerate(tested)
+                f = fits[k, :]
+                merge!(rows[j], Dict{String,Any}(
+                    "status" => f.status, "note" => f.note, "estimate" => f.estimate,
+                    "log2_fold_change" => isnothing(f.estimate) ? nothing : f.estimate / log(2),
+                    "standard_error" => f.se, "statistic" => f.statistic,
+                    "pvalue" => f.pvalue, "dispersion_theta" => f.theta,
+                    "zero_process_log_odds" => f.zero_estimate))
+            end
+        end
+        model_label = config.method == "zinb" ? "zero-inflated negative-binomial" :
+                      "hurdle negative-binomial"
+        method_text = "$model_label two-part model per taxon with TSS/RLE offset in the count process; " *
+                      "joint 2-df likelihood-ratio test of group effects in the zero and count processes; " *
+                      "Benjamini-Hochberg adjustment over fitted taxa. The reported fold change is the " *
+                      "count-process effect; zero_process_log_odds is the group effect in the zero process."
+        effect = Dict("key" => "log2_fold_change", "label" => "count-process log2 fold change")
+        size_factor_rows = [Dict("sample" => String(samples[i]), "group" => String(groups[i]),
+                                 "factor" => factors[i]) for i in eachindex(samples)]
+        config_echo = Dict{String,Any}("method" => config.method, "offset" => config.offset,
+                                       "min_prevalence" => config.min_prevalence)
+        diagnostics["model_warning"] = config.method == "zinb" ?
+            "Use ZINB only when diagnostics support excess zeros beyond NB; zero inflation and NB overdispersion can be difficult to distinguish." :
+            "Hurdle NB assumes every zero comes from the occurrence process and positive counts follow a zero-truncated NB."
     end
 
     family = findall(r -> r["status"] in ("ok", "boundary") && !isnothing(r["pvalue"]), rows)
diff --git a/src/server/routes/config.jl b/src/server/routes/config.jl
index f4f2ba0b..6fa9faf5 100644
--- a/src/server/routes/config.jl
+++ b/src/server/routes/config.jl
@@ -197,8 +197,8 @@ function _differential_pair_error(resolved::AbstractDict, changes::AbstractDict)
     method = value(keys_[1], "nb_glm")
     dispersion = value(keys_[2], "ml")
     if method isa AbstractString && dispersion isa AbstractString &&
-       lowercase(strip(method)) == "clr_welch" && lowercase(strip(dispersion)) == "glmgampoi"
-        return "analysis.differential.dispersion: glmGamPoi requires method: nb_glm; use dispersion: ml with clr_welch"
+       lowercase(strip(method)) != "nb_glm" && lowercase(strip(dispersion)) == "glmgampoi"
+        return "analysis.differential.dispersion: glmGamPoi requires method: nb_glm; use dispersion: ml with $(strip(method))"
     end
     nothing
 end
diff --git a/test/unit/test_differential.jl b/test/unit/test_differential.jl
index dcab0549..09c31703 100644
--- a/test/unit/test_differential.jl
+++ b/test/unit/test_differential.jl
@@ -10,6 +10,7 @@ using MetaManifold
 using MetaManifold.Differential
 using RCall
 using Test
+using Random
 using DuckDB, DBInterface, JSON3, HTTP, YAML
 using OrderedCollections: OrderedDict
 
@@ -163,6 +164,7 @@ end
     end
 
     @testset "Configuration" begin
+        @test Differential.METHODS == ("nb_glm", "clr_welch", "zinb", "hurdle_nb")
         @test DifferentialConfig().offset == "tss"
         @test DifferentialConfig().min_prevalence == 0.0
         @test DifferentialConfig("rle", 0.5).offset == "rle"
@@ -172,6 +174,9 @@ end
         @test DifferentialConfig().method == "nb_glm"
         @test DifferentialConfig().replacement_delta == 0.65
         @test DifferentialConfig(; method="clr_welch", replacement_delta=0.3).replacement_delta == 0.3
+        @test DifferentialConfig(; method="zinb").method == "zinb"
+        @test DifferentialConfig(; method="hurdle_nb").method == "hurdle_nb"
+        @test_throws ArgumentError DifferentialConfig(; method="occupancy") # repeated visits use fit_occupancy
         @test_throws ArgumentError DifferentialConfig(; method="ancombc")
         # The pre-review name, never released, is not accepted.
         @test_throws ArgumentError DifferentialConfig(; method="clr_lm")
@@ -182,6 +187,7 @@ end
         @test DifferentialConfig(dispersion="  GLMGAMPOI ").dispersion == "glmGamPoi"
         @test_throws ArgumentError DifferentialConfig(dispersion="shrinkage")
         @test_throws ArgumentError DifferentialConfig(method="clr_welch", dispersion="glmGamPoi")
+        @test_throws ArgumentError DifferentialConfig(method="zinb", dispersion="glmGamPoi")
         @test DifferentialConfig(method="clr_welch").dispersion == "ml"
         for d in (0, 1, -0.5, NaN)
             @test_throws ArgumentError DifferentialConfig(; method="clr_welch", replacement_delta=d)
@@ -647,6 +653,57 @@ end
         end
     end
 
+    @testset "two-part negative-binomial models" begin
+        if !HAVE_R
+            @test_skip false
+        else
+            rng = MersenneTwister(18)
+            g = vcat(fill("A", 60), fill("B", 60))
+            x = zeros(Int, 120, 2)
+            for i in eachindex(g), j in 1:2
+                pzero = j == 1 ? (g[i] == "A" ? 0.55 : 0.25) : 0.4
+                x[i,j] = rand(rng) < pzero ? 0 : 1 + rand(rng, 0:5) + (j == 1 && g[i] == "B" ? 3 : 0)
+            end
+            for kind in ("zinb", "hurdle_nb")
+                fitted = fit_two_part_nb(x, g, zeros(120); levels=("A", "B"), kind)
+                @test nrow(fitted) == 2
+                @test all(fitted.status .== "ok")
+                @test all(isfinite, fitted.pvalue)
+                @test all((0 .<= fitted.pvalue) .& (fitted.pvalue .<= 1))
+                @test all(isfinite, fitted.zero_estimate)
+                @test fitted.estimate[1] > 0
+            end
+        end
+    end
+
+    @testset "occupancy requires repeated detection histories" begin
+        rng = MersenneTwister(29)
+        history = falses(100, 4, 100)
+        g = vcat(fill("A", 50), fill("B", 50))
+        for i in eachindex(g), j in axes(history,3)
+            occupied = rand(rng) < 0.7
+            for v in 1:4
+                history[i,v,j] = occupied && rand(rng) < 0.55
+            end
+        end
+        taxa = ["taxon_$j" for j in axes(history,3)]
+        @test_throws ArgumentError fit_occupancy(history[:,1:1,:], g, taxa; levels=("A","B"))
+        @test_throws ArgumentError fit_occupancy(history, g[1:end-1], taxa; levels=("A","B"))
+        if !HAVE_R
+            @test_skip false
+        else
+            fitted = fit_occupancy(history, g, taxa; levels=("A","B"))
+            @test nrow(fitted) == 2
+            @test all(fitted.status .== "ok")
+            @test all((0 .<= fitted.detection_probability) .& (fitted.detection_probability .<= 1))
+            @test all((0 .<= fitted.psi_reference) .& (fitted.psi_reference .<= 1))
+            @test all((0 .<= fitted.psi_contrast) .& (fitted.psi_contrast .<= 1))
+            @test sum(fitted.detection_probability) / length(fitted.detection_probability) ≈ 0.55 atol=0.2
+            @test sum(fitted.psi_reference) / length(fitted.psi_reference) ≈ 0.7 atol=0.1
+            @test all((0 .<= fitted.padj) .& (fitted.padj .<= 1))
+        end
+    end
+
     @testset "CLR transform" begin
         @test clr_transform([1.0 exp(1) exp(2)]) ≈ [-1.0 0.0 1.0]
         z = clr_transform([1.0 2 3 4; 5 1 1 9; 0.5 0.25 8 2])

From ce961fc2f0f57e34dd470c058546c463b254409e Mon Sep 17 00:00:00 2001
From: "Jonathan D.A. Jewell" <6759885+hyperpolymath@users.noreply.github.com>
Date: Thu, 22 Oct 2026 23:34:09 +0100
Subject: [PATCH 2/6] Update test/unit/test_differential.jl

Co-authored-by: coderabbitai[bot] <136622811+coderabbitai[bot]@users.noreply.github.com>
Signed-off-by: Jonathan D.A. Jewell <6759885+hyperpolymath@users.noreply.github.com>
---
 test/unit/test_differential.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/unit/test_differential.jl b/test/unit/test_differential.jl
index 09c31703..5c260c2d 100644
--- a/test/unit/test_differential.jl
+++ b/test/unit/test_differential.jl
@@ -693,7 +693,7 @@ end
             @test_skip false
         else
             fitted = fit_occupancy(history, g, taxa; levels=("A","B"))
-            @test nrow(fitted) == 2
+            @test nrow(fitted) == length(taxa)
             @test all(fitted.status .== "ok")
             @test all((0 .<= fitted.detection_probability) .& (fitted.detection_probability .<= 1))
             @test all((0 .<= fitted.psi_reference) .& (fitted.psi_reference .<= 1))

From 2a0ddccbdf1e877e3035a4e91b038e86b0cdfdde Mon Sep 17 00:00:00 2001
From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com>
Date: Thu, 22 Oct 2026 22:34:19 +0000
Subject: [PATCH 3/6] fix(analysis): unblock model fitting and Julia CI

Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com>
---
 docs/occupancy-models.md     |  8 +++++---
 src/analysis/differential.jl | 25 +++++++++++++------------
 src/core/doi.jl              |  2 ++
 3 files changed, 20 insertions(+), 15 deletions(-)

diff --git a/docs/occupancy-models.md b/docs/occupancy-models.md
index 7e22880c..86e42a8c 100644
--- a/docs/occupancy-models.md
+++ b/docs/occupancy-models.md
@@ -42,9 +42,11 @@ 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. A per-taxon likelihood-ratio
-test assesses the group effect in occupancy; p-values are BH-adjusted. The
-returned table includes `ψ` in each group and the estimated `p`.
+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
diff --git a/src/analysis/differential.jl b/src/analysis/differential.jl
index 7d761dcf..6015b0d0 100644
--- a/src/analysis/differential.jl
+++ b/src/analysis/differential.jl
@@ -527,13 +527,13 @@ da_result <- data.frame(status=rep("ok", da_n), note=rep("", da_n),
   estimate=rep(NA_real_, da_n), se=rep(NA_real_, da_n), statistic=rep(NA_real_, da_n),
   pvalue=rep(NA_real_, da_n), theta=rep(NA_real_, da_n), zero_estimate=rep(NA_real_, da_n))
 da_tp_logspace_add <- function(a,b) { m <- pmax(a,b); m + log(exp(a-m)+exp(b-m)) }
-da_tp_nll <- function(par, y, g, off, kind, null=FALSE) {
+da_tp_nll <- function(par, y, grp, off, kind, null=FALSE) {
   if (null) { b0 <- par[1]; b1 <- 0; z0 <- par[2]; z1 <- 0; lt <- par[3] }
   else { b0 <- par[1]; b1 <- par[2]; z0 <- par[3]; z1 <- par[4]; lt <- par[5] }
   if (!is.finite(lt) || lt < -15 || lt > 15) return(1e100)
   theta <- exp(lt)
-  mu <- exp(pmax(-30, pmin(30, b0 + b1*g + off)))
-  q <- plogis(pmax(-30, pmin(30, z0 + z1*g)))
+  mu <- exp(pmax(-30, pmin(30, b0 + b1*grp + off)))
+  q <- plogis(pmax(-30, pmin(30, z0 + z1*grp)))
   nb <- dnbinom(y, size=theta, mu=mu, log=TRUE)
   nb0 <- dnbinom(0, size=theta, mu=mu)
   if (kind == "zinb") {
@@ -549,11 +549,11 @@ for (da_j in seq_len(da_n)) {
   if (length(unique(da_tp_y)) < 2L) { da_result[da_j,"status"] <- "failed"; da_result[da_j,"note"] <- "counts are constant"; next }
   da_tp_init <- c(log((sum(da_tp_y)+0.5)/sum(exp(da_offset))), 0,
                   qlogis(min(.95,max(.05,mean(da_tp_y==0)*.8))), 0, 0)
-  da_tp_fit <- tryCatch(optim(da_tp_init, da_tp_nll, y=da_tp_y, g=da_g01,
+  da_tp_fit <- tryCatch(optim(da_tp_init, da_tp_nll, y=da_tp_y, grp=da_g01,
     off=da_offset, kind=da_kind, method="BFGS", hessian=TRUE,
     control=list(maxit=1000, reltol=1e-9)), error=function(e)e)
   da_tp_null_init <- da_tp_init[c(1,3,5)]
-  da_tp_nullfit <- tryCatch(optim(da_tp_null_init, da_tp_nll, y=da_tp_y, g=da_g01,
+  da_tp_nullfit <- tryCatch(optim(da_tp_null_init, da_tp_nll, y=da_tp_y, grp=da_g01,
     off=da_offset, kind=da_kind, null=TRUE, method="BFGS",
     control=list(maxit=1000, reltol=1e-9)), error=function(e)e)
   if (inherits(da_tp_fit,"error") || inherits(da_tp_nullfit,"error") ||
@@ -629,13 +629,13 @@ da_occ_result <- data.frame(status=rep("ok",da_n), note=rep("",da_n),
   estimate=rep(NA_real_,da_n), se=rep(NA_real_,da_n), statistic=rep(NA_real_,da_n),
   pvalue=rep(NA_real_,da_n), psi_reference=rep(NA_real_,da_n),
   psi_contrast=rep(NA_real_,da_n), detection_probability=rep(NA_real_,da_n))
-da_occ_nll <- function(par, h, g, null=FALSE) {
+da_occ_nll <- function(par, history, grp, null=FALSE) {
   if (null) { b0 <- par[1]; b1 <- 0; a0 <- par[2] }
   else { b0 <- par[1]; b1 <- par[2]; a0 <- par[3] }
-  psi <- plogis(pmax(-30,pmin(30,b0+b1*g)))
+  psi <- plogis(pmax(-30,pmin(30,b0+b1*grp)))
   p <- plogis(pmax(-30,pmin(30,a0)))
-  k <- ncol(h); detected <- rowSums(h)>0
-  lp_detect <- rowSums(ifelse(h==1,log(p),log1p(-p)))
+  k <- ncol(history); detected <- rowSums(history)>0
+  lp_detect <- rowSums(ifelse(history==1,log(p),log1p(-p)))
   lp_none <- da_occ_logspace_add(log1p(-psi), log(psi)+k*log1p(-p))
   ll <- ifelse(detected,log(psi)+lp_detect,lp_none)
   if(any(!is.finite(ll))) return(1e100)
@@ -648,9 +648,9 @@ for (da_j in seq_len(da_n)) {
   }
   da_occ_init <- c(qlogis(min(.95,max(.05,mean(rowSums(da_occ_h)>0)*1.5))), 0,
                    qlogis(min(.95,max(.05,mean(da_occ_h)))))
-  da_occ_fit <- tryCatch(optim(da_occ_init,da_occ_nll,h=da_occ_h,g=da_g01,
+  da_occ_fit <- tryCatch(optim(da_occ_init,da_occ_nll,history=da_occ_h,grp=da_g01,
     method="BFGS",hessian=TRUE,control=list(maxit=1000,reltol=1e-10)),error=function(e)e)
-  da_occ_nullfit <- tryCatch(optim(da_occ_init[c(1,3)],da_occ_nll,h=da_occ_h,g=da_g01,
+  da_occ_nullfit <- tryCatch(optim(da_occ_init[c(1,3)],da_occ_nll,history=da_occ_h,grp=da_g01,
     null=TRUE,method="BFGS",control=list(maxit=1000,reltol=1e-10)),error=function(e)e)
   if (inherits(da_occ_fit,"error") || inherits(da_occ_nullfit,"error") ||
       da_occ_fit$convergence!=0 || da_occ_nullfit$convergence!=0) {
@@ -683,7 +683,8 @@ Single-species occupancy likelihood with imperfect detection and repeated visits
 required: single-visit histories do not separate occupancy from detection
 without strong external assumptions. Detection probability is constant across
 visits in this initial implementation; occupancy has an intercept and a
-binary reference/contrast effect. Results include ψ estimates, p, the occupancy
+binary reference/contrast effect. Results include `psi_reference` for `levels[1]`,
+`psi_contrast` for `levels[2]`, the shared detection estimate `p`, the occupancy
 log-odds contrast and an LRT p-value. This is a lower-level API; it is not
 silently applied to ordinary one-observation-per-sample count tables.
 """
diff --git a/src/core/doi.jl b/src/core/doi.jl
index 4db53f69..684ecf36 100644
--- a/src/core/doi.jl
+++ b/src/core/doi.jl
@@ -1,6 +1,8 @@
 # © 2026 Joshua Benjamin Jewell. All rights reserved.
 # Licensed under the GNU Affero General Public License version 3 (AGPLv3).
 
+module DOI
+
 # Zenodo deposition-v1 client for automated DOI minting.
 #
 # Publication is DISABLED by default. Opt-in requires DOI_ENABLED=true and a

From 5428506f7ce8474bdd2b8889b1da0306132fffa6 Mon Sep 17 00:00:00 2001
From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com>
Date: Thu, 22 Oct 2026 22:44:26 +0000
Subject: [PATCH 4/6] docs: clarify occupancy delivery and repair DOI route
 syntax

Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com>
---
 EXPLAINME.adoc           | 52 ++++++++++++++++++++++++++++++++++++++++
 README.md                |  9 ++++---
 ROADMAP.md               | 37 ++++++++++++++++++++++++++++
 docs/occupancy-models.md |  6 +++++
 src/server/routes/doi.jl |  4 ++--
 5 files changed, 103 insertions(+), 5 deletions(-)
 create mode 100644 EXPLAINME.adoc
 create mode 100644 ROADMAP.md

diff --git a/EXPLAINME.adoc b/EXPLAINME.adoc
new file mode 100644
index 00000000..47e051f1
--- /dev/null
+++ b/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.
diff --git a/README.md b/README.md
index 482ee321..63b5790b 100644
--- a/README.md
+++ b/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**
 
@@ -273,8 +273,11 @@ process fold change and the zero-process coefficient separately. The public
 `MetaManifold.Differential.fit_occupancy` lower-level estimator is available for
 explicit repeated-visit detection histories (sites × visits × taxa). It refuses
 single-visit histories because occupancy and detection are not separately
-identified from one observation without strong external assumptions. See
-[occupancy model notes](docs/occupancy-models.md).
+identified from one observation without strong external assumptions. Occupancy
+is not yet selectable in the standard analysis configuration or UI, and the
+current fitter has no detection covariates; do not use it as a substitute for
+the full occupancy workflow tracked in [issue #205](https://github.com/hyperpolymath/MetaManifold-WebUI/issues/205).
+See [occupancy model notes](docs/occupancy-models.md).
 
 - `clr_welch`: Welch's t-test on centred log-ratios after multiplicative zero
   replacement with `zCompositions::cmultRepl`. The effect is a CLR difference,
diff --git a/ROADMAP.md b/ROADMAP.md
new file mode 100644
index 00000000..bc98c01c
--- /dev/null
+++ b/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.
diff --git a/docs/occupancy-models.md b/docs/occupancy-models.md
index 86e42a8c..46a28661 100644
--- a/docs/occupancy-models.md
+++ b/docs/occupancy-models.md
@@ -1,5 +1,11 @@
 # 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.
diff --git a/src/server/routes/doi.jl b/src/server/routes/doi.jl
index 76321e85..a54672cc 100644
--- a/src/server/routes/doi.jl
+++ b/src/server/routes/doi.jl
@@ -55,7 +55,7 @@ function _doi_publish_body(req)::Dict{String,Any}
     Dict{String,Any}(string(k) => v for (k, v) in body)
 end
 
-@post \"/api/v1/studies/{study}/doi/publish\" function(req, study::String)
+@post "/api/v1/studies/{study}/doi/publish" function(req, study::String)
     study in _study_names() || return json_error(404, "study_not_found",
                                                   "Study '$study' not found")
     _valid_name(study) || return json_error(400, "invalid_study", "Invalid study name")
@@ -87,7 +87,7 @@ end
             published=receipt.published, created_at=receipt.created_at))
 end
 
-@get \"/api/v1/studies/{study}/doi/status\" function(req, study::String)
+@get "/api/v1/studies/{study}/doi/status" function(req, study::String)
     study in _study_names() || return json_error(404, "study_not_found",
                                                   "Study '$study' not found")
     _valid_name(study) || return json_error(400, "invalid_study", "Invalid study name")

From fab7e4ad612cc10456ee6fbe59069688496e94d4 Mon Sep 17 00:00:00 2001
From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com>
Date: Thu, 22 Oct 2026 23:14:58 +0000
Subject: [PATCH 5/6] fix(stats): canonicalize records and align numeric policy
 tests

Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com>
---
 src/canonical.jl                 | 7 ++++++-
 test/unit/test_high_precision.jl | 2 +-
 test/unit/test_numeric_policy.jl | 2 +-
 3 files changed, 8 insertions(+), 3 deletions(-)

diff --git a/src/canonical.jl b/src/canonical.jl
index b287b7fb..3da806a1 100644
--- a/src/canonical.jl
+++ b/src/canonical.jl
@@ -40,7 +40,7 @@ const SAFE_INT_MAX = 2^53 - 1
     jcs_bytes(x) -> Vector{UInt8}
 
 The RFC 8785 canonical form of `x`, as bytes. Accepts the JSON-shaped
-values the server produces: `Dict` and `JSON3.Object` (key order is
+values the server produces: `Dict`, `NamedTuple` and `JSON3.Object` (key order is
 irrelevant; it is sorted away), `AbstractVector` and `JSON3.Array` (order
 is preserved, as JCS preserves array order), `AbstractString`, `Integer`,
 `Float16`, `Float32`, `Float64`, `Bool` and `nothing`. Floats serialize per ECMAScript's
@@ -120,6 +120,11 @@ function _write_char!(io, c::Char)
     end
 end
 
+# NamedTuple is convenient for internal records. Convert to string-keyed
+# entries so it follows the same canonical sorting and duplicate checks as
+# dictionary objects.
+_write!(io, nt::NamedTuple) = _write!(io, Dict{String,Any}(string(k) => v for (k, v) in pairs(nt)))
+
 # AbstractDict (Dict and JSON3.Object alike): sorted keys, then each pair.
 function _write!(io, d::AbstractDict)
     print(io, '{')
diff --git a/test/unit/test_high_precision.jl b/test/unit/test_high_precision.jl
index d5a29933..477f7b9a 100644
--- a/test/unit/test_high_precision.jl
+++ b/test/unit/test_high_precision.jl
@@ -65,7 +65,7 @@ end
         sample_labels=HP_LABELS, feature_labels=HP_FEATURES, policy=policy)
     @test_throws ArgumentError high_precision_summary([1 -2; 3 4];
         sample_labels=HP_LABELS, feature_labels=["x", "y"], policy=policy)
-    @test_throws ArgumentError high_precision_summary(HP_FIXTURE;
+    @test_throws UnsupportedRepresentationError high_precision_summary(HP_FIXTURE;
         sample_labels=HP_LABELS, feature_labels=HP_FEATURES,
         policy=numeric_policy(:exact_counts))
 
diff --git a/test/unit/test_numeric_policy.jl b/test/unit/test_numeric_policy.jl
index 01b74131..164a389d 100644
--- a/test/unit/test_numeric_policy.jl
+++ b/test/unit/test_numeric_policy.jl
@@ -51,7 +51,7 @@ end
     @test_throws ResourceLimitError exact_fraction(1, big(2)^40; policy=tight)
     @test_throws ResourceLimitError exact_fraction(big(2)^40, 1; policy=tight)
     @test_throws ResourceLimitError exact_rational_sum(
-        [big(1)//big(2)^31, big(1)//big(2)^31 + big(1)//3]; policy=tight)
+        [big(1)//big(2)^31, big(1)//(big(2)^31 - 1)]; policy=tight)
 end
 
 @testset "Precision scopes and restoration" begin

From a5dca81f6e0e2498b0761136efdaf1b20317145c Mon Sep 17 00:00:00 2001
From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com>
Date: Thu, 22 Oct 2026 23:36:54 +0000
Subject: [PATCH 6/6] fix(stats): correct nullable config defaults and
 edge-case tests

Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com>
---
 src/server/routes/analysis.jl     | 5 ++++-
 test/unit/test_exact_summaries.jl | 8 ++++----
 test/unit/test_high_precision.jl  | 6 +++---
 3 files changed, 11 insertions(+), 8 deletions(-)

diff --git a/src/server/routes/analysis.jl b/src/server/routes/analysis.jl
index cb865d0d..8bcfc1d9 100644
--- a/src/server/routes/analysis.jl
+++ b/src/server/routes/analysis.jl
@@ -480,7 +480,10 @@ function _alpha_normalisation(study::String;
                               group::Union{String,Nothing}=nothing,
                               config_cache::Union{Dict,Nothing}=nothing)
     resolved = _resolved_run_config(study, run, group, config_cache)
-    val(k, d) = something(get(get(resolved, k, (; value=d)), :value, d), d)
+    val(k, d) = begin
+        value = get(get(resolved, k, (; value=d)), :value, d)
+        isnothing(value) ? d : value
+    end
     (; method = String(val("analysis.alpha.normalisation", "none")),
        depth = Int(val("analysis.alpha.depth", 0)),
        iterations = Int(val("analysis.alpha.rarefaction_iterations", 1)))
diff --git a/test/unit/test_exact_summaries.jl b/test/unit/test_exact_summaries.jl
index 28ae4087..b5bdbe34 100644
--- a/test/unit/test_exact_summaries.jl
+++ b/test/unit/test_exact_summaries.jl
@@ -13,7 +13,7 @@ const EXACT_FEATURES = ["f1", "f2", "f3"]
     @test !s.approximate
     @test s.samples[1].total == 3
     @test s.samples[1].proportions == fill(big(1)//3, 3)
-    @test s.samples[2].proportions == Union{Rational{BigInt},Nothing}[big(1), nothing, nothing]
+    @test s.samples[2].proportions == Union{Rational{BigInt},Nothing}[big(1), big(0), big(0)]
     @test length(s.groups) == 2
     @test s.groups[1].label == "g1"
     @test s.groups[2].counts == BigInt[2, 0, 0]
@@ -27,7 +27,7 @@ const EXACT_FEATURES = ["f1", "f2", "f3"]
     @test storage["input_status"] == "exact"
     @test storage["samples"][1]["proportions"] == ["1/3", "1/3", "1/3"]
     @test storage["samples"][1]["total"] == 3
-    @test storage["samples"][2]["proportions"][2] === nothing
+    @test storage["samples"][2]["proportions"][2] == "0/1"
     @test length(jcs_sha256(storage)) == 64
     display = summary_to_display(s; digits=8)
     @test occursin("exact", display)
@@ -75,8 +75,8 @@ end
     s = exact_summary(Float64.(EXACT_FIXTURE); sample_labels=EXACT_LABELS,
         feature_labels=EXACT_FEATURES)
     @test s.approximate
-    @test all(x.approximate, s.samples)
-    @test all(any(w -> occursin("approximations", w), x.warnings), s.samples)
+    @test all(x -> x.approximate, s.samples)
+    @test all(x -> any(w -> occursin("approximations", w), x.warnings), s.samples)
     @test summary_to_storage(s)["input_status"] == "approximate"
     @test occursin("remain approximate", summary_to_display(s))
     @test_throws UnsupportedRepresentationError exact_summary(Float64.(EXACT_FIXTURE);
diff --git a/test/unit/test_high_precision.jl b/test/unit/test_high_precision.jl
index 477f7b9a..24901109 100644
--- a/test/unit/test_high_precision.jl
+++ b/test/unit/test_high_precision.jl
@@ -3,13 +3,13 @@ using MetaManifold.NumericPolicy
 using MetaManifold.HighPrecision
 using MetaManifold.Canonical
 
-# Independently produced with Python decimal at 140 digits. They are references,
-# not values copied from Julia/MPFR, and are long enough to check a 512-bit run.
+# Independently produced with Python decimal. These are reference values, not
+# values copied from Julia/MPFR; REF_LN2 has enough digits for the 1024-bit test.
 const HP_FIXTURE = [1 2; 1 0; 1 0]
 const HP_LABELS = ["A", "B"]
 const HP_FEATURES = ["f1", "f2", "f3"]
 
-const REF_LN2 = "0.6931471805599453094172321214581765680755001343602552541206800094933936219696947156058633269964186875420014810205706857336855202357581305570326707516350759619307275708283714351903070386238916734711233501153644979552391205"
+const REF_LN2 = "0.6931471805599453094172321214581765680755001343602552541206800094933936219696947156058633269964186875420014810205706857336855202357581305570326707516350759619307275708283714351903070386238916734711233501153644979552391204751726815749320651555247341395258829504530070953263666426541042391578149520437404303855008019441706416715186447128399681717845469570262716310645461502572074024816377733896385506952606683411372738737229289564935470257626520988596932019650586"
 const REF_LN3 = "1.098612288668109691395245236922525704647490557822749451734694333637494293218608966873615754813732088787970029065957865742368004225930519821052801870767277410603162769183381367179373698844360959903742570316795911521145592"
 const REF_H_1234 = "1.279854225833667467182185762149429928131354186921692886392239591382930699742075074568049151649059984163575294301928203177013603961545520253317776056233522504371230806742168835499680746549033237086147384098562106321046048"