Exercise 13.1

bmd <- read.csv("https://www.richardpmann.com/MATH3701/Datasets/bmd.csv")

(a)

plot(bmd$age, bmd$Tscore, pch=16,
     xlab="age", ylab="Tscore",
     main="Bone mineral density T-score vs age")

The relationship is clearly decreasing overall (older women tend to have lower T-scores, i.e. lower bone density), but it does not look linear. Through the 40s the T-scores sit close to 0 with very little scatter; from about age 50 onwards the scores decline more steeply and the scatter increases substantially. This suggests the rate of decline may accelerate with age, though the increasing noise at older ages makes this hard to judge from the scatter plot alone – a smoothing spline should help clarify this.

(b)

fit.small <- smooth.spline(bmd$age, bmd$Tscore, lambda=1e-6)
fit.large <- smooth.spline(bmd$age, bmd$Tscore, lambda=10)
fit.small$lambda; fit.small$df
## [1] 1e-06
## [1] 25.18682
fit.large$lambda; fit.large$df
## [1] 10
## [1] 2.013645
plot(bmd$age, bmd$Tscore, pch=16, xlab="age", ylab="Tscore",
     main="lambda=1e-6 vs lambda=10")
lines(fit.small, col="red", lwd=2)
lines(fit.large, col="blue", lwd=2)
legend("topright", legend=c("lambda=1e-6","lambda=10"), col=c("red","blue"), lwd=2)

With \(\lambda=10^{-6}\) we get df \(=25.2\): this is close to interpolating the data (note there are only 29 distinct ages among the 66 women, which is why the maximum achievable df here is around 29 rather than 66) and the fitted curve is very wiggly, chasing individual points and small clusters of points rather than showing a clear trend – classic overfitting.

With \(\lambda=10\) we get df \(=2.01\), essentially indistinguishable from the straight-line OLS solution (\(\nu=2\)). This fit is too rigid: it forces a single constant rate of decline across the whole age range and misses the possibility that decline accelerates with age that we noted in part (a) – underfitting.

(c)

fit.gcv <- smooth.spline(bmd$age, bmd$Tscore)
fit.gcv$lambda; fit.gcv$df
## [1] 2.522923e-07
## [1] 27.70183
plot(bmd$age, bmd$Tscore, pch=16, xlab="age", ylab="Tscore",
     main="GCV-selected smoothing spline")
lines(fit.gcv, col="darkgreen", lwd=2)

GCV selects \(\lambda \approx 2.5\times10^{-7}\), giving df \(=27.7\). Does this look reasonable? Only partially. The curve does capture the overall accelerating decline suggested in part (a), but it is still quite wiggly – almost as wiggly as the \(\lambda=10^{-6}\) fit in part (b) – and zig-zags noticeably between individual points and small clusters in the 55-65 age range rather than showing a clean trend. With only \(n=66\) observations (and only 29 distinct ages), automatic GCV selection can be unstable and tends to slightly overfit here; a human eye might reasonably choose a somewhat larger \(\lambda\) than GCV does.

(d)

lambdas.eye <- 10^seq(-6, -1, by=1)
plot(bmd$age, bmd$Tscore, pch=16, col="grey50", xlab="age", ylab="Tscore",
     main="Exploring lambda by eye")
cols <- c("red","orange","darkgreen","blue","purple","brown")
for (i in seq_along(lambdas.eye)) {
  f <- smooth.spline(bmd$age, bmd$Tscore, lambda=lambdas.eye[i])
  lines(f, col=cols[i], lwd=2)
}
legend("topright", legend=paste0("lambda=",formatC(lambdas.eye, format="e", digits=0)),
       col=cols, lwd=2)

fit.eye <- smooth.spline(bmd$age, bmd$Tscore, lambda=1e-3)
fit.eye$df; fit.eye$lambda
## [1] 6.631113
## [1] 0.001
r <- residuals(fit.eye)
plot(bmd$age, r, pch=16, col="grey30", xlab="age", ylab="residual",
     main="Residuals from by-eye fit (lambda=1e-3, df=6.6)")
abline(h=0, lty=2, col="red")

By eye, \(\lambda=10^{-6}\) is clearly too wiggly – it bends to chase individual points and small local clusters, just like the fit in (b) – while \(\lambda=10^{-1}\) looks nearly straight again and starts to lose the visible acceleration in the decline after age 55. Something around \(\lambda=10^{-3}\) (df\(\approx6.6\)) looks like a reasonable compromise: smooth enough not to chase noise, but still clearly curving rather than a straight line, tracking the flat T-scores through the 40s and the steeper decline from the mid-50s onward that we saw in (a).

The residual plot, though, reveals something the fitted curve alone did not: the residuals are small and tightly clustered around zero for ages under about 50 (sd \(\approx0.11\)), but noticeably more spread out from age 55 onward (sd \(\approx0.64\), roughly six times larger), including several residuals beyond \(\pm1\). This is a classic fan/funnel pattern – not a sign that the mean curve is wrong, but a sign that the variance of Tscore around that curve is not constant with age, which is an assumption smoothing splines (and GCV, which is built on a constant-variance argument) take for granted. In other words, the by-eye fit looks visually sensible, but the residual plot shows the data are more heteroscedastic than a single smoothing-spline fit accounts for.

(e)

Recall \(\nu=2\) (cubic spline) and, strictly, \(n=66\), though as noted above the practical upper limit on df for this particular dataset is closer to 29 because of the repeated ages. GCV’s edf of 27.7 sits almost at that practical upper limit, and a long way from \(\nu=2\) – in the language of @sec-edf, GCV has chosen a fit that is close to interpolation rather than close to the OLS line. This supports the visual impression from (a) that the age-density relationship is genuinely non-linear (a straight line would clearly underfit the steepening decline at older ages), but the fact that GCV’s chosen edf is so close to its achievable maximum is also a sign that some of that flexibility may be fitting sampling noise rather than a real population trend, especially given the small, unevenly-spaced sample – and the heteroscedasticity uncovered in (d) is a plausible reason why: GCV may be inflating the fitted edf partly to chase the extra noise variance at older ages, rather than genuine extra curvature. A slightly smoother fit than the literal GCV minimiser, closer to the by-eye choice in (d), might be a more defensible summary of the age effect.

Exercise 13.2

iron <- read.csv("https://www.richardpmann.com/MATH3701/Datasets/iron.csv")

(a)

plot(iron$Year, iron$iron, pch=16, type="b",
     xlab="Year", ylab="iron", main="Iron measurement vs Year")

There is a clear, strong downward trend (correlation \(-0.88\)) and, by eye, a straight line looks like a fairly good description of it – there is no obvious curvature such as a plateau or a change of direction, just a roughly constant year-on-year decline with a fair amount of noise superimposed. Whether a smoothing spline is actually needed (as opposed to just permitted) is exactly what parts (b)-(d) investigate.

(b)

test <- iron[iron$Year >= 1990, ]
train <- iron[iron$Year < 1990, ]

lambdas <- 10^seq(-4, 4, by=0.5)
Q <- numeric(length(lambdas))
for (i in seq_along(lambdas)) {
  fit <- smooth.spline(train$Year, train$iron, lambda=lambdas[i])
  pred <- predict(fit, test$Year)$y
  Q[i] <- sum((test$iron - pred)^2)
}

plot(lambdas, Q, log="x", type="b",
     xlab=expression(lambda), ylab=expression(Q[I[1]:I[2]](lambda)),
     main="Train/test criterion vs lambda")

(c)

best.lambda <- lambdas[which.min(Q)]
best.lambda; min(Q)
## [1] 1e-04
## [1] 105.567
fit.best <- smooth.spline(iron$Year, iron$iron, lambda=best.lambda)
fit.best$df
## [1] 10.51614
plot(iron$Year, iron$iron, pch=16, xlab="Year", ylab="iron",
     main="Full-data fit at train/test-selected lambda")
lines(fit.best, col="purple", lwd=2)

Over the specified grid, \(Q_{I_1:I_2}(\lambda)\) is minimised at \(\lambda = 10^{-4} \approx 0.0001\) (with \(Q\approx 105.6\)), which is the smallest value in the grid. Refitting the full dataset at this \(\lambda\) gives df \(\approx10.5\) – a fairly flexible fit that follows the medium-term wiggles in the series (e.g. the rise around 1968 and the dip around 1990) rather than reducing to a straight line.

(d)

fit.gcv2 <- smooth.spline(iron$Year, iron$iron)
fit.gcv2$lambda; fit.gcv2$df
## [1] 424.4194
## [1] 2.000297

GCV on the full dataset chooses \(\lambda \approx 424\), giving df \(\approx 2.00\) – essentially indistinguishable from the OLS straight line (\(\nu=2\)). This is a striking disagreement with the train/test choice in (c), which preferred a much more flexible fit (df \(\approx10.5\)).

Which is more reliable here? I would trust the GCV answer more in this case, for two reasons. First, GCV uses the whole 50-year series and an averaged leave-one-out-type criterion, whereas the train/test approach depends on a single hold-out block of only 10 points at the very end of the series – a high-variance estimate that could easily be swayed by the particular pattern of those 10 years (for example the unusually low value in 1990 right at the start of the test window). Second, the train/test criterion is really measuring extrapolation accuracy just past the edge of the training data, and a flexible spline’s behaviour right at the boundary of its training range can be quite different from its behaviour on new data further out – so a \(\lambda\) that happens to extrapolate well into exactly this one test block need not reflect the smoothness of the underlying trend. The GCV answer is also more consistent with the visual impression from part (a), where a straight line already looked like a reasonable description of the data.


End of Solutions to Exercises 13