8
$\begingroup$

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 
# Link function: identity 
# 
# 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 
# Link function: identity 
# 
# 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 ***

EDIT: added model from comments

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

#Family: gaussian 
#Link function: identity 

#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
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Roland Feb 23 '16 at 15:15
  • $\begingroup$ 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. $\endgroup$ – user2957945 Feb 23 '16 at 15:36
  • $\begingroup$ @user2957945 can you edit your question to include the output from the model that used select = TRUE? $\endgroup$ – Gavin Simpson Feb 23 '16 at 19:01
  • $\begingroup$ Hello @GavinSimpson , I have added the model output, thanks $\endgroup$ – user2957945 Feb 23 '16 at 20:49
2
$\begingroup$

tl;dr: AIC is predictive whereas p-values are for inference. Also, your test of significance may simply lack power.

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.

$\endgroup$
  • $\begingroup$ 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) $\endgroup$ – user2957945 Dec 9 '16 at 8:35
  • $\begingroup$ Oooh, good to know! $\endgroup$ – eric_kernfeld Dec 9 '16 at 15:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.