Is there a hypothesis test that tells us whether we should use GAM vs GLM?

Question

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?

Answer 1

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.