The output from anova()
is a series of likelihood ratio tests. The lines in the output are:
- The first line in the output corresponds to the simplest model with only a smooth of
x1
(I'm ignoring the factor x0
as it isn't up for consideration in your example) — this is not tested against anything simpler hence the last few column entries are empty.
The second line is a likelihood ratio test between the model in line 1 and the model in line 2. At the cost of 0.97695
extra degrees of freedom, the residual deviance is decreased by 1180.2
. This reduction in deviance (or conversely, increase in deviance explained), at the cost of <1 degree of freedom, is highly unlikely if the true effect of x2
were 0.
Why 0.97695
degrees of freedom increase? Well, the linear function of x2
would add 1 df to the model but the smoother for x1
will be penalised back a little bit more than before and hence use slightly fewer effective degrees of freedom, hence the <1 change in overall degrees of freedom.
- The third line is exactly the same as I described above but for a comparison between the model in the second line and the model in the third line: i.e. the third line is evaluating the improvement in moving from modelling
x2
as a linear term to modelling x2
as a smooth function. Again, this improvement in model fit (change in deviance is now 2211.8
at the cost of 7.37288
more degrees of freedom) is unlikely if the extra parameters associated with s(x2)
were all equal to 0.
In summary, line 2 says Model 2 fits better than Model 1, so a linear function of x2
is better than no effect of x1
. But line 3 says that Model 3 fits the data better than Model 2, so a smooth function of x2
is preferred over a linear function of x2
. This is a sequential analysis of models, not a series of comparisons against the simplest model.
However…
What they're showing is not the best way to do this — recent theory would suggest that the output from summary(m3)
would have the most "correct" coverage properties. Furthermore, to select between models, one should probably use select = TRUE
when fitting the full model (the one with two smooths), which would allow for shrinkage of terms that would include the model with linear x2
or even no effect of this variable. They're also not fitting using REML or ML smoothness selection, which many of us mgcv users would consider the default option (even though it isn't the actual default in gam()
)
What I would do is:
library("mgcv")
gam_data <- gamSim(eg=5)
m3 <- gam(y ~ x0 + s(x1) + s(x2), data = gam_data, select = TRUE,
method = "REML")
summary(m3)
The final line produces the following:
> summary(m3)
Family: gaussian
Link function: identity
Formula:
y ~ x0 + s(x1) + s(x2)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.4097 0.2153 39.053 < 2e-16 ***
x02 1.9311 0.3073 6.284 8.93e-10 ***
x03 4.4241 0.3052 14.493 < 2e-16 ***
x04 5.7639 0.3042 18.948 < 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(x1) 2.487 9 25.85 <2e-16 ***
s(x2) 7.627 9 76.03 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.769 Deviance explained = 77.7%
-REML = 892.61 Scale est. = 4.5057 n = 400
We can see that both smooth terms are significantly different from null functions.
What select = TRUE
is doing is putting an extra penalty on the null space of the penalty (this is the part of the spline that is perfectly smooth). If you don't have this, smoothness selection can only penalise a smooth back to a linear function (because the penalty that's doing smoothness selection only works on the non-smooth (the wiggly) parts of the basis). To perform selection we need to be able to penalise the null space (the smooth parts of the basis) as well.
select = TRUE
achieves this through the use of a second penalty added to all smooth terms in the model (Marra and Wood, 2011). This acts as a kinds of shrinkage, pulling all smooth terms somewhat towards 0, but it will pull superfluous terms towards 0 much more quickly, hence selecting them out of the model if they don't have any explanatory power. We pay a price for this when evaluating the significance of the smooths; note the Ref.df
column above (the 9 comes from the default value of k = 10
, which for thin plate splines with centring constraints means 9 basis functions), instead of paying something like 2.5 and 7.7 degrees of freedom for the splines, we're paying 9 degrees of freedom each. This reflects that fact that we've done the selection, that we weren't sure which terms should be in the model.
Note: it is important that you don't use anova(m1, m2, m3)
type calls on models using select = TRUE
. As noted in ?mgcv:::anova.gam
, the approximation used can be very bad for smooths with penalties on their null spaces.
In the comments, @BillyJean mentioned using AIC for selection. Recent work by Simon Wood and colleagues (Wood et al, 2016) derived an AIC that accounts for the extra uncertainty due to us having estimated the smoothness parameters in the model. This AIC works reasonably well, but there is some discussion as to the behaviour of their derivation of AIC when IIRC smooths are close to linear functions. Anyway, AIC would give us:
m1 <- gam(y ~ x0 + s(x1), data = gam_data, method = "ML")
m2 <- gam(y ~ x0 + s(x1) + x2, data = gam_data, method = "ML")
m3 <- gam(y ~ x0 + s(x1) + s(x2), data = gam_data, method = "ML")
AIC(m1, m2, m3)
> AIC(m1, m2, m3)
df AIC
m1 7.307712 2149.046
m2 8.608444 2055.651
m3 16.589330 1756.890
Note I refitted all of these with ML smoothness selection as I'm not certain what the AIC does when select = TRUE
and you have to be careful comparing models with different fixed effects, that aren't fully penalised, using REML.
Again the inference is clear; the model with smooths of x1
and x2
has substantially better fit than either of the other two models.
Marra, G. & Wood, S. N. Practical variable selection for generalized additive models. Comput. Stat. Data Anal. 55, 2372–2387 (2011).
Wood, S. N., Pya, N. & Säfken, B. Smoothing Parameter and Model Selection for General Smooth Models. J. Am. Stat. Assoc. 111, 1548–1563 (2016).
s
is defined. (I presume it's some kind of spline, but I am reluctant to presume any more than that.) We can tell from the output that it's fairly complicated--the move fromx2
tos(x2)
adds $7.37288$ degrees of freedom--but that's all we can determine about it from this output. $\endgroup$AIC(model1, model2, model3)
reveals that model 3 has a lowerAIC
. This could be a further proof that it is the optimal model among the three $\endgroup$