Fix Wald thresholds for multi-coefficient and reversed contrasts

src/pydeseq2/glm.py

@@ -246,18 +246,24 @@ def wald_test(
     # Evaluate standard error and Wald statistic
     wald_se: float = np.sqrt(Hc.T @ M @ Hc)
 
+    # Apply the null and directional truncation to the tested scalar effect.
+    # Applying them coefficient-wise depends on the design parameterization.
+    contrast_lfc = contrast @ lfc
+
     def greater(lfc_null):
-        stat = contrast @ np.fmax((lfc - lfc_null) / wald_se, 0)
+        stat = np.fmax((contrast_lfc - lfc_null) / wald_se, 0)
         pval = norm.sf(stat)
         return stat, pval
 
     def less(lfc_null):
-        stat = contrast @ np.fmin((lfc - lfc_null) / wald_se, 0)
+        stat = np.fmin((contrast_lfc - lfc_null) / wald_se, 0)
         pval = norm.sf(np.abs(stat))
         return stat, pval
 
     def greater_abs(lfc_null):
-        stat = contrast @ (np.sign(lfc) * np.fmax((np.abs(lfc) - lfc_null) / wald_se, 0))
+        stat = np.sign(contrast_lfc) * np.fmax(
+            (np.abs(contrast_lfc) - lfc_null) / wald_se, 0
+        )
         pval = 2 * norm.sf(np.abs(stat))  # Only case where the test is two-tailed
         return stat, pval
 
@@ -276,7 +282,7 @@ def less_abs(lfc_null):
             "less": less(lfc_null),
         }[alt_hypothesis]
     else:
-        wald_statistic = float(contrast @ (lfc - lfc_null) / wald_se)
+        wald_statistic = float((contrast_lfc - lfc_null) / wald_se)
         wald_p_value = 2 * norm.sf(np.abs(wald_statistic))
 
     return wald_p_value, wald_statistic, wald_se

tests/test_wald_contrasts.py

@@ -0,0 +1,53 @@
+import numpy as np
+import pytest
+from scipy.stats import norm
+
+from pydeseq2.glm import wald_test
+
+
+@pytest.mark.parametrize(
+    "alternative", [None, "greater", "less", "greaterAbs", "lessAbs"]
+)
+@pytest.mark.parametrize("null", [0.0, 1.5])
+@pytest.mark.parametrize(
+    "contrast",
+    [[0.0, -1.0, 1.0], [0.0, 1.0, -1.0], [0.0, -1.0, 0.0], [0.0, 0.0, 1.0]],
+)
+def test_wald_tests_the_scalar_contrast(alternative, null, contrast):
+    # All tests concern the scalar effect c @ beta, including multi-coefficient
+    # and reversed contrasts. Thresholding individual coefficients is not invariant
+    # to the chosen model parameterization.
+    design = np.column_stack([np.ones(12), np.tile([0, 1, 0], 4), np.tile([0, 0, 1], 4)])
+    beta = np.array([4.0, 2.0, 3.0])
+    contrast = np.asarray(contrast)
+    mu = np.exp(design @ beta)
+    dispersion = 0.1
+    ridge = np.eye(3) * 1e-6
+    pvalue, statistic, se = wald_test(
+        design, dispersion, beta, mu, ridge, contrast, null, alternative
+    )
+    weight = mu / (1 + mu * dispersion)
+    information = design.T @ (weight[:, None] * design)
+    inverse = np.linalg.inv(information + ridge)
+    expected_se = np.sqrt(contrast @ inverse @ information @ inverse @ contrast)
+    effect = contrast @ beta
+    if alternative is None:
+        expected_stat = (effect - null) / expected_se
+        expected_p = 2 * norm.sf(abs(expected_stat))
+    elif alternative == "greater":
+        expected_stat = max((effect - null) / expected_se, 0)
+        expected_p = norm.sf(expected_stat)
+    elif alternative == "less":
+        expected_stat = min((effect - null) / expected_se, 0)
+        expected_p = norm.sf(abs(expected_stat))
+    elif alternative == "greaterAbs":
+        expected_stat = np.sign(effect) * max((abs(effect) - null) / expected_se, 0)
+        expected_p = 2 * norm.sf(abs(expected_stat))
+    else:
+        above = max((effect + abs(null)) / expected_se, 0)
+        below = min((effect - abs(null)) / expected_se, 0)
+        expected_stat = min(above, below, key=abs)
+        expected_p = max(norm.sf(abs(above)), norm.sf(abs(below)))
+    np.testing.assert_allclose(se, expected_se)
+    np.testing.assert_allclose(statistic, expected_stat)
+    np.testing.assert_allclose(pvalue, expected_p, atol=0)