Null-calibration bench + small-sample caveat (#193), with cross-PR integration repair

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.

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/

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.

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.

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.