Linear Model Selection and Regularization
Justin Gould
2026-07-06
Introduction This lab will explore regularization techniques such as OLS, Lasso and Ridge regression, Subset Selection, Principal Componant Regression, and Partial Least Squares. This lab explores questions 2, 9, and 11 from the ISLR textbook applying data from the College and Boston dataset from the ISLR2 library.
Limitations:
- On both problems 9 and 11 for the Ordinary Least Squares and Ridge sections you will notice some of the predictions are negative. This is a limitation of linear regularization methods - these regularization methods do not account for the fact that some variables cannot be negative. For example, you cannot have negative college applications or negative crime rates. The amount of negative predictions that are produces within these methods should be gaged to measure the severity of the limitation. OLS and Ridge regression are most likely to produce negative predictions because they don’t penalize coeffecients.
2 Comparing Methods
For parts (a) through (c), indicate which of i. through iv. is correct. Justify your answer.
(a) Lasso vs OLS
The lasso, relative to least squares, is:
- Less flexible and hence will give improved prediction accuracy when its increase in bias is less than its decrease in variance.
True Lasso regression is a less flexible model than OLS because it restrains the coeffecients unlike OLS. As the coeffiencents become more restrained, bias decreases because the assumptions are more complex and variance increases.
(b) Ridge vs OLS
Repeat (a) for ridge regression relative to least squares.
- Less flexible and hence will give improved prediction accuracy when its increase in bias is less than its decrease in variance.
True Ridge regression is a less flexible model because it restrains the coefficients unlike OLS. As the coefficients become more restrained, bias decreases because the assumptions are more complex and variance increases.
(c) Non-linear methods vs OLS
Repeat (a) for non-linear methods relative to least squares. i. More flexible and hence will give improved prediction accuracy when its increase in bias is less than its decrease in variance. True Non-linear methods tend to be more flexible than linear methods. This flexibility causes a more accurate prediction to occur when the bias is greatest than the variance Due to its flexibility, the non linear method would have more complex assumptions.
9. College dataset
In this exercise, we will predict the number of applications received
using the other variables in the College data set.
Collegeattach(College)(a) Test/Train split
Split the data set into a training set and a test set.
# Step 1. Setting the random sample at 1, anyone who runs this code will get the same random sample.
set.seed(1)
# Step 2. Create a 30/70 split
train <- sample(777, 0.7 * 777)
# Step 3. Training data: 70% of the total observations
train_data <- College[train, ]
# Step 4. Test data: 30% of total
test_data <- College[-train, ](b) OLS model
Fit a linear model using least squares on the training set, and report the test error obtained.
# Step 1. Fit the linear model using OLS
lm_fit <- lm(Apps ~ ., data = train_data)
# Step 2. Predict how many application each college have
prediction <- predict(lm_fit, newdata = test_data)
prediction |> head(20)## Albright College Alderson-Broaddus College
## 933.02355 601.08836
## Allegheny College Allentown Coll. of St. Francis de Sales
## 2968.77127 1815.84026
## Alma College Amherst College
## 1939.42215 3802.65706
## Aquinas College Arizona State University Main campus
## 289.25505 15319.12774
## Arkansas Tech University Augustana College IL
## 2040.11507 2380.98350
## Barat College Barnard College
## 130.38897 2895.20627
## Bemidji State University Benedictine College
## 1181.71002 478.39857
## Berry College Bethel College KS
## 2080.75592 -38.70755
## Bethel College Bloomsburg Univ. of Pennsylvania
## 39.37242 4296.55287
## Bluffton College Bradford College
## 250.49734 -138.65982# Step 3. Calculate MSE
error_lm <- mean((test_data$Apps - prediction)^2)
error_lm## [1] 1261630OLS Test MSE 1,261,630
# Due Diligence: Determine how many of the predictions are negative.
sum(prediction < 0)## [1] 13sum(prediction > 0)## [1] 221Limitation: When performed on the College test dataset OLS produces 13 negative predictions and 221 positive. Since this is only a small fraction of the total this method should not be ruled out entirely.
(c) Ridge model
Fit a ridge regression model on the training set, with \(\lambda\) chosen by cross-validation. Report the test error obtained.
# Step 1. Define x (the predictors)
x = model.matrix(Apps ~ ., data = train_data)[, -1]
# Step 2. Define y (the response)
y = train_data$Apps
# Step 3. Fit the ridge model
ridge_fit <- cv.glmnet(x, y, alpha =0 )
# Step 4. Run the dimensions of the model.
dim(coef(ridge_fit))## [1] 18 1# Evaluate Lambda on training data
ridge_fit$lambda.min## [1] 367.5286# Step 5. Calculate test x and y
x_test = model.matrix(Apps ~ ., data = test_data)[, -1]
y_test = test_data$Apps
# Step 6. Predict on the test set
prediction_ridge <- predict(ridge_fit, newx = x_test)
# Step 7. Calculate test MSE
error_ridge <- mean((test_data$Apps - prediction_ridge)^2)
error_ridge## [1] 1453474# Due Diligence: Determine how many of the predictions are negative.
sum(prediction_ridge < 0)## [1] 3sum(prediction_ridge > 0)## [1] 231Ridge Test MSE 1,453,474. With only 3 negative predictions the Ridge method is a strong regularization method.
(d) Lasso model
Fit a lasso model on the training set, with \(\lambda\) chosen by cross-validation. Report the test error obtained, along with the number of non-zero coefficient estimates.
# Step 1. Define x (the predictors)
x = model.matrix(Apps ~ ., data = train_data)[, -1]
# Step 2. Define y (the response)
y = train_data$Apps
# Step 3. Fit the lasso model
lasso_fit <- cv.glmnet(x, y, alpha =1 )
# Step 4. Run the dimensions of the model.
dim(coef(lasso_fit))## [1] 18 1# Evaluate Lambda on training data
lasso_fit$lambda.min## [1] 1.96139# Step 5. Calculate test x and y
x_test = model.matrix(Apps ~ ., data = test_data)[, -1]
y_test = test_data$Apps
# Step 6. Predict on the test set
prediction_lasso <- predict(lasso_fit, newx = x_test)
# Step 7. Calculate test MSE
error_lasso <- mean((test_data$Apps - prediction_lasso)^2)
error_lasso## [1] 1314717# Due Diligence: Determine how many of the predictions are negative.
sum(prediction_lasso < 0)## [1] 0sum(prediction_lasso > 0)## [1] 234Lasso Test MSE 1,314,717. With 0 negative predictions Lasso is a great reguarization method.
(e) PCR model
Fit a PCR model on the training set, with \(M\) chosen by cross-validation. Report the test error obtained, along with the value of \(M\) selected by cross-validation.
# Step 1. Set seed for the CV
set.seed(2)
# Step 2. Fit the PCR model
pcr.fit = pcr(Apps ~ .,
data = train_data,
scale = TRUE,
validation = "CV")
summary(pcr.fit)## Data: X dimension: 543 17
## Y dimension: 543 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 3895 3811 2119 2127 1799 1704 1707
## adjCV 3895 3813 2116 2124 1779 1696 1702
## 7 comps 8 comps 9 comps 10 comps 11 comps 12 comps 13 comps
## CV 1707 1660 1613 1609 1613 1613 1616
## adjCV 1704 1649 1607 1603 1607 1608 1611
## 14 comps 15 comps 16 comps 17 comps
## CV 1617 1549 1154 1120
## adjCV 1612 1531 1145 1113
##
## TRAINING: % variance explained
## 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps 7 comps 8 comps
## X 32.051 57.00 64.42 70.27 75.65 80.65 84.26 87.61
## Apps 5.788 71.69 71.70 80.97 82.60 82.60 82.69 84.06
## 9 comps 10 comps 11 comps 12 comps 13 comps 14 comps 15 comps
## X 90.58 92.84 94.93 96.74 97.82 98.72 99.39
## Apps 84.55 84.82 84.86 84.86 85.01 85.05 89.81
## 16 comps 17 comps
## X 99.85 100.00
## Apps 93.03 93.32# Step 3. Principal component graph
validationplot(pcr.fit, val.type = "MSEP")
prediction_pcr <- predict(pcr.fit, newdata = test_data, ncomp = 17)
# Step 7. Calculate test MSE
error_pcr <- mean((test_data$Apps - prediction_pcr)^2)
error_pcr## [1] 1261630# Due Diligence: Determine how many of the predictions are negative.
sum(prediction_pcr < 0)## [1] 13sum(prediction_pcr > 0)## [1] 221PCR Test MSE 1,261,630, PCR produced 13 negative values, identifcal to OLS.
(f) PLS model
Fit a PLS model on the training set, with \(M\) chosen by cross-validation. Report the test error obtained, along with the value of \(M\) selected by cross-validation.
# Step 1. Set seed for the CV
set.seed(3)
# Step 2. Fit the PLS model
pls.fit = plsr(Apps ~ .,
data = train_data,
scale = TRUE,
validation = "CV")
summary(pls.fit)## Data: X dimension: 543 17
## Y dimension: 543 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 3895 1942 1748 1533 1454 1279 1189
## adjCV 3895 1937 1746 1525 1428 1244 1177
## 7 comps 8 comps 9 comps 10 comps 11 comps 12 comps 13 comps
## CV 1172 1162 1153 1154 1151 1151 1151
## adjCV 1162 1153 1144 1145 1142 1142 1142
## 14 comps 15 comps 16 comps 17 comps
## CV 1150 1150 1150 1150
## adjCV 1141 1141 1141 1141
##
## TRAINING: % variance explained
## 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps 7 comps 8 comps
## X 25.68 47.43 62.46 64.88 67.34 72.68 77.20 80.92
## Apps 76.62 82.39 86.93 90.76 92.82 93.05 93.13 93.20
## 9 comps 10 comps 11 comps 12 comps 13 comps 14 comps 15 comps
## X 82.69 85.16 87.35 90.73 92.49 95.10 97.09
## Apps 93.26 93.28 93.30 93.31 93.32 93.32 93.32
## 16 comps 17 comps
## X 98.40 100.00
## Apps 93.32 93.32# Step 3. Principal component graph
validationplot(pls.fit, val.type = "MSEP")
prediction_pls <- predict(pls.fit, newdata = test_data, ncomp = 14)
# Step 7. Calculate test MSE
error_pls <- mean((test_data$Apps - prediction_pls)^2)
error_pls## [1] 1262592# Due Diligence: Determine how many of the predictions are negative.
sum(prediction_pls< 0)## [1] 13sum(prediction_pls > 0)## [1] 221PLS Test MSE 1,262,592. Identifcal to PCR and OLS this method has a 13/221 split.
(g) Conclusion
Comment on the results obtained. How accurately can we predict the number of college applications received? Is there much difference among the test errors resulting from these five approaches?
Conclusion: The lowest test MSE of the five models was Ordinary Least Squares containing a value of 1,261,630.41. This indicates that OLS was the best regularization method for predicting applications on the College test data. By taking the square root of the MSE we get a root mean squared error (RMSE) of 1123.223. This value shows that a OLS model’s predictions would be off by around 1,123.22 college applications on average when making predictions. If prescribing a regularization method with no negative values it the goal than Lasso would be the best approach.
11 Boston dataset
We will now try to predict per capita crime rate in the
Boston data set.
Bostonattach(Boston)(a) Exploring Regularization Methods
In this section I explore OLS and Stepwise Selection to select the best model.
i. Ordinary Least Squares
Try out some of the regression methods explored in this chapter, such as best subset selection, the lasso, ridge regression, and PCR. Present and discuss results for the approaches that you consider.
# Step 1. Setting the random sample at 1, anyone who runs this code will get the same random sample.
set.seed(5)
# Step 2. Create a 30/70 split
training <- sample(506, 0.7 * 506)
# Step 3. Training data: 70% of the total observations
boston_train <- Boston[training, ]
# Step 4. Test data: 30% of total
boston_test <- Boston[-train, ]# Step 1. Fit the linear model using OLS
boston_lm <- lm(crim ~ ., data = boston_train)
# Step 2. Predict how many application each college have
boston_prediction <- predict(boston_lm, newdata = boston_test)
boston_prediction |> head(20)## 9 10 12 13 14 17 23
## 3.6607881 1.1798532 1.1704771 1.7674738 -1.8887969 -2.4532991 1.9532504
## 24 26 30 36 38 46 47
## 1.7289145 0.8759291 -0.3390880 2.5550312 1.1728961 0.2009310 0.6200915
## 50 52 53 57 59 63
## -0.5752299 0.3007753 -0.8462041 -1.8123780 -0.4158228 0.7350519# Step 3. Calculate MSE
error_lm <- mean((test_data$Apps - boston_prediction)^2)
error_lm## [1] 24184666# Due Diligence: Determine how many of the predictions are negative.
sum(boston_prediction < 0)## [1] 43sum(boston_prediction > 0)## [1] 110Takeaway: Despite the MSE of the OLS model the limitations prove that this model should NOT be used. With 43 out of 153 negative values this OLS model is predicting far to many negative values for this to be a trusted method.
ii. Subset Selection
# Step 1: Run subset selection on Boston dataset
subset <- regsubsets(crim ~ ., Boston, nvmax = 12)
# Step 2: View Results
reg.summary <- summary(subset)
reg.summary$rsq## [1] 0.3912567 0.4207965 0.4244920 0.4334892 0.4378328 0.4421077 0.4451514
## [8] 0.4470236 0.4482891 0.4487917 0.4493353 0.4493378# Step 3: Plot the results
plot(reg.summary$rsq,
main = "Variable increase vs R squared",
xlab = "Number of Variables",
ylab = "R squared",
type = "o",
col = "darkblue",
las = 1)
Takeaway: This graph displays the increase in R squared while adding variables. Specifically, when variables change from 1-7 we see the highest increase in R squared and from 8-12 there is minor increase. Next, we examine which variables specifically will have the highest impact on the response.
# Step 1: View the subset selection plot to see which variables are impactfull
reg.summary$outmat## zn indus chas nox rm age dis rad tax ptratio lstat medv
## 1 ( 1 ) " " " " " " " " " " " " " " "*" " " " " " " " "
## 2 ( 1 ) " " " " " " " " " " " " " " "*" " " " " "*" " "
## 3 ( 1 ) " " " " " " " " " " " " " " "*" " " " " "*" "*"
## 4 ( 1 ) "*" " " " " " " " " " " "*" "*" " " " " " " "*"
## 5 ( 1 ) "*" "*" " " " " " " " " "*" "*" " " " " " " "*"
## 6 ( 1 ) "*" " " " " "*" " " " " "*" "*" " " "*" " " "*"
## 7 ( 1 ) "*" " " " " "*" " " " " "*" "*" " " "*" "*" "*"
## 8 ( 1 ) "*" "*" " " "*" " " " " "*" "*" " " "*" "*" "*"
## 9 ( 1 ) "*" "*" " " "*" "*" " " "*" "*" " " "*" "*" "*"
## 10 ( 1 ) "*" " " "*" "*" "*" " " "*" "*" "*" "*" "*" "*"
## 11 ( 1 ) "*" "*" "*" "*" "*" " " "*" "*" "*" "*" "*" "*"
## 12 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*" "*" "*" "*" "*"# Step 2: Plot this on a graph (AI was used on this step)
out <- as.data.frame(reg.summary$outmat)
out$model <- 1:nrow(out)
out_long <- pivot_longer(
out,
cols = -model,
names_to = "Variable",
values_to = "Included"
)
out_long$Included <- out_long$Included == "*"
ggplot(out_long,
aes(x = Variable,
y = factor(model),
fill = Included)) +
geom_tile(color = "white") +
scale_fill_manual(values = c("lightgrey", "steelblue")) +
labs(title = "Best Subset Selection",
x = "Predictor",
y = "Model Size") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
Takeaway: The variables that appear in almost every model are Rad, Zn, Medv, and Dis. This shows that these predictors have a major impact on the Crime Rates per capita in Boston and should be considered in the model. Next, I use K-fold cross validation to determine which subset to select to form the best model.
(b) Cross Validation on Subsets
Propose a model (or set of models) that seem to perform well on this data set, and justify your answer. Make sure that you are evaluating model performance using validation set error, cross-validation, or some other reasonable alternative, as opposed to using training error.
# Prediction Function
predict.regsubsets <- function(object, newdata, id, ...) {
form <- as.formula(object$call[[2]])
mat <- model.matrix(form, newdata)
coefi <- coef(object, id = id)
xvars <- names(coefi)
mat[, xvars] %*% coefi
}
# Step 1: Create a matrix to store the result of each fold and CV
k = 10
set.seed(1)
folds = sample(1:k,
nrow(Boston),
replace = TRUE)
cv.error = matrix(NA, k, 12)
# Step 2. Write a for-loop to perform cross-validation
for(j in 1:k){
best.fit = regsubsets(crim ~ .,
data = Boston[folds != j, ],
nvmax = 12)
for(i in 1:12){
pred = predict(best.fit,
Boston[folds == j, ],
id = i)
cv.error[j, i] =
mean((Boston$crim[folds == j] - pred)^2)
}
}
# Step 3. Average the CV errors
mean.cv.error = apply(cv.error, 2, mean)plot(mean.cv.error,
main = "Mean CV Error",
xlab = "Number of Predictors",
ylab = "Mean CV MSE",
type = "o",
col = "darkblue",
las = 1)
Interpretation: An 11 variable model would perform the best on this dataset. The MSE of 42.39886 plotted on the graph should that his model would perform the best between all the subsets. Next, I determine which 11 variables to use
# Add the predictors to a data table
coef_table <- data.frame(
Predictor = names(coef(best.fit, id = 11))[-1],
Coefficient = round(as.numeric(coef(best.fit, id = 11))[-1], 4)
)
datatable(
coef_table,
caption = "Coefficients for the Best 11-Variable Model",
rownames = FALSE,
options = list(
dom = "t",
pageLength = 11
)
)Interpretation: The 11 Variables to predict crime rate are zn, indus, chas, nox, rm, dis, rad, tax, ptratio, lstat, and medv. Exactly as I suspected earlier, the variable that only showed up in 1 model (age) was NOT a useful factor in predicting crime rates per capita.
(c) Conclusion
Does your chosen model involve all of the features in the data set? Why or why not?
Response: No, the model does not contain all the variables from the Boston dataset. Age was ruled out as an insignificant predictor in the final model after using cross-validation to select the best subset of predictors. The boxchart showed the age was only a useful predictor in the 12 variable model so it is no surprise that it does not show up in the final result.