Exercise 14.1

library(mgcv)
## Loading required package: nlme
## This is mgcv 1.9-4. For overview type '?mgcv'.
hr <- read.table("https://www.richardpmann.com/MATH3701/Datasets/SAheart.txt",
                  sep=",", head=TRUE, row.names=1)

gam.full <- gam(chd ~ s(tobacco, k=10) + s(alcohol, k=10) + s(obesity, k=10) +
                   s(sbp, k=10) + s(ldl, k=10) + s(adiposity, k=10) + s(typea, k=10) +
                   famhist,
                 family="binomial", data=hr)
summary(gam.full)
## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## chd ~ s(tobacco, k = 10) + s(alcohol, k = 10) + s(obesity, k = 10) + 
##     s(sbp, k = 10) + s(ldl, k = 10) + s(adiposity, k = 10) + 
##     s(typea, k = 10) + famhist
## 
## Parametric coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -1.2703     0.1668  -7.618 2.58e-14 ***
## famhistPresent   1.0370     0.2306   4.497 6.89e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                edf Ref.df Chi.sq  p-value    
## s(tobacco)   6.145  7.246 28.049 0.000284 ***
## s(alcohol)   1.000  1.000  0.232 0.629933    
## s(obesity)   1.905  2.444  9.346 0.014852 *  
## s(sbp)       1.445  1.771  5.270 0.102268    
## s(ldl)       1.000  1.000  8.855 0.002924 ** 
## s(adiposity) 2.729  3.448 11.334 0.013465 *  
## s(typea)     2.328  2.982 10.623 0.015104 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.251   Deviance explained = 23.4%
## UBRE = 0.06928  Scale est. = 1         n = 462

(a)

The edf for each smooth term (from the summary table) are:

Term edf
tobacco 6.14
alcohol 1.00
obesity 1.90
sbp 1.45
ldl 1.00
adiposity 2.73
typea 2.33

An edf close to 1 means the penalty has shrunk that smooth term all the way to a straight line, i.e. no evidence of non-linearity is being used. By that standard: alcohol and ldl (edf \(=1.00\) exactly) show no non-linearity at all; sbp (edf \(=1.45\)) shows only a little; obesity (edf \(=1.90\)) and typea (edf \(=2.33\)) show moderate non-linearity; and tobacco (edf \(=6.14\)) shows the strongest non-linearity by a clear margin, and adiposity (edf \(=2.73\)) also shows a fair amount.

(b)

gam.noalc <- gam(chd ~ s(tobacco, k=10) + s(obesity, k=10) + s(sbp, k=10) +
                    s(ldl, k=10) + s(adiposity, k=10) + s(typea, k=10) + famhist,
                  family="binomial", data=hr)

D.full <- deviance(gam.full); D.noalc <- deviance(gam.noalc)
edf.full <- sum(gam.full$edf); edf.noalc <- sum(gam.noalc$edf)

dev.diff <- D.noalc - D.full
df.diff  <- edf.full - edf.noalc
dev.diff; df.diff
## [1] 0.2602646
## [1] 0.999577
pchisq(dev.diff, df=df.diff, lower.tail=FALSE)
## [1] 0.6097624

Removing alcohol increases the deviance only slightly, from \(456.91\) to \(457.17\) (a change of \(0.260\)), while the total edf falls from \(18.55\) to \(17.55\) – a change of almost exactly \(1\) degree of freedom (unsurprising, since alcohol’s own edf in the full model was \(1.00\), i.e. it was already using no more flexibility than a single linear term). Referring \(D_1-D_2=0.260\) to \(\chi^2_{1}\) gives \(p=0.610\), so there is no evidence that alcohol improves the fit – consistent with the very high \(p\)-value for s(alcohol) in the full-model summary table. We would drop alcohol from the model.

(c)

gam.full$gcv.ubre
##     GCV.Cp 
## 0.06928044
gam.noalc$gcv.ubre
##     GCV.Cp 
## 0.06551661

The GCV score is \(0.06928\) for the full model and \(0.06552\) for the model without alcohol. This suggests that removing alcohol actually improves predictions. This agrees with the conclusion from the deviance test in (b): alcohol was not contributing any real explanatory power. Both routes point to the same answer: alcohol can be dropped from the model.

(d)

gam.notob <- gam(chd ~ s(alcohol, k=10) + s(obesity, k=10) + s(sbp, k=10) +
                    s(ldl, k=10) + s(adiposity, k=10) + s(typea, k=10) + famhist,
                  family="binomial", data=hr)

D.notob <- deviance(gam.notob)
edf.notob <- sum(gam.notob$edf)

dev.diff2 <- D.notob - D.full
df.diff2  <- edf.full - edf.notob
dev.diff2; df.diff2
## [1] 38.65358
## [1] 6.730893
pchisq(dev.diff2, df=df.diff2, lower.tail=FALSE)
## [1] 1.753606e-06
gam.full$gcv.ubre
##     GCV.Cp 
## 0.06928044
gam.notob$gcv.ubre
##    GCV.Cp 
## 0.1238081

This time the picture is completely different. Removing tobacco increases the deviance from \(456.91\) to \(495.56\), a change of \(38.65\), against a change in edf of about \(6.73\) (close to, though not identical to, tobacco‘s own edf of \(6.14\) in the full model, because refitting without tobacco lets the other smooth terms’ penalties shift slightly too, so the edf change isn’t exactly equal to the removed term’s own edf). Referring \(38.65\) to \(\chi^2_{6.73}\) gives \(p\approx1.8\times10^{-6}\): overwhelming evidence that tobacco matters. The GCV score confirms this, nearly doubling from \(0.0693\) to \(0.1238\) once tobacco is removed. This shows that removing tobacco makes predictions much worse.

(e)

The two predictors behave in opposite ways, and the deviance test and GCV comparison agree with each other in both cases: for alcohol, a tiny, non-significant change in deviance (\(p=0.61\)) is matched by an improved GCV score, telling us alcohol does not improve predictions. For tobacco a highly significant change in deviance (\(p\approx1.8\times10^{-6}\)) is matched by a much worse GCV score, telling us tobacco is one of the most important predictors in the model. This is a reassuring result: two quite different ways of asking “does this term matter?” – a formal hypothesis test based on the sampling distribution of the deviance, and an out-of-sample-flavoured criterion based on prediction accuracy adjusted for complexity – give the same answer for both terms.

Exercise 14.2

First fit the model again to get the model summary:

library(mgcv)
gam2 = gam(chd ~ age + s(tobacco, k=20) + famhist,
           family = 'binomial', data=hr)
summary.gam(gam2)
## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## chd ~ age + s(tobacco, k = 20) + famhist
## 
## Parametric coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -3.273522   0.471089  -6.949 3.68e-12 ***
## age             0.047544   0.009684   4.910 9.12e-07 ***
## famhistPresent  0.962834   0.223318   4.311 1.62e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##             edf Ref.df Chi.sq p-value  
## s(tobacco) 6.08  7.573  17.89  0.0178 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.212   Deviance explained = 19.1%
## UBRE = 0.083267  Scale est. = 1         n = 462

(a)

We know from the model summary that if FamilyHist is Absent: \[\log(p/(1-p)) = -3.273522 + 0.047544 \times \text{Age} + s(\text{Tobacco})\] where \(s(\cdot)\) is the smooth function for tobacco consumption.

We are told \(p=0.20\) when \(\text{Age}=22\) and \(\text{Tobacco} = 10\). So we can calculate: \[\log(0.2/0.8) = -\log(4) = -3.273522 + 0.047544 \times 22 + s(10)\] Rearranging and simplifying gives: \[s(10) = 3.273522-\log(4) - 0.047544 \times 22 = 0.8412596\]

To calculate the desired probability we substitute this value for \(s(10\) back in. For Age = 55, we want: \[p = \frac{\exp(-3.273522 + 0.047544 \times 55 + s(10))}{1 + \exp(-3.273522 + 0.047544 \times 55 + s(10))} = \frac{\exp(-3.273522 + 0.047544 \times 55 + 0.8412596}{1 + \exp(-3.273522 + 0.047544 \times 55 + 0.8412596} = 0.55\]

(b) Now we need to add in the effect from FamilyHist=Present. Here, the model summary tells us: \[\log(p/(1-p)) = -3.273522 + 0.962834 + 0.047544 \times \text{Age} + s(\text{Tobacco})\] The tobacco consumption remains 10kg, we can substitute in the value of \(s(10) = 0.8412596\) that we calculated before. We want to find the value of Age that gives \(p=0.5\), so we have: \[\log(0.5/0.5) = 0 = -3.273522 + 0.962834 + 0.047544 \times \text{Age} + 0.8412596\] Rearranging and simplifying gives: \[\text{Age} = \frac{3.273522 - 0.962834 -0.8412596}{0.047544} = 30.9\]

Exercise 14.3

library(mgcv); library(MASS)
data(mcycle)

out <- gam(accel ~ s(times, k=20), data=mcycle)
summary(out)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## accel ~ s(times, k = 20)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -25.546      1.966     -13   <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(times) 10.9  13.22 34.9  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =   0.78   Deviance explained = 79.8%
## GCV = 564.33  Scale est. = 513.84    n = 133
out$sp
##    s(times) 
## 0.002207576

(a)

The GCV-selected smoothing parameter is sp \(\approx 0.0022\), and the smooth term for times has edf \(\approx10.9\) (giving a total model edf, including the intercept, of about \(11.9\)). This is a fairly flexible fit – well above the \(\nu=2\) “essentially linear” end of the scale – reflecting the sharp dip-and-rebound pattern visible in the data around the crash.

(b)

newdat <- data.frame(times = seq(min(mcycle$times), max(mcycle$times), length=200))
pred <- predict(out, newdata=newdat, se.fit=TRUE)

plot(mcycle$times, mcycle$accel, pch=16, col="grey40",
     xlab="time (ms)", ylab="acceleration (g)",
     main="GAM fit with 95% confidence band")
lines(newdat$times, pred$fit, col="blue", lwd=2)
lines(newdat$times, pred$fit + 1.96*pred$se.fit, col="blue", lty=2)
lines(newdat$times, pred$fit - 1.96*pred$se.fit, col="blue", lty=2)

The fitted curve traces the characteristic pattern of the crash test clearly: near-zero acceleration before impact, a sharp dip to around \(-110g\) at roughly 20ms, a rebound to a smaller positive peak around 30-35ms, before settling back towards zero with some residual oscillation. The confidence band is narrow where the data are dense (e.g. before 15ms) and widens noticeably in the sparser regions and near the extremes of the time range.

(c)

out.under <- gam(accel ~ s(times, k=20, sp=1e-4), data=mcycle)
out.over  <- gam(accel ~ s(times, k=20, sp=1e2), data=mcycle)
sum(out.under$edf); sum(out.over$edf)
## [1] 18.40118
## [1] 2.048793
par(mfrow=c(1,3))
for (fit in list(out.under, out, out.over)) {
  pr <- predict(fit, newdata=newdat)
  plot(mcycle$times, mcycle$accel, pch=16, col="grey60", xlab="time", ylab="accel")
  lines(newdat$times, pr, col="red", lwd=2)
}

With sp=\(10^{-4}\) the fit is noticeably more flexible (edf \(\approx18.4\)) and starts to chase individual noisy observations, particularly visible as small spurious wiggles superimposed on the main dip-and-rebound shape. With sp=\(10^2\) the fit is heavily oversmoothed (edf \(\approx2.05\), essentially a straight line) and completely misses the crash dip and rebound – it captures only the (slight) overall downward trend across the full time range and none of the interesting structure. The GCV-optimal fit from (a)/(b) sits sensibly between these two extremes: flexible enough to capture the real dip-and-rebound shape, but without the extraneous wiggle seen at sp=\(10^{-4}\).


End of Solutions to Exercises 14