Fix over-rejecting p-values for penalized terms (Wood 2013 Tr statistic, fixes #163)

pygam/pygam.py

@@ -1252,47 +1252,60 @@ def _compute_p_value(self, term_i):
 
         Notes
         -----
-        Wood 2006, section 4.8.5:
-            The p-values, calculated in this manner, behave correctly for un-penalized
-            models, or models with known smoothing parameters, but when smoothing
-            parameters have been estimated, the p-values are typically lower than they
-            should be, meaning that the tests reject the null too readily.
+        Implements the Tr test statistic from Wood (2013), "On p-values for
+        smooth components of an extended generalized additive model"
+        (Biometrika 100(1), 221-228), as used by `mgcv::summary.gam`.
 
-                (...)
+        For a fitted GAM term with coefficient subvector beta and Bayesian
+        posterior covariance V, the rank-r pseudoinverse of V is built from
+        the top-r eigencomponents, where r is the integer closest to the
+        term's effective degrees of freedom. The statistic
 
-            In practical terms, if these p-values suggest that a term is not needed in
-            a model, then this is probably true, but if a term is deemed ‘significant’
-            it is important to be aware that this significance may be overstated.
+            T_r = beta^T V^{-r} beta
 
-        based on equations from Wood 2006 section 4.8.5 page 191
-        and errata https://people.maths.bris.ac.uk/~sw15190/igam/iGAMerrata-12.pdf
-
-        the errata show a correction for the f-statistic.
+        is referenced against a chi-square distribution with r degrees of
+        freedom. Using r = round(edof) corrects the over-rejection seen
+        when smoothing parameters are estimated (issue #163).
         """
         if not self._is_fitted:
             raise AttributeError("GAM has not been fitted. Call fit first.")
 
         idxs = self.terms.get_coef_indices(term_i)
         cov = self.statistics_["cov"][idxs][:, idxs]
-        coef = self.coef_[idxs]
+        beta = self.coef_[idxs]
+
+        # term edof; fall back to nominal length when edof_per_coef is not
+        # available for every coefficient (e.g. for the intercept term, or
+        # when there are more splines than samples).
+        edof_per_coef = self.statistics_["edof_per_coef"]
+        if len(edof_per_coef) >= max(idxs) + 1:
+            edof_term = edof_per_coef[idxs].sum()
+        else:
+            edof_term = len(idxs)
 
-        # center non-intercept term functions
-        if isinstance(self.terms[term_i], SplineTerm):
-            coef -= coef.mean()
+        # eigendecomposition of the (symmetric) covariance
+        eig_vals, eig_vecs = np.linalg.eigh(cov)
+        order = eig_vals.argsort()[::-1]
+        eig_vals = eig_vals[order]
+        eig_vecs = eig_vecs[:, order]
 
-        inv_cov, rank = sp.linalg.pinv(cov, return_rank=True)
-        score = coef.T.dot(inv_cov).dot(coef)
+        # Wood (2013): truncate to rank = round(edof) for the term
+        rank = int(round(edof_term))
+        rank = max(1, min(rank, len(idxs)))
 
-        # compute p-values
-        if self.distribution._known_scale:
-            # for known scale use chi-squared statistic
-            return 1 - sp.stats.chi2.cdf(x=score, df=rank)
-        else:
-            # if scale has been estimated, prefer to use f-statistic
-            score = score / rank
-            return 1 - sp.stats.f.cdf(
-                score, rank, self.statistics_["n_samples"] - self.statistics_["edof"]
-            )
+        # numerical floor: only keep eigenvalues that are clearly nonzero
+        tol = np.max(np.abs(eig_vals)) * np.finfo(float).eps * len(eig_vals) * 100
+        rank = min(rank, int(np.sum(eig_vals > tol)))
+
+        if rank == 0:
+            return 1.0
+
+        inv_eig = np.zeros_like(eig_vals)
+        inv_eig[:rank] = 1.0 / eig_vals[:rank]
+        cov_inv = eig_vecs @ np.diag(inv_eig) @ eig_vecs.T
+
+        stat = beta.T @ cov_inv @ beta
+        return stats.chi2.sf(stat, df=rank)
 
     def confidence_intervals(self, X, width=0.95, quantiles=None):
         """Estimate confidence intervals for the model.

pygam/tests/test_pvalue.py

@@ -0,0 +1,160 @@
+import warnings
+
+import numpy as np
+
+from pygam import LinearGAM, f, s
+
+# --- DATA GENERATORS ---
+
+
+def gen_noise_univariate(rng):
+    """Standard Univariate Noise"""
+    n = 200
+    X = rng.standard_normal((n, 1))
+    y = rng.standard_normal(n)
+    return X, y
+
+
+def gen_strong_signal(rng):
+    """Strong Sine Wave"""
+    n = 200
+    X = rng.standard_normal((n, 1))
+    y = np.sin(X * 2).flatten() + rng.standard_normal(n) * 0.5
+    return X, y
+
+
+def gen_multivariate_noise(rng):
+    """Two uncorrelated noise predictors"""
+    n = 200
+    X = rng.standard_normal((n, 2))
+    y = rng.standard_normal(n)
+    return X, y
+
+
+def gen_mixed_data(rng):
+    """Continuous + Categorical (3 levels)"""
+    n = 200
+    x_cont = rng.standard_normal((n, 1))
+    x_cat = rng.choice([0, 1, 2], size=(n, 1))
+    X = np.hstack([x_cont, x_cat])
+    y = rng.standard_normal(n)  # Pure noise response
+    return X, y
+
+
+def gen_collinear_data(rng):
+    """Two highly correlated predictors"""
+    n = 200
+    x1 = rng.standard_normal((n, 1))
+    x2 = x1 + rng.standard_normal((n, 1)) * 0.001  # Almost identical
+    X = np.hstack([x1, x2])
+    y = rng.standard_normal(n)
+    return X, y
+
+
+# --- TEST ENGINE ---
+
+N_SIMS = 100
+SEED = 12345
+FPR_BOUNDS = (1.0, 10.0)  # window around the 5% nominal level
+
+
+def run_test_scenario(data_gen_func, gam_factory, term_idx=0, n_sims=N_SIMS, seed=SEED):
+    """Run a simulation loop and return the rejection rate (%) for the term.
+
+    The RNG is seeded so that the loop is reproducible across runs.
+    Fits that fail to converge are skipped; if every fit fails the test
+    will raise rather than silently report 0%.
+    """
+    rng = np.random.default_rng(seed)
+    rejections = 0
+    n_fitted = 0
+
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        for _ in range(n_sims):
+            X, y = data_gen_func(rng)
+            try:
+                gam = gam_factory().fit(X, y)
+            except Exception:  # noqa: S112, BLE001
+                continue
+            n_fitted += 1
+            if gam.statistics_["p_values"][term_idx] < 0.05:
+                rejections += 1
+
+    if n_fitted == 0:
+        raise RuntimeError("no GAM fits succeeded in this scenario")
+
+    return (rejections / n_fitted) * 100
+
+
+# --- PYTEST CASES ---
+
+
+def test_fpr_noise_unpenalized():
+    """Univariate noise, lam=0. FPR should be near 5%."""
+    factory = lambda: LinearGAM(s(0, n_splines=10), lam=0)
+    rate = run_test_scenario(gen_noise_univariate, factory)
+    assert FPR_BOUNDS[0] <= rate <= FPR_BOUNDS[1], (
+        f"FPR {rate:.1f}% outside {FPR_BOUNDS}"
+    )
+
+
+def test_fpr_noise_smoothed():
+    """Univariate noise, lam=0.6. FPR should be near 5%."""
+    factory = lambda: LinearGAM(s(0, n_splines=10), lam=0.6)
+    rate = run_test_scenario(gen_noise_univariate, factory)
+    assert FPR_BOUNDS[0] <= rate <= FPR_BOUNDS[1], (
+        f"FPR {rate:.1f}% outside {FPR_BOUNDS}"
+    )
+
+
+def test_power_strong_signal():
+    """Strong signal. Power should be near 100%."""
+    factory = lambda: LinearGAM(s(0, n_splines=10), lam=0.6)
+    rate = run_test_scenario(gen_strong_signal, factory)
+    assert rate >= 95.0, f"Power {rate:.1f}% is too low (Target >=95%)"
+
+
+def test_multivariate_term0():
+    """Multivariate s(0) + s(1), check term 0."""
+    factory = lambda: LinearGAM(s(0) + s(1), lam=0.6)
+    rate = run_test_scenario(gen_multivariate_noise, factory, term_idx=0)
+    assert FPR_BOUNDS[0] <= rate <= FPR_BOUNDS[1], (
+        f"FPR {rate:.1f}% outside {FPR_BOUNDS}"
+    )
+
+
+def test_multivariate_term1():
+    """Multivariate s(0) + s(1), check term 1."""
+    factory = lambda: LinearGAM(s(0) + s(1), lam=0.6)
+    rate = run_test_scenario(gen_multivariate_noise, factory, term_idx=1)
+    assert FPR_BOUNDS[0] <= rate <= FPR_BOUNDS[1], (
+        f"FPR {rate:.1f}% outside {FPR_BOUNDS}"
+    )
+
+
+def test_mixed_types():
+    """Spline + factor, check the spline term."""
+    factory = lambda: LinearGAM(s(0) + f(1), lam=0.6)
+    rate = run_test_scenario(gen_mixed_data, factory, term_idx=0)
+    assert FPR_BOUNDS[0] <= rate <= FPR_BOUNDS[1], (
+        f"FPR {rate:.1f}% outside {FPR_BOUNDS}"
+    )
+
+
+def test_high_complexity():
+    """High spline count (n_splines=25), univariate noise."""
+    factory = lambda: LinearGAM(s(0, n_splines=25), lam=0.6)
+    rate = run_test_scenario(gen_noise_univariate, factory)
+    assert FPR_BOUNDS[0] <= rate <= FPR_BOUNDS[1], (
+        f"FPR {rate:.1f}% outside {FPR_BOUNDS}"
+    )
+
+
+def test_collinearity():
+    """Two near-identical predictors. Convergence is allowed to fail; if it fits, FPR should be calibrated."""
+    factory = lambda: LinearGAM(s(0) + s(1), lam=0.6)
+    rate = run_test_scenario(gen_collinear_data, factory, term_idx=0)
+    assert FPR_BOUNDS[0] <= rate <= FPR_BOUNDS[1], (
+        f"FPR {rate:.1f}% outside {FPR_BOUNDS}"
+    )