@@ -305,6 +305,79 @@ not exclude samples from this analysis.
   multiplicative (CZM) replacement; `replacement_delta` sets the replacement
   fraction and optional `bayesian_alpha` overrides GBM's concentration. The
   offset is not used. See [zero-handling notes](docs/statistics/zero-handling.md).
+- `ancom_bc`: the bias-corrected log-linear model of Lin & Peddada (2020), which
+  estimates within-group sample-specific offsets and a between-group sampling
+  fraction from a 3-component Gaussian mixture fitted by safeguarded SQUAREM.
+  The effect is a bias-corrected log2 fold change; `offset` and the zero
+  replacement settings are not used.
+- `aldex2`: ALDEx2's Dirichlet Monte-Carlo CLR instances (Fernandes et al. 2014),
+  with a Welch test per instance and a per-instance sign-consistency check; the
+  effect is the standardized CLR difference. `mc_samples` sets the number of
+  Dirichlet draws (whole number in [16, 4096], seeded by the study seed).
+- `songbird`: Songbird's baseline-category compositional multinomial regression
+  (Morton et al. 2019; Zhang et al. 2017), fitted in R by damped Newton-Raphson
+  with exact Sherman-Morrison information inversion. It reports CLR-centred
+  differentials, which sum to zero over the fitted taxa. `songbird` and
+  `multinomial` are two implementations of the same model, in R and in Julia:
+  `songbird` reports centred differentials, `multinomial` reports log ratios
+  against one reference taxon and has the Dirichlet-multinomial sibling.
+- `multinomial`: one multinomial logistic model over the whole composition
+  (Songbird's model, Morton et al. 2019), fitted in Julia. The effect is a log2
+  fold change in relative abundance against a reference taxon. `offset` is not
+  used and no zero is replaced: a zero count is part of the likelihood.
+  `mn_penalty` sets an L2 (ridge) penalty on the log-odds (0 is maximum
+  likelihood); `mn_reference_taxon` names the reference (empty picks the most
+  abundant fitted taxon).
+- `dirichlet_multinomial`: the same model with a Dirichlet prior over the
+  composition, whose concentration absorbs between-sample overdispersion (La
+  Rosa et al. 2012). The result reports the concentration `c` and
+  `theta = 1/(1+c)`; `theta` near 0 is the multinomial.
+  `dm_overdispersion_method` picks the estimator: `mle` (default) profiles the
+  likelihood alternately with the coefficients, `mom` takes one
+  method-of-moments estimate.
+
+#### Why a count model rather than a t-test
+
+Sequencing measures a sample of the library, not an absolute quantity: a sample
+of twice the depth has twice every count, and one taxon can only increase at
+another's expense. Comparisons of raw or rarefied abundances therefore carry the
+compositionality problem (Gloor et al. 2017). The multinomial models fit the whole
+composition at once and condition on each sample's total, so each coefficient is
+the log ratio of two taxa. That ratio is unchanged when every count in a sample
+is multiplied by a constant, which is why no size factor and no pseudocount is
+needed. Because a log ratio is the sum of log ratios through any third taxon,
+every reported effect is coherent under compositionality and can be
+re-expressed against any other taxon.
+
+#### Choosing the reference taxon
+
+Every reported effect is anchored on the reference taxon, so the choice matters
+for interpretability and for stability, but not for the fit: the likelihood and
+the log ratio between two non-reference taxa are the same whichever taxon is the
+reference. An empty `mn_reference_taxon` takes the most abundant taxon that
+passes the prevalence filter. An explicit one must be a taxon at the analysed
+rank with reads in at least half the samples, and reads in both groups. Taxa with
+no reads in one group are removed from the composition and listed as `failed`
+with their note rather than fitted: in a joint model a single diverging
+coefficient would otherwise pull the whole fit with it.
+
+#### Overdispersion
+
+The multinomial treats every read as independent, so its variance is the sampling
+variance and nothing else. Real samples vary more than that: the same community,
+sequenced twice, differs by more than resampling reads from one library. The
+Dirichlet-multinomial adds that between-sample variation as a concentration `c`,
+with `theta = 1/(1+c)` the model's intra-sample correlation. When the data are
+overdispersed the multinomial's standard errors are too small, and the
+Dirichlet-multinomial corrects them; when they are not, `c` runs to its upper
+bound, the fit collapses to the multinomial exactly, and the affected rows are
+marked `boundary` with a note saying so. The two models are otherwise the same
+model on the same scale, so their coefficients are directly comparable.
+
+`max_iter`, `tolerance` and `max_seconds` bound one multinomial fit; a fit that
+runs out of time fails with an explicit error instead of holding the request
+open. Runtime and memory at 100, 1000 and 10000 taxa are measured by
+`bench/compositional_scaling/scaling.jl`, which writes `results/scaling.csv`.
 
 ZINB and hurdle models can overfit sparse data; do not select ZINB solely
 because counts are sparse. Compare against NB diagnostics. The methods report
@@ -316,15 +389,19 @@ 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`,
+Across the differential methods, `min_prevalence` (default 0) filters taxa by the
+fraction of samples with nonzero counts. Benjamini–Hochberg adjustment is applied
+to every fitted taxon and cannot be turned off: it 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/compositional.jl` (the two
+multinomial models), `src/analysis/differential.jl`,
 `src/server/routes/analysis.jl`, `bench/compositional_differential/benchmark.jl`,
-and `test/unit/test_differential.jl`.
+and `test/unit/test_differential.jl`. The reference values in the tests come from
+an independent implementation; the special functions are cross-checked against
+mpmath at 40 digits.
 
 ### Figures and report
 
@@ -662,13 +739,20 @@ analysis:
     normalisation: hellinger     # hellinger | none | rarefy | srs
     depth: 0                     # rarefy/srs reads per sample; 0 = smallest non-empty sample
   differential:
-    method: nb_glm               # nb_glm | clr_welch | ancom_bc | aldex2 | songbird
+    method: nb_glm               # nb_glm | clr_welch | zinb | hurdle_nb | ancom_bc
+                                 #   | aldex2 | songbird | multinomial | dirichlet_multinomial
     dispersion: ml               # nb_glm: ml | glmGamPoi
     replacement_method: bayesian_multiplicative
                                  # clr_welch: bayesian_multiplicative | multiplicative_replacement
     replacement_delta: 0.65      # clr_welch: strictly between 0 and 1
     bayesian_alpha: ~            # clr_welch: null (GBM estimate) or a number > 0
     mc_samples: 128              # aldex2: Dirichlet Monte-Carlo instances in [16, 4096]
+    mn_reference_taxon: ~        # multinomial/DM: reference taxon, ~ picks the most abundant
+    mn_penalty: 0.0              # multinomial/DM: L2 (ridge) penalty, 0 = maximum likelihood
+    dm_overdispersion_method: mle  # DM: mle (profiled) | mom (method of moments)
+    max_iter: 500                # multinomial/DM: Newton iteration cap
+    tolerance: 1.0e-8            # multinomial/DM: convergence tolerance
+    max_seconds: 300.0           # multinomial/DM: wall-clock budget for one fit
   nmds:
     max_stress: 0.2              # warn if NMDS stress exceeds this value
 ```
 
@@ -6,7 +6,12 @@ occupancy work. The GitHub Wiki roadmap still needs a synchronized update; see
 
 ## Available now
 
-- `nb_glm` and `clr_welch` remain the existing differential methods.
+- The differential-abundance config offers four per-taxon count models
+  (`nb_glm`, `clr_welch`, `zinb`, `hurdle_nb`) and five compositional models
+  (`ancom_bc`, `aldex2`, `songbird`, `multinomial`, `dirichlet_multinomial`);
+  see [the differential-abundance section](README.md#differential-abundance).
+  `songbird` and `multinomial` are two implementations of the same
+  baseline-category multinomial model.
 - `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.
 
@@ -40,6 +40,11 @@ log-ratio units, not a fold change. Replacement choice and parameters are
 returned in `config`; zero counts, imputed-mass summaries, warnings, and (for
 GBM) concentration information are returned in `diagnostics.zero_replacement`.
 
+The multinomial models (`multinomial`, `dirichlet_multinomial`) replace no
+zero: a zero count is part of the likelihood, and because the model conditions on
+each sample's total every effect is a log ratio against a reference taxon rather
+than an absolute abundance.
+
 Every zero-replacement method introduces assumptions and can bias results,
 especially when zeros are common or sample sizes are small. Do not choose a
 replacement or `delta` to obtain a preferred p-value. Consider the negative
 
@@ -44,7 +44,9 @@ include("canonical.jl")
 include("analysis/clade_cumulus.jl")
 include("analysis/diversity.jl")
 include("analysis/analysis.jl")
-include("analysis/clade_cumulus.jl")
+# The compositional (multinomial / Dirichlet-multinomial) fits are included
+# before differential.jl, which imports CompositionalModels for those methods.
+include("analysis/compositional.jl")
 include("analysis/differential.jl")
 # Opt-in validated descriptive statistics (issue #1). Ordinary workflows above
 # remain unchanged; these modes must be selected explicitly.
 
@@ -18,7 +18,6 @@ module DOI
 #   * Durable write-ahead journal; immutable DOI/provenance receipts
 #   * DOI insertion never rewrites scientific objects
 
-module DOI
 
 using Dates, SHA, HTTP, JSON3, Logging