Assignment #6
Lisa Smith
2022-11-05
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.
library(ISLR2)## Warning: package 'ISLR2' was built under R version 4.1.3attach(Wage)library(boot)
set.seed(1)
degree <- 10
cr.va.err = rep(NA,10)
for (i in 1:10) {
glm.fit = glm(wage~poly(age, i), data=Wage)
cr.va.err[i] = cv.glm(Wage, glm.fit, K=10)$delta[2]
}Plot:
plot(1:10, cr.va.err, xlab="Degree", ylab="Cross Validation Error", type = "l", pch=20, lwd=2)
min.deg = min(cr.va.err)
std.dev.points = sd(cr.va.err)
abline(h=min.deg + 0.2 * std.dev.points, col="blue", lty="dashed")
abline(h=min.deg - 0.2 * std.dev.points, col="blue", lty="dashed") 
fit.1 = lm(wage~poly(age, 1), 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)
fit.6 = lm(wage~poly(age, 6), data=Wage)
fit.7 = lm(wage~poly(age, 7), data=Wage)
fit.8 = lm(wage~poly(age, 8), data=Wage)
fit.9 = lm(wage~poly(age, 9), data=Wage)
fit.10 = lm(wage~poly(age, 10), data=Wage)
anova(fit.1, fit.2, fit.3, fit.4, fit.5, fit.6, fit.7, fit.8, fit.9, fit.10)## Analysis of Variance Table
##
## Model 1: wage ~ poly(age, 1)
## Model 2: wage ~ poly(age, 2)
## Model 3: wage ~ poly(age, 3)
## Model 4: wage ~ poly(age, 4)
## Model 5: wage ~ poly(age, 5)
## Model 6: wage ~ poly(age, 6)
## Model 7: wage ~ poly(age, 7)
## Model 8: wage ~ poly(age, 8)
## Model 9: wage ~ poly(age, 9)
## Model 10: wage ~ poly(age, 10)
## 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 ' ' 1Solution: Our Cross Validation plot showing the standard deviation suggests that d=3 is the degree with the smallest cv error. Our ANOVA shows that the polynomials with degree 2 and 3 are significant at the alpha = 0.01 level. Polynomials over 3 degrees are not significant at alpha = 0.01.
Make a plot of the resulting polynomial fit to the data:
plot(wage~age, data=Wage, col="darkgrey")
age.range = range(Wage$age)
age.grid = seq(from=age.range[1], to=age.range[2])
lm.fit = lm(wage~poly(age, 3), data=Wage)
lm.pred = predict(lm.fit, data.frame(age=age.grid))
lines(age.grid, lm.pred, col="blue", lwd=2)
- 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.
cr.va.err = rep(NA, 10)
for (i in 2:10) {
Wage$age.cut = cut(Wage$age, i)
lm.fit = glm(wage~age.cut, data=Wage)
cr.va.err[i] = cv.glm(Wage, lm.fit, K=10)$delta[2]
}
plot(2:10, cr.va.err[-1], xlab="# of cuts", ylab="Cross-Validation Error", type="l", pch=20, lwd=2)
We can see on the plot above that the Cross-Validation Error is the lowest for 8 cuts. Now, we apply 8 cuts to the entire data:
lm.fit = glm(wage~cut(age, 8), data=Wage)
age.range = range(Wage$age)
age.grid = seq(from=age.range[1], to=age.range[2])
lm.pred = predict(lm.fit, data.frame(age=age.grid))
plot(wage~age, data=Wage, col="darkgrey")
lines(age.grid, lm.pred, col="blue", lwd=2)
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.
library(leaps)## Warning: package 'leaps' was built under R version 4.1.3attach(College)Col.train <- sample(1: nrow(College), nrow(College)/2)
Col.test <- -Col.train
reg.fit <- regsubsets(Outstate ~ ., data = College, subset = Col.train, method = 'forward')
reg.summary <- summary(reg.fit)
reg.summary## Subset selection object
## Call: regsubsets.formula(Outstate ~ ., data = College, subset = Col.train,
## 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 8
## 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 ) "*" " " " " " " " " "*" " "
## 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 ) " " "*" " " "*" "*" " " " "
## perc.alumni Expend Grad.Rate
## 1 ( 1 ) " " " " " "
## 2 ( 1 ) "*" " " " "
## 3 ( 1 ) "*" "*" " "
## 4 ( 1 ) "*" "*" " "
## 5 ( 1 ) "*" "*" " "
## 6 ( 1 ) "*" "*" "*"
## 7 ( 1 ) "*" "*" "*"
## 8 ( 1 ) "*" "*" "*"The scores show that 6 is the minimum size for the subset. We select 6 as the best subset size and now apply it to the entire data.
reg.fit = regsubsets(Outstate ~ ., data = College, method = "forward")
coef = coef(reg.fit, id = 6)
names(coef)## [1] "(Intercept)" "PrivateYes" "Room.Board" "PhD" "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(gam)## Warning: package 'gam' was built under R version 4.1.3## Loading required package: splines## Loading required package: foreach## Warning: package 'foreach' was built under R version 4.1.3## Loaded gam 1.20.2gam.fit = gam(Outstate ~ Private + s(Room.Board, 5) + s(Terminal, 5) + s(perc.alumni, 5) + s(Expend, 5) + s(Grad.Rate, 5), data = College, subset = Col.train)
par(mfrow = c(2, 3))
plot(gam.fit, se = T, col = "blue")
- Evaluate the model obtained on the test set, and explain the results obtained.
gam.pred <- predict(gam.fit, College[Col.test , ])
RSS <- sum((College[Col.test, ]$Outstate - gam.pred)^2)
TSS <- sum((College[Col.test, ]$Outstate - mean(College[Col.test, ]$Outstate)) ^ 2)
1 - (RSS / TSS)## [1] 0.7649038Solution: The R-squared on the College test set is: 0.76
- 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, 5) + s(Terminal,
## 5) + s(perc.alumni, 5) + s(Expend, 5) + s(Grad.Rate, 5),
## data = College, subset = Col.train)
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -7289.5 -1004.3 18.3 1123.6 4218.8
##
## (Dispersion Parameter for gaussian family taken to be 3138798)
##
## Null Deviance: 6139053909 on 387 degrees of freedom
## Residual Deviance: 1133105994 on 361 degrees of freedom
## AIC: 6933.339
##
## Number of Local Scoring Iterations: NA
##
## Anova for Parametric Effects
## Df Sum Sq Mean Sq F value Pr(>F)
## Private 1 1658551575 1658551575 528.404 < 2.2e-16 ***
## s(Room.Board, 5) 1 1093958629 1093958629 348.528 < 2.2e-16 ***
## s(Terminal, 5) 1 239592419 239592419 76.332 < 2.2e-16 ***
## s(perc.alumni, 5) 1 189302589 189302589 60.310 8.461e-14 ***
## s(Expend, 5) 1 671008681 671008681 213.779 < 2.2e-16 ***
## s(Grad.Rate, 5) 1 87504239 87504239 27.878 2.236e-07 ***
## Residuals 361 1133105994 3138798
## ---
## 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.6201 0.006576 **
## s(Terminal, 5) 4 2.3018 0.058243 .
## s(perc.alumni, 5) 4 0.8690 0.482600
## s(Expend, 5) 4 28.0768 < 2.2e-16 ***
## s(Grad.Rate, 5) 4 2.7848 0.026556 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Solution: Our results in the non-parametric ANOVA indicate evidence of a non-linear relationshio between “Expend” and the response.