# ANOVA for natural splines

I want to use ANOVA for selecting the best model with two independent variables, where either one of them (or both) can be fitted with natural splines. Hence, I try the following models:

library("MASS") # Includes the Boston dataset
gam.1 <- gam(medv ~ lstat + rm, data=Boston)
gam.2 <- gam(medv ~ ns(lstat, df=2) + rm, data=Boston)
gam.3 <- gam(medv ~ ns(lstat, df=3) + rm, data=Boston)
gam.4 <- gam(medv ~ ns(lstat, df=4) + rm, data=Boston)
gam.5 <- gam(medv ~ ns(lstat, df=5) + rm, data=Boston)
gam.6 <- gam(medv ~ lstat + ns(rm, df=2), data=Boston)
gam.7 <- gam(medv ~ lstat + ns(rm, df=3), data=Boston)
gam.8 <- gam(medv ~ lstat + ns(rm, df=4), data=Boston)
gam.9 <- gam(medv ~ lstat + ns(rm, df=5), data=Boston)
gam.10 <- gam(medv ~ ns(lstat, df=2) + ns(rm, df=2), data=Boston)
gam.11 <- gam(medv ~ ns(lstat, df=3) + ns(rm, df=3), data=Boston)
gam.12 <- gam(medv ~ ns(lstat, df=4) + ns(rm, df=4), data=Boston)
gam.13 <- gam(medv ~ ns(lstat, df=5) + ns(rm, df=5), data=Boston)


I am new to ANOVA, so I am not sure, how to interpret the following table (i.e. which model would be the best):

   Resid. Df Resid. Dev Df Deviance        F    Pr(>F)
1        503    15439.3
2        502    12801.2  1  2638.11 141.0841 < 2.2e-16 ***
3        501    12555.7  1   245.51  13.1294 0.0003207 ***
4        500    11642.0  1   913.69  48.8634 8.918e-12 ***
5        499    11543.8  1    98.16   5.2494 0.0223741 *
6        502    11034.6 -3   509.22
7        501    11038.5  1    -3.84
8        500    10504.8  1   533.69  28.5410 1.402e-07 ***
9        499    10252.3  1   252.46  13.5011 0.0002644 ***
10       501    10425.3 -2  -173.02   4.6264 0.0102174 *
11       499    10241.5  2   183.82   4.9153 0.0076951 **
12       497     9569.0  2   672.48  17.9819 2.892e-08 ***
13       495     9255.9  2   313.08   8.3717 0.0002657 ***


Actually, I am not sure if it makes sense at all to use ANOVA on all these models together, or if I should break them into 3 groups (i.e. 1. natural spline on lstat, 2. natural spline on rm, 3. natural spline on both) and use ANOVA separately for each group.

What particularly confuses me on this result are the missing values in F & Pr(>F) columns for some rows (except the first one).

So: should I split these models in separate groups for using ANOVA? And if I do, how would I choose the best one among the groups? By performing ANOVA again or by, for example, just checking train/test score?

Assuming this is the gam package, I would suggest that you rethink this and fit GAMs using penalised splines, where the selection of the complexity of the smooths becomes part of the model fitting. This sort of model is easily fitted using the mgcv package, which ships with R.

library('mgcv')
data('Boston', package = 'MASS')

m <- gam(medv ~ s(lstat) + s(rm), data = Boston, method = 'REML')


This results in

r$> summary(m) Family: gaussian Link function: identity Formula: medv ~ s(lstat) + s(rm) Parametric coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 22.5328 0.1912 117.9 <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(lstat) 5.458 6.634 68.28 <2e-16 *** s(rm) 6.731 7.863 32.47 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 R-sq.(adj) = 0.781 Deviance explained = 78.7% -REML = 1470.5 Scale est. = 18.491 n = 506  The output strongly suggests that limiting the smooths to only 5 degrees of freedom was an underestimate of the complexity. If you really want to specifically test for extra wiggliness over and above a linear effect you can do that with some appropriate manipulation of the bases used to fit the GAM: m <- gam(medv ~ lstat + s(lstat, bs = 'tp', m = c(2,0)) + rm + s(rm, bs = 'tp', m = c(2,0)), data = Boston, method = 'REML') summary(m)  which gives: r$> summary(m)

Family: gaussian

Formula:
medv ~ lstat + s(lstat, bs = "tp", m = c(2, 0)) + rm + s(rm,
bs = "tp", m = c(2, 0))

Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  63.0186    21.9730   2.868  0.00431 **
lstat        -1.0374     0.3184  -3.258  0.00120 **
rm           -4.3535     3.4925  -1.247  0.21317
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
edf Ref.df     F  p-value
s(lstat) 4.458      8  4.26 1.12e-07 ***
s(rm)    5.731      8 18.38  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.781   Deviance explained = 78.7%
-REML = 1472.1  Scale est. = 18.491    n = 506


Again, there is clear support for a non-linear effect as the smooth terms are both unlikely under the null hypothesis tested here. This model works for thin plate splines (TPRS: bs = 'tp') because we can specify that we want a second-order wiggliness penalty but no null space (m = c(2, 0)). No null space means there are no perfectly smooth functions in the span of the basis; for the TPRS this means that the linear function has been removed from the basis, so the resulting basis only contains non-linear functions.

If you want to potentially exclude one of the two terms from the model then we need a slightly different approach; we need to penalise the null space as well as the range space (the wiggly bits). In the first model we fitted, the simplest model you could end up with would be linear effects of both covariates; we couldn't get any simpler a model because the wiggliness penalty that is used to select how wiggly (complex) the fitted functions are doesn't affect the linear part of the basis as it is already perfectly smooth from the point of view of the penalty (its second derivative is 0).

If we want to admit the possibility that one or both covariates might have no effect on the response, then we need to be able to penalise functions in the null space of the penalty, those that are perfectly smooth. We have two ways in mgcv to achieve this:

1. by per-smooth using the bs = 'ts' basis, for shrinkage TPRS (or the bs = 'cs' basis, for cubic regression splines with shrinkage) or
2. for all smooths in the model using the select = TRUE argument to gam().

Option 1 works by tweaking the usual wiggliness penalty so that the smooth bits of the basis are no longer exactly smooth from the view point of the penalty. Option 2 works by adding a second penalty that only works on the null space to each smooth in the model.

Option 1 assumes that you want to shrink the wiggly parts of the basis more than shrink the linear bits. Option 2 assumes you a prior want to shrink both the wiggly (range space) and smooth (null space) parts of the basis the same. Option 2 tends to work slightly better than Option 1, but it does require an additional smoothness parameter to be estimated for every smooth in the model, whereas Option 1 doesn't involve any additional smoothness parameters being estimated.

### Option 1:

m <- gam(medv ~ s(lstat, bs = 'ts') + s(rm, bs = 'ts'),
data = Boston, method = 'REML')
summary(m)


giving:

r$> summary(m) Family: gaussian Link function: identity Formula: medv ~ s(lstat, bs = "ts") + s(rm, bs = "ts") Parametric coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 22.5328 0.1914 117.7 <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(lstat) 4.903 9 50.15 <2e-16 *** s(rm) 6.522 9 29.09 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 R-sq.(adj) = 0.781 Deviance explained = 78.6% -REML = 1480 Scale est. = 18.543 n = 506  ### Option 2 m <- gam(medv ~ s(lstat) + s(rm), data = Boston, method = 'REML', select = TRUE) summary(m)  giving r$> summary(m)

Family: gaussian

Formula:
medv ~ s(lstat) + s(rm)

Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  22.5328     0.1917   117.5   <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(lstat) 5.345      9 51.17  <2e-16 ***
s(rm)    5.744      9 28.25  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =   0.78   Deviance explained = 78.5%
-REML = 1476.4  Scale est. = 18.598    n = 506
`

It doesn't make much difference which you use in this example.