Question 6

Further Analysis of the Wage Data Set

a.) Polynomial regression with cross‑validation

Perform a polynomial regression of wage on age and use cross‑validation to determine the optimal polynomial degree 𝑑. Report which degree is selected and compare this result with the conclusions obtained from an ANOVA hypothesis test comparing nested polynomial models. Finally, produce a plot showing the fitted polynomial curve along with the observed data.

library(ISLR2)
library(caret)
library(ggplot2)

data(Wage)

set.seed(1)

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

poly.models <- list()
cv.errors <- numeric(10)

for(i in 1:10) {
  poly.formula <- as.formula(
    paste("wage~poly(age,", i, ")")
  )
  poly.models[[i]] <- train(
    poly.formula, 
    data = Wage, 
    method = "lm",
    trControl = ctrl
  )
  cv.errors[i] <- min(poly.models[[i]]$results$RMSE^2)
}
cv.errors

best.degree <- which.min(poly.models[[i]]$results$RMSE^2)
cat("Degree Chosen:", best.degree, "\n")

poly.fit <- lm(wage ~ poly(age, best.degree), data = Wage)
summary(poly.fit)

#Compare to anova
fit1 <- lm(wage ~ age, 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)

#plot
age.grid <- data.frame(age = seq(min(Wage$age), 
                                 max(Wage$age), 
                                 length = 100))
pred <- predict(poly.fit, newdata = age.grid)

ggplot(Wage, aes(age, wage))+
  geom_point(alpha = .4)+
  geom_line(data = age.grid,
            aes(age, pred),
            color = "blue",
            linewidth = 1.2)+
  labs(title = paste("Polynomial Regression (Degree", best.degree, ")"))
##  [1] 1670.277 1596.717 1582.212 1581.502 1584.812 1578.691 1591.044 1585.468
##  [9] 1586.835 1583.362
## Degree Chosen: 1 
## 
## Call:
## lm(formula = wage ~ poly(age, best.degree), data = Wage)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -100.265  -25.115   -6.063   16.601  205.748 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            111.7036     0.7473  149.48   <2e-16 ***
## poly(age, best.degree) 447.0679    40.9291   10.92   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 40.93 on 2998 degrees of freedom
## Multiple R-squared:  0.03827,    Adjusted R-squared:  0.03795 
## F-statistic: 119.3 on 1 and 2998 DF,  p-value: < 2.2e-16

b.) Step function with cross‑validation

Fit a step function model for predicting wage from age, and use cross‑validation to determine the optimal number of age intervals (cuts). Then, generate a plot showing the fitted step function alongside the observed data.

cuts <- 2:10
step.models <- list()
step.errors <- numeric(length(cuts))

for(i in seq_along(cuts)){
  step.formula <- as.formula(
    paste("wage ~ cut(age,", cuts[i], ")")
  )
  step.models[[i]] <- train(
    step.formula,
    data = Wage,
    method = "lm", 
    trControl = ctrl
  )
  step.errors[i] <- min(step.models[[i]]$results$RMSE^2)
}
step.errors

best.cuts <- cuts[which.min(step.errors)]
cat("Best cuts:", best.cuts, "\n")

step.fit <- lm(
  wage ~ cut(age, best.cuts),
  data = Wage
)

summary(step.fit)

#plot

step.pred <-  predict(
  step.fit, 
  newdata = data.frame(age = age.grid)
)

step.plot <- data.frame(
  age = age.grid, 
  wage = step.pred
)

ggplot(Wage, aes(age, wage)) +
  geom_point(alpha = .4) +
  geom_line(
    data = step.plot,
    aes(age, wage), 
    color = "blue",
    linewidth = 1.2
  ) +
  labs(title = paste("Step Function (", best.cuts, "cuts)", sep = ""),
       xlab = "Age",
       ylab = "Wage")
## [1] 1729.420 1677.387 1628.351 1615.789 1619.801 1597.982 1599.407 1603.159
## [9] 1602.212
## Best cuts: 7 
## 
## Call:
## lm(formula = wage ~ cut(age, best.cuts), data = Wage)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -98.571 -24.584  -4.956  15.602 207.239 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      80.064      2.405  33.288  < 2e-16 ***
## cut(age, best.cuts)(26.9,35.7]   24.516      2.893   8.476  < 2e-16 ***
## cut(age, best.cuts)(35.7,44.6]   38.593      2.792  13.820  < 2e-16 ***
## cut(age, best.cuts)(44.6,53.4]   38.055      2.812  13.531  < 2e-16 ***
## cut(age, best.cuts)(53.4,62.3]   39.040      3.082  12.667  < 2e-16 ***
## cut(age, best.cuts)(62.3,71.1]   31.039      4.884   6.355  2.4e-10 ***
## cut(age, best.cuts)(71.1,80.1]   15.994      9.076   1.762   0.0781 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 40.1 on 2993 degrees of freedom
## Multiple R-squared:  0.07824,    Adjusted R-squared:  0.07639 
## F-statistic: 42.34 on 6 and 2993 DF,  p-value: < 2.2e-16

Question 10

This Question relates to the College data set

a.) Stepwise Selection

Split the College data set into a training set and a test set. Using Outstate (out‑of‑state tuition) as the response variable and all remaining variables as predictors, perform forward stepwise selection on the training set to identify a parsimonious model that includes only a subset of the available predictors.

library(ISLR2)
library(caret)
library(MASS)
library(gam)

data(College)

set.seed(1)

train.index <- createDataPartition(College$Outstate,
                                   p = 0.7,
                                   list = FALSE)

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

full.fit <- lm(Outstate ~ ., data = College.train)

step.fit <- stepAIC(full.fit,
                    direction = "both",
                    trace = FALSE)

summary(step.fit)

b.) GAM

Fit a generalized additive model (GAM) to the training data, using Outstate (out‑of‑state tuition) as the response and the subset of predictors selected in the previous step as the explanatory variables. Produce plots of the fitted smooth functions and provide an explanation of the patterns revealed in the model.

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

summary(gam.fit)

par(mfrow = c(2,3))

plot(gam.fit,
     se = TRUE,
     col = "blue")
## 
## Call: gam(formula = Outstate ~ Private + lo(Room.Board) + lo(Expend) + 
##     lo(PhD) + lo(perc.alumni) + lo(Grad.Rate), data = College.train)
## Deviance Residuals:
##    Min     1Q Median     3Q    Max 
##  -6785  -1139     62   1229   4768 
## 
## (Dispersion Parameter for gaussian family taken to be 3242650)
## 
##     Null Deviance: 8648835880 on 545 degrees of freedom
## Residual Deviance: 1700199687 on 524.3242 degrees of freedom
## AIC: 9758.293 
## 
## Number of Local Scoring Iterations: NA 
## 
## Anova for Parametric Effects
##                     Df     Sum Sq    Mean Sq F value    Pr(>F)    
## Private           1.00 2619605560 2619605560 807.860 < 2.2e-16 ***
## lo(Room.Board)    1.00 1684465323 1684465323 519.472 < 2.2e-16 ***
## lo(Expend)        1.00 1284709139 1284709139 396.191 < 2.2e-16 ***
## lo(PhD)           1.00  174755812  174755812  53.893 8.128e-13 ***
## lo(perc.alumni)   1.00  160221174  160221174  49.411 6.494e-12 ***
## lo(Grad.Rate)     1.00   78079935   78079935  24.079 1.236e-06 ***
## Residuals       524.32 1700199687    3242650                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Anova for Nonparametric Effects
##                 Npar Df  Npar F    Pr(F)    
## (Intercept)                                 
## Private                                     
## lo(Room.Board)      3.0  3.4006  0.01816 *  
## lo(Expend)          4.0 20.3418 1.11e-15 ***
## lo(PhD)             2.8  0.9328  0.41989    
## lo(perc.alumni)     2.3  0.7208  0.50664    
## lo(Grad.Rate)       2.6  2.0513  0.11582    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

c.) Results

Evaluate the fitted model on the test set and discuss the results, commenting on its predictive performance and any insights gained from comparing test‑set error to the training‑set behavior.

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

gam.mse <- mean((gam.pred - College.test$Outstate)^2)

gam.rmse <- sqrt(gam.mse)

gam.mse
gam.rmse
## [1] 4193772
## [1] 2047.87

The fitted model produced a mean squared error (MSE) of 4,193,773, corresponding to a root mean squared error (RMSE) of 2,047. This indicates that, on average, the model’s predicted out‑of‑state tuition values differ from the actual tuition values by approximately $2,047.

d.) Non-linearity

Identify which predictors, if any, exhibit evidence of a non‑linear relationship with the response variable.

summary(gam.fit)
## 
## Call: gam(formula = Outstate ~ Private + lo(Room.Board) + lo(Expend) + 
##     lo(PhD) + lo(perc.alumni) + lo(Grad.Rate), data = College.train)
## Deviance Residuals:
##    Min     1Q Median     3Q    Max 
##  -6785  -1139     62   1229   4768 
## 
## (Dispersion Parameter for gaussian family taken to be 3242650)
## 
##     Null Deviance: 8648835880 on 545 degrees of freedom
## Residual Deviance: 1700199687 on 524.3242 degrees of freedom
## AIC: 9758.293 
## 
## Number of Local Scoring Iterations: NA 
## 
## Anova for Parametric Effects
##                     Df     Sum Sq    Mean Sq F value    Pr(>F)    
## Private           1.00 2619605560 2619605560 807.860 < 2.2e-16 ***
## lo(Room.Board)    1.00 1684465323 1684465323 519.472 < 2.2e-16 ***
## lo(Expend)        1.00 1284709139 1284709139 396.191 < 2.2e-16 ***
## lo(PhD)           1.00  174755812  174755812  53.893 8.128e-13 ***
## lo(perc.alumni)   1.00  160221174  160221174  49.411 6.494e-12 ***
## lo(Grad.Rate)     1.00   78079935   78079935  24.079 1.236e-06 ***
## Residuals       524.32 1700199687    3242650                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Anova for Nonparametric Effects
##                 Npar Df  Npar F    Pr(F)    
## (Intercept)                                 
## Private                                     
## lo(Room.Board)      3.0  3.4006  0.01816 *  
## lo(Expend)          4.0 20.3418 1.11e-15 ***
## lo(PhD)             2.8  0.9328  0.41989    
## lo(perc.alumni)     2.3  0.7208  0.50664    
## lo(Grad.Rate)       2.6  2.0513  0.11582    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Based on the ANOVA test for nonparametric effects, Expend (p = 1.1e‑15) and Room.Board (p = 0.0182) exhibit evidence of non‑linear relationships with the response. Among these, Expend shows particularly strong evidence of non‑linearity.