Conceptual

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 felible 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 flecible and hence will give improved prediction accuracy when its increase in variance is less than its decrease in bias.
  
  The correct option is iii. This is because the lasso shrinks the coefficient estimates towards zero or sets some to zero. With this, the lasso has less flexibility than a least squares. Also, since the lasso is a less flexible model it also has higher bias but lower variance of the variables in a model.

(b) Repear (a) for ridge regression relative to least squares.
  
  i. More flexible and hence will give improved prediction accuracy when its increase in bias is less than its decrease in variance.
  ii. More felible 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 flecible and hence will give improved prediction accuracy when its increase in variance is less than its decrease in bias.
  
  The correct option is iii. This is because ridge regression's advantage over least squares is in the fact that ridge regression has a bias-variance trade-off. As λ increase, the flexibility of the ridge regression fit decreases, leading to decreased variance but increased bias.
  
(c) Repeat (a) for non-linear methods relative to least squares.

  i. More flexible and hence will give improved prediction accuracy when its increase in bias is less than its decrease in variance.
  ii. More felible 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 flecible and hence will give improved prediction accuracy when its increase in variance is less than its decrease in bias.
  
  The correct option is ii. This is because non-linear models allow for more flexibility in relationships between predictors. Because of this increase flexibility, that means the model will also have an increase in variance and a decrease in bias.

Applied

Question 9

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

# install.packages("glmnet")
# install.packages("pls")
# install.packages("leaps")
# install.packages("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
library(leaps)
library(ISLR2)
college_data <- College
  1. Split the data set into a training set and a test set.
set.seed(1)
train <- sample(1:nrow(college_data), nrow(college_data) / 2)
test <- (-train)
  1. Fit a linear model using least squares on the training set, and report the test error obtained.
lm_college <- lm(Apps ~ .,
                 data = college_data,
                 subset = train)

lm_college_preds <- predict(lm_college, college_data[test, ])

mean((college_data$Apps[test] - lm_college_preds) ^ 2)
## [1] 1135758
  1. Fit a ridge regression model on the training set, with λ chosen by cross-validation. Report the test error obtained.
x <- model.matrix(Apps ~ ., college_data)[, -1]
y <- college_data$Apps

set.seed(1)
cv.out <- cv.glmnet(x[train, ], y[train], alpha = 0)
bestlam <- cv.out$lambda.min
bestlam
## [1] 405.8404
ridge.mod <- glmnet(x[train, ], y[train], alpha = 0, lambda = bestlam)

ridge.prediction <- predict(ridge.mod, s = bestlam, newx = x[test, ])

mean((ridge.prediction - y[test]) ^ 2)
## [1] 976897.6
  1. Fit a lasso model on the training set, with λ chosen by cross-validation. Report the test error obtained, along with the number of non-zero coefficient estimates.
set.seed(1)
lasso_cv_out <- cv.glmnet(x[train, ], y[train], alpha = 1)
lasso_bestlam <- lasso_cv_out$lambda.min
lasso_bestlam
## [1] 1.97344
lasso.mod <- glmnet(x[train, ], y[train], alpha = 1, lambda = lasso_bestlam)

lasso.predictions <- predict(lasso.mod, s = lasso_bestlam, newx = x[test, ])

mean((lasso.predictions - y[test]) ^ 2)
## [1] 1116402
lasso.coef <- predict(lasso.mod, s = lasso_bestlam, type = "coefficients")
sum(lasso.coef != 0) - 1
## [1] 17
  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(2)

pcr.fit <- pcr(Apps ~ ., data = college_data, subset = 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     4032     2359     2362     2055     1924     1881
## adjCV         4288     4035     2355     2358     1971     1911     1875
##        7 comps  8 comps  9 comps  10 comps  11 comps  12 comps  13 comps
## CV        1882     1890     1834      1799      1798      1799      1794
## adjCV     1877     1887     1829      1790      1791      1792      1787
##        14 comps  15 comps  16 comps  17 comps
## CV         1795      1764      1348      1278
## adjCV      1790      1737      1332      1265
## 
## 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.89
## (389 observations deleted due to missingness)
validationplot(pcr.fit, val.type = "MSEP")

pcr.pred <- predict( pcr.fit, college_data[test, ], ncomp = 15)

mean((pcr.pred - y[test]) ^ 2)
## [1] 1268629
  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 = college_data, subset = 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
## (389 observations deleted due to missingness)
validationplot(pls.fit, val.type = "MSEP")

pls.pred <- predict(pls.fit, college_data[test, ], ncomp = 6)

mean((pls.pred - y[test]) ^ 2)
## [1] 1066991
  1. 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?
Approach MSE
Linear model using least squares 1,135,758
Ridge regression model 976,897.6
Lasso model 1,116,402
PCR model 1,268,629
PLS model 1,066,991

From the MSE of the five different approaches we see that the ridge regression model performs the best. However, to see how accurately predict the number of applications we have to take the root mean squared error (RMSE).

sqrt(1135758)
## [1] 1065.719
sqrt(976897.6)
## [1] 988.3813
sqrt(1116402)
## [1] 1056.599
sqrt(1268629)
## [1] 1126.334
sqrt(1066991)
## [1] 1032.953
Approach MSE RMSE
Linear model using least squares 1,135,758 1065.72
Ridge regression model 976,897.6 988.38
Lasso model 1,116,402 1056.60
PCR model 1,268,629 1126.33
PLS model 1,066,991 1032.95

Based on the RSME, it reconfirms that ridge regression performs the best. We can say this because the ridge regression model has the lower difference in the number of applications predicted.

While ridge is the best performing model, on the other end of the spectrum we see that the PCR model is the one performing the worse in predicting using unseen data.

Question 11

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

boston_data <- Boston
  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.
# Ridge Regression

boston_x <- model.matrix(crim ~ ., boston_data)[, -1]
boston_y <- boston_data$crim

boston_train <- sample(1:nrow(boston_x), nrow(boston_x) / 2)
boston_test <- (-boston_train)
y.test <- boston_y[boston_test]

set.seed(1)
boston_cv_out <- cv.glmnet(boston_x[boston_train, ], boston_y[boston_train], alpha = 0)
boston_best_lam <- boston_cv_out$lambda.min
boston_best_lam
## [1] 0.4700022
boston_ridge <- glmnet(boston_x[boston_train, ], boston_y[boston_train], alpha = 0,
                       lambda = boston_best_lam)

boston_ridge.prediction <- predict(boston_ridge, s = boston_best_lam, newx = boston_x[boston_test, ])

boston_mse <- mean((boston_ridge.prediction - boston_y[boston_test]) ^ 2)
boston_mse
## [1] 72.15903
sqrt(boston_mse)
## [1] 8.494647

From the ridge regression model on the Boston data, we can see that the model is off by about 4.74 crime rate units when predicting on unseen Boston data.

# PCR

set.seed(2)
boston_pcr <- pcr(crim ~ ., data = boston_data, subset = boston_train, scale = TRUE, validation = "CV")
summary(boston_pcr)
## 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  6 comps
## CV           6.282    4.768    4.754    4.381    4.313    4.238    4.205
## adjCV        6.282    4.765    4.751    4.352    4.307    4.234    4.200
##        7 comps  8 comps  9 comps  10 comps  11 comps  12 comps
## CV       4.067    4.022    4.005     4.004     4.005     3.943
## adjCV    4.059    4.013    3.992     3.997     3.997     3.935
## 
## TRAINING: % variance explained
##       1 comps  2 comps  3 comps  4 comps  5 comps  6 comps  7 comps  8 comps
## X       51.11    65.01    73.40    80.91    86.82    90.21    92.94    95.00
## crim    43.20    43.72    52.96    54.30    55.80    56.82    59.61    60.96
##       9 comps  10 comps  11 comps  12 comps
## X       96.83     98.51     99.51    100.00
## crim    61.81     62.01     62.11     63.34
## (253 observations deleted due to missingness)
validationplot(boston_pcr, val.type = "MSEP")

#PCR continued

boston_pcr_pred <- predict(boston_pcr, boston_x[boston_test, ], ncomp = 7)
pcr_mse <- mean((boston_pcr_pred - y.test) ^ 2)
pcr_mse
## [1] 72.67164
sqrt(pcr_mse)
## [1] 8.524766

From the PCR model on the Boston data, we can see that the model is off by about 4.85 crime rate units when predicting on unseen Boston data.

# PLS

set.seed(1)
boston_plsr <- plsr(crim ~ ., data = boston_data, subset = boston_train, scale = TRUE, validation = "CV")
summary(boston_plsr)
## Data:    X dimension: 253 12 
##  Y dimension: 253 1
## Fit method: kernelpls
## 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           6.282    4.573    4.056    3.973    3.939    3.916    3.914
## adjCV        6.282    4.571    4.053    3.969    3.933    3.911    3.908
##        7 comps  8 comps  9 comps  10 comps  11 comps  12 comps
## CV       3.903    3.901    3.900     3.899     3.899     3.899
## adjCV    3.897    3.895    3.894     3.893     3.893     3.893
## 
## TRAINING: % variance explained
##       1 comps  2 comps  3 comps  4 comps  5 comps  6 comps  7 comps  8 comps
## X       50.77    59.39    68.34    75.92    81.06    84.34    89.23    92.49
## crim    47.64    59.38    61.42    62.54    62.85    63.14    63.26    63.33
##       9 comps  10 comps  11 comps  12 comps
## X       95.12     96.61     98.26    100.00
## crim    63.34     63.34     63.34     63.34
## (253 observations deleted due to missingness)
boston_plsr_pred <- predict(boston_plsr, boston_x[boston_test, ], ncomp = 10)
plsr_mse <- mean((boston_plsr_pred - y.test) ^ 2)
plsr_mse
## [1] 70.86988
sqrt(plsr_mse)
## [1] 8.418425

From the PLS model on the Boston data, we can see that the model is off by about 4.87 crime rate units when predicting on unseen Boston data.

Out of the three models, we see that the ridge regression model is the most accurate when predicting the per capita crime rate in the Boston data set.

  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.

From the previous section, the model that I would suggest perform best on the data set are ridge regression.

I would recommend this model because the Boston data has multiple predictors that have collinearity but ridge handles the risk of collinearity by shrinking the coefficients of each predictor and lowering the variance of the model.

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

Yes, the ridge regression model does involve all the features of the data set. All features are included because a ridge regression shrinks the coefficients close to zero but never zero so nothing will be removed from the model.