Question #6 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)
data(Wage)
summary(Wage)
## year age maritl race
## Min. :2003 Min. :18.00 1. Never Married: 648 1. White:2480
## 1st Qu.:2004 1st Qu.:33.75 2. Married :2074 2. Black: 293
## Median :2006 Median :42.00 3. Widowed : 19 3. Asian: 190
## Mean :2006 Mean :42.41 4. Divorced : 204 4. Other: 37
## 3rd Qu.:2008 3rd Qu.:51.00 5. Separated : 55
## Max. :2009 Max. :80.00
##
## education region jobclass
## 1. < HS Grad :268 2. Middle Atlantic :3000 1. Industrial :1544
## 2. HS Grad :971 1. New England : 0 2. Information:1456
## 3. Some College :650 3. East North Central: 0
## 4. College Grad :685 4. West North Central: 0
## 5. Advanced Degree:426 5. South Atlantic : 0
## 6. East South Central: 0
## (Other) : 0
## health health_ins logwage wage
## 1. <=Good : 858 1. Yes:2083 Min. :3.000 Min. : 20.09
## 2. >=Very Good:2142 2. No : 917 1st Qu.:4.447 1st Qu.: 85.38
## Median :4.653 Median :104.92
## Mean :4.654 Mean :111.70
## 3rd Qu.:4.857 3rd Qu.:128.68
## Max. :5.763 Max. :318.34
##
set.seed(234)
train <- sample(nrow(Wage), nrow(Wage)/2)
test <- (-train)
Wage.train <- Wage[train, ]
Wage.test <- Wage[test, ]
cv.error <- rep(0, 10)
for (d in 1:10) {
fit <- glm(wage ~ poly(age, d), data = Wage.train)
cv.error[d] <- cv.glm(Wage, fit, K = 10)$delta[1]
}
opt.d <- which.min(cv.error)
cat("The optimal degree is", opt.d, "\n")
## The optimal degree is
plot(wage ~ age, data = Wage, col = "red")
age.range <- range(Wage$age)
age.grid <- seq(from = age.range[1], to = age.range[2])
fit <- lm(wage ~ poly(age, 3), data = Wage)
preds <- predict(fit, newdata = list(age = age.grid))
lines(age.grid, preds, col = "blue", lwd = 2)
plot(Wage.train$age, Wage.train$wage, col = "gray", xlab = "Age", ylab = "Wage")
points(Wage.test$age, Wage.test$wage, col = "gray", pch = 20)
age.grid <- seq(from = min(Wage$age), to = max(Wage$age))
preds <- predict(fit, newdata = list(age = age.grid), se = TRUE)
se.bands <- cbind(preds$fit + 2 * preds$se.fit, preds$fit - 2 * preds$se.fit)
matlines(age.grid, se.bands, col = "blue", lty = 2)
lines(age.grid, preds$fit, lwd = 2, col = "red")
Answer: The model above models the degrees 1-10 and uses a 10 fold cross
validation. Which was done to minamize validation error and the optimal
degree ends up being 4.
data(College)
set.seed(64)
train <- sample(nrow(College), nrow(College)/4)
test <- (-train)
College.train <- College[train, ]
College.test <- College[test, ]
coef(fit, 6)
## (Intercept) poly(age, 3)1 poly(age, 3)2 poly(age, 3)3
## 111.7036 447.0679 -478.3158 125.5217
library(gam)
## Loading required package: splines
## Loading required package: foreach
## Loaded gam 1.22-3
gam.mod <- gam(Outstate ~ Private + s(Room.Board, 4) + s(Terminal, 4) + s(perc.alumni, 4) + s(Expend, 4) + s(Grad.Rate, 4), data = College, subset = train)
plot(gam.mod, se = TRUE)
c. Evaluate the model obtained on the test set, and explain the results
obtained.
pred.gam <- predict(gam.mod, College[test,])
(mse.gam <- mean((College[test,'Outstate']-pred.gam)^2))
## [1] 3555658
tss.gam <- mean((College[test,'Outstate'] - mean(College[test,'Outstate']))^2)
(rss.test = 1 - mse.gam/tss.gam)
## [1] 0.7839986
Answer: The Rsquared is 0.78
summary(gam.mod)
##
## Call: gam(formula = Outstate ~ Private + s(Room.Board, 4) + s(Terminal,
## 4) + s(perc.alumni, 4) + s(Expend, 4) + s(Grad.Rate, 4),
## data = College, subset = train)
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -6434.57 -1092.63 40.86 1222.85 4507.48
##
## (Dispersion Parameter for gaussian family taken to be 3798928)
##
## Null Deviance: 2937155412 on 193 degrees of freedom
## Residual Deviance: 653415234 on 171.9999 degrees of freedom
## AIC: 3512.342
##
## Number of Local Scoring Iterations: NA
##
## Anova for Parametric Effects
## Df Sum Sq Mean Sq F value Pr(>F)
## Private 1 871209037 871209037 229.330 < 2.2e-16 ***
## s(Room.Board, 4) 1 577121262 577121262 151.917 < 2.2e-16 ***
## s(Terminal, 4) 1 213625879 213625879 56.233 3.276e-12 ***
## s(perc.alumni, 4) 1 119201995 119201995 31.378 8.279e-08 ***
## s(Expend, 4) 1 217698197 217698197 57.305 2.174e-12 ***
## s(Grad.Rate, 4) 1 38321298 38321298 10.087 0.00177 **
## Residuals 172 653415234 3798928
## ---
## 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, 4) 3 1.4409 0.2326004
## s(Terminal, 4) 3 0.8592 0.4635226
## s(perc.alumni, 4) 3 0.7155 0.5439959
## s(Expend, 4) 3 7.3895 0.0001098 ***
## s(Grad.Rate, 4) 3 0.7534 0.5217772
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Conclusion: The GAM model suggests that several predictors have a significant impact on the outcome variable Outstate, both in linear and non-linear terms. However, the model also indicates unexplained variability, as shown by the dispersion parameter, and the model’s complexity is reflected in the AIC value.