bench(null_calibration): replicate #193 at 1,000 sims and measure the alternatives

bench/null_calibration/README.md

@@ -21,6 +21,7 @@ test gets an easier dataset):
 | `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 |
+| `ql_f` | quasi-likelihood F test on the same full `glm.nb` fit as `lrt`: `φ = Pearson X² / residual df`, `F = (β/se)² / φ` against `F(1, residual df)`. A tagwise dispersion scale, **not** empirical-Bayes shrinkage and **not** glmGamPoi's `test_de` | 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 |
 
@@ -52,9 +53,13 @@ file), never silently replaced.
   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).
+- **Seeds**: a per-cell key hashed from the cell's identity string (truth,
+  shape, offset, samples per group, taxa, log2FC), then
+  `seed0 + key × 1000003 + sim`. One seed per dataset, so a cell reproduces
+  exactly whatever else is requested, in whatever order or process. (It
+  previously used the cell's *index* in the invocation, which changed every
+  dataset when a grid was split across processes.) The key is recorded in each
+  provenance file. Default `seed0 = 20261020` (the #193 seed).
 
 ### Reading a cell
 
@@ -76,21 +81,31 @@ 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
+  --tests=wald_ml,lrt,ql_f
 
 # the full #193 grid — 1,000 sims per cell; expect hours, use several
-# processes by splitting --n-per / --shapes across terminals
+# processes by splitting --n-per / --shapes across terminals (with the
+# identity-hashed seeds, splitting the grid does not change any dataset)
 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).
+taxon × test. Re-running a cell **resumes** it: simulations whose number is
+already in the file are skipped, so an interrupted long cell costs only the
+current simulation. Because the per-sim seed is a function of cell identity and
+sim number, the resumed rows are the same rows the uninterrupted run would have
+written. A cell whose provenance records a different seed key is refused rather
+than appended to; delete that cell's two files to start it over.
+
+Each cell also writes `provenance_<cell>.txt` with the arguments, seed and cell
+key, executor (webR/native R), package versions, MD5s of `calibration.R` and of
+the extracted `_FIT_R` block, and `sessionInfo()`. `summarise.R` writes
+`results/calibration_summary.csv`, `power_summary.csv` when power cells exist,
+and `provenance_summary.csv` — the run record, generated from the per-cell
+provenance files rather than assembled by hand.
 
 ### No native R? The WebR fallback
 

bench/null_calibration/UPSTREAM_DRAFT.md

@@ -1,51 +1,88 @@
 # 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.*
+upstream. Every number is from `REPORT.md` in this directory, which defines
+exactly what ran, where, and what did not run. This replaces the earlier draft
+based on the reduced 2026-10-22 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).
+(`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, including 20 samples per group. At α = 0.05 the rejection rate is 0.172
+at 3/group, 0.140 at 4, 0.104 at 6, 0.090 at 8, 0.075 at 12 and 0.064 at 20
+(1,000 sims × 40 taxa per cell; every Wilson 95% lower bound is above 0.05).
 
 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.
+theta and ignores its sampling uncertainty.
 
-## What was and was not measured
+Two cheaper alternatives, measured on identical datasets:
 
-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.
+- a likelihood-ratio test on the same two `glm.nb` fits removes about a third of
+  the excess at 4/group (0.140 → 0.101) but is still anti-conservative in every
+  cell;
+- a quasi-likelihood F test (`φ` = Pearson X² / residual df, `F = (β/se)²/φ`
+  against `F(1, df)`; a tagwise scale, **not** glmGamPoi's or edgeR's
+  empirical-Bayes version) is much closer: 0.044 at 3/group (conservative),
+  0.055 at 4, 0.058–0.060 through 8, 0.056 at 20. It remains ~1 point above
+  nominal, not at nominal.
 
-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.
+The inflation is not specific to one data shape: constant dispersion, sparse
+zero-heavy tables, unequal library sizes and RLE offsets all reproduce it
+(4/group Wald 0.112–0.142; 12/group 0.070–0.081). RLE is not a remedy; at
+12/group it was the worst configuration measured.
+
+## What this means for Benjamini–Hochberg
+
+Under a power design (25% of taxa at |log2FC| = 1, 400 sims), BH at 0.05:
+
+| n/group | test | BH power | BH FDR |
+|---:|---|---:|---:|
+| 4 | Wald (production) | 0.275 | **0.362** |
+| 4 | LRT | 0.115 | 0.192 |
+| 4 | QL F (`ql_f`) | 0.022 | 0.037 |
+| 12 | Wald (production) | 0.536 | **0.118** |
+| 12 | LRT | 0.468 | 0.078 |
+| 12 | QL F (`ql_f`) | 0.433 | 0.063 |
+
+The default's apparent power advantage at small n is an artefact of an
+uncontrolled FDR. The QL F test is the only measured option whose FDR is near
+nominal, but at 4/group it is nearly powerless after BH because its F reference
+has only six residual degrees of freedom (`F(1,6)` ≈ 5.99, |t| > 2.45).
 
 ## Options for the upstream maintainer's consideration
 
-1. **Result-notes caveat for small samples** (adopted downstream): record a
+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.
+   Changes no statistic or default. The downstream text is justified unchanged
+   by this run: Wald is above nominal at every measured size through 20.
+2. **Likelihood-ratio test for the `glm_nb` path.** Measurably closer to nominal
+   but still anti-conservative everywhere, slower, and it fails more often at
+   tiny n. It is a mitigation, not a calibration fix; the numbers above quantify
+   the trade-off.
+3. **Quasi-likelihood F test below a group-size threshold.** Closest to nominal
+   at every size and the only option with near-nominal BH FDR, but it changes
+   critical values substantially at 3–4 samples per group and costs almost all
+   BH power there. glmGamPoi's own QL F-test (the calibrated path in the
+   downstream issue) is the version worth considering rather than a tagwise
+   reimplementation, and it was not measured in this sandbox (see `REPORT.md`
+   §"What was not run" for the native-R command).
+
+The downstream recommendation is **1 now, 2 and 3 as an upstream decision**.
+A default change would need the glmGamPoi QL F row measured on native R first;
+this bench provides the code path and the protocol for it.
+
+## Reproducing
+
+```sh
+Rscript -e 'renv::restore()'
+# the 16 cells actually run are listed in REPORT.md §"Cells run"
+Rscript bench/null_calibration/calibration.R --shapes=original --offsets=tss \
+  --n-per=3,4,6,8,12,20 --sims=1000 --taxa=40 --tests=wald_ml,lrt,ql_f
+Rscript bench/null_calibration/summarise.R
+```
+
+The committed `results/` CSVs carry every rate and interval quoted above;
+`results/provenance_summary.csv` carries the executor, the package versions and
+the input digests for each cell.