Chris Serrano - Assignment 6

Question 6

a

library(ISLR2)
library(boot)

set.seed(1)
cv.errors <- numeric(10)

for (i in 1:10) {
  glm.fit <- glm(wage ~ poly(age, i), data = Wage)
  cv.errors[i] <- cv.glm(Wage, glm.fit, K = 10)$delta[1]
}

best_d <- which.min(cv.errors)
cat("Optimal Polynomial Degree selected by CV:", best_d, "\n")

Optimal Polynomial Degree selected by CV: 9 

plot(1:10, cv.errors, type = "b", xlab = "Polynomial Degree", ylab = "10-Fold CV MSE",
     main = "10-Fold CV Error by Polynomial Degree", col = "royalblue", pch = 19, lwd = 2)
points(best_d, cv.errors[best_d], col = "red", cex = 2, pch = 20)

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)

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)
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

Cross validation selects degree 9. However, the ANOVA test shows that a degree 3 or degree 4 polynomial is good, as there is no statistically significant improvement for degree 5.

fit.poly <- lm(wage ~ poly(age, 4), data = Wage)
age.range <- range(Wage$age)
age.grid <- seq(from = age.range[1], to = age.range[2])
preds <- predict(fit.poly, newdata = list(age = age.grid), se = TRUE)
se.bands <- cbind(preds$fit + 2 * preds$se.fit, preds$fit - 2 * preds$se.fit)

plot(Wage$age, Wage$wage, xlim = age.range, cex = 0.5, col = "darkgrey",
     main = "Degree-4 Polynomial Fit of Wage vs. Age", xlab = "Age", ylab = "Wage")
lines(age.grid, preds$fit, lwd = 3, col = "royalblue")

matlines(age.grid, se.bands, lwd = 1.5, col = "royalblue", lty = 3)

b

set.seed(1)
cv.errors.step <- numeric(10)

for (i in 2:10) {
  Wage$age.cut <- cut(Wage$age, i)
  glm.fit <- glm(wage ~ age.cut, data = Wage)
  cv.errors.step[i] <- cv.glm(Wage, glm.fit, K = 10)$delta[1]
}

best_cuts <- which.min(cv.errors.step[-1]) + 1
cat("Optimal Number of Cuts selected by CV:", best_cuts, "\n")

Optimal Number of Cuts selected by CV: 8 

plot(2:10, cv.errors.step[-1], type = "b", xlab = "Number of Cuts", ylab = "10-Fold CV MSE",
     main = "10-Fold CV Error by Step Function Cuts", col = "forestgreen", pch = 19, lwd = 2)
points(best_cuts, cv.errors.step[best_cuts], col = "red", cex = 2, pch = 20)

Cross-validation selects 8 cuts.

fit.step <- lm(wage ~ cut(age, best_cuts), data = Wage)
preds.step <- predict(fit.step, newdata = list(age = age.grid), se = TRUE)
se.bands.step <- cbind(preds.step$fit + 2 * preds.step$se.fit, preds.step$fit - 2 * preds.step$se.fit)

plot(Wage$age, Wage$wage, xlim = age.range, cex = 0.5, col = "darkgrey",
     main = paste(best_cuts, "-Cut Step Function Fit of Wage vs. Age"), xlab = "Age", ylab = "Wage")
lines(age.grid, preds.step$fit, lwd = 3, col = "forestgreen")

matlines(age.grid, se.bands.step, lwd = 1.5, col = "forestgreen", lty = 3)

---

Question 10

a

library(leaps)
library(gam)

set.seed(1)
train_idx <- sample(nrow(College), nrow(College) / 2)
train_set <- College[train_idx, ]
test_set <- College[-train_idx, ]

fit.fwd <- regsubsets(Outstate ~ ., data = train_set, method = "forward", nvmax = 17)
fwd.summary <- summary(fit.fwd)

par(mfrow = c(1, 3))

# Cp
plot(fwd.summary$cp, xlab = "Number of variables", ylab = "Cp", type = "l")
min.cp <- min(fwd.summary$cp)
std.cp <- sd(fwd.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)

# BIC
plot(fwd.summary$bic, xlab = "Number of variables", ylab = "BIC", type = "l")

min.bic <- min(fwd.summary$bic)
std.bic <- sd(fwd.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)

# Adjusted R2
plot(fwd.summary$adjr2, xlab = "Number of variables", ylab = "Adjusted R2", type = "l", ylim = c(0.4, 0.84))
max.adjr2 <- max(fwd.summary$adjr2)
std.adjr2 <- sd(fwd.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)

Plots show that size 6 is the minimum size for which the scores are within 0.2 standard deviations of the optimum.

coeffs <- coef(fit.fwd, id = 6)
print(coeffs)

  (Intercept)    PrivateYes    Room.Board      Terminal   perc.alumni        Expend 
-4726.8810613  2717.7019276     1.1032433    36.9990286    59.0863753     0.1930814 
    Grad.Rate 
   33.8303314 

The selected features are: PrivateYes, Room.Board, PhD, perc.alumni, Expend, and Grad.Rate.

b

gam.fit <- gam(Outstate ~ Private + 
                 s(Room.Board, df = 2) +
                 s(PhD, df = 2) +
                 s(perc.alumni, df = 2) +
                 s(Expend, df = 5) +
                 s(Grad.Rate, df = 2),
               data = train_set)

par(mfrow = c(2, 3))
plot(gam.fit, se = TRUE, col = "blue")

c

gam.pred <- predict(gam.fit, newdata = test_set)
err <- mean((test_set$Outstate - gam.pred)^2)
cat("Test MSE:", err, "\n")

Test MSE: 3349290 

tss <- mean((test_set$Outstate - mean(test_set$Outstate))^2)
r2 <- 1 - err / tss
cat("Test R2:", r2, "\n")

Test R2: 0.7660016 

We obtain a test R2 of 0.77 using a GAM with 6 predictors.

d

summary(gam.fit)

Call: gam(formula = Outstate ~ Private + s(Room.Board, df = 2) + s(PhD, 
    df = 2) + s(perc.alumni, df = 2) + s(Expend, df = 5) + s(Grad.Rate, 
    df = 2), data = train_set)
Deviance Residuals:
     Min       1Q   Median       3Q      Max 
-7402.89 -1114.45   -12.67  1282.69  7470.60 

(Dispersion Parameter for gaussian family taken to be 3711182)

    Null Deviance: 6989966760 on 387 degrees of freedom
Residual Deviance: 1384271126 on 373 degrees of freedom
AIC: 6987.021 

Number of Local Scoring Iterations: NA 

Anova for Parametric Effects
                        Df     Sum Sq    Mean Sq F value    Pr(>F)    
Private                  1 1778718277 1778718277 479.286 < 2.2e-16 ***
s(Room.Board, df = 2)    1 1577115244 1577115244 424.963 < 2.2e-16 ***
s(PhD, df = 2)           1  322431195  322431195  86.881 < 2.2e-16 ***
s(perc.alumni, df = 2)   1  336869281  336869281  90.771 < 2.2e-16 ***
s(Expend, df = 5)        1  530538753  530538753 142.957 < 2.2e-16 ***
s(Grad.Rate, df = 2)     1   86504998   86504998  23.309 2.016e-06 ***
Residuals              373 1384271126    3711182                      
---
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 = 2)        1  1.9157    0.1672    
s(PhD, df = 2)               1  0.9699    0.3253    
s(perc.alumni, df = 2)       1  0.1859    0.6666    
s(Expend, df = 5)            4 20.5075 2.665e-15 ***
s(Grad.Rate, df = 2)         1  0.5702    0.4506    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

ANOVA shows evidence of a non linear relationship between Outstate and Expend and moderately strong evidence of a non linear relationship between Outstate and Grad.Rate.