Question 6
# (6a)
set.seed(1)
cv_errors <- sapply(1:5, function(d) {
fit <- glm(wage ~ poly(age, d), data = Wage)
cv.glm(Wage, fit, K = 10)$delta[1]
})
print(cv_errors)
models <- lapply(1:5, function(d) lm(wage ~ poly(age, d), data = Wage))
do.call(anova, models)
plot(wage ~ age, data = Wage, col = "black")
age_grid <- seq(min(Wage$age), max(Wage$age))
preds <- predict(models[[4]], newdata = list(age = age_grid))
lines(age_grid, preds, col = "blue", lwd = 2)
invisible("The 10-fold cross-validation strictly minimizes test error at degree 5 (1594.977), though the improvements stop after degree 4. The ANOVA table shows that adding the 4th degree is only marginally significant (p = 0.051) and the 5th degree is insignificant (p = 0.370).")
# (6b)
set.seed(1)
folds <- sample(1:10, nrow(Wage), replace = TRUE)
cv_errors <- sapply(2:8, function(c) {
Wage$age_cut <- cut(Wage$age, breaks = c)
errors <- sapply(1:10, function(f) {
fit <- lm(wage ~ age_cut, data = Wage[folds != f, ])
pred <- predict(fit, newdata = Wage[folds == f, ])
mean((Wage[folds == f, ]$wage - pred)^2, na.rm = TRUE)
})
mean(errors)
})
names(cv_errors) <- 2:8
print(cv_errors)
plot(wage ~ age, data = Wage, col = "black")
step_fit <- lm(wage ~ cut(age, 8), data = Wage)
age_grid <- data.frame(age = seq(min(Wage$age), max(Wage$age)))
lines(age_grid$age, predict(step_fit, age_grid), col = "red", lwd = 2)
Question 10
# (10a)
set.seed(1)
train_idx <- sample(1:nrow(College), 0.7 * nrow(College))
train <- College[train_idx, ]
test <- College[-train_idx, ]
forward_fit <- regsubsets(Outstate ~ ., data = train, nvmax = 17, method = "forward")
forward_summary <- summary(forward_fit)
best_adjr2 <- which.max(forward_summary$adjr2)
best_bic <- which.min(forward_summary$bic)
cat("Optimal variables by Adjusted R-squared:", best_adjr2, "\n")
## Optimal variables by Adjusted R-squared: 13
cat("Optimal variables by BIC:", best_bic, "\n")
## Optimal variables by BIC: 13
coef(forward_fit, id = best_bic)
## (Intercept) PrivateYes Apps Accept Top10perc
## -1739.5725417 2276.7996721 -0.3358567 0.7814587 28.9687655
## F.Undergrad Room.Board Personal PhD Terminal
## -0.1559550 0.9134134 -0.3484815 11.8113175 24.9233138
## S.F.Ratio perc.alumni Expend Grad.Rate
## -55.0649149 48.6046652 0.1744677 20.9498491
# (10b)
gam_fit <- gam(Outstate ~ Private + s(Apps) + s(Accept) + s(Top10perc) +
s(F.Undergrad) + s(Room.Board) + s(Personal) + s(PhD) +
s(Terminal) + s(S.F.Ratio) + s(perc.alumni) + s(Expend) +
s(Grad.Rate), data = train)
par(mfrow = c(3, 5))
plot(gam_fit, se = TRUE, col = "blue")
invisible("Private status, higher student quality (Top10perc), and higher institutional investments (Room.Board, perc.alumni) steadily drive up out-of-state tuition. Non-linear threshold effect for instructional spending (Expend), where tuition sharply increases before flattening out past $20,000.")
# (10c)
gam_preds <- predict(gam_fit, newdata = test)
test_mse <- mean((test$Outstate - gam_preds)^2)
cat("Test MSE:", test_mse, "\n")
## Test MSE: 3146155
tss <- mean((test$Outstate - mean(test$Outstate))^2)
test_r2 <- 1 - (test_mse / tss)
cat("Test R-squared:", test_r2, "\n")
## Test R-squared: 0.7720737
invisible("The GAM achieves a Test R square of 0.7721, explains 77.21% of the variance in out-of-state tuition on unseen data. The model's predictions deviate from the true tuition values by approximately $1,774 (3,146,155), demonstrating a reliable fit.")
# (10d)
summary(gam_fit)
##
## Call: gam(formula = Outstate ~ Private + s(Apps) + s(Accept) + s(Top10perc) +
## s(F.Undergrad) + s(Room.Board) + s(Personal) + s(PhD) + s(Terminal) +
## s(S.F.Ratio) + s(perc.alumni) + s(Expend) + s(Grad.Rate),
## data = train)
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5910.22 -1088.99 73.01 1113.20 7156.97
##
## (Dispersion Parameter for gaussian family taken to be 3195196)
##
## Null Deviance: 9260683704 on 542 degrees of freedom
## Residual Deviance: 1575229760 on 492.9994 degrees of freedom
## AIC: 9723.111
##
## Number of Local Scoring Iterations: NA
##
## Anova for Parametric Effects
## Df Sum Sq Mean Sq F value Pr(>F)
## Private 1 2531384624 2531384624 792.2471 < 2.2e-16 ***
## s(Apps) 1 920829741 920829741 288.1920 < 2.2e-16 ***
## s(Accept) 1 111647663 111647663 34.9424 6.346e-09 ***
## s(Top10perc) 1 1058267443 1058267443 331.2058 < 2.2e-16 ***
## s(F.Undergrad) 1 273900647 273900647 85.7226 < 2.2e-16 ***
## s(Room.Board) 1 513712100 513712100 160.7764 < 2.2e-16 ***
## s(Personal) 1 54852808 54852808 17.1673 4.028e-05 ***
## s(PhD) 1 61155825 61155825 19.1399 1.483e-05 ***
## s(Terminal) 1 22024889 22024889 6.8931 0.0089216 **
## s(S.F.Ratio) 1 89787883 89787883 28.1009 1.741e-07 ***
## s(perc.alumni) 1 139078026 139078026 43.5272 1.081e-10 ***
## s(Expend) 1 468470487 468470487 146.6171 < 2.2e-16 ***
## s(Grad.Rate) 1 37708950 37708950 11.8018 0.0006417 ***
## Residuals 493 1575229760 3195196
## ---
## 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(Apps) 3 2.1126 0.097725 .
## s(Accept) 3 7.7040 4.868e-05 ***
## s(Top10perc) 3 1.9302 0.123747
## s(F.Undergrad) 3 2.0602 0.104615
## s(Room.Board) 3 1.8472 0.137676
## s(Personal) 3 2.7447 0.042544 *
## s(PhD) 3 1.9813 0.115856
## s(Terminal) 3 1.8216 0.142254
## s(S.F.Ratio) 3 4.8727 0.002386 **
## s(perc.alumni) 3 1.3279 0.264519
## s(Expend) 3 25.8833 1.443e-15 ***
## s(Grad.Rate) 3 1.3980 0.242634
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
invisible("The ANOVA test for nonparametric effects provides evidence of a non-linear relationship for Expend, Accept, S.F.Ratio, and Personal (p<0.05). All other variables are using a standard linear relationship.")
