Assignment 6
Jason Villa
2026-07-28
Problem 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.
library(ISLR2)## Warning: package 'ISLR2' was built under R version 4.3.3library(boot)
data(Wage)
head(Wage)## year age maritl race education region
## 231655 2006 18 1. Never Married 1. White 1. < HS Grad 2. Middle Atlantic
## 86582 2004 24 1. Never Married 1. White 4. College Grad 2. Middle Atlantic
## 161300 2003 45 2. Married 1. White 3. Some College 2. Middle Atlantic
## 155159 2003 43 2. Married 3. Asian 4. College Grad 2. Middle Atlantic
## 11443 2005 50 4. Divorced 1. White 2. HS Grad 2. Middle Atlantic
## 376662 2008 54 2. Married 1. White 4. College Grad 2. Middle Atlantic
## jobclass health health_ins logwage wage
## 231655 1. Industrial 1. <=Good 2. No 4.318063 75.04315
## 86582 2. Information 2. >=Very Good 2. No 4.255273 70.47602
## 161300 1. Industrial 1. <=Good 1. Yes 4.875061 130.98218
## 155159 2. Information 2. >=Very Good 1. Yes 5.041393 154.68529
## 11443 2. Information 1. <=Good 1. Yes 4.318063 75.04315
## 376662 2. Information 2. >=Very Good 1. Yes 4.845098 127.11574set.seed(42)
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]
}
degree <- which.min(cv.error)
degree## [1] 9bestfit1 <- lm(wage ~ age, data = Wage)
bestfit2 <- lm(wage ~ poly(age, 2), data = Wage)
bestfit3 <- lm(wage ~ poly(age, 3), data = Wage)
bestfit4 <- lm(wage ~ poly(age, 4), data = Wage)
anova(bestfit1, bestfit2, bestfit3, bestfit4)## 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)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2998 5022216
## 2 2997 4793430 1 228786 143.6025 < 2.2e-16 ***
## 3 2996 4777674 1 15756 9.8894 0.001679 **
## 4 2995 4771604 1 6070 3.8101 0.051039 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(Wage$age, Wage$wage, col = "darkgrey", xlab = "Age", ylab = "Wage")
- 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(42)
cv.error.b <- numeric(10)
for (i in 2:10) {
Wage$age.cut <- cut(Wage$age, i)
fit <- glm(wage ~ age.cut, data = Wage)
cv.error[i] <- cv.glm(Wage, fit, K = 10)$delta[1]
}
print(cv.error.b[2:10])## [1] 0 0 0 0 0 0 0 0 0fit.step <- lm(wage ~ cut(age, 4), data = Wage)
summary(fit.step)##
## Call:
## lm(formula = wage ~ cut(age, 4), data = Wage)
##
## Residuals:
## Min 1Q Median 3Q Max
## -98.126 -24.803 -6.177 16.493 200.519
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 94.158 1.476 63.790 <2e-16 ***
## cut(age, 4)(33.5,49] 24.053 1.829 13.148 <2e-16 ***
## cut(age, 4)(49,64.5] 23.665 2.068 11.443 <2e-16 ***
## cut(age, 4)(64.5,80.1] 7.641 4.987 1.532 0.126
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 40.42 on 2996 degrees of freedom
## Multiple R-squared: 0.0625, Adjusted R-squared: 0.06156
## F-statistic: 66.58 on 3 and 2996 DF, p-value: < 2.2e-16age.grid <- seq(from = min(Wage$age), to = max(Wage$age), by = 1)
preds1 <- predict(fit.step, newdata = data.frame(age = age.grid), se = TRUE)
plot(Wage$age, Wage$wage, col = "darkgrey", xlab = "Age", ylab = "Wage ($1000s)")
lines(age.grid, preds1$fit, col = "red", lwd = 3)
lines(age.grid, preds1$fit + 2 * preds1$se.fit, col = "red", lty = 2)
lines(age.grid, preds1$fit - 2 * preds1$se.fit, col = "red", lty = 2)
Problem 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.
library(leaps)## Warning: package 'leaps' was built under R version 4.3.3set.seed(1)
attach(College)
train <- sample(length(Outstate), length(Outstate) / 2)
test <- -train
College.train <- College[train, ]
College.test <- College[test, ]
fit <- regsubsets(Outstate ~ ., data = College.train, nvmax = 17, method = "forward")
fit.summary <- summary(fit)
par(mfrow = c(1, 3))
plot(fit.summary$cp, xlab = "Number of variables", ylab = "Cp", type = "l")
min.cp <- min(fit.summary$cp)
std.cp <- sd(fit.summary$cp)
abline(h = min.cp + 0.2 * std.cp, col = "red", lty = 2)
abline(h = min.cp - 0.2 * std.cp, col = "red", lty = 2)
plot(fit.summary$bic, xlab = "Number of variables", ylab = "BIC", type='l')
min.bic <- min(fit.summary$bic)
std.bic <- sd(fit.summary$bic)
abline(h = min.bic + 0.2 * std.bic, col = "red", lty = 2)
abline(h = min.bic - 0.2 * std.bic, col = "red", lty = 2)
plot(fit.summary$adjr2, xlab = "Number of variables", ylab = "Adjusted R2", type = "l", ylim = c(0.4, 0.84))
max.adjr2 <- max(fit.summary$adjr2)
std.adjr2 <- sd(fit.summary$adjr2)
abline(h = max.adjr2 + 0.2 * std.adjr2, col = "red", lty = 2)
abline(h = max.adjr2 - 0.2 * std.adjr2, col = "red", lty = 2)
#Cp, BIC and adjr2 show that size 6 is the minimum size for the subset for which the scores are within 0.2 standard devitations of optimum.
lm1 <- regsubsets(Outstate ~ ., data = College, method = "forward")
coeffs <- coef(fit, id = 6)
names(coeffs)## [1] "(Intercept)" "PrivateYes" "Room.Board" "Terminal" "perc.alumni"
## [6] "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(mgcv)## Warning: package 'mgcv' was built under R version 4.3.3## Loading required package: nlme## Warning: package 'nlme' was built under R version 4.3.3## This is mgcv 1.9-1. For overview type 'help("mgcv-package")'.gam1 <- gam(Outstate ~ Private + s(Room.Board, k = 3) + s(PhD, k = 3) + s(perc.alumni, k = 3) + s(Expend, k = 6) + s(Grad.Rate, k = 3), data=College.train)
par(mfrow = c(2, 3))
plot(gam1, se = T, col = "blue")
- Evaluate the model obtained on the test set, and explain the results obtained.
preds <- predict(gam1, College.test)
err <- mean((College.test$Outstate - preds)^2)
err## [1] 3412194tss <- mean((College.test$Outstate - mean(College.test$Outstate))^2)
rss <- 1 - err / tss
rss## [1] 0.7616068We obtain a test R^2 of 0.77 using GAM with 6 predictors.
- For which variables, if any, is there evidence of a non-linear relationship with the response?
summary(gam1)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Outstate ~ Private + s(Room.Board, k = 3) + s(PhD, k = 3) + s(perc.alumni,
## k = 3) + s(Expend, k = 6) + s(Grad.Rate, k = 3)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8740.0 225.6 38.737 < 2e-16 ***
## PrivateYes 2408.1 280.9 8.571 2.67e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Room.Board) 1.416 1.658 38.604 < 2e-16 ***
## s(PhD) 1.000 1.000 1.837 0.176
## s(perc.alumni) 1.000 1.000 21.085 6.43e-06 ***
## s(Expend) 4.829 4.985 32.819 < 2e-16 ***
## s(Grad.Rate) 1.000 1.000 24.237 1.95e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.796 Deviance explained = 80.1%
## GCV = 3.7993e+06 Scale est. = 3.6892e+06 n = 388ANOVA shows a strong evidence of non-linear relationship between “Outstate” and “Expend”“, and a moderately strong non-linear relationship (using p-value of 0.05) between”Outstate” and “Grad.Rate”” or “PhD”.