Assignment #5

Question 2

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

(a)

The lasso, relative to least squares, is: i. More flexible and hence will give improved prediction accuracy when its increase in bias is less than its decrease in variance. ii. More flexible and hence will give improved prediction accuracy when its increase in variance is less than its decrease in bias. iii. Less flexible and hence will give improved prediction accuracy when its increase in bias is less than its decrease in variance. iv. Less flexible and hence will give improved prediction accuracy when its increase in variance is less than its decrease in bias.

The correct answer is III. Lasso allows for a reduction in variance in return for a small increase in bias which leads to more accurate predictions.

(b)

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

The correct answer is III. Ridge regression also performs well when trading a small increase in bias for a large decrease in variance. Works best when least squares estimates have high variance.

##(c) >Repeat (a) for non-linear methods relative to least squares.

The correct answer is II. A non linear method sacrifices low variaince of least squares to better capture the curvature of the model. This lenience towards linearity necessary for least squares makes the model more flexible.

Question 9

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

(a)

Split the data set into a training set and a test set.

library(ISLR2)
set.seed(40)

data(College)

train <- sample(1:nrow(College), nrow(College) / 2)
test <- (-train)

College_train <- College[train, ]
College_test <- College[test, ]

(b)

Fit a linear model using least squares on the training set, and report the test error obtained.

lm_fit <- lm(Apps ~ ., data = College_train)

# predicting on test set
lm_pred <- predict(lm_fit, newdata = College_test)

# computing test error
lm_test_mse <- mean((College_test$Apps - lm_pred)^2)
lm_test_mse

[1] 1202848

(c)

Fit a ridge regression model on the training set, with λ chosen by cross-validation. Report the test error obtained.

library(glmnet)

# build model matrices
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

# fit ridge regression
grid <- 10^seq(10, -2, length = 100)
ridge_mod <- glmnet(x_train, y_train, alpha = 0, lambda = grid, thresh = 1e-12)

# choose lambda by cross-validation
set.seed(40)
cv_ridge <- cv.glmnet(x_train, y_train, alpha = 0)
best_lambda_ridge <- cv_ridge$lambda.min
best_lambda_ridge

[1] 393.1805

# predict on the test set
ridge_pred <- predict(ridge_mod, s = best_lambda_ridge, newx = x_test)

# test error
ridge_test_mse <- mean((ridge_pred - y_test)^2)
ridge_test_mse

[1] 1191184

(d)

Fit a lasso model on the training set, with λ chosen by crossvalidation. Report the test error obtained, along with the number of non-zero coefficient estimates.

# fit lasso across
lasso_mod <- glmnet(x_train, y_train, alpha = 1, lambda = grid, thresh = 1e-12)

# choose lambda by cross-validation
set.seed(1)
cv_lasso <- cv.glmnet(x_train, y_train, alpha = 1)
best_lambda_lasso <- cv_lasso$lambda.min
best_lambda_lasso

[1] 28.39081

# predict on the test set
lasso_pred <- predict(lasso_mod, s = best_lambda_lasso, newx = x_test)

# test error
lasso_test_mse <- mean((lasso_pred - y_test)^2)
lasso_test_mse

[1] 1193408

# number of non-zero coefficient estimates
lasso_coef <- predict(lasso_mod, type = "coefficients", s = best_lambda_lasso)
lasso_coef

18 x 1 sparse Matrix of class "dgCMatrix"
               s=28.39081
(Intercept) -378.87942212
PrivateYes  -446.76851885
Accept         1.57591196
Enroll        -0.55989641
Top10perc     32.07488689
Top25perc      .         
F.Undergrad    .         
P.Undergrad    .         
Outstate      -0.06546089
Room.Board     0.04003995
Books          0.12506277
Personal       0.07108963
PhD           -2.94902010
Terminal      -2.12899196
S.F.Ratio      2.25909544
perc.alumni    .         
Expend         0.05124086
Grad.Rate      4.24681410

sum(lasso_coef != 0) - 1  # subtract 1 to exclude the intercept

[1] 13

(e)

Fit a PCR model on the training set, with M chosen by crossvalidation. Report the test error obtained, along with the value of M selected by cross-validation.

library(pls)

set.seed(40)
pcr_fit <- pcr(Apps ~ ., data = College_train, scale = TRUE, validation = "CV")

# cross validation to find M
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
CV            4137     4015     2269     2280     2125     1918
adjCV         4137     4017     2263     2277     2112     1906
       6 comps  7 comps  8 comps  9 comps  10 comps  11 comps  12 comps
CV        1896     1888     1883     1838      1825      1829      1834
adjCV     1886     1878     1876     1828      1814      1819      1825
       13 comps  14 comps  15 comps  16 comps  17 comps
CV         1854      1871      1441      1301      1209
adjCV      1845      1868      1413      1286      1197

TRAINING: % variance explained
      1 comps  2 comps  3 comps  4 comps  5 comps  6 comps  7 comps
X      31.166    55.85    63.03    68.88    74.52    79.72    83.52
Apps    6.236    71.84    71.90    75.82    80.56    81.95    82.24
      8 comps  9 comps  10 comps  11 comps  12 comps  13 comps
X       86.84    89.97     92.51     94.84     96.76     97.85
Apps    82.24    83.57     84.01     84.06     84.06     84.06
      14 comps  15 comps  16 comps  17 comps
X        98.66     99.35     99.83    100.00
Apps     84.06     93.05     93.66     94.14

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

# predict on the test set using M = 1
pcr_pred <- predict(pcr_fit, College_test, ncomp = 17)

# test error
pcr_test_mse <- mean((pcr_pred - College_test$Apps)^2)
pcr_test_mse

[1] 1202848

(f)

Fit a PLS model on the training set, with M chosen by crossvalidation. Report the test error obtained, along with the value of M selected by cross-validation.

set.seed(40)
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
CV            4137     2110     1932     1717     1556     1305
adjCV         4137     2103     1924     1704     1523     1283
       6 comps  7 comps  8 comps  9 comps  10 comps  11 comps  12 comps
CV        1262     1253     1242     1219      1213      1211      1210
adjCV     1248     1239     1228     1206      1200      1199      1198
       13 comps  14 comps  15 comps  16 comps  17 comps
CV         1210      1209      1209      1209      1209
adjCV      1198      1197      1197      1197      1197

TRAINING: % variance explained
      1 comps  2 comps  3 comps  4 comps  5 comps  6 comps  7 comps
X       25.39    46.27    61.15    63.42    66.68    71.95    76.75
Apps    76.54    82.36    87.02    91.86    93.61    93.89    93.93
      8 comps  9 comps  10 comps  11 comps  12 comps  13 comps
X       79.82    81.84     85.78     89.32     90.59     92.63
Apps    94.00    94.09     94.11     94.12     94.14     94.14
      14 comps  15 comps  16 comps  17 comps
X        94.17     96.28     99.12    100.00
Apps     94.14     94.14     94.14     94.14

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

pls_pred <- predict(pls_fit, College_test, ncomp = 1)

pls_test_mse <- mean((pls_pred - College_test$Apps)^2)
pls_test_mse

[1] 2735408

(g)

Comment on the results obtained. How accurately can we predict the number of college applications received? Is there much difference among the test errors resulting from these five approaches?

sqrt(lm_test_mse)

[1] 1096.744

sqrt(ridge_test_mse)

[1] 1091.414

sqrt(lasso_test_mse)

[1] 1092.432

sqrt(pcr_test_mse)

[1] 1096.744

sqrt(pls_test_mse)

[1] 1653.907

Across all 5 approaches, the test errors are very close to one another, such as PCR and PLS. Taking the square root of each mean square error, there is a range of prediction error between 1000-1700 applications. Overall, there is not much difference in the 5 approaches, with only ridge and lasso having a bit of improvement.

Question 11

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

(a)

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.

set.seed(40)
train <- sample(1:nrow(Boston), nrow(Boston) / 2)
test <- (-train)

Boston_train <- Boston[train, ]
Boston_test <- Boston[test, ]

Ridge Regression

x_train <- model.matrix(crim ~ ., data = Boston_train)[, -1]
y_train <- Boston_train$crim
x_test <- model.matrix(crim ~ ., data = Boston_test)[, -1]
y_test <- Boston_test$crim

grid <- 10^seq(10, -2, length = 100)
ridge_mod <- glmnet(x_train, y_train, alpha = 0, lambda = grid, thresh = 1e-12)

set.seed(40)
cv_ridge <- cv.glmnet(x_train, y_train, alpha = 0)
best_lambda_ridge <- cv_ridge$lambda.min

ridge_pred <- predict(ridge_mod, s = best_lambda_ridge, newx = x_test)
ridge_test_mse <- mean((ridge_pred - y_test)^2)
ridge_test_mse

[1] 55.75216

Lasso

lasso_mod <- glmnet(x_train, y_train, alpha = 1, lambda = grid, thresh = 1e-12)

set.seed(40)
cv_lasso <- cv.glmnet(x_train, y_train, alpha = 1)
best_lambda_lasso <- cv_lasso$lambda.min

lasso_pred <- predict(lasso_mod, s = best_lambda_lasso, newx = x_test)
lasso_test_mse <- mean((lasso_pred - y_test)^2)
lasso_test_mse

[1] 55.16232

lasso_coef <- predict(lasso_mod, type = "coefficients", s = best_lambda_lasso)
lasso_coef

13 x 1 sparse Matrix of class "dgCMatrix"
            s=0.05825362
(Intercept)   4.36505654
zn            0.03464220
indus        -0.06950000
chas         -0.54476018
nox           .         
rm           -0.02524042
age           .         
dis          -0.55421158
rad           0.50916802
tax           .         
ptratio      -0.11561034
lstat         0.14589860
medv         -0.11465293

sum(lasso_coef != 0) - 1

[1] 9

PCR

set.seed(40)
pcr_fit <- pcr(crim ~ ., data = Boston_train, scale = TRUE, validation = "CV")
summary(pcr_fit)

Data:   X dimension: 253 12 
    Y dimension: 253 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
CV           7.902    6.124    6.138    5.799    5.663    5.621
adjCV        7.902    6.120    6.135    5.794    5.657    5.616
       6 comps  7 comps  8 comps  9 comps  10 comps  11 comps  12 comps
CV       5.627    5.497    5.441    5.427     5.419     5.430     5.357
adjCV    5.623    5.480    5.434    5.419     5.413     5.423     5.350

TRAINING: % variance explained
      1 comps  2 comps  3 comps  4 comps  5 comps  6 comps  7 comps
X       52.84    66.94    74.96    82.32    88.39    91.69    93.80
crim    40.60    40.67    47.63    49.51    50.58    50.58    52.93
      8 comps  9 comps  10 comps  11 comps  12 comps
X       95.72    97.25     98.64     99.59    100.00
crim    53.84    54.19     54.43     54.73     55.85

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

cv_rmsep_pcr <- RMSEP(pcr_fit)$val[1, 1, ]
best_M_pcr <- which.min(cv_rmsep_pcr) - 1
best_M_pcr

12 comps 
      12 

pcr_pred <- predict(pcr_fit, Boston_test, ncomp = best_M_pcr)
pcr_test_mse <- mean((pcr_pred - Boston_test$crim)^2)
pcr_test_mse

[1] 54.68814

The PCR model performed the lowest of the three with a MSE of 54.59. Since Boston has 12 predictors, using all 12 means no actual dimension reduction took place.

(b)

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, crossvalidation, or some other reasonable alternative, as opposed to using training error.

Comparing the three approaches, ridge produced 55.75, the lasso 55.16, and PCR 54.69. However, the result from PCR can be misleading since M=12 is using all 12 predictors. Using the lasso instead, it produced almost the same MSE will only using 3 of the 12 predictors, which is why it is the best model to use in this case.

(c)

Does your chosen model involve all of the features in the data set? Why or why not?

No, the lasso model does not involve all features of the data set as the model only retained 9 coefficients that weren’t zero and excluded the other 3 entirely. This is due to the lasso penalty forcing the coefficient estimates to be exactly zero when the tuning parameter is large enough.