Assignment#6
Ashley Windham
7/19/2021
Chapter 7: 6, 10
Question 6
In this exercise, you will furthur analyze the Wage dataset 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, how does this compare to the results of hypothesis testing using ANOVA? Make a plot of the resulting polynomial fit to the data.
library(ISLR)## Warning: package 'ISLR' was built under R version 4.0.5library(boot)## Warning: package 'boot' was built under R version 4.0.5set.seed(0)
data(Wage)
attach(Wage)
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[2]
}
plot(1:10, cv.error, xlab="Degree of Polynomial", ylab="CV error", type="b", pch=20, lwd=2, ylim=c(1590, 1700))
min.point = min(cv.error)
sd.points = sd(cv.error)
abline(h=min.point + 0.2 * sd.points, col="red", lty="dashed")
abline(h=min.point - 0.2 * sd.points, col="red", lty="dashed")
legend("topright", "0.2-standard deviation lines", lty="dashed", col="red")
m0 <- lm(wage ~ 1, data = Wage)
m1 <- lm(wage ~ poly(age, 1), data = Wage)
m2 <- lm(wage ~ poly(age, 2), data = Wage)
m3 <- lm(wage ~ poly(age, 3), data = Wage)
m4 <- lm(wage ~ poly(age, 4), data = Wage)
m5 <- lm(wage ~ poly(age, 5), data = Wage)
m6 <- lm(wage ~ poly(age, 6), data = Wage)
m7 <- lm(wage ~ poly(age, 7), data = Wage)
m8 <- lm(wage ~ poly(age, 8), data = Wage)
m9 <- lm(wage ~ poly(age, 9), data = Wage)
m10 <- lm(wage ~ poly(age, 10), data = Wage)
anova(m0, m1, m2, m3, m4, m5, m6, m7, m8, m9, m10)## Analysis of Variance Table
##
## Model 1: wage ~ 1
## Model 2: wage ~ poly(age, 1)
## Model 3: wage ~ poly(age, 2)
## Model 4: wage ~ poly(age, 3)
## Model 5: wage ~ poly(age, 4)
## Model 6: wage ~ poly(age, 5)
## Model 7: wage ~ poly(age, 6)
## Model 8: wage ~ poly(age, 7)
## Model 9: wage ~ poly(age, 8)
## Model 10: wage ~ poly(age, 9)
## Model 11: wage ~ poly(age, 10)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2999 5222086
## 2 2998 5022216 1 199870 125.5934 < 2.2e-16 ***
## 3 2997 4793430 1 228786 143.7638 < 2.2e-16 ***
## 4 2996 4777674 1 15756 9.9005 0.001669 **
## 5 2995 4771604 1 6070 3.8143 0.050909 .
## 6 2994 4770322 1 1283 0.8059 0.369398
## 7 2993 4766389 1 3932 2.4709 0.116074
## 8 2992 4763834 1 2555 1.6057 0.205199
## 9 2991 4763707 1 127 0.0796 0.777865
## 10 2990 4756703 1 7004 4.4014 0.035994 *
## 11 2989 4756701 1 3 0.0017 0.967529
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(wage~age, data=Wage, col="darkgrey")
agelims = range(Wage$age)
age.grid = seq(from=agelims[1], to=agelims[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)
- ANOVA shows that the polynomials above 4 are insigificant with p-values above 0.05.
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.
for (i in 2:10) {
Wage$age.cut = cut(Wage$age, i)
lm.fit = glm(wage~age.cut, data=Wage)
cv.error[i] = cv.glm(Wage, lm.fit, K=10)$delta[2]
}
plot(2:10, cv.error[-1], xlab="Number of cuts", ylab="CV error", type="l", pch=20, lwd=2)
lm.fit = glm(wage~cut(age, 8), data=Wage)
agelims = range(Wage$age)
age.grid = seq(from=agelims[1], to=agelims[2])
lm.pred = predict(lm.fit, data.frame(age=age.grid))
plot(wage~age, data=Wage, col="darkgrey")
lines(age.grid, lm.pred, col="red", lwd=2)
Question 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.0.5library(ISLR)
data(College)
set.seed(1)
train = sample(length(College$Outstate), length(College$Outstate)/2)
test = -train
college.train = College[train,]
college.test = College[test, ]
p = ncol(College) - 1
reg.fit = regsubsets(Outstate~., data = college.train, nvmax = p, method = "forward")
reg.summary = summary(reg.fit)
par(mfrow = c(1, 3))
plot(reg.summary$cp, xlab = "Number of Variables", ylab = "Cp", type = "l")
min.cp = min(reg.summary$cp)
std.cp = sd(reg.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(reg.summary$bic, xlab = "Number of Variables", ylab = "BIC", type = "l")
min.bic = min(reg.summary$bic)
std.bic = sd(reg.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(reg.summary$adjr2, xlab = "Number of Variables", ylab = "Adjusted R2",
type = "l", ylim = c(0.4, 0.84))
max.adjr2 = max(reg.summary$adjr2)
std.adjr2 = sd(reg.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)
which.min(reg.summary$cp)## [1] 14which.min(reg.summary$bic)## [1] 6which.max(reg.summary$adjr2)## [1] 14reg.fit = regsubsets(Outstate ~., data = College, method ="forward")
coefi = coef(reg.fit, id=6)
names(coefi)## [1] "(Intercept)" "PrivateYes" "Room.Board" "PhD" "perc.alumni"
## [6] "Expend" "Grad.Rate"- For the lowest BIC, the number of variables is 6. Additionally, 6 variables is the lowest number of variables between the 0.2 standard deviations for the optimal CP and adjusted R2.
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 predictors. Plot the results, and explain your findings.
library(gam)## Warning: package 'gam' was built under R version 4.0.5## Loading required package: splines## Loading required package: foreach## Warning: package 'foreach' was built under R version 4.0.5## Loaded gam 1.20fit = gam(Outstate ~ Private + s(Room.Board, df=2) + s(Personal, df=2) + s(PhD, df=2) + s(Terminal, df=2) + s(perc.alumni, df=2) + s(Expend, df=2) + s(Grad.Rate, df=2), data = college.train)
par(mfrow = c(2, 3))
plot(fit, se = T, col = "blue")
- The R squared is 0.77 using GAM with 6 predictors.
c) Evaluate the model obtained on the test set, and explain the results obtained.
gam.pred = predict(fit, college.test)
gam.err = mean((college.test$Outstate - gam.pred)^2)
gam.err## [1] 3419856gam.tss = mean((college.test$Outstate - mean(college.test$Outstate))^2)
test.rss = 1-gam.err/gam.tss
test.rss## [1] 0.7610715d)For which variables, if any, is there evidence of a non-linear relationship with the response?
summary(fit)##
## Call: gam(formula = Outstate ~ Private + s(Room.Board, df = 2) + s(Personal,
## df = 2) + s(PhD, df = 2) + s(Terminal, df = 2) + s(perc.alumni,
## df = 2) + s(Expend, df = 2) + s(Grad.Rate, df = 2), data = college.train)
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -6287.87 -1297.72 -45.08 1311.20 8263.91
##
## (Dispersion Parameter for gaussian family taken to be 3876712)
##
## Null Deviance: 6989966760 on 387 degrees of freedom
## Residual Deviance: 1442137582 on 372.0002 degrees of freedom
## AIC: 7004.91
##
## Number of Local Scoring Iterations: NA
##
## Anova for Parametric Effects
## Df Sum Sq Mean Sq F value Pr(>F)
## Private 1 1838502725 1838502725 474.2428 < 2.2e-16 ***
## s(Room.Board, df = 2) 1 1696989483 1696989483 437.7393 < 2.2e-16 ***
## s(Personal, df = 2) 1 86755162 86755162 22.3785 3.184e-06 ***
## s(PhD, df = 2) 1 378603646 378603646 97.6610 < 2.2e-16 ***
## s(Terminal, df = 2) 1 24898611 24898611 6.4226 0.01168 *
## s(perc.alumni, df = 2) 1 291408295 291408295 75.1689 < 2.2e-16 ***
## s(Expend, df = 2) 1 442273794 442273794 114.0848 < 2.2e-16 ***
## s(Grad.Rate, df = 2) 1 63695744 63695744 16.4304 6.149e-05 ***
## Residuals 372 1442137582 3876712
## ---
## 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, df = 2) 1 1.989 0.15930
## s(Personal, df = 2) 1 1.825 0.17753
## s(PhD, df = 2) 1 3.844 0.05068 .
## s(Terminal, df = 2) 1 0.858 0.35483
## s(perc.alumni, df = 2) 1 0.379 0.53837
## s(Expend, df = 2) 1 47.723 2.139e-11 ***
## s(Grad.Rate, df = 2) 1 0.880 0.34890
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- The Expend variable is showing a strong non-linear relationship with the response. Additionally, there is a moderate non-linear relationship between the PhD variable and the response.