Assignment 6
STA-6543 Predictive Modeling — Chapter 7: An Introduction to Statistical Learning with Applications in R
Jerimey Simons
July 15, 2026
Question 6.
In this exercise, you will further analyze the Wage data set considered throughout this chapter.
library(ISLR2)
library(caret)
library(ggplot2)
attach(Wage)
set.seed(42)(a) 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.
library(ISLR2)
library(caret)
set.seed(42)
cv_folds <- createFolds(
Wage$wage,
k = 5,
returnTrain = TRUE
)
cv_control <- trainControl(
method = "cv",
number = 5,
index = cv_folds
)
cv_results <- data.frame(
Degree = 0:10,
CV_MSE = NA_real_
)
for (d in 0:10) {
if (d == 0) {
# Intercept-only model
fold_mse <- numeric(5)
for (j in seq_along(cv_folds)) {
train_rows <- cv_folds[[j]]
validation_rows <- setdiff(seq_len(nrow(Wage)), train_rows)
mean_prediction <- mean(Wage$wage[train_rows])
fold_mse[j] <- mean(
(Wage$wage[validation_rows] - mean_prediction)^2
)
}
cv_results$CV_MSE[d + 1] <- mean(fold_mse)
} else {
x_poly <- as.data.frame(
poly(
Wage$age,
degree = d,
raw = TRUE
)
)
names(x_poly) <- paste0("age_degree_", seq_len(d))
cv_model <- train(
x = x_poly,
y = Wage$wage,
method = "lm",
trControl = cv_control
)
cv_results$CV_MSE[d + 1] <- cv_model$results$RMSE^2
}
}
best_degree <- cv_results$Degree[
which.min(cv_results$CV_MSE)
]
cat("Optimal polynomial degree:", best_degree, "\n")## Optimal polynomial degree: 7knitr::kable(
cv_results,
digits = 2,
caption = "Five-Fold Cross-Validation MSE by Polynomial Degree"
)| Degree | CV_MSE |
|---|---|
| 0 | 1741.01 |
| 1 | 1672.16 |
| 2 | 1596.95 |
| 3 | 1592.26 |
| 4 | 1590.61 |
| 5 | 1590.47 |
| 6 | 1589.36 |
| 7 | 1588.76 |
| 8 | 1589.29 |
| 9 | 1589.11 |
| 10 | 1589.76 |
Five-fold cross-validation was used to compare polynomial regression models with degrees ranging from 0 to 10. The degree-0 model, which predicts the mean wage for all observations, had the largest cross-validation MSE (1738.98). As additional polynomial terms were added, the prediction error decreased substantially through degree 4.
The lowest cross-validation MSE (1588.76) was obtained using a 7th-degree polynomial. However, the improvement beyond a 4th-degree polynomial was very small, indicating that higher-order polynomial terms provide only marginal gains in predictive accuracy. Therefore, although cross-validation selected degree 7, a lower-degree polynomial may be preferred because it provides a simpler and more interpretable model with nearly the same predictive performance.
fit1 <- lm(wage ~ poly(age, 1), data = Wage)
fit2 <- lm(wage ~ poly(age, 2), data = Wage)
fit3 <- lm(wage ~ poly(age, 3), data = Wage)
fit4 <- lm(wage ~ poly(age, 4), data = Wage)
fit5 <- lm(wage ~ poly(age, 5), data = Wage)
anova(fit1, fit2, fit3, fit4, fit5)## Analysis of Variance Table
##
## Model 1: wage ~ poly(age, 1)
## Model 2: wage ~ poly(age, 2)
## Model 3: wage ~ poly(age, 3)
## Model 4: wage ~ poly(age, 4)
## Model 5: wage ~ poly(age, 5)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2998 5022216
## 2 2997 4793430 1 228786 143.5931 < 2.2e-16 ***
## 3 2996 4777674 1 15756 9.8888 0.001679 **
## 4 2995 4771604 1 6070 3.8098 0.051046 .
## 5 2994 4770322 1 1283 0.8050 0.369682
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Final Model
best_model <- lm(wage ~ poly(age, 7), data = Wage)
# Create a sequence of age values
age_grid <- data.frame(
age = seq(
min(Wage$age),
max(Wage$age),
length.out = 100
)
)
# Predicted wages
age_grid$wage_hat <- predict(
best_model,
newdata = age_grid
)
# Plot
ggplot(Wage, aes(x = age, y = wage)) +
geom_point(alpha = 0.30) +
geom_line(
data = age_grid,
aes(y = wage_hat),
color = "blue",
linewidth = 1.2
) +
labs(
title = "Polynomial Regression of Wage on Age",
subtitle = "Degree 7 Polynomial Selected by Cross-Validation",
x = "Age",
y = "Wage"
) +
theme_minimal()
Five-fold cross-validation selected a 7th-degree
polynomial because it achieved the lowest cross-validation mean
squared error (MSE = 1588.76), making it the best model
for prediction. In contrast, the ANOVA results indicated that adding
polynomial terms beyond the 3rd degree did not provide
statistically significant improvements in model fit. Specifically, the
quadratic and cubic terms significantly improved the model, while the
fourth-degree term was only marginally significant (p = 0.051) and the
fifth-degree term was not significant (p = 0.370).
Although the 7th-degree polynomial achieved the lowest prediction error, the improvement over lower-degree models was very small. Therefore, if the goal is prediction, the 7th-degree polynomial is preferred. If the goal is model simplicity and statistical inference, the 3rd-degree polynomial is the more appropriate choice because it is more parsimonious while providing nearly the same predictive performance.
(b) Fit a step function to predict wage using age, and perform crossvalidation to choose the optimal number of cuts. Make a plot of the fit obtained.
# 5-fold cross-validation
ctrl <- trainControl(
method = "cv",
number = 5
)
cv_results <- data.frame(
Cuts = 2:10,
CV_MSE = NA
)
for (k in 2:10) {
Wage$age_cut <- cut(
Wage$age,
breaks = k
)
model <- train(
wage ~ age_cut,
data = Wage,
method = "lm",
trControl = ctrl
)
cv_results$CV_MSE[k - 1] <- model$results$RMSE^2
}
knitr::kable(
cv_results,
digits = 2,
caption = "Five-Fold Cross-Validation Results for Step Function"
)| Cuts | CV_MSE |
|---|---|
| 2 | 1732.83 |
| 3 | 1678.10 |
| 4 | 1632.46 |
| 5 | 1629.05 |
| 6 | 1618.40 |
| 7 | 1607.81 |
| 8 | 1598.10 |
| 9 | 1605.82 |
| 10 | 1603.99 |
best_cuts <- cv_results$Cuts[
which.min(cv_results$CV_MSE)
]
cat("Optimal Number of Cuts:", best_cuts)## Optimal Number of Cuts: 8Wage$age_cut <- cut(
Wage$age,
breaks = best_cuts
)
step_fit <- lm(
wage ~ age_cut,
data = Wage
)
summary(step_fit)##
## Call:
## lm(formula = wage ~ age_cut, data = Wage)
##
## Residuals:
## Min 1Q Median 3Q Max
## -99.697 -24.552 -5.307 15.417 198.560
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 76.282 2.630 29.007 < 2e-16 ***
## age_cut(25.8,33.5] 25.833 3.161 8.172 4.44e-16 ***
## age_cut(33.5,41.2] 40.226 3.049 13.193 < 2e-16 ***
## age_cut(41.2,49] 43.501 3.018 14.412 < 2e-16 ***
## age_cut(49,56.8] 40.136 3.177 12.634 < 2e-16 ***
## age_cut(56.8,64.5] 44.102 3.564 12.373 < 2e-16 ***
## age_cut(64.5,72.2] 28.948 6.042 4.792 1.74e-06 ***
## age_cut(72.2,80.1] 15.224 9.781 1.556 0.12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 39.97 on 2992 degrees of freedom
## Multiple R-squared: 0.08467, Adjusted R-squared: 0.08253
## F-statistic: 39.54 on 7 and 2992 DF, p-value: < 2.2e-16# Prediction grid
age_grid <- data.frame(
age = seq(
min(Wage$age),
max(Wage$age),
length.out = 1000
)
)
age_grid$age_cut <- cut(
age_grid$age,
breaks = best_cuts
)
age_grid$pred <- predict(
step_fit,
newdata = age_grid
)
ggplot(Wage, aes(age, wage)) +
geom_point(alpha = 0.30) +
geom_step(
data = age_grid,
aes(age, pred),
color = "blue",
linewidth = 1.2
) +
labs(
title = paste(
"Step Function Fit (",
best_cuts,
" Cuts)",
sep = ""
),
x = "Age",
y = "Wage"
) +
theme_minimal()
Five-fold cross-validation was used to compare step function models with
the number of cuts ranging from 2 to 10. The model with 8
cuts produced the lowest cross-validation mean squared error
(MSE) of 1595.45, indicating that it provided the best
predictive performance among the models considered.
The cross-validation results show that prediction error decreased substantially as the number of cuts increased from 2 to 8. After 8 cuts, the prediction error increased slightly, suggesting that additional cuts did not improve predictive performance and may have introduced unnecessary model complexity.
The fitted step function shows that average wage generally increases from early adulthood through middle age, reaching its highest levels between approximately 40 and 65 years of age. After age 65, the average wage begins to decline. Most of the age intervals were statistically significant (p < 0.001), indicating meaningful differences in average wage across these age groups. However, the oldest age group (72.2–80.1 years) was not statistically significant (p = 0.12), suggesting that its average wage was not significantly different from that of the reference age group.
Compared with the polynomial regression model, the step function provides a simpler and more interpretable representation of the relationship between age and wage by modeling wage as constant within age intervals. In contrast, polynomial regression produces a smooth nonlinear curve that captures gradual changes in wage with age. Both approaches indicate that wages tend to increase during early and middle adulthood before leveling off and declining at older ages.
Question 10.
This question relates to the College data set.
library(leaps)
library(ISLR2)
attach(College)
set.seed(42)(a) 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.
train <- sample(length(Outstate), length(Outstate) / 2)
test <- -train
College.train <- College[train, ]
College.test <- College[test, ]
fit <- regsubsets(Outstate ~ ., data = College.train, nvmax = 17, method = "forward")
fit.summary <- summary(fit)
par(mfrow = c(1, 3))
plot(fit, scale = "adjr2", main = "Adjusted R²")
plot(fit, scale = "Cp", main = "Mallows' Cp")
plot(fit, scale = "bic", main = "BIC")
best_size <- which.min(fit.summary$bic)
best_size## [1] 10coef(fit, best_size)## (Intercept) PrivateYes Accept Top10perc F.Undergrad
## -2621.3803473 2513.8383008 0.3753840 31.7605086 -0.2069239
## Room.Board Books Terminal perc.alumni Expend
## 0.9932388 -1.3038896 38.7475364 51.2338704 0.1456424
## Grad.Rate
## 17.9418689par(mfrow = c(1, 3))
plot(fit.summary$cp, xlab = "Number of variables", ylab = "Cp", type = "l")
min.cp <- min(fit.summary$cp)
std.cp <- sd(fit.summary$cp)
abline(h = min.cp + 0.2 * std.cp, col = "red", lty = 2)
abline(h = min.cp - 0.2 * std.cp, col = "red", lty = 2)
plot(fit.summary$bic, xlab = "Number of variables", ylab = "BIC", type='l')
min.bic <- min(fit.summary$bic)
std.bic <- sd(fit.summary$bic)
abline(h = min.bic + 0.2 * std.bic, col = "red", lty = 2)
abline(h = min.bic - 0.2 * std.bic, col = "red", lty = 2)
plot(fit.summary$adjr2, xlab = "Number of variables", ylab = "Adjusted R2", type = "l", ylim = c(0.4, 0.84))
max.adjr2 <- max(fit.summary$adjr2)
std.adjr2 <- sd(fit.summary$adjr2)
abline(h = max.adjr2 + 0.2 * std.adjr2, col = "red", lty = 2)
abline(h = max.adjr2 - 0.2 * std.adjr2, col = "red", lty = 2)
The College data set was randomly divided into equal-sized training and
test sets . Forward stepwise selection was then performed on the
training data using Outstate as the response variable
and the remaining variables as predictors. Models containing up to 17
predictors were evaluated.
Model selection was based on Mallows’ Cp, the Bayesian Information Criterion (BIC), and Adjusted \(R^2\). Both Cp and BIC selected a model containing 10 predictors, while the Adjusted \(R^2\) increased rapidly before leveling off after approximately 10 predictors, indicating that additional variables provided little improvement in model fit. Therefore, the 10-predictor model was selected because it provides an excellent balance between predictive accuracy and model complexity.
The final forward stepwise model included the following predictors: Private, Accept, Top10perc, F.Undergrad, Room.Board, Books, Terminal, perc.alumni, Expend, and Grad.Rate. These variables were identified as the most informative for predicting out-of-state tuition while excluding less influential predictors.
(b) 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.
library(mgcv)
gam_fit <- gam(
Outstate ~
Private +
s(Accept) +
s(Top10perc) +
s(F.Undergrad) +
s(Room.Board) +
s(Books) +
s(Terminal) +
s(perc.alumni) +
s(Expend) +
s(Grad.Rate),
data = College.train
)
summary(gam_fit)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Outstate ~ Private + s(Accept) + s(Top10perc) + s(F.Undergrad) +
## s(Room.Board) + s(Books) + s(Terminal) + s(perc.alumni) +
## s(Expend) + s(Grad.Rate)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9242.7 272.7 33.898 <2e-16 ***
## PrivateYes 1839.4 356.2 5.164 4e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Accept) 4.355 5.459 5.190 8.37e-05 ***
## s(Top10perc) 1.000 1.000 5.018 0.0257 *
## s(F.Undergrad) 8.507 8.921 5.129 2.88e-06 ***
## s(Room.Board) 1.000 1.000 59.081 < 2e-16 ***
## s(Books) 1.954 2.486 2.845 0.0391 *
## s(Terminal) 1.000 1.000 6.088 0.0141 *
## s(perc.alumni) 1.000 1.000 8.347 0.0041 **
## s(Expend) 5.192 6.268 14.553 < 2e-16 ***
## s(Grad.Rate) 1.175 1.329 5.659 0.0091 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.828 Deviance explained = 84%
## GCV = 3.0582e+06 Scale est. = 2.8439e+06 n = 388A generalized additive model (GAM) was fitted to the training data using the 10 predictors selected by forward stepwise selection. The categorical variable Private was included as a linear term, while smoothing splines were used for the remaining continuous predictors. The model explained approximately 84% of the variation in out-of-state tuition (Adjusted \(R^2 = 0.828\); Deviance Explained = 84%), indicating an excellent overall fit.
The GAM results show that all selected predictors were statistically significant at the 5% significance level. The variable Private had a significant positive effect (\(p < 0.001\)), indicating that private colleges generally charge higher out-of-state tuition than public institutions. Several predictors exhibited nonlinear relationships with out-of-state tuition, as indicated by estimated degrees of freedom (edf) greater than one. In particular, Accept (edf = 4.36), F.Undergrad (edf = 8.51), Books (edf = 1.95), Expend (edf = 5.19), and Grad.Rate (edf = 1.18) showed nonlinear effects. In contrast, Top10perc, Room.Board, Terminal, and perc.alumni all had estimated degrees of freedom close to one, suggesting that their relationships with out-of-state tuition are approximately linear.
The GAM plots provide insight into the shape of each predictor’s relationship with out-of-state tuition. The smooth for Accept shows a nonlinear relationship, with tuition increasing more rapidly as the number of accepted students becomes larger before leveling off. F.Undergrad displays the strongest nonlinear pattern, suggesting that undergraduate enrollment influences tuition differently across small and large institutions. Expend also exhibits a pronounced nonlinear positive relationship, indicating that institutions with higher instructional expenditures tend to charge substantially higher out-of-state tuition.
The remaining predictors display nearly linear effects. Room.Board has a strong positive relationship with tuition, while Books shows a slight negative trend. Top10perc, Terminal, perc.alumni, and Grad.Rate exhibit relatively modest but statistically significant effects. Overall, the GAM reveals that only a subset of predictors require nonlinear smoothing, while the remaining variables are adequately modeled with approximately linear relationships.
(c) Evaluate the model obtained on the test set, and explain the results obtained.
gam_pred <- predict(gam_fit, newdata = College.test)
gam_mse <- mean( (College.test$Outstate - gam_pred)^2)
gam_rmse <- sqrt(gam_mse)
cat("Test Error:", gam_mse, "\n")## Test Error: 4393719cat("Root Mean Squared Error:", gam_rmse, "\n")## Root Mean Squared Error: 2096.12The generalized additive model (GAM) was evaluated on the test set by predicting out-of-state tuition for the held-out observations. The model achieved a test mean squared error (MSE) of 4,393,719 and a root mean squared error (RMSE) of 2,096.12.
The relatively low prediction error indicates that the GAM generalized well to the test data and was able to accurately predict out-of-state tuition for unseen observations. Combined with the training performance (Adjusted \(R^2 = 0.828\) and 84% deviance explained), these results suggest that the model fits the data well without substantial overfitting.
The GAM benefited from modeling nonlinear relationships for variables such as Accept, F.Undergrad, and Expend, while treating variables such as Top10perc, Room.Board, Terminal, and perc.alumni as approximately linear. This flexibility allowed the model to capture complex relationships in the data while maintaining good predictive accuracy. Overall, the GAM provides an effective model for predicting out-of-state tuition by balancing model flexibility with generalization to new data.
(d) For which variables, if any, is there evidence of a non-linear relationship with the response?
The estimated degrees of freedom (edf) from the generalized additive model (GAM) were used to determine whether each predictor exhibited a linear or nonlinear relationship with out-of-state tuition. Predictors with an estimated degrees of freedom close to one were considered approximately linear, while predictors with estimated degrees of freedom greater than one indicated evidence of nonlinearity.
The strongest evidence of nonlinear relationships was observed for F.Undergrad (edf = 8.51), Expend (edf = 5.19), and Accept (edf = 4.36), indicating that these variables have complex relationships with out-of-state tuition that are better modeled with smoothing splines than with simple linear terms. Books also showed mild evidence of nonlinearity (edf = 1.95). In contrast, Top10perc, Room.Board, Terminal, and perc.alumni all had estimated degrees of freedom equal to one, suggesting approximately linear relationships with the response. Although Grad.Rate had an estimated degrees of freedom of 1.18, its relationship was only slightly nonlinear and is close enough to one that it can be considered approximately linear.
Overall, the GAM indicates that the primary nonlinear effects are associated with Accept, F.Undergrad, Books, and Expend, while the remaining predictors have relationships that are largely linear.