Assignment #6

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(ISLR2)
library(boot)

attach(Wage)

# cross validation
set.seed(1)
cv.error.10 <- rep(0, 10)
for (i in 1:10) {
  glm.fit <- glm(wage ~ poly(age, i), data = Wage)
  cv.error.10[i] <- cv.glm(Wage, glm.fit, K = 10)$delta[1]
}
cv.error.10

 [1] 1676.826 1600.763 1598.399 1595.651
 [5] 1594.977 1596.061 1594.298 1598.134
 [9] 1593.913 1595.950

plot(1:10, cv.error.10, type = "b", pch = 19,
     xlab = "Degree of Polynomial", ylab = "CV Error")

d <- which.min(cv.error.10)
d

[1] 9

# ANOVA testing
fit.1 <- lm(wage ~ age, data = Wage)
fit.2 <- lm(wage ~ poly(age, 2), data = Wage)
fit.3 <- lm(wage ~ poly(age, 3), data = Wage)
fit.4 <- lm(wage ~ poly(age, 4), data = Wage)
fit.5 <- lm(wage ~ poly(age, 5), data = Wage)

anova(fit.1, fit.2, fit.3, fit.4, fit.5)

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
1   2998 5022216                      
2   2997 4793430  1    228786 143.5931
3   2996 4777674  1     15756   9.8888
4   2995 4771604  1      6070   3.8098
5   2994 4770322  1      1283   0.8050
     Pr(>F)    
1              
2 < 2.2e-16 ***
3  0.001679 ** 
4  0.051046 .  
5  0.369682    
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Fitting the cross validation with the degree
fit <- lm(wage ~ poly(age, d), data = Wage)
coef(summary(fit))

                Estimate Std. Error     t value
(Intercept)    111.70361  0.7282103 153.3947154
poly(age, d)1  447.06785 39.8857195  11.2087198
poly(age, d)2 -478.31581 39.8857195 -11.9921569
poly(age, d)3  125.52169 39.8857195   3.1470333
poly(age, d)4  -77.91118 39.8857195  -1.9533603
poly(age, d)5  -35.81289 39.8857195  -0.8978875
poly(age, d)6   62.70772 39.8857195   1.5721847
poly(age, d)7   50.54979 39.8857195   1.2673656
poly(age, d)8  -11.25473 39.8857195  -0.2821745
poly(age, d)9  -83.69180 39.8857195  -2.0982898
                  Pr(>|t|)
(Intercept)   0.000000e+00
poly(age, d)1 1.361797e-28
poly(age, d)2 2.136281e-32
poly(age, d)3 1.665592e-03
poly(age, d)4 5.087006e-02
poly(age, d)5 3.693178e-01
poly(age, d)6 1.160135e-01
poly(age, d)7 2.051233e-01
poly(age, d)8 7.778293e-01
poly(age, d)9 3.596326e-02

agelims <- range(age)
age.grid <- seq(from = agelims[1], to = agelims[2])
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)

par(mfrow = c(1, 1), mar = c(4.5, 4.5, 1, 1), oma = c(0, 0, 4, 0))

plot(age, wage, xlim = agelims, cex = .5, col = "grey")

lines(age.grid, preds$fit, lwd = 2, col = "blue")

matlines(age.grid, se.bands, lwd = 1, col = "blue", lty = 3)

The cross validation uses the 9th degree but the cross validation error curve shows that is it flat after the 3rd degree and coefficient at the 9th degree is a bit weak because of the p value at 0.036. ANOVA testing also shows that the degree of 3 and 4 were the most meaningful in improvements.

(b)

Fit a step function to predict wage using age, and perform crossvalidation to choose the optimal number of cuts. Make a plot of the fit obtained.

# CV to select  number of cuts
set.seed(1)
cv.error.cuts <- rep(NA, 9)
for (i in 2:10) {
  Wage$age.cut <- cut(age, i)
  glm.fit <- glm(wage ~ age.cut, data = Wage)
  cv.error.cuts[i - 1] <- cv.glm(Wage, glm.fit, K = 10)$delta[1]
}

cv.error.cuts

[1] 1734.489 1684.271 1635.552 1632.080 1623.415
[6] 1614.996 1601.318 1613.954 1606.331

plot(2:10, cv.error.cuts, type = "b", pch = 19,
     xlab = "Number of Cuts", ylab = "CV Error")

best.cuts <- which.min(cv.error.cuts) + 1
best.cuts

[1] 8

# fitting the step function
fit <- lm(wage ~ cut(age, best.cuts), data = Wage)
coef(summary(fit))

                               Estimate
(Intercept)                    76.28175
cut(age, best.cuts)(25.8,33.5] 25.83329
cut(age, best.cuts)(33.5,41.2] 40.22568
cut(age, best.cuts)(41.2,49]   43.50112
cut(age, best.cuts)(49,56.8]   40.13583
cut(age, best.cuts)(56.8,64.5] 44.10243
cut(age, best.cuts)(64.5,72.2] 28.94825
cut(age, best.cuts)(72.2,80.1] 15.22418
                               Std. Error
(Intercept)                      2.629812
cut(age, best.cuts)(25.8,33.5]   3.161343
cut(age, best.cuts)(33.5,41.2]   3.049065
cut(age, best.cuts)(41.2,49]     3.018341
cut(age, best.cuts)(49,56.8]     3.176792
cut(age, best.cuts)(56.8,64.5]   3.564299
cut(age, best.cuts)(64.5,72.2]   6.041576
cut(age, best.cuts)(72.2,80.1]   9.781110
                                 t value
(Intercept)                    29.006542
cut(age, best.cuts)(25.8,33.5]  8.171618
cut(age, best.cuts)(33.5,41.2] 13.192791
cut(age, best.cuts)(41.2,49]   14.412262
cut(age, best.cuts)(49,56.8]   12.634076
cut(age, best.cuts)(56.8,64.5] 12.373380
cut(age, best.cuts)(64.5,72.2]  4.791505
cut(age, best.cuts)(72.2,80.1]  1.556488
                                    Pr(>|t|)
(Intercept)                    3.110596e-163
cut(age, best.cuts)(25.8,33.5]  4.440913e-16
cut(age, best.cuts)(33.5,41.2]  1.136044e-38
cut(age, best.cuts)(41.2,49]    1.406253e-45
cut(age, best.cuts)(49,56.8]    1.098741e-35
cut(age, best.cuts)(56.8,64.5]  2.481643e-34
cut(age, best.cuts)(64.5,72.2]  1.736008e-06
cut(age, best.cuts)(72.2,80.1]  1.196978e-01

# plotting the fit
agelims <- range(age)
age.grid <- seq(from = agelims[1], to = agelims[2])
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)

plot(age, wage, xlim = agelims, cex = .5, col = "grey")
lines(age.grid, preds$fit, lwd = 2, col = "blue")

matlines(age.grid, se.bands, lwd = 1, col = "blue", lty = 3)

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)
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 = "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)

coeffs <- coef(fit, id = 6)
names(coeffs)

[1] "(Intercept)" "PrivateYes"  "Room.Board" 
[4] "Terminal"    "perc.alumni" "Expend"     
[7] "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)

gam1 <- 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(gam1, se = T, col = "blue")

Private schools charge more than public schools. Room.board, perc.alumni, and PhD show a linear positive relationship with tuition. Expend is clearly non linear with tuition, as tuition rises at first and then levels out. Grad.Rate tends positive but has a wider confidence band at high graduation rates.

(c)

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

preds <- predict(gam1, College.test)
err <- mean((College.test$Outstate - preds)^2)
err

[1] 3349290

tss <- mean((College.test$Outstate - mean(College.test$Outstate))^2)
rss <- 1 - err / tss
rss

[1] 0.7660016

The GAM achieves a test MSE of about 3,349,290 (RMSE ≈ $1,830) and an R square of 0.766, meaning it explains about 76.6% of the variance in out-of-state tuition. This suggests that Private, Room.Board, PhD, perc.alumni, Expend, and Grad.Rate captures the meaningful variation in tuition.

(d)

For which variables, if any, is there evidence of a non-linear relationship with the response?

summary(gam1)

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 
-7402.89 -1114.45   -12.67  1282.69  7470.60 

(Dispersion Parameter for gaussian family taken to be 3711182)

    Null Deviance: 6989966760 on 387 degrees of freedom
Residual Deviance: 1384271126 on 373 degrees of freedom
AIC: 6987.021 

Number of Local Scoring Iterations: NA 

Anova for Parametric Effects
                        Df     Sum Sq
Private                  1 1778718277
s(Room.Board, df = 2)    1 1577115244
s(PhD, df = 2)           1  322431195
s(perc.alumni, df = 2)   1  336869281
s(Expend, df = 5)        1  530538753
s(Grad.Rate, df = 2)     1   86504998
Residuals              373 1384271126
                          Mean Sq F value
Private                1778718277 479.286
s(Room.Board, df = 2)  1577115244 424.963
s(PhD, df = 2)          322431195  86.881
s(perc.alumni, df = 2)  336869281  90.771
s(Expend, df = 5)       530538753 142.957
s(Grad.Rate, df = 2)     86504998  23.309
Residuals                 3711182        
                          Pr(>F)    
Private                < 2.2e-16 ***
s(Room.Board, df = 2)  < 2.2e-16 ***
s(PhD, df = 2)         < 2.2e-16 ***
s(perc.alumni, df = 2) < 2.2e-16 ***
s(Expend, df = 5)      < 2.2e-16 ***
s(Grad.Rate, df = 2)   2.016e-06 ***
Residuals                           
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Anova for Nonparametric Effects
                       Npar Df  Npar F
(Intercept)                           
Private                               
s(Room.Board, df = 2)        1  1.9157
s(PhD, df = 2)               1  0.9699
s(perc.alumni, df = 2)       1  0.1859
s(Expend, df = 5)            4 20.5075
s(Grad.Rate, df = 2)         1  0.5702
                           Pr(F)    
(Intercept)                         
Private                             
s(Room.Board, df = 2)     0.1672    
s(PhD, df = 2)            0.3253    
s(perc.alumni, df = 2)    0.6666    
s(Expend, df = 5)      2.665e-15 ***
s(Grad.Rate, df = 2)      0.4506    
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Expend shows evidence of non linear relationship with out of state tuition with a very small p value < 0.05 whereas Room.Board, PhD, perc.alumni, and Grad.Rate are linear as their p values are > 0.05.