HW6
Bella Lin
2024-04-12
6. In this exercise, you will further analyze the Wage data set considered throughout this chapter.
(a) 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.
set.seed(1)
degree_range <- 1:10
cv_errors <- rep(NA, length(degree_range))
for (i in degree_range) {
fit <- glm(wage ~ poly(age, i), data = Wage)
cv_errors[i] <- cv.glm(Wage, fit)$delta[1]
}
# Plot cross-validation errors
plot(degree_range, cv_errors, xlab = 'Degree', ylab = 'Test MSE', type = 'l')
degree_min <- which.min(cv_errors)
points(degree_min, cv_errors[degree_min], col = 'red', cex = 2, pch = 19)
The cross-validation suggested that the optimal degree d for the polynomial is 9. However, according to the plot, we can see that 4 is small enough for the polynomial, and the result from ANOVA also suggested degree 4 for polynomial. The graph below showed prediction with 9 degree model.
chosen_degree <- degree_min
fit <- glm(wage ~ poly(age, chosen_degree), data = Wage)
plot(wage ~ age, data = Wage, col = "salmon")
age_range <- range(Wage$age)
age_grid <- seq(from = age_range[1], to = age_range[2], length.out = 100)
preds <- predict(fit, newdata = list(age = age_grid))
lines(age_grid, preds, col = "purple", lwd = 2)
(b) 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.
cuts_range <- 2:10
cut_errors <- rep(NA, length(cuts_range))
for (i in cuts_range) {
Wage$age.cut <- cut(Wage$age, i)
fit <- glm(wage ~ age.cut, data = Wage)
cut_errors[i] <- cv.glm(Wage, fit)$delta[1]
}
plot(cuts_range, cut_errors[-1], xlab = 'Number of Cuts', ylab = 'Test MSE', type = 'l')
cuts_min <- which.min(cut_errors)
points(cuts_min, cut_errors[cuts_min], col = 'red', cex = 2, pch = 19)
> The cross validation suggested 8 cuts for the optimal test MSE
result, and below is the plot for 8-cuts step prediction.
plot(wage ~ age, data = Wage, col = "salmon")
fit <- glm(wage ~ cut(age, 8), data = Wage)
preds <- predict(fit, list(age = age_grid))
lines(age_grid, preds, col = "purple", lwd = 2)
10. This question relates to the College data set.
(a) 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.
set.seed(123)
train_index <- sample(1:nrow(College), nrow(College) * 0.7)
train_data <- College[train_index, ]
test_data <- College[-train_index, ]
m1 <- regsubsets(Outstate ~ ., data = train_data, nvmax = ncol(train_data) - 1, method = 'forward')
summary(m1)## Subset selection object
## Call: regsubsets.formula(Outstate ~ ., data = train_data, nvmax = ncol(train_data) -
## 1, method = "forward")
## 17 Variables (and intercept)
## Forced in Forced out
## PrivateYes FALSE FALSE
## Apps FALSE FALSE
## Accept FALSE FALSE
## Enroll FALSE FALSE
## Top10perc FALSE FALSE
## Top25perc FALSE FALSE
## F.Undergrad FALSE FALSE
## P.Undergrad FALSE FALSE
## Room.Board FALSE FALSE
## Books FALSE FALSE
## Personal FALSE FALSE
## PhD FALSE FALSE
## Terminal FALSE FALSE
## S.F.Ratio FALSE FALSE
## perc.alumni FALSE FALSE
## Expend FALSE FALSE
## Grad.Rate FALSE FALSE
## 1 subsets of each size up to 17
## Selection Algorithm: forward
## PrivateYes Apps Accept Enroll Top10perc Top25perc F.Undergrad
## 1 ( 1 ) " " " " " " " " " " " " " "
## 2 ( 1 ) "*" " " " " " " " " " " " "
## 3 ( 1 ) "*" " " " " " " " " " " " "
## 4 ( 1 ) "*" " " " " " " " " " " " "
## 5 ( 1 ) "*" " " " " " " " " " " " "
## 6 ( 1 ) "*" " " " " " " " " " " " "
## 7 ( 1 ) "*" " " "*" " " " " " " " "
## 8 ( 1 ) "*" " " "*" " " " " " " "*"
## 9 ( 1 ) "*" " " "*" " " " " "*" "*"
## 10 ( 1 ) "*" "*" "*" " " " " "*" "*"
## 11 ( 1 ) "*" "*" "*" " " " " "*" "*"
## 12 ( 1 ) "*" "*" "*" " " " " "*" "*"
## 13 ( 1 ) "*" "*" "*" " " "*" "*" "*"
## 14 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## 15 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## 16 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## 17 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## P.Undergrad Room.Board Books Personal PhD Terminal S.F.Ratio
## 1 ( 1 ) " " " " " " " " " " " " " "
## 2 ( 1 ) " " " " " " " " " " " " " "
## 3 ( 1 ) " " "*" " " " " " " " " " "
## 4 ( 1 ) " " "*" " " " " " " " " " "
## 5 ( 1 ) " " "*" " " " " "*" " " " "
## 6 ( 1 ) " " "*" " " " " "*" " " " "
## 7 ( 1 ) " " "*" " " " " "*" " " " "
## 8 ( 1 ) " " "*" " " " " "*" " " " "
## 9 ( 1 ) " " "*" " " " " "*" " " " "
## 10 ( 1 ) " " "*" " " " " "*" " " " "
## 11 ( 1 ) " " "*" " " " " "*" "*" " "
## 12 ( 1 ) " " "*" " " "*" "*" "*" " "
## 13 ( 1 ) " " "*" " " "*" "*" "*" " "
## 14 ( 1 ) " " "*" " " "*" "*" "*" " "
## 15 ( 1 ) " " "*" " " "*" "*" "*" "*"
## 16 ( 1 ) " " "*" "*" "*" "*" "*" "*"
## 17 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## perc.alumni Expend Grad.Rate
## 1 ( 1 ) " " "*" " "
## 2 ( 1 ) " " "*" " "
## 3 ( 1 ) " " "*" " "
## 4 ( 1 ) "*" "*" " "
## 5 ( 1 ) "*" "*" " "
## 6 ( 1 ) "*" "*" "*"
## 7 ( 1 ) "*" "*" "*"
## 8 ( 1 ) "*" "*" "*"
## 9 ( 1 ) "*" "*" "*"
## 10 ( 1 ) "*" "*" "*"
## 11 ( 1 ) "*" "*" "*"
## 12 ( 1 ) "*" "*" "*"
## 13 ( 1 ) "*" "*" "*"
## 14 ( 1 ) "*" "*" "*"
## 15 ( 1 ) "*" "*" "*"
## 16 ( 1 ) "*" "*" "*"
## 17 ( 1 ) "*" "*" "*"coef(m1, id = 6)## (Intercept) PrivateYes Room.Board PhD perc.alumni
## -3413.8786263 3074.9826173 0.8694753 39.2943681 48.0930260
## Expend Grad.Rate
## 0.2118468 27.7508305The final predictors selected from forward stepwise selection are: PrivateYes, Room.Board, PhD, perc.alumni, Expend, and Grad.Rate.
(b) 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.
gm <- gam(Outstate ~ Private + s(Room.Board, 5) + s(Terminal, 5) + s(perc.alumni, 5) + s(Expend, 5) + s(Grad.Rate, 5), data = train_data)
par(mfrow = c(2,3))
plot(gm, se = TRUE, col = 'purple')
> According to the plots above, Expend and Grad.Rate are strong
non-linear with Outstate.
(c) Evaluate the model obtained on the test set, and explain the results obtained.
test_pred <- predict(gm, newdata = test_data)
test_mse <- mean((test_data$Outstate - test_pred)^2)
print(test_mse)## [1] 3139704SST <- sum((test_data$Outstate - mean(test_data$Outstate))^2)
SSE <- sum((test_data$Outstate - test_pred)^2)
R_squared <- 1 - SSE / SST
print(R_squared)## [1] 0.7913189The mean squared error on the test set is 3139704 with r square value of 0.7913189, which indicates that around 79.13% of the variance in the Outstate tuition is explained by the model. This suggests a reasonably good fit of the model to the test data,
(d) For which variables, if any, is there evidence of a non-linear relationship with the response?
summary(gm)##
## Call: gam(formula = Outstate ~ Private + s(Room.Board, 5) + s(Terminal,
## 5) + s(perc.alumni, 5) + s(Expend, 5) + s(Grad.Rate, 5),
## data = train_data)
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -7430.0 -1168.7 63.9 1272.5 7366.9
##
## (Dispersion Parameter for gaussian family taken to be 3602037)
##
## Null Deviance: 9035173364 on 542 degrees of freedom
## Residual Deviance: 1858649299 on 515.9996 degrees of freedom
## AIC: 9766.95
##
## Number of Local Scoring Iterations: NA
##
## Anova for Parametric Effects
## Df Sum Sq Mean Sq F value Pr(>F)
## Private 1 2574098716 2574098716 714.623 < 2.2e-16 ***
## s(Room.Board, 5) 1 1617075141 1617075141 448.934 < 2.2e-16 ***
## s(Terminal, 5) 1 503631746 503631746 139.819 < 2.2e-16 ***
## s(perc.alumni, 5) 1 352592770 352592770 97.887 < 2.2e-16 ***
## s(Expend, 5) 1 818326211 818326211 227.184 < 2.2e-16 ***
## s(Grad.Rate, 5) 1 94541716 94541716 26.247 4.252e-07 ***
## Residuals 516 1858649299 3602037
## ---
## 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, 5) 4 3.2551 0.01186 *
## s(Terminal, 5) 4 0.5660 0.68741
## s(perc.alumni, 5) 4 0.7176 0.58013
## s(Expend, 5) 4 24.8793 < 2e-16 ***
## s(Grad.Rate, 5) 4 1.9779 0.09658 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Based on the summary output of the GAM model, Room.Board, Expend, and potentially Grad.Rate have non-linear associations with Outstate in the GAM model.