Introduction

This lab will cover polynomial regression and GAMs using linear and non-linear methods. The data will be pulled from the ISLR2 Wage data set. Data analysis methods such as cross-validation, prediction sampling, plots, and outlier analysis will be performed on the data throughout this lab.

6. Wage Dataset

To begin, I start with loading the necessary packages and performing an explortatory data analysis on the Wage data

library(ISLR2)
library(tidyverse)
library(plotly)
library(caret)
library(GGally)
library(leaps)

Wage 
attach(Wage)

# Run a quick visual analysis on the variables that are being used
wage1 <- Wage |> select(age, wage)
ggpairs(wage1)

w <- Wage |> ggplot(
    aes(
        x=age,
        y=wage)
    ) +
    geom_point(
        color = "darkred",
        size = 1.5,
        alpha = .6
    ) +
    labs(
        title = "The Affect of Age on Wages"
    ) +
    theme_light()

ggplotly(w)

**Initial Thoughts: This dataset seems to be heavily populated with information about individuals with ages ranging from 25 - 60 years old. The wage variable seems to have a noticeable upward trend towards individuals making $105,000 to $120k.000

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.

# Fit the basic polynomial model
fit <- lm(wage ~ poly(age,4), data=Wage)

# Fit the amplified model
fit.raw <- lm(wage ~ poly(age,4,raw=TRUE), data=Wage)

# Fit another polynomal model
fit.manual <- lm(wage ~ I(age^2) + I(age^3) + I(age^4), data=Wage)


coef(summary(fit))
##                 Estimate Std. Error    t value     Pr(>|t|)
## (Intercept)    111.70361  0.7287409 153.283015 0.000000e+00
## poly(age, 4)1  447.06785 39.9147851  11.200558 1.484604e-28
## poly(age, 4)2 -478.31581 39.9147851 -11.983424 2.355831e-32
## poly(age, 4)3  125.52169 39.9147851   3.144742 1.678622e-03
## poly(age, 4)4  -77.91118 39.9147851  -1.951938 5.103865e-02
coef(summary(fit.raw))
##                                Estimate   Std. Error   t value     Pr(>|t|)
## (Intercept)               -1.841542e+02 6.004038e+01 -3.067172 0.0021802539
## poly(age, 4, raw = TRUE)1  2.124552e+01 5.886748e+00  3.609042 0.0003123618
## poly(age, 4, raw = TRUE)2 -5.638593e-01 2.061083e-01 -2.735743 0.0062606446
## poly(age, 4, raw = TRUE)3  6.810688e-03 3.065931e-03  2.221409 0.0263977518
## poly(age, 4, raw = TRUE)4 -3.203830e-05 1.641359e-05 -1.951938 0.0510386498
coef(summary(fit.manual))
##                  Estimate   Std. Error   t value     Pr(>|t|)
## (Intercept)  3.117328e+01 6.731628e+00  4.630868 3.796406e-06
## I(age^2)     1.762538e-01 2.068405e-02  8.521243 2.462922e-17
## I(age^3)    -4.051882e-03 5.850462e-04 -6.925747 5.285491e-12
## I(age^4)     2.495627e-05 4.482611e-06  5.567351 2.815024e-08

10-Fold Cross-validation

# 10 fold Cross Validation
ctrl <- trainControl(method = "cv", number = 10)

rmse <- numeric(10)

# Fit the model using Caret and perform a For loop to find the best fit
for (d in 1:10) {

  form <- as.formula(
    paste0("wage ~ poly(age,", d, ")")
  )

  model <- train(
    form,
    data = Wage,
    method = "lm",
    trControl = ctrl
  )

  rmse[d] <- model$results$RMSE
}

# Plot the result
plot(rmse,
     col = "navyblue",
     pch = 16,
     main = "Best CV Model")

Initial Thoughts: The graph shows the result of the best polynomial model using CV. The dot with the lowest RMSE is the best fitting model to the Wage data. This model will be stored inside of the variable called best.degree and ran in the next section as a polynomial model.

Best Fitting Model

# fit the best model
best.degree <- which.min(rmse)

#Fit the best model
fit <- lm(wage ~ poly(age, best.degree), data = Wage)
# Create an age grid
age.grid <- seq(
  min(Wage$age),
  max(Wage$age),
  by = 1
  )

# run predictions
preds <- predict(
  fit,
  newdata = data.frame(
    age = age.grid),
  se = TRUE)


se.bands <- cbind(
  preds $fit + 2 * preds$se.fit,
  preds$fit,
  preds$fit - 2 * preds$se.fit
  )

Predictions

plot(
  Wage$age,
  Wage$wage,
  col = "grey",
  pch = 16,
  xlab = "age",
  ylab = "Wage",
  main = "Prediction of Age Affect on Wage"
)

lines(age.grid,
      preds$fit,
      col="blue",
      lwd=2)

Initial Thoughts: The grey dots in the graph above are the actual data points of the Wage datasets. The blue line is what the polynomial regression has predicted. According to the prediction, there is a gradual increase in wage from individuals aging from 20 - 40 and a steady hold from 40 -65. As individuals age past 65 they begin to slowly taper off with the exception of a slight increase in the 80s.

Choosing Cuts

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.

# 10 fold Cross Validation
ctrl <- trainControl(method = "cv", number = 10)

rmse1 <- numeric(9)

# Fit the model using Caret and perform a For loop to find the best fit
for (i in 2:10) {

   Wage$age.group <- cut(Wage$age, breaks =i)

  model <- train(
    wage ~ age.group,
    data = Wage,
    method = "lm",
    trControl = ctrl
  )

  rmse1[i - 1] <- model$results$RMSE
}

# Plot the result
plot(rmse1,
     col = "navyblue",
     pch = 16,
     main = "Best CV Model")

best.cuts <- which.min(rmse1) + 1

# use the step function to predict the best breaks
Wage$age.group <- cut(Wage$age, breaks = best.cuts)

# fit the best cuts
fit <- lm(wage ~ age.group, data = Wage)

v <- seq(
  min(Wage$age),
  max(Wage$age),
  by = 1
)

age.grid.group <- cut(
  age.grid,
  breaks = best.cuts
)

# Predict wages
preds <- predict(
  fit, newdata = data.frame(age.group = age.grid.group)
)

# Plot the best step fit
z<-Wage |> ggplot(
  aes(
    age,
    wage)) + 
  geom_point(
      col="gray", 
      pch = 16,) +
    geom_line(
      data = data.frame(age = age.grid, wage = preds),
      aes(age, wage),
      color = "blue",
      linewidth = 1
  ) +
    labs(x = "Age", 
         y = "Wage", 
         title = "Step Function Best fit")

ggplotly(z)

Final Thoughts: The step function breaks out the data in this graph into portions that experience the largest spikes in wages. These are plotted above. The breaks in the data where wage spikes occur the most frequently appear to be from (17.9,25.8], (25.8,33.5], (33.5,41.2], (41.2,49], (49,56.8], (56.8,64.5], (64.5,72.2], (72.2,80.1]. The spike from ages 18 - 25 and ages 25 - 33 appear to be where the largest shifts occour.

10. College Dataset

This question uses the College data set.

Test/train split

set.seed(1)

train <- sample(nrow(College), nrow(College) / 2)
College.train <- College[train, ]
College.test  <- College[-train, ]

fit <- regsubsets(
  Outstate ~ .,
  data = College.train,
  method = "forward",
  nvmax = 17
)

fit.summary <- summary(fit)

Model Selection Criteria

x <- data.frame(
  nvars = 1:length(fit.summary$rsq),
  rsq = fit.summary$rsq
)

g<-ggplot(x, aes(x = nvars, y = rsq)) +
  geom_point() +
  geom_line() + 
  labs(title = "Best Variable Amount using Subset Selection",
       x = "r squared",
       y = "number of variables")

ggplotly(g)

Thoughts: The highest R squared is .7794 showing that the best fitting model contain all 17 variables in the model. The plot also shows that as more variables are added, out of state tuition predictions become more accurate.

Could not continue this homework question due to the gam package being able to load in my version of R