The answer for part a would be iii. The lasso has a penalty tearm that limits the coefficient estimates which can shrinl some of them to 0. Since it is less flexible it will introduce a small amount of bias. The MSE will decrease as long as the increase in bias is smaller than the corresponding decrease in variance.
The answer for part b is iii. This also introduces a penalty term that shrinks the coefficients toward 0 limiting the model’s freedom which makes it less flexible.
The answer for part c is ii. Non-linear methods can show complex relationships. This gives it a wider range of shapes to fit the data making them more flexible. Fitting the data more closely reduces bias.
library(ISLR2)
## Warning: package 'ISLR2' was built under R version 4.3.3
data(College)
head(College)
## Private Apps Accept Enroll Top10perc Top25perc
## Abilene Christian University Yes 1660 1232 721 23 52
## Adelphi University Yes 2186 1924 512 16 29
## Adrian College Yes 1428 1097 336 22 50
## Agnes Scott College Yes 417 349 137 60 89
## Alaska Pacific University Yes 193 146 55 16 44
## Albertson College Yes 587 479 158 38 62
## F.Undergrad P.Undergrad Outstate Room.Board Books
## Abilene Christian University 2885 537 7440 3300 450
## Adelphi University 2683 1227 12280 6450 750
## Adrian College 1036 99 11250 3750 400
## Agnes Scott College 510 63 12960 5450 450
## Alaska Pacific University 249 869 7560 4120 800
## Albertson College 678 41 13500 3335 500
## Personal PhD Terminal S.F.Ratio perc.alumni Expend
## Abilene Christian University 2200 70 78 18.1 12 7041
## Adelphi University 1500 29 30 12.2 16 10527
## Adrian College 1165 53 66 12.9 30 8735
## Agnes Scott College 875 92 97 7.7 37 19016
## Alaska Pacific University 1500 76 72 11.9 2 10922
## Albertson College 675 67 73 9.4 11 9727
## Grad.Rate
## Abilene Christian University 60
## Adelphi University 56
## Adrian College 54
## Agnes Scott College 59
## Alaska Pacific University 15
## Albertson College 55
set.seed(42)
train_indices <- sample(nrow(College), 0.8 * nrow(College))
college_train <- College[train_indices, ]
college_test <- College[-train_indices, ]
lm_fit <- lm(Apps ~ ., data = college_train)
lm_pred <- predict(lm_fit, newdata = college_test)
lm_testerror <- mean((college_test$Apps - lm_pred)^2)
print(lm_testerror)
## [1] 1941715
#RMSE
RMSE <- sqrt(lm_testerror)
print(RMSE)
## [1] 1393.454
library(glmnet)
## Warning: package 'glmnet' was built under R version 4.3.3
## Loading required package: Matrix
## Loaded glmnet 4.1-8
library(caret)
## Warning: package 'caret' was built under R version 4.3.3
## Loading required package: ggplot2
## Loading required package: lattice
## Warning: package 'lattice' was built under R version 4.3.3
library(Matrix)
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.3.3
##
## 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(ggplot2)
X_train <- model.matrix(Apps ~ ., data = college_train)[, -1]
y_train <- college_train$Apps
X_test <- model.matrix(Apps ~ ., data = college_test)[, -1]
y_test <- college_test$Apps
set.seed(42)
cv_ridge <- cv.glmnet(X_train, y_train, alpha = 0)
lambda <- cv_ridge$lambda.min
print(lambda)
## [1] 337.0816
plot(cv_ridge)
ridge.pred <- predict(cv_ridge, , newx = X_test)
ridge_testerror <- mean((y_test - ridge.pred)^2)
print(ridge_testerror)
## [1] 4683980
fitControl <- trainControl(
method = "repeatedcv",
number = 10,
repeats = 3,
verboseIter = TRUE
)
ridgeGrid <- expand.grid(alpha = 0,
lambda = 10^seq(-3, 1, length = 100))
ridge_caret <- train(
Apps ~ .,
data = college_test,
method = "glmnet",
tuneGrid = ridgeGrid,
trControl = fitControl,
preProcess = c("center", "scale")
)
## + Fold01.Rep1: alpha=0, lambda=10
## - Fold01.Rep1: alpha=0, lambda=10
## + Fold02.Rep1: alpha=0, lambda=10
## - Fold02.Rep1: alpha=0, lambda=10
## + Fold03.Rep1: alpha=0, lambda=10
## - Fold03.Rep1: alpha=0, lambda=10
## + Fold04.Rep1: alpha=0, lambda=10
## - Fold04.Rep1: alpha=0, lambda=10
## + Fold05.Rep1: alpha=0, lambda=10
## - Fold05.Rep1: alpha=0, lambda=10
## + Fold06.Rep1: alpha=0, lambda=10
## - Fold06.Rep1: alpha=0, lambda=10
## + Fold07.Rep1: alpha=0, lambda=10
## - Fold07.Rep1: alpha=0, lambda=10
## + Fold08.Rep1: alpha=0, lambda=10
## - Fold08.Rep1: alpha=0, lambda=10
## + Fold09.Rep1: alpha=0, lambda=10
## - Fold09.Rep1: alpha=0, lambda=10
## + Fold10.Rep1: alpha=0, lambda=10
## - Fold10.Rep1: alpha=0, lambda=10
## + Fold01.Rep2: alpha=0, lambda=10
## - Fold01.Rep2: alpha=0, lambda=10
## + Fold02.Rep2: alpha=0, lambda=10
## - Fold02.Rep2: alpha=0, lambda=10
## + Fold03.Rep2: alpha=0, lambda=10
## - Fold03.Rep2: alpha=0, lambda=10
## + Fold04.Rep2: alpha=0, lambda=10
## - Fold04.Rep2: alpha=0, lambda=10
## + Fold05.Rep2: alpha=0, lambda=10
## - Fold05.Rep2: alpha=0, lambda=10
## + Fold06.Rep2: alpha=0, lambda=10
## - Fold06.Rep2: alpha=0, lambda=10
## + Fold07.Rep2: alpha=0, lambda=10
## - Fold07.Rep2: alpha=0, lambda=10
## + Fold08.Rep2: alpha=0, lambda=10
## - Fold08.Rep2: alpha=0, lambda=10
## + Fold09.Rep2: alpha=0, lambda=10
## - Fold09.Rep2: alpha=0, lambda=10
## + Fold10.Rep2: alpha=0, lambda=10
## - Fold10.Rep2: alpha=0, lambda=10
## + Fold01.Rep3: alpha=0, lambda=10
## - Fold01.Rep3: alpha=0, lambda=10
## + Fold02.Rep3: alpha=0, lambda=10
## - Fold02.Rep3: alpha=0, lambda=10
## + Fold03.Rep3: alpha=0, lambda=10
## - Fold03.Rep3: alpha=0, lambda=10
## + Fold04.Rep3: alpha=0, lambda=10
## - Fold04.Rep3: alpha=0, lambda=10
## + Fold05.Rep3: alpha=0, lambda=10
## - Fold05.Rep3: alpha=0, lambda=10
## + Fold06.Rep3: alpha=0, lambda=10
## - Fold06.Rep3: alpha=0, lambda=10
## + Fold07.Rep3: alpha=0, lambda=10
## - Fold07.Rep3: alpha=0, lambda=10
## + Fold08.Rep3: alpha=0, lambda=10
## - Fold08.Rep3: alpha=0, lambda=10
## + Fold09.Rep3: alpha=0, lambda=10
## - Fold09.Rep3: alpha=0, lambda=10
## + Fold10.Rep3: alpha=0, lambda=10
## - Fold10.Rep3: alpha=0, lambda=10
## Aggregating results
## Selecting tuning parameters
## Fitting alpha = 0, lambda = 10 on full training set
print(ridge_caret)
## glmnet
##
## 156 samples
## 17 predictor
##
## Pre-processing: centered (17), scaled (17)
## Resampling: Cross-Validated (10 fold, repeated 3 times)
## Summary of sample sizes: 141, 140, 140, 140, 141, 140, ...
## Resampling results across tuning parameters:
##
## lambda RMSE Rsquared MAE
## 0.001000000 1524.014 0.9182428 800.2012
## 0.001097499 1524.014 0.9182428 800.2012
## 0.001204504 1524.014 0.9182428 800.2012
## 0.001321941 1524.014 0.9182428 800.2012
## 0.001450829 1524.014 0.9182428 800.2012
## 0.001592283 1524.014 0.9182428 800.2012
## 0.001747528 1524.014 0.9182428 800.2012
## 0.001917910 1524.014 0.9182428 800.2012
## 0.002104904 1524.014 0.9182428 800.2012
## 0.002310130 1524.014 0.9182428 800.2012
## 0.002535364 1524.014 0.9182428 800.2012
## 0.002782559 1524.014 0.9182428 800.2012
## 0.003053856 1524.014 0.9182428 800.2012
## 0.003351603 1524.014 0.9182428 800.2012
## 0.003678380 1524.014 0.9182428 800.2012
## 0.004037017 1524.014 0.9182428 800.2012
## 0.004430621 1524.014 0.9182428 800.2012
## 0.004862602 1524.014 0.9182428 800.2012
## 0.005336699 1524.014 0.9182428 800.2012
## 0.005857021 1524.014 0.9182428 800.2012
## 0.006428073 1524.014 0.9182428 800.2012
## 0.007054802 1524.014 0.9182428 800.2012
## 0.007742637 1524.014 0.9182428 800.2012
## 0.008497534 1524.014 0.9182428 800.2012
## 0.009326033 1524.014 0.9182428 800.2012
## 0.010235310 1524.014 0.9182428 800.2012
## 0.011233240 1524.014 0.9182428 800.2012
## 0.012328467 1524.014 0.9182428 800.2012
## 0.013530478 1524.014 0.9182428 800.2012
## 0.014849683 1524.014 0.9182428 800.2012
## 0.016297508 1524.014 0.9182428 800.2012
## 0.017886495 1524.014 0.9182428 800.2012
## 0.019630407 1524.014 0.9182428 800.2012
## 0.021544347 1524.014 0.9182428 800.2012
## 0.023644894 1524.014 0.9182428 800.2012
## 0.025950242 1524.014 0.9182428 800.2012
## 0.028480359 1524.014 0.9182428 800.2012
## 0.031257158 1524.014 0.9182428 800.2012
## 0.034304693 1524.014 0.9182428 800.2012
## 0.037649358 1524.014 0.9182428 800.2012
## 0.041320124 1524.014 0.9182428 800.2012
## 0.045348785 1524.014 0.9182428 800.2012
## 0.049770236 1524.014 0.9182428 800.2012
## 0.054622772 1524.014 0.9182428 800.2012
## 0.059948425 1524.014 0.9182428 800.2012
## 0.065793322 1524.014 0.9182428 800.2012
## 0.072208090 1524.014 0.9182428 800.2012
## 0.079248290 1524.014 0.9182428 800.2012
## 0.086974900 1524.014 0.9182428 800.2012
## 0.095454846 1524.014 0.9182428 800.2012
## 0.104761575 1524.014 0.9182428 800.2012
## 0.114975700 1524.014 0.9182428 800.2012
## 0.126185688 1524.014 0.9182428 800.2012
## 0.138488637 1524.014 0.9182428 800.2012
## 0.151991108 1524.014 0.9182428 800.2012
## 0.166810054 1524.014 0.9182428 800.2012
## 0.183073828 1524.014 0.9182428 800.2012
## 0.200923300 1524.014 0.9182428 800.2012
## 0.220513074 1524.014 0.9182428 800.2012
## 0.242012826 1524.014 0.9182428 800.2012
## 0.265608778 1524.014 0.9182428 800.2012
## 0.291505306 1524.014 0.9182428 800.2012
## 0.319926714 1524.014 0.9182428 800.2012
## 0.351119173 1524.014 0.9182428 800.2012
## 0.385352859 1524.014 0.9182428 800.2012
## 0.422924287 1524.014 0.9182428 800.2012
## 0.464158883 1524.014 0.9182428 800.2012
## 0.509413801 1524.014 0.9182428 800.2012
## 0.559081018 1524.014 0.9182428 800.2012
## 0.613590727 1524.014 0.9182428 800.2012
## 0.673415066 1524.014 0.9182428 800.2012
## 0.739072203 1524.014 0.9182428 800.2012
## 0.811130831 1524.014 0.9182428 800.2012
## 0.890215085 1524.014 0.9182428 800.2012
## 0.977009957 1524.014 0.9182428 800.2012
## 1.072267222 1524.014 0.9182428 800.2012
## 1.176811952 1524.014 0.9182428 800.2012
## 1.291549665 1524.014 0.9182428 800.2012
## 1.417474163 1524.014 0.9182428 800.2012
## 1.555676144 1524.014 0.9182428 800.2012
## 1.707352647 1524.014 0.9182428 800.2012
## 1.873817423 1524.014 0.9182428 800.2012
## 2.056512308 1524.014 0.9182428 800.2012
## 2.257019720 1524.014 0.9182428 800.2012
## 2.477076356 1524.014 0.9182428 800.2012
## 2.718588243 1524.014 0.9182428 800.2012
## 2.983647240 1524.014 0.9182428 800.2012
## 3.274549163 1524.014 0.9182428 800.2012
## 3.593813664 1524.014 0.9182428 800.2012
## 3.944206059 1524.014 0.9182428 800.2012
## 4.328761281 1524.014 0.9182428 800.2012
## 4.750810162 1524.014 0.9182428 800.2012
## 5.214008288 1524.014 0.9182428 800.2012
## 5.722367659 1524.014 0.9182428 800.2012
## 6.280291442 1524.014 0.9182428 800.2012
## 6.892612104 1524.014 0.9182428 800.2012
## 7.564633276 1524.014 0.9182428 800.2012
## 8.302175681 1524.014 0.9182428 800.2012
## 9.111627561 1524.014 0.9182428 800.2012
## 10.000000000 1524.014 0.9182428 800.2012
##
## Tuning parameter 'alpha' was held constant at a value of 0
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were alpha = 0 and lambda = 10.
lassoGrid <- expand.grid(alpha = 1,
lambda = 10^seq(-3, 1, length = 100))
lasso_caret <- train(
Apps ~ .,
data = college_test,
method = "glmnet",
tuneGrid = lassoGrid,
trControl = fitControl,
preProcess = c("center", "scale")
)
## + Fold01.Rep1: alpha=1, lambda=10
## - Fold01.Rep1: alpha=1, lambda=10
## + Fold02.Rep1: alpha=1, lambda=10
## - Fold02.Rep1: alpha=1, lambda=10
## + Fold03.Rep1: alpha=1, lambda=10
## - Fold03.Rep1: alpha=1, lambda=10
## + Fold04.Rep1: alpha=1, lambda=10
## - Fold04.Rep1: alpha=1, lambda=10
## + Fold05.Rep1: alpha=1, lambda=10
## - Fold05.Rep1: alpha=1, lambda=10
## + Fold06.Rep1: alpha=1, lambda=10
## - Fold06.Rep1: alpha=1, lambda=10
## + Fold07.Rep1: alpha=1, lambda=10
## - Fold07.Rep1: alpha=1, lambda=10
## + Fold08.Rep1: alpha=1, lambda=10
## - Fold08.Rep1: alpha=1, lambda=10
## + Fold09.Rep1: alpha=1, lambda=10
## - Fold09.Rep1: alpha=1, lambda=10
## + Fold10.Rep1: alpha=1, lambda=10
## - Fold10.Rep1: alpha=1, lambda=10
## + Fold01.Rep2: alpha=1, lambda=10
## - Fold01.Rep2: alpha=1, lambda=10
## + Fold02.Rep2: alpha=1, lambda=10
## - Fold02.Rep2: alpha=1, lambda=10
## + Fold03.Rep2: alpha=1, lambda=10
## - Fold03.Rep2: alpha=1, lambda=10
## + Fold04.Rep2: alpha=1, lambda=10
## - Fold04.Rep2: alpha=1, lambda=10
## + Fold05.Rep2: alpha=1, lambda=10
## - Fold05.Rep2: alpha=1, lambda=10
## + Fold06.Rep2: alpha=1, lambda=10
## - Fold06.Rep2: alpha=1, lambda=10
## + Fold07.Rep2: alpha=1, lambda=10
## - Fold07.Rep2: alpha=1, lambda=10
## + Fold08.Rep2: alpha=1, lambda=10
## - Fold08.Rep2: alpha=1, lambda=10
## + Fold09.Rep2: alpha=1, lambda=10
## - Fold09.Rep2: alpha=1, lambda=10
## + Fold10.Rep2: alpha=1, lambda=10
## - Fold10.Rep2: alpha=1, lambda=10
## + Fold01.Rep3: alpha=1, lambda=10
## - Fold01.Rep3: alpha=1, lambda=10
## + Fold02.Rep3: alpha=1, lambda=10
## - Fold02.Rep3: alpha=1, lambda=10
## + Fold03.Rep3: alpha=1, lambda=10
## - Fold03.Rep3: alpha=1, lambda=10
## + Fold04.Rep3: alpha=1, lambda=10
## - Fold04.Rep3: alpha=1, lambda=10
## + Fold05.Rep3: alpha=1, lambda=10
## - Fold05.Rep3: alpha=1, lambda=10
## + Fold06.Rep3: alpha=1, lambda=10
## - Fold06.Rep3: alpha=1, lambda=10
## + Fold07.Rep3: alpha=1, lambda=10
## - Fold07.Rep3: alpha=1, lambda=10
## + Fold08.Rep3: alpha=1, lambda=10
## - Fold08.Rep3: alpha=1, lambda=10
## + Fold09.Rep3: alpha=1, lambda=10
## - Fold09.Rep3: alpha=1, lambda=10
## + Fold10.Rep3: alpha=1, lambda=10
## - Fold10.Rep3: alpha=1, lambda=10
## Aggregating results
## Selecting tuning parameters
## Fitting alpha = 1, lambda = 10 on full training set
print(lasso_caret)
## glmnet
##
## 156 samples
## 17 predictor
##
## Pre-processing: centered (17), scaled (17)
## Resampling: Cross-Validated (10 fold, repeated 3 times)
## Summary of sample sizes: 140, 140, 140, 142, 140, 142, ...
## Resampling results across tuning parameters:
##
## lambda RMSE Rsquared MAE
## 0.001000000 1435.880 0.8898335 810.0747
## 0.001097499 1435.880 0.8898335 810.0747
## 0.001204504 1435.880 0.8898335 810.0747
## 0.001321941 1435.880 0.8898335 810.0747
## 0.001450829 1435.880 0.8898335 810.0747
## 0.001592283 1435.880 0.8898335 810.0747
## 0.001747528 1435.880 0.8898335 810.0747
## 0.001917910 1435.880 0.8898335 810.0747
## 0.002104904 1435.880 0.8898335 810.0747
## 0.002310130 1435.880 0.8898335 810.0747
## 0.002535364 1435.880 0.8898335 810.0747
## 0.002782559 1435.880 0.8898335 810.0747
## 0.003053856 1435.880 0.8898335 810.0747
## 0.003351603 1435.880 0.8898335 810.0747
## 0.003678380 1435.880 0.8898335 810.0747
## 0.004037017 1435.880 0.8898335 810.0747
## 0.004430621 1435.880 0.8898335 810.0747
## 0.004862602 1435.880 0.8898335 810.0747
## 0.005336699 1435.880 0.8898335 810.0747
## 0.005857021 1435.880 0.8898335 810.0747
## 0.006428073 1435.880 0.8898335 810.0747
## 0.007054802 1435.880 0.8898335 810.0747
## 0.007742637 1435.880 0.8898335 810.0747
## 0.008497534 1435.880 0.8898335 810.0747
## 0.009326033 1435.880 0.8898335 810.0747
## 0.010235310 1435.880 0.8898335 810.0747
## 0.011233240 1435.880 0.8898335 810.0747
## 0.012328467 1435.880 0.8898335 810.0747
## 0.013530478 1435.880 0.8898335 810.0747
## 0.014849683 1435.880 0.8898335 810.0747
## 0.016297508 1435.880 0.8898335 810.0747
## 0.017886495 1435.880 0.8898335 810.0747
## 0.019630407 1435.880 0.8898335 810.0747
## 0.021544347 1435.880 0.8898335 810.0747
## 0.023644894 1435.880 0.8898335 810.0747
## 0.025950242 1435.880 0.8898335 810.0747
## 0.028480359 1435.880 0.8898335 810.0747
## 0.031257158 1435.880 0.8898335 810.0747
## 0.034304693 1435.880 0.8898335 810.0747
## 0.037649358 1435.880 0.8898335 810.0747
## 0.041320124 1435.880 0.8898335 810.0747
## 0.045348785 1435.880 0.8898335 810.0747
## 0.049770236 1435.880 0.8898335 810.0747
## 0.054622772 1435.880 0.8898335 810.0747
## 0.059948425 1435.880 0.8898335 810.0747
## 0.065793322 1435.880 0.8898335 810.0747
## 0.072208090 1435.880 0.8898335 810.0747
## 0.079248290 1435.880 0.8898335 810.0747
## 0.086974900 1435.880 0.8898335 810.0747
## 0.095454846 1435.880 0.8898335 810.0747
## 0.104761575 1435.880 0.8898335 810.0747
## 0.114975700 1435.880 0.8898335 810.0747
## 0.126185688 1435.880 0.8898335 810.0747
## 0.138488637 1435.880 0.8898335 810.0747
## 0.151991108 1435.880 0.8898335 810.0747
## 0.166810054 1435.880 0.8898335 810.0747
## 0.183073828 1435.880 0.8898335 810.0747
## 0.200923300 1435.880 0.8898335 810.0747
## 0.220513074 1435.880 0.8898335 810.0747
## 0.242012826 1435.880 0.8898335 810.0747
## 0.265608778 1435.880 0.8898335 810.0747
## 0.291505306 1435.880 0.8898335 810.0747
## 0.319926714 1435.880 0.8898335 810.0747
## 0.351119173 1435.880 0.8898335 810.0747
## 0.385352859 1435.880 0.8898335 810.0747
## 0.422924287 1435.880 0.8898335 810.0747
## 0.464158883 1435.880 0.8898335 810.0747
## 0.509413801 1435.880 0.8898335 810.0747
## 0.559081018 1435.880 0.8898335 810.0747
## 0.613590727 1435.880 0.8898335 810.0747
## 0.673415066 1435.880 0.8898335 810.0747
## 0.739072203 1435.880 0.8898335 810.0747
## 0.811130831 1435.880 0.8898335 810.0747
## 0.890215085 1435.880 0.8898335 810.0747
## 0.977009957 1435.880 0.8898335 810.0747
## 1.072267222 1436.200 0.8898322 810.1355
## 1.176811952 1436.765 0.8898299 810.2436
## 1.291549665 1437.383 0.8898272 810.3616
## 1.417474163 1437.907 0.8898399 810.3722
## 1.555676144 1438.444 0.8898592 810.3456
## 1.707352647 1438.904 0.8898813 810.2721
## 1.873817423 1439.174 0.8899838 809.9645
## 2.056512308 1439.330 0.8901385 809.4955
## 2.257019720 1439.551 0.8902953 808.9607
## 2.477076356 1439.770 0.8904724 808.3715
## 2.718588243 1440.048 0.8906695 807.7372
## 2.983647240 1440.370 0.8908791 807.0418
## 3.274549163 1440.783 0.8910966 806.2986
## 3.593813664 1441.293 0.8913173 805.4984
## 3.944206059 1441.782 0.8915507 804.6095
## 4.328761281 1442.168 0.8917991 803.5512
## 4.750810162 1442.287 0.8920643 802.5242
## 5.214008288 1441.667 0.8924080 801.2134
## 5.722367659 1440.314 0.8928492 799.5583
## 6.280291442 1438.681 0.8933613 797.7385
## 6.892612104 1437.016 0.8939027 795.8179
## 7.564633276 1435.171 0.8945058 793.7226
## 8.302175681 1432.811 0.8951903 791.4643
## 9.111627561 1429.734 0.8961292 788.8419
## 10.000000000 1426.235 0.8972572 785.9028
##
## Tuning parameter 'alpha' was held constant at a value of 1
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were alpha = 1 and lambda = 10.
library(pls)
##
## Attaching package: 'pls'
## The following object is masked from 'package:caret':
##
## R2
## The following object is masked from 'package:stats':
##
## loadings
set.seed(42)
pcr_fit <- pcr(Apps ~., data = College, scale = TRUE, validation = "CV")
summary(pcr_fit)
## Data: X dimension: 777 17
## Y dimension: 777 1
## Fit method: svdpc
## Number of components considered: 17
##
## VALIDATION: RMSEP
## Cross-validated using 10 random segments.
## (Intercept) 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps
## CV 3873 3833 2039 2049 1841 1585 1576
## adjCV 3873 3833 2036 2049 1739 1577 1572
## 7 comps 8 comps 9 comps 10 comps 11 comps 12 comps 13 comps
## CV 1565 1533 1491 1491 1497 1497 1500
## adjCV 1563 1527 1489 1488 1494 1494 1497
## 14 comps 15 comps 16 comps 17 comps
## CV 1502 1423 1146 1114
## adjCV 1500 1410 1141 1109
##
## TRAINING: % variance explained
## 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps 7 comps 8 comps
## X 31.670 57.30 64.30 69.90 75.39 80.38 83.99 87.40
## Apps 2.316 73.06 73.07 82.08 84.08 84.11 84.32 85.18
## 9 comps 10 comps 11 comps 12 comps 13 comps 14 comps 15 comps
## X 90.50 92.91 95.01 96.81 97.9 98.75 99.36
## Apps 85.88 86.06 86.06 86.10 86.1 86.13 90.32
## 16 comps 17 comps
## X 99.84 100.00
## Apps 92.52 92.92
validationplot(pcr_fit,val.type="MSEP")
pcr_pred=predict(pcr_fit, newdata = college_test, ncomp = 17)
mean((college_test$Apps - pcr_pred)^2)
## [1] 1307328
set.seed(42)
pls.fit=plsr(Apps ~ ., data = college_train, scale = TRUE, validation ="CV")
summary(pls.fit)
## Data: X dimension: 621 17
## Y dimension: 621 1
## Fit method: kernelpls
## Number of components considered: 17
##
## VALIDATION: RMSEP
## Cross-validated using 10 random segments.
## (Intercept) 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps
## CV 3606 1546 1324 1154 1130 1096 1058
## adjCV 3606 1544 1324 1152 1125 1083 1051
## 7 comps 8 comps 9 comps 10 comps 11 comps 12 comps 13 comps
## CV 1041 1039 1037 1034 1036 1035 1035
## adjCV 1037 1036 1034 1031 1032 1031 1032
## 14 comps 15 comps 16 comps 17 comps
## CV 1035 1035 1035 1035
## adjCV 1032 1032 1032 1032
##
## TRAINING: % variance explained
## 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps 7 comps 8 comps
## X 26.34 46.91 62.99 66.01 68.24 71.98 75.40 81.19
## Apps 82.05 87.06 90.33 91.21 92.20 92.54 92.63 92.64
## 9 comps 10 comps 11 comps 12 comps 13 comps 14 comps 15 comps
## X 83.51 85.23 87.18 88.80 91.34 93.31 97.11
## Apps 92.67 92.69 92.70 92.71 92.71 92.71 92.71
## 16 comps 17 comps
## X 99.34 100.00
## Apps 92.71 92.71
validationplot(pls.fit,val.type="MSEP")
pls.pred=predict(pls.fit, newdata = college_test, ncomp = 17)
mean((college_test$Apps - pls.pred)^2)
## [1] 1941715
To comment on the accuracy we look at the RMSE. We see the RMSE is 1393 which is pretty low compared to universities that receive the higher end of applications. If the school is a small college only receiving ~800 applications, this number is pretty large. We found that OLS, PCR and PLS will have similar test errors. Ridge and lasso regression offer only a slight improvement compared to the OLS model since it introduce regularization to the model.
library(ISLR2)
library(leaps)
## Warning: package 'leaps' was built under R version 4.3.3
data(Boston)
head(Boston)
## crim zn indus chas nox rm age dis rad tax ptratio lstat medv
## 1 0.00632 18 2.31 0 0.538 6.575 65.2 4.0900 1 296 15.3 4.98 24.0
## 2 0.02731 0 7.07 0 0.469 6.421 78.9 4.9671 2 242 17.8 9.14 21.6
## 3 0.02729 0 7.07 0 0.469 7.185 61.1 4.9671 2 242 17.8 4.03 34.7
## 4 0.03237 0 2.18 0 0.458 6.998 45.8 6.0622 3 222 18.7 2.94 33.4
## 5 0.06905 0 2.18 0 0.458 7.147 54.2 6.0622 3 222 18.7 5.33 36.2
## 6 0.02985 0 2.18 0 0.458 6.430 58.7 6.0622 3 222 18.7 5.21 28.7
train_idx <- sample(1:nrow(Boston), 0.7 * nrow(Boston))
train_data <- Boston[train_idx, ]
test_data <- Boston[-train_idx, ]
x.train <- model.matrix(crim ~ ., data = train_data)[, -1]
y.train <- train_data$crim
x.test <- model.matrix(crim ~ ., data = test_data)[, -1]
y.test <- test_data$crim
best_subset <- regsubsets(crim ~ ., data = train_data, nvmax = 13)
subset_summary <- summary(best_subset)
best_size <- which.min(subset_summary$bic)
print(coef(best_subset, best_size))
## (Intercept) rad lstat
## -4.3591153 0.4300705 0.2817968
coef_best <- coef(best_subset, id = best_size)
subset_preds <- model.matrix(crim ~ ., data = test_data)[, names(coef_best)] %*% coef_best
a <- print(mean((y.test - subset_preds)^2))
## [1] 81.54036
cv_ridge_boston <- cv.glmnet(x.train, y.train, alpha = 0)
best_lambda_ridge <- cv_ridge_boston$lambda.min
print(best_lambda_ridge)
## [1] 0.4709069
ridge_preds_boston <- predict(cv_ridge_boston, s = best_lambda_ridge, newx = x.test)
b <- print(mean((y.test - ridge_preds_boston)^2))
## [1] 80.5106
cv_lasso_boston <- cv.glmnet(x.train, y.train, alpha = 1)
best_lambda_lasso <- cv_lasso_boston$lambda.min
lasso_coef <- predict(cv_lasso_boston, s = best_lambda_lasso, type = "coefficients")
print(lasso_coef)
## 13 x 1 sparse Matrix of class "dgCMatrix"
## s1
## (Intercept) 7.910636e+00
## zn 3.306823e-02
## indus -7.704398e-02
## chas -1.081364e+00
## nox -5.285668e+00
## rm -3.064599e-01
## age -4.962980e-03
## dis -6.245458e-01
## rad 4.735539e-01
## tax -3.658554e-05
## ptratio -1.759717e-01
## lstat 2.685792e-01
## medv -4.880134e-02
lasso_preds <- predict(cv_lasso_boston, s = best_lambda_lasso, newx = x.test)
c <- print(mean((y.test - lasso_preds)^2))
## [1] 79.5121
pcr_fit_boston <- pcr(crim ~ ., data = train_data, scale = TRUE, validation = "CV")
summary(pcr_fit_boston)
## Data: X dimension: 354 12
## Y dimension: 354 1
## Fit method: svdpc
## Number of components considered: 12
##
## VALIDATION: RMSEP
## Cross-validated using 10 random segments.
## (Intercept) 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps
## CV 7.224 5.968 5.980 5.618 5.608 5.516 5.494
## adjCV 7.224 5.965 5.976 5.610 5.604 5.510 5.489
## 7 comps 8 comps 9 comps 10 comps 11 comps 12 comps
## CV 5.361 5.362 5.343 5.310 5.327 5.273
## adjCV 5.354 5.355 5.335 5.301 5.317 5.263
##
## TRAINING: % variance explained
## 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps 7 comps 8 comps
## X 49.94 64.05 73.38 80.84 87.28 90.61 93.10 95.10
## crim 32.78 32.97 41.13 41.46 43.58 43.93 47.06 47.67
## 9 comps 10 comps 11 comps 12 comps
## X 96.90 98.46 99.54 100.00
## crim 48.35 49.15 49.32 50.29
best_msep <- which.min(MSEP(pcr_fit_boston)$val[1,1,-1])
pcr.preds <- predict(pcr_fit_boston, test_data, ncomp = best_msep)
d <- print(mean((y.test - pcr.preds)^2))
## [1] 79.11225
print(c(a,b,c,d))
## [1] 81.54036 80.51060 79.51210 79.11225
I propose the Lasso Regression model as the optimal choice. It utilizes 10-fold cross-validation on the training data to dictate the bias-variance tradeoff via lambda. When tested on a completely independent validation set, it consistently yields the lowest Test MSE by eliminating unhelpful noise variables.
The chosen models does not involve all of the data. The lasso model performs automatic variable selection. Since lambda increase, penalty forces coefficient estimates of lesser variables to become 0.