feat: add likelihood ratio test (LRT) to DeseqStats

[pankev-in]

Summary

Adds the likelihood ratio test (LRT) to DeseqStats, the Python equivalent of R DESeq2's DESeq(dds, test="LRT", reduced=...). This resolves the long-standing request in #176 (and supersedes the stalled #178).

DeseqStats gains two arguments:

• test: "Wald" (default, unchanged behaviour) or "LRT";
• reduced: the nested reduced model for the LRT, as a formula string (e.g. "~group") or an explicit design matrix (numpy array / DataFrame).

dds = DeseqDataSet(counts, metadata, design="~group + condition")
dds.deseq2()

ds = DeseqStats(dds, contrast=["condition", "B", "A"], test="LRT", reduced="~group")
ds.summary()   # LRT p-values in ds.results_df

Unlike #178 (which only ever removed a single coefficient, hard-coded df=1, and dropped lfcSE), this supports arbitrary nested reduced models, including df > 1 — e.g. jointly testing a multi-level factor or comparing nested models in general — which is the main reason to use the LRT over the Wald test.

Closes #176.

The test

For each gene the statistic is

Λ = 2 · ( ℓ_full − ℓ_reduced )

where ℓ is the negative-binomial log-likelihood at the maximum-likelihood fit of each model, both using the same gene-wise dispersions (estimated once from the full model, exactly as in DESeq2). Under H₀ that the extra full-model coefficients are zero, Λ ~ χ² with

df = (# full-model coefficients) − (# reduced-model coefficients)

and pvalue = chi2.sf(Λ, df). As in R DESeq2, the reported log2FoldChange/lfcSE still correspond to the requested contrast of the full model, while the p-value reflects the model comparison. Cooks outlier replacement is honoured by refitting the reduced model on the same imputed counts as the full model (R's counts(dds, replaced=TRUE)).

Validation against R DESeq2 (the math, proven)

The implementation is validated numerically against R DESeq2 1.46.0 on the package's synthetic dataset, mirroring how the existing Wald reference data is produced. Reference data is regenerable via tests/data/generate_lrt_reference.R.

Example — full ~group + condition vs reduced ~group (df=1), PyDESeq2 vs R:

gene	R log2FC	Py log2FC	R stat	Py stat	R pvalue	Py pvalue
gene1	0.7315	0.7315	6.3158	6.3170	1.197e-02	1.196e-02
gene2	0.5353	0.5353	12.6914	12.6914	3.673e-04	3.673e-04
gene3	-0.6737	-0.6737	5.3890	5.3889	2.026e-02	2.027e-02
gene4	-0.4235	-0.4235	15.8217	15.8210	6.960e-05	6.963e-05
gene5	0.5880	0.5880	14.7313	14.7332	1.240e-04	1.238e-04
gene6	-0.0195	-0.0195	0.0040	0.0040	9.498e-01	9.498e-01
gene7	0.1351	0.1351	0.8134	0.8132	3.671e-01	3.672e-01
gene8	-0.2716	-0.2716	4.2516	4.2507	3.921e-02	3.923e-02
gene9	-0.2140	-0.2140	2.5794	2.5801	1.083e-01	1.082e-01
gene10	0.3895	0.3895	2.5095	2.5095	1.132e-01	1.132e-01

Log fold changes agree to ~1e-6 relative error; statistics and p-values for genes with signal agree to ~1e-3. Across all validated scenarios, LFCs match R to ≤1e-5 and, for genes with signal, statistics/p-values match well within the suite's existing tolerance.

Tests

Added to tests/test_pydeseq2.py (all passing; full suite 74/74 green):

• test_lrt_matches_r — single-factor df=1, multi-factor df=1, and joint df=2, vs R;
• test_lrt_no_independent_filtering — isolates the raw LRT statistic/p-value vs R;
• test_lrt_with_outliers — Cooks outlier replacement + reduced-model refit vs R (and asserts the raw counts are not mutated);
• test_lrt_reduced_as_design_matrix — a reduced model given as a matrix equals the formula;
• test_lrt_statistic_definition — R-independent checks: statistic ≥ 0, pvalue == chi2.sf(stat, df), and Λ ≈ (Wald stat)² at df=1;
• test_lrt_lfc_matches_wald — LRT reports the same LFC as the Wald test;
• test_lrt_input_validation — informative errors for missing/invalid reduced, non-nested reduced, unknown test, etc.

Docs: an LRT subsection is added to the minimal pipeline example, and the reference data is documented in tests/data/README.md.

Notes

• Fully backwards compatible: test defaults to "Wald", and the Wald path is unchanged (reduced must be None for it).
• The reduced-model fit reuses the existing Inference.irls routine, so custom Inference backends work with no changes to the ABC.

🤖 Generated with Claude Code

[pankev-in]

Quick note on the failing CI: it's the mypy -p pydeseq2 step that's red, not the tests. The same step is currently failing on main (its scheduled runs are red on it, and the pytest step gets skipped because mypy exits non-zero first), so this looks like a stubs/dependency regression that hit the whole repo around Jul 24 rather than anything in this PR.

The LRT tests pass locally (74/74 for the full suite). I also tidied up the typing in the code this PR touches, so the branch actually brings the total mypy count down from 228 to 188 instead of adding to it.

I didn't want to fold a repo-wide typing/deps fix into a feature PR, but happy to open a separate one for the mypy step if that's helpful. Let me know how you'd prefer to handle it.