HW5
NastaranTaabodi -dew050
2026-07-21
For parts (a) through (c), indicate which of i. through iv. is correct.Justify your answer.
#(a) The lasso, relative to least squares, is: Less fexible and hence will give improved prediction accuracy when its increase in bias is less than its decrease in Variance. ##The lasso shrinks coefficients toward zero (some exactly to zero), constraining the model relative to unconstrained least squares making it less flexible. This adds some bias but substantially reduces variance, improving prediction accuracy whenever that variance reduction outweighs the added bias.
b) 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. #Ridge regression also constrains the coefficient estimates (shrinking them toward zero via the L2 penalty), making it less flexible than least squares, which is unconstrained. Like the lasso, this introduces bias but reduces variance particularly useful when predictors are correlated or when p is large relative to n, situations where least squares estimates tend to have high variance. Ridge improves on least squares prediction accuracy precisely when this variance reduction exceeds the bias introduced
#(c) Repeat (a) for non-linear methods relative to least squares: More flexible and hence will give improved prediction accuracy when its increase in variance is less than its decrease in bias. #Non-linear methods relax the linearity assumption of least squares, allowing the fitted model to follow more complex, curved relationships in the data making them more flexible. This flexibility reduces bias (the model can better capture the true underlying relationship, especially if it’s non-linear), but at the cost of increased variance (more flexible models fit the noise in training data more closely).
Question 9 :
library(ISLR2)
data(College)
set.seed(1)
train <- sample(1:nrow(College), nrow(College) / 2)
test <- -train
College.train <- College[train, ]
College.test <- College[test, ]
dim(College.train)## [1] 388 18dim(College.test)## [1] 389 18lm.fit <- lm(Apps ~ ., data = College.train)
summary(lm.fit)##
## Call:
## lm(formula = Apps ~ ., data = College.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5741.2 -479.5 15.3 359.6 7258.0
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.902e+02 6.381e+02 -1.238 0.216410
## PrivateYes -3.070e+02 2.006e+02 -1.531 0.126736
## Accept 1.779e+00 5.420e-02 32.830 < 2e-16 ***
## Enroll -1.470e+00 3.115e-01 -4.720 3.35e-06 ***
## Top10perc 6.673e+01 8.310e+00 8.030 1.31e-14 ***
## Top25perc -2.231e+01 6.533e+00 -3.415 0.000708 ***
## F.Undergrad 9.269e-02 5.529e-02 1.676 0.094538 .
## P.Undergrad 9.397e-03 5.493e-02 0.171 0.864275
## Outstate -1.084e-01 2.700e-02 -4.014 7.22e-05 ***
## Room.Board 2.115e-01 7.224e-02 2.928 0.003622 **
## Books 2.912e-01 3.985e-01 0.731 0.465399
## Personal 6.133e-03 8.803e-02 0.070 0.944497
## PhD -1.548e+01 6.681e+00 -2.316 0.021082 *
## Terminal 6.415e+00 7.290e+00 0.880 0.379470
## S.F.Ratio 2.283e+01 2.047e+01 1.115 0.265526
## perc.alumni 1.134e+00 6.083e+00 0.186 0.852274
## Expend 4.857e-02 1.619e-02 2.999 0.002890 **
## Grad.Rate 7.490e+00 4.397e+00 1.703 0.089324 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1083 on 370 degrees of freedom
## Multiple R-squared: 0.9389, Adjusted R-squared: 0.9361
## F-statistic: 334.3 on 17 and 370 DF, p-value: < 2.2e-16lm.pred <- predict(lm.fit, College.test)
lm.err <- mean((lm.pred - College.test$Apps)^2)
lm.err## [1] 1135758library(glmnet)## Loading required package: Matrix## Loaded glmnet 4.1-10train.mat <- model.matrix(Apps ~ ., data = College.train)[, -1]
test.mat <- model.matrix(Apps ~ ., data = College.test)[, -1]
set.seed(1)
grid <- 10^seq(10, -2, length = 100)
ridge.mod <- cv.glmnet(train.mat, College.train$Apps, alpha = 0, lambda = grid)
best.lambda.ridge <- ridge.mod$lambda.min
best.lambda.ridge## [1] 0.01ridge.pred <- predict(ridge.mod, s = best.lambda.ridge, newx = test.mat)
ridge.err <- mean((ridge.pred - College.test$Apps)^2)
ridge.err## [1] 1134677##Ridge regression: chosen λ = 0.01 essentially the smallest value in our grid, meaning very little shrinkage was preferred, giving test MSE = 1,134,677 nearly identical to least squares (1,135,758), which makes sense since λ is so close to 0 that ridge is behaving almost like ordinary least squares here.
set.seed(1)
lasso.mod <- cv.glmnet(train.mat, College.train$Apps, alpha = 1, lambda = grid)
best.lambda.lasso <- lasso.mod$lambda.min
best.lambda.lasso## [1] 0.01lasso.pred <- predict(lasso.mod, s = best.lambda.lasso, newx = test.mat)
lasso.err <- mean((lasso.pred - College.test$Apps)^2)
lasso.err## [1] 1133422lasso.coef <- predict(lasso.mod, s = best.lambda.lasso, type = "coefficients")
lasso.coef## 18 x 1 sparse Matrix of class "dgCMatrix"
## s=0.01
## (Intercept) -7.931498e+02
## PrivateYes -3.078903e+02
## Accept 1.777242e+00
## Enroll -1.450532e+00
## Top10perc 6.659456e+01
## Top25perc -2.221506e+01
## F.Undergrad 8.983869e-02
## P.Undergrad 1.005260e-02
## Outstate -1.082871e-01
## Room.Board 2.118762e-01
## Books 2.922508e-01
## Personal 6.234085e-03
## PhD -1.542914e+01
## Terminal 6.364841e+00
## S.F.Ratio 2.284667e+01
## perc.alumni 1.114025e+00
## Expend 4.861825e-02
## Grad.Rate 7.466015e+00sum(lasso.coef != 0)## [1] 18library(pls)## Warning: package 'pls' was built under R version 4.4.3##
## Attaching package: 'pls'## The following object is masked from 'package:stats':
##
## loadingsset.seed(1)
pcr.fit <- pcr(Apps ~ ., data = College.train, scale = TRUE, validation = "CV")
summary(pcr.fit)## Data: X dimension: 388 17
## Y dimension: 388 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 4288 4006 2373 2372 2069 1961 1919
## adjCV 4288 4007 2368 2369 1999 1948 1911
## 7 comps 8 comps 9 comps 10 comps 11 comps 12 comps 13 comps
## CV 1919 1921 1876 1832 1832 1836 1837
## adjCV 1912 1915 1868 1821 1823 1827 1827
## 14 comps 15 comps 16 comps 17 comps
## CV 1853 1759 1341 1270
## adjCV 1850 1733 1326 1257
##
## TRAINING: % variance explained
## 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps 7 comps 8 comps
## X 32.20 57.78 65.31 70.99 76.37 81.27 84.8 87.85
## Apps 13.44 70.93 71.07 79.87 81.15 82.25 82.3 82.33
## 9 comps 10 comps 11 comps 12 comps 13 comps 14 comps 15 comps
## X 90.62 92.91 94.98 96.74 97.79 98.72 99.42
## Apps 83.38 84.76 84.80 84.84 85.11 85.14 90.55
## 16 comps 17 comps
## X 99.88 100.00
## Apps 93.42 93.89validationplot(pcr.fit, val.type = "MSEP")
pcr.pred <- predict(pcr.fit, College.test, ncomp = 17)
pcr.err <- mean((pcr.pred - College.test$Apps)^2)
pcr.err## [1] 1135758set.seed(1)
pls.fit <- plsr(Apps ~ ., data = College.train, scale = TRUE, validation = "CV")
summary(pls.fit)## Data: X dimension: 388 17
## Y dimension: 388 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 4288 2217 2019 1761 1630 1533 1347
## adjCV 4288 2211 2012 1749 1605 1510 1331
## 7 comps 8 comps 9 comps 10 comps 11 comps 12 comps 13 comps
## CV 1309 1303 1286 1283 1283 1277 1271
## adjCV 1296 1289 1273 1270 1270 1264 1258
## 14 comps 15 comps 16 comps 17 comps
## CV 1270 1270 1270 1270
## adjCV 1258 1257 1257 1257
##
## TRAINING: % variance explained
## 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps 7 comps 8 comps
## X 27.21 50.73 63.06 65.52 70.20 74.20 78.62 80.81
## Apps 75.39 81.24 86.97 91.14 92.62 93.43 93.56 93.68
## 9 comps 10 comps 11 comps 12 comps 13 comps 14 comps 15 comps
## X 83.29 87.17 89.15 91.37 92.58 94.42 96.98
## Apps 93.76 93.79 93.83 93.86 93.88 93.89 93.89
## 16 comps 17 comps
## X 98.78 100.00
## Apps 93.89 93.89#plateaus around M = 9 or 10 components (CV ≈ 1283–1286), with only marginal further improvement out to M = 17 (RMSEP = 1270). The minimum technically occurs at M = 14–17 (tied at 1270), but since the improvement beyond M ≈ 10 is negligible, a reasonable practical choice is M = 10, which achieves nearly the same performance with far fewer components (much more parsimonious than PCR’s M=17).
pls.pred <- predict(pls.fit, College.test, ncomp = 10)
pls.err <- mean((pls.pred - College.test$Apps)^2)
pls.err## [1] 1131661#PLS with M = 10 components gives test MSE = 1,131,661 — the best (lowest) test MSE among all five methods so far, and it achieves this using only 10 components instead of the full 17 predictors.
test.avg <- mean(College.test$Apps)
tss <- sum((College.test$Apps - test.avg)^2)
lm.r2 <- 1 - sum((lm.pred - College.test$Apps)^2) / tss
ridge.r2 <- 1 - sum((ridge.pred - College.test$Apps)^2) / tss
lasso.r2 <- 1 - sum((lasso.pred - College.test$Apps)^2) / tss
pcr.r2 <- 1 - sum((pcr.pred - College.test$Apps)^2) / tss
pls.r2 <- 1 - sum((pls.pred - College.test$Apps)^2) / tss
c(lm.r2, ridge.r2, lasso.r2, pcr.r2, pls.r2)## [1] 0.9015413 0.9016351 0.9017438 0.9015413 0.9018965##We can predict the number of college applications received quite accurately — all five methods explain roughly 90% of the variance in Apps on the held-out test data (R² ≈ 0.90–0.90). This is a strong result.
#the test MSEs range only from about 1,131,661 to 1,135,758, a difference of less than 0.4%, and the R² values are essentially indistinguishable (all between 0.9015 and 0.9019). Practically speaking, all five methods perform almost identically on this dataset.
#The overall takeaway is that for this particular dataset, the predictors don’t appear to be so collinear or numerous relative to n that regularization/dimension reduction provides a substantial edge over ordinary least squares College has only 17 predictors and 388 training observations, a comfortable n >> p regime where OLS is already well-behaved and doesn’t overfit much. Regularization methods tend to shine more when p is large relative to n or predictors are highly collinear, neither of which is a major issue here.
library(ISLR2)
library(leaps)
data(Boston)
set.seed(1)
dim(Boston)## [1] 506 13sum(is.na(Boston))## [1] 0predict.regsubsets <- function(object, newdata, id, formula, ...) {
mat <- model.matrix(formula, newdata)
coefi <- coef(object, id = id)
xvars <- names(coefi)
mat[, xvars] %*% coefi
}
k <- 10
n <- nrow(Boston)
set.seed(1)
folds <- sample(1:k, n, replace = TRUE)
cv.errors <- matrix(NA, k, 13, dimnames = list(NULL, paste(1:13)))
for (j in 1:k) {
best.fit <- regsubsets(crim ~ ., data = Boston[folds != j, ], nvmax = 13)
n.sizes <- nrow(summary(best.fit)$which) # how many sizes actually available
for (i in 1:min(13, n.sizes)) {
pred <- predict.regsubsets(best.fit, Boston[folds == j, ], id = i, formula = crim ~ .)
cv.errors[j, i] <- mean((Boston$crim[folds == j] - pred)^2)
}
}
mean.cv.errors <- apply(cv.errors, 2, mean, na.rm = TRUE)
mean.cv.errors## 1 2 3 4 5 6 7 8
## 45.44573 43.87260 43.95794 43.80360 43.23984 43.06580 42.52569 42.69789
## 9 10 11 12 13
## 42.56403 42.47280 42.39886 42.48245 NaNwhich.min(mean.cv.errors)## 11
## 11reg.best <- regsubsets(crim ~ ., data = Boston, nvmax = 13)
coef(reg.best, 12)## (Intercept) zn indus chas nox
## 13.7783937863 0.0457100386 -0.0583501107 -0.8253775522 -9.9575865471
## rm age dis rad tax
## 0.6289106622 -0.0008482791 -1.0122467382 0.6124653115 -0.0037756465
## ptratio lstat medv
## -0.3040727572 0.1388005968 -0.2200563590library(glmnet)
x <- model.matrix(crim ~ ., Boston)[, -1]
y <- Boston$crim
set.seed(1)
cv.ridge <- cv.glmnet(x, y, alpha = 0)
best.lambda.ridge <- cv.ridge$lambda.min
best.lambda.ridge## [1] 0.5374992min(cv.ridge$cvm)## [1] 42.64358set.seed(1)
cv.lasso <- cv.glmnet(x, y, alpha = 1)
best.lambda.lasso <- cv.lasso$lambda.min
best.lambda.lasso## [1] 0.02675146min(cv.lasso$cvm)## [1] 42.33945lasso.coef <- predict(cv.lasso, s = best.lambda.lasso, type = "coefficients")
lasso.coef## 13 x 1 sparse Matrix of class "dgCMatrix"
## s=0.02675146
## (Intercept) 11.455698159
## zn 0.041001602
## indus -0.065035444
## chas -0.734769185
## nox -8.380135083
## rm 0.516753238
## age .
## dis -0.900912511
## rad 0.568102801
## tax -0.001360719
## ptratio -0.262679725
## lstat 0.136742891
## medv -0.198157635#Lasso: chosen λ = 0.0563, CV error = 42.52 — essentially tied with best subset (42.46) and slightly better than ridge (42.71). Lasso zeroed out 2 predictors (age and tax), keeping 11 non-zero coefficients notably, it also dropped age, consistent with best subset selection
library(pls)
set.seed(1)
pcr.fit <- pcr(crim ~ ., data = Boston, scale = TRUE, validation = "CV")
summary(pcr.fit)## Data: X dimension: 506 12
## Y dimension: 506 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 8.61 7.315 7.318 6.917 6.887 6.849 6.840
## adjCV 8.61 7.310 7.312 6.908 6.882 6.842 6.833
## 7 comps 8 comps 9 comps 10 comps 11 comps 12 comps
## CV 6.705 6.733 6.715 6.702 6.684 6.607
## adjCV 6.694 6.724 6.705 6.693 6.672 6.594
##
## TRAINING: % variance explained
## 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps 7 comps 8 comps
## X 49.93 63.64 72.94 80.21 86.83 90.26 92.79 94.99
## crim 29.39 29.55 37.39 37.85 38.85 39.23 41.73 41.82
## 9 comps 10 comps 11 comps 12 comps
## X 96.78 98.33 99.48 100.00
## crim 42.12 42.43 43.58 44.93M = 4: RMSEP = 6.815 → MSE ≈ 46.44
M = 13: RMSEP = 6.624 → MSE ≈ 43.88
##for (a) : Best subset selection, ridge, and lasso all perform very similarly, with CV MSEs clustered tightly between about 42.46 and 42.71. PCR performs slightly worse even at its best (M=13, MSE≈43.88), and considerably worse if we restrict to a small number of components (M=4, MSE≈46.44) suggesting that reducing dimensionality via principal components (which are formed without reference to the response) loses some predictive signal relative to methods that select or shrink coefficients using their relationship to crim directl.
#Given these results, the lasso stands out as the most attractive choice. Its cross-validated MSE (42.52) is essentially tied with the best CV error achieved by any method — best subset selection at size 12 (42.46) — and is slightly better than ridge (42.71). Unlike ridge and PCR, however, the lasso also performs automatic variable selection, dropping age and tax and reducing the model to 11 predictors. This yields a more interpretable, parsimonious model with no real cost in predictive accuracy. I propose the lasso model with λ = 0.056 as the final choice. Its CV error nearly matches the best result of any method tested, it naturally performs variable selection to improve interpretability, and it is more computationally efficient and stable than an exhaustive best subset search, particularly as the number of predictors grows. Importantly, all comparisons were based on 10-fold cross-validation error rather than training error, ensuring a fair evaluation of out-of-sample performance. ## the model model doesnt involve all of the features The lasso model excludes age and tax, retaining only 11 of the 13 available predictors. This is consistent with best subset selection, which also excluded age from its top-performing model. This makes intuitive sense: age (proportion of owner-occupied units built before 1940) likely carries information that is largely redundant with other retained predictors, such as dis (distance to employment centers, correlated with neighborhood development) and medv (median home value) so once these are included, age adds little unique predictive value for crime rate. Similarly, tax may be collinear with rad (accessibility to radial highways), as both relate to urban infrastructure and proximity, making one largely redundant once the other is in the model.