Questions 6, 10

0.0.1 Analyzing the Wage Data Set

  1. In this exercise, you will further analyze the Wage data set considered throughout this chapter.
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6      ✔ purrr   0.3.5 
## ✔ tibble  3.1.8      ✔ dplyr   1.0.10
## ✔ tidyr   1.2.1      ✔ stringr 1.4.1 
## ✔ readr   2.1.3      ✔ forcats 0.5.2 
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(ISLR2)
library(boot)
library(leaps)
library(gam)
## Loading required package: splines
## Loading required package: foreach
## 
## Attaching package: 'foreach'
## 
## The following objects are masked from 'package:purrr':
## 
##     accumulate, when
## 
## Loaded gam 1.20.2
str(Wage)
## 'data.frame':    3000 obs. of  11 variables:
##  $ year      : int  2006 2004 2003 2003 2005 2008 2009 2008 2006 2004 ...
##  $ age       : int  18 24 45 43 50 54 44 30 41 52 ...
##  $ maritl    : Factor w/ 5 levels "1. Never Married",..: 1 1 2 2 4 2 2 1 1 2 ...
##  $ race      : Factor w/ 4 levels "1. White","2. Black",..: 1 1 1 3 1 1 4 3 2 1 ...
##  $ education : Factor w/ 5 levels "1. < HS Grad",..: 1 4 3 4 2 4 3 3 3 2 ...
##  $ region    : Factor w/ 9 levels "1. New England",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ jobclass  : Factor w/ 2 levels "1. Industrial",..: 1 2 1 2 2 2 1 2 2 2 ...
##  $ health    : Factor w/ 2 levels "1. <=Good","2. >=Very Good": 1 2 1 2 1 2 2 1 2 2 ...
##  $ health_ins: Factor w/ 2 levels "1. Yes","2. No": 2 2 1 1 1 1 1 1 1 1 ...
##  $ logwage   : num  4.32 4.26 4.88 5.04 4.32 ...
##  $ wage      : num  75 70.5 131 154.7 75 ...
attach(Wage)
  1. 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.
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 = 'red',
       cex = 2,
       pch = 19)

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

  1. 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 = 'red',
       cex = 2,
       pch = 19)

The plot shows 8 as the minimum test MSE.

plot(wage ~ age,
     data = Wage,
     col = "darkgrey")
fit = glm(wage ~ cut(age, 8),
          data = Wage)
preds = predict(fit, list(age = age.grid))
lines(age.grid, preds,
      col = "red",
      lwd = 2)

0.0.2 Analyzing the College Data Set

  1. This question relates to the College data set.
  1. 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.
train = sample(1: nrow(College),
               nrow(College)/2)
test = -train
fit = regsubsets(Outstate ~ .,
                 data = College,
                 subset = train,
                 method = 'forward')
fit.summary = summary(fit)
fit.summary
## Subset selection object
## Call: regsubsets.formula(Outstate ~ ., data = College, subset = 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 ) "*"         "*"    "*"
coef(fit, id = 6)
##   (Intercept)    PrivateYes    Room.Board           PhD   perc.alumni 
## -3815.6574509  2880.3858979     0.9861841    43.6735045    40.4602197 
##        Expend     Grad.Rate 
##     0.1770944    30.8363935
  1. 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.
gam.mod = gam(Outstate ~ Private + s(Room.Board, 5) + s(Terminal, 5) + s(perc.alumni, 5) + s(Expend, 5) + s(Grad.Rate, 5),
              data = College,
              subset = train)
par(mfrow = c(2,3))
plot(gam.mod,
     se = TRUE,
     col = "blue")

  1. Evaluate the model obtained on the test set, and explain the results obtained.

The R2 statistic on test set is .78.

preds = predict(gam.mod,
                College[test, ])
RSS = sum((College[test, ]$Outstate - preds)^2)
TSS = sum((College[test, ]$Outstate - mean(College[test, ]$Outstate)) ^ 2)
1 - (RSS / TSS)
## [1] 0.7649038
  1. For which variables, if any, is there evidence of a non-linear relationship with the response?

ANOVA for Nonparametric Effects shows Expend has strong non-linear relationship with Outstate.

summary(gam.mod)
## 
## 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 = 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 ' ' 1