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.

data(Wage)
data1=Wage

set.seed(100)
deltas <- rep(NA, 10)
for (i in 1:10) {
    fit <- glm(wage ~ poly(age, i), data = data1)
    deltas[i] <- cv.glm(data1, fit, K = 10)$delta[1]
}
plot(1:10, deltas, xlab = "Degree", ylab = "Test MSE", type = "l")
d.min <- which.min(deltas)
points(which.min(deltas), deltas[which.min(deltas)], col = "red", cex = 2, pch = 20)

We note that using a degree equal to 9 will give us the lowest test MSE based on our K-fold cross validation (K=10).

fit1 <- lm(wage ~ age, data = data1)
fit2 <- lm(wage ~ poly(age, 2), data = data1)
fit3 <- lm(wage ~ poly(age, 3), data = data1)
fit4 <- lm(wage ~ poly(age, 4), data = data1)
fit5 <- lm(wage ~ poly(age, 5), data = data1)
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

Based on our ANOVA, we see that models with polynomials higher than cubic are not significant, which is contradictory to what our cv showed. Based on the ANOVA and the K-fold cv, I would recommend using a polynomial of 3 for the predictor age when predicting on wage.

plot(wage ~ age, data = data1, col = "darkgrey")
agelims <- range(data1$age)
age.grid <- seq(from = agelims[1], to = agelims[2])
fit <- lm(wage ~ poly(age, 3), data = data1)
preds <- predict(fit, newdata = list(age = age.grid))
lines(age.grid, preds, col = "green", lwd = 2)

Above is the plot for the polynomial fit to the data.

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

cvs <- rep(NA, 10)
for (i in 2:10) {
    data1$age.cut <- cut(data1$age, i)
    fit <- glm(wage ~ age.cut, data = data1)
    cvs[i] <- cv.glm(data1, fit, K = 10)$delta[1]
}
plot(2:10, cvs[-1], xlab = "Cuts", ylab = "Test MSE", type = "l")
d.min <- which.min(cvs)
points(which.min(cvs), cvs[which.min(cvs)], col = "red", cex = 2, pch = 20)

Based on the above plot, we can see that the lowest Test MSE value occurs when using 8 cuts.

plot(wage ~ age, data = data1, col = "darkgrey")
agelims <- range(data1$age)
age.grid <- seq(from = agelims[1], to = agelims[2])
fit <- glm(wage ~ cut(age, 8), data = data1)
preds <- predict(fit, data.frame(age = age.grid))
lines(age.grid, preds, col = "green", lwd = 2)

Above is the plot for the step function fit to the data.

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 stepwise selection on the training set in order to identify a satisfactory model that uses just a subset of the predictors.

data2 =College
set.seed(100)

index = sample(nrow(data2), nrow(data2)*0.8)
train = data2[index,]
test = data2[-index,]



fit <- regsubsets(Outstate ~ ., data = 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)

Based on the Cp, BIC, and adj \(R^2\) we can see that 6 is the optimal number of variables to include in our model.

fit <- regsubsets(Outstate ~ ., data = data2, 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.

fit <- gam(Outstate ~ Private + s(Room.Board, df=2) + s(perc.alumni, df=2) + 
                   s(PhD, df=2) + s(Expend, df=2) + s(Grad.Rate, df=2), data=train)
par(mfrow=c(2, 3))
plot(fit, se=TRUE, col="blue")

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

preds <- predict(fit, test)
MSE <- mean((test$Outstate - preds)^2)
MSE
## [1] 3503217
tss <- mean((test$Outstate - mean(test$Outstate))^2)
rss <- 1 - MSE / tss
rss
## [1] 0.77415

We obtain a test \(R^2\) of 0.77415 using GAM with 6 predictors.

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

summary(fit)
## 
## Call: gam(formula = Outstate ~ Private + s(Room.Board, df = 2) + s(perc.alumni, 
##     df = 2) + s(PhD, df = 2) + s(Expend, df = 2) + s(Grad.Rate, 
##     df = 2), data = train)
## Deviance Residuals:
##      Min       1Q   Median       3Q      Max 
## -7577.35 -1158.57    10.38  1280.72  8896.24 
## 
## (Dispersion Parameter for gaussian family taken to be 3673569)
## 
##     Null Deviance: 10130673051 on 620 degrees of freedom
## Residual Deviance: 2237202639 on 608.9998 degrees of freedom
## AIC: 11163.66 
## 
## Number of Local Scoring Iterations: NA 
## 
## Anova for Parametric Effects
##                         Df     Sum Sq    Mean Sq F value    Pr(>F)    
## Private                  1 2967958313 2967958313 807.922 < 2.2e-16 ***
## s(Room.Board, df = 2)    1 2181592923 2181592923 593.862 < 2.2e-16 ***
## s(perc.alumni, df = 2)   1  787193426  787193426 214.286 < 2.2e-16 ***
## s(PhD, df = 2)           1  380388969  380388969 103.547 < 2.2e-16 ***
## s(Expend, df = 2)        1  667015965  667015965 181.572 < 2.2e-16 ***
## s(Grad.Rate, df = 2)     1   73074653   73074653  19.892 9.759e-06 ***
## Residuals              609 2237202639    3673569                      
## ---
## 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  2.239   0.13506    
## s(perc.alumni, df = 2)       1  0.432   0.51148    
## s(PhD, df = 2)               1  3.726   0.05405 .  
## s(Expend, df = 2)            1 60.290 3.464e-14 ***
## s(Grad.Rate, df = 2)         1  3.776   0.05245 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Based on our ANOVA analysis, we can see that Expend, Phd, and Grad.Rate have evidence of a non-linear relationship if we assume a p-value of 0.06.