nb_glm's default glm_nb Wald test is anti-conservative at small n: replicate before raising upstream #193

What was measured

A null calibration found that the default nb_glm path's p-values reject too often in small studies. The default path is analysis.differential.dispersion: glm_nb: MASS::glm.nb per taxon, a Wald test at fixed scale 1, from upstream PR 24. The glmGamPoi quasi-likelihood F-test, added as an option in upstream PR 62, stayed at nominal.

The setup: the real fit_nb R block from src/analysis/differential.jl, branch claude/elegant-rubin-z92rx4. There were no true effects, TSS offsets, and library sizes drawn from lognormal(0, 0.3). Taxon means ran from 2 to 500, and NB dispersion was 0.05 + 2 / sqrt(mu), falling with the mean.

The binomial SE at 0.05 is about 0.003. At 4 vs 4, 154 glm_nb fits also failed, against 9 for glmGamPoi. Environment: R 4.3.3, MASS 7.3-60.0.1, glmGamPoi 1.14.3 with 1.22.0's pure-R test_de/predict loaded over it.

Why this is expected: a Wald test that plugs in the per-taxon ML theta ignores the uncertainty in theta. That uncertainty is large when n is small. The QL F-test carries it in its denominator degrees of freedom (Lund et al., 2012).

What this issue is, and is not

This is a fork-side tracking issue, kept here to replicate the finding before anything is raised with JoshuaJewell/MetaManifold-WebUI. The default is his decision. The glm_nb default is unchanged by PR 62, and changing it is out of scope for issue 37.

Proposed replication, before raising upstream

• More runs: n per group in {3, 4, 6, 8, 12, 20}, with 1,000 simulations per cell. Report Wilson intervals at alpha = 0.01, 0.05 and 0.10, plus a KS test of the null p-values against Uniform(0, 1).
• More data shapes: constant dispersion; sparse, amplicon-like tables (many zeros, a few dominant taxa); library sizes that differ between the groups; RLE as well as TSS offsets.
• More tests: glm_nb Wald (current), the glm_nb likelihood-ratio test against the intercept-only fit, glmGamPoi QL F, and edgeR's QL F as an outside reference only, not a dependency.
• Power: at fixed fold changes, for whichever tests stay calibrated, and the BH false discovery rate at 0.05.
• Where it lives: as a bench (e.g. bench/null_calibration/, next to bench/layer1_mock_recovery/), not a unit test. It is slow and stochastic.

Options to put to the maintainer once it replicates, rated:

1. Use a likelihood-ratio test for glm_nb. Most correct while keeping the model. It still ignores theta uncertainty, but is usually better calibrated than Wald. The bench above should measure it.
2. Warn about small n in the result's notes. Most dependable, and cheap. It changes no numbers.
3. Make glmGamPoi the default below a sample threshold. A different path, but a default change and a new dependency on every path. His call.

Reproduction script

# Null calibration of the two nb_glm dispersion settings, through the _FIT_R block.
# Usage: Rscript calibration.R <fit_block.R> <n_per_group> <n_sims> <n_taxa>
# <fit_block.R> is the text of `const _FIT_R = raw"""..."""` from src/analysis/differential.jl.
# With glmGamPoi >= 1.22 (the renv.lock pin) no overlay is needed; drop the source() line.
args <- commandArgs(trailingOnly = TRUE)
fit_code <- paste(readLines(args[[1]]), collapse = "\n")
n_per <- as.integer(args[[2]]); n_sims <- as.integer(args[[3]]); n_taxa <- as.integer(args[[4]])
source(file.path(dirname(normalizePath(args[[1]])), "overlay_profile.R"))
DA_VARS <- c("da_counts", "da_group", "da_levels", "da_offset", "da_theta_upper", "da_theta_lower",
             "da_use_glp", "da_g", "da_term", "da_n", "da_result", "da_glp", "da_glp_failed", "da_gp",
             "da_j", "da_y", "da_absent", "da_od", "da_est", "da_se", "da_stat", "da_p",
             "da_warns", "da_fit", "da_co", "da_row", "da_theta", "da_design", "da_contrast")
run_block <- function(counts, use_glp) {
  g <- globalenv(); lib <- rowSums(counts); lv <- c("A", "B")
  assign("da_counts", counts, g); assign("da_group", rep(lv, each = n_per), g)
  assign("da_levels", lv, g); assign("da_offset", log(lib / exp(mean(log(lib)))), g)
  assign("da_theta_upper", 1e7, g); assign("da_theta_lower", 1e-8, g); assign("da_use_glp", use_glp, g)
  on.exit(rm(list = intersect(ls(g), DA_VARS), envir = g))
  suppressWarnings(eval(parse(text = fit_code), envir = g))
  get("da_result", envir = g)
}
set.seed(20261020)
p_nb <- p_gp <- c(); fails <- c(nb = 0, gp = 0)
for (s in seq_len(n_sims)) {
  mu <- exp(runif(n_taxa, log(2), log(500)))
  disp <- 0.05 + 2 / sqrt(mu)
  lib <- exp(rnorm(2 * n_per, 0, 0.3))
  y <- sapply(seq_len(n_taxa), function(j) rnbinom(2 * n_per, mu = mu[j] * lib, size = 1 / disp[j]))
  storage.mode(y) <- "integer"
  y <- y[, colSums(y > 0) >= 2 & apply(y, 2, function(v) length(unique(v)) > 1), drop = FALSE]
  if (ncol(y) < 5) next
  a <- run_block(y, FALSE); b <- run_block(y, TRUE)
  p_nb <- c(p_nb, a$pvalue[a$status %in% c("ok", "boundary")]); fails["nb"] <- fails["nb"] + sum(a$status == "failed")
  p_gp <- c(p_gp, b$pvalue[b$status %in% c("ok", "boundary")]); fails["gp"] <- fails["gp"] + sum(b$status == "failed")
}
rate <- function(p, a) sprintf("%.4f", mean(p < a))
cat(sprintf("glm_nb    %d tests  %s %s %s  failed %d\n", length(p_nb), rate(p_nb, .01), rate(p_nb, .05), rate(p_nb, .10), fails[["nb"]]))
cat(sprintf("glmGamPoi %d tests  %s %s %s  failed %d\n", length(p_gp), rate(p_gp, .01), rate(p_gp, .05), rate(p_gp, .10), fails[["gp"]]))

Related: #192, which prepares the upstream PR 62 fixes. This issue is filed under #1, the validated statistics layer.