6a

library(ISLR2)
## Warning: package 'ISLR2' was built under R version 4.5.3
library(caret)
## Warning: package 'caret' was built under R version 4.5.3
## Loading required package: ggplot2
## Loading required package: lattice
data(Wage)
set.seed(1)


ctrl <- trainControl(method = "cv", number = 10)

cv.errors <- rep(NA, 10)

for (d in 1:10) {
  poly.data <- data.frame(wage = Wage$wage, poly(Wage$age, d))
  
  set.seed(1)
  fit <- train(wage ~ ., data = poly.data,
               method = "lm", trControl = ctrl)
  cv.errors[d] <- fit$results$RMSE^2
}

Plot polynomial fit.

plot(1:10, cv.errors, type = "b", xlab = "Degree", ylab = "CV MSE")

best.d <- which.min(cv.errors)
best.d
## [1] 6

Look at anova hypothesis testing for comparison.

fit.1 <- lm(wage ~ poly(age, 1), 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 ~ 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)
## 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

6b Stepwise and Cross Validation

cv.errors2 <- rep(NA, 9)
for (c in 2:10) {
  Wage$age.cut <- cut(Wage$age, c)
  set.seed(1)
  fit <- train(wage ~ age.cut, data = Wage,
               method = "lm", trControl = ctrl)
  cv.errors2[c - 1] <- fit$results$RMSE^2
}

plot(2:10, cv.errors2, type = "b", xlab = "Number of cuts", ylab = "CV MSE")

best.cuts <- which.min(cv.errors2) + 1
best.cuts
## [1] 8

10a & b

library(mgcv)
## Loading required package: nlme
## This is mgcv 1.9-3. For overview type 'help("mgcv-package")'.
train.idx <- createDataPartition(College$Outstate, p = 0.5, list = FALSE)
College.train <- College[train.idx, ]
College.test  <- College[-train.idx, ]

gam.vars <- c("Private", "Room.Board", "PhD", "perc.alumni", "Expend", "Grad.Rate")

train.sub <- College.train[, c("Outstate", gam.vars)]
test.sub  <- College.test[, c("Outstate", gam.vars)]

ctrl <- trainControl(method = "cv", number = 10)

gam.grid <- data.frame(select = FALSE, method = "GCV.Cp")

set.seed(1)

gam.caret <- train(Outstate ~ Private + Room.Board + PhD +
                     perc.alumni + Expend + Grad.Rate,
                   data = train.sub,
                   method = "gam",
                   trControl = ctrl,
                   tuneGrid = gam.grid)

gam.caret
## Generalized Additive Model using Splines 
## 
## 389 samples
##   6 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 351, 349, 349, 352, 349, 351, ... 
## Resampling results:
## 
##   RMSE      Rsquared   MAE     
##   1889.452  0.7788583  1469.138
## 
## Tuning parameter 'select' was held constant at a value of FALSE
## 
## Tuning parameter 'method' was held constant at a value of GCV.Cp
plot(gam.caret$finalModel, se = TRUE, pages = 1)

10c

The GAM achieves a test R squared of 0.76 (MSE ≈ 3.89 million), which shows good predictive performance and validates that the non-linear terms are adding value beyond what a linear model would do.

preds <- predict(gam.caret, newdata = test.sub)
gam.mse <- mean((test.sub$Outstate - preds)^2)
gam.mse
## [1] 3768141
tss <- mean((test.sub$Outstate - mean(test.sub$Outstate))^2)
test.r2 <- 1 - gam.mse / tss
test.r2
## [1] 0.7699125

10d

Phd, grad.rate, room.board, and expend all have nonlinear relationships.

summary(gam.caret$finalModel)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## .outcome ~ PrivateYes + s(perc.alumni) + s(PhD) + s(Grad.Rate) + 
##     s(Room.Board) + s(Expend)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   8617.0      225.6  38.204   <2e-16 ***
## PrivateYes    2571.4      286.2   8.984   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                  edf Ref.df      F  p-value    
## s(perc.alumni) 2.547  3.220  6.708 0.000141 ***
## s(PhD)         4.631  5.674  2.716 0.014447 *  
## s(Grad.Rate)   1.000  1.000 12.479 0.000463 ***
## s(Room.Board)  1.331  1.595 21.651 2.13e-07 ***
## s(Expend)      5.414  6.572 19.948  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.795   Deviance explained = 80.3%
## GCV = 3.429e+06  Scale est. = 3.2798e+06  n = 389