Question 2

For parts (a) through (c), indicate which of i. through iv. is correct. Justify your answer.

  1. The lasso, relative to least squares, is:
  1. More flexible and hence will give improved prediction accuracy when its increase in bias is less than its decrease in variance.
  2. More flexible and hence will give improved prediction accuracy when its increase in variance is less than its decrease in bias.
  3. Less flexible and hence will give improved prediction accuracy when its increase in bias is less than its decrease in variance.
  4. Less flexible and hence will give improved prediction accuracy when its increase in variance is less than its decrease in bias.

–Answer: iii is the correct answer: Lasso regression is less flexible than least squares. It improves prediction accuracy because the reduction in variance is greater than the increase in bias.

  1. Repeat (a) for ridge regression relative to least squares.

–Answer: iii is the correct statement: being less flexible than least squares, Ridge regression improves prediction accuracy when variance decreases more than bias increases.

  1. Repeat (a) for non-linear methods relative to least squares.

–Answer: ii is the correct stmnt: being more flexible than lease squares, non-linear method improves prediction accuracy when variance increases less than bias decreases.

Question 9

In this exercise, we will predict the number of applications received using the other variables in the College data set.

  1. Split the data set into a training set and a test set.
library(ISLR2)
library(glmnet)
## Loading required package: Matrix
## Loaded glmnet 5.0
library(pls)
## 
## Attaching package: 'pls'
## The following object is masked from 'package:stats':
## 
##     loadings
set.seed(1)

# split
train_idx <-sample(nrow(College), nrow(College) /2)
train_set <-College[train_idx, ]
test_set <-College[-train_idx, ]

# model for glmnet
x_train <- model.matrix(Apps ~ ., data = train_set) [, -1]
y_train <- train_set$Apps
x_test <- model.matrix(Apps ~ ., data = test_set)[, -1]
y_test <- test_set$Apps
  1. Fit a linear model using least squares on the training set, and report the test error obtained.
lm.fit <- lm(Apps ~ ., data = train_set)
lm.pred <- predict(lm.fit, test_set)
lm.err <- mean((test_set$Apps - lm.pred)^2)
cat("OLS Linear Model Test MSE:", lm.err, "\n")
## OLS Linear Model Test MSE: 1135758
  1. Fit a ridge regression model on the training set, with λ chosen by cross-validation. Report the test error obtained.
set.seed(1)
cv.ridge <- cv.glmnet(x_train, y_train, alpha = 0)
bestlambda.ridge <- cv.ridge$lambda.min
cat("Ridge Regression Best Lambda:", bestlambda.ridge, "\n")
## Ridge Regression Best Lambda: 405.8404
ridge.pred <- predict(cv.ridge, s = bestlambda.ridge, newx = x_test)
ridge.err <- mean((y_test - ridge.pred)^2)
cat("Ridge Regression Test MSE:", ridge.err, "\n")
## Ridge Regression Test MSE: 976261.5
  1. Fit a lasso model on the training set, with λ chosen by cross- validation. Report the test error obtained, along with the num- ber of non-zero coefficient estimates.
  set.seed(1)
cv.lasso <- cv.glmnet(x_train, y_train, alpha = 1)
bestlambda.lasso <- cv.lasso$lambda.min
cat("Lasso Best Lambda:", bestlambda.lasso, "\n")
## Lasso Best Lambda: 1.97344
lasso.pred <- predict(cv.lasso, s = bestlambda.lasso, newx = x_test)
lasso.err <- mean((y_test - lasso.pred)^2)
cat("Lasso Test MSE:", lasso.err, "\n")
## Lasso Test MSE: 1115901
# Ext coefficients 
lasso.coef <- predict(cv.lasso, type = "coefficients", s = bestlambda.lasso)

# Count non-zero coefficients excluding intercept
non_zero_coefs <- sum(lasso.coef[-1, 1] != 0)
cat("Number of Non-zero Coefficients (excluding intercept):", non_zero_coefs, "\n")
## Number of Non-zero Coefficients (excluding intercept): 17
print(lasso.coef)
## 18 x 1 sparse Matrix of class "dgCMatrix"
##                 s=1.97344
## (Intercept) -7.688896e+02
## PrivateYes  -3.127034e+02
## Accept       1.762718e+00
## Enroll      -1.318195e+00
## Top10perc    6.482356e+01
## Top25perc   -2.081406e+01
## F.Undergrad  7.119149e-02
## P.Undergrad  1.246161e-02
## Outstate    -1.049091e-01
## Room.Board   2.088305e-01
## Books        2.926466e-01
## Personal     3.955068e-03
## PhD         -1.455463e+01
## Terminal     5.395858e+00
## S.F.Ratio    2.171398e+01
## perc.alumni  5.088260e-01
## Expend       4.824455e-02
## Grad.Rate    7.036148e+00
  1. 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.
set.seed(1)
pcr.fit <- pcr(Apps ~ ., data = train_set, scale = TRUE, validation = "CV")
validationplot(pcr.fit, val.type = "MSEP")

# Min CV error
pcr.cv.err <- RMSEP(pcr.fit)$val[1,, -1] 
best_M_pcr <- which.min(pcr.cv.err)
cat("PCR Best Number of Components (M):", best_M_pcr, "\n")
## PCR Best Number of Components (M): 17
pcr.pred <- predict(pcr.fit, test_set, ncomp = best_M_pcr)
pcr.err <- mean((test_set$Apps - pcr.pred)^2)
cat("PCR Test MSE:", pcr.err, "\n")
## PCR Test MSE: 1135758
  1. 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.
set.seed(1)
pls.fit <- plsr(Apps ~ ., data = train_set, scale = TRUE, validation = "CV")
validationplot(pls.fit, val.type = "MSEP")

# Find component size M that minimizes CV error
pls.cv.err <- RMSEP(pls.fit)$val[1,, -1] # Exclude 0 components
best_M_pls <- which.min(pls.cv.err)
cat("PLS Best Number of Components (M):", best_M_pls, "\n")
## PLS Best Number of Components (M): 17
pls.pred <- predict(pls.fit, test_set, ncomp = best_M_pls)
pls.err <- mean((test_set$Apps - pls.pred)^2)
cat("PLS Test MSE:", pls.err, "\n")
## PLS Test MSE: 1135758
  1. Comment on the results obtained. How accurately can we pre- dict the number of college applications received? Is there much difference among the test errors resulting from these five ap- proaches?
test_avg <- mean(test_set$Apps)
test_var <- sum((test_set$Apps - test_avg)^2)

# Compute R-squared values
lm.r2 <- 1 - sum((test_set$Apps - lm.pred)^2) / test_var
ridge.r2 <- 1 - sum((y_test - ridge.pred)^2) / test_var
lasso.r2 <- 1 - sum((y_test - lasso.pred)^2) / test_var
pcr.r2 <- 1 - sum((test_set$Apps - pcr.pred)^2) / test_var
pls.r2 <- 1 - sum((test_set$Apps - pls.pred)^2) / test_var

# Construct comparison table
results <- data.frame(
  Method = c("OLS Linear Model", "Ridge Regression", "Lasso Regression", "PCR", "PLS"),
  Test_MSE = c(lm.err, ridge.err, lasso.err, pcr.err, pls.err),
  Test_R2 = c(lm.r2, ridge.r2, lasso.r2, pcr.r2, pls.r2)
)
print(results)
##             Method  Test_MSE   Test_R2
## 1 OLS Linear Model 1135758.3 0.9015413
## 2 Ridge Regression  976261.5 0.9153681
## 3 Lasso Regression 1115900.6 0.9032628
## 4              PCR 1135758.3 0.9015413
## 5              PLS 1135758.3 0.9015413

Question 11

We will now try to predict per capita crime rate in the Boston data set.

  1. 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.
library(leaps)

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}
set.seed(1)
k <- 10

# Boston Partition into 10 folds
folds <- sample(1:k, nrow(Boston), replace = TRUE)

err.best <- numeric(k)
err.ridge <- numeric(k)
err.lasso <- numeric(k)
err.pcr <- numeric(k)

# matrices for glmnet
x_boston <- model.matrix(crim ~ ., data = Boston)[, -1]
y_boston <- Boston$crim

for (j in 1:k) {
  train_fold <- (folds != j)
  test_fold <- (folds == j)
  
  # 1. Best Selection
  fit.best <- regsubsets(crim ~ ., data = Boston[train_fold, ], nvmax = 13)
  best_size_bic <- which.min(summary(fit.best)$bic)
  pred.best <- predict(fit.best, Boston[test_fold, ], id = best_size_bic)
  err.best[j] <- mean((y_boston[test_fold] - pred.best)^2)
  
  # 2. Ridge 
  cv.r <- cv.glmnet(x_boston[train_fold, ], y_boston[train_fold], alpha = 0)
  pred.r <- predict(cv.r, s = cv.r$lambda.min, newx = x_boston[test_fold, ])
  err.ridge[j] <- mean((y_boston[test_fold] - pred.r)^2)
  
  # 3. Lasso 
  cv.l <- cv.glmnet(x_boston[train_fold, ], y_boston[train_fold], alpha = 1)
  pred.l <- predict(cv.l, s = cv.l$lambda.min, newx = x_boston[test_fold, ])
  err.lasso[j] <- mean((y_boston[test_fold] - pred.l)^2)
  
  # 4. PCR
  fit.pcr <- pcr(crim ~ ., data = Boston[train_fold, ], scale = TRUE, validation = "CV")
  pcr.cv.err <- RMSEP(fit.pcr)$val[1,, -1]
  best_M <- which.min(pcr.cv.err)
  pred.pcr <- predict(fit.pcr, Boston[test_fold, ], ncomp = best_M)
  err.pcr[j] <- mean((y_boston[test_fold] - pred.pcr)^2)
}

# avg CV error rate 
cat("Best Subset CV MSE:", mean(err.best), "\n")
## Best Subset CV MSE: 43.8726
cat("Ridge CV MSE:", mean(err.ridge), "\n")
## Ridge CV MSE: 42.68628
cat("Lasso CV MSE:", mean(err.lasso), "\n")
## Lasso CV MSE: 42.23101
cat("PCR CV MSE:", mean(err.pcr), "\n")
## PCR CV MSE: 42.48245
  1. 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.

Answer: The Lasso Regression model is best for this dataset, achieving the lowest cross-validation test MSE and best predictive accuracy

  1. Does your chosen model involve all of the features in the data set? Why or why not?
lasso.full <- cv.glmnet(x_boston, y_boston, alpha = 1)
best_lambda <- lasso.full$lambda.min
coef(lasso.full, s = best_lambda)
## 13 x 1 sparse Matrix of class "dgCMatrix"
##              s=0.03881179
## (Intercept) 10.4134569616
## zn           0.0389011405
## indus       -0.0676069618
## chas        -0.6935740430
## nox         -7.6342999354
## rm           0.4679060817
## age          .           
## dis         -0.8526255420
## rad          0.5485729634
## tax         -0.0003100068
## ptratio     -0.2437901821
## lstat        0.1360221554
## medv        -0.1883947640

Answer: No, the Lasso excludes some features. By performing automatic variable selection, it drops redundant or non-informative variables to simplify interpretation and lower prediction variance without adding significant bias.