Two different GAMs result in the same AIC - How is this possible?

Question

I have two different GAM models as follows, generated via 'mgcv' in R.

I run checks to ensure that the basis dimension is specified correctly via gam.check().

I also run summary() to get a basic report.

Now I want to select the best model via AIC. However, R finds that the two models are equally parsimonious - how can this be?

#### Linear effect of distance ####

nz_gam_linear <- gam(reducedhap ~ s(latitude, k = 20) + distances, data = final2) # here, distance is a linear effect
gam.check(nz_gam_linear) # based on p-value being large (not significant), k = 20 is sufficient 

Method: GCV   Optimizer: magic
Smoothing parameter selection converged after 6 iterations.
The RMS GCV score gradient at convergence was 1.206796e-06 .
The Hessian was positive definite.
Model rank =  21 / 21 

Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.

              k'  edf k-index p-value
s(latitude) 19.0 17.4    1.19    0.94


#### Smooth effect of distance ####

nz_gam_smooth <- gam(reducedhap ~ s(latitude, k = 20) + s(distances, k = 20), data = final2) # here, distance is a smmoth effect
gam.check(nz_gam_smooth) # based on p-value being large (not significant), k = 20 is sufficient

Method: GCV   Optimizer: magic
Smoothing parameter selection converged after 14 iterations.
The RMS GCV score gradient at convergence was 6.818248e-08 .
The Hessian was positive definite.
Model rank =  29 / 29 

Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.

               k'  edf k-index p-value
s(latitude)  19.0 17.4    1.19    0.94
s(distances)  9.0  1.0    1.00    0.36


#### Summarizing and choosing models ####

summary(nz_gam_linear) # linear effect of distance on genetic diversity is not significant
summary(nz_gam_smooth) # smooth effect of distance on genetic diversity is significant

nz_gam_linear$aic # AIC = -15.90926
nz_gam_smooth$aic # AIC = -15.90926

Any idea what is going on here?

Specifying a third model that excludes distance would be appropriate. Trying this gives me an AIC = -12.0363. Clearly, in comparison to the other models, it's not the best.

NB - I can't share the actual code (sorry) as it confidential until we get the manuscript submitted.

Can anyone weigh in here?

Answer 1

Thomas already provided the answer - I just wanted to add some further comments.

When we fit a gam(y ~ s(x1) + s(x2)) model, we assume that the effect of x1 on y (adjusted for the effect of x2) is smooth, possibly nonlinear and the effect of x2 on y (adjusted for the effect of x1) is smooth, possibly nonlinear.

The data will help determine the actual shape of each of these effects - but, before we even see the data, we keep an open mind towards the possibility that the effects may be nonlinear. There is no guarantee that they will be nonlinear though (e.g., it is possible that x1 has a linear effect on y after adjusting for the effect of x2; it is also possible that x2 has no effect on y after adjusting for the effect of x1).

The edf (effective degrees of freedom) value reported for each estimated effect can be compared against 1 to get an immediate indication of the shape of the estimated effects obtained after fitting the model to the data.

An edf of 1 would correspond to a linear effect. (An edf close to 1 would correspond to an approximately linear effect.)

An edf greater than 1 would correspond to a smooth, nonlinear effect - the farther the edf would be from 1, the more complex the shape of the smooth, nonlinear effect.

In your case, x1 = latitude and x2 = distance. The corresponding edf values for the estimated effects of x1 and x2 are:

s(x1)    edf = 17.4
s(x2)    edf = 1.0 

This suggests that the estimated effect of x1 = latitude is smooth and highly non-linear, whereas the estimated effect of x2 = distance is linear. You can visualize the shape of these effects by applying the plot() command to your model object:

plot(nz_gam_smooth)

So you can simplify your model from gam(y ~ s(x1) + s(x2)) to gam(y ~ s(x1) + x2). Whether or not you simplify your model further from gam(y ~ s(x1) + x2) to gam(y ~ s(x1)) depends on your research question. If you are interested in describing how both x1 and x2 affect y, then you should keep x2 in the model even if its p-value is not statistically significant at your chosen significance level.

Answer 2

In

               k'  edf k-index p-value
s(latitude)  19.0 17.4    1.19    0.94
s(distances)  9.0  1.0    1.00    0.36

the s(distances) term has 1 degree of freedom, so it looks like the optimal spline fit chosen by gam() was actually a straight line and the two models really are the same.