Assignment 6
Valeria De La Rocha
today
library(ISLR2)
library(boot)
library(leaps)
library(gam)
set.seed(123)Problem 6
data("Wage")
head(Wage)6(a) Polynomial Regression
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.
# Create a vector to store the cross-validation errors
cv_errors <- rep(NA, 10)
# Test polynomial degrees from 1 through 10
for (degree in 1:10) {
polynomial_model <- glm(
wage ~ poly(age, degree),
data = Wage
)
cv_result <- cv.glm(
Wage,
polynomial_model,
K = 10
)
cv_errors[degree] <- cv_result$delta[1]
}
# Display the cross-validation errors
polynomial_results <- data.frame(
Degree = 1:10,
CV_Error = cv_errors
)
polynomial_results# Select the degree with the smallest cross-validation error
best_degree <- which.min(cv_errors)
best_degree## [1] 10The polynomial degree selected by cross-validation is:
## Degree 10ANOVA
fit_degree_1 <- lm(wage ~ age, data = Wage)
fit_degree_2 <- lm(wage ~ poly(age, 2), data = Wage)
fit_degree_3 <- lm(wage ~ poly(age, 3), data = Wage)
fit_degree_4 <- lm(wage ~ poly(age, 4), data = Wage)
fit_degree_5 <- lm(wage ~ poly(age, 5), data = Wage)
anova(
fit_degree_1,
fit_degree_2,
fit_degree_3,
fit_degree_4,
fit_degree_5
)Plot of the Polynomial Fit
# Fit the model using the degree selected by cross-validation
best_polynomial_model <- lm(
wage ~ poly(age, best_degree),
data = Wage
)
# Create a sequence of age values for the prediction line
age_grid <- seq(
from = min(Wage$age),
to = max(Wage$age),
length.out = 200
)
# Predict wage for the age grid
polynomial_predictions <- predict(
best_polynomial_model,
newdata = data.frame(age = age_grid),
interval = "confidence"
)
# Plot the original observations
plot(
Wage$age,
Wage$wage,
col = "lightgray",
pch = 16,
xlab = "Age",
ylab = "Wage",
main = paste("Polynomial Regression: Degree", best_degree)
)
# Add the fitted polynomial line
lines(
age_grid,
polynomial_predictions[, "fit"],
lwd = 3
)
# Add confidence interval lines
lines(
age_grid,
polynomial_predictions[, "lwr"],
lty = 2
)
lines(
age_grid,
polynomial_predictions[, "upr"],
lty = 2
)
6(b) Step Function
Fit a step function to predict wage using age, and perform crossvalidation to choose the optimal number of cuts. Make a plot of the fit obtained.
# Test models with 2 through 10 age groups
number_of_groups <- 2:10
step_cv_errors <- rep(NA, length(number_of_groups))
for (i in seq_along(number_of_groups)) {
groups <- number_of_groups[i]
# Create age groups using equally spaced breaks
age_breaks <- seq(
min(Wage$age),
max(Wage$age),
length.out = groups + 1
)
Wage$age_group <- cut(
Wage$age,
breaks = age_breaks,
include.lowest = TRUE
)
step_model <- glm(
wage ~ age_group,
data = Wage
)
step_cv <- cv.glm(
Wage,
step_model,
K = 10
)
step_cv_errors[i] <- step_cv$delta[1]
}
step_results <- data.frame(
Number_of_Groups = number_of_groups,
CV_Error = step_cv_errors
)
step_resultsbest_groups <- number_of_groups[which.min(step_cv_errors)]
best_groups## [1] 8The step function selected by cross-validation uses:
## 8 age groupsPlot of the Step Function
# Create the best set of age groups
best_breaks <- seq(
min(Wage$age),
max(Wage$age),
length.out = best_groups + 1
)
Wage$best_age_group <- cut(
Wage$age,
breaks = best_breaks,
include.lowest = TRUE
)
# Fit the final step-function model
best_step_model <- lm(
wage ~ best_age_group,
data = Wage
)
# Create prediction data
step_age_grid <- seq(
min(Wage$age),
max(Wage$age),
length.out = 200
)
step_prediction_data <- data.frame(
age = step_age_grid
)
step_prediction_data$best_age_group <- cut(
step_prediction_data$age,
breaks = best_breaks,
include.lowest = TRUE
)
step_predictions <- predict(
best_step_model,
newdata = step_prediction_data
)
# Plot observations
plot(
Wage$age,
Wage$wage,
col = "lightgray",
pch = 16,
xlab = "Age",
ylab = "Wage",
main = paste("Step Function with", best_groups, "Age Groups")
)
# Add the step-function prediction
lines(
step_age_grid,
step_predictions,
type = "s",
lwd = 3
)
Problem 10
data("College")
head(College)10(a) Training and Test Sets and Forward Stepwise Selection
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)
training_rows <- sample(
1:nrow(College),
size = round(0.70 * nrow(College))
)
college_train <- College[training_rows, ]
college_test <- College[-training_rows, ]
nrow(college_train)## [1] 544nrow(college_test)## [1] 233# Model with no predictors
empty_model <- lm(
Outstate ~ 1,
data = college_train
)
# Model with all available predictors
full_model <- lm(
Outstate ~ .,
data = college_train
)
# Forward stepwise selection using AIC
forward_model <- step(
empty_model,
scope = list(
lower = formula(empty_model),
upper = formula(full_model)
),
direction = "forward",
trace = 0
)
summary(forward_model)##
## Call:
## lm(formula = Outstate ~ Expend + Private + Room.Board + perc.alumni +
## PhD + Grad.Rate + Top25perc + Personal + Accept + F.Undergrad +
## Apps + Terminal + Top10perc + Enroll, data = college_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6498.6 -1327.5 -1.8 1287.5 10012.9
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.362e+03 7.027e+02 -3.361 0.000832 ***
## Expend 1.968e-01 2.492e-02 7.899 1.64e-14 ***
## PrivateYes 2.739e+03 3.075e+02 8.907 < 2e-16 ***
## Room.Board 7.492e-01 1.031e-01 7.268 1.32e-12 ***
## perc.alumni 4.247e+01 9.307e+00 4.564 6.26e-06 ***
## PhD 1.440e+01 1.077e+01 1.337 0.181662
## Grad.Rate 1.890e+01 6.830e+00 2.767 0.005848 **
## Top25perc 2.888e+00 1.031e+01 0.280 0.779458
## Personal -2.592e-01 1.436e-01 -1.805 0.071608 .
## Accept 8.419e-01 1.567e-01 5.372 1.17e-07 ***
## F.Undergrad -1.221e-01 7.667e-02 -1.593 0.111826
## Apps -2.324e-01 8.801e-02 -2.640 0.008529 **
## Terminal 2.327e+01 1.159e+01 2.008 0.045162 *
## Top10perc 2.613e+01 1.313e+01 1.990 0.047138 *
## Enroll -6.679e-01 4.407e-01 -1.516 0.130183
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2019 on 529 degrees of freedom
## Multiple R-squared: 0.7613, Adjusted R-squared: 0.755
## F-statistic: 120.5 on 14 and 529 DF, p-value: < 2.2e-16The predictors selected by forward stepwise selection are:
selected_predictors <- attr(
terms(forward_model),
"term.labels"
)
selected_predictors## [1] "Expend" "Private" "Room.Board" "perc.alumni" "PhD"
## [6] "Grad.Rate" "Top25perc" "Personal" "Accept" "F.Undergrad"
## [11] "Apps" "Terminal" "Top10perc" "Enroll"The forward selection procedure chooses a smaller group of predictors instead of using every variable in the data set. This creates a model that is easier to interpret.
10(b) Generalized Additive Model
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.
# Check whether each selected variable is numeric or categorical
numeric_selected <- selected_predictors[
sapply(
college_train[selected_predictors],
is.numeric
)
]
factor_selected <- selected_predictors[
!sapply(
college_train[selected_predictors],
is.numeric
)
]
numeric_selected## [1] "Expend" "Room.Board" "perc.alumni" "PhD" "Grad.Rate"
## [6] "Top25perc" "Personal" "Accept" "F.Undergrad" "Apps"
## [11] "Terminal" "Top10perc" "Enroll"factor_selected## [1] "Private"# Create smooth terms for numeric predictors
smooth_terms <- paste0(
"s(",
numeric_selected,
", df = 4)"
)
# Keep factor predictors as regular terms
gam_terms <- c(
smooth_terms,
factor_selected
)
# Create the GAM formula
gam_formula <- as.formula(
paste(
"Outstate ~",
paste(gam_terms, collapse = " + ")
)
)
gam_formula## Outstate ~ s(Expend, df = 4) + s(Room.Board, df = 4) + s(perc.alumni,
## df = 4) + s(PhD, df = 4) + s(Grad.Rate, df = 4) + s(Top25perc,
## df = 4) + s(Personal, df = 4) + s(Accept, df = 4) + s(F.Undergrad,
## df = 4) + s(Apps, df = 4) + s(Terminal, df = 4) + s(Top10perc,
## df = 4) + s(Enroll, df = 4) + Private# Fit the GAM
gam_model <- gam(
gam_formula,
data = college_train
)
summary(gam_model)##
## Call: gam(formula = gam_formula, data = college_train)
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -6192.47 -1088.43 90.76 1183.24 7314.86
##
## (Dispersion Parameter for gaussian family taken to be 3389677)
##
## Null Deviance: 9035479394 on 543 degrees of freedom
## Residual Deviance: 1660944196 on 490.0007 degrees of freedom
## AIC: 9776.65
##
## Number of Local Scoring Iterations: NA
##
## Anova for Parametric Effects
## Df Sum Sq Mean Sq F value Pr(>F)
## s(Expend, df = 4) 1 3902947905 3902947905 1151.4217 < 2.2e-16 ***
## s(Room.Board, df = 4) 1 775586461 775586461 228.8083 < 2.2e-16 ***
## s(perc.alumni, df = 4) 1 515550905 515550905 152.0944 < 2.2e-16 ***
## s(PhD, df = 4) 1 17176366 17176366 5.0673 0.024824 *
## s(Grad.Rate, df = 4) 1 165303525 165303525 48.7667 9.433e-12 ***
## s(Top25perc, df = 4) 1 920588 920588 0.2716 0.602505
## s(Personal, df = 4) 1 63117577 63117577 18.6205 1.930e-05 ***
## s(Accept, df = 4) 1 38797560 38797560 11.4458 0.000774 ***
## s(F.Undergrad, df = 4) 1 229343143 229343143 67.6593 1.762e-15 ***
## s(Apps, df = 4) 1 24726149 24726149 7.2945 0.007156 **
## s(Terminal, df = 4) 1 243005 243005 0.0717 0.789005
## s(Top10perc, df = 4) 1 23173664 23173664 6.8365 0.009206 **
## s(Enroll, df = 4) 1 7581960 7581960 2.2368 0.135405
## Private 1 263474936 263474936 77.7286 < 2.2e-16 ***
## Residuals 490 1660944196 3389677
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Anova for Nonparametric Effects
## Npar Df Npar F Pr(F)
## (Intercept)
## s(Expend, df = 4) 3 24.2166 1.255e-14 ***
## s(Room.Board, df = 4) 3 1.5865 0.19175
## s(perc.alumni, df = 4) 3 0.6311 0.59518
## s(PhD, df = 4) 3 2.2005 0.08718 .
## s(Grad.Rate, df = 4) 3 2.8843 0.03533 *
## s(Top25perc, df = 4) 3 1.6570 0.17541
## s(Personal, df = 4) 3 2.5274 0.05675 .
## s(Accept, df = 4) 3 11.5677 2.452e-07 ***
## s(F.Undergrad, df = 4) 3 2.8680 0.03611 *
## s(Apps, df = 4) 3 0.5370 0.65713
## s(Terminal, df = 4) 3 0.6671 0.57258
## s(Top10perc, df = 4) 3 0.7914 0.49911
## s(Enroll, df = 4) 3 1.6684 0.17290
## Private
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1GAM Plots
# Plot the smooth effects
par(mfrow = c(2, 2))
plot(
gam_model,
se = TRUE
)
par(mfrow = c(1, 1))
10(c) Test Set Evaluation
Evaluate the model obtained on the test set, and explain the results obtained.
gam_test_predictions <- predict(
gam_model,
newdata = college_test
)
# Calculate test mean squared error
gam_test_mse <- mean(
(college_test$Outstate - gam_test_predictions)^2
)
# Calculate test root mean squared error
gam_test_rmse <- sqrt(gam_test_mse)
gam_test_mse## [1] 3000587gam_test_rmse## [1] 1732.22The test mean squared error is:
## [1] 3000587The test root mean squared error is:
## [1] 1732.2210(d) Evidence of Non-Linear Relationships
For which variables, if any, is there evidence of a non-linear relationship with the response?
summary(gam_model)##
## Call: gam(formula = gam_formula, data = college_train)
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -6192.47 -1088.43 90.76 1183.24 7314.86
##
## (Dispersion Parameter for gaussian family taken to be 3389677)
##
## Null Deviance: 9035479394 on 543 degrees of freedom
## Residual Deviance: 1660944196 on 490.0007 degrees of freedom
## AIC: 9776.65
##
## Number of Local Scoring Iterations: NA
##
## Anova for Parametric Effects
## Df Sum Sq Mean Sq F value Pr(>F)
## s(Expend, df = 4) 1 3902947905 3902947905 1151.4217 < 2.2e-16 ***
## s(Room.Board, df = 4) 1 775586461 775586461 228.8083 < 2.2e-16 ***
## s(perc.alumni, df = 4) 1 515550905 515550905 152.0944 < 2.2e-16 ***
## s(PhD, df = 4) 1 17176366 17176366 5.0673 0.024824 *
## s(Grad.Rate, df = 4) 1 165303525 165303525 48.7667 9.433e-12 ***
## s(Top25perc, df = 4) 1 920588 920588 0.2716 0.602505
## s(Personal, df = 4) 1 63117577 63117577 18.6205 1.930e-05 ***
## s(Accept, df = 4) 1 38797560 38797560 11.4458 0.000774 ***
## s(F.Undergrad, df = 4) 1 229343143 229343143 67.6593 1.762e-15 ***
## s(Apps, df = 4) 1 24726149 24726149 7.2945 0.007156 **
## s(Terminal, df = 4) 1 243005 243005 0.0717 0.789005
## s(Top10perc, df = 4) 1 23173664 23173664 6.8365 0.009206 **
## s(Enroll, df = 4) 1 7581960 7581960 2.2368 0.135405
## Private 1 263474936 263474936 77.7286 < 2.2e-16 ***
## Residuals 490 1660944196 3389677
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Anova for Nonparametric Effects
## Npar Df Npar F Pr(F)
## (Intercept)
## s(Expend, df = 4) 3 24.2166 1.255e-14 ***
## s(Room.Board, df = 4) 3 1.5865 0.19175
## s(perc.alumni, df = 4) 3 0.6311 0.59518
## s(PhD, df = 4) 3 2.2005 0.08718 .
## s(Grad.Rate, df = 4) 3 2.8843 0.03533 *
## s(Top25perc, df = 4) 3 1.6570 0.17541
## s(Personal, df = 4) 3 2.5274 0.05675 .
## s(Accept, df = 4) 3 11.5677 2.452e-07 ***
## s(F.Undergrad, df = 4) 3 2.8680 0.03611 *
## s(Apps, df = 4) 3 0.5370 0.65713
## s(Terminal, df = 4) 3 0.6671 0.57258
## s(Top10perc, df = 4) 3 0.7914 0.49911
## s(Enroll, df = 4) 3 1.6684 0.17290
## Private
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1