STA 6543 Predictive Modeling
Homework 5
Haniya Farhat
July 16, 2026
Problem 2: Ridge, Lasso, and Non-linear Methods vs. Least Squares
(a) The lasso, relative to least squares, is:
iii. Less flexible and hence will give improved prediction accuracy when its increase in bias is less than its decrease in variance.
The lasso constrains the size of the coefficients by adding an \(\ell_1\) penalty (\(\sum |\beta_j| \le s\)), which shrinks coefficients toward zero and can set some of them exactly to zero. That’s a less flexible fit than ordinary least squares, which places no constraint on the coefficients at all. Because it’s less flexible, the lasso trades a bit of bias for a often-large reduction in variance, and it beats least squares on prediction accuracy exactly when that variance reduction outweighs the added bias — which tends to happen when \(p\) is large relative to \(n\) or predictors are correlated.
(b) Ridge regression relative to least squares is:
iii. Less flexible and hence will give improved prediction accuracy when its increase in bias is less than its decrease in variance.
Same logic as the lasso: ridge regression adds an \(\ell_2\) penalty (\(\sum \beta_j^2 \le s\)) that shrinks the coefficient estimates toward zero, making it less flexible than least squares. As the penalty (via \(\lambda\)) increases, variance drops off quickly while bias increases slowly, so ridge improves on least squares whenever that variance reduction dominates. The difference from the lasso is only that ridge shrinks coefficients smoothly without ever forcing them to exactly zero.
(c) Non-linear methods relative to least squares:
ii. More flexible and hence will give improved prediction accuracy when its increase in variance is less than its decrease in bias.
Non-linear methods can adapt to curvature in the true relationship between predictors and response, so they are strictly more flexible than a linear least squares fit. That extra flexibility lets them fit the training data more closely (lower bias) but makes the fitted function more sensitive to the particular training sample (higher variance). They improve on least squares when the true relationship really is non-linear enough that the bias reduction from the extra flexibility outweighs the added variance.
Problem 9: Predicting College Applications
data(College, package = "ISLR2")
cat("Dimensions:", nrow(College), "rows x", ncol(College), "columns\n")## Dimensions: 777 rows x 18 columns## 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(a) Split the data into training and test sets
set.seed(1)
n_col <- nrow(College)
train9 <- sample(n_col, n_col * 0.7)
col_tr <- College[train9, ]
col_te <- College[-train9, ]
cat("Training:", nrow(col_tr), " Test:", nrow(col_te), "\n")## Training: 543 Test: 234(b) Linear model fit by least squares
lm_fit <- lm(Apps ~ ., data = col_tr)
lm_pred <- predict(lm_fit, col_te)
lm_mse <- mean((lm_pred - col_te$Apps)^2)
cat("Least squares test MSE:", round(lm_mse, 2), "\n")## Least squares test MSE: 1261630The least squares test MSE is 1.26163^{6} (RMSE of about 1123 applications).
(c) Ridge regression
x_tr9 <- model.matrix(Apps ~ ., col_tr)[, -1]
y_tr9 <- col_tr$Apps
x_te9 <- model.matrix(Apps ~ ., col_te)[, -1]
y_te9 <- col_te$Apps
set.seed(1)
cv_ridge <- cv.glmnet(x_tr9, y_tr9, alpha = 0)
lam_ridge <- cv_ridge$lambda.min
cat("Ridge lambda chosen by CV:", round(lam_ridge, 3), "\n")## Ridge lambda chosen by CV: 367.529ridge_fit <- glmnet(x_tr9, y_tr9, alpha = 0, lambda = lam_ridge)
ridge_pred <- predict(ridge_fit, s = lam_ridge, newx = x_te9)
ridge_mse <- mean((ridge_pred - y_te9)^2)
cat("Ridge test MSE:", round(ridge_mse, 2), "\n")## Ridge test MSE: 1121808Cross-validation picks \(\lambda =\) 367.53, giving a test MSE of 1.121808^{6} — a clear improvement over least squares.
(d) Lasso
set.seed(1)
cv_lasso <- cv.glmnet(x_tr9, y_tr9, alpha = 1)
lam_lasso <- cv_lasso$lambda.min
cat("Lasso lambda chosen by CV:", round(lam_lasso, 3), "\n")## Lasso lambda chosen by CV: 8.69lasso_fit <- glmnet(x_tr9, y_tr9, alpha = 1, lambda = lam_lasso)
lasso_pred <- predict(lasso_fit, s = lam_lasso, newx = x_te9)
lasso_mse <- mean((lasso_pred - y_te9)^2)
cat("Lasso test MSE:", round(lasso_mse, 2), "\n")## Lasso test MSE: 1232781lasso_coef <- predict(lasso_fit, s = lam_lasso, type = "coefficients")
n_nonzero9 <- sum(lasso_coef[-1, 1] != 0)
print(lasso_coef)## 18 x 1 sparse Matrix of class "dgCMatrix"
## s=8.690175
## (Intercept) -587.60604731
## PrivateYes -467.57943263
## Accept 1.66583967
## Enroll -0.73836466
## Top10perc 47.66628066
## Top25perc -11.91753233
## F.Undergrad .
## P.Undergrad 0.06044763
## Outstate -0.07706177
## Room.Board 0.15118171
## Books 0.22000324
## Personal .
## PhD -8.47525968
## Terminal -0.16546727
## S.F.Ratio 11.14075510
## perc.alumni .
## Expend 0.05730903
## Grad.Rate 6.19785260cat("\nNumber of non-zero coefficients (excluding intercept):", n_nonzero9,
"out of", nrow(lasso_coef) - 1, "\n")##
## Number of non-zero coefficients (excluding intercept): 14 out of 17The lasso’s test MSE is 1.232781^{6}, and it keeps
14 of the 17 predictors (it zeroes out F.Undergrad,
Personal, and perc.alumni). It lands between
least squares and ridge here — it does some real variable selection, but
at a small cost in test error relative to ridge.
(e) Principal components regression (PCR)
set.seed(1)
pcr_fit9 <- pcr(Apps ~ ., data = col_tr, scale = TRUE, validation = "CV")
validationplot(pcr_fit9, val.type = "MSEP", main = "PCR: 10-fold CV MSEP by Number of Components")
val_pcr9 <- MSEP(pcr_fit9)
best_M_pcr <- which.min(val_pcr9$val[1, 1, ]) - 1
cat("Number of components chosen by CV:", best_M_pcr, "\n")## Number of components chosen by CV: 17pcr_pred9 <- predict(pcr_fit9, col_te, ncomp = best_M_pcr)
pcr_mse9 <- mean((pcr_pred9 - col_te$Apps)^2)
cat("PCR test MSE:", round(pcr_mse9, 2), "\n")## PCR test MSE: 1261630Cross-validation error keeps dropping as more components are added and never turns back up, so the CV-optimal choice is \(M =\) 17, i.e. essentially the full set of 17 predictors. With all components retained, PCR reduces to ordinary least squares in a rotated coordinate system, so its test MSE (1.26163^{6}) matches part (b) exactly. Looking at the validation plot, though, the CV curve is already close to flat by around 8-9 components — most of the predictive signal is captured well before the technical minimum, it just never quite stops improving enough to justify trimming components on this particular split.
(f) Partial least squares (PLS)
set.seed(1)
pls_fit9 <- plsr(Apps ~ ., data = col_tr, scale = TRUE, validation = "CV")
validationplot(pls_fit9, val.type = "MSEP", main = "PLS: 10-fold CV MSEP by Number of Components")
val_pls9 <- MSEP(pls_fit9)
best_M_pls <- which.min(val_pls9$val[1, 1, ]) - 1
cat("Number of components chosen by CV:", best_M_pls, "\n")## Number of components chosen by CV: 17pls_pred9 <- predict(pls_fit9, col_te, ncomp = best_M_pls)
pls_mse9 <- mean((pls_pred9 - col_te$Apps)^2)
cat("PLS test MSE:", round(pls_mse9, 2), "\n")## PLS test MSE: 1261630PLS behaves the same way here: CV selects \(M =\) 17 components, so it also reduces to the least-squares fit and produces the identical test MSE of 1.26163^{6}. Unlike PCR, PLS builds its components using the response as well as the predictors, so in principle it should reach a low error with fewer components — and looking at the validation curve, it does flatten out slightly earlier than PCR’s, even though the strict minimum still falls at the maximum number of components on this split.
(g) Comparing the five approaches
results9 <- data.frame(
Method = c("Least Squares", "Ridge", "Lasso", "PCR", "PLS"),
Test_MSE = round(c(lm_mse, ridge_mse, lasso_mse, pcr_mse9, pls_mse9), 1)
)
tss9 <- sum((col_te$Apps - mean(col_te$Apps))^2)
results9$Test_R2 <- round(1 - results9$Test_MSE * nrow(col_te) / tss9, 4)
print(results9)## Method Test_MSE Test_R2
## 1 Least Squares 1261630 0.9134
## 2 Ridge 1121808 0.9230
## 3 Lasso 1232781 0.9154
## 4 PCR 1261630 0.9134
## 5 PLS 1261630 0.9134All five methods land in a fairly tight band, explaining roughly
91-92% of the variance in Apps on the test set
(91.3%-92.3%). Ridge does the best here, followed by lasso, with least
squares, PCR, and PLS essentially tied since the latter two ended up
equivalent to least squares once CV kept all their components. The
takeaway is that Apps is quite predictable from the other
variables in College — the number of applications a school
receives is strongly and fairly linearly tied to things like
acceptances, enrollment, and cost, so a simple least-squares fit is
already close to as good as it gets, and the shrinkage methods only buy
a modest improvement by reining in a few noisy coefficients.
Problem 11: Predicting Per-Capita Crime Rate in the Boston Data Set
data(Boston, package = "ISLR2")
cat("Dimensions:", nrow(Boston), "rows x", ncol(Boston), "columns\n")## Dimensions: 506 rows x 13 columns## 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(a) Fitting the models
I’ll hold out a test set, then use 10-fold cross-validation on the training data to tune each method (best subset size, ridge/lasso \(\lambda\), and PCR’s number of components), and finally check each tuned model against the untouched test set.
set.seed(1)
n_bos <- nrow(Boston)
train11 <- sample(n_bos, n_bos * 0.7)
bos_tr <- Boston[train11, ]
bos_te <- Boston[-train11, ]
cat("Training:", nrow(bos_tr), " Test:", nrow(bos_te), "\n")## Training: 354 Test: 152Best subset selection, choosing model size by 10-fold CV on the training set (following the approach from Section 6.5.3 of the text):
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
}
p_bos <- ncol(Boston) - 1
k11 <- 10
set.seed(1)
folds11 <- sample(1:k11, nrow(bos_tr), replace = TRUE)
cv_errors11 <- matrix(NA, k11, p_bos, dimnames = list(NULL, paste(1:p_bos)))
for (j in 1:k11) {
fit_j <- regsubsets(crim ~ ., data = bos_tr[folds11 != j, ], nvmax = p_bos)
for (i in 1:p_bos) {
pred_j <- predict.regsubsets(fit_j, bos_tr[folds11 == j, ], id = i)
cv_errors11[j, i] <- mean((bos_tr$crim[folds11 == j] - pred_j)^2)
}
}
mean_cv11 <- apply(cv_errors11, 2, mean)
plot(1:p_bos, mean_cv11, type = "b", pch = 19,
xlab = "Number of Predictors", ylab = "10-fold CV MSE",
main = "Best Subset Selection: CV Error by Model Size")
## Best subset size chosen by CV: 3subset_full <- regsubsets(crim ~ ., data = bos_tr, nvmax = p_bos)
coef_subset <- coef(subset_full, best_size11)
print(coef_subset)## (Intercept) zn rad lstat
## -5.34841662 0.02589714 0.49660246 0.30694586x_te_subset <- model.matrix(crim ~ ., bos_te)
subset_pred <- x_te_subset[, names(coef_subset)] %*% coef_subset
subset_mse11 <- mean((subset_pred - bos_te$crim)^2)
cat("Best subset test MSE:", round(subset_mse11, 3), "\n")## Best subset test MSE: 60.847Ridge regression:
x_tr11 <- model.matrix(crim ~ ., bos_tr)[, -1]
y_tr11 <- bos_tr$crim
x_te11 <- model.matrix(crim ~ ., bos_te)[, -1]
y_te11 <- bos_te$crim
set.seed(1)
cv_ridge11 <- cv.glmnet(x_tr11, y_tr11, alpha = 0)
lam_ridge11 <- cv_ridge11$lambda.min
ridge_fit11 <- glmnet(x_tr11, y_tr11, alpha = 0, lambda = lam_ridge11)
ridge_pred11 <- predict(ridge_fit11, s = lam_ridge11, newx = x_te11)
ridge_mse11 <- mean((ridge_pred11 - y_te11)^2)
cat("Ridge lambda (CV):", round(lam_ridge11, 3), "\n")## Ridge lambda (CV): 0.519## Ridge test MSE : 58.746Lasso:
set.seed(1)
cv_lasso11 <- cv.glmnet(x_tr11, y_tr11, alpha = 1)
lam_lasso11 <- cv_lasso11$lambda.min
lasso_fit11 <- glmnet(x_tr11, y_tr11, alpha = 1, lambda = lam_lasso11)
lasso_pred11 <- predict(lasso_fit11, s = lam_lasso11, newx = x_te11)
lasso_mse11 <- mean((lasso_pred11 - y_te11)^2)
cat("Lasso lambda (CV):", round(lam_lasso11, 4), "\n")## Lasso lambda (CV): 0.0284## Lasso test MSE : 58.026## 13 x 1 sparse Matrix of class "dgCMatrix"
## s=0.02835162
## (Intercept) 7.67712712
## zn 0.03229332
## indus -0.07627009
## chas -0.65254550
## nox -5.07810152
## rm 0.33334857
## age .
## dis -0.59808124
## rad 0.52240377
## tax .
## ptratio -0.29555492
## lstat 0.22551015
## medv -0.13005115n_nonzero11 <- sum(lasso_coef11[-1, 1] != 0)
cat("\nNon-zero coefficients:", n_nonzero11, "out of", nrow(lasso_coef11) - 1, "\n")##
## Non-zero coefficients: 10 out of 12Principal components regression:
set.seed(1)
pcr_fit11 <- pcr(crim ~ ., data = bos_tr, scale = TRUE, validation = "CV")
validationplot(pcr_fit11, val.type = "MSEP", main = "PCR: 10-fold CV MSEP by Number of Components")
val_pcr11 <- MSEP(pcr_fit11)
best_M_pcr11 <- which.min(val_pcr11$val[1, 1, ]) - 1
cat("Number of components chosen by CV:", best_M_pcr11, "\n")## Number of components chosen by CV: 12pcr_pred11 <- predict(pcr_fit11, bos_te, ncomp = best_M_pcr11)
pcr_mse11 <- mean((pcr_pred11 - bos_te$crim)^2)
cat("PCR test MSE:", round(pcr_mse11, 3), "\n")## PCR test MSE: 57.613Comparison:
results11 <- data.frame(
Method = c("Best Subset", "Ridge", "Lasso", "PCR"),
Predictors = c(best_size11, p_bos, n_nonzero11, best_M_pcr11),
Test_MSE = round(c(subset_mse11, ridge_mse11, lasso_mse11, pcr_mse11), 3)
)
print(results11)## Method Predictors Test_MSE
## 1 Best Subset 3 60.847
## 2 Ridge 12 58.746
## 3 Lasso 10 58.026
## 4 PCR 12 57.613All four methods land in a similar range, roughly 57.6 to 60.8 test
MSE. PCR edges out the rest slightly, with lasso close behind, and both
beat plain best subset selection. Ridge and PCR use all 12 predictors
(PCR through its components rather than the original variables
directly), while the lasso does real variable selection and drops down
to 10 predictors (removing age and tax), and
CV-tuned best subset goes even further, keeping only 3 (zn, rad, lstat).
Shrinking or restricting the model doesn’t hurt performance relative to
using every predictor — a sign that a good chunk of the variables in
Boston carry limited independent information about
crim once the strongest ones (rad,
nox, lstat, and similar) are already in the
model.
(b) Which model would I propose?
I’d propose the lasso model. It essentially ties PCR for the lowest test MSE (58.03 vs. 57.61), all of it measured on a held-out test set rather than training error, so it’s not just an artifact of overfitting. Compared to PCR, the lasso has a real practical advantage: PCR’s components are linear combinations of all the original predictors, so the fitted model isn’t directly interpretable in terms of the original variables, whereas the lasso model uses a genuine subset of the original predictors with directly interpretable coefficients. Compared to the CV-selected best subset model, the lasso achieves meaningfully lower test error (58.03 vs. 60.85) while still cutting the predictor set down noticeably (from 12 to 10). So the lasso gives close to the best predictive accuracy of any method tried, along with a sparse, interpretable set of predictors — the best balance of the two goals.
(c) Does the chosen model involve all of the features?
No. The lasso model uses 10 of the 12 predictors in
Boston, dropping age (proportion of
owner-occupied units built before 1940) and tax (property
tax rate) entirely — their coefficients were shrunk all the way to zero
by the \(\ell_1\) penalty. That’s
exactly the behavior that makes the lasso attractive over ridge or PCR
here: age and tax apparently don’t carry
independent predictive value for crim once variables like
rad (highway accessibility), nox (pollution),
dis (distance to employment centers), and
lstat (socioeconomic status) are already in the model —
likely because they’re correlated with those stronger predictors rather
than because they’re unrelated to crime on their own. Automatically
zeroing out that redundancy is exactly what the lasso is designed to do,
and it’s the reason its fitted model ends up both simpler and slightly
more accurate on the test set than a model that keeps every
predictor.