2
$\begingroup$

Is there a hypothesis test that's ideally uniformly most powerful or metric that tells us whether we should use GAM vs GLM?

Does there exists some kind of metric i.e. AIC/BIC or loglikelihood to compare whether GAM vs GLM should be used?

May Cross-Validation be used to determine whether GAM or GLM should have been used?

$\endgroup$
1
  • 2
    $\begingroup$ A GAM is a GLM, just with a sophisticated predictor-generator (using base-splines). And yes, AIC/BIC, logLik and cross-validation are all fine to compare GAM and GLM. Note that mgcv::gam uses generalised cross-validation for fitting (as default), not maximum likelihood, but that shouldn't keep you. Just don't be surprised if the same model will be estimated slightly differently. $\endgroup$
    – Carsten
    Commented Feb 15, 2020 at 22:10

1 Answer 1

6
$\begingroup$

Simon Wood, the author of the mgcv package for R and a statistician who has made significant contributions to GAM theory and methods, has developed well-performing Wald-like tests for smooths. As such, you could fit a model with a linear term and a smooth with no null space and then use the test on the smooth term as a test for non-linearity (on the link scale).

m <- gam(y ~ x + s(x, m=c(2,0)), method="REML")

In the model above we ask for a linear term of x, plus a smooth of x. The m argument specifies that we want the usual penalty on the 2nd derivative of the smooth (measuring wiggliness) but we want a basis with no (0) null space. The null space of the basis is all the functions for which the wiggliness penalty doesn't apply, which includes the flat constant function and the linear function. By removing these two functions from the basis (usually only the flat, constant function is removed as it is confounded with the model intercept) the smooth represents the non-linear aspects of the model fit. Hence if it is significant that suggest non-linearity on the link scale which would not be modelled by the corresponding GLM.

Here's an example from ?mgcv::tprs

library("mgcv")
n <- 100
set.seed(2)
x <- runif(n)
y <- x + x^2*.2 + rnorm(n) *.1
df <- data.frame(x = x, y = y)

We're simulating data with a linear effect of x plus a weak quadratic effect of x, to simulate some small non-linearity. So, is the smooth significantly different from straight line?

m1 <- gam(y ~ x+ s(x, m=c(2, 0)), data = df, method = "REML")
summary(m1) ## not quite

The final output is:

Family: gaussian 
Link function: identity 

Formula:
y ~ x + s(x, m = c(2, 0))

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.02133    0.05315  -0.401    0.689    
x            1.18249    0.10564  11.193   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
        edf Ref.df     F p-value  
s(x) 0.9334      8 0.304  0.0781 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =   0.91   Deviance explained = 91.1%
-REML = -70.567  Scale est. = 0.012767  n = 100

We see that the smooth isn't quite significant at the 0.05 level (note the p values are approximate anyway, just like in a GLM, just a little more so as we chose the smoothness parameter too).

This compares with the model containing a smooth that has a null space (i.e. it contains the linear term). Now we see that the full smooth is significantly different from 0, but we cannot separate the linear and non-linear parts of the smooth:

## is smooth significantly different from zero?
m2 <- gam(y ~ s(x), data = df, method = "REML")
summary(m2) ## yes!

Family: gaussian 
Link function: identity 

Formula:
y ~ s(x)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.5600     0.0113   49.56   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
       edf Ref.df     F p-value    
s(x) 1.933  2.417 413.1  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =   0.91   Deviance explained = 91.1%
-REML = -69.353  Scale est. = 0.012767  n = 100

Notice that these models are essentially the same; the first is just separating out the linear basis function from the smooth basis which we add back in as an explicit linear parametric term.

In this case there isn't much to be gained from the GAM, but the true data generating process would be non-linearity, i.e. the GAM. With more data we'd be able to detect the non-linearity more easily, but this shows the problem of selecting between a GLM and a GAM. As such, unless you have a strong a prior guide that the effect is linear (on the link scale), then you should use a GAM and trust in the wiggliness penalty to decide if we need the non-linear bits.

As such you could just use the second model form (m2). You only need to decompose (m1) if you want an explicit test on the non-linear component.

Notice also that I use method = "REML"; choosing smoothness parameters using REML or ML selection is preferred today as GCV has been shown to undersmooth in circumstances where the GCV criterion is flat around the minimum.

$\endgroup$
2
  • $\begingroup$ This is a very informative post. Thank you. I am wondering if the ideas here can be generalized, however. For example, how might you specify testing the signficance of the nonparametric component above say a quadratic parametric model? Would you increase the order of the derivative penalty while retaining a basis with no null space in an analogous way? $\endgroup$
    – aphe
    Commented Oct 6, 2023 at 13:42
  • $\begingroup$ I am away from keyboard just now, but I’m pretty sure I’ve answered that somewhere before… I think it you were to set m=c(3,0) then that would mean constant, linear, & quadratic functions are in the span of the null space and thus removes from that basis. But I’d need to check and the devil might be in the detail. $\endgroup$ Commented Oct 7, 2023 at 8:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.