Prepare upstream PR 62 review fixes: theta, label-safe glmGamPoi contrast

README.md

@@ -249,7 +249,7 @@ is the reference. Counts are aggregated by taxon after applicable
 read-count bounds to the differential-abundance request, so those bounds do
 not exclude samples from this analysis.
 
-All `analysis.differential` options below (`method`, `offset`,
+All `analysis.differential` options below (`method`, `offset`, `dispersion`,
 `min_prevalence` and `replacement_delta`) are resolved at study scope.
 Set them in the study configuration; group and run overrides are not used
 by this endpoint, including when called from a group Analysis workspace.
@@ -259,7 +259,15 @@ by this endpoint, including when called from a group Analysis workspace.
 - `nb_glm` (default): a negative-binomial GLM per taxon using `MASS::glm.nb`,
   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` scaling.
+  (default) or `rle` scaling. `dispersion` selects the fits and the test:
+  `glm_nb` (default) fits `MASS::glm.nb` per taxon with its own
+  maximum-likelihood theta and Wald-tests the group coefficient; `glmGamPoi`
+  fits one `glmGamPoi::glm_gp` across the taxa passing `min_prevalence`,
+  sharing a dispersion trend across them, and reports `glmGamPoi::test_de`'s
+  quasi-likelihood F-test. Under `glmGamPoi` the p-value comes from that
+  F-test, the standard error is `test_de`'s `lfc_se`, and theta is each
+  taxon's own estimate before shrinkage. It requires the pinned glmGamPoi R
+  package and is refused for `clr_welch`.
 - `clr_welch`: Welch's t-test on centred log-ratios after multiplicative zero
   replacement with `zCompositions::cmultRepl`. The effect is a CLR difference,
   not a fold change. `replacement_delta` controls the replacement fraction

config/defaults/pipeline.yml

@@ -363,9 +363,10 @@ analysis:
     # the group coefficient. glmGamPoi fits one glmGamPoi::glm_gp across the
     # taxa passing min_prevalence, with overdispersion shrinkage toward the
     # cross-taxon trend, and reports glmGamPoi::test_de's quasi-likelihood
-    # F-test of the group contrast with its log2 fold changes. A zero shrunken
-    # dispersion withholds inference for its taxon. Requires the pinned
-    # glmGamPoi R package. Refused with clr_welch.
+    # F-test of the group contrast with its log2 fold changes; the reported
+    # theta inverts glm_gp's negative-binomial overdispersion, and a taxon
+    # whose overdispersion is zero is marked at the Poisson boundary, as under
+    # glm_nb. Requires the pinned glmGamPoi R package. Refused with clr_welch.
     dispersion: glm_nb
     # 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. For

frontend/src/api/types.ts

@@ -439,7 +439,7 @@ export interface DifferentialResult {
     min_prevalence: number
     /** nb_glm only. */
     offset?: 'tss' | 'rle'
-    /** nb_glm only: where each taxon's theta comes from. */
+    /** nb_glm only: which fits and test produced the rows (per-taxon MASS::glm.nb with a Wald test, or one glmGamPoi::glm_gp with test_de's quasi-likelihood F-test). */
     dispersion?: 'glm_nb' | 'glmGamPoi'
     /** clr_welch only. */
     replacement_delta?: number

src/analysis/differential.jl

@@ -291,27 +291,32 @@ da_result <- data.frame(status = rep("ok", da_n), note = rep("", da_n),
                         estimate = rep(NA_real_, da_n), se = rep(NA_real_, da_n),
                         statistic = rep(NA_real_, da_n), pvalue = rep(NA_real_, da_n),
                         theta = rep(NA_real_, da_n), stringsAsFactors = FALSE)
-# "glmGamPoi": one glmGamPoi::glm_gp fit across all taxa supplies each taxon's
-# "glmGamPoi": one glmGamPoi::glm_gp fit across all taxa, with each taxon's
-# overdispersion shrunk toward the cross-taxon trend, supplies glmGamPoi's own
-# published inference -- glmGamPoi::test_de's quasi-likelihood F-test of the
-# group coefficient. Its genes x samples convention, so
-# the counts are transposed, and size_factors takes the linear size factor
-# (the per-taxon fits below use its log as their offset). A fit- or test-level
-# error fails every taxon. No fallback to another dispersion estimator or
-# significance test is permitted.
+# "glmGamPoi": one glmGamPoi::glm_gp fit across all taxa supplies glmGamPoi's
+# own published inference -- glmGamPoi::test_de's quasi-likelihood F-test of
+# the group coefficient. glm_gp refits the coefficients at its cross-taxon
+# dispersion trend and shrinks each taxon's quasi-likelihood dispersion by
+# empirical Bayes. glm_gp works genes x samples, so the counts are transposed,
+# and size_factors takes the linear size factor (the per-taxon fits below use
+# its log as their offset). The contrast reaches test_de as numeric weights on
+# the design's columns, never as text: test_de parses and evaluates a text
+# contrast as R code, so a condition label such as "run-B" or "grpB/run_2"
+# would stop it, and a label that reads as R code would run. A fit- or
+# test-level error fails every taxon. No fallback to another dispersion
+# estimator or significance test is permitted.
 da_glp <- NULL
 da_glp_failed <- ""
 if (da_use_glp) {
+  da_design <- model.matrix(~ da_g)
+  da_contrast <- as.numeric(colnames(da_design) == da_term)
   da_gp <- tryCatch(
-    glmGamPoi::glm_gp(t(da_counts), design = model.matrix(~ da_g),
+    glmGamPoi::glm_gp(t(da_counts), design = da_design,
                       size_factors = exp(da_offset), overdispersion_shrinkage = TRUE),
     error = function(e) e)
   if (inherits(da_gp, "error")) {
     da_glp_failed <- paste("glm_gp stopped with an error:", conditionMessage(da_gp))
   } else {
     da_glp <- tryCatch(
-      glmGamPoi::test_de(da_gp, contrast = da_term, max_lfc = Inf, compute_lfc_se = TRUE),
+      glmGamPoi::test_de(da_gp, contrast = da_contrast, max_lfc = Inf, compute_lfc_se = TRUE),
       error = function(e) e)
     if (inherits(da_glp, "error")) {
       da_glp_failed <- paste("test_de stopped with an error:", conditionMessage(da_glp))
@@ -349,20 +354,20 @@ for (da_j in seq_len(da_n)) {
     # glmGamPoi's own test on its own fit. `lfc` is the log2 fold change of
     # the group contrast and `lfc_se` its standard error; both are converted
     # to the natural-log scale so every method reports one effect scale. The
-    # per-taxon scale the test used is the shrunken quasi-likelihood
-    # dispersion `ql_disp_shrunken`, whose inverse is reported as the
-    # dispersion theta. Zero dispersion is the Poisson boundary, where the F
-    # ratio would report p = 0 for the wrong reason, so inference is withheld;
-    # an unusable dispersion fails only its own taxon.
-    da_ql <- da_gp[["overdispersion_shrinkage_list"]][["ql_disp_shrunken"]][da_j]
-    if (!is.finite(da_ql) || da_ql < 0) {
+    # reported theta is the inverse of glm_gp's negative-binomial
+    # overdispersion for the taxon (`overdispersions`: variance = mu + mu^2 *
+    # overdispersion), the same parameter glm.nb's theta inverts. It is the
+    # taxon's own estimate before shrinkage: the coefficients and the F-test's
+    # reduced fit use the cross-taxon trend, and the F ratio divides by the
+    # shrunken quasi-likelihood dispersion; only lfc_se, from predict(), is
+    # weighted by this estimate. Zero overdispersion is the Poisson boundary:
+    # the row is marked, as under glm_nb, and keeps its statistics, since the
+    # test does not read it. An unusable overdispersion fails only its own
+    # taxon.
+    da_od <- da_gp[["overdispersions"]][da_j]
+    if (!is.finite(da_od) || da_od < 0) {
       da_result[da_j, "status"] <- "failed"
-      da_result[da_j, "note"] <- "the shrunken dispersion for this taxon is not usable"
-      next
-    }
-    if (da_ql == 0) {
-      da_result[da_j, "status"] <- "boundary"
-      da_result[da_j, "note"] <- "zero dispersion: Poisson boundary; quasi-likelihood inference withheld"
+      da_result[da_j, "note"] <- "the negative-binomial overdispersion for this taxon is not usable"
       next
     }
     da_est <- da_glp[["lfc"]][da_j] * log(2)
@@ -376,7 +381,11 @@ for (da_j in seq_len(da_n)) {
       next
     }
     da_result[da_j, c("estimate", "se", "statistic", "pvalue", "theta")] <-
-      c(da_est, da_se, da_stat, da_p, 1 / da_ql)
+      c(da_est, da_se, da_stat, da_p, 1 / da_od)
+    if (da_od == 0) {
+      da_result[da_j, "status"] <- "boundary"
+      da_result[da_j, "note"] <- "the taxon's negative-binomial overdispersion estimate is zero: these counts show no overdispersion beyond Poisson, so theta is unbounded; the test used the cross-taxon trend, so its statistics stand"
+    }
   } else {
     # "glm_nb" estimates theta per taxon inside MASS::glm.nb's alternating
     # fit, then Wald-tests the group coefficient at fixed scale 1, which is
@@ -453,17 +462,21 @@ toward the cross-taxon trend and reports `glmGamPoi::test_de`'s
 quasi-likelihood F-test of the group contrast — its `lfc` (log2 fold change)
 and `lfc_se`, both converted to the natural-log scale so the columns carry one
 effect scale. In that path `statistic` is the quasi-likelihood F, `pvalue`
-its p-value, and `dispersion_theta` the inverse of the shrunken
-quasi-likelihood dispersion the test used for the taxon. Names are stripped
-and matched without regard to case.
+its p-value, and `theta` the inverse of `glm_gp`'s negative-binomial
+overdispersion for the taxon (its `overdispersions` element, the taxon's own
+estimate before shrinkage), the same parameter `MASS::glm.nb`'s theta inverts;
+it is `nothing` where that overdispersion is zero. The contrast is passed to
+`test_de` as numeric weights, so condition labels are never parsed as R code.
+Names are stripped and matched without regard to case.
 
 Returns one row per taxon with `status` (`ok`, `boundary` or `failed`), `note`,
 `estimate`, `se`, `statistic`, `pvalue` and `theta`; a failed taxon has
 `nothing` for every statistic. Boundary fits retain their statistics: under
 `"glm_nb"` theta at or above 1e7 marks the Poisson end and theta at or below
-1e-8 marks extreme variance; under `"glmGamPoi"` a zero shrunken dispersion
-marks the Poisson boundary, where the F ratio would report p = 0 for the wrong
-reason, so inference is withheld. Constant counts, complete separation and
+1e-8 marks extreme variance; under `"glmGamPoi"` a zero negative-binomial
+overdispersion marks the Poisson boundary, and the row keeps its statistics,
+since the F-test divides by the shrunken quasi-likelihood dispersion, which
+stays defined there. Constant counts, complete separation and
 unusable coefficients produce failed rows. Errors caught from the model-fitting
 calls also produce failed rows; a caught `glm_gp` or `test_de` error fails every
 taxon, with no fallback to `"glm_nb"`.
@@ -500,8 +513,9 @@ function fit_nb(counts::AbstractMatrix{<:Integer}, groups::AbstractVector{<:Abst
         finally
             RCall.reval("rm(list = intersect(ls(), 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_ql', 'da_est', 'da_se', 'da_stat', " *
+                        "'da_term', 'da_n', 'da_result', 'da_design', 'da_contrast', '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'))); invisible(gc())")
         end
     end

test/unit/test_differential.jl

@@ -295,7 +295,9 @@ end
             # quasi-likelihood F-test, computed independently. lfc is log2 and
             # lfc_se its standard error, so both are scaled by log(2) to meet
             # the natural-log convention of every method's estimate and se.
-            # One matrix column per quantity, one row per taxon.
+            # theta inverts glm_gp's negative-binomial overdispersion, not the
+            # quasi-likelihood scale the F-test divides by. One matrix column
+            # per quantity, one row per taxon.
             direct = rcopy(R"""
                 local({
                   g <- factor($g, levels = c("A", "B"))
@@ -304,7 +306,7 @@ end
                   td <- glmGamPoi::test_de(gp, contrast = "gB", max_lfc = Inf,
                                            compute_lfc_se = TRUE)
                   cbind(td$pval, td$f_statistic, td$lfc * log(2), td$lfc_se * log(2),
-                        1 / gp$overdispersion_shrinkage_list$ql_disp_shrunken)
+                        1 / gp$overdispersions)
                 })
             """)
             for (j, t) in enumerate(effect_taxa)
@@ -320,43 +322,86 @@ end
         end
     end
 
-    @testset "glmGamPoi unusable dispersions are isolated per taxon" begin
+    @testset "glmGamPoi unusable overdispersions are isolated per taxon" begin
         if !HAVE_GLMGAMPOI
             # Issue #37's contract: these comparisons fail, not skip, in CI.
             @test !_GLMGAMPOI_IN_CI
             @test_skip false
         else
+            # test_de refits the reduced model from the shrinkage list's
+            # dispersion trend and divides by ql_disp_shrunken, so altering
+            # `overdispersions` after the fit changes theta, the boundary check
+            # and lfc_se (predict() weights by it), never the F statistic or
+            # its p-value.
             MetaManifold.RRuntime.with_r_lock() do
                 reval(raw"""
                     test_original_gp <- get("glm_gp", asNamespace("glmGamPoi"))
-                    test_ql <- 0
+                    test_od <- 0
                     assignInNamespace("glm_gp", function(...) {
                         fit <- test_original_gp(...)
-                        fit$overdispersion_shrinkage_list$ql_disp_shrunken[1] <- test_ql
+                        fit$overdispersions[1] <- test_od
                         fit
                     }, ns="glmGamPoi")
                 """)
                 try
                     off = log.(tss_factors(effect, effect_samples))
-                    for (ql, status) in ((0.0, "boundary"), (-1.0, "failed"), (NaN, "failed"))
-                        RCall.globalEnv[:test_ql] = ql
+                    # A zero overdispersion is the Poisson boundary: marked, with
+                    # the statistics kept and no finite theta, as under glm_nb.
+                    # A negative or non-finite one fails its taxon alone.
+                    for (od, status, tested) in ((0.0, "boundary", true), (-1.0, "failed", false),
+                                                 (NaN, "failed", false))
+                        RCall.globalEnv[:test_od] = od
                         rows = fit_nb(effect, effect_groups, off; levels=("A", "B"),
                                       dispersion="glmGamPoi")
                         @test rows.status[1] == status
-                        @test isnothing(rows.pvalue[1])
-                        @test occursin("dispersion", rows.note[1])
+                        @test isnothing(rows.pvalue[1]) == !tested
+                        @test isnothing(rows.theta[1])
+                        @test occursin("overdispersion", rows.note[1])
                         @test all(rows.status[2:end] .== "ok")
                     end
                 finally
                     reval(raw"""
                         assignInNamespace("glm_gp", test_original_gp, ns="glmGamPoi")
-                        rm(test_original_gp, test_ql)
+                        rm(test_original_gp, test_od)
                     """)
                 end
             end
         end
     end
 
+    @testset "Condition labels never change the fits or reach R as code" begin
+        if !HAVE_GLMGAMPOI
+            # Issue #37's contract: these comparisons fail, not skip, in CI.
+            @test !_GLMGAMPOI_IN_CI
+            @test_skip false
+        else
+            # Condition labels join group, run and sub-group names with '/', and
+            # those names may carry '-' and '.'. test_de parses and evaluates a
+            # text contrast as R code, so the contrast goes to it as numeric
+            # weights. These labels, and one that reads as R code, must give the
+            # A/B results under both dispersions, and that code must never run.
+            # The coded label opens with '[' so that the text the old path built,
+            # "da_g[...]", indexed the live factor da_g and so ran the code.
+            off = log.(tss_factors(effect, effect_samples))
+            evaluated = "mm_contrast_label_evaluated"
+            labels = (("run-A", "run-B"), ("grpA/run_1", "grpB/run_2"), ("2024-01", "2024-02"),
+                      ("A", "[assign(\"$evaluated\", TRUE, envir = globalenv())]"))
+            for dispersion in ("glm_nb", "glmGamPoi")
+                base = fit_nb(effect, effect_groups, off; levels=("A", "B"), dispersion)
+                @test all(base.status .== "ok")
+                for (ref, con) in labels
+                    groups = [g == "A" ? ref : con for g in effect_groups]
+                    rows = fit_nb(effect, groups, off; levels=(ref, con), dispersion)
+                    @test rows.status == base.status
+                    for col in (:estimate, :se, :statistic, :pvalue, :theta)
+                        @test all(isapprox.(rows[!, col], base[!, col]; rtol=1e-10))
+                    end
+                end
+            end
+            @test !rcopy(R"exists($evaluated, envir = globalenv())")
+        end
+    end
+
     @testset "Degenerate and filtered taxa are reported, not dropped" begin
         if !HAVE_MASS
             @test_skip false