1 Exercise 6 — The Wage data set

1.1 (a) Polynomial regression: wage ~ age

We use 10-fold cross-validation to select the optimal polynomial degree d, then compare the result to nested-model hypothesis testing via anova().

data(Wage)

set.seed(1)
max.degree <- 10
cv.errors  <- rep(NA, max.degree)

for (d in 1:max.degree) {
  glm.fit      <- glm(wage ~ poly(age, d), data = Wage)
  cv.errors[d] <- cv.glm(Wage, glm.fit, K = 10)$delta[1]
}

plot(1:max.degree, cv.errors, type = "b", pch = 19,
     xlab = "Polynomial Degree", ylab = "10-fold CV Error",
     main = "CV Error vs. Polynomial Degree (wage ~ age)")

best.degree.cv <- which.min(cv.errors)
best.degree.cv
## [1] 9

The CV curve typically drops sharply from degree 1 to 2, flattens out, and the minimum is often achieved somewhere between degree 3 and degree 9, but the improvement past degree ~4 is usually negligible (“one-standard-error” style reasoning would pick a lower degree).

ANOVA / nested F-test comparison:

fit.1 <- lm(wage ~ age, data = Wage)
fit.2 <- lm(wage ~ poly(age, 2), data = Wage)
fit.3 <- lm(wage ~ poly(age, 3), data = Wage)
fit.4 <- lm(wage ~ poly(age, 4), data = Wage)
fit.5 <- lm(wage ~ poly(age, 5), data = Wage)
fit.6 <- lm(wage ~ poly(age, 6), data = Wage)

anova(fit.1, fit.2, fit.3, fit.4, fit.5, fit.6)
## Analysis of Variance Table
## 
## Model 1: wage ~ age
## Model 2: wage ~ poly(age, 2)
## Model 3: wage ~ poly(age, 3)
## Model 4: wage ~ poly(age, 4)
## Model 5: wage ~ poly(age, 5)
## Model 6: wage ~ poly(age, 6)
##   Res.Df     RSS Df Sum of Sq        F    Pr(>F)    
## 1   2998 5022216                                    
## 2   2997 4793430  1    228786 143.6636 < 2.2e-16 ***
## 3   2996 4777674  1     15756   9.8936  0.001675 ** 
## 4   2995 4771604  1      6070   3.8117  0.050989 .  
## 5   2994 4770322  1      1283   0.8054  0.369565    
## 6   2993 4766389  1      3932   2.4692  0.116201    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The ANOVA sequence usually shows that going from degree 1 → 2 and 2 → 3 is highly significant (very small p-values), degree 3 → 4 is borderline significant, and degree 4 → 5 (and beyond) is not significant. This points to degree 3 or 4 as a reasonable choice — broadly consistent with the CV result, which often selects a similar or slightly higher degree, though CV can sometimes pick a higher degree because it optimizes prediction error rather than testing nested significance.

Plot of the fitted polynomial (using the CV-selected degree):

agelims  <- range(Wage$age)
age.grid <- seq(agelims[1], agelims[2])

fit.best <- lm(wage ~ poly(age, best.degree.cv), data = Wage)
preds    <- predict(fit.best, newdata = list(age = age.grid), se = TRUE)
se.bands <- cbind(preds$fit + 2 * preds$se.fit,
                   preds$fit - 2 * preds$se.fit)

plot(Wage$age, Wage$wage, xlim = agelims, cex = 0.5, col = "darkgrey",
     xlab = "Age", ylab = "Wage",
     main = paste("Polynomial Fit, Degree =", best.degree.cv))
lines(age.grid, preds$fit, lwd = 2, col = "blue")
matlines(age.grid, se.bands, lwd = 1, col = "blue", lty = 3)

1.2 (b) Step function: wage ~ age

We use 10-fold cross-validation to select the optimal number of cuts (i.e., the number of levels created by cut(age, k)).

set.seed(1)
max.cuts       <- 10
cv.errors.step <- rep(NA, max.cuts)


for (k in 2:max.cuts) {
  Wage$age.cut  <- cut(Wage$age, k)
  fit.step      <- glm(wage ~ age.cut, data = Wage)
  cv.errors.step[k] <- cv.glm(Wage, fit.step, K = 10)$delta[1]
}

plot(2:max.cuts, cv.errors.step[2:max.cuts], type = "b", pch = 19,
     xlab = "Number of Cuts", ylab = "10-fold CV Error",
     main = "CV Error vs. Number of Cuts (Step Function)")

best.cuts <- which.min(cv.errors.step)
best.cuts
## [1] 8

CV error for step functions tends to decrease quickly and then level off; somewhere around 7–8 cuts is a typical optimum, after which additional cuts add little predictive value.

Plot of the resulting step-function fit:

fit.step.best <- glm(wage ~ cut(age, best.cuts), data = Wage)
preds.step    <- predict(fit.step.best, newdata = list(age = age.grid), se = TRUE)
se.bands.step <- cbind(preds.step$fit + 2 * preds.step$se.fit,
                        preds.step$fit - 2 * preds.step$se.fit)

plot(Wage$age, Wage$wage, xlim = agelims, cex = 0.5, col = "darkgrey",
     xlab = "Age", ylab = "Wage",
     main = paste("Step Function Fit,", best.cuts, "Cuts"))
lines(age.grid, preds.step$fit, lwd = 2, col = "red")
matlines(age.grid, se.bands.step, lwd = 1, col = "red", lty = 3)

2 Exercise 10 — The College data set

2.1 (a) Train/test split + forward stepwise selection

We predict Outstate (out-of-state tuition) using all other variables as candidate predictors, and perform forward stepwise selection on the training set.

data(College)

set.seed(1)
train <- sample(1:nrow(College), nrow(College) / 2)
test  <- setdiff(1:nrow(College), train)

College.train <- College[train, ]
College.test  <- College[test, ]

regfit.fwd  <- regsubsets(Outstate ~ ., data = College.train,
                           nvmax = 17, method = "forward")
reg.summary <- summary(regfit.fwd)
par(mfrow = c(1, 3))
plot(reg.summary$cp, type = "b", pch = 19,
     xlab = "Number of Variables", ylab = "Cp")
points(which.min(reg.summary$cp), min(reg.summary$cp), col = "red", pch = 19, cex = 1.5)

plot(reg.summary$bic, type = "b", pch = 19,
     xlab = "Number of Variables", ylab = "BIC")
points(which.min(reg.summary$bic), min(reg.summary$bic), col = "red", pch = 19, cex = 1.5)

plot(reg.summary$adjr2, type = "b", pch = 19,
     xlab = "Number of Variables", ylab = "Adjusted R2")
points(which.max(reg.summary$adjr2), max(reg.summary$adjr2), col = "red", pch = 19, cex = 1.5)

par(mfrow = c(1, 1))
which.min(reg.summary$cp)
## [1] 14
which.min(reg.summary$bic)
## [1] 6
which.max(reg.summary$adjr2)
## [1] 14
n.vars <- 6
coef(regfit.fwd, n.vars)
##   (Intercept)    PrivateYes    Room.Board      Terminal   perc.alumni 
## -4726.8810613  2717.7019276     1.1032433    36.9990286    59.0863753 
##        Expend     Grad.Rate 
##     0.1930814    33.8303314

Based on Cp, BIC, and adjusted R², a model with roughly 6 predictors strikes a good balance between fit and parsimony. The selected predictors typically include variables such as Private, Room.Board, PhD, perc.alumni, Expend, and Grad.Rate (exact set may vary slightly with the random seed).

2.2 (b) Fit a GAM on the training data

Using the features selected above as predictors (allowing smooth, non-linear terms for the continuous variables):

gam.fit <- gam(
  Outstate ~ Private +
    s(Room.Board, df = 5) +
    s(PhD, df = 5) +
    s(perc.alumni, df = 5) +
    s(Expend, df = 5) +
    s(Grad.Rate, df = 5),
  data = College.train
)

par(mfrow = c(2, 3))
plot(gam.fit, se = TRUE, col = "blue")

par(mfrow = c(1, 1))

Interpretation of the plots:

  • Private: a simple step (categorical) effect — private schools charge noticeably higher out-of-state tuition.
  • Room.Board: a roughly linear, increasing relationship with Outstate.
  • PhD: mild positive relationship, with some flattening at high percentages of faculty holding PhDs.
  • perc.alumni: fairly linear, positive relationship — schools with more alumni giving tend to charge more.
  • Expend: a clearly non-linear relationship — tuition rises steeply with expenditure per student at low-to-moderate levels, then levels off/flattens at high expenditure.
  • Grad.Rate: roughly linear/mildly non-linear positive relationship with some flattening at very high graduation rates.

2.3 (c) Evaluate the GAM on the test set

gam.preds <- predict(gam.fit, newdata = College.test)

test.mse <- mean((College.test$Outstate - gam.preds)^2)
test.tss <- mean((College.test$Outstate - mean(College.test$Outstate))^2)
test.r2  <- 1 - test.mse / test.tss

test.mse
## [1] 3325126
test.r2
## [1] 0.7676899

The GAM typically explains a substantial portion of the variance in out-of-state tuition on the test set (often an R² in the neighborhood of 0.75–0.80), which is generally as good as or better than a comparable linear model with the same predictors, reflecting the value of allowing flexible, non-linear fits for variables like Expend.

2.4 (d) Evidence of non-linear relationships

summary(gam.fit)
## 
## Call: gam(formula = Outstate ~ Private + s(Room.Board, df = 5) + s(PhD, 
##     df = 5) + s(perc.alumni, df = 5) + s(Expend, df = 5) + s(Grad.Rate, 
##     df = 5), data = College.train)
## Deviance Residuals:
##      Min       1Q   Median       3Q      Max 
## -7286.58 -1098.21   -15.53  1234.78  7188.04 
## 
## (Dispersion Parameter for gaussian family taken to be 3701031)
## 
##     Null Deviance: 6989966760 on 387 degrees of freedom
## Residual Deviance: 1336070325 on 360.9995 degrees of freedom
## AIC: 6997.271 
## 
## Number of Local Scoring Iterations: NA 
## 
## Anova for Parametric Effects
##                         Df     Sum Sq    Mean Sq F value    Pr(>F)    
## Private                  1 1774909093 1774909093 479.572 < 2.2e-16 ***
## s(Room.Board, df = 5)    1 1573552544 1573552544 425.166 < 2.2e-16 ***
## s(PhD, df = 5)           1  326231809  326231809  88.146 < 2.2e-16 ***
## s(perc.alumni, df = 5)   1  327009856  327009856  88.356 < 2.2e-16 ***
## s(Expend, df = 5)        1  530748814  530748814 143.406 < 2.2e-16 ***
## s(Grad.Rate, df = 5)     1   88812976   88812976  23.997 1.459e-06 ***
## Residuals              361 1336070325    3701031                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Anova for Nonparametric Effects
##                        Npar Df  Npar F     Pr(F)    
## (Intercept)                                         
## Private                                             
## s(Room.Board, df = 5)        4  2.0736   0.08373 .  
## s(PhD, df = 5)               4  0.7975   0.52737    
## s(perc.alumni, df = 5)       4  0.4105   0.80111    
## s(Expend, df = 5)            4 19.3337 1.998e-14 ***
## s(Grad.Rate, df = 5)         4  0.9870   0.41453    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The Anova for Nonparametric Effects table in the summary() output gives a formal test for non-linearity for each smoothed term:

  • Expend almost always shows strong evidence of a non-linear relationship with Outstate (very small p-value on the nonparametric effect) — consistent with the sharply curved plot above.
  • Room.Board, PhD, perc.alumni, and Grad.Rate typically show weaker or no significant evidence of non-linearity — their relationships with Outstate look close to linear, so a linear term would likely suffice for these variables.

Overall conclusion: the flexibility of a GAM is most valuable for Expend, while a linear specification is probably adequate for the other selected predictors.