p-values are too small

[moptis]

Hello,

Is there a plan to output p-values or f-stats to the individual predictors?

Thanks!

[dswah]

Hi @moptis . Yes! the plan is to add these asap!

Also, I happily welcome collaborators. :wink: Do you have experience with these tests?

[dswah]

hi @moptis i've just added p-values to the model summary!

you should be able to see the following for a fitted GAM.

plz let me know how this works for you.

[moptis]

@dswah This is great! And sorry for the late reply.

I gave it a test run and I'm getting the same p-values for 3 different features (I assure you they are different features!). See the snapshot below. Any ideas on what's gone wrong? Note that I normalized 'y' to get p-values that weren't equal to one, but normalizing X didn't make a difference.

Thanks for all the work here!

[dswah]

@moptis 😒 that does seem weird. also strange that the intercept is the same. it looks like all the pvalues are getting as small as possible... im going to dig in a little more.

what happens when you dont normalize y?

and could i ask you to add some noise (maybe a lot?) to one of the features and show me the pvalues again?

lets fix this.

[moptis]

@dswah if I don't normalize y, all the p-values are equal to 1. If i use different combinations of features (e.g. x1, x2, x3; x1, x2; x1; etc.) I get the same result: the same low p-value for all features if 'y' is normalized and 1 if 'y' is not normalized.

I wonder if I'm using some package dependency that is out of date?

[moptis]

Also for reference is a snapshot of the result in R using the 'mgcv' library.

[dswah]

@moptis thanks for this info. it is really helpful.

I was able to reproduce your findings.

When i divide y by its standard deviation, then all my pvalues become tiny and equal, just like yours. When i multiply y by a large number, then i get all the pvalues are equal to 1.

I hope this is some obvious bug in the math.

[dswah]

there is at least one obvious bug in the math :) https://github.com/dswah/pyGAM/pull/171

[dswah]

@moptis I've pushed a new version of the code that fixes a bug in the F-statistic (used in the p-value).

could you update the library and show me the results again?

[moptis]

@dswah Sorry for the delay. I was out of office over the last 10 days. I'll update the library soon and will reach out when I have results.

[dswah]

@moptis have you had a chance to try this out again?

[moptis]

@dswah I'm still getting the same result with version 0.5.2:

[dswah]

@moptis wow Im so sorry i havent answered yet! let me take another look!

[dswah]

@moptis grrr... hmm the p-values certainly look wrong, but are they still scaling wrt the data scale?

[alexgain]

@dswah Would just like to mention that the p-values are still broken I believe, as I am also getting p_values of 1.11e-16 myself in a totally different context.

[dswah]

@alexgain i believe it :/ I frequently see these values as well.

I'm not sure what to do until we fix this bug.

• Should we remove the p-values from the model summary?
• should we leave them in, but add a warning that there is a known bug in the computation?
• something else?

your ideas are very welcome

[moptis]

@dswah @alexgain I'm not so sure these are bugs. I think 1e-16 is just the lowest possible value. I see the same in the R mgcv package, except the value is 2e-16. So if I see 1e-16, I just know that the feature is VERY statistically significant.

[alexgain]

@moptis @dswah I think that is the case and it is reaching machine precision, however I believe the p-values are much smaller than they should be. For a test case, I tried an expression where I would expect a significantly larger p-value. Indeed, it was much larger, but was still on the order of 1e-05 and again decreased over time until it was 1e-16.

I will be testing this against the R implementation and report back the results in a couple weeks or so.

I'm not sure what to do until we fix this bug.

I would go with the second bullet-point (just give a warning) until it is fixed.

[canthonyscott]

Hi all, Is there any new information regarding this? I am about to perform a large analysis and would love to use pyGAM but I do need the p-values for each feature in the model. Are they trust-worthy? Is it save to move forward?

Thanks for the guidance!

[alittlesliceoftom]

I don't think so but just a fellow user atm, easiest test is to run the same model through MGCV in R.

[shyamcody]

@dswah probably a reason the f-statistics are shooting higher values is that maybe the wald statistics being calculated is wrong according to this cross-validated answer. To summarize, the wald statistics we calculate is beta inv(V) Transpose(beta). But with rank r in the penalized model, Simon wood used beta((V)^-rank)Transpose(beta) which is lesser value in case of rank>1; and therefore will return higher p-value as the corresponding CDF will be less. One more thing: we should probably look into how the covariance matrix is getting calculated. Cause it is having near-zero determinant in a number of my test cases.

[dswah]

@shyamcody nice find. I took a look at the cross-validated answer, and it looks like we accounted for the rank of the covariance matrix in the F-statistic.

(also, I believe we are correctly following the updated formula from the errata as shown in equation for p.195.)

Then again, I am not sure what Wood means by "Vr- is the rank-r pseudoinverse of the covariance matrix." In the code I have simply calculated the pseudo-inverse and collected the rank returned by the scipy's pinv, which I pass to the corresponding statistic...

Perhaps like you suggest we need to somehow scale the pseudoinverse by the computed rank? I am looking at mgcv source code to try to get an idea of how they compute their p-values, but I haven't understood much yet.

With respect to poor conditioning of the determinant matrix, let's look into this! Perhaps the un-centered features are causing problems in the covariance estimates?

[shyamcody]

hey @dswah, my bad in reading. I misunderstood the r-- notation as the power of (-r) from the StackOverflow answer. Can you exactly mention how are you creating the parameter covariance matrix or refer to where in Simon's book it is written? my idea about the covariance matrix calculation is not clear and therefore comments may be insignificant. And yes, I too believe the unpenalized model related implementation is correct.

[dswah]

@shyamcody here is the seemingly innocent line where the covariance matrix is computed.

The B matrix is computed in the PIRLS loop on https://github.com/dswah/pyGAM/blob/b57b4cf8783a90976031e1857e748ca3e6ec650b/pygam/pygam.py#L756 which corresponds to equation (4.32) on pg 184

The action happens between pg 184 and 185 of the 1st edition book.

[uditrana]

Has this been fixed? If the P scores are wrong than the significance codes are also extremely wrong, right?

[shyamcody]

The issue is that, me and dswah were not able to find out what exactly is wrong in the code. Can you figure it out?

On Wed, Aug 11, 2021 at 3:45 AM Udit Ranasaria @.***> wrote:

Has this been fixed? If the P scores are wrong than the significance codes are also extremely wrong, right?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/dswah/pyGAM/issues/163#issuecomment-896347926, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALED6RLMLRE45FXNQ7TB6M3T4GQGJANCNFSM4EXFWVGA . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&utm_campaign=notification-email .

[uditrana]

I investigated slightly but unfortunately do not have the domain expertise nor the time to investigate. I was just experimenting lightly with GAMs anyways. Might have to switch to R and mgcv if I end up wanting to use a GAM for my model...

I also noticed that the ordering in which you passed in features seemed to change results (not sure if this is related or not).

[jonathan-taylor]

Hmm... this line is a bit curious -- why center here? https://github.com/dswah/pyGAM/blob/b57b4cf8783a90976031e1857e748ca3e6ec650b/pygam/pygam.py#L1259

Perhaps this is to try to separate into linear and non-parametric parts but I don't think this centering does the trick.

[jonathan-taylor]

For models with linear terms only and lam=0, the p-values do agree with statsmodels.OLS.

Including a factor the p-values do not agree.

One thing I do see is that you have a rank-deficient covariance matrix and are trying to compute \chi^2 for inestimable contrasts. E.g. for a factor the rows for the indicator encoding are not in the rowspace of the design... Can't say for sure this would make p-values too small always but inference will be affected.

This reference also probably useful: https://www.jstor.org/stable/43304548?seq=1#metadata_info_tab_contents

[Sorkanius]

any plan for updating these issues? Looking for a solid python GAM library, such as R's mgcv