Assignment 6

ISLR Chapter 07 (page 297): 6, 10

Q6. 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(ISLR)
## Warning: package 'ISLR' was built under R version 4.5.3
attach(Wage)
library(boot)
library(ggplot2)
set.seed(35)

cv_error <- rep(0,5)

for (i in 1:5) {
  fit_poly <- glm(wage~poly(age,i),data=Wage)
  cv_error[i] <- cv.glm(Wage, fit_poly, K = 5)$delta[1]
}

print(cv_error)
## [1] 1674.609 1602.414 1595.699 1596.547 1597.229
best_degree <- which.min(cv_error)
best_degree
## [1] 3

Our 5 fold cross validation output shows us the error for each polynomial degree. The error decreases from 1674.609 (degree 1) down to 1595.699 (degree 3) then goes upward trend on 4 amd 5. Lowest degree is 3. Now we can compare this to ANOVA with different polynomoial degrees.

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)

ANOVA results further confirms that 3 degree polynominal is the best choice as far as the lowest error rate. Now we can fit the final model and plot the results.

final_model <- lm(wage ~ poly(age, 3), data = Wage)
summary(final_model)
## 
## Call:
## lm(formula = wage ~ poly(age, 3), data = Wage)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -99.693 -24.562  -5.222  15.096 206.119 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    111.7036     0.7291 153.211  < 2e-16 ***
## poly(age, 3)1  447.0679    39.9335  11.195  < 2e-16 ***
## poly(age, 3)2 -478.3158    39.9335 -11.978  < 2e-16 ***
## poly(age, 3)3  125.5217    39.9335   3.143  0.00169 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 39.93 on 2996 degrees of freedom
## Multiple R-squared:  0.0851, Adjusted R-squared:  0.08419 
## F-statistic: 92.89 on 3 and 2996 DF,  p-value: < 2.2e-16
agelims <- range(Wage$age)
age.grid <- seq(from = agelims[1], to = agelims[2])
preds <- predict(final_model, newdata = list(age = age.grid), se.fit = 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,2,0))
plot(age,wage,xlim=agelims,cex =.5,col="darkgrey")
title("Degree-3 Polynomial",outer=T)

lines(age.grid,preds$fit,lwd=2,col="darkblue")
matlines(age.grid,se.bands,lwd=2,col="pink",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.

set.seed(35)

max_cuts <- 10
cv_errors <- rep(NA, max_cuts)

for(i in 2:max_cuts){
  Wage$age_group <- cut(Wage$age, breaks = i)
  glm_fit <- glm(wage ~ age_group, data = Wage)
  cv_errors[i] <- cv.glm(Wage, glm_fit, K = 10)$delta[1]
}

best_cuts <- which.min(cv_errors)

cat("Best number of cuts:", best_cuts)
## Best number of cuts: 8

Our cross-validation results shows that optimal number of cuts is 8, thereforwe we will fit our final model and plot the results.

final_cut_model <- glm(wage ~ cut(age, best_cuts), data = Wage)
summary(final_cut_model)
## 
## Call:
## glm(formula = wage ~ cut(age, best_cuts), data = Wage)
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      76.282      2.630  29.007  < 2e-16 ***
## cut(age, best_cuts)(25.8,33.5]   25.833      3.161   8.172 4.44e-16 ***
## cut(age, best_cuts)(33.5,41.2]   40.226      3.049  13.193  < 2e-16 ***
## cut(age, best_cuts)(41.2,49]     43.501      3.018  14.412  < 2e-16 ***
## cut(age, best_cuts)(49,56.8]     40.136      3.177  12.634  < 2e-16 ***
## cut(age, best_cuts)(56.8,64.5]   44.102      3.564  12.373  < 2e-16 ***
## cut(age, best_cuts)(64.5,72.2]   28.948      6.042   4.792 1.74e-06 ***
## cut(age, best_cuts)(72.2,80.1]   15.224      9.781   1.556     0.12    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 1597.576)
## 
##     Null deviance: 5222086  on 2999  degrees of freedom
## Residual deviance: 4779946  on 2992  degrees of freedom
## AIC: 30652
## 
## Number of Fisher Scoring iterations: 2
age_grid <- data.frame(
  age = seq(min(Wage$age), max(Wage$age), length = 200)
)

pred <- predict(final_cut_model, newdata = age_grid)

ggplot(Wage, aes(age, wage)) +
  geom_point(alpha = 0.3) +
  geom_line(data = age_grid,aes(age, pred), color = "blue", linewidth = 1) +
  labs(title = paste("Step Function with 8 Cuts"), x = "Age",y = "Wage")

Q10. 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)
attach(College)

set.seed(35)

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)

Cp, BIC and adjr2 plots suggest that size 6 is the minimum size for the subsets.

lm1 <- regsubsets(Outstate ~ ., data = College, method = "forward")
coeffs <- coef(fit, id = 6)
names(coeffs)
## [1] "(Intercept)" "PrivateYes"  "Room.Board"  "Terminal"    "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.

library(mgcv)
## Warning: package 'mgcv' was built under R version 4.5.3
## Loading required package: nlme
## This is mgcv 1.9-4. For overview type '?mgcv'.
fit_gam <- gam(Outstate ~ Private + s(Room.Board) + s(Terminal) + s(perc.alumni) + s(Expend) + s(Grad.Rate), data = College.train)

par(mfrow = c(2, 3))
plot(fit_gam, se = TRUE, col = "blue")

> Out-of-state tuition increases steadily and linearly with higher values of Room.Board, perc.alumni, and Grad.Rate. Terminal displays a flat line, showing almost no effect on tuition rates. The most complex predictor is Expend, which exhibits a strong, wavy, non-linear relationship where tuition climbs sharply at first, dips slightly in the middle range, and then rises again at the highest spending levels.

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

preds <- predict(fit_gam, College.test)
err <- mean((College.test$Outstate - preds)^2)
cat("Error Rate:", err, "\n")
## Error Rate: 4021645
tss <- mean((College.test$Outstate - mean(College.test$Outstate))^2)
rss <- 1 - err / tss
cat("Test R-Squared:", rss)
## Test R-Squared: 0.7546697

The model performs well on test data, explaining 75.5% of the total variation in out-of-state tuition. 4,021,645 mean squared error rate.

  1. For which variables, if any, is there evidence of a non-linear relationship with the response?
summary(fit_gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Outstate ~ Private + s(Room.Board) + s(Terminal) + s(perc.alumni) + 
##     s(Expend) + s(Grad.Rate)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   8565.1      227.2  37.691   <2e-16 ***
## PrivateYes    2578.6      280.1   9.207   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                  edf Ref.df      F  p-value    
## s(Room.Board)  2.242  2.867 12.066 1.01e-06 ***
## s(Terminal)    1.913  2.432  5.619 0.002274 ** 
## s(perc.alumni) 1.179  1.337  9.326 0.000808 ***
## s(Expend)      8.141  8.805 14.356  < 2e-16 ***
## s(Grad.Rate)   5.356  6.516  4.230 0.000281 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.781   Deviance explained = 79.2%
## GCV = 3.6974e+06  Scale est. = 3.4989e+06  n = 388

Summary of our model shows that there are non-linear relationship exists for Expend, Grad.Rate, Room.Board, and Terminal with highly significant p-values (\(p < 0.01\)).