Assignment #6: Moving Beyond Linearity

Question 6. In this exercise, you will further analyze the Wage data set considered throughout this chapter.

library(ISLR2)
attach(Wage)

(a) Perform polynomial regression to predict wage using age. Use cross-validation to select the optimal degree d for the polyno- mial. 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(boot)
set.seed(15)
cv.error = rep(0,5)

for (i in 1:5){
glm.fit = glm(wage ~ poly(age,i),data=Wage)
cv.error[i]= cv.glm(Wage,glm.fit,K=10)$delta[1]
}
cv.error

[1] 1676.192 1599.786 1596.256 1593.755 1595.106

plot(cv.error, type="b", xlab="Degree", ylab="TestMSE",lwd = 2)
points(which.min(cv.error), cv.error[4], col="blue", pch="X", cex=2)

• Based on cross-validation, the optimal degree is 4. The comparison by ANOVA suggest degree 4 is enough

age.range <- range(age)
age.grid <- seq(from=age.range[1], to=age.range[2])
lm.fit <- lm(wage ~ poly(age, 3), data = Wage)
pred1 <- predict(lm.fit, newdata = list(age = age.grid))

plot(wage ~ age, data = Wage, xlim=age.range, cex=.5, col="darkgrey", xlab = "Age", ylab = "Wage")
lines(age.grid, pred1, 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.

set.seed(2)
cv.errors = rep(NA, 10)
for(i in 2:10){
  Wage$age.cut <- cut(Wage$age,i)
  glm.fit <- glm(wage ~ age.cut, data=Wage)
  cv.errors[i] <- cv.glm(Wage, glm.fit, K=10)$delta[1]
}

cv.errors

set.seed(2)
cv.errors = rep(NA, 10)
for(i in 2:10){
  Wage$age.cut <- cut(Wage$age,i)
  glm.fit <- glm(wage ~ age.cut, data=Wage)
  cv.errors[i] <- cv.glm(Wage, glm.fit, K=10)$delta[1]
}

cv.errors

 [1]       NA 1734.448 1681.623 1635.869 1632.864 1621.643 1613.148 1601.850 1610.852 1603.607

• 8 is the number of optimal cuts.

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

---

Question 10. This question relates to the College data set.

library(leaps)
attach(College)

(a) Split the data into a training set and a test set.Using out-of-state tuition as the response and the other variabgout-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.

set.seed(5)
train = sample(nrow(College) * 0.7)
train.college = College[train, ]
test.college = College[-train, ]

forward.subset = regsubsets(Outstate ~ ., data = train.college, nvmax = ncol(College)-1, method = "forward")
model.summary = summary(forward.subset)

plot.metric = function(metric, yaxis_label, reverse = FALSE) {
  plot(metric, xlab = "# of Variables", ylab = yaxis_label, xaxt = "n", type = "l")
  axis(side = 1, at = 1:length(metric))
  
  if (reverse) {
    metric_1se = max(metric) - (sd(metric) / sqrt(length(metric)))
    min_subset = which(metric > metric_1se)
  } else {
    metric_1se = min(metric) + (sd(metric) / sqrt(length(metric)))
    min_subset = which(metric < metric_1se)
  }
  
  abline(h = metric_1se, col = "red", lty = 2)
  abline(v = min_subset[1], col = "blue", lty = 2)
}

par(mfrow=c(1, 3))
  
plot.metric(model.summary$cp, "Cp")
plot.metric(model.summary$bic, "BIC")
plot.metric(model.summary$adjr2, "Adjusted R2", reverse = TRUE)

*Based on the models, Cp and Adjusted R2 have a subset with 7 variables and BIC has 6.

coef(forward.subset, 6)

(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.model = 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.college)
par(mfrow=c(2, 3))
plot(gam.model, se=TRUE, col="blue")

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

calc_mse = function(y, y_hat) {
  return(mean((y - y_hat)^2))
}

calc_rmse = function(y, y_hat) {
  return(sqrt(calc_mse(y, y_hat)))
}

calc_r2 = function(y, y_hat) {
  y_bar = mean(y)
  rss = sum((y - y_hat)^2)
  tss = sum((y - y_bar)^2)
  return(1 - (rss / tss))
}

gam.pred = predict(gam.model, test.college)

gam.rmse = calc_rmse(test.college$Outstate, gam.pred)
cat("RMSE:", gam.rmse, "\n")

gam_r2 = calc_r2(test.college$Outstate, gam.pred)
cat(" R^2:", gam_r2, "\n")

• We get a test RMSE of 1984.38 and R^2 of 0.7614 using GAM with 6 predictors.

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

summary(gam.model)

• The results of the non-parametric ANOVA test indicate that there is compelling evidence of a non-linear relationship between the response variable and both Expend and PhD, and a moderately strong non-linear relationship between the response variable and both Room.Board and Grad.Rate.