# Interpretation of gam model results

I am fitting some generalized additive models using the mgcv package in R, and I am wanting to test between two models; whether I can remove a term or not. I am, however, getting conflicting results (as far as I can tell).

A model, m1, with a smooth term for x added, appears to give a better fit in terms of $R^{2}_{adj}$, AIC, deviance explained, and when comparing the models using an F-test. However, the significance of the smooth term is not significant (nor is it when I added to the model as a linear covariate, instead of a spline).

Is my interpretation of the smooth terms tests in correct? As much as I could understand the help page, was that the tests are approximate, but there is quite a large difference here.

The model outputs

m1 <- gam(out ~ s(x) + s(y) + s(z), data=dat)
> summary(m1)
#
# Family: gaussian
#
# Formula:
# out ~ s(x) + s(y) + s(z)
#
# Parametric coefficients:
#               Estimate Std. Error t value Pr(>|t|)
# (Intercept) -7.502e-16  1.209e-01       0        1
#
# Approximate significance of smooth terms:
#        edf Ref.df     F  p-value
# s(x) 4.005  4.716 1.810    0.136
# s(y) 8.799  8.951 4.032 4.01e-05 ***
# s(z) 7.612  8.526 5.649 4.83e-07 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# R-sq.(adj) =  0.213   Deviance explained = 24.8%
# GCV = 6.9741  Scale est. = 6.6459    n = 455

> AIC(m1)
#[1] 2175.898

> m2 <- gam(out ~ s(y) + s(z), data=dat)
> summary(m2)
#
# Family: gaussian
#
# Formula:
# out ~ s(y) + s(z)
#
# Parametric coefficients:
#              Estimate Std. Error t value Pr(>|t|)
# (Intercept) 1.705e-15  1.228e-01       0        1
#
# Approximate significance of smooth terms:
#        edf Ref.df     F  p-value
# s(y) 8.726  8.968 5.137 6.78e-07 ***
# s(z) 8.110  8.793 5.827 1.55e-07 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# R-sq.(adj) =  0.187   Deviance explained = 21.7%
# GCV =  7.144  Scale est. = 6.8639    n = 455

> AIC(m2)
#[1] 2187.168

> anova(m1, m2, test="F")
# Analysis of Deviance Table
#
# Model 1: out ~ s(x) + s(y) + s(z)
# Model 2: out ~ s(y) + s(z)
#   Resid. Df Resid. Dev      Df Deviance      F    Pr(>F)
# 1    433.58     2881.6
# 2    437.16     3000.7 -3.5791   -119.1 5.0073 0.0009864 ***


> summary(m3 <- gam(out ~ s(x) + s(y) + s(z), data=dat, select=TRUE))

#Family: gaussian

#Formula:
#out ~ s(x) + s(y) + s(z)

#Parametric coefficients:
#              Estimate Std. Error t value Pr(>|t|)
#(Intercept) -1.588e-14  1.209e-01       0        1

#Approximate significance of smooth terms:
#       edf Ref.df     F  p-value
#s(x) 4.424      9 1.750  0.00161 **
#s(y) 8.260      9 3.623 5.56e-06 ***
#s(z) 7.150      9 5.329 4.19e-09 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#R-sq.(adj) =  0.212   Deviance explained = 24.7%
#GCV = 6.9694  Scale est. = 6.6502    n = 455

• mgcv::gam does penalized regression. Set select = TRUE and then terms can be removed from the model during fitting. However, if your goal is prediction I'd recommend using independent training and testing datasets or at least cross-validation. – Roland Feb 23 '16 at 15:15
• Thanks Roland. I have added select , but I'm even more unsure of how to interpret this. It result in a model with fit statistics almost the same (slightly worse) in terms of r2, aic etc but the p-value for the s(x) term is now much lower. So if the parameter is not getting shrunk what is changing. – user2957945 Feb 23 '16 at 15:36
• @user2957945 can you edit your question to include the output from the model that used select = TRUE? – Gavin Simpson Feb 23 '16 at 19:01
• Hello @GavinSimpson , I have added the model output, thanks – user2957945 Feb 23 '16 at 20:49

One possible explanation is that the null hypothesis $s(x) = 0$ is false, but you have low power and so your p-value is not very impressive. Just because an effect is present doesn't mean it is easy to detect. That's why clinical trials must be designed with a certain effect size in mind (usually the MCID).
Another way to resolve this: different measures should give different results because they encode different priorities. AIC is a predictive criterion, and it behaves similarly to cross-validation. It may result in overly complex models that happen to have strong predictive performance. By contrast, mgcv's p-values are used to determine the presence or absence of a given effect$^*$, and predictive performance is a secondary concern.
$^*$substitute "association" for "effect" unless you're working with data from a controlled trial where $x$ was assigned randomly or you have other reasons to believe the observed association is causal.
• Thanks. Actually I was mistaken in my understanding. What I should of tested is the model with the s(x) term and then model with x (rather than no x. The significance of the smooth term (or otherwise) defaults to testing whether a smooth curve is needed compared to a linear (rather than omitting completely) – user2957945 Dec 9 '16 at 8:35