diff --git a/README.md b/README.md
index 3fde4c68..49f6d86d 100644
--- a/README.md
+++ b/README.md
@@ -293,7 +293,12 @@ 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`; `dispersion` selects maximum-likelihood theta or a
-  cross-taxon `glmGamPoi` trend.
+  cross-taxon `glmGamPoi` trend. Both options Wald-test the group coefficient
+  with the dispersion treated as known, so every `ml` result — and a
+  `glmGamPoi` result whose smaller group has at most 12 samples — carries a
+  calibration caveat in its diagnostics, shown with the result and written
+  into its CSV export. The caveat changes no p-value; the null-simulation
+  evidence and full replication protocol live in `bench/null_calibration`.
 - `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.
diff --git a/bench/null_calibration/.gitignore b/bench/null_calibration/.gitignore
new file mode 100644
index 00000000..088565ad
--- /dev/null
+++ b/bench/null_calibration/.gitignore
@@ -0,0 +1,4 @@
+# Raw per-taxon p-values are large, append-only working state. The seeds in
+# calibration.R regenerate them dataset-for-dataset; committed evidence is
+# the summarised CSVs and REPORT.md under results/.
+results/raw/
diff --git a/bench/null_calibration/README.md b/bench/null_calibration/README.md
new file mode 100644
index 00000000..f9598d2a
--- /dev/null
+++ b/bench/null_calibration/README.md
@@ -0,0 +1,124 @@
+# Null-calibration bench for `nb_glm` (issue #193)
+
+Issue #193 measured that the default `nb_glm` path — `MASS::glm.nb` per taxon
+with a fixed-scale Wald test, `analysis.differential.dispersion: ml` — rejects
+true nulls too often in small studies (nominal 0.05 rejecting at 0.132 with 4
+samples per group and 0.069 with 12). The mechanism is textbook: the Wald test
+treats the per-taxon maximum-likelihood theta as known, ignoring its sampling
+uncertainty, which is large when groups are small.
+
+This bench replicates that finding at the scale #193 asked for, measures the
+candidate remedies on identical datasets, and records what was actually run.
+It is a bench, not a unit test: slow, stochastic, never run by CI.
+
+## What is compared
+
+Every simulated dataset is analysed by every enabled test (paired design — no
+test gets an easier dataset):
+
+| test | what it is | status |
+|---|---|---|
+| `wald_ml` | the production default, through the real `_FIT_R` block extracted from `src/analysis/differential.jl` (`da_use_gp = FALSE`) | production |
+| `wald_trend` | the production `dispersion: glmGamPoi` option on this branch, through the same block (`da_use_gp = TRUE`): dispersion trend plugged into a fixed-theta Wald fit. **Not** a quasi-likelihood test | production option |
+| `lrt` | `MASS::glm.nb` likelihood-ratio test, full vs intercept-only, each fit estimating its own theta; `2 × ΔlogLik` against χ²(1) | bench-only candidate |
+| `glmgampoi_ql` | glmGamPoi's own quasi-likelihood F-test (`glm_gp` + `test_de`, numeric contrast weights) | reference (the calibrated path in #193) |
+| `edger_ql` | edgeR's quasi-likelihood F-test, fed the same offsets | outside reference only, not a dependency |
+
+A test whose package is unavailable is skipped loudly (stdout + provenance
+file), never silently replaced.
+
+## Protocol (mapped to #193's checklist)
+
+- **Group sizes** `n ∈ {3, 4, 6, 8, 12, 20}` per group, equal groups.
+- **Summaries**: rejection rates at α = 0.01 / 0.05 / 0.10 with Wilson 95%
+  intervals, Kolmogorov–Smirnov D and p of the null p-values against
+  Uniform(0, 1), and per-test fit-failure shares (`summarise.R`).
+- **Data shapes** (`--shapes`):
+  - `original` — the #193 generator: taxon means exp(U(log 2, log 500)),
+    dispersion `0.05 + 2/√μ`, library sizes lognormal(0, 0.3);
+  - `constant` — same means, dispersion 0.5 for every taxon;
+  - `sparse` — amplicon-like: means exp(N(log 3, 1.6)), so most taxa are rare
+    and a few dominate; many zeros;
+  - `unequal_lib` — `original` with group B library sizes ×2.5.
+- **Offsets** (`--offsets`): `tss` and `rle`, implemented to mirror
+  `tss_factors` / `rle_factors` in `differential.jl`. An RLE refusal (no taxon
+  present in every sample) is counted, never replaced by TSS.
+- **Power** (`--truth=power --log2fc=L`): 25% of taxa truly change by ±L (log2),
+  the rest are null. Reports raw and BH power on the changed taxa, the BH
+  false discovery rate at 0.05 (averaged per dataset, BH applied within each
+  dataset over its tested family, as production does), and null-taxon
+  calibration within the same cells.
+- **Taxa filter**: the #193 script's filter is kept for fidelity — a taxon
+  needs reads in ≥2 samples and non-constant counts. Taxon parameters are
+  redrawn for every dataset, as in #193. Datasets with <5 testable taxa are
+  skipped and counted.
+- **Seeds**: `seed0 + cell_index × 1000003 + sim`, one per dataset, so cells
+  can run in any order or in parallel processes and reproduce exactly.
+  Default `seed0 = 20261020` (the #193 seed).
+
+### Reading a cell
+
+A test is **flagged anti-conservative** in a cell when the lower bound of the
+Wilson 95% interval for its rejection rate at α = 0.05 exceeds 0.05, and
+**flagged conservative** when the upper bound is below 0.05. The KS statistic
+is supporting evidence, not the criterion — with ~10⁴ p-values it detects
+departures too small to matter at any single α.
+
+## Running it
+
+Native R with the pinned library (the supported route; `Rscript -e
+'renv::restore()'` first):
+
+```sh
+# smoke: 10 sims, one cell, a few seconds
+Rscript bench/null_calibration/calibration.R --mode=quick
+
+# one direct-replication cell
+Rscript bench/null_calibration/calibration.R \
+  --shapes=original --offsets=tss --n-per=4 --sims=200 --taxa=40 \
+  --tests=wald_ml,lrt
+
+# the full #193 grid — 1,000 sims per cell; expect hours, use several
+# processes by splitting --n-per / --shapes across terminals
+Rscript bench/null_calibration/calibration.R --mode=full
+
+# aggregate whatever raw cells exist
+Rscript bench/null_calibration/summarise.R
+```
+
+Outputs append to `results/raw/pvalues_<cell>.csv` (gitignored), one row per
+taxon × test, so an interrupted run keeps finished simulations; each cell also
+writes `provenance_<cell>.txt` with arguments, package versions and
+`sessionInfo()`. `summarise.R` writes `results/calibration_summary.csv` (and
+`power_summary.csv` when power cells exist).
+
+### No native R? The WebR fallback
+
+`webr_run.mjs` runs the same `calibration.R` / `summarise.R` under WebR (R in
+WebAssembly) for environments with Node but no R. `glm.nb` is pure R (MASS
+7.3-65's `negbin.R` makes no `.C`/`.Call`/`.Fortran` call), so sourcing that
+one file provides genuine MASS code; `calibration.R` detects this mode,
+strips the block's `MASS::` prefixes as a single mechanical substitution, and
+records the fact in every provenance file. glmGamPoi and edgeR are compiled
+Bioconductor packages and cannot run this way — their tests skip loudly.
+Setup and usage are in the header of `webr_run.mjs`. Prefer native Rscript
+whenever it exists.
+
+## What has actually been run
+
+`REPORT.md` records each committed run: cells, sizes, environment and results,
+with the corresponding summary CSVs in `results/`. Anything not in `REPORT.md`
+has not been run — in particular, a committed bench is not a result.
+
+## Limitations (read before quoting rates)
+
+- Taxa are simulated independently. Real amplicon tables are compositional and
+  correlated; independence makes BH behaviour here cleaner than in practice.
+- Groups are equal-sized. Unequal groups change both theta estimation and
+  offset behaviour; #193's grid did not cross that factor either.
+- The ≥2-sample presence filter mirrors #193, not production's default
+  `min_prevalence = 0` (under which untestable taxa instead surface as failed
+  rows outside the BH family). Failure shares are reported so the tested
+  family is never a hidden denominator.
+- `original`-shape library sizes are mild (lognormal σ = 0.3). Harsher depth
+  variation is covered only through `unequal_lib`'s group shift.
diff --git a/bench/null_calibration/REPORT.md b/bench/null_calibration/REPORT.md
new file mode 100644
index 00000000..68c1f34a
--- /dev/null
+++ b/bench/null_calibration/REPORT.md
@@ -0,0 +1,112 @@
+# Null-calibration bench report (issue #193)
+
+This report records the one reduced run that has actually been executed.
+Anything not listed here has not been run. The committed summary CSVs in
+`results/` (`calibration_summary.csv`, `power_summary.csv`) are the source of
+truth for every number below; per-cell run records are summarised in
+`results/provenance_summary.csv` (the raw p-value files under `results/raw/`
+are gitignored and reproducible from the recorded seeds and arguments).
+
+## Run identity
+
+- **Date:** 2026-10-22 (UTC, per-cell provenance records).
+- **Executor:** WebR 0.6.0 (R 4.6.0, wasm32) through `webr_run.mjs`, because the
+  sandbox used for this run had Node but no native R. `MASS::glm.nb` was
+  provided by sourcing MASS 7.3-65's pure-R `negbin.R` (the exact MASS version
+  pinned in `renv.lock`); in that mode `calibration.R` strips the `MASS::`
+  prefixes from the extracted production block as a single mechanical
+  substitution and records the fact in every provenance file. Native
+  `Rscript` runs use the extracted block byte-identically. The sandbox R
+  version (4.6.0) differs from the lockfile's R 4.5.0.
+- **Production path under test:** the real `_FIT_R` block extracted from
+  `src/analysis/differential.jl` (`wald_ml`, `dispersion: ml`), paired on
+  identical datasets with a bench-only likelihood-ratio test (`lrt`: full vs
+  intercept-only `glm.nb`, each fit estimating its own theta, `2 × ΔlogLik`
+  against χ²(1)).
+- **Not run in this sandbox:** `wald_trend`, `glmgampoi_ql`, `edger_ql`
+  (glmGamPoi/edgeR are compiled Bioconductor packages and were absent — the
+  provenance records say so). Their code paths in `calibration.R` are
+  unexecuted here and need a native R/Bioconductor run.
+- **Cells:** direct replication (original/TSS null) at 3/4/6/8/12/20 samples
+  per group, 200 sims × 40 taxa; shape/offset variations (constant, sparse,
+  unequal library sizes, original/RLE) at 4 and 12 per group, 150 × 40; power
+  cells (original/TSS, 25% of taxa at |log2FC| = 1) at 4 and 12 per group,
+  150 × 40. This is a *reduced* run: the full 1,000-sim grid in `README.md`
+  has not been run.
+
+## Results
+
+A cell is flagged anti-conservative when the lower bound of the Wilson 95%
+interval for its α = 0.05 rejection rate exceeds 0.05. **All 32 summarised
+α = 0.05 test-cells were flagged**, including 20 samples per group; the
+closest cell was LRT at 20/group on original/TSS (lower bound 0.0534).
+
+### Direct replication (original shape, TSS offset, null)
+
+Rejection rates at α = 0.01 / 0.05 / 0.10. Issue #193 reported Wald at
+4/group: 0.060 / 0.132 / 0.198 and 12/group: 0.019 / 0.069 / 0.124 — the
+replication lands essentially on those values.
+
+| n/group | sims | `wald_ml` (production) | `lrt` (bench-only) | Wald fit failures | LRT fit failures |
+|---:|---:|---|---|---:|---:|
+| 3 | 200 | 0.092 / 0.179 / 0.251 | 0.046 / 0.131 / 0.203 | 6.41% | 8.43% |
+| 4 | 200 | 0.063 / 0.138 / 0.204 | 0.030 / 0.103 / 0.169 | 2.68% | 3.61% |
+| 6 | 200 | 0.045 / 0.111 / 0.171 | 0.023 / 0.087 / 0.149 | 0.79% | 1.20% |
+| 8 | 200 | 0.029 / 0.089 / 0.145 | 0.017 / 0.070 / 0.129 | 0.20% | 0.33% |
+| 12 | 200 | 0.025 / 0.076 / 0.130 | 0.017 / 0.065 / 0.119 | 0.01% | 0.04% |
+| 20 | 200 | 0.019 / 0.065 / 0.116 | 0.014 / 0.058 / 0.110 | 0.00% | 0.00% |
+
+### Other shapes and the RLE offset (α = 0.05, Wilson 95% interval)
+
+| shape / offset | n/group | sims | `wald_ml` | `lrt` | Wald failures | LRT failures |
+|---|---:|---:|---|---|---:|---:|
+| `constant` / `tss` | 4 | 150 | 0.135 (0.126–0.144) | 0.104 (0.096–0.112) | 4.03% | 5.35% |
+| `constant` / `tss` | 12 | 150 | 0.073 (0.066–0.080) | 0.065 (0.059–0.072) | 0.18% | 0.23% |
+| `sparse` / `tss` | 4 | 150 | 0.107 (0.098–0.116) | 0.082 (0.074–0.090) | 15.12% | 18.36% |
+| `sparse` / `tss` | 12 | 150 | 0.068 (0.062–0.075) | 0.060 (0.054–0.067) | 4.52% | 5.07% |
+| `unequal_lib` / `tss` | 4 | 150 | 0.149 (0.140–0.158) | 0.114 (0.106–0.122) | 1.68% | 3.17% |
+| `unequal_lib` / `tss` | 12 | 150 | 0.074 (0.068–0.081) | 0.062 (0.056–0.068) | 0.00% | 0.08% |
+| `original` / `rle` | 4 | 150 | 0.138 (0.130–0.147) | 0.100 (0.093–0.108) | 2.57% | 3.37% |
+| `original` / `rle` | 12 | 150 | 0.080 (0.074–0.087) | 0.068 (0.061–0.074) | 0.02% | 0.02% |
+
+The inflation is not specific to the original mean-varying dispersion shape,
+and the RLE offset does not rescue calibration.
+
+### Power and BH behaviour (25% of taxa at |log2FC| = 1)
+
+BH is applied within each simulated dataset over its tested family, as
+production does. “Null rate” is the α = 0.05 rejection rate among the truly
+null taxa in the same cells.
+
+| n/group | test | sims | null rate @0.05 | raw power | BH power | BH FDR | mean discoveries |
+|---:|---|---:|---:|---:|---:|---:|---:|
+| 4 | `lrt` | 150 | 0.114 | 0.414 | 0.147 | 0.184 | 1.89 |
+| 4 | `wald_ml` | 150 | 0.150 | 0.480 | 0.290 | 0.379 | 4.70 |
+| 12 | `lrt` | 150 | 0.083 | 0.668 | 0.422 | 0.101 | 4.82 |
+| 12 | `wald_ml` | 150 | 0.096 | 0.688 | 0.509 | 0.129 | 5.99 |
+
+## Interpretation
+
+- #193 replicates. The production fixed-scale Wald test on per-taxon ML theta
+  rejects true nulls well above nominal at every measured group size,
+  including 20/group (0.065 at α = 0.05).
+- The LRT is consistently closer to nominal but remains anti-conservative in
+  every cell, and fails more often than Wald at tiny n (it needs two fits).
+  It is a mitigation, not a calibration fix.
+- Wald's extra BH “power” is bought with uncontrolled FDR (0.38 at 4/group
+  against a nominal 0.05). Nothing measured in this sandbox “stays
+  calibrated”; the calibrated comparator from #193 (glmGamPoi’s
+  quasi-likelihood F-test) still needs the native run.
+- Because even 20/group remained flagged, the production calibration caveat
+  (see `_calibration_caveat` in `src/analysis/differential.jl`) uses “caution”
+  severity through the largest measured size and only drops to “info” beyond
+  the measured range. That caveat is fork-side only; any change to upstream
+  defaults remains an upstream decision (see `UPSTREAM_DRAFT.md`).
+
+## Checklist mapping (#193)
+
+Done in this reduced run: direct replication at larger sim counts and more
+group sizes; additional shapes; RLE offset at 4/12 per group; an additional
+test (LRT) on identical datasets; power and BH-FDR at |log2FC| = 1.
+Not done: the full 1,000-sim grid; `wald_trend`; glmGamPoi QL F; edgeR QL;
+native Rscript/Bioconductor execution; power at other effect/group sizes.
diff --git a/bench/null_calibration/UPSTREAM_DRAFT.md b/bench/null_calibration/UPSTREAM_DRAFT.md
new file mode 100644
index 00000000..5de5b181
--- /dev/null
+++ b/bench/null_calibration/UPSTREAM_DRAFT.md
@@ -0,0 +1,51 @@
+# Draft upstream report: null calibration of the `glm.nb` Wald test
+
+*Status: draft prepared downstream for review before anything is filed
+upstream. All numbers are from `REPORT.md` in this directory, which defines
+exactly what ran, where, and what did not run.*
+
+## Summary
+
+In null simulations that drive the production negative-binomial fitting code
+(`MASS::glm.nb` per taxon, fixed-scale Wald test, per-taxon ML theta treated
+as known), true nulls are rejected more often than nominal at every measured
+group size. At α = 0.05 the rejection rate was 0.138 with 4 samples per group
+and 0.076 with 12 (original/TSS shape, 200 sims × 40 taxa), and still 0.065
+at 20 per group. A likelihood-ratio test on the same fits/datasets was closer
+to nominal (0.103 / 0.065 / 0.058) but still anti-conservative in every cell
+(Wilson 95% lower bound > 0.05 in all 32 summarised test-cells). Under a
+power design (25% of taxa at |log2FC| = 1), Benjamini–Hochberg at 0.05
+delivered FDR 0.38 (Wald, 4/group) and 0.13 (12/group).
+
+The mechanism is the standard one: the Wald test plugs in the per-taxon ML
+theta and ignores its sampling uncertainty. The bench that produces these
+numbers extracts the production fitting block directly from
+`src/analysis/differential.jl`, so the `wald_ml` column is the shipped code
+path, not a reimplementation.
+
+## What was and was not measured
+
+Measured (reduced WebR run, genuine MASS 7.3-65 `glm.nb`; see `REPORT.md` for
+provenance limits): Wald vs LRT, paired on identical datasets; group sizes
+3–20; constant/sparse/unequal-library shapes; RLE offsets; one power design.
+
+Not yet measured: the glmGamPoi quasi-likelihood F-test as an in-bench
+comparator on this branch (the calibrated path reported in the downstream
+issue), the dispersion-trend Wald option, edgeR, and the full 1,000-sim grid.
+The bench (`bench/null_calibration/`, native `Rscript` route) is ready to run
+those where R + Bioconductor are available.
+
+## Options for the upstream maintainer's consideration
+
+1. **Result-notes caveat for small samples** (adopted downstream): record a
+   calibration caveat with `ml`-dispersion results, scaled by group size.
+   Changes no statistics or defaults.
+2. **LRT instead of Wald for the `glm.nb` path:** measurably closer to
+   nominal in this bench, but still anti-conservative, slower, and fails
+   more often at tiny n. The bench quantifies the trade-off; it is not a
+   calibration fix.
+3. **Prefer a quasi-likelihood path below a group-size threshold:** the
+   downstream measurements that motivated this bench found the QL F-test at
+   nominal rates where Wald was inflated. Any such default change is
+   upstream's call; this draft only supplies the replication bench and
+   numbers to decide with.
diff --git a/bench/null_calibration/calibration.R b/bench/null_calibration/calibration.R
new file mode 100644
index 00000000..a424dc4c
--- /dev/null
+++ b/bench/null_calibration/calibration.R
@@ -0,0 +1,395 @@
+#!/usr/bin/env Rscript
+# SPDX-License-Identifier: AGPL-3.0-only
+# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) <j.d.a.jewell@open.ac.uk>
+#
+# Null-calibration and power bench for the nb_glm tests in
+# src/analysis/differential.jl (fork issue #193). This is a bench, not a unit
+# test: it is slow and stochastic, and it never runs in CI by default.
+#
+# What it measures, on the SAME simulated datasets (paired comparison):
+#   wald_ml      - the production default: the real `_FIT_R` block extracted
+#                  from differential.jl, da_use_gp = FALSE (MASS::glm.nb per
+#                  taxon, Wald test at fixed scale 1).
+#   wald_trend   - the production glmGamPoi option ON THIS BRANCH: the same
+#                  `_FIT_R` block with da_use_gp = TRUE (glm_gp dispersion
+#                  trend plugged into a fixed-theta Wald fit). This is NOT
+#                  glmGamPoi's quasi-likelihood test. Needs glmGamPoi.
+#   lrt          - MASS::glm.nb likelihood-ratio test, full (~ group) against
+#                  intercept-only, each fit estimating its own theta;
+#                  2 * (logLik(full) - logLik(null)) against chisq(1).
+#                  Bench-only candidate fix; production does not use it.
+#   glmgampoi_ql - glmGamPoi's own quasi-likelihood F-test (glm_gp + test_de),
+#                  as on the upstream PR 62 branch. Needs glmGamPoi.
+#   edger_ql     - edgeR's quasi-likelihood F-test as an outside reference
+#                  only. Needs edgeR. Not a dependency of the app.
+#
+# Tests whose package is unavailable are skipped LOUDLY (a line on stdout and
+# in the provenance file), never silently replaced by another test.
+#
+# Usage (from the repo root, with the pinned R library restored):
+#   Rscript bench/null_calibration/calibration.R --mode=quick
+#   Rscript bench/null_calibration/calibration.R --mode=full     # hours; see README
+#   Rscript bench/null_calibration/calibration.R --truth=power --log2fc=1 \
+#       --n-per=4,12 --sims=200
+# Any --key=value below overrides the mode defaults. Under WebR (no Rscript
+# argv), settings may instead come from the NB_CALIB_ARGS environment
+# variable as one space-separated string; see webr_run.mjs.
+#
+# Output: one appended CSV per cell under --out (default
+# bench/null_calibration/results/raw), so an interrupted long run keeps the
+# simulations it finished. summarise.R aggregates them.
+
+## ---------------------------------------------------------------- settings
+
+parse_args <- function(argv) {
+  out <- list()
+  for (a in argv) {
+    if (!grepl("^--[a-z0-9-]+=", a)) stop("unrecognised argument '", a, "' (want --key=value)")
+    kv <- strsplit(sub("^--", "", a), "=", fixed = TRUE)[[1]]
+    out[[kv[1]]] <- paste(kv[-1], collapse = "=")
+  }
+  out
+}
+
+argv <- commandArgs(trailingOnly = TRUE)
+if (length(argv) == 0 && nzchar(Sys.getenv("NB_CALIB_ARGS"))) {
+  argv <- strsplit(trimws(Sys.getenv("NB_CALIB_ARGS")), "\\s+")[[1]]
+}
+args <- parse_args(argv)
+
+mode <- if (!is.null(args$mode)) args$mode else "quick"
+defaults <- switch(mode,
+  quick = list(`n-per` = "4", sims = "10", taxa = "20", shapes = "original",
+               offsets = "tss", tests = "wald_ml,lrt", truth = "null",
+               log2fc = "1", seed = "20261020"),
+  full  = list(`n-per` = "3,4,6,8,12,20", sims = "1000", taxa = "40",
+               shapes = "original,constant,sparse,unequal_lib",
+               offsets = "tss,rle",
+               tests = "wald_ml,wald_trend,lrt,glmgampoi_ql,edger_ql",
+               truth = "null", log2fc = "1", seed = "20261020"),
+  stop("unknown --mode '", mode, "' (want quick or full, plus overrides)")
+)
+for (k in names(defaults)) if (is.null(args[[k]])) args[[k]] <- defaults[[k]]
+if (is.null(args$out)) args$out <- file.path("bench", "null_calibration", "results", "raw")
+if (is.null(args[["differential-jl"]])) {
+  args[["differential-jl"]] <- file.path("src", "analysis", "differential.jl")
+}
+
+csv_ints <- function(s) as.integer(strsplit(s, ",", fixed = TRUE)[[1]])
+csv_strs <- function(s) trimws(strsplit(s, ",", fixed = TRUE)[[1]])
+
+n_pers   <- csv_ints(args[["n-per"]])
+n_sims   <- as.integer(args$sims)
+n_taxa   <- as.integer(args$taxa)
+shapes   <- csv_strs(args$shapes)
+offsets  <- csv_strs(args$offsets)
+want     <- csv_strs(args$tests)
+truth    <- args$truth
+lfc      <- as.numeric(args$log2fc)
+seed0    <- as.integer(args$seed)
+out_dir  <- args$out
+stopifnot(truth %in% c("null", "power"))
+dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)
+
+## ------------------------------------------------------------ MASS access
+# Normally MASS is installed and the extracted production block runs
+# byte-identical. In a minimal R with no MASS but with MASS's pure-R
+# negbin.R already sourced (glm.nb is pure R; MASS 7.3-65 negbin.R makes no
+# .C/.Call/.Fortran call), the block's `MASS::` prefixes are stripped as a
+# single mechanical substitution and that fact is recorded in the
+# provenance. The functions executed are still MASS's own source.
+
+have_mass_ns <- requireNamespace("MASS", quietly = TRUE)
+mass_mode <- if (have_mass_ns) {
+  paste0("installed MASS ", as.character(packageVersion("MASS")))
+} else if (exists("glm.nb", envir = globalenv())) {
+  "sourced negbin.R (MASS:: stripped from the extracted block; see header)"
+} else {
+  stop("MASS is not available and glm.nb is not sourced: restore renv.lock first")
+}
+glm_nb <- if (have_mass_ns) MASS::glm.nb else get("glm.nb", envir = globalenv())
+
+have_gp <- requireNamespace("glmGamPoi", quietly = TRUE)
+have_edger <- requireNamespace("edgeR", quietly = TRUE)
+
+## ------------------------------------------------------- production block
+
+extract_fit_r <- function(path) {
+  lines <- readLines(path, warn = FALSE)
+  start <- grep('const _FIT_R = raw"""', lines, fixed = TRUE)
+  if (length(start) != 1) stop("could not find exactly one _FIT_R block in ", path)
+  rest <- lines[(start + 1):length(lines)]
+  end <- which(rest == '"""')[1]
+  if (is.na(end)) stop("_FIT_R block in ", path, " is not terminated")
+  block <- paste(rest[seq_len(end - 1)], collapse = "\n")
+  if (!grepl("da_result", block, fixed = TRUE) || !grepl("glm.nb", block, fixed = TRUE)) {
+    stop("extracted _FIT_R does not look like the production fit loop")
+  }
+  block
+}
+
+fit_code <- extract_fit_r(args[["differential-jl"]])
+if (!have_mass_ns) fit_code <- gsub("MASS::", "", fit_code, fixed = TRUE)
+
+# Every da_* global fit_nb sets or the block creates on the current branch
+# (mirrors the cleanup list in differential.jl's fit_nb).
+DA_VARS <- c("da_counts", "da_group", "da_levels", "da_offset", "da_theta_upper",
+             "da_theta_lower", "da_use_gp", "da_g", "da_term", "da_n", "da_result",
+             "da_phi", "da_gp_failed", "da_gp", "da_j", "da_y", "da_warns",
+             "da_fit", "da_co", "da_row", "da_theta", "da_absent")
+
+run_block <- function(counts, offset, use_gp, n_per) {
+  g <- globalenv()
+  assign("da_counts", counts, g)
+  assign("da_group", rep(c("A", "B"), each = n_per), g)
+  assign("da_levels", c("A", "B"), g)
+  assign("da_offset", offset, g)
+  assign("da_theta_upper", 1e7, g)
+  assign("da_theta_lower", 1e-8, g)
+  assign("da_use_gp", use_gp, g)
+  on.exit(rm(list = intersect(ls(g), DA_VARS), envir = g))
+  suppressWarnings(eval(parse(text = fit_code), envir = g))
+  res <- get("da_result", envir = g)
+  res
+}
+
+## ---------------------------------------------------------------- offsets
+# Mirror src/analysis/differential.jl exactly: factors centred to geometric
+# mean 1, offset = log(factor). RLE refuses (returns NULL) when no taxon is
+# present in every sample, as rle_factors does; the refusal is recorded,
+# never replaced by TSS.
+
+geomean <- function(x) exp(mean(log(x)))
+
+make_offset <- function(counts, method) {
+  if (method == "tss") {
+    lib <- rowSums(counts)
+    if (any(lib <= 0)) return(NULL)
+    return(log(lib / geomean(lib)))
+  }
+  if (method == "rle") {
+    shared <- which(colSums(counts > 0) == nrow(counts))
+    if (length(shared) == 0) return(NULL)
+    sub <- counts[, shared, drop = FALSE]
+    gm <- apply(sub, 2, geomean)
+    raw <- apply(sub, 1, function(row) median(row / gm))
+    return(log(raw / geomean(raw)))
+  }
+  stop("unknown offset method '", method, "'")
+}
+
+## ------------------------------------------------------------- simulation
+# One dataset: samples x taxa integer counts plus the true log2 fold change
+# per taxon (all zero under the null). Taxon parameters are redrawn per
+# dataset, as in the #193 script. The data-generating library sizes are
+# realisations, not observed library sizes: offsets are computed from the
+# simulated counts afterwards, as production computes them from data.
+
+simulate <- function(shape, n_per, n_taxa, truth, lfc) {
+  n <- 2 * n_per
+  mu <- switch(shape,
+    original    = exp(runif(n_taxa, log(2), log(500))),
+    constant    = exp(runif(n_taxa, log(2), log(500))),
+    sparse      = exp(rnorm(n_taxa, mean = log(3), sd = 1.6)),
+    unequal_lib = exp(runif(n_taxa, log(2), log(500))),
+    stop("unknown shape '", shape, "'")
+  )
+  disp <- if (shape == "constant") rep(0.5, n_taxa) else 0.05 + 2 / sqrt(mu)
+  lib <- exp(rnorm(n, 0, 0.3))
+  if (shape == "unequal_lib") lib[(n_per + 1):n] <- lib[(n_per + 1):n] * 2.5
+  true_lfc <- rep(0, n_taxa)
+  if (truth == "power") {
+    n_de <- max(1L, round(0.25 * n_taxa))
+    de <- seq_len(n_de)
+    true_lfc[de] <- rep(c(lfc, -lfc), length.out = n_de)
+  }
+  y <- vapply(seq_len(n_taxa), function(j) {
+    m <- mu[j] * lib
+    if (true_lfc[j] != 0) m[(n_per + 1):n] <- m[(n_per + 1):n] * 2^true_lfc[j]
+    as.integer(rnbinom(n, mu = m, size = 1 / disp[j]))
+  }, integer(n))
+  storage.mode(y) <- "integer"
+  # The #193 script's filter, kept so the replication is faithful: a taxon
+  # needs reads in at least 2 samples and non-constant counts to be tested.
+  keep <- colSums(y > 0) >= 2 & apply(y, 2, function(v) length(unique(v)) > 1)
+  list(counts = y[, keep, drop = FALSE], true_lfc = true_lfc[keep],
+       zero_frac = mean(y == 0))
+}
+
+## ------------------------------------------------------------ bench tests
+
+# Production's per-taxon screens, applied to the bench-only tests so every
+# test faces the same family: constant counts, and complete separation (no
+# reads in any sample of one group). Returns "" when fittable.
+screen <- function(yv, g) {
+  if (length(unique(yv)) < 2L) return("constant counts")
+  if (all(yv[g == "A"] == 0) || all(yv[g == "B"] == 0)) return("complete separation")
+  ""
+}
+
+fit_lrt <- function(counts, offset, n_per) {
+  g <- factor(rep(c("A", "B"), each = n_per), levels = c("A", "B"))
+  out <- data.frame(status = "failed", pvalue = NA_real_)
+  out <- out[rep(1, ncol(counts)), , drop = FALSE]
+  rownames(out) <- NULL
+  one <- function(formula, yv) {
+    warns <- character(0)
+    fit <- tryCatch(
+      withCallingHandlers(glm_nb(formula, control = glm.control(maxit = 100)),
+        warning = function(w) { warns <<- c(warns, conditionMessage(w)); invokeRestart("muffleWarning") }),
+      error = function(e) e)
+    list(fit = fit, warns = warns)
+  }
+  bad_fit <- function(r) {
+    inherits(r$fit, "error") || !isTRUE(r$fit[["converged"]]) ||
+      !is.null(r$fit[["th.warn"]]) ||
+      any(grepl("iteration limit|alternation limit|did not converge|NaNs produced", r$warns))
+  }
+  off <- offset
+  for (j in seq_len(ncol(counts))) {
+    yv <- counts[, j]
+    if (nzchar(screen(yv, g))) next
+    # glm.nb looks its variables up in the formula's environment; yv, g and
+    # off all live in this frame (one()'s enclosing environment), which is
+    # where these formulas are created.
+    full <- one(yv ~ g + offset(off), yv)
+    if (bad_fit(full)) next
+    null <- one(yv ~ 1 + offset(off), yv)
+    if (bad_fit(null)) next
+    lr <- 2 * (as.numeric(logLik(full$fit)) - as.numeric(logLik(null$fit)))
+    if (!is.finite(lr) || lr < 0) next
+    p <- pchisq(lr, df = 1, lower.tail = FALSE)
+    if (!is.finite(p)) next
+    out$status[j] <- "ok"
+    out$pvalue[j] <- p
+  }
+  out
+}
+
+fit_glmgampoi_ql <- function(counts, offset, n_per) {
+  g <- factor(rep(c("A", "B"), each = n_per), levels = c("A", "B"))
+  out <- data.frame(status = rep("failed", ncol(counts)), pvalue = NA_real_)
+  screens <- vapply(seq_len(ncol(counts)),
+                    function(j) nzchar(screen(counts[, j], g)), logical(1))
+  res <- tryCatch(suppressWarnings({
+    gp <- glmGamPoi::glm_gp(t(counts), design = model.matrix(~ g),
+                            size_factors = exp(offset),
+                            overdispersion_shrinkage = TRUE)
+    # Numeric contrast weights on the design columns ((Intercept), gB):
+    # test_de parses a TEXT contrast as R code, so bench and production
+    # alike must never hand it a label.
+    glmGamPoi::test_de(gp, contrast = c(0, 1), sort_by = "none")
+  }), error = function(e) e)
+  if (inherits(res, "error")) return(out)
+  pcol <- intersect(c("pval", "pvalue", "p_value"), names(res))[1]
+  if (is.na(pcol) || length(res[[pcol]]) != ncol(counts)) return(out)
+  p <- as.numeric(res[[pcol]])
+  ok <- !screens & is.finite(p) & p >= 0 & p <= 1
+  out$status[ok] <- "ok"
+  out$pvalue[ok] <- p[ok]
+  out
+}
+
+fit_edger_ql <- function(counts, offset, n_per) {
+  g <- factor(rep(c("A", "B"), each = n_per), levels = c("A", "B"))
+  out <- data.frame(status = rep("failed", ncol(counts)), pvalue = NA_real_)
+  screens <- vapply(seq_len(ncol(counts)),
+                    function(j) nzchar(screen(counts[, j], g)), logical(1))
+  res <- tryCatch(suppressWarnings({
+    y <- edgeR::DGEList(counts = t(counts), group = g)
+    # Feed edgeR the bench offset rather than its own normalisation, so all
+    # tests see the same library-size handling: effective library size
+    # proportional to exp(offset).
+    y$samples$lib.size <- exp(offset) * geomean(rowSums(counts))
+    y$samples$norm.factors <- 1
+    design <- model.matrix(~ g)
+    y <- edgeR::estimateDisp(y, design)
+    fit <- edgeR::glmQLFit(y, design)
+    qlf <- edgeR::glmQLFTest(fit, coef = 2)
+    edgeR::topTags(qlf, n = Inf, sort.by = "none")$table
+  }), error = function(e) e)
+  if (inherits(res, "error")) return(out)
+  pcol <- intersect(c("pvalue", "p.value", "pvalue"), names(res))[1]
+  if (is.na(pcol) || length(res[[pcol]]) != ncol(counts)) return(out)
+  p <- as.numeric(res[[pcol]])
+  ok <- !screens & is.finite(p) & p >= 0 & p <= 1
+  out$status[ok] <- "ok"
+  out$pvalue[ok] <- p[ok]
+  out
+}
+
+## ------------------------------------------------------------------- main
+
+enabled <- want
+for (t in want) {
+  if (t %in% c("wald_trend", "glmgampoi_ql") && !have_gp) {
+    cat("SKIPPED test", t, ": glmGamPoi is not available in this R library\n")
+    enabled <- setdiff(enabled, t)
+  }
+  if (t == "edger_ql" && !have_edger) {
+    cat("SKIPPED test edger_ql: edgeR is not available in this R library\n")
+    enabled <- setdiff(enabled, t)
+  }
+}
+if (length(enabled) == 0) stop("no requested test can run in this R library")
+
+cell_index <- 0L
+for (shape in shapes) for (off_m in offsets) for (n_per in n_pers) {
+  cell_index <- cell_index + 1L
+  # The taxa count is part of the cell identity: reruns append to a cell's
+  # CSV, so two settings must never share a file by accident.
+  cell <- sprintf("%s_%s_%s_n%d_t%d%s", truth, shape, off_m, n_per, n_taxa,
+                  if (truth == "power") sprintf("_lfc%s", lfc) else "")
+  f <- file.path(out_dir, paste0("pvalues_", cell, ".csv"))
+  prov <- file.path(out_dir, paste0("provenance_", cell, ".txt"))
+  writeLines(c(
+    paste("cell:", cell), paste("date:", format(Sys.time(), tz = "UTC", usetz = TRUE)),
+    paste("args:", paste(argv, collapse = " ")),
+    paste("MASS:", mass_mode),
+    paste("glmGamPoi:", if (have_gp) as.character(packageVersion("glmGamPoi")) else "absent"),
+    paste("edgeR:", if (have_edger) as.character(packageVersion("edgeR")) else "absent"),
+    paste("tests run:", paste(enabled, collapse = ",")),
+    capture.output(sessionInfo())
+  ), prov)
+  cat(sprintf("cell %s: %d sims x %d taxa, tests %s\n",
+              cell, n_sims, n_taxa, paste(enabled, collapse = "+")))
+  header_done <- file.exists(f)
+  skipped_sims <- 0L
+  refused_offsets <- 0L
+  for (s in seq_len(n_sims)) {
+    # One seed per (cell, sim): cells can run in any order, in parallel
+    # processes, and reproduce dataset-for-dataset.
+    set.seed(seed0 + cell_index * 1000003L + s)
+    sim <- simulate(shape, n_per, n_taxa, truth, lfc)
+    counts <- sim$counts
+    if (ncol(counts) < 5) { skipped_sims <- skipped_sims + 1L; next }
+    offset <- make_offset(counts, off_m)
+    if (is.null(offset)) { refused_offsets <- refused_offsets + 1L; next }
+    fits <- list()
+    if ("wald_ml" %in% enabled) {
+      r <- run_block(counts, offset, FALSE, n_per)
+      fits$wald_ml <- data.frame(status = r$status, pvalue = r$pvalue)
+    }
+    if ("wald_trend" %in% enabled) {
+      r <- run_block(counts, offset, TRUE, n_per)
+      fits$wald_trend <- data.frame(status = r$status, pvalue = r$pvalue)
+    }
+    if ("lrt" %in% enabled) fits$lrt <- fit_lrt(counts, offset, n_per)
+    if ("glmgampoi_ql" %in% enabled) fits$glmgampoi_ql <- fit_glmgampoi_ql(counts, offset, n_per)
+    if ("edger_ql" %in% enabled) fits$edger_ql <- fit_edger_ql(counts, offset, n_per)
+    rows <- do.call(rbind, lapply(names(fits), function(t) {
+      data.frame(truth = truth, shape = shape, offset = off_m, n_per = n_per,
+                 taxa = n_taxa, sim = s, taxon = seq_len(ncol(counts)),
+                 true_log2fc = sim$true_lfc, test = t,
+                 status = fits[[t]]$status, pvalue = fits[[t]]$pvalue,
+                 stringsAsFactors = FALSE)
+    }))
+    write.table(rows, f, sep = ",", row.names = FALSE,
+                col.names = !header_done, append = header_done, quote = FALSE)
+    header_done <- TRUE
+    if (s %% 25 == 0) { cat(sprintf("  %s sim %d/%d\n", cell, s, n_sims)); flush.console() }
+  }
+  cat(sprintf("cell %s done: %d sims skipped (<5 testable taxa), %d offsets refused\n",
+              cell, skipped_sims, refused_offsets))
+}
+cat("calibration run complete\n")
diff --git a/bench/null_calibration/results/calibration_summary.csv b/bench/null_calibration/results/calibration_summary.csv
new file mode 100644
index 00000000..a73478c8
--- /dev/null
+++ b/bench/null_calibration/results/calibration_summary.csv
@@ -0,0 +1,33 @@
+"truth","shape","offset","n_per_group","n_taxa_setting","test","n_sims","n_attempted","n_pvalues","n_failed","n_boundary","fail_rate","rate_01","rate_01_lo","rate_01_hi","rate_05","rate_05_lo","rate_05_hi","rate_10","rate_10_lo","rate_10_hi","ks_D","ks_p","mean_p_null"
+"null","constant","tss",4,40,"lrt",150,6000,5679,321,0,0.0535,0.0343370311674591,0.0299070144118022,0.0393966009845745,0.103539355520338,0.0958817540765541,0.111732951941575,0.167635147032928,0.158145307258128,0.177574327478367,0.0843830805075691,1.50521819708956e-35,0.444054181891207
+"null","constant","tss",4,40,"wald_ml",150,6000,5758,242,147,0.0403333333333333,0.067210837096214,0.0610277913962478,0.0739709697428515,0.134769017019799,0.126192001748904,0.14383303624288,0.195206668982286,0.185173587708648,0.205646165822299,0.104078202508359,1.3341509762296e-54,0.433539982517614
+"null","constant","tss",12,40,"lrt",150,6000,5986,14,0,0.00233333333333333,0.0143668559973271,0.0116487013673261,0.0177079118301315,0.0649849649181423,0.0590152557624449,0.0715126495578783,0.115603073838958,0.107748388407623,0.123950809035566,0.0376376568717833,8.62249043869314e-08,0.485492376920367
+"null","constant","tss",12,40,"wald_ml",150,6000,5989,11,6,0.00183333333333333,0.0213725162798464,0.0180049117092347,0.0253537288632171,0.0728001335782267,0.066490416983351,0.0796575271661756,0.127233260978461,0.119031959024314,0.135912455794454,0.0445719596399606,9.25695055432286e-11,0.481410606985188
+"null","original","rle",4,40,"lrt",150,5997,5795,202,0,0.0336835084208771,0.0291630716134599,0.0251329702386332,0.0338169874810617,0.100431406384815,0.0929553445016146,0.108436858943485,0.167385677308024,0.15799492693621,0.177217110265444,0.0999852967297406,9.5741484666021e-51,0.437673327317577
+"null","original","rle",4,40,"wald_ml",150,5997,5843,154,79,0.0256795064198766,0.0641793599178504,0.0581774162795938,0.0707539843971647,0.138113982543214,0.129504908831245,0.147198584878162,0.199897313024131,0.189841646759709,0.210347322840407,0.121629858990801,1.65933724634618e-75,0.426544757195022
+"null","original","rle",12,40,"lrt",150,6000,5999,1,0,0.000166666666666667,0.0170028338056343,0.0140269075883063,0.0205969385096592,0.0675112518753126,0.0614348261387838,0.0741412113767197,0.126187697949658,0.118023356072251,0.134830474773197,0.0374388440525218,9.94059353574812e-08,0.476768173372017
+"null","original","rle",12,40,"wald_ml",150,6000,5999,1,0,0.000166666666666667,0.024837472912152,0.021192855359116,0.029090241552791,0.0801800300050008,0.0735734789459446,0.0873239003112516,0.138523087181197,0.130012534390866,0.14749628709781,0.0425443113010688,7.40678934295594e-10,0.472785456890374
+"null","original","tss",3,40,"lrt",200,7988,7315,673,0,0.0842513770655984,0.045933014354067,0.0413694275154052,0.0509732557662729,0.131237183868763,0.12369250293294,0.139168971747875,0.203417634996582,0.19434974613276,0.21279685973194,0.135182796664952,1.55071172931089e-116,0.415587740716814
+"null","original","tss",3,40,"wald_ml",200,7988,7476,512,273,0.0640961442163245,0.0918940609951846,0.0855537288027265,0.0986535791203389,0.179240235420011,0.170711249407119,0.188098689760206,0.250668806848582,0.240974270777598,0.260619443404176,0.161656443920632,4.03457010651865e-170,0.401366722446981
+"null","original","tss",4,40,"lrt",200,7997,7708,289,0,0.0361385519569839,0.0303580695381422,0.026755637175536,0.0344283824052878,0.102880124545926,0.0962945813035816,0.109861298236999,0.169434353918007,0.161224866817917,0.177973166853889,0.0933013909161196,1.0455020069296e-58,0.443086204545911
+"null","original","tss",4,40,"wald_ml",200,7997,7783,214,106,0.0267600350131299,0.0628292432224078,0.0576509835270075,0.0684388391582219,0.13825003212129,0.130760030253546,0.146096956012513,0.203520493383014,0.194723061730016,0.212610447703097,0.113643078400003,9.87146482984532e-88,0.431436851581758
+"null","original","tss",6,40,"lrt",200,8000,7904,96,0,0.012,0.0232793522267206,0.0201794197046358,0.0268424459659169,0.086918016194332,0.0809063372182911,0.09333102780675,0.148658906882591,0.14098681509377,0.156672346521709,0.0645951152356356,4.52120174669666e-29,0.461303882402388
+"null","original","tss",6,40,"wald_ml",200,8000,7937,63,28,0.007875,0.0447272269119315,0.0403957707152817,0.0494991683188881,0.111251102431649,0.104520582486221,0.118357744451632,0.170719415396245,0.162601431411705,0.179155984712773,0.0787526584971869,3.50452954348045e-43,0.453506610955192
+"null","original","tss",8,40,"lrt",200,8000,7974,26,0,0.00325,0.0171808377225984,0.0145524202088425,0.0202742256019111,0.0701028342111863,0.0647034005936096,0.0759162726113939,0.128668171557562,0.121497414397098,0.136196533199172,0.0453296175635426,1.17327512113714e-14,0.473696316973658
+"null","original","tss",8,40,"wald_ml",200,8000,7984,16,6,0.002,0.0293086172344689,0.0258291673085302,0.0332407906085308,0.0886773547094188,0.0826378982356264,0.0951124322059991,0.145415831663327,0.137853814201029,0.153318897577363,0.0552258902999433,1.41431816022389e-21,0.467708032138111
+"null","original","tss",12,40,"lrt",200,8000,7997,3,0,0.000375,0.0171314242841065,0.0145105151972706,0.0202160145372216,0.0645241965737151,0.059345808818977,0.0701207559357476,0.119294735525822,0.112372755059871,0.126582293438552,0.0329734074224069,5.6092770626337e-08,0.480625097706312
+"null","original","tss",12,40,"wald_ml",200,8000,7999,1,1,0.000125,0.0251281410176272,0.0219194190982992,0.028792751186341,0.0757594699337417,0.0701620857850788,0.0817641350549703,0.130016252031504,0.122823186861121,0.137564510376576,0.0381036842306565,1.63501153735572e-10,0.476689665536424
+"null","original","tss",20,40,"lrt",200,8000,8000,0,0,0,0.014125,0.0117624417845792,0.0169539514627402,0.058375,0.0534462764601491,0.0637276410450682,0.11025,0.103572963173447,0.117301159323374,0.020738749813529,0.00205314046936761,0.487703822259244
+"null","original","tss",20,40,"wald_ml",200,8000,8000,0,0,0,0.018625,0.0158851798715362,0.0218268938087898,0.0645,0.0593234953117845,0.0700945427828028,0.116,0.109166456166853,0.123202146883254,0.026145410693358,3.55641366754023e-05,0.48530796517549
+"null","sparse","tss",4,40,"lrt",150,5365,4380,985,0,0.183597390493942,0.0237442922374429,0.0196354002458669,0.0286878479396827,0.0817351598173516,0.0739836169914801,0.0902197342433841,0.144063926940639,0.133976279706903,0.154775371262067,0.0691867976694867,1.2302251751396e-18,0.455338729993458
+"null","sparse","tss",4,40,"wald_ml",150,5365,4554,811,236,0.15116495806151,0.0408432147562582,0.0354710856191992,0.0469893207412142,0.106719367588933,0.0980810744949493,0.116020593506725,0.170838823012736,0.160186222147154,0.182046273852014,0.0823244156861178,3.1118285502842e-27,0.44737008911558
+"null","sparse","tss",12,40,"lrt",150,5859,5562,297,0,0.0506912442396313,0.015282272563826,0.0123767061844007,0.0188569284972326,0.0604099244875944,0.054446958150994,0.0669796875450717,0.110931319669184,0.102945034249156,0.119454663618446,0.0326604118895657,1.40505153549206e-05,0.480162876414684
+"null","sparse","tss",12,40,"wald_ml",150,5859,5594,265,33,0.0452295613585936,0.0196639256346085,0.0163415019144068,0.0236456003405133,0.0682874508401859,0.0619693876083744,0.075198030309613,0.119413657490168,0.111176063604179,0.128173597988946,0.0345742276076426,3.11048126202616e-06,0.476462685961144
+"null","unequal_lib","tss",4,40,"lrt",150,5999,5809,190,0,0.0316719453242207,0.0395937338612498,0.0348758076860468,0.0449201857375891,0.113616801514891,0.105710109456138,0.122034181902751,0.17868824238251,0.169050117338897,0.188751049879355,0.0933208396676188,2.28941310680147e-44,0.440256995667013
+"null","unequal_lib","tss",4,40,"wald_ml",150,5999,5898,101,56,0.0168361393565594,0.0764665988470668,0.0699568766474143,0.0835276696829271,0.14852492370295,0.139678040812683,0.157829350917821,0.211766700576467,0.20152919358122,0.222379424818132,0.115545521434757,8.05155216532249e-69,0.429486140222383
+"null","unequal_lib","tss",12,40,"lrt",150,6000,5995,5,0,0.000833333333333333,0.0128440366972477,0.0102895388904349,0.0160224515004104,0.0618849040867389,0.0560617457461293,0.0682691711384984,0.116263552960801,0.108394130866014,0.124624439205597,0.023602226833183,0.00251353912720786,0.48526062352396
+"null","unequal_lib","tss",12,40,"wald_ml",150,6000,6000,0,1,0,0.0196666666666667,0.0164483267012244,0.0234996733360393,0.0741666666666667,0.0678051993175382,0.0810730588695571,0.13,0.121726673024468,0.138746803756492,0.0305808238880169,2.67468380464776e-05,0.481384434000172
+"power","original","tss",4,40,"lrt",150,5997,5768,229,0,0.0381857595464399,0.0325109522711552,0.0276325349035269,0.038216785218464,0.114134194143417,0.105010338234011,0.123940999484436,0.183306433018215,0.172073169738517,0.195100216790487,0.105618632297901,1.89778614361325e-42,0.434784684139481
+"power","original","tss",4,40,"wald_ml",150,5997,5841,156,88,0.0260130065032516,0.0721602549510585,0.0648766537321604,0.0801914513597557,0.150466651491008,0.140199738589643,0.161344329177655,0.218757113589802,0.206780906784307,0.23122475633943,0.127157614018003,4.02317867444891e-62,0.422783825089475
+"power","original","tss",12,40,"lrt",150,6000,5997,3,0,5e-04,0.0237777777777778,0.0197159062402261,0.0286520171896841,0.0831111111111111,0.0753968197212776,0.0915365560857739,0.150888888888889,0.140728778500017,0.161644533543095,0.0677456389208428,2.30351108889158e-18,0.456624861778725
+"power","original","tss",12,40,"wald_ml",150,6000,5998,2,1,0.000333333333333333,0.0328888888888889,0.0280635063862089,0.0385110969967455,0.096,0.0877341616927958,0.104955005267343,0.16,0.1495793559677,0.171000636028456,0.0737015111182361,1.17376844197725e-21,0.452433541509862
diff --git a/bench/null_calibration/results/power_summary.csv b/bench/null_calibration/results/power_summary.csv
new file mode 100644
index 00000000..b9685f56
--- /dev/null
+++ b/bench/null_calibration/results/power_summary.csv
@@ -0,0 +1,5 @@
+"shape","offset","n_per_group","test","n_sims","power_raw_05","power_bh_05","bh_fdr_05","mean_discoveries"
+"original","tss",4,"lrt",150,0.413696715583508,0.147449336128581,0.183793650793651,1.89333333333333
+"original","tss",4,"wald_ml",150,0.479972375690608,0.290055248618785,0.379349206349206,4.7
+"original","tss",12,"lrt",150,0.668002672010688,0.421509686038744,0.100953102453102,4.82
+"original","tss",12,"wald_ml",150,0.688251001335114,0.509345794392523,0.128748011248011,5.99333333333333
diff --git a/bench/null_calibration/results/provenance_summary.csv b/bench/null_calibration/results/provenance_summary.csv
new file mode 100644
index 00000000..5cc78144
--- /dev/null
+++ b/bench/null_calibration/results/provenance_summary.csv
@@ -0,0 +1,17 @@
+cell,date_utc,args,mass_mode,glmGamPoi,edgeR,tests_run,r_version,platform,webr
+null_constant_tss_n12_t40,2026-10-22 22:42:23 UTC,"--shapes=constant --offsets=tss --n-per=12 --sims=150 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0
+null_constant_tss_n4_t40,2026-10-22 22:43:30 UTC,"--shapes=constant --offsets=tss --n-per=4 --sims=150 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0
+null_original_rle_n12_t40,2026-10-22 22:54:06 UTC,"--shapes=original --offsets=rle --n-per=12 --sims=150 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0
+null_original_rle_n4_t40,2026-10-22 22:55:33 UTC,"--shapes=original --offsets=rle --n-per=4 --sims=150 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0
+null_original_tss_n12_t40,2026-10-22 22:32:15 UTC,"--shapes=original --offsets=tss --n-per=12 --sims=200 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0
+null_original_tss_n20_t40,2026-10-22 22:37:21 UTC,"--shapes=original --offsets=tss --n-per=20 --sims=200 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0
+null_original_tss_n3_t40,2026-10-22 22:27:17 UTC,"--shapes=original --offsets=tss --n-per=3 --sims=200 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0
+null_original_tss_n4_t40,2026-10-22 22:32:41 UTC,"--shapes=original --offsets=tss --n-per=4 --sims=200 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0
+null_original_tss_n6_t40,2026-10-22 22:38:09 UTC,"--shapes=original --offsets=tss --n-per=6 --sims=200 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0
+null_original_tss_n8_t40,2026-10-22 22:27:17 UTC,"--shapes=original --offsets=tss --n-per=8 --sims=200 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0
+null_sparse_tss_n12_t40,2026-10-22 22:46:13 UTC,"--shapes=sparse --offsets=tss --n-per=12 --sims=150 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0
+null_sparse_tss_n4_t40,2026-10-22 22:47:35 UTC,"--shapes=sparse --offsets=tss --n-per=4 --sims=150 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0
+null_unequal_lib_tss_n12_t40,2026-10-22 22:50:11 UTC,"--shapes=unequal_lib --offsets=tss --n-per=12 --sims=150 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0
+null_unequal_lib_tss_n4_t40,2026-10-22 22:51:15 UTC,"--shapes=unequal_lib --offsets=tss --n-per=4 --sims=150 --taxa=40 --tests=wald_ml,lrt --truth=null --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0
+power_original_tss_n12_t40_lfc1,2026-10-22 22:57:57 UTC,"--shapes=original --offsets=tss --n-per=12 --sims=150 --taxa=40 --tests=wald_ml,lrt --truth=power --log2fc=1 --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0
+power_original_tss_n4_t40_lfc1,2026-10-22 22:59:41 UTC,"--shapes=original --offsets=tss --n-per=4 --sims=150 --taxa=40 --tests=wald_ml,lrt --truth=power --log2fc=1 --differential-jl=/work/differential.jl --out=/work/results/raw --results=/work/results/raw",sourced negbin.R (MASS:: stripped from the extracted block; see header),absent,absent,"wald_ml,lrt",R version 4.6.0 (2026-04-24),wasm32-unknown-emscripten (32-bit),webr_0.6.0
diff --git a/bench/null_calibration/summarise.R b/bench/null_calibration/summarise.R
new file mode 100644
index 00000000..21231182
--- /dev/null
+++ b/bench/null_calibration/summarise.R
@@ -0,0 +1,131 @@
+#!/usr/bin/env Rscript
+# SPDX-License-Identifier: AGPL-3.0-only
+# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) <j.d.a.jewell@open.ac.uk>
+#
+# Aggregate calibration.R's per-cell p-value CSVs into the #193 summary
+# tables. Base R only, so it runs anywhere calibration.R's output lands.
+#
+#   Rscript bench/null_calibration/summarise.R \
+#       --results=bench/null_calibration/results/raw \
+#       --out=bench/null_calibration/results
+#
+# Writes calibration_summary.csv (one row per cell x test: rejection rates
+# at alpha = 0.01/0.05/0.10 with Wilson 95% intervals, Kolmogorov-Smirnov D
+# and p against Uniform(0,1), fit-failure share) and, when power cells are
+# present, power_summary.csv (raw and BH power on truly-changed taxa, and
+# the Benjamini-Hochberg false discovery rate at 0.05 averaged per dataset).
+#
+# A "test" here is one taxon in one simulated dataset that produced a
+# p-value (production's BH family: status ok or boundary). Failures are
+# counted separately and never enter a rate's denominator silently: every
+# row reports both n_pvalues and n_failed.
+
+argv <- commandArgs(trailingOnly = TRUE)
+if (length(argv) == 0 && nzchar(Sys.getenv("NB_CALIB_ARGS"))) {
+  argv <- strsplit(trimws(Sys.getenv("NB_CALIB_ARGS")), "\\s+")[[1]]
+}
+args <- list()
+for (a in argv) {
+  if (!grepl("^--[a-z0-9-]+=", a)) stop("unrecognised argument '", a, "' (want --key=value)")
+  kv <- strsplit(sub("^--", "", a), "=", fixed = TRUE)[[1]]
+  args[[kv[1]]] <- paste(kv[-1], collapse = "=")
+}
+results_dir <- if (!is.null(args$results)) args$results else
+  file.path("bench", "null_calibration", "results", "raw")
+out_dir <- if (!is.null(args$out)) args$out else
+  file.path("bench", "null_calibration", "results")
+dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)
+
+files <- list.files(results_dir, pattern = "^pvalues_.*[.]csv$", full.names = TRUE)
+if (length(files) == 0) stop("no pvalues_*.csv under ", results_dir)
+dat <- do.call(rbind, lapply(files, function(f) {
+  d <- read.csv(f, stringsAsFactors = FALSE)
+  d$file <- basename(f)
+  d
+}))
+if (!"taxa" %in% names(dat)) dat$taxa <- NA_integer_
+dat$has_p <- !is.na(dat$pvalue) & dat$status %in% c("ok", "boundary")
+
+# Wilson score interval for k/n at 95% (the same interval prop.test uses
+# without continuity correction, written out so the bench needs nothing
+# beyond base R distributions).
+wilson <- function(k, n) {
+  if (n == 0) return(c(NA_real_, NA_real_, NA_real_))
+  z <- qnorm(0.975)
+  p <- k / n
+  den <- 1 + z^2 / n
+  mid <- (p + z^2 / (2 * n)) / den
+  half <- z * sqrt(p * (1 - p) / n + z^2 / (4 * n^2)) / den
+  c(p, max(0, mid - half), min(1, mid + half))
+}
+
+key_of <- function(d) paste(d$truth, d$shape, d$offset, d$n_per, d$taxa, d$test, sep = "|")
+groups <- split(dat, key_of(dat))
+
+summ_rows <- lapply(groups, function(d) {
+  p <- d$pvalue[d$has_p]
+  # Under a power cell only the truly-null taxa calibrate the test.
+  p_null <- d$pvalue[d$has_p & d$true_log2fc == 0]
+  ks <- if (length(p_null) >= 10) {
+    k <- suppressWarnings(ks.test(p_null, "punif"))
+    c(unname(k$statistic), k$p.value)
+  } else c(NA_real_, NA_real_)
+  w <- function(alpha) wilson(sum(p_null < alpha), length(p_null))
+  w01 <- w(0.01); w05 <- w(0.05); w10 <- w(0.10)
+  data.frame(
+    truth = d$truth[1], shape = d$shape[1], offset = d$offset[1],
+    n_per_group = d$n_per[1], n_taxa_setting = d$taxa[1], test = d$test[1],
+    n_sims = length(unique(d$sim)), n_attempted = nrow(d),
+    n_pvalues = length(p), n_failed = sum(d$status == "failed"),
+    n_boundary = sum(d$status == "boundary"),
+    fail_rate = sum(d$status == "failed") / nrow(d),
+    rate_01 = w01[1], rate_01_lo = w01[2], rate_01_hi = w01[3],
+    rate_05 = w05[1], rate_05_lo = w05[2], rate_05_hi = w05[3],
+    rate_10 = w10[1], rate_10_lo = w10[2], rate_10_hi = w10[3],
+    ks_D = ks[1], ks_p = ks[2], mean_p_null = if (length(p_null)) mean(p_null) else NA_real_,
+    stringsAsFactors = FALSE
+  )
+})
+summary_df <- do.call(rbind, summ_rows)
+summary_df <- summary_df[order(summary_df$truth, summary_df$shape, summary_df$offset,
+                               summary_df$n_per_group, summary_df$test), ]
+write.csv(summary_df, file.path(out_dir, "calibration_summary.csv"), row.names = FALSE)
+cat(sprintf("wrote %s (%d rows)\n", file.path(out_dir, "calibration_summary.csv"), nrow(summary_df)))
+
+pow <- dat[dat$truth == "power", , drop = FALSE]
+if (nrow(pow) > 0) {
+  pow$padj <- NA_real_
+  # BH within each simulated dataset and test, over that dataset's tested
+  # family - the same family production adjusts over.
+  fams <- split(seq_len(nrow(pow)),
+                paste(pow$file, pow$sim, pow$test, sep = "|"))
+  for (idx in fams) {
+    fam <- idx[pow$has_p[idx]]
+    if (length(fam) > 0) pow$padj[fam] <- p.adjust(pow$pvalue[fam], method = "BH")
+  }
+  pg <- split(pow, key_of(pow))
+  prows <- lapply(pg, function(d) {
+    de <- d[d$true_log2fc != 0 & d$has_p, , drop = FALSE]
+    by_sim <- split(d, d$sim)
+    fdp <- vapply(by_sim, function(s) {
+      disc <- s$has_p & !is.na(s$padj) & s$padj < 0.05
+      if (!any(disc)) 0 else mean(s$true_log2fc[disc] == 0)
+    }, numeric(1))
+    data.frame(
+      shape = d$shape[1], offset = d$offset[1], n_per_group = d$n_per[1],
+      test = d$test[1], n_sims = length(by_sim),
+      power_raw_05 = if (nrow(de)) mean(de$pvalue < 0.05) else NA_real_,
+      power_bh_05 = if (nrow(de)) mean(!is.na(de$padj) & de$padj < 0.05) else NA_real_,
+      bh_fdr_05 = mean(fdp),
+      mean_discoveries = mean(vapply(by_sim, function(s)
+        sum(s$has_p & !is.na(s$padj) & s$padj < 0.05), numeric(1))),
+      stringsAsFactors = FALSE
+    )
+  })
+  power_df <- do.call(rbind, prows)
+  power_df <- power_df[order(power_df$shape, power_df$offset,
+                             power_df$n_per_group, power_df$test), ]
+  write.csv(power_df, file.path(out_dir, "power_summary.csv"), row.names = FALSE)
+  cat(sprintf("wrote %s (%d rows)\n", file.path(out_dir, "power_summary.csv"), nrow(power_df)))
+}
+cat("summarise complete\n")
diff --git a/bench/null_calibration/webr_run.mjs b/bench/null_calibration/webr_run.mjs
new file mode 100644
index 00000000..50eedd29
--- /dev/null
+++ b/bench/null_calibration/webr_run.mjs
@@ -0,0 +1,104 @@
+#!/usr/bin/env node
+// SPDX-License-Identifier: AGPL-3.0-only
+// SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) <j.d.a.jewell@open.ac.uk>
+//
+// Sandbox-only driver: run calibration.R / summarise.R under WebR (R in
+// WebAssembly) when no native R is available. The R code executed is the
+// bench's own, unmodified; MASS comes from sourcing MASS's pure-R negbin.R
+// (glm.nb makes no compiled calls), which calibration.R detects and records
+// in its provenance output. Prefer native Rscript with the renv.lock pin
+// whenever it exists - it is faster and carries glmGamPoi/edgeR.
+//
+// One-time setup (outside the repo, so no repo lockfile changes):
+//   mkdir ~/webr-probe && cd ~/webr-probe && npm init -y && npm i webr@0.6.0
+//   git clone https://github.com/cran/MASS.git mass-src
+//   git -C mass-src checkout b6a2274   # MASS 7.3-65, the renv.lock pin
+//
+// Run, from the repo root:
+//   WEBR_PATH=~/webr-probe/node_modules/webr \
+//   node bench/null_calibration/webr_run.mjs \
+//     --mass-r ~/webr-probe/mass-src/R/negbin.R \
+//     --script calibration -- --n-per=4 --sims=200 --tests=wald_ml,lrt
+// Files land in bench/null_calibration/results/raw (override with
+// --host-out). For summarise, existing raw CSVs in --host-out are mounted
+// in first, and the summary CSVs are copied back to --host-summary-out
+// (default bench/null_calibration/results).
+
+import { existsSync, mkdirSync, readFileSync, readdirSync, writeFileSync } from 'node:fs';
+import { dirname, join, resolve } from 'node:path';
+import { fileURLToPath, pathToFileURL } from 'node:url';
+
+const here = dirname(fileURLToPath(import.meta.url));
+const argv = process.argv.slice(2);
+const sep = argv.indexOf('--');
+const driverArgs = sep === -1 ? argv : argv.slice(0, sep);
+const rArgs = sep === -1 ? [] : argv.slice(sep + 1);
+
+function opt(name, fallback) {
+  const hit = driverArgs.find((a) => a.startsWith(`${name}=`));
+  if (hit) return hit.slice(name.length + 1);
+  const i = driverArgs.indexOf(name);
+  return i === -1 ? fallback : driverArgs[i + 1];
+}
+
+const script = opt('--script', 'calibration');
+const massR = opt('--mass-r', '');
+const hostOut = resolve(opt('--host-out', join(here, 'results', 'raw')));
+const hostSummaryOut = resolve(opt('--host-summary-out', join(here, 'results')));
+const webrPath = process.env.WEBR_PATH || 'webr';
+// WEBR_PATH may name the package directory or its ESM entry directly.
+const webrEntry = webrPath === 'webr' ? 'webr'
+  : pathToFileURL(resolve(webrPath.endsWith('.mjs') ? webrPath : join(webrPath, 'dist', 'webr.mjs'))).href;
+const { WebR } = await import(webrEntry);
+
+// No stdout/stderr callbacks here: under Node the WebR worker cannot clone
+// function options, and the defaults already stream R output to this
+// process's stdout/stderr line by line.
+const webR = new WebR({ interactive: false });
+await webR.init();
+
+const put = async (vfsPath, hostPath) => {
+  const dir = vfsPath.slice(0, vfsPath.lastIndexOf('/'));
+  if (dir) await webR.evalRVoid(`dir.create("${dir}", recursive = TRUE, showWarnings = FALSE)`);
+  await webR.FS.writeFile(vfsPath, readFileSync(hostPath));
+};
+
+await put('/work/calibration.R', join(here, 'calibration.R'));
+await put('/work/summarise.R', join(here, 'summarise.R'));
+await put('/work/differential.jl', join(here, '..', '..', 'src', 'analysis', 'differential.jl'));
+if (massR) await put('/work/negbin.R', resolve(massR.replace(/^~/, process.env.HOME || '')));
+
+// Mount any existing raw results so summarise (and resumed calibration
+// cells appending to the same CSV) see them.
+if (existsSync(hostOut)) {
+  for (const f of readdirSync(hostOut)) {
+    if (/\.(csv|txt)$/.test(f)) await put(`/work/results/raw/${f}`, join(hostOut, f));
+  }
+}
+
+const benchArgs = [...rArgs.filter((a) => !a.startsWith('--differential-jl=') && !a.startsWith('--out=') && !a.startsWith('--results=')),
+  '--differential-jl=/work/differential.jl',
+  // calibration.R writes raw cells to --out; summarise.R writes summaries to
+  // --out. They are different directories.
+  script === 'summarise' ? '--out=/work/results' : '--out=/work/results/raw',
+  '--results=/work/results/raw'];
+const quote = (s) => `'${s.replace(/'/g, "'\\''")}'`;
+await webR.evalRVoid(`Sys.setenv(NB_CALIB_ARGS = ${quote(benchArgs.join(' '))})`);
+if (massR) await webR.evalRVoid('source("/work/negbin.R")');
+await webR.evalRVoid('setwd("/work")');
+await webR.evalRVoid(`source("/work/${script}.R")`);
+
+mkdirSync(hostOut, { recursive: true });
+mkdirSync(hostSummaryOut, { recursive: true });
+for (const [vfsDir, hostDir] of [['/work/results/raw', hostOut], ['/work/results', hostSummaryOut]]) {
+  const listing = await (await webR.evalR(
+    `list.files("${vfsDir}", pattern = "\\\\.(csv|txt)$")`)).toJs();
+  for (const name of listing.values) {
+    if (vfsDir.endsWith('/raw') === false && name.startsWith('pvalues_')) continue;
+    const bytes = await webR.FS.readFile(`${vfsDir}/${name}`);
+    writeFileSync(join(hostDir, name), bytes);
+    console.log(`[webr_run] ${name} -> ${join(hostDir, name)}`);
+  }
+}
+await webR.close();
+process.exit(0);
diff --git a/config/defaults/pipeline.yml b/config/defaults/pipeline.yml
index 0be643e4..d9930ae9 100644
--- a/config/defaults/pipeline.yml
+++ b/config/defaults/pipeline.yml
@@ -365,6 +365,12 @@ analysis:
     # steadies low-count taxa -- and fits each taxon through
     # MASS::negative.binomial(theta = 1/dispersion), so the Wald test is
     # unchanged. Requires the glmGamPoi R package; refused for non-nb_glm methods.
+    # Both fixed-theta Wald tests treat the dispersion they use as known rather
+    # than propagating its estimation uncertainty: every ml result therefore
+    # carries a calibration caveat in its diagnostics, strongest when a group
+    # is small, as does a glmGamPoi result whose smaller group has at most 12
+    # samples. The caveat records the limitation; it changes no p-value.
+    # Evidence and the replication protocol: bench/null_calibration.
     dispersion: ml
     # Fraction of samples, 0 to 1, in which a taxon must have reads to be
     # tested. Taxa below it are listed as filtered, without a p-value. Applied
diff --git a/frontend/src/api/types.ts b/frontend/src/api/types.ts
index 8c4edf97..48df92c7 100644
--- a/frontend/src/api/types.ts
+++ b/frontend/src/api/types.ts
@@ -474,6 +474,18 @@ export interface DifferentialResult {
   }
   diagnostics: {
     n_taxa: number; n_tested: number; n_failed: number; n_boundary: number; n_filtered: number
+    /**
+     * nb_glm only: interpretation caveat for the fixed-theta Wald test.
+     * Present for every `ml` result, scaled by the smaller group, and for
+     * `glmGamPoi` results whose smaller group has at most 12 samples. It
+     * changes no statistic.
+     */
+    calibration_caveat?: {
+      severity: 'info' | 'caution' | 'strong'
+      dispersion: 'ml' | 'glmGamPoi'
+      min_group_n: number
+      message: string
+    }
     /** Safety note for two-part models. */
     model_warning?: string
     /** clr_welch only. */
diff --git a/frontend/src/components/DifferentialPanel.tsx b/frontend/src/components/DifferentialPanel.tsx
index 0f61f6e9..545aa875 100644
--- a/frontend/src/components/DifferentialPanel.tsx
+++ b/frontend/src/components/DifferentialPanel.tsx
@@ -237,6 +237,15 @@ export function DifferentialPanel({ study, runs, option, aggregate }: {
               return ` Reference frame taxon: ${s.reference_taxon ?? 'none'} over ${s.n_identifiable_taxa} identifiable taxa${sep}.`
             })()}
           </p>
+          {result.diagnostics.calibration_caveat && (
+            <p role="note" style={{
+              fontSize: '.82rem', margin: '0 0 8px',
+              color: result.diagnostics.calibration_caveat.severity === 'info'
+                ? 'var(--color-muted-fg)' : 'inherit',
+            }}>
+              <strong>Calibration caveat:</strong> {result.diagnostics.calibration_caveat.message}
+            </p>
+          )}
 
           {result.diagnostics.zero_replacement && result.diagnostics.zero_replacement.zeros_replaced > 0 && (() => {
             const z = result.diagnostics.zero_replacement!
diff --git a/frontend/src/components/PipelineStages.tsx b/frontend/src/components/PipelineStages.tsx
index 86ad5d24..9d84a40e 100644
--- a/frontend/src/components/PipelineStages.tsx
+++ b/frontend/src/components/PipelineStages.tsx
@@ -152,7 +152,7 @@ export const CONFIG_DESCRIPTIONS: Record<string, string> = {
   'analysis.beta.depth':          'Reads per sample for rarefy or srs; 0 uses the smallest non-empty sample.',
   'analysis.differential.method': 'Differential model: nb_glm (negative-binomial GLM); clr_welch (Welch test on zero-replaced CLR values); zinb or hurdle_nb (two-part count models with a joint zero/count-process test); ancom_bc (bias-corrected log-linear model); aldex2 (Monte-Carlo CLR); songbird (CLR-centred multinomial regression); multinomial or dirichlet_multinomial (joint compositional models with log ratios against a reference taxon). Compositional effects are relative log-ratios, not absolute abundance changes.',
   '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.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. Both Wald tests treat the dispersion as known, so ml results always carry a calibration caveat (strongest in small groups; see bench/null_calibration).',
   '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. Compositional models filter before constructing their subcomposition.',
   'analysis.differential.replacement_method': 'clr_welch only: bayesian_multiplicative uses cmultRepl GBM posterior estimates (Martín-Fernández et al. 2015); multiplicative_replacement uses CZM (Martín-Fernández et al. 2003). Both preserve observed-part ratios, but every replacement is biased. For sparse data consider NB on raw counts, prevalence filtering, occupancy or zero-inflated models instead.',
   'analysis.differential.replacement_delta': "clr_welch only: cmultRepl frac, strictly between 0 and 1. CZM starts at delta × 0.5/sample total; either method caps high imputations at delta × the taxon's smallest observed proportion. Below 0.01 or at least 0.9 warns. Do not tune delta for significance; every replacement is biased.",
diff --git a/frontend/src/components/differentialTable.test.ts b/frontend/src/components/differentialTable.test.ts
index 523e893b..07e3f1dd 100644
--- a/frontend/src/components/differentialTable.test.ts
+++ b/frontend/src/components/differentialTable.test.ts
@@ -193,3 +193,34 @@ describe('advanced compositional methods (ancom_bc, aldex2, songbird)', () => {
     expect(configLine(sb)).toBe('method songbird, min_prevalence 0.2; estimate is a CLR-centred multinomial differential, not an absolute fold change')
   })
 })
+
+describe('calibration caveat', () => {
+  const caveated: DifferentialResult = {
+    ...result,
+    config: { method: 'nb_glm', offset: 'tss', dispersion: 'ml', min_prevalence: 0 },
+    diagnostics: {
+      ...result.diagnostics,
+      calibration_caveat: {
+        severity: 'strong', dispersion: 'ml', min_group_n: 2,
+        message: 'Wald p-values treat the per-taxon theta as known.',
+      },
+    },
+  }
+
+  test('the settings line records the dispersion the result carries', () => {
+    expect(configLine(caveated)).toBe('method nb_glm, offset tss, dispersion ml, min_prevalence 0')
+  })
+
+  test('the caveat travels into the CSV header, after the settings line', () => {
+    const lines = differentialCsv(caveated).trimEnd().split('\r\n')
+    expect(lines[2]).toBe(`# ${configLine(caveated)}`)
+    expect(lines[3]).toBe('# Calibration caveat (strong): Wald p-values treat the per-taxon theta as known.')
+    expect(lines[4]).toBe('taxon,status,log2_fold_change,estimate,standard_error,statistic,pvalue,padj,dispersion_theta,zero_process_log_odds,prevalence,note')
+    expect(lines).toHaveLength(7)
+  })
+
+  test('a result without a caveat keeps the three-line header', () => {
+    const lines = differentialCsv(result).trimEnd().split('\r\n')
+    expect(lines[3]).toBe('taxon,status,log2_fold_change,estimate,standard_error,statistic,pvalue,padj,dispersion_theta,zero_process_log_odds,prevalence,note')
+  })
+})
diff --git a/frontend/src/components/differentialTable.ts b/frontend/src/components/differentialTable.ts
index f173314f..70dfae03 100644
--- a/frontend/src/components/differentialTable.ts
+++ b/frontend/src/components/differentialTable.ts
@@ -22,10 +22,14 @@ export function csvField(value: string | number | null): string {
  */
 export function differentialCsv(result: DifferentialResult): string {
   const { reference, contrast } = result.groups
+  const caveat = result.diagnostics.calibration_caveat
   const header = [
     `# Differential abundance: ${contrast} vs ${reference} (reference), rank ${result.rank}, table ${result.table}`,
     `# ${result.method}`,
     `# ${configLine(result)}`,
+    // The caveat travels with the numbers it qualifies: a downloaded table
+    // must not look better-calibrated than the result it came from.
+    ...(caveat ? [`# Calibration caveat (${caveat.severity}): ${caveat.message}`] : []),
   ].map(l => l.replace(/[\r\n]+/g, ' '))
   const rows = result.rows.map(r => COLUMNS.map(c => csvField(r[c] ?? null)).join(','))
   return [...header, COLUMNS.join(','), ...rows].join('\r\n') + '\r\n'
diff --git a/src/analysis/differential.jl b/src/analysis/differential.jl
index be010b84..aaec5ba7 100644
--- a/src/analysis/differential.jl
+++ b/src/analysis/differential.jl
@@ -1130,6 +1130,66 @@ _row(taxon, prevalence, note) = Dict{String,Any}(
     "pvalue" => nothing, "padj" => nothing, "dispersion_theta" => nothing,
     "prevalence" => prevalence)
 
+"""
+    _calibration_caveat(dispersion, n_ref, n_con) -> Union{Dict{String,Any},Nothing}
+
+A result-level calibration caveat for the negative-binomial Wald tests, or
+`nothing` when none is recorded. Every `"ml"` result carries one, scaled by
+the smaller group size: `MASS::glm.nb` estimates theta taxon by taxon and the
+fixed-scale Wald test then treats that estimate as known, ignoring its
+sampling uncertainty (issue #193). The `"glmGamPoi"` dispersion trend shares
+information across taxa, so its caveat is recorded only through 12 samples in
+the smaller group; even there, this implementation plugs the trend into the
+same fixed-theta Wald test rather than propagating trend uncertainty or using
+glmGamPoi's quasi-likelihood F-test. The caveat changes no statistic: it is
+provenance and interpretation guidance recorded with the result.
+"""
+function _calibration_caveat(dispersion::AbstractString, n_ref::Integer, n_con::Integer)
+    min_group_n = min(n_ref, n_con)
+    # "caution" runs through 20, the largest group size the bench measured:
+    # the ml Wald test was still anti-conservative there. Beyond the measured
+    # range the caveat drops to "info" but, for "ml", never disappears.
+    severity = min_group_n <= 4 ? "strong" : min_group_n <= 20 ? "caution" : "info"
+    if dispersion == "ml"
+        message = if severity == "strong"
+            "The smallest group has $min_group_n sample(s). With dispersion \"ml\", " *
+            "theta is estimated separately for each taxon and the fixed-scale Wald test " *
+            "then treats that estimate as known. Null simulations for issue #193 rejected " *
+            "substantially too often at 4 samples per group; treat small p-values and " *
+            "Benjamini-Hochberg adjusted p-values at this size as exploratory rather than " *
+            "confirmatory. See bench/null_calibration."
+        elseif severity == "caution"
+            "The smallest group has $min_group_n samples. With dispersion \"ml\", theta " *
+            "is estimated separately for each taxon and the fixed-scale Wald test then " *
+            "treats that estimate as known. Null simulations in bench/null_calibration " *
+            "(issue #193) found above-nominal rejection at every measured group size " *
+            "through 20 per group, so interpret borderline p-values and " *
+            "Benjamini-Hochberg adjusted p-values cautiously."
+        else
+            "With dispersion \"ml\", theta is estimated separately for each taxon and " *
+            "the fixed-scale Wald test treats that estimate as known; theta uncertainty is " *
+            "not propagated. Null simulations in bench/null_calibration (issue #193) " *
+            "found above-nominal rejection through 20 samples per group, the largest " *
+            "size measured; beyond that range the approximation is expected to improve, " *
+            "but it has not been measured there."
+        end
+        return Dict{String,Any}("severity" => severity, "dispersion" => "ml",
+                                 "min_group_n" => min_group_n, "message" => message)
+    elseif dispersion == "glmGamPoi" && min_group_n <= 12
+        message = "The smallest group has $min_group_n sample(s). The glmGamPoi " *
+                  "dispersion trend shares information across taxa and steadies theta " *
+                  "relative to a separate maximum-likelihood estimate for every taxon, " *
+                  "but this implementation still plugs the trend into a fixed-theta, " *
+                  "fixed-scale Wald test. It does not propagate trend uncertainty and " *
+                  "does not implement glmGamPoi's quasi-likelihood F-test, so interpret " *
+                  "borderline p-values cautiously at this sample size. " *
+                  "See bench/null_calibration and issue #193."
+        return Dict{String,Any}("severity" => severity, "dispersion" => "glmGamPoi",
+                                 "min_group_n" => min_group_n, "message" => message)
+    end
+    nothing
+end
+
 # Two-group ANCOM-BC (Lin & Peddada 2020, Nat. Commun. 11:3514):
 # (1) Estimate within-group sample-specific sampling-fraction offsets on
 #     log(counts + 1) by inverse-variance-weighted least squares, centred within
@@ -1993,7 +2053,12 @@ enter each negative-binomial model as `log(factor)`; the estimate is a
 natural-log fold change. `config.dispersion` selects per-taxon maximum-likelihood
 theta (`"ml"`) or fixed theta from `glmGamPoi::glm_gp`'s cross-taxon dispersion trend
 (`"glmGamPoi"`), fitted in one call on the taxa passing the prevalence filter.
-The result's `config` includes the selected dispersion name.
+The result's `config` includes the selected dispersion name. Its `diagnostics`
+also records a `calibration_caveat` for every `"ml"` result (scaled by the
+smaller group) and for a `"glmGamPoi"` result whose smaller group has at most
+12 samples: both fixed-theta Wald paths treat the dispersion used by the test
+as known rather than propagating its estimation uncertainty. The caveat alters
+no estimate, p-value or adjusted p-value.
 
 `clr_welch`, in this order: (1) `min_prevalence` chooses the taxa, and taxa with
 no reads in any sample are always left out; (2) the zeros of those taxa only
@@ -2090,6 +2155,8 @@ 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)
+        caveat = _calibration_caveat(config.dispersion, n_ref, n_con)
+        isnothing(caveat) || (diagnostics["calibration_caveat"] = caveat)
     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; " *
diff --git a/src/server/server.jl b/src/server/server.jl
index 44f916f4..ebf5abb6 100644
--- a/src/server/server.jl
+++ b/src/server/server.jl
@@ -116,6 +116,7 @@ module Server
     include(joinpath(@__DIR__, "routes", "doi.jl"))
     include(joinpath(@__DIR__, "routes", "figures.jl"))
     include(joinpath(@__DIR__, "routes", "phylogeny.jl"))
+    include(joinpath(@__DIR__, "routes", "cladistics.jl"))
     include(joinpath(@__DIR__, "routes", "clade_cumulus.jl"))
     include(joinpath(@__DIR__, "routes", "cladistics.jl"))
 
diff --git a/test/unit/test_differential.jl b/test/unit/test_differential.jl
index 234437dc..b8543eff 100644
--- a/test/unit/test_differential.jl
+++ b/test/unit/test_differential.jl
@@ -220,6 +220,23 @@ end
         @test_throws ArgumentError replace_zeros([0 1; 2 3]; method="pseudocount")
     end
 
+    @testset "Calibration caveat scales with the smaller group" begin
+        # Pure thresholds, no R needed: every "ml" result is caveated, the
+        # glmGamPoi trend only through 12 in the smaller group.
+        caveat = Differential._calibration_caveat
+        @test caveat("ml", 4, 4)["severity"] == "strong"
+        @test caveat("ml", 30, 3)["severity"] == "strong"
+        @test caveat("ml", 5, 5)["severity"] == "caution"
+        @test caveat("ml", 12, 12)["severity"] == "caution"
+        @test caveat("ml", 13, 20)["severity"] == "caution"
+        @test caveat("ml", 13, 20)["min_group_n"] == 13
+        @test caveat("ml", 21, 25)["severity"] == "info"
+        @test caveat("glmGamPoi", 4, 9)["severity"] == "strong"
+        @test caveat("glmGamPoi", 12, 12)["severity"] == "caution"
+        @test isnothing(caveat("glmGamPoi", 13, 13))
+        @test isnothing(caveat("glmGamPoi", 20, 13))
+    end
+
     @testset "Input validation" begin
         taxa = ["a", "b", "c"]
         groups = ["A", "A", "B"]
@@ -258,6 +275,9 @@ end
                                    "TSS size-factor offset; Wald test of the group coefficient; " *
                                    "Benjamini-Hochberg adjustment over the fitted taxa"
             @test res["diagnostics"]["n_tested"] == 3
+            caveat = res["diagnostics"]["calibration_caveat"]
+            @test caveat["severity"] == "caution" && caveat["dispersion"] == "ml"
+            @test caveat["min_group_n"] == 12 && occursin("theta", caveat["message"])
             @test all(r -> r["status"] == "ok", res["rows"])
             @test isapprox(byt["t_up"]["estimate"], log(4); atol=0.5)
             @test isapprox(byt["t_up"]["log2_fold_change"], 2; atol=0.8)
@@ -271,6 +291,13 @@ end
                                              reference="B", contrast="A")
             sb = Dict(r["taxon"] => r for r in swapped["rows"])
             @test sb["t_up"]["estimate"] ≈ -byt["t_up"]["estimate"] rtol=1e-6
+            # The caveat strengthens as the smaller group shrinks; the p-values
+            # themselves are untouched by it (checked against MASS below).
+            keep = [1:4; 13:16]
+            small = differential_abundance(effect[keep, :], effect_samples[keep], effect_taxa,
+                                           effect_groups[keep]; reference="A", contrast="B")
+            @test small["diagnostics"]["calibration_caveat"]["severity"] == "strong"
+            @test small["diagnostics"]["calibration_caveat"]["min_group_n"] == 4
 
             @testset "matches a direct MASS::glm.nb call" begin
                 off = log.(tss_factors(effect, effect_samples))
@@ -308,6 +335,10 @@ end
                                    "with dispersion trend from glmGamPoi::glm_gp, log link) with a " *
                                    "TSS size-factor offset; Wald test of the group coefficient; " *
                                    "Benjamini-Hochberg adjustment over the fitted taxa"
+            caveat = res["diagnostics"]["calibration_caveat"]
+            @test caveat["severity"] == "caution" && caveat["dispersion"] == "glmGamPoi"
+            @test caveat["min_group_n"] == 12 &&
+                  occursin("quasi-likelihood F-test", caveat["message"])
             byt = Dict(r["taxon"] => r for r in res["rows"])
             off = log.(tss_factors(effect, effect_samples))
             g = effect_groups
@@ -876,6 +907,7 @@ end
             @test isempty(res["size_factors"])
             @test res["config"]["method"] == "clr_welch" && !haskey(res["config"], "offset")
             @test res["diagnostics"]["zero_replacement"]["zeros_replaced"] == 0
+            @test !haskey(res["diagnostics"], "calibration_caveat")
             @test occursin("Welch", res["method"]) && occursin("cmultRepl", res["method"])
             # Known answer, independent of R: the estimate is mean(B) - mean(A)
             # of each taxon's CLR.