From ed0605c4b072aeb2b11639125fe2cbb64c07fb0b Mon Sep 17 00:00:00 2001 From: "Jonathan D.A. Jewell" Date: Mon, 19 Oct 2026 12:43:59 +0100 Subject: [PATCH 1/3] Report glm_gp's negative-binomial overdispersion as theta The glmGamPoi path reported 1 / ql_disp_shrunken as the dispersion theta, but that is the quasi-likelihood scale test_de divides by, not the negative-binomial parameter MASS::glm.nb's theta inverts. Read glm_gp's overdispersions element instead (variance = mu + mu^2 * overdispersion) and report its inverse, so dispersion_theta means one thing under both dispersion settings. The Poisson boundary moves with it: a taxon whose negative-binomial overdispersion is zero is marked boundary, as under glm_nb, and keeps its statistics, because the F-test's denominator is the shrunken quasi-likelihood dispersion and stays defined there. The earlier withholding was justified only by a zero quasi-likelihood scale, which this check no longer reads. A negative or non-finite overdispersion still fails only its own taxon. theta is nothing where it is unbounded. The direct R reference in the tests inverts gp$overdispersions, and the isolation test overrides that element, which test_de reads neither for its refit nor for its F ratio, so the test under comparison is unchanged. Requested in review on #62. --- config/defaults/pipeline.yml | 7 ++--- src/analysis/differential.jl | 47 +++++++++++++++++++--------------- test/unit/test_differential.jl | 31 ++++++++++++++-------- 3 files changed, 51 insertions(+), 34 deletions(-) diff --git a/config/defaults/pipeline.yml b/config/defaults/pipeline.yml index 2181a23..debcff0 100644 --- a/config/defaults/pipeline.yml +++ b/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 diff --git a/src/analysis/differential.jl b/src/analysis/differential.jl index 5e46c3f..0e47773 100644 --- a/src/analysis/differential.jl +++ b/src/analysis/differential.jl @@ -349,20 +349,18 @@ 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 not + # the quasi-likelihood scale the F-test divides by. Zero overdispersion + # is the Poisson boundary: the row is marked, as under glm_nb, and keeps + # its statistics, since the F-test's denominator is the shrunken + # quasi-likelihood dispersion and stays defined. 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 +374,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 negative-binomial overdispersion is zero: these counts show no overdispersion, so the fit is effectively Poisson and theta is unbounded" + } } 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 +455,20 @@ 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 same +parameter `MASS::glm.nb`'s theta inverts; it is `nothing` where that +overdispersion is zero. 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"`. @@ -501,7 +506,7 @@ function fit_nb(counts::AbstractMatrix{<:Integer}, groups::AbstractVector{<:Abst 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_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 diff --git a/test/unit/test_differential.jl b/test/unit/test_differential.jl index 46a6e79..6c4e801 100644 --- a/test/unit/test_differential.jl +++ b/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,37 +322,46 @@ 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 only what theta and the + # boundary check read, never the test itself. 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 From 8dac711a8eb8c33d61eba980f70bf74f60f3e000 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 20 Oct 2026 12:41:56 +0000 Subject: [PATCH 2/3] Pass the glmGamPoi contrast to test_de as numeric weights glmGamPoi::test_de parses a text contrast with rlang::parse_expr and evaluates it, so the coefficient name built from the contrast label reached R as code. Condition labels join group, run and sub-group names with "/", and those names may carry "-", so comparing "grpA/run_1" with "grpB/run_2", or "run-A" with "run-B", stopped test_de and failed every taxon, while glm_nb fitted the same labels. A label that reads as R code ran: with glmGamPoi 1.14.3 carrying the 1.22.0 test_de, a label of [file.create("")] created the file before test_de stopped. The design matrix is now built once, and test_de receives the numeric contrast that selects the contrast column, which it uses without parsing. For labels that were already syntactic the rows are identical. A test fits both dispersions under labels with "-", "/" and R code, and requires the A/B results and that the code never ran. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HtyPfZ1MWBB94RawgzQgz7 --- src/analysis/differential.jl | 35 ++++++++++++++++++++-------------- test/unit/test_differential.jl | 33 ++++++++++++++++++++++++++++++++ 2 files changed, 54 insertions(+), 14 deletions(-) diff --git a/src/analysis/differential.jl b/src/analysis/differential.jl index 0e47773..836d2b1 100644 --- a/src/analysis/differential.jl +++ b/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)) @@ -458,8 +463,9 @@ effect scale. In that path `statistic` is the quasi-likelihood F, `pvalue` its p-value, and `theta` the inverse of `glm_gp`'s negative-binomial overdispersion for the taxon (its `overdispersions` element), the same parameter `MASS::glm.nb`'s theta inverts; it is `nothing` where that -overdispersion is zero. Names are stripped and matched without regard to -case. +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 @@ -505,7 +511,8 @@ 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_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 diff --git a/test/unit/test_differential.jl b/test/unit/test_differential.jl index 6c4e801..783a03d 100644 --- a/test/unit/test_differential.jl +++ b/test/unit/test_differential.jl @@ -368,6 +368,39 @@ 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 From 10e099b6ba59cae3ab225d6b118b0da5bdfce374 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 20 Oct 2026 12:41:57 +0000 Subject: [PATCH 3/3] State what glmGamPoi's theta, standard error and boundary mean With overdispersion_shrinkage = TRUE, glm_gp returns each taxon's own overdispersion estimate in `overdispersions` but refits the coefficients at the cross-taxon trend, and test_de refits the reduced model at that trend and divides by the shrunken quasi-likelihood dispersion; only predict(), which supplies lfc_se, reads `overdispersions`. The comments, the boundary note, the isolation test's comment, the README and the result type now say so: theta is the taxon's estimate before shrinkage, and a zero estimate leaves the test's statistics standing because the test used the trend, so the fit there is not Poisson. The README lists dispersion with the other study-scoped analysis.differential settings. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HtyPfZ1MWBB94RawgzQgz7 --- README.md | 12 ++++++++++-- frontend/src/api/types.ts | 2 +- src/analysis/differential.jl | 26 ++++++++++++++------------ test/unit/test_differential.jl | 5 +++-- 4 files changed, 28 insertions(+), 17 deletions(-) diff --git a/README.md b/README.md index ecffd7d..d77f7c8 100644 --- a/README.md +++ b/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 diff --git a/frontend/src/api/types.ts b/frontend/src/api/types.ts index e6ee284..fff182b 100644 --- a/frontend/src/api/types.ts +++ b/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 diff --git a/src/analysis/differential.jl b/src/analysis/differential.jl index 836d2b1..f87d7bc 100644 --- a/src/analysis/differential.jl +++ b/src/analysis/differential.jl @@ -356,12 +356,14 @@ for (da_j in seq_len(da_n)) { # to the natural-log scale so every method reports one effect scale. The # 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 not - # the quasi-likelihood scale the F-test divides by. Zero overdispersion - # is the Poisson boundary: the row is marked, as under glm_nb, and keeps - # its statistics, since the F-test's denominator is the shrunken - # quasi-likelihood dispersion and stays defined. An unusable - # overdispersion fails only its own taxon. + # 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" @@ -382,7 +384,7 @@ for (da_j in seq_len(da_n)) { 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 negative-binomial overdispersion is zero: these counts show no overdispersion, so the fit is effectively Poisson and theta is unbounded" + 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 @@ -461,11 +463,11 @@ 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 `theta` the inverse of `glm_gp`'s negative-binomial -overdispersion for the taxon (its `overdispersions` element), 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. +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 diff --git a/test/unit/test_differential.jl b/test/unit/test_differential.jl index 783a03d..2c686a2 100644 --- a/test/unit/test_differential.jl +++ b/test/unit/test_differential.jl @@ -330,8 +330,9 @@ end 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 only what theta and the - # boundary check read, never the test itself. + # `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"))