Chapter 7

Question 6

In this exercise, you will further analyze the Wage data set considered throughout this chapter.

library(ISLR2)
library(boot)
  1. Perform polynomial regression to predict wage using age. Use cross-validation to select the optimal degree d for the polynomial. What degree was chosen, and how does this compare to the results of hypothesis testing using ANOVA? Make a plot of the resulting polynomial fit to the data.
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.

  1. Fit a step function to predict wage using age, and perform cross validation to choose the optimal number of cuts. Make a plot of the fit obtained.
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.

Question 10

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)
  1. Split the data into a training set and a test set. Using out-of-state tuition as the response and the other variables as the predictors, perform forward stepwise selection on the training set in order to identify a satisfactory model that uses just a subset of the predictors.
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
  1. Fit a GAM on the training data, using out-of-state tuition as the response and the features selected in the previous step as the predictors. Plot the results, and explain your findings.
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.

  1. Evaluate the model obtained on the test set, and explain the results obtained.
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.

  1. For which variables, if any, is there evidence of a non-linear relationship with the response?
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.