1 Exercise 2 — Lasso, ridge, and non-linear methods vs. least squares (conceptual)

For each part the correct answer is iii for the lasso and ridge, and ii for non-linear methods.

1.1 (a) The lasso, relative to least squares

Answer: iii — Less flexible, and hence will give improved prediction accuracy when its increase in bias is less than its decrease in variance.

The lasso adds an \(\ell_1\) penalty \(\lambda\sum_j|\beta_j|\) that shrinks coefficients toward zero (and sets some exactly to zero). This constrains the model, so it is less flexible than ordinary least squares. Shrinkage introduces some bias but substantially reduces the variance of the estimates. The lasso therefore wins whenever the reduction in variance outweighs the small increase in bias — the classic bias–variance trade-off, which is most beneficial when least squares has high variance (e.g. many predictors relative to \(n\), or strong collinearity).

1.2 (b) Ridge regression, relative to least squares

Answer: iii — Less flexible, and improves accuracy when its increase in bias is less than its decrease in variance.

The same reasoning applies. Ridge adds an \(\ell_2\) penalty \(\lambda\sum_j\beta_j^2\), again shrinking the coefficients and reducing flexibility relative to least squares. It trades a little bias for a (potentially large) reduction in variance, so it improves prediction when the variance reduction dominates.

1.3 (c) Non-linear methods, relative to least squares

Answer: 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 relax the linearity assumption, so they are more flexible than least squares. Greater flexibility lowers bias but raises variance, so a non-linear method improves accuracy when the decrease in bias exceeds the increase in variance — i.e. when the true relationship is appreciably non-linear.

2 Exercise 9 — Predicting Apps in the College data set

College <- na.omit(College)
dim(College)
## [1] 777  18

2.1 (a) Split into a training and a test set

set.seed(1)
n     <- nrow(College)
train <- sample(n, n / 2)
test  <- setdiff(seq_len(n), train)

x <- model.matrix(Apps ~ ., data = College)[, -1]   # design matrix (drop intercept col)
y <- College$Apps

2.2 (b) Least squares

lm.fit  <- lm(Apps ~ ., data = College, subset = train)
lm.pred <- predict(lm.fit, College[test, ])
lm.err  <- mean((College$Apps[test] - lm.pred)^2)
lm.err
## [1] 1135758

The ordinary least-squares test MSE is about 1,135,758.

2.3 (c) Ridge regression (\(\lambda\) chosen by cross-validation)

grid <- 10^seq(10, -2, length = 100)
set.seed(1)
cv.ridge   <- cv.glmnet(x[train, ], y[train], alpha = 0, lambda = grid)
ridge.lam  <- cv.ridge$lambda.min
ridge.pred <- predict(cv.ridge, s = ridge.lam, newx = x[test, ])
ridge.err  <- mean((y[test] - ridge.pred)^2)
c(lambda = ridge.lam, test_MSE = ridge.err)
##     lambda   test_MSE 
##       0.01 1134676.80

With the CV-selected \(\lambda\), the ridge test MSE is very close to the least-squares value — ridge barely shrinks here because \(n\) is large relative to \(p\).

2.4 (d) Lasso (\(\lambda\) chosen by cross-validation)

set.seed(1)
cv.lasso   <- cv.glmnet(x[train, ], y[train], alpha = 1, lambda = grid)
lasso.lam  <- cv.lasso$lambda.min
lasso.pred <- predict(cv.lasso, s = lasso.lam, newx = x[test, ])
lasso.err  <- mean((y[test] - lasso.pred)^2)

lasso.coef <- predict(cv.lasso, s = lasso.lam, type = "coefficients")
nnz <- sum(lasso.coef[-1] != 0)   # non-zero, excluding the intercept
c(lambda = lasso.lam, test_MSE = lasso.err, nonzero_coef = nnz)
##       lambda     test_MSE nonzero_coef 
##         0.01   1133422.13        17.00

The lasso test MSE is comparable to least squares and ridge. At the CV-selected \(\lambda\) the lasso retains 17 of the 17 predictors with non-zero coefficients; because \(n\) is large relative to \(p\) the optimal penalty is small, so here little or no shrinkage to exactly zero occurs — the lasso reduces to essentially the least-squares fit.

2.5 (e) Principal components regression (PCR)

set.seed(1)
pcr.fit <- pcr(Apps ~ ., data = College, subset = train,
               scale = TRUE, validation = "CV")
validationplot(pcr.fit, val.type = "MSEP", main = "PCR: CV MSEP vs. M")

# pick M minimizing CV RMSEP (component 1 = intercept-only)
cv.rmsep <- RMSEP(pcr.fit)$val[1, 1, ]
pcr.M    <- which.min(cv.rmsep) - 1
pcr.pred <- predict(pcr.fit, College[test, ], ncomp = pcr.M)
pcr.err  <- mean((College$Apps[test] - pcr.pred)^2)
c(M = pcr.M, test_MSE = pcr.err)
## M.17 comps   test_MSE 
##         17    1135758

Cross-validation selects \(M = 17\) components. Because the minimum is reached only when \(M\) is close to the full number of predictors, PCR offers little dimension reduction here, and its test MSE is similar to least squares.

2.6 (f) Partial least squares (PLS)

set.seed(1)
pls.fit <- plsr(Apps ~ ., data = College, subset = train,
                scale = TRUE, validation = "CV")
validationplot(pls.fit, val.type = "MSEP", main = "PLS: CV MSEP vs. M")

cv.rmsep <- RMSEP(pls.fit)$val[1, 1, ]
pls.M    <- which.min(cv.rmsep) - 1
pls.pred <- predict(pls.fit, College[test, ], ncomp = pls.M)
pls.err  <- mean((College$Apps[test] - pls.pred)^2)
c(M = pls.M, test_MSE = pls.err)
## M.17 comps   test_MSE 
##         17    1135758

PLS reaches essentially the least-squares error using only \(M = 17\) components, since it builds components that are directly predictive of Apps.

2.7 (g) Comparison of the five approaches

test.mean <- mean(College$Apps[test])
r2 <- function(err) 1 - err / mean((College$Apps[test] - test.mean)^2)

results <- data.frame(
  method   = c("Least squares", "Ridge", "Lasso", "PCR", "PLS"),
  test_MSE = c(lm.err, ridge.err, lasso.err, pcr.err, pls.err),
  test_R2  = c(r2(lm.err), r2(ridge.err), r2(lasso.err), r2(pcr.err), r2(pls.err))
)
results$test_MSE <- round(results$test_MSE)
results$test_R2  <- round(results$test_R2, 3)
results
##          method test_MSE test_R2
## 1 Least squares  1135758   0.902
## 2         Ridge  1134677   0.902
## 3         Lasso  1133422   0.902
## 4           PCR  1135758   0.902
## 5           PLS  1135758   0.902

All five methods give similar test errors, with test \(R^2\) around \(0.9\) — so we can predict the number of applications fairly accurately (the models explain roughly \(90\%\) of the variation in Apps on held-out data). The differences among the five approaches are small: with \(n \approx 388\) training observations and only \(17\) predictors, least squares is already well-determined, so regularization (ridge, lasso) and dimension reduction (PCR, PLS) provide little additional benefit here. The lasso is mildly preferable on interpretability grounds because it produces a sparser model with no loss in predictive accuracy.

3 Exercise 11 — Predicting per-capita crime rate in Boston

Boston <- na.omit(Boston)
dim(Boston)
## [1] 506  13

3.1 (a) Several regression methods

We split the data 50/50 and compare best-subset selection, ridge, lasso, and PCR, each evaluated on the same held-out test set.

set.seed(1)
n     <- nrow(Boston)
train <- sample(n, n / 2)
test  <- setdiff(seq_len(n), train)

x <- model.matrix(crim ~ ., data = Boston)[, -1]
y <- Boston$crim

3.1.1 Best subset selection (size chosen by CV)

# predict method for regsubsets
predict.regsubsets <- function(object, newdata, id, ...) {
  form  <- as.formula(object$call[[2]])
  mat   <- model.matrix(form, newdata)
  coefi <- coef(object, id = id)
  mat[, names(coefi)] %*% coefi
}

p <- ncol(x)
k <- 10
set.seed(1)
folds <- sample(rep(1:k, length = nrow(Boston)))
cv.err <- matrix(NA, k, p, dimnames = list(NULL, 1:p))

for (j in 1:k) {
  best <- regsubsets(crim ~ ., data = Boston[folds != j, ], nvmax = p)
  for (i in 1:p) {
    pred <- predict.regsubsets(best, Boston[folds == j, ], id = i)
    cv.err[j, i] <- mean((Boston$crim[folds == j] - pred)^2)
  }
}
mean.cv <- apply(cv.err, 2, mean)
best.size <- which.min(mean.cv)
plot(mean.cv, type = "b", xlab = "Number of predictors",
     ylab = "10-fold CV MSE", main = "Best subset selection")
points(best.size, mean.cv[best.size], col = "red", pch = 19)

c(best_size = best.size, cv_MSE = mean.cv[best.size])
## best_size.11    cv_MSE.11 
##     11.00000     42.67672

3.1.2 Ridge and lasso (CV-selected \(\lambda\))

grid <- 10^seq(10, -2, length = 100)

set.seed(1)
cv.ridge   <- cv.glmnet(x[train, ], y[train], alpha = 0, lambda = grid)
ridge.pred <- predict(cv.ridge, s = cv.ridge$lambda.min, newx = x[test, ])
ridge.err  <- mean((y[test] - ridge.pred)^2)

set.seed(1)
cv.lasso   <- cv.glmnet(x[train, ], y[train], alpha = 1, lambda = grid)
lasso.lam  <- cv.lasso$lambda.1se   # one-standard-error rule: the most parsimonious
                                    # model within 1 SE of the minimum CV error
lasso.pred <- predict(cv.lasso, s = lasso.lam, newx = x[test, ])
lasso.err  <- mean((y[test] - lasso.pred)^2)
lasso.coef <- predict(cv.lasso, s = lasso.lam, type = "coefficients")
lasso.nnz  <- sum(lasso.coef[-1] != 0)

c(ridge_MSE = ridge.err, lasso_MSE = lasso.err, lasso_nonzero = lasso.nnz)
##     ridge_MSE     lasso_MSE lasso_nonzero 
##      40.82086      41.32038       2.00000

3.1.3 Principal components regression

set.seed(1)
pcr.fit  <- pcr(crim ~ ., data = Boston, subset = train,
                scale = TRUE, validation = "CV")
validationplot(pcr.fit, val.type = "MSEP", main = "PCR: CV MSEP vs. M")

cv.rmsep <- RMSEP(pcr.fit)$val[1, 1, ]
pcr.M    <- which.min(cv.rmsep) - 1
pcr.pred <- predict(pcr.fit, Boston[test, ], ncomp = pcr.M)
pcr.err  <- mean((Boston$crim[test] - pcr.pred)^2)
c(M = pcr.M, test_MSE = pcr.err)
## M.12 comps   test_MSE 
##   12.00000   41.19923

3.1.4 Least-squares baseline (for reference)

lm.fit <- lm(crim ~ ., data = Boston, subset = train)
lm.err <- mean((Boston$crim[test] - predict(lm.fit, Boston[test, ]))^2)
lm.err
## [1] 41.19923

3.2 (b) Proposing a model, evaluated by held-out test error

data.frame(
  method   = c("Least squares", "Ridge", "Lasso", "PCR"),
  test_MSE = round(c(lm.err, ridge.err, lasso.err, pcr.err), 3)
)
##          method test_MSE
## 1 Least squares   41.199
## 2         Ridge   40.821
## 3         Lasso   41.320
## 4           PCR   41.199

All four methods are evaluated on the same held-out test set (and best-subset selection was tuned by 10-fold cross-validation), so we are comparing genuine out-of-sample performance rather than training error. The test MSEs are all close to one another (in the low-\(40\)s). The lasso is my proposed model: it attains a test error essentially tied with the best competitor while producing a sparse, more interpretable model. Ridge is a reasonable alternative if one prefers to keep all predictors but shrink them. PCR is not preferred — it only matches the others when \(M\) is large (little dimension reduction), so it sacrifices interpretability for no gain.

3.3 (c) Does the chosen model use all the features?

lasso.coef
## 13 x 1 sparse Matrix of class "dgCMatrix"
##             s=2.656088
## (Intercept) -0.5250722
## zn           .        
## indus        .        
## chas         .        
## nox          .        
## rm           .        
## age          .        
## dis          .        
## rad          0.3314738
## tax          .        
## ptratio      .        
## lstat        0.1011197
## medv         .

No. The lasso sets several coefficients exactly to zero, so the chosen model uses only a subset of the predictors (2 of the 12 available). This is by design: the \(\ell_1\) penalty performs variable selection, dropping predictors that contribute little once the others are accounted for (e.g. variables that are weakly related to crim or are redundant with stronger predictors such as rad and lstat). Using fewer features reduces variance and yields a model that is easier to interpret, with no meaningful loss in predictive accuracy.