Assignment #6
Alliah Kelley
2022-05-11
Question 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(ISLR)
library(boot)set.seed(1)
Q6delta=rep(NA, 10)
for (i in 1:10) {
Q6fit=glm(wage ~ poly(age, i), data=Wage)
Q6delta[i]=cv.glm(Wage, Q6fit, K = 10)$delta[1]
}
plot(1:10, Q6delta, xlab = "Optimal Degree", ylab = "Cross Val", type = "l")
Q6min=which.min(Q6delta)
points(which.min(Q6delta), Q6delta[which.min(Q6delta)], cex = 2, pch = 20)
This shows that the optimal degree d = 9 using Cross Validation.
Q6fit.1=lm(wage~poly(age, 1), data=Wage)
Q6fit.2=lm(wage~poly(age, 2), data=Wage)
Q6fit.3=lm(wage~poly(age, 3), data=Wage)
Q6fit.4=lm(wage~poly(age, 4), data=Wage)
Q6fit.5=lm(wage~poly(age, 5), data=Wage)
Q6fit.6=lm(wage~poly(age, 6), data=Wage)
Q6fit.7=lm(wage~poly(age, 7), data=Wage)
Q6fit.8=lm(wage~poly(age, 8), data=Wage)
Q6fit.9=lm(wage~poly(age, 9), data=Wage)
Q6fit.10=lm(wage~poly(age, 10), data=Wage)
anova(Q6fit.1, Q6fit.2, Q6fit.3, Q6fit.4, Q6fit.5, Q6fit.6, Q6fit.7, Q6fit.8, Q6fit.9, Q6fit.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 ' ' 1Using ANOVA shows that almost all polynomials above degree 3 show a higher significance.
plot(wage~age, data=Wage, col="black")
Q6range=range(Wage$age)
Q6age=seq(from=Q6range[1], to=Q6range[2])
Q6fit=lm(wage~poly(age, 3), data=Wage)
Q6pred=predict(Q6fit, data.frame(age=Q6age))
lines(Q6age, Q6pred, col="red", lwd=2)
(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.
#K=10
Crossval=rep(NA, 10)
for (i in 2:10) {
Wage$age.cut=cut(Wage$age, i)
CVfit=glm(wage ~ age.cut, data = Wage)
Crossval[i]=cv.glm(Wage, CVfit, K = 10)$delta[1]
}
plot(2:10, Crossval[-1], xlab = "Optimal Cuts", ylab = "Cross Val", type = "l")
Q6bmin=which.min(Crossval)
points(which.min(Crossval), Crossval[which.min(Crossval)], col = "blue", cex = 2, pch = 20)
This shows that the optimal number of cuts is 8.
plot(wage ~ age, data = Wage, col = "black")
Q6range=range(Wage$age)
Q6age=seq(from = Q6range[1], to = Q6range[2])
Q6fit=glm(wage ~ cut(age, 8), data = Wage)
Q6pred=predict(Q6fit, data.frame(age = Q6age))
lines(Q6age, Q6pred, 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(ISLR)
library(leaps)
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 Adjusted R2 seem to show that a size of 6 is the minimum size with a 0.2 standard deviation.
Q10fit=regsubsets(Outstate ~ ., data = College, method = "forward")
Q10coeffs=coef(fit, id = 6)
names(Q10coeffs)## [1] "(Intercept)" "PrivateYes" "Room.Board" "PhD" "perc.alumni"
## [6] "Expend" "Grad.Rate"(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 the predictors. Plot the results, and explain your findings.
library(gam)Q10gamfit=gam(Outstate ~ Private + s(Room.Board, df = 2) + s(PhD, df = 2) +
s(perc.alumni, df = 2) + s(Expend, df = 5) + s(Grad.Rate, df = 2), data = College.train)
par(mfrow = c(2, 3))
plot(Q10gamfit, se = T, col = "red")
(c) Evaluate the model obtained on the test set, and explain the results obtained.
Q10gampred = predict(Q10gamfit, College.test)
Q10gamerr = mean((College.test$Outstate - Q10gampred)^2)
Q10gamerr## [1] 3756005Q10gamtss = mean((College.test$Outstate - mean(College.test$Outstate))^2)
Q10testrss = 1 - Q10gamerr/Q10gamtss
Q10testrss## [1] 0.7722687This showed a test R-squared of 0.77 using GAM and is a slight improvement over test RSS of 0.74.
(d) For which variables, if any, is there evidence of a non-linear relationship with the response?
summary(Q10gamfit)##
## Call: gam(formula = Outstate ~ Private + s(Room.Board, df = 2) + s(PhD,
## df = 2) + s(perc.alumni, df = 2) + s(Expend, df = 5) + s(Grad.Rate,
## df = 2), data = College.train)
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -7758.60 -1070.79 32.98 1220.39 4232.54
##
## (Dispersion Parameter for gaussian family taken to be 3275358)
##
## Null Deviance: 6139053909 on 387 degrees of freedom
## Residual Deviance: 1221709840 on 373.0004 degrees of freedom
## AIC: 6938.55
##
## Number of Local Scoring Iterations: NA
##
## Anova for Parametric Effects
## Df Sum Sq Mean Sq F value Pr(>F)
## Private 1 1676778467 1676778467 511.938 < 2.2e-16 ***
## s(Room.Board, df = 2) 1 1077471252 1077471252 328.963 < 2.2e-16 ***
## s(PhD, df = 2) 1 314826907 314826907 96.120 < 2.2e-16 ***
## s(perc.alumni, df = 2) 1 183001813 183001813 55.872 5.565e-13 ***
## s(Expend, df = 5) 1 579739180 579739180 177.000 < 2.2e-16 ***
## s(Grad.Rate, df = 2) 1 80903699 80903699 24.701 1.022e-06 ***
## Residuals 373 1221709840 3275358
## ---
## 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 4.5861 0.03288 *
## s(PhD, df = 2) 1 0.8275 0.36357
## s(perc.alumni, df = 2) 1 0.5100 0.47554
## s(Expend, df = 5) 4 26.2683 < 2e-16 ***
## s(Grad.Rate, df = 2) 1 2.0718 0.15088
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1This shows a significant relation between Outstate and
Expend, and semi-significance between Outstate
and Grad.Rate or PhD.