Assignment 6

Problem 3

(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)
attach(Wage)

cv_errors <- rep(NA, 10)
for (d in 1:10) {
  glm_fit <- glm(wage ~ poly(age, d), data = Wage)
  cv_errors[d] <- cv.glm(Wage, glm_fit, K = 10)$delta[1]
}
which.min(cv_errors) #7

[1] 5

fit1 <- lm(wage ~ poly(age, 1), 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) #3

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

bestfit <- lm(wage~poly(age,7), data = Wage)

agelims=range(age)
age.grid=seq(from=agelims[1],to=agelims[2])
preds=predict(bestfit,newdata = list(age=age.grid),se=T)
se.bands=cbind(preds$fit+2*preds$se.fit,preds$fit-2*preds$se.fit)
plot(age,wage,xlim=agelims,cex=.5,col='darkgrey')
title('Degree-7 Polynomial')

lines(age.grid, preds$fit, lwd=2, col='steelblue')
matlines(age.grid,se.bands,lwd = 1, col = 'lightblue',lty = 3)

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

library(boot)
attach(Wage)

set.seed(1)

cv_errors_step <- rep(NA, 10)

for (k in 2:10) {
  Wage$age_cut <- cut(Wage$age, k)
  step_fit <- glm(wage ~ age_cut, data = Wage)
  cv_errors_step[k] <- cv.glm(Wage, step_fit, K = 10)$delta[1]
}
cv_errors_step

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

which.min(cv_errors_step)

[1] 8

step_best <- lm(wage~cut(age,8), data = Wage)
plot(wage ~ age, data = Wage, col = "darkgrey",main = paste("Step Function Fit with 8 Cuts"))
preds_step=predict(step_best,newdata = list(age=age.grid),se=T)
se.bands_step=cbind(preds_step$fit+2*preds_step$se.fit,preds_step$fit-2*preds_step$se.fit)
lines(age.grid, preds_step$fit, col = "red", lwd = 2)

matlines(age.grid,se.bands_step,lwd = 1, col = 'hotpink',lty = 3)

Problem 10

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

names(coef(fit,6))

[1] "(Intercept)" "PrivateYes"  "Room.Board"  "Terminal"    "perc.alumni" "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,4)+s(PhD,4)+s(perc.alumni,4)+s(Expend,4)+s(Grad.Rate,4),data = College.train)
par(mfrow=c(2,3))
plot(gam_fit,se=T,col="steelblue")

• By evaluating the charts you can see that some are actually pretty close to having a linear relationship with Outstate. The one that stands out to be non-linear is the Expend variable. This suggests that a GAM may better capture the relationship between Outstate and Expend.

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

gam_pred <- predict(gam_fit,newdata = College.test)
gam_mse <- mean((College.test$Outstate - gam_pred)^2)
gam_mse

[1] 3324814

• The model has an MSE of 3324814

(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, 4) + s(PhD, 
    4) + s(perc.alumni, 4) + s(Expend, 4) + s(Grad.Rate, 4), 
    data = College.train)
Deviance Residuals:
     Min       1Q   Median       3Q      Max 
-7174.12 -1141.26   -86.67  1264.75  7506.52 

(Dispersion Parameter for gaussian family taken to be 3750826)

    Null Deviance: 6989966760 on 387 degrees of freedom
Residual Deviance: 1372800367 on 365.9994 degrees of freedom
AIC: 6997.793 

Number of Local Scoring Iterations: NA 

Anova for Parametric Effects
                   Df     Sum Sq    Mean Sq F value    Pr(>F)    
Private             1 1761501708 1761501708 469.630 < 2.2e-16 ***
s(Room.Board, 4)    1 1563405241 1563405241 416.816 < 2.2e-16 ***
s(PhD, 4)           1  334516516  334516516  89.185 < 2.2e-16 ***
s(perc.alumni, 4)   1  331488702  331488702  88.377 < 2.2e-16 ***
s(Expend, 4)        1  525593541  525593541 140.127 < 2.2e-16 ***
s(Grad.Rate, 4)     1   87873241   87873241  23.428 1.916e-06 ***
Residuals         366 1372800367    3750826                      
---
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.9158    0.1266    
s(PhD, 4)               3  0.8504    0.4671    
s(perc.alumni, 4)       3  0.3520    0.7877    
s(Expend, 4)            3 23.6176 5.373e-14 ***
s(Grad.Rate, 4)         3  1.0701    0.3617    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

• The ANOVA results show a strong non-linear relationship between Outstate and Expend. The other predictors do not show significant evidence of a non-linear relationship.