Assignment #6
2024-04-12
QUESTION 6
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.
set.seed(20)
k_folds <- 10
cv_errors <- vapply(
1:10,
function(degree) {
# Fit a generalized linear model with polynomial terms
model <- glm(wage ~ poly(age, degree), data = Wage)
# Calculate the cross-validation error
cv.glm(Wage, model, K = k_folds)$delta[1]
},
numeric(1)
)
cv_error_std_error <- sd(cv_errors) / sqrt(k_folds)
ggplot(data = data.frame(Degree = 1:10, CV_Error = cv_errors), aes(x = Degree, y = CV_Error)) +
geom_line() +
geom_hline(yintercept = min(cv_errors) + cv_error_std_error, color = "red", linetype = "dashed") +
geom_point(aes(x = which.min(cv_errors), y = min(cv_errors)), color = "red", size = 3) +
ylab("CV Error") +
xlab("Degree") +
scale_x_continuous(breaks = 1:10) +
theme_minimal() 
polynomial_models <- lapply(1:10, function(degree) {
lm(wage ~ poly(age, degree), data = Wage)
})
model_comparison <- do.call(anova, polynomial_models)
model_comparison## Analysis of Variance Table
##
## Model 1: wage ~ poly(age, degree)
## Model 2: wage ~ poly(age, degree)
## Model 3: wage ~ poly(age, degree)
## Model 4: wage ~ poly(age, degree)
## Model 5: wage ~ poly(age, degree)
## Model 6: wage ~ poly(age, degree)
## Model 7: wage ~ poly(age, degree)
## Model 8: wage ~ poly(age, degree)
## Model 9: wage ~ poly(age, degree)
## Model 10: wage ~ poly(age, degree)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2998 5022216
## 2 2997 4793430 1 228786 143.7638 < 2.2e-16 ***
## 3 2996 4777674 1 15756 9.9005 0.001669 **
## 4 2995 4771604 1 6070 3.8143 0.050909 .
## 5 2994 4770322 1 1283 0.8059 0.369398
## 6 2993 4766389 1 3932 2.4709 0.116074
## 7 2992 4763834 1 2555 1.6057 0.205199
## 8 2991 4763707 1 127 0.0796 0.777865
## 9 2990 4756703 1 7004 4.4014 0.035994 *
## 10 2989 4756701 1 3 0.0017 0.967529
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ninth degree had the lowest error, but third degree was chosen in accordance with the one standard error rule. ANOVA results support this because the P-value was just above 0.05 at the third degree.
fit <- lapply(1:10, function(degree) {
lm(wage ~ poly(age, degree), data = Wage)
})
selected_model <- fit[[3]]
age_range <- range(Wage$age)
age_grid <- seq(from = age_range[1], to = age_range[2], length.out = 100)
predictions <- predict(selected_model, newdata = list(age = age_grid), se = TRUE)
predictions_df <- data.frame(
Age = age_grid,
Predicted_Wage = predictions$fit,
Lower_CI = predictions$fit - 2 * predictions$se.fit,
Upper_CI = predictions$fit + 2 * predictions$se.fit
)
ggplot() +
geom_point(data = Wage, aes(x = age, y = wage), alpha = 0.3) +
geom_line(data = predictions_df, aes(x = Age, y = Predicted_Wage), colour = "blue") +
geom_ribbon(data = predictions_df, aes(x = Age, ymin = Lower_CI, ymax = Upper_CI),
fill = "blue", alpha = 0.2) +
theme_minimal() +
labs(x = "Age", y = "Wage", title = "Fitted Model with Confidence Intervals")
- 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(20)
cuts <- 1:10
k_folds <- 10
cv_errors <- sapply(cuts, function(num_breaks) {
Wage$age_cut <- cut(Wage$age, breaks = num_breaks + 1, include.lowest = TRUE)
model <- glm(wage ~ age_cut, data = Wage)
cv.glm(Wage, model, K = k_folds)$delta[1]
})
cv_stderr <- sd(cv_errors) / sqrt(k_folds)
ggplot(data = data.frame(Cuts = cuts, CV_Error = cv_errors), aes(x = Cuts, y = CV_Error)) +
geom_line() +
geom_hline(yintercept = min(cv_errors) + cv_stderr, color = "red", linetype = "dashed") +
geom_point(aes(x = which.min(cv_errors), y = min(cv_errors)), color = "red", size = 3) +
labs(x = "Breaks", y = "CV Error",
title = "Effect of Categorizing Age on Model Error") +
scale_x_continuous(breaks = cuts) +
theme_minimal() 
fit <- lm(wage ~ cut(age, 7), data = Wage)
preds <- predict(fit, newdata = list(age = age_grid), se = TRUE)
preds_df <- data.frame(
predictors = age_grid,
predictions = preds$fit,
lower = preds$fit - 2 * preds$se.fit,
upper = preds$fit + 2 * preds$se.fit)
ggplot() +
geom_point(data = Wage, aes(x = age, y = wage), alpha = 0.3) +
geom_line(data = preds_df, aes(x = predictors, y = predictions), colour = "blue", size = 1) +
geom_ribbon(data = preds_df, aes(ymin = lower, ymax = upper, x = predictors), fill = "blue", alpha = 0.2) +
theme_minimal() +
labs(x = "Age", y = "Wage", title = "Fitted Model with Confidence Intervals")## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
QUESTION 10
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.
set.seed(1)
train_indices <- sample(nrow(College), nrow(College) / 2)
test_indices <- setdiff(1:nrow(College), train_indices)
College_train <- College[train_indices, ]
College_test <- College[test_indices, ]
fit <- regsubsets(Outstate ~ ., data = College_train, nvmax = 17, method = "forward")
fit_summary <- summary(fit)
par(mfrow = c(1, 3))
plot_criteria <- function(statistics, label, higher_is_better = FALSE) {
plot(statistics, xlab = "Number of variables", ylab = label, type = "l")
optimal_value <- if(higher_is_better) max(statistics) else min(statistics)
stat_std <- sd(statistics)
abline(h = optimal_value + 0.2 * stat_std, col = "red", lty = 2)
abline(h = optimal_value - 0.2 * stat_std, col = "red", lty = 2)
}
plot_criteria(fit_summary$cp, "Cp")
plot_criteria(fit_summary$bic, "BIC")
plot_criteria(fit_summary$adjr2, "Adjusted R2", higher_is_better = TRUE)
coeffs <- coef(fit, id = 10)
names(coeffs)## [1] "(Intercept)" "PrivateYes" "Top10perc" "Room.Board" "Books"
## [6] "Personal" "Terminal" "S.F.Ratio" "perc.alumni" "Expend"
## [11] "Grad.Rate"10 Variable model containing “Private” “Top10perc” “Room.Board” “Books” “Personal” “Terminal” “S.F.Ratio” “perc.alumni” “Expend” “Grad.Rate” is suggested by forward selection.
- 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.
fit <- gam(Outstate ~ Private + s(Top10perc) + s(Room.Board) + s(Books) + s(Personal) + s(Terminal) + s(S.F.Ratio) + s(perc.alumni) + s(Expend) + s(Grad.Rate), data=College_train)
par(mfrow = c(3, 4))
plot(fit, se = T, col = "blue")
Some non-linear relationship between Outstate and Books and Expend (presence of curved lines). Some linear relationship between out-of-state tuition and Room.Board, perc.alumni (presence of straight diagonal lines).
- Evaluate the model obtained on the test set, and explain the results obtained.
preds <- predict(fit, newdata = College_test)
mse_nonlinear <- mean((preds - College_test$Outstate)^2)
fit_lm <- lm(Outstate ~ Private + Top10perc + Room.Board + Books + Personal + Terminal + S.F.Ratio + perc.alumni + Expend + Grad.Rate, data = College_train)
preds_linear <- predict(fit_lm, newdata = College_test)
mse_linear <- mean((preds_linear - College_test$Outstate)^2)
paste("Non-linear RMSE:", sqrt(mse_nonlinear))## [1] "Non-linear RMSE: 1865.17865088501"paste("Linear RMSE:", sqrt(mse_linear))## [1] "Linear RMSE: 1965.64187456604"Non-linear model is performing better on the college_test dataset because its RMSE of 1865.18 is lower than the linear model’s RMSE of 1965.64.
- For which variables, if any, is there evidence of a non-linear relationship with the response?
summary(fit)##
## Call: gam(formula = Outstate ~ Private + s(Top10perc) + s(Room.Board) +
## s(Books) + s(Personal) + s(Terminal) + s(S.F.Ratio) + s(perc.alumni) +
## s(Expend) + s(Grad.Rate), data = College_train)
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -5873.54 -1055.69 -57.15 1231.86 7083.04
##
## (Dispersion Parameter for gaussian family taken to be 3606488)
##
## Null Deviance: 6989966760 on 387 degrees of freedom
## Residual Deviance: 1262269084 on 349.9995 degrees of freedom
## AIC: 6997.224
##
## Number of Local Scoring Iterations: NA
##
## Anova for Parametric Effects
## Df Sum Sq Mean Sq F value Pr(>F)
## Private 1 1671815897 1671815897 463.558 < 2.2e-16 ***
## s(Top10perc) 1 1738332759 1738332759 482.002 < 2.2e-16 ***
## s(Room.Board) 1 709692318 709692318 196.782 < 2.2e-16 ***
## s(Books) 1 55141092 55141092 15.289 0.0001108 ***
## s(Personal) 1 38773902 38773902 10.751 0.0011467 **
## s(Terminal) 1 45671307 45671307 12.664 0.0004243 ***
## s(S.F.Ratio) 1 66142478 66142478 18.340 2.390e-05 ***
## s(perc.alumni) 1 79848520 79848520 22.140 3.655e-06 ***
## s(Expend) 1 463112678 463112678 128.411 < 2.2e-16 ***
## s(Grad.Rate) 1 61349309 61349309 17.011 4.647e-05 ***
## Residuals 350 1262269084 3606488
## ---
## 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(Top10perc) 3 0.6531 0.58148
## s(Room.Board) 3 1.7642 0.15370
## s(Books) 3 1.7482 0.15683
## s(Personal) 3 1.0351 0.37705
## s(Terminal) 3 1.4086 0.24002
## s(S.F.Ratio) 3 2.2068 0.08702 .
## s(perc.alumni) 3 0.3706 0.77429
## s(Expend) 3 20.1035 4.811e-12 ***
## s(Grad.Rate) 3 0.5756 0.63138
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1P-values for Anova Nonparametric effects indicate that Expend has a non-linear relationship with Outstate. However, as opposed to previous plots, Books does not seem to have a non-linear relationship.