Assignment 5

Question 2

2(a) Lasso compared with least squares

Correct answer: iii.

The lasso is less flexible than least squares because it shrinks some coefficients, and some can become zero. It can improve prediction accuracy when the decrease in variance is greater than the increase in bias.

2(b) Ridge regression compared with least squares

Correct answer: iii.

Ridge regression is less flexible than least squares because it shrinks the coefficients toward zero. It can improve prediction accuracy when the decrease in variance is greater than the increase in bias.

2(c) Non-linear methods compared with least squares

Correct answer: ii.

Non-linear methods are more flexible than least squares because they can fit more complex patterns. They can improve prediction accuracy if the increase in variance is smaller than the decrease in bias.Question 9: College Applications

The goal is to predict Apps, the number of applications received by each college.

data(College)
college <- College

# Look at the data.
dim(college)

[1] 777  18

head(college)

                             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

9(a) Split the data into training and test sets

We will randomly place approximately 70% of the observations in the training set and 30% in the test set.

set.seed(123)

train_id <- sample(
  1:nrow(college),
  size = round(0.70 * nrow(college))
)

college_train <- college[train_id, ]
college_test  <- college[-train_id, ]

nrow(college_train)

[1] 544

nrow(college_test)

[1] 233

9(b) Least squares regression

college_lm <- lm(Apps ~ ., data = college_train)

lm_predictions <- predict(
  college_lm,
  newdata = college_test
)

lm_test_mse <- mean(
  (college_test$Apps - lm_predictions)^2
)

lm_test_rmse <- sqrt(lm_test_mse)

lm_test_mse

[1] 1733856

lm_test_rmse

[1] 1316.76

The least-squares test mean squared error is 1.733856^{6}.
The test root mean squared error is 1316.76 applications.

Prepare matrices for ridge regression and lasso

glmnet() requires a numeric predictor matrix. The model.matrix() function automatically converts the Private variable into a dummy variable.

x_train <- model.matrix(Apps ~ ., data = college_train)[, -1]
y_train <- college_train$Apps

x_test <- model.matrix(Apps ~ ., data = college_test)[, -1]
y_test <- college_test$Apps

9(c) Ridge regression

For ridge regression, alpha = 0. Cross-validation is used to choose lambda.

set.seed(123)

ridge_cv <- cv.glmnet(
  x_train,
  y_train,
  alpha = 0
)

best_ridge_lambda <- ridge_cv$lambda.min

ridge_predictions <- predict(
  ridge_cv,
  s = "lambda.min",
  newx = x_test
)

ridge_test_mse <- mean(
  (y_test - ridge_predictions)^2
)

ridge_test_rmse <- sqrt(ridge_test_mse)

best_ridge_lambda

[1] 314.0251

ridge_test_mse

[1] 2981466

ridge_test_rmse

[1] 1726.692

plot(ridge_cv)

The cross-validated ridge lambda is 314.0251.
The ridge test MSE is 2.9814658^{6}.
The ridge test RMSE is 1726.69 applications.

9(d) Lasso regression

For the lasso, alpha = 1.

set.seed(123)

lasso_cv <- cv.glmnet(
  x_train,
  y_train,
  alpha = 1
)

best_lasso_lambda <- lasso_cv$lambda.min

lasso_predictions <- predict(
  lasso_cv,
  s = "lambda.min",
  newx = x_test
)

lasso_test_mse <- mean(
  (y_test - lasso_predictions)^2
)

lasso_test_rmse <- sqrt(lasso_test_mse)

lasso_coefficients <- coef(
  lasso_cv,
  s = "lambda.min"
)

# Subtract 1 so the intercept is not counted.
lasso_nonzero <- sum(lasso_coefficients != 0) - 1

best_lasso_lambda

[1] 8.149028

lasso_test_mse

[1] 1730230

lasso_test_rmse

[1] 1315.382

lasso_nonzero

[1] 16

plot(lasso_cv)

The cross-validated lasso lambda is 8.149.
The lasso test MSE is 1.73023^{6}.
The lasso test RMSE is 1315.38 applications.
The lasso selected 16 non-zero predictor coefficients.

9(e) Principal components regression

PCR selects the number of components using cross-validation.

set.seed(123)

college_pcr <- pcr(
  Apps ~ .,
  data = college_train,
  scale = TRUE,
  validation = "CV"
)

validationplot(
  college_pcr,
  val.type = "MSEP"
)

# RMSEP() includes a 0-component model first.
pcr_cv_rmse <- RMSEP(college_pcr, estimate = "CV")$val[1, 1, ]

best_pcr_m <- which.min(pcr_cv_rmse[-1])
best_pcr_m

16 comps 
      16 

pcr_predictions <- predict(
  college_pcr,
  newdata = college_test,
  ncomp = best_pcr_m
)

pcr_test_mse <- mean(
  (college_test$Apps - pcr_predictions)^2
)

pcr_test_rmse <- sqrt(pcr_test_mse)

pcr_test_mse

[1] 1855093

pcr_test_rmse

[1] 1362.018

Cross-validation selected 16 principal components.
The PCR test MSE is 1.8550934^{6}.
The PCR test RMSE is 1362.02 applications.

9(f) Partial least squares

PLS also selects the number of components using cross-validation.

set.seed(123)

college_pls <- plsr(
  Apps ~ .,
  data = college_train,
  scale = TRUE,
  validation = "CV"
)

validationplot(
  college_pls,
  val.type = "MSEP"
)

# RMSEP() includes a 0-component model first.
pls_cv_rmse <- RMSEP(college_pls, estimate = "CV")$val[1, 1, ]

best_pls_m <- which.min(pls_cv_rmse[-1])
best_pls_m

9 comps 
      9 

pls_predictions <- predict(
  college_pls,
  newdata = college_test,
  ncomp = best_pls_m
)

pls_test_mse <- mean(
  (college_test$Apps - pls_predictions)^2
)

pls_test_rmse <- sqrt(pls_test_mse)

pls_test_mse

[1] 1762688

pls_test_rmse

[1] 1327.663

Cross-validation selected 9 PLS components.
The PLS test MSE is 1.7626882^{6}.
The PLS test RMSE is 1327.66 applications.

9(g) Compare the five methods

college_results <- data.frame(
  Method = c(
    "Least Squares",
    "Ridge",
    "Lasso",
    "PCR",
    "PLS"
  ),
  Test_MSE = c(
    lm_test_mse,
    ridge_test_mse,
    lasso_test_mse,
    pcr_test_mse,
    pls_test_mse
  ),
  Test_RMSE = c(
    lm_test_rmse,
    ridge_test_rmse,
    lasso_test_rmse,
    pcr_test_rmse,
    pls_test_rmse
  )
)

college_results <- college_results[
  order(college_results$Test_MSE),
]

college_results

         Method Test_MSE Test_RMSE
3         Lasso  1730230  1315.382
1 Least Squares  1733856  1316.760
5           PLS  1762688  1327.663
4           PCR  1855093  1362.018
2         Ridge  2981466  1726.692

The model with the lowest test MSE is Lasso. Its test RMSE is about 1315.38 applications.

I compared the models using their test errors because they show how well each model predicts new data. If the test errors are very similar, I would choose the simpler model.

Question 11: Boston Crime Rate

The response variable is crim, the per capita crime rate by town.

data(Boston)
boston <- Boston

dim(boston)

[1] 506  13

head(boston)

     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

Split the Boston data

We will use the same type of 70% training and 30% test split.

set.seed(123)

boston_train_id <- sample(
  1:nrow(boston),
  size = round(0.70 * nrow(boston))
)

boston_train <- boston[boston_train_id, ]
boston_test  <- boston[-boston_train_id, ]

11(a) Fit and compare several regression methods

Least squares

boston_lm <- lm(crim ~ ., data = boston_train)

boston_lm_predictions <- predict(
  boston_lm,
  newdata = boston_test
)

boston_lm_mse <- mean(
  (boston_test$crim - boston_lm_predictions)^2
)

boston_lm_rmse <- sqrt(boston_lm_mse)

boston_lm_mse

[1] 18.67968

boston_lm_rmse

[1] 4.322

Prepare matrices for ridge and lasso

bx_train <- model.matrix(crim ~ ., data = boston_train)[, -1]
by_train <- boston_train$crim

bx_test <- model.matrix(crim ~ ., data = boston_test)[, -1]
by_test <- boston_test$crim

Ridge regression

set.seed(123)

boston_ridge_cv <- cv.glmnet(
  bx_train,
  by_train,
  alpha = 0
)

boston_ridge_predictions <- predict(
  boston_ridge_cv,
  s = "lambda.min",
  newx = bx_test
)

boston_ridge_mse <- mean(
  (by_test - boston_ridge_predictions)^2
)

boston_ridge_rmse <- sqrt(boston_ridge_mse)

boston_ridge_cv$lambda.min

[1] 0.5863068

boston_ridge_mse

[1] 17.74969

boston_ridge_rmse

[1] 4.213038

Lasso regression

set.seed(123)

boston_lasso_cv <- cv.glmnet(
  bx_train,
  by_train,
  alpha = 1
)

boston_lasso_predictions <- predict(
  boston_lasso_cv,
  s = "lambda.min",
  newx = bx_test
)

boston_lasso_mse <- mean(
  (by_test - boston_lasso_predictions)^2
)

boston_lasso_rmse <- sqrt(boston_lasso_mse)

boston_lasso_coef <- coef(
  boston_lasso_cv,
  s = "lambda.min"
)

boston_lasso_nonzero <- sum(
  boston_lasso_coef != 0
) - 1

boston_lasso_cv$lambda.min

[1] 0.06741104

boston_lasso_mse

[1] 18.31403

boston_lasso_rmse

[1] 4.27949

boston_lasso_nonzero

[1] 10

PCR

set.seed(123)

boston_pcr <- pcr(
  crim ~ .,
  data = boston_train,
  scale = TRUE,
  validation = "CV"
)

boston_pcr_cv_rmse <- RMSEP(
  boston_pcr,
  estimate = "CV"
)$val[1, 1, ]

boston_best_pcr_m <- which.min(
  boston_pcr_cv_rmse[-1]
)

boston_pcr_predictions <- predict(
  boston_pcr,
  newdata = boston_test,
  ncomp = boston_best_pcr_m
)

boston_pcr_mse <- mean(
  (boston_test$crim - boston_pcr_predictions)^2
)

boston_pcr_rmse <- sqrt(boston_pcr_mse)

boston_best_pcr_m

12 comps 
      12 

boston_pcr_mse

[1] 18.67968

boston_pcr_rmse

[1] 4.322

Best subset selection

First, I fit the best subset selection model using the training data. Then I used 10-fold cross-validation to choose the best model size. After that, I tested the final model using the test data.

predict_regsubsets <- function(object, newdata, id) {
  formula_used <- as.formula(object$call[[2]])
  new_matrix <- model.matrix(formula_used, newdata)
  coefficients <- coef(object, id = id)
  variables <- names(coefficients)

  new_matrix[, variables, drop = FALSE] %*% coefficients
}

set.seed(123)

number_of_folds <- 10

fold_id <- sample(
  rep(1:number_of_folds, length.out = nrow(boston_train))
)

max_predictors <- ncol(boston_train) - 1

cv_errors <- matrix(
  NA,
  nrow = number_of_folds,
  ncol = max_predictors
)

for (fold in 1:number_of_folds) {

  fold_training <- boston_train[fold_id != fold, ]
  fold_validation <- boston_train[fold_id == fold, ]

  subset_fit <- regsubsets(
    crim ~ .,
    data = fold_training,
    nvmax = max_predictors
  )

  for (model_size in 1:max_predictors) {

    subset_predictions <- predict_regsubsets(
      subset_fit,
      newdata = fold_validation,
      id = model_size
    )

    cv_errors[fold, model_size] <- mean(
      (fold_validation$crim - subset_predictions)^2
    )
  }
}

mean_cv_errors <- colMeans(cv_errors)

best_subset_size <- which.min(mean_cv_errors)

best_subset_size

[1] 9

mean_cv_errors

 [1] 57.46905 56.68955 56.30207 55.36947 54.84814 54.95260 55.02720 54.69178
 [9] 54.38770 54.50736 54.58297 54.61265

Fit the selected subset size using the complete training set, then evaluate it on the test set.

final_subset_fit <- regsubsets(
  crim ~ .,
  data = boston_train,
  nvmax = max_predictors
)

subset_test_predictions <- predict_regsubsets(
  final_subset_fit,
  newdata = boston_test,
  id = best_subset_size
)

boston_subset_mse <- mean(
  (boston_test$crim - subset_test_predictions)^2
)

boston_subset_rmse <- sqrt(boston_subset_mse)

selected_subset_coefficients <- coef(
  final_subset_fit,
  id = best_subset_size
)

boston_subset_mse

[1] 18.83233

boston_subset_rmse

[1] 4.339624

selected_subset_coefficients

  (Intercept)            zn           nox            rm           dis 
 15.351793874   0.053958377 -11.656963767   0.781683691  -0.988097250 
          rad           tax       ptratio         lstat          medv 
  0.699489646  -0.006010871  -0.362964361   0.128422085  -0.266541031 

11(b) Select a model using test performance

boston_results <- data.frame(
  Method = c(
    "Least Squares",
    "Ridge",
    "Lasso",
    "PCR",
    "Best Subset"
  ),
  Test_MSE = c(
    boston_lm_mse,
    boston_ridge_mse,
    boston_lasso_mse,
    boston_pcr_mse,
    boston_subset_mse
  ),
  Test_RMSE = c(
    boston_lm_rmse,
    boston_ridge_rmse,
    boston_lasso_rmse,
    boston_pcr_rmse,
    boston_subset_rmse
  )
)

boston_results <- boston_results[
  order(boston_results$Test_MSE),
]

boston_results

         Method Test_MSE Test_RMSE
2         Ridge 17.74969  4.213038
3         Lasso 18.31403  4.279490
1 Least Squares 18.67968  4.322000
4           PCR 18.67968  4.322000
5   Best Subset 18.83233  4.339624

Based on this training and test split, Ridge had the lowest test MSE, so I would choose this model.

I used the test error to compare the models because it shows how well they predict new data. If two models have very similar test errors, I would choose the simpler one because it is easier to understand and gives similar results.

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

# Variables selected by the lasso
lasso_selected_names <- rownames(boston_lasso_coef)[
  as.vector(boston_lasso_coef != 0)
]

# Remove the intercept from the displayed list.
lasso_selected_names <- setdiff(
  lasso_selected_names,
  "(Intercept)"
)

lasso_selected_names

 [1] "zn"      "indus"   "chas"    "nox"     "rm"      "dis"     "rad"    
 [8] "ptratio" "lstat"   "medv"   

length(lasso_selected_names)

[1] 10

# Variables selected by best subset selection
subset_selected_names <- setdiff(
  names(selected_subset_coefficients),
  "(Intercept)"
)

subset_selected_names

[1] "zn"      "nox"     "rm"      "dis"     "rad"     "tax"     "ptratio"
[8] "lstat"   "medv"   

length(subset_selected_names)

[1] 9

The final model does not have to use every predictor. In this analysis, the lasso selected 10 predictors, and the best subset model selected 9 predictors. Using fewer predictors can make the model simpler while still giving good prediction results.