Exercise 5

library(rsample)

## Warning: package 'rsample' was built under R version 4.6.1

library(glmnet)

## Warning: package 'glmnet' was built under R version 4.6.1

## Loading required package: Matrix

## Loaded glmnet 5.0

library(leaps)

## Warning: package 'leaps' was built under R version 4.6.1

library(pls)

## Warning: package 'pls' was built under R version 4.6.1

## 
## Attaching package: 'pls'

## The following object is masked from 'package:stats':
## 
##     loadings

library(ggplot2)
library(plotly)

## Warning: package 'plotly' was built under R version 4.6.1

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ISLR2)

## Warning: package 'ISLR2' was built under R version 4.6.1

get_metrics <- function(actual, predicted) {
  
  actual <- as.numeric(actual)
  predicted <- as.numeric(predicted)
  
  mse <- mean((actual - predicted)^2)
  
  data.frame(
    MSE = mse,
    RMSE = sqrt(mse),
    MAE = mean(abs(actual - predicted)),
    R_Squared = 1 -
      sum((actual - predicted)^2) /
      sum((actual - mean(actual))^2)
  )
}

get_test_results <- function(model, test_data, outcome) {
  
  predictions <- predict(
    model,
    newdata = test_data
  )
  
  actual <- test_data[[outcome]]
  
  metrics <- caret::postResample(
    pred = predictions,
    obs = actual
  )
  
  data.frame(
    MSE = mean((actual - predictions)^2),
    RMSE = unname(metrics["RMSE"]),
    MAE = unname(metrics["MAE"]),
    R_Squared = unname(metrics["Rsquared"])
  )
}

Q2a

iii. With Lasso it’s the best choice when a small number of predictors are associated with the response. The L1 penalty shrinks many coefficients to zero which automatically performs the variable selection. This will produce a simpler and more interpertable model that reduces overfitting.

Q2b

iii. Ridge works well when many or most of the predictors have a moderate or smaller impact on the response. It shrinks the coefficients toward zeo which will reduce variance while keeping infromation from every variable in the model.

Q2c

ii. When the relationship between the predictors and outcome is highly nonlinear squares, ridge, lasso and the other linear methods aren’t able to capture the true pattern. Nonlinear methods are designed to model more complex relationships.

Q9a

data(College)
set.seed(5)

college_split <- initial_split(
  College,
  prop = 0.80,
  strata = Apps
)

college_train <- training(college_split)
college_test <- testing(college_split)

Q9b

college_lm <- lm(
  Apps ~ .,
  data = college_train
)

summary(college_lm)

## 
## Call:
## lm(formula = Apps ~ ., data = college_train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5051.8  -440.4   -31.3   337.3  7047.5 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -730.14682  453.19962  -1.611  0.10768    
## PrivateYes  -458.33505  152.68876  -3.002  0.00280 ** 
## Accept         1.63555    0.04394  37.222  < 2e-16 ***
## Enroll        -1.05569    0.20776  -5.081 5.00e-07 ***
## Top10perc     51.92719    6.21339   8.357 4.43e-16 ***
## Top25perc    -15.46667    5.04137  -3.068  0.00225 ** 
## F.Undergrad    0.06488    0.03615   1.794  0.07324 .  
## P.Undergrad    0.05654    0.03533   1.600  0.11009    
## Outstate      -0.08950    0.02126  -4.210 2.94e-05 ***
## Room.Board     0.10965    0.05422   2.022  0.04358 *  
## Books         -0.03267    0.25687  -0.127  0.89884    
## Personal       0.08671    0.07383   1.174  0.24067    
## PhD          -10.78025    5.27253  -2.045  0.04133 *  
## Terminal      -1.72083    5.84066  -0.295  0.76838    
## S.F.Ratio     29.68427   14.82947   2.002  0.04576 *  
## perc.alumni    3.58905    4.58041   0.784  0.43360    
## Expend         0.08969    0.01466   6.116 1.73e-09 ***
## Grad.Rate     10.08881    3.38278   2.982  0.00298 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1048 on 603 degrees of freedom
## Multiple R-squared:  0.9322, Adjusted R-squared:  0.9303 
## F-statistic: 487.7 on 17 and 603 DF,  p-value: < 2.2e-16

college_lm_pred <- predict(
  college_lm,
  newdata = college_test
)

college_lm_results <- get_metrics(
  actual = college_test$Apps,
  predicted = college_lm_pred
)

college_lm_results

Q9c

x_college_train <- model.matrix(
  Apps ~ .,
  data = college_train
)[, -1]

x_college_test <- model.matrix(
  Apps ~ .,
  data = college_test
)[, -1]

y_college_train <- college_train$Apps
y_college_test <- college_test$Apps

college_ridge <- cv.glmnet(
  x = x_college_train,
  y = y_college_train,
  alpha = 0,
  nfolds = 10,
  standardize = TRUE
)

college_ridge$lambda.min

## [1] 374.5783

college_ridge_pred <- predict(
  college_ridge,
  newx = x_college_test,
  s = "lambda.min"
)

college_ridge_pred <- as.numeric(college_ridge_pred)

college_ridge_results <- get_metrics(
  actual = y_college_test,
  predicted = college_ridge_pred
)

data.frame(
  Model = "Ridge",
  Lambda = college_ridge$lambda.min,
  college_ridge_results
)

Q9d

college_lasso <- cv.glmnet(
  x = x_college_train,
  y = y_college_train,
  alpha = 1,
  nfolds = 10,
  standardize = TRUE
)

college_lasso$lambda.min

## [1] 1.999011

college_lasso_pred <- predict(
  college_lasso,
  newx = x_college_test,
  s = "lambda.min"
)

college_lasso_pred <- as.numeric(college_lasso_pred)

college_lasso_results <- get_metrics(
  actual = y_college_test,
  predicted = college_lasso_pred
)

college_lasso_results

college_lasso_coef <- as.matrix(
  coef(
    college_lasso,
    s = "lambda.min"
  )
)

college_lasso_selected <- college_lasso_coef[
  college_lasso_coef[, 1] != 0,
  ,
  drop = FALSE
]

college_lasso_count <- sum(
  rownames(college_lasso_selected) != "(Intercept)"
)

college_lasso_count

## [1] 17

data.frame(
  Model = "Lasso",
  Lambda = college_lasso$lambda.min,
  Nonzero_Coefficients = college_lasso_count,
  college_lasso_results
)

Q9e

college_pcr <- pcr(
  Apps ~ .,
  data = college_train,
  scale = TRUE,
  validation = "CV"
)

validationplot(
  college_pcr,
  val.type = "MSEP"
)

pcr_cv <- RMSEP(
  college_pcr,
  estimate = "CV"
)

pcr_cv_values <- drop(
  pcr_cv$val[1, 1, ]
)

best_pcr_M <- which.min(pcr_cv_values) - 1

best_pcr_M

## 17 comps 
##       17

college_pcr_pred <- predict(
  college_pcr,
  newdata = college_test,
  ncomp = best_pcr_M
)

college_pcr_results <- get_metrics(
  actual = college_test$Apps,
  predicted = college_pcr_pred
)

data.frame(
  Model = "PCR",
  M_Selected = best_pcr_M,
  college_pcr_results
)

Q9f

college_pls <- plsr(
  Apps ~ .,
  data = college_train,
  scale = TRUE,
  validation = "CV"
)

validationplot(
  college_pls,
  val.type = "MSEP"
)

pls_cv <- RMSEP(
  college_pls,
  estimate = "CV"
)

pls_cv_values <- drop(
  pls_cv$val[1, 1, ]
)

best_pls_M <- which.min(pls_cv_values) - 1

best_pls_M

## 16 comps 
##       16

college_pls_pred <- predict(
  college_pls,
  newdata = college_test,
  ncomp = best_pls_M
)

college_pls_results <- get_metrics(
  actual = college_test$Apps,
  predicted = college_pls_pred
)

data.frame(
  Model = "PLS",
  M_Selected = best_pls_M,
  college_pls_results
)

Q9g

college_results <- bind_rows(
  data.frame(
    Model = "Least Squares",
    college_lm_results
  ),
  data.frame(
    Model = "Ridge",
    college_ridge_results
  ),
  data.frame(
    Model = "Lasso",
    college_lasso_results
  ),
  data.frame(
    Model = "PCR",
    college_pcr_results
  ),
  data.frame(
    Model = "PLS",
    college_pls_results
  )
) |>
  arrange(RMSE)

college_results

college_rmse_plot <- ggplot(
  college_results,
  aes(
    x = reorder(Model, RMSE),
    y = RMSE
  )
) +
  geom_col() +
  coord_flip() +
  labs(
    title = "College Application Model Performance",
    x = NULL,
    y = "Test RMSE"
  ) +
  theme_minimal()

college_rmse_plot

So It looks like overall that Ridge ended up being the best with a RMSE of ~995 so it’s predictions were off by about that many applicants on average. It’s R2 was one of the best with ~.92 so it accounts for ~92% of variation in the test set. The differences aren’t huge though as you can see the R2 are all above .90 and no RMSE is above 1035.

Q11a

library(caret)

## Warning: package 'caret' was built under R version 4.6.1

## Loading required package: lattice

## 
## Attaching package: 'caret'

## The following object is masked from 'package:pls':
## 
##     R2

## The following object is masked from 'package:rsample':
## 
##     calibration

data(Boston)

cv_control <- trainControl(
  method = "repeatedcv",
  number = 10,
  repeats = 3,
  savePredictions = "final"
)

## Least Squares

boston_lm <- train(
  crim ~ .,
  data = Boston,
  method = "lm",
  trControl = cv_control,
  metric = "RMSE"
)

## Best Subset

boston_subset <- train(
  crim ~ .,
  data = Boston,
  method = "leapSeq",
  tuneGrid = data.frame(
    nvmax = 1:(ncol(Boston) - 1)
  ),
  trControl = cv_control,
  metric = "RMSE"
)

boston_subset$bestTune

## Ridge Regression:
ridge_grid <- expand.grid(
  alpha = 0,
  lambda = 10^seq(-4, 4, length.out = 100)
)

boston_ridge <- train(
  crim ~ .,
  data = Boston,
  method = "glmnet",
  preProcess = c("center", "scale"),
  tuneGrid = ridge_grid,
  trControl = cv_control,
  metric = "RMSE"
)

## Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo,
## : There were missing values in resampled performance measures.

## Lasso

lasso_grid <- expand.grid(
  alpha = 1,
  lambda = 10^seq(-4, 4, length.out = 100)
)


boston_lasso <- train(
  crim ~ .,
  data = Boston,
  method = "glmnet",
  preProcess = c("center", "scale"),
  tuneGrid = lasso_grid,
  trControl = cv_control,
  metric = "RMSE"
)

## Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo,
## : There were missing values in resampled performance measures.

## PCR

boston_pcr <- train(
  crim ~ .,
  data = Boston,
  method = "pcr",
  preProcess = c("center", "scale"),
  tuneGrid = data.frame(
    ncomp = 1:(ncol(Boston) - 1)
  ),
  trControl = cv_control,
  metric = "RMSE"
)

## Model Comparison
boston_resamples <- resamples(
  list(
    Least_Squares = boston_lm,
    Best_Subset = boston_subset,
    Ridge = boston_ridge,
    Lasso = boston_lasso,
    PCR = boston_pcr
  )
)



bwplot(
  boston_resamples,
  metric = "RMSE"
)

boston_results <- data.frame(
  Method = c("Least Squares", "Best Subset", "Ridge", "Lasso", "PCR"),
  RMSE = round(
    c(
      min(boston_lm$results$RMSE),
      min(boston_subset$results$RMSE),
      min(boston_ridge$results$RMSE),
      min(boston_lasso$results$RMSE),
      min(boston_pcr$results$RMSE)
    ),
    3
  )
)


boston_results

Q11b

The best model looks to be Best Subset which has a 5.562 RMSE, it also sbustantial has a smaller spread that leans towards a lower RMSE in the box-whisker chart.Overall there’s not a ton of advantage for any of the methods over the others.

Q11c

This just shows that the regularization and dimension reduction in methods didn’t substantially improve the predictiveness over the OLS for the Boston data set. Since all 5 methods got pretty similar RMSE it shows that the predictors contained enough information for least squares to perform comparably to the other methods.