In this exercise, you will further analyze the Wage data set considered throughout this chapter.
library(ISLR2)
library(boot)
set.seed(1)
cv.error <- rep(0, 5)
for (i in 1:5) {
glm.fit <- glm(wage ~ poly(age, i), data = Wage)
cv.error[i] <- cv.glm(Wage, glm.fit, K = 10)$delta[1]
}
print(cv.error)
## [1] 1676.826 1600.763 1598.399 1595.651 1594.977
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)
anova(fit.1, fit.2, fit.3, fit.4)
## 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)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2998 5022216
## 2 2997 4793430 1 228786 143.6025 < 2.2e-16 ***
## 3 2996 4777674 1 15756 9.8894 0.001679 **
## 4 2995 4771604 1 6070 3.8101 0.051039 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(wage ~ age, data = Wage, col = "darkgrey")
age.grid <- seq(min(Wage$age), max(Wage$age))
preds <- predict(fit.4, newdata = list(age = age.grid))
lines(age.grid, preds, col = "blue", lwd = 2)
The degree-5 polynomial produced the lowest cross-validation error, indicating the best predictive performance. However, the improvement beyond degree 3 was very small, suggesting that higher-degree models provide only marginal gains.
The ANOVA results showed that adding quadratic and cubic terms significantly improved the model (p-values < 0.001), while the fourth-degree term was only marginally significant (p ≈ 0.05). Therefore, ANOVA favors a third-degree polynomial, whereas cross-validation slightly favors a fifth-degree polynomial. Overall, a third-degree polynomial offers a good balance between prediction accuracy, simplicity, and interpretability.
cv.error.step <- rep(0, 10)
for (i in 2:10) {
Wage$age.cut <- cut(Wage$age, i)
glm.fit <- glm(wage ~ age.cut, data = Wage)
cv.error.step[i] <- cv.glm(Wage, glm.fit, K = 10)$delta[1]
}
print(cv.error.step)
## [1] 0.000 1734.999 1682.289 1636.768 1632.149 1624.180 1611.908 1601.480
## [9] 1611.079 1604.125
A step-function regression model was fitted by dividing age into several intervals and using 10-fold cross-validation to determine the optimal number of groups. The cross-validation error decreased as the number of intervals increased, reaching its minimum at approximately six age groups. Beyond six groups, the prediction error increased, indicating that additional intervals caused overfitting.
Unlike polynomial regression, which produces a smooth curve, the step-function model predicts a constant wage within each age interval and changes only at the interval boundaries. The resulting plot therefore appears as a series of horizontal steps rather than a continuous curve.
Although the step-function model captures broad differences in wages across age groups, polynomial regression generally provides a better fit because wages change gradually with age rather than abruptly. Overall, the polynomial model offers slightly better predictive performance and a more realistic representation of the relationship between age and wage.
This question relates to the College data set.
library(leaps)
## Warning: package 'leaps' was built under R version 4.6.1
library(gam)
## Warning: package 'gam' was built under R version 4.6.1
## Loading required package: splines
## Loading required package: foreach
## Warning: package 'foreach' was built under R version 4.6.1
## Loaded gam 1.22-7
data(College)
set.seed(123)
train_idx <- sample(1:nrow(College), 0.75 * nrow(College))
college_train <- College[train_idx, ]
college_test <- College[-train_idx, ]
regfit_fwd <- regsubsets(Outstate ~ ., data = college_train, nvmax = 17, method = "forward")
reg_summary <- summary(regfit_fwd)
best_size <- which.min(reg_summary$bic)
cat("Optimal number of predictors chosen by BIC:", best_size, "\n")
## Optimal number of predictors chosen by BIC: 10
coef(regfit_fwd, best_size)
## (Intercept) PrivateYes Accept Enroll Top25perc
## -2082.6361943 2978.7433045 0.4798842 -1.0568538 15.7230189
## Room.Board Personal PhD perc.alumni Expend
## 0.7556423 -0.3128766 31.4883824 48.9685736 0.2009429
## Grad.Rate
## 17.0082561
gam_model <- gam(Outstate ~ Private + s(Room.Board, df = 4) + s(PhD, df = 4) +
s(perc.alumni, df = 4) + s(Expend, df = 4) + s(Grad.Rate, df = 4),
data = college_train)
par(mfrow = c(2, 3))
plot(gam_model, se = TRUE, col = "blue")
Continuous predictors were modeled using smoothing splines with four degrees of freedom, while Private was included as a parametric factor.
The resulting partial dependence plots reveal several nonlinear relationships. Room.Board, Expend, and Grad.Rate exhibit clear curvature, while PhD and perc.alumni show more modest nonlinear patterns.
gam_preds <- predict(gam_model, newdata = college_test)
test_mse <- mean((college_test$Outstate - gam_preds)^2)
cat("GAM Test MSE:", test_mse, "\n")
## GAM Test MSE: 3316231
rss <- sum((college_test$Outstate - gam_preds)^2)
tss <- sum((college_test$Outstate - mean(college_test$Outstate))^2)
test_r2 <- 1 - (rss / tss)
cat("GAM Test R-squared:", test_r2, "\n")
## GAM Test R-squared: 0.7949882
The high test R2 indicates that the model explains approximately 79% of the variation in out‑of‑state tuition for unseen data, while the MSE reflects strong predictive accuracy. Together, these metrics demonstrate that the GAM generalizes well without evidence of overfitting.
summary(gam_model)
##
## Call: gam(formula = Outstate ~ Private + s(Room.Board, df = 4) + s(PhD,
## df = 4) + s(perc.alumni, df = 4) + s(Expend, df = 4) + s(Grad.Rate,
## df = 4), data = college_train)
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -7458.30 -1192.99 63.72 1253.03 7586.55
##
## (Dispersion Parameter for gaussian family taken to be 3548309)
##
## Null Deviance: 9382434662 on 581 degrees of freedom
## Residual Deviance: 1987054487 on 560.0003 degrees of freedom
## AIC: 10452.93
##
## Number of Local Scoring Iterations: NA
##
## Anova for Parametric Effects
## Df Sum Sq Mean Sq F value Pr(>F)
## Private 1 2574993073 2574993073 725.696 < 2.2e-16 ***
## s(Room.Board, df = 4) 1 1699351821 1699351821 478.919 < 2.2e-16 ***
## s(PhD, df = 4) 1 626090053 626090053 176.447 < 2.2e-16 ***
## s(perc.alumni, df = 4) 1 376899217 376899217 106.219 < 2.2e-16 ***
## s(Expend, df = 4) 1 720605349 720605349 203.084 < 2.2e-16 ***
## s(Grad.Rate, df = 4) 1 89023079 89023079 25.089 7.349e-07 ***
## Residuals 560 1987054487 3548309
## ---
## 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 = 4) 3 3.1708 0.02397 *
## s(PhD, df = 4) 3 2.5325 0.05621 .
## s(perc.alumni, df = 4) 3 1.0622 0.36464
## s(Expend, df = 4) 3 29.6028 < 2e-16 ***
## s(Grad.Rate, df = 4) 3 2.9289 0.03316 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
There is strong statistical evidence of a non-linear relationship with out-of-state tuition for Expend, Room.Board, and Grad.Rate.