From 2304a119c1c1e4a4b0384653472cad446d2ecea8 Mon Sep 17 00:00:00 2001 From: Vasanthkumar GR Date: Tue, 31 Mar 2026 19:31:13 +0530 Subject: [PATCH 1/5] fix: clip linear predictor in LogLink and LogitLink to prevent overflow --- pygam/links.py | 4 ++-- pygam/tests/test_links.py | 19 +++++++++++++++++++ 2 files changed, 21 insertions(+), 2 deletions(-) create mode 100644 pygam/tests/test_links.py diff --git a/pygam/links.py b/pygam/links.py index cc976eec..328fe726 100644 --- a/pygam/links.py +++ b/pygam/links.py @@ -118,7 +118,7 @@ def mu(self, lp, dist): ------- mu : np.array of length n """ - elp = np.exp(lp) + elp = np.exp(np.clip(lp, -500, 500)) return dist.levels * elp / (elp + 1) def gradient(self, mu, dist): @@ -178,7 +178,7 @@ def mu(self, lp, dist): ------- mu : np.array of length n """ - return np.exp(lp) + return np.exp(np.clip(lp, -500, 500)) def gradient(self, mu, dist): """ diff --git a/pygam/tests/test_links.py b/pygam/tests/test_links.py new file mode 100644 index 00000000..b3e2c10b --- /dev/null +++ b/pygam/tests/test_links.py @@ -0,0 +1,19 @@ +def test_log_link_no_overflow(): + """LogLink.mu should not overflow for large linear predictor values.""" + from pygam.links import LogLink + import numpy as np + link = LogLink() + lp = np.array([-1000.0, 0.0, 1000.0]) + result = link.mu(lp, dist=None) + assert np.all(np.isfinite(result)), "LogLink.mu produced inf or nan" + +def test_logit_link_no_overflow(): + """LogitLink.mu should not overflow for large linear predictor values.""" + from pygam.links import LogitLink + from pygam.distributions import BinomialDist + import numpy as np + link = LogitLink() + dist = BinomialDist() + lp = np.array([-1000.0, 0.0, 1000.0]) + result = link.mu(lp, dist) + assert np.all(np.isfinite(result)), "LogitLink.mu produced inf or nan" From ade5b27fcd3e2a2bd96705d5b143d8472604e1c7 Mon Sep 17 00:00:00 2001 From: Vasanthkumar GR Date: Tue, 31 Mar 2026 19:37:17 +0530 Subject: [PATCH 2/5] style: fix ruff formatting and imports in test_links.py --- pygam/tests/test_links.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/pygam/tests/test_links.py b/pygam/tests/test_links.py index b3e2c10b..ae5fc40f 100644 --- a/pygam/tests/test_links.py +++ b/pygam/tests/test_links.py @@ -1,17 +1,19 @@ +import numpy as np + +from pygam.distributions import BinomialDist +from pygam.links import LogitLink, LogLink + + def test_log_link_no_overflow(): """LogLink.mu should not overflow for large linear predictor values.""" - from pygam.links import LogLink - import numpy as np link = LogLink() lp = np.array([-1000.0, 0.0, 1000.0]) result = link.mu(lp, dist=None) assert np.all(np.isfinite(result)), "LogLink.mu produced inf or nan" + def test_logit_link_no_overflow(): """LogitLink.mu should not overflow for large linear predictor values.""" - from pygam.links import LogitLink - from pygam.distributions import BinomialDist - import numpy as np link = LogitLink() dist = BinomialDist() lp = np.array([-1000.0, 0.0, 1000.0]) From 434fff46d97df28e02636f0ad9eda8bacc13ed8f Mon Sep 17 00:00:00 2001 From: Vasanthkumar GR Date: Tue, 31 Mar 2026 19:49:34 +0530 Subject: [PATCH 3/5] test: set MPLBACKEND in conftest.py to prevent Tkinter Initialization --- pygam/tests/conftest.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/pygam/tests/conftest.py b/pygam/tests/conftest.py index 9a49a78a..60599dca 100644 --- a/pygam/tests/conftest.py +++ b/pygam/tests/conftest.py @@ -1,3 +1,6 @@ +import os +os.environ["MPLBACKEND"] = "Agg" + import pytest from pygam import ( From 09a7e3dd6d31abe5b51aa5097606f6e7836dd0a0 Mon Sep 17 00:00:00 2001 From: Vasanthkumar GR Date: Tue, 31 Mar 2026 19:53:29 +0530 Subject: [PATCH 4/5] style: run ruff format and revert conftest.py --- gen_imgs.py | 1 - pygam/tests/conftest.py | 3 --- test_output.txt | Bin 0 -> 41234 bytes 3 files changed, 4 deletions(-) create mode 100644 test_output.txt diff --git a/gen_imgs.py b/gen_imgs.py index 4e9d3ca0..6fc8694d 100644 --- a/gen_imgs.py +++ b/gen_imgs.py @@ -323,7 +323,6 @@ def chicago_tensor(): def expectiles(): - """Generate expectiles visualization.""" X, y = mcycle(return_X_y=True) diff --git a/pygam/tests/conftest.py b/pygam/tests/conftest.py index 60599dca..9a49a78a 100644 --- a/pygam/tests/conftest.py +++ b/pygam/tests/conftest.py @@ -1,6 +1,3 @@ -import os -os.environ["MPLBACKEND"] = "Agg" - import pytest from pygam import ( diff --git a/test_output.txt b/test_output.txt new file mode 100644 index 0000000000000000000000000000000000000000..f2aa924a96dc14d66b0294910fc90dea64b4f26d GIT binary patch literal 41234 zcmeI5eQz8`a>nQH0QnA^GXnE2CrYH=jpzVdmJQ*2vLI2(0jI-4A}NU{ie$KyWZflS zJxTrg*<$z1^vq80?ohO0ESB8edFiLCtE-->ewqLK-+R@vuSs=Q?N>AVJF&lfnpFqY zadi~FpH(laQ#&@ZZ_Ac&ulncepQ@AU(C&3sy|R0sR`2Y7AFZ6D>Uy=SWqPO_SJ;`@ar<53MxH;fj-RzrWcpyLQh_dmme!4$oT%3DPRFaK)}2{b16(ZD%6U zYv=brv}<$^&VC(gv1w;a?7qi#%#TvfiCufk+IeDkM%R!bHIr6+va6Bn)o=_woF$*9 zR!8bJwKMkZ7!uvHGbFRG45e%D)Y=Rz`*yd}^Q)+@WKDftuOHhzpW5qzwfLm*`b+qJ zt@_vMuOlEtj}8n!Sa#1*q7%EOy=d63wGusH5yGue5I#2ff^TO3kIwO-zbE#~EMVMU zfNxSsh8KzO(9pd$tvjgRT6_1)j}{lg?NMZiojWx&`=+(h-o<#DTIdL`kAjePd}l9k zF4DF(GXF&~90%b$YfpOa)N~Hdp`Jan@7MrrHb<}<9K|E#WGsX?Mv>v3mHNSOfn_%! z1uwpH&!mgz$YXNp>HnArvBda+X^zY_cykmP_9OhS zAjhL%E7iu+z!~;}4y3ihqtXuaxfj*Tpb<@dP(5yahI*jMXp1~i6IaxSF%i~#e7*J2 zhAm@J)1;r}Q3~PhQDlhMnFhXcX6=j5o(8rOy>}i{lpmiyFG5(#`>Z93D24FOg~?EM zN-Zwyv`CIZc=y5(_E~{E){91X68i7Sf|iYKEIlN48G}XUdP=>2=cSh7bA;n?G~&p! z(0}nOikf7*@!nHA`cvrRr)IP9l=O5zvx%Ly-m=#il3axIB*Ys><`ows9EVFThY&l8 z*MXPh>g^Hh^U~@-9tb_9Rwu3cc%5ffSDpJNFn_W$ z^?btRtLaOt5sa^{fzd0pV=2{P96q1R8l~|DnM1r5z5T7dkIgnHB2~_?%j%B9oQI;l z3`ixp_o0Up4|ygM?*HCeA3KSyRxCg?rx>+MqT}#ra;NyQ)B=JF9RS&`p~PeG%YJBo z@}_S(`q?;4x*Ox8>c?v5(xb$NAx&C)TKz}pQRu}3vu=3P)2uYN*UG8gK0mQ?$?r?^ zlqDtuo^pVbdG18IOFjx+hA-e@&=zv(Sg<}D@cHm>2m1<>RrZ-2J$WDgW4~hTt&G0x z-866`+9}Nl=lsd`S~+@2%V2-8zTl=En%QdQ5ny5O10%-IaE-3;{8+xGx)TB>PtPbPiP)(C~kXTULR|TO*##pNj8+86lHU4*)jT@GJ!rKlJ~BkuZAF=&8@Gp zW**Vda{EDB_d;zM!6W9#YvdR)#1LD74M6U+MLKSE`B0)k@?x?I(xLn~8H{Ut6uv)1 zru7jXMg$lOq*v#qx(ts``DSi}k&zwZ2)8)sBJp6%o-2pCXbeh{90MV-k8&#XZTVJm zlC)4Ow7gZ6Kd z@)!u|*>FQf@Fh_+0R{kS<~G#htu{9OB81CC1&d_pBW7r(a=%S}GYq2;F2k}HAw>Te zy+gZ@9U6-N!fP?Ih92X$@Xz!&d3HXDYjaHU?($w^pe7=j&GUiGH-W{9M=|2ds4BQM zU#3q7B-dpkgS=dKtXw2G#FC2ph$W#XyFrpsdllixW+-nWyZ6#;?TeOfQ?6-c`;G4l3{uNIWI%arGx*j1DjB;7vcRNJI&Zz7e*m0XQFu- zJ_z%~d|n`tg>WiUx?FD-!g3Z@iqJFA@{4I57TND|gdQE{VvPIdH^eH*GQu_%l^=UV z*C!z?;hB91@g8x;6WQlgZH-%-62?}BQ2tnQg=I2^QYKf)umoW}Ly^}<7{sDzv8z>^ zzYAdrryPdx*f6$pV9nBpd>;oPqpPnBqd3EAx(LD>MQnD63_(Y<^2Tt2autk1Si%p7 zAY?{5d>qcjyR2{2+kPB`@66_WvNv9375RZx2ut|k5E)`q_U*Nb>|B!LAY`5-w)=&} zs_l%9mp$yi5SH-6VKO`p(RIuLF(V$0B*$OY?m}3?4~HQ97^Jtl4n{p(H$wWdM$?Jx z$)b^$cxAt$pPH>aRFaZm2|pYnLpUrWQIX-0=aFWXwzUwJ@Iz)PDYx{_`a#|=`#G>e zW=(1=g>^AQB^P_qI0&6NOhZUUGja)SQ?k%`VmHMlgbaz z0&23vi)E*NBKCE&Xtux*8TvUfcNmcl#}@nY&7$#rFTKwm`ZUy zV7>G{gmr$X5dlY9hR<8MgdYw;Se~76t!w|649gHQ--q?5Ff5qWo9s$d#@d?aX%x;^ zT)`*eEP<}RE`%lgFjsOe%al8>jSy+eCZ%Y7R8BL{(B}LY{7tiz!)M6k5DiC84uuF- zW8B!JMfxh)UK~RD9l0`ARBNP)xi~)m@gV%##4L?UU`Q`DRVExXoNiZT#B=BHWG-CSInRKjF`E;VIeoBA9e z)wP6iUZ%Q|L8Kot>xmg#tg6Dl#aY+b=C)`1dDoi#TAZt>|4iXx`VB#CA_ zlRID)0(OEtFnf;hIVXv9Zom9xnjvY&#z!OJ&k6!yH1i z^70^5?jMHnQ?bi~VAcb1(<`gT^4Nq`2uq@gei`~vjerx& zsh40yJqBVGtyWxM5#&$1>oEpGACZ!;NmpUhvtgd;I&69QiDifWX!)VX_Nvc^5W#gqaR*ixu)nnoU98DRLr;tKwSZ@x z*=37;g2_qdvZbm`30IAhDQTB_=0&^5lDR{d@@?5A`LVD@`KGoQ8T(pP<{;vG=@Bts zAREa#RkFL~!gKo(rz_9m zE>zeMb^vRZ-t=&}@6zFQB7}WKE3Mr*yR(O~>eK4ao#;V$r-@?U#^c_BU;gPc3^CSe4DS2 zfnzVpZ^i65ypfgLz7LuA?G$(0Gc_f;K7I>#W#4No5|RTka~+5kn$n3V`aY@6*M~5* z=f2)Sc@^#V7S}ot#VdsC?QU*)2<7OVzE8-Ijf%GR|J#h1{4!CB|I+fXc5b#nf7=gT$N+) zn7u((Pwam)=hWP3OrFX2U4bc)g+`W5=?Uk>BTV%H{ywywS7ntYQfy)zrUA(E?#kS< z6lbmirrokLOIIBC-WUMI6DF zMaz))utHzLIftXi!)s(hoP)^NzrV8jv!=B@Qz7m5d=|nIMpuRq`C=6$EwWaMe~EsP zEz)kbc^%`}2wi4yPjV#b{;Yk>Eus{{5|&veLq{ll+UeK#v=R?~J4lA(Akb}s_q zDtcJLGKb3$gkO7?)Z z@5jAvgo+rL%Mf#a%AGDdIwyoBY_ea5SeI!SA$n@#Ll12%`LT^CKe1hu7i(PFw(;{{ zY?S+2!XtUs-SW>*j9120z0c;eo>zG;W4y1P(Kt0JEc1%bkQ8TRHBM@CwEW2i4z$#I%NJGmEub^05|RIkP&v zsWz;Ya_a5EV_8_=t zZ|2+hKJDx^3deNL?q_ZUo%uP-h*hJh-VCXADmkR$t8D2xv=UJltcZ4_$1x`ri?C<* z5-0v33y<1|7TjV>FwcC^XrsI0DMcP(% zYqU{o5n4iR@=I%pEt{))$9h5ZVXnJ%c`TN2{d_pG@9el~pRrhq%hH)JFb!{sjc|qR zPTCuO8M4QRw(awoSy#73H*JUMv>k@XENvB-q@DuyaqdfX_Z=%kvV|qiZ0)G{?mc6f zv{}=u!zj7$<>g~|JP`7Rksn%L6?^!1s||DD$Vw4uY+7_nW;9)6&9mtL5Ns=ZQ=pf4 zOP7l7^T@8CAE0~gHR@-48Q6}q+%E%`sPWBnb8`MN+?|qq!ab8mw*=}(+qJ$|eHm&*4Z-Ai zy8lNh+{;Sfp6&D1u6@MOAMrjhc}BLx?(VX2M2^hhnV4L)8~>qwrXTIi85NoXuTM3F zS+@5>UgCjC>5=U}_9vTdyJ=^AZEwHxPOxT*3#5yBf|)+66k^_-e3%(XED>$yUbz1^ zR>D&&%gc;+K})$ta%86xk70@DIYOP3ti?+LOd&nW(r{~hmp%z}3>5<{lULke))io&gkIj&O$!B2Pj%@~;%y=30 ztKGkweJEZ>=0N-N(-LSc_T(Sc|77JKuVZB7t+h0xh3E*{!;DVP5XvNP)lM{m?!uZk80D`?1r=k>RFKP2WFmcu3;Z2J=7@GJ9q zJ8yICu0xKC8Zmy5E6LApmjybZQhoF2Pzj>GTJU!J7tB^+6;WKGbY4f8(YRPwO?ZP$M9@)L{dqgpg;IjR2b zyq>w1h&+U-oGhuO)l$#pI~OAXK|0_+BPw}O6HVt zHT(`N_y_y{((Z-F_+-avAKsrQhOV-$<@ ztMM_|T_i563m=&*Spvws5AZ~I`r$SED;v(U&B;^$x4lbWT=T$?fuDW>_8n_o>nbRlb%-dhTaH{PZ}gya zv|rdc{jYIQ-VD#nLzcG&p>(`mSjRzmJIt#G<((pwj63fw=-d*^+fk05{PR4Nspr+`bN|?#KiMW<{++p}dKT>R zY4x9Gqod5Qz<)GO{Gyalrdd+lg}c!nUK$ZG=y=L$hHE10C?~Q$yCXw;DYUh$o z6!xoTt!`QVGTzJYi505v?=6jEgg@<1KJv+R>+BbdfT1b%`9qqEDvfB)_jKgr_8n~K z?I8VMd(W@Er?2;1Y=5cu%Jy!_TgIZa!jep9%s`CYTke>3!B+j+TYl{=eZ6JQa1{2; z3ME!6JG&-_#hj>#c|mlDH<_L6N%ZoD@^1Q=6jR#MpOcwRQ~Mb4XMBi;Vi`E7MMkGS8bPOO5!XwIg2`N20Jm#98 Date: Tue, 31 Mar 2026 22:45:24 +0530 Subject: [PATCH 5/5] feat: replace pinv p-value rank with Wood 2013 eigendecomposition of Vb Use np.linalg.eigh to decompose the Bayesian posterior covariance Vb, then project the test statistic into directions with eigenvalue > tol, where tol = max_eigval * p * eps. This fixes rank ambiguity for heavily-penalized smooth terms and matches mgcv's summary.gam approach. Reference: Wood, S.N. (2013) A simple test for random effects in regression models. Biometrika, 100(4), pp.1005-1010. Tests: - test_pvalue_rejects_useless_feature (existing, passes) - test_pvalue_invariant_to_scale (existing, passes) - test_pvalue_wood2013_range (new: all p-values in [0,1]) --- pygam/pygam.py | 43 +++++++++++++++++++++++---------- pygam/tests/test_GAM_methods.py | 12 +++++++++ 2 files changed, 42 insertions(+), 13 deletions(-) diff --git a/pygam/pygam.py b/pygam/pygam.py index 6f86ae43..18c0385c 100644 --- a/pygam/pygam.py +++ b/pygam/pygam.py @@ -1252,6 +1252,16 @@ def _compute_p_value(self, term_i): Notes ----- + Uses the eigendecomposition-based approach from Wood (2013), which projects + the test statistic into the effective range space of the Bayesian covariance + Vb. This avoids rank ambiguity from pseudo-inverse tolerances for penalized + smooth terms, matching the behaviour of mgcv's summary.gam. + + References + ---------- + Wood, S.N. (2013) A simple test for random effects in regression models. + Biometrika, 100(4), pp.1005-1010. + 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 @@ -1261,36 +1271,43 @@ def _compute_p_value(self, term_i): (...) 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’ + 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. - - 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. """ 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] + Vb = self.statistics_["cov"][idxs][:, idxs] # Bayesian posterior covariance + coef = self.coef_[idxs].copy() - # center non-intercept term functions + # center non-intercept smooth term functions if isinstance(self.terms[term_i], SplineTerm): coef -= coef.mean() - inv_cov, rank = sp.linalg.pinv(cov, return_rank=True) - score = coef.T.dot(inv_cov).dot(coef) + # Wood (2013) eigendecomposition: project into effective range space of Vb. + # Only eigenvectors with eigenvalue > tol contribute a non-trivial direction; + # directions with near-zero eigenvalues are fully penalized to zero. + eigvals, eigvecs = np.linalg.eigh(Vb) + tol = eigvals.max() * len(eigvals) * np.finfo(float).eps + keep = eigvals > tol + rank = int(keep.sum()) + if rank == 0: + return 1.0 # fully penalized term: not significant + + # Score = coef^T Vb^{-1} coef restricted to the kept subspace + # Equivalent to ||V^{-1/2} coef||^2 in the range of Vb + sqrt_inv = eigvecs[:, keep] / np.sqrt(eigvals[keep]) + score = float(np.sum((sqrt_inv.T.dot(coef)) ** 2)) # 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) + return 1.0 - 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( + return 1.0 - sp.stats.f.cdf( score, rank, self.statistics_["n_samples"] - self.statistics_["edof"] ) diff --git a/pygam/tests/test_GAM_methods.py b/pygam/tests/test_GAM_methods.py index ee29de14..58c36c3e 100644 --- a/pygam/tests/test_GAM_methods.py +++ b/pygam/tests/test_GAM_methods.py @@ -433,6 +433,18 @@ def test_pvalue_rejects_useless_feature(wage_X_y): assert p_values[-2] > 0.5 # because -1 is intercept +def test_pvalue_wood2013_range(wage_X_y): + """ + Wood 2013: all p-values must be in [0, 1] for every term, + and the eigendecomposition rank must always be at least 1. + """ + X, y = wage_X_y + gam = LinearGAM(s(0) + s(1) + f(2)).fit(X, y) + for term_i in range(len(gam.terms)): + p = gam._compute_p_value(term_i) + assert 0.0 <= p <= 1.0, f"p-value out of [0,1] for term {term_i}: {p}" + + def test_fit_quantile_is_close_enough(head_circumference_X_y): """see that we get close to the desired quantile