# ANOVA to compare models

I'm looking at this site for a workshop on GAM in R: http://qcbs.ca/wiki/r_workshop8

In the end of the section 2. Multiple smooth terms they show an example, where they use anova to compare three different models to determine the best fit model. The output is

  Analysis of Deviance Table
Model 1: y ~ x0 + s(x1)
Model 2: y ~ x0 + s(x1) + x2
Model 3: y ~ x0 + s(x1) + s(x2)
Resid. Df Resid. Dev      Df Deviance  Pr(>Chi)
1    394.08     5231.6
2    393.10     4051.3 0.97695   1180.2 < 2.2e-16 ***
3    385.73     1839.5 7.37288   2211.8 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Based on this they conclude that model 3 is best. My question is how they see that?

My current understanding is: The Pr(>Chi)-value is small for both model 2 and 3, so these are better than model 1. However, what other variable are they using to determine that 3 is better than 2?

• Hint: what variable(s) appear in Model 3 that do not appear in Model 2? To answer this, you need to know what "s(x2)" means--and that depends on how the function 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 from x2 to s(x2) adds $7.37288$ degrees of freedom--but that's all we can determine about it from this output. – whuber Apr 17 '17 at 14:53
• Doing AIC(model1, model2, model3) reveals that model 3 has a lower AIC. This could be a further proof that it is the optimal model among the three – BillyJean Apr 17 '17 at 15:59
• The huge difference in deviance between models (2) and (3) is compelling. – whuber Apr 17 '17 at 16:02

The output from anova() is a series of likelihood ratio tests. The lines in the output are:

1. 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.
2. 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.

3. 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

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).

• +1 for the detailed answer with code and reference. I will read and learn this answer in detail later. One question, you said it is likelihood ratio test. but in ?anova.lm there is no such option, can be F chisq or CP – Haitao Du Apr 19 '17 at 18:22
• @hxd1011 this is using mgcv:::anova.gam not the method for lm models. These are analysis of deviance tests but that's the same thing as likelihood ratios. – Reinstate Monica - G. Simpson Apr 19 '17 at 18:46
• Thanks. Could you answer my question here with some high level summary? Or your answer is here is already covered. – Haitao Du Apr 19 '17 at 18:49
• @hxd1011 You'd need to be more specific. What kinds of models? There are a lot of assumptions behind anova() but which ones depends on what the model is. Often for non-Gaussian models these are performing likelihood ratio tests or similar tests, but the assumptions will vary; they do vary even for GLMs and GAMs. anova() is a convenience function but it isn't doing ANOVA sensu the general linear model for anything but a general linear model (fitted via lm() say). – Reinstate Monica - G. Simpson Apr 19 '17 at 18:54
• @DeltaIV Sorry, that was poor wording on my part. What select = TRUE does is fully penalize all smooth terms, which AFAIU makes comparisons using REML OK. I haven't looked at the detail of the new AIC for GAMS to see what it would do with the extra penalties added when using select = TRUE. So if we drop select = TRUE to be on the safe side, we have the problem that REML is not a true likelihood and is not used in AIC comparisons because it is dependent upon the fixed effects in the model. Accounting for both concerns means I used method = "ML" (not method = "REML") when refitting. – Reinstate Monica - G. Simpson Apr 20 '17 at 17:10

You may want to test the two models with lrest.

lrtest(two_term_model, two_smooth_model)

Model 1: y ~ x0 + s(x1) + x2
Model 2: y ~ x0 + s(x1) + s(x2)
#Df  LogLik    Df  Chisq Pr(>Chisq)
1  8.1107 -995.22
2 15.0658 -848.95 6.955 292.55  < 2.2e-16 ***


While adding a smooth function to both terms indeed complicate the model, the improvement in the log-likelihood is significant. This shouldn't be surprising because the data was generated by a GAM simulator.

You may also want to print out the summary statistics:

Link function: identity

Formula:
y ~ x0 + s(x1) + x2

Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  11.6234     0.3950  29.429  < 2e-16 ***
x02           2.1147     0.4180   5.059 6.48e-07 ***
x03           4.3813     0.4172  10.501  < 2e-16 ***
x04           6.2644     0.4173  15.010  < 2e-16 ***
x2           -6.4110     0.5212 -12.300  < 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.111  2.626 64.92  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.583   Deviance explained = 58.9%
GCV = 8.7944  Scale est. = 8.6381    n = 400


and

Family: gaussian

Formula:
y ~ x0 + s(x1) + s(x2)

Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   8.3328     0.2074  40.185  < 2e-16 ***
x02           2.1057     0.2955   7.125 5.15e-12 ***
x03           4.3715     0.2934  14.901  < 2e-16 ***
x04           6.1197     0.2935  20.853  < 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.691  3.343 95.00  <2e-16 ***
s(x2) 7.375  8.356 85.07  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.796   Deviance explained = 80.2%
GCV = 4.3862  Scale est. = 4.232     n = 400


Note the difference in deviance explained (it's huge). The more complicated model also has better R-sq.(adj). The second smoothing term is highly significant and fits nicely with the data.

• Doesn't this just produce another example like that in the question? Could you indicate more explicitly how it responds to the question? – whuber Apr 18 '17 at 3:27