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(ISLR)
library(boot)
set.seed(1)
degree = 10
cv.errs = rep(NA, degree)
for (i in 1:degree) {
  fit = glm(wage ~ poly(age, i), data = Wage)
  cv.errs[i] = cv.glm(Wage, fit)$delta[1]
}
plot(1:degree, cv.errs, xlab = 'Degree', ylab = 'Test MSE', type = 'l')
deg.min = which.min(cv.errs)
points(deg.min, cv.errs[deg.min], col = 'blue', cex = 2, pch = 19)

fit1 <- lm(wage ~ age, data = Wage)
fit2 <- lm(wage ~ poly(age, 2), data = Wage)
fit3 <- lm(wage ~ poly(age, 3), data = Wage)
fit4 <- lm(wage ~ poly(age, 4), data = Wage)
fit5 <- lm(wage ~ poly(age, 5), data = Wage)
anova(fit1, fit2, fit3, fit4, fit5)
## 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)
## Model 5: wage ~ poly(age, 5)
##   Res.Df     RSS Df Sum of Sq        F    Pr(>F)    
## 1   2998 5022216                                    
## 2   2997 4793430  1    228786 143.5931 < 2.2e-16 ***
## 3   2996 4777674  1     15756   9.8888  0.001679 ** 
## 4   2995 4771604  1      6070   3.8098  0.051046 .  
## 5   2994 4770322  1      1283   0.8050  0.369682    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The minimum of test MSE is at the degree 9. But test MSE of degree 4 is also small and with comparison to ANOVA degree 4 should be used.

plot(wage ~ age, data = Wage, col = "lightblue")
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)

(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.

cv.errs = rep(NA, degree)
for (i in 2:degree) {
  Wage$age.cut = cut(Wage$age, i)
  fit = glm(wage ~ age.cut, data = Wage)
  cv.errs[i] = cv.glm(Wage, fit)$delta[1]
}
plot(2:degree, cv.errs[-1], xlab = 'Cuts', ylab = 'Test MSE', type = 'l')
deg.min = which.min(cv.errs)
points(deg.min, cv.errs[deg.min], col = 'blue', cex = 2, pch = 19)

Minimum Test MSE = 8 cuts.

plot(wage ~ age, data = Wage, col = "lightblue")
fit <- glm(wage ~ cut(age, 8), data = Wage)
preds <- predict(fit, list(age = age.grid))
lines(age.grid, preds, col = "blue", lwd = 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 step-wise selection on the training set in order to identify a satisfactory model that uses just a subset of the predictors.

library(leaps)
library(ISLR)
set.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 = "blue", lty = 2)
abline(h = min.cp - 0.2 * std.cp, col = "blue", 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 = "blue", lty = 2)
abline(h = min.bic - 0.2 * std.bic, col = "blue", 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 = "blue", lty = 2)
abline(h = max.adjr2 - 0.2 * std.adjr2, col = "blue", 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.

fit = regsubsets(Outstate ~ ., data = College, method = "forward")
coeffs = coef(fit, id = 6)
names(coeffs)
## [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)

gam.fit = gam(Outstate ~ Private + s(Room.Board, df = 2) + s(Terminal, 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(gam.fit, se=TRUE, col="red")

(c) Evaluate the model obtained on the test set, and explain the results obtained.

gam.pred = predict(gam.fit, College.test)
gam.err = mean((College.test$Outstate - gam.pred)^2)
gam.err
## [1] 3378340
tss = mean((College.test$Outstate - mean(College.test$Outstate))^2)
rss = 1 - gam.err / tss
rss
## [1] 0.763972

We obtain a test R^2 of 0.77 using GAM with 6 predictors.

(d) 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, df = 2) + s(Terminal, 
##     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 
## -7318.32 -1123.66   -28.31  1259.62  7338.91 
## 
## (Dispersion Parameter for gaussian family taken to be 3701607)
## 
##     Null Deviance: 6989966760 on 387 degrees of freedom
## Residual Deviance: 1380699954 on 373.0002 degrees of freedom
## AIC: 6986.018 
## 
## Number of Local Scoring Iterations: 2 
## 
## Anova for Parametric Effects
##                         Df     Sum Sq    Mean Sq F value    Pr(>F)    
## Private                  1 1782594667 1782594667 481.573 < 2.2e-16 ***
## s(Room.Board, df = 2)    1 1604874368 1604874368 433.562 < 2.2e-16 ***
## s(Terminal, df = 2)      1  289618280  289618280  78.241 < 2.2e-16 ***
## s(perc.alumni, df = 2)   1  349447569  349447569  94.404 < 2.2e-16 ***
## s(Expend, df = 5)        1  578389738  578389738 156.254 < 2.2e-16 ***
## s(Grad.Rate, df = 2)     1   90976435   90976435  24.578 1.086e-06 ***
## Residuals              373 1380699954    3701607                      
## ---
## 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.7567    0.1858    
## s(Terminal, df = 2)          1  1.2035    0.2733    
## s(perc.alumni, df = 2)       1  0.1715    0.6790    
## s(Expend, df = 5)            4 21.6541 4.441e-16 ***
## s(Grad.Rate, df = 2)         1  0.4668    0.4948    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ANOVA shows a strong non-linear relationship between “Outstate” and “Expend”, Using p-value of 0.05, ANOVA shows a moderately strong non-linear relationship between ”Outstate" , “Grad.Rate”" and “PhD”.