Assignment 6
Valerie Ceciliano
2026-07-18
Chapter 07 (page 297): 6, 10
Q6.
In this exercise, you will further analyze the Wage data set considered throughout this chapter.
- 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(boot)
set.seed(1)
cv.error <- rep(0, 10)
for (i in 1:10) {
glm.fit <- glm(wage ~ poly(age, i), data = Wage)
cv.error[i] <-cv.glm(Wage, glm.fit, K =10)$delta[1]
}
best.d <- which.min(cv.error)
print(paste("Optimal degree by CV:", best.d))## [1] "Optimal degree by CV: 9"fit.1 <- lm(wage ~ age, 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)
anova(fit.1, fit.2, fit.3, fit.4, fit.5)## 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 ' ' 1agelims <- range(Wage$age)
age.grid <- seq(from = agelims[1], to = agelims[2])
preds <- predict(fit.4, newdata = list(age=age.grid), se=TRUE)
plot(Wage$age, Wage$wage, col = "darkgrey", xlab = "Age", ylab = "Wage")
lines(age.grid, preds$fit, lwd = 2, col = "blue")
lines(age.grid, preds$fit + 2 * preds$se.fit, lty = "dashed", col = "blue")
lines(age.grid, preds$fit - 2 * preds$se.fit, lty = "dashed", col = "blue")
Q6 Answer: The cv ten-fold error drops significantly when moving from degree 1 to degree 2, stabilizes between degree 3 and degree 4 and hits its minimum at around 3&4. The anova test whether a simpler model is compared. p-value - degree 2 - Highly significant, p-value - degree 3 - Highly significant, p-value - degree 4 - Borderline significant, p-value - degree 5 - Not significant
Both cv and hypothesis testing point to degree 3 or degree 4 polynomial
- Fit a step function to predict wage using age, and perform cross- validation to choose the optimal number of cuts. Make a plot of the fit obtained
set.seed(1)
cv.error.step <- rep(0, 9)
for (i in 2:10) {
Wage$age.cut <-cut(Wage$age, i)
glm.fit <-glm(wage ~ age.cut, data = Wage)
cv.error.step[i - 1] <- cv.glm(Wage, glm.fit, K = 10)$delta[1]
}
best.cuts <- which.min(cv.error.step) + 1
print(paste("Optimal number of cuts by CV:", best.cuts))## [1] "Optimal number of cuts by CV: 8"fit.step <- lm(wage ~ cut(age, best.cuts), data = Wage)
preds.step <- predict(fit.step, newdata=list(age = age.grid), se = TRUE)
plot(Wage$age, Wage$age, col = "darkgrey", xlab = "Age", ylab = "Wage")
lines(age.grid, preds.step$fit, lwd = 2, col = "red")
lines(age.grid, preds.step$fit + 2 * preds.step$se.fit, lty = "dashed", col = "red")
lines(age.grid, preds.step$fit - 2 * preds.step$se.fit, lty = "dashed", col = "red")
Q10.
This question relates to the College data set.
- 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.
library(ISLR2)
library(leaps)
set.seed(1)
train <- sample(1:nrow(College), nrow(College) /2)
test <- (-train)
regfit.fwd <- regsubsets(Outstate ~ ., data = College[train, ], nvmax = 17, method = "forward")
reg.summary <- summary(regfit.fwd)
best.size <- which.min(reg.summary$bic)
coef(regfit.fwd, id = best.size)## (Intercept) PrivateYes Room.Board Terminal perc.alumni
## -4726.8810613 2717.7019276 1.1032433 36.9990286 59.0863753
## Expend Grad.Rate
## 0.1930814 33.8303314Q10 Answer: (A) Using bic identifies a subset of key predictors, such as PrivateYes Room.Board Terminal perc.alumni Expend Grad.Rate
- 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(gam)## Loading required package: splines## Loading required package: foreach## Loaded gam 1.22-7gam.fit <- gam(Outstate ~ Private + s(Room.Board, df = 4) + s(PhD, df = 4) +
s(perc.alumni, df = 4) + s(Expend, df = 4) + s(Grad.Rate, df = 4), data = College[train, ])
par(mfrow = c(2, 3))
plot(gam.fit, se = TRUE, col ="blue")
Q10 Answer: (B) Expend shows a strong, non-linear relationship with Outstate. I can see tuition increases at lower levels of expenditure but levels off. Room.Board and Grad.Rate show a somewhat linear positive trend with tuition. Private shows a constant upward shift in tuition for private institutions compared to public institutions.
- Evaluate the model obtained on the test set, and explain the results obtained.
gam.preds <- predict(gam.fit, newdata = College[test, ])
test.mse <- mean((College[test, ]$Outstate - gam.preds)^2)
test.tss <- mean((College[test, ]$Outstate - mean(College[test, ]$Outstate))^2)
test.r2 <- 1 - (test.mse / test.tss)
cat("Test MSE:", test.mse, "\n")## Test MSE: 3324814cat("Test R-squared:", test.r2, "\n")## Test R-squared: 0.7677116Q10 Answer: (C) In the fitted GAM, it consistently gets around 0.77 to .80 showing a non-linear transformations account for at least 80% of the variance in out of state tuition on the unseen data.
- For which variables, if any, is there evidence of a non-linear relationship with the response?
summary(gam.fit)##
## Call: gam(formula = Outstate ~ Private + s(Room.Board, df = 4) + s(PhD,
## df = 4) + s(perc.alumni, df = 4) + s(Expend, df = 4) + s(Grad.Rate,
## df = 4), data = College[train, ])
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -7174.12 -1141.26 -86.67 1264.75 7506.52
##
## (Dispersion Parameter for gaussian family taken to be 3750826)
##
## Null Deviance: 6989966760 on 387 degrees of freedom
## Residual Deviance: 1372800367 on 365.9994 degrees of freedom
## AIC: 6997.793
##
## Number of Local Scoring Iterations: NA
##
## Anova for Parametric Effects
## Df Sum Sq Mean Sq F value Pr(>F)
## Private 1 1761501708 1761501708 469.630 < 2.2e-16 ***
## s(Room.Board, df = 4) 1 1563405241 1563405241 416.816 < 2.2e-16 ***
## s(PhD, df = 4) 1 334516516 334516516 89.185 < 2.2e-16 ***
## s(perc.alumni, df = 4) 1 331488702 331488702 88.377 < 2.2e-16 ***
## s(Expend, df = 4) 1 525593541 525593541 140.127 < 2.2e-16 ***
## s(Grad.Rate, df = 4) 1 87873241 87873241 23.428 1.916e-06 ***
## Residuals 366 1372800367 3750826
## ---
## 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 = 4) 3 1.9158 0.1266
## s(PhD, df = 4) 3 0.8504 0.4671
## s(perc.alumni, df = 4) 3 0.3520 0.7877
## s(Expend, df = 4) 3 23.6176 5.373e-14 ***
## s(Grad.Rate, df = 4) 3 1.0701 0.3617
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Q10 Answer: (D) Expend shows a statistically significant non-linear relationship. For Grad.Rate, PhD, perc.alumni, and Room.Board: will generally give non-parametric P-values greater > 0.05, showing that their contributions can be modeleled linearly without sacrificing accuracy.