Assignment 6

This assignment follows questions modeled in chapter 7 of the ISLR textbook for data analysis.

Question 6

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

a.) Polynomial Regression

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(caret)
library(ggplot2)

data(Wage)

set.seed(1)

ctrl <- trainControl(method = "cv", number = 10)

poly.models <- list()
cv.errors <- numeric(10)

for(i in 1:10) {
  poly.formula <- as.formula(
    paste("wage~poly(age,", i, ")")
  )
  poly.models[[i]] <- train(
    poly.formula, 
    data = Wage, 
    method = "lm",
    trControl = ctrl
  )
  cv.errors[i] <- min(poly.models[[i]]$results$RMSE^2)
}
cv.errors

 [1] 1670.277 1596.717 1582.212 1581.502 1584.812 1578.691
 [7] 1591.044 1585.468 1586.835 1583.362

best.degree <- which.min(poly.models[[i]]$results$RMSE^2)
cat("Degree Chosen:", best.degree, "\n")

Degree Chosen: 1 

poly.fit <- lm(wage ~ poly(age, best.degree), data = Wage)
summary(poly.fit)

Call:
lm(formula = wage ~ poly(age, best.degree), data = Wage)

Residuals:
     Min       1Q   Median       3Q      Max 
-100.265  -25.115   -6.063   16.601  205.748 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)            111.7036     0.7473  149.48   <2e-16 ***
poly(age, best.degree) 447.0679    40.9291   10.92   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 40.93 on 2998 degrees of freedom
Multiple R-squared:  0.03827,   Adjusted R-squared:  0.03795 
F-statistic: 119.3 on 1 and 2998 DF,  p-value: < 2.2e-16

#Compare to anova
fit1 <- lm(wage ~ age, 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)

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

#plot
age.grid <- data.frame(age = seq(min(Wage$age), 
                                 max(Wage$age), 
                                 length = 100))
pred <- predict(poly.fit, newdata = age.grid)

ggplot(Wage, aes(age, wage))+
  geom_point(alpha = .4)+
  geom_line(data = age.grid,
            aes(age, pred),
            color = "blue",
            linewidth = 1.2)+
  labs(title = paste("Polynomial Regression (Degree", best.degree, ")"))

b.) Fit a step function

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.

cuts <- 2:10
step.models <- list()
step.errors <- numeric(length(cuts))

for(i in seq_along(cuts)){
  step.formula <- as.formula(
    paste("wage ~ cut(age,", cuts[i], ")")
  )
  step.models[[i]] <- train(
    step.formula,
    data = Wage,
    method = "lm", 
    trControl = ctrl
  )
  step.errors[i] <- min(step.models[[i]]$results$RMSE^2)
}
step.errors

[1] 1724.581 1674.529 1613.872 1623.198 1614.032 1603.570 1594.231
[8] 1600.528 1602.152

best.cuts <- cuts[which.min(step.errors)]
cat("Best cuts:", best.cuts, "\n")

Best cuts: 8 

step.fit <- lm(
  wage ~ cut(age, best.cuts),
  data = Wage
)

summary(step.fit)

Call:
lm(formula = wage ~ cut(age, best.cuts), data = Wage)

Residuals:
    Min      1Q  Median      3Q     Max 
-99.697 -24.552  -5.307  15.417 198.560 

Coefficients:
                               Estimate Std. Error t value
(Intercept)                      76.282      2.630  29.007
cut(age, best.cuts)(25.8,33.5]   25.833      3.161   8.172
cut(age, best.cuts)(33.5,41.2]   40.226      3.049  13.193
cut(age, best.cuts)(41.2,49]     43.501      3.018  14.412
cut(age, best.cuts)(49,56.8]     40.136      3.177  12.634
cut(age, best.cuts)(56.8,64.5]   44.102      3.564  12.373
cut(age, best.cuts)(64.5,72.2]   28.948      6.042   4.792
cut(age, best.cuts)(72.2,80.1]   15.224      9.781   1.556
                               Pr(>|t|)    
(Intercept)                     < 2e-16 ***
cut(age, best.cuts)(25.8,33.5] 4.44e-16 ***
cut(age, best.cuts)(33.5,41.2]  < 2e-16 ***
cut(age, best.cuts)(41.2,49]    < 2e-16 ***
cut(age, best.cuts)(49,56.8]    < 2e-16 ***
cut(age, best.cuts)(56.8,64.5]  < 2e-16 ***
cut(age, best.cuts)(64.5,72.2] 1.74e-06 ***
cut(age, best.cuts)(72.2,80.1]     0.12    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 39.97 on 2992 degrees of freedom
Multiple R-squared:  0.08467,   Adjusted R-squared:  0.08253 
F-statistic: 39.54 on 7 and 2992 DF,  p-value: < 2.2e-16

#plot

step.pred <-  predict(
  step.fit, 
  newdata = data.frame(age = age.grid)
)

step.plot <- data.frame(
  age = age.grid, 
  wage = step.pred
)

ggplot(Wage, aes(age, wage)) +
  geom_point(alpha = .4) +
  geom_line(
    data = step.plot,
    aes(age, wage), 
    color = "blue",
    linewidth = 1.2
  ) +
  labs(title = paste("Step Function (", best.cuts, "cuts)", sep = ""),
       xlab = "Age",
       ylab = "Wage")

Question 10

This Question relates to the College data set

a.) Stepwise Selection

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

library(ISLR2)
library(caret)
library(MASS)
library(gam)

data(College)

set.seed(1)

train.index <- createDataPartition(College$Outstate,
                                   p = 0.7,
                                   list = FALSE)

College.train <- College[train.index, ]
College.test  <- College[-train.index, ]

full.fit <- lm(Outstate ~ ., data = College.train)

step.fit <- stepAIC(full.fit,
                    direction = "both",
                    trace = FALSE)

summary(step.fit)

b.) GAM

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.

gam.fit <- gam(
  Outstate ~
    Private +
    lo(Room.Board) +
    lo(Expend) +
    lo(PhD) +
    lo(perc.alumni) +
    lo(Grad.Rate),
  data = College.train
)

summary(gam.fit)

Call: gam(formula = Outstate ~ Private + lo(Room.Board) + lo(Expend) + 
    lo(PhD) + lo(perc.alumni) + lo(Grad.Rate), data = College.train)
Deviance Residuals:
   Min     1Q Median     3Q    Max 
 -6785  -1139     62   1229   4768 

(Dispersion Parameter for gaussian family taken to be 3242650)

    Null Deviance: 8648835880 on 545 degrees of freedom
Residual Deviance: 1700199687 on 524.3242 degrees of freedom
AIC: 9758.293 

Number of Local Scoring Iterations: NA 

Anova for Parametric Effects
                    Df     Sum Sq    Mean Sq F value    Pr(>F)    
Private           1.00 2619605560 2619605560 807.860 < 2.2e-16 ***
lo(Room.Board)    1.00 1684465323 1684465323 519.472 < 2.2e-16 ***
lo(Expend)        1.00 1284709139 1284709139 396.191 < 2.2e-16 ***
lo(PhD)           1.00  174755812  174755812  53.893 8.128e-13 ***
lo(perc.alumni)   1.00  160221174  160221174  49.411 6.494e-12 ***
lo(Grad.Rate)     1.00   78079935   78079935  24.079 1.236e-06 ***
Residuals       524.32 1700199687    3242650                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Anova for Nonparametric Effects
                Npar Df  Npar F    Pr(F)    
(Intercept)                                 
Private                                     
lo(Room.Board)      3.0  3.4006  0.01816 *  
lo(Expend)          4.0 20.3418 1.11e-15 ***
lo(PhD)             2.8  0.9328  0.41989    
lo(perc.alumni)     2.3  0.7208  0.50664    
lo(Grad.Rate)       2.6  2.0513  0.11582    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

par(mfrow = c(2,3))

plot(gam.fit,
     se = TRUE,
     col = "blue")

c.) Results

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

gam.pred <- predict(gam.fit,
                    newdata = College.test)

gam.mse <- mean((gam.pred - College.test$Outstate)^2)

gam.rmse <- sqrt(gam.mse)

gam.mse

[1] 4193772

gam.rmse

[1] 2047.87

The model used produced a mean squared error (MSE) of 4,193,773, while the root mean squared error (RMSE) produced was 2,047. Meaning the predicted tuition differs $2,047 on average when compared to the actual tuition values.

d.) Non-linearity

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

summary(gam.fit)

Call: gam(formula = Outstate ~ Private + lo(Room.Board) + lo(Expend) + 
    lo(PhD) + lo(perc.alumni) + lo(Grad.Rate), data = College.train)
Deviance Residuals:
   Min     1Q Median     3Q    Max 
 -6785  -1139     62   1229   4768 

(Dispersion Parameter for gaussian family taken to be 3242650)

    Null Deviance: 8648835880 on 545 degrees of freedom
Residual Deviance: 1700199687 on 524.3242 degrees of freedom
AIC: 9758.293 

Number of Local Scoring Iterations: NA 

Anova for Parametric Effects
                    Df     Sum Sq    Mean Sq F value    Pr(>F)    
Private           1.00 2619605560 2619605560 807.860 < 2.2e-16 ***
lo(Room.Board)    1.00 1684465323 1684465323 519.472 < 2.2e-16 ***
lo(Expend)        1.00 1284709139 1284709139 396.191 < 2.2e-16 ***
lo(PhD)           1.00  174755812  174755812  53.893 8.128e-13 ***
lo(perc.alumni)   1.00  160221174  160221174  49.411 6.494e-12 ***
lo(Grad.Rate)     1.00   78079935   78079935  24.079 1.236e-06 ***
Residuals       524.32 1700199687    3242650                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Anova for Nonparametric Effects
                Npar Df  Npar F    Pr(F)    
(Intercept)                                 
Private                                     
lo(Room.Board)      3.0  3.4006  0.01816 *  
lo(Expend)          4.0 20.3418 1.11e-15 ***
lo(PhD)             2.8  0.9328  0.41989    
lo(perc.alumni)     2.3  0.7208  0.50664    
lo(Grad.Rate)       2.6  2.0513  0.11582    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Based on Anova for Nonparametric Effects, variables that show evidence of non-linearity include Expend (p = 1.1e-15) and Room.Board (p = .0182). With Expend showing significant non-linearity.