Question 2

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

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.

Solution: iii. Less flexible and hence will give improved prediction accuracy when its increase in bias is less than its decrease in variance. The Lasso reduces flexibility and is therefore a more restrictive model, it can potentially decrease the variance and overfitting.

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

    Solution: iii. Less flexible and hence will give improved prediction accuracy when its increase in bias is less than its decrease in variance.The ridge regression is a little more flexible than the lasso, however, it is more restrictive than least squares.

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

    Solution: 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 tend to be more flexible than least squares method.

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)
## Warning: package 'ISLR2' was built under R version 4.1.3
library(glmnet)
## Warning: package 'glmnet' was built under R version 4.1.3
## Loading required package: Matrix
## Warning: package 'Matrix' was built under R version 4.1.3
## Loaded glmnet 4.1-3
library(pls)
## Warning: package 'pls' was built under R version 4.1.3
## 
## Attaching package: 'pls'
## The following object is masked from 'package:stats':
## 
##     loadings
library(leaps)
## Warning: package 'leaps' was built under R version 4.1.3
set.seed(1)
train_sub=sample(c(TRUE,FALSE),nrow(College),rep=TRUE)
test_sub=(!train_sub)
Col.train=College[train_sub,,drop=F]
Col.test=College[test_sub,,drop=F]
  1. Fit a linear model using least squares on the training set, and report the test error obtained.
lm.fit = lm(Apps~., data=Col.train)
lm.pred = predict(lm.fit, Col.test, type="response")
mean((Col.test[, "Apps"] - lm.pred)^2)
## [1] 984743.1

Solution: The test error obtained is 984743.

  1. Fit a ridge regression model on the training set, with λ chosen by cross-validation. Report the test error obtained.
train.mat = model.matrix(Apps~., data=Col.train)
test.mat = model.matrix(Apps~., data=Col.test)

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

ridge.mod = cv.glmnet(train.mat, Col.train[, "Apps"], alpha=0, lambda=grid, thresh=1e-12)

ridge.pred = predict(ridge.mod, newx=test.mat, s=ridge.mod$lambda.min)
mean((Col.test[, "Apps"] - ridge.pred)^2)
## [1] 984731.2

Solution: The test error obtained is 984731.

  1. 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.
lasso.mod = cv.glmnet(train.mat, Col.train[, "Apps"],
alpha=1, lambda=grid, thresh=1e-12)
lambda.best = lasso.mod$lambda.min
lambda.best
## [1] 0.01
lasso.pred = predict(lasso.mod, newx=test.mat, s=lasso.mod$lambda.min)
mean((Col.test[, "Apps"] - lasso.pred)^2)
## [1] 984715.2
lasso.mod = glmnet(model.matrix(Apps~., data=College),
College[, "Apps"], alpha=1)
predict(lasso.mod, s=lambda.best, type="coefficients")
## 19 x 1 sparse Matrix of class "dgCMatrix"
##                        s1
## (Intercept) -471.39372052
## (Intercept)    .         
## PrivateYes  -491.04485137
## Accept         1.57033288
## Enroll        -0.75961467
## Top10perc     48.14698892
## Top25perc    -12.84690695
## F.Undergrad    0.04149116
## P.Undergrad    0.04438973
## Outstate      -0.08328388
## Room.Board     0.14943472
## Books          0.01532293
## Personal       0.02909954
## PhD           -8.39597537
## Terminal      -3.26800340
## S.F.Ratio     14.59298267
## perc.alumni   -0.04404771
## Expend         0.07712632
## Grad.Rate      8.28950241

Solution: The test error obtained is 998042.

  1. 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.
pcr.fit = pcr(Apps~., data=Col.train, scale=T, validation="CV")
validationplot(pcr.fit, val.type="MSEP")

pcr.pred = predict(pcr.fit, Col.test, ncomp=10)
mean((Col.test[, "Apps"] - pcr.pred)^2)
## [1] 1682909

Solution: The test error obtained is 1682909.

  1. 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.
pls.fit = plsr(Apps~., data=Col.train, scale=T, validation="CV")
validationplot(pls.fit, val.type="MSEP")

pls.pred = predict(pls.fit, Col.test, ncomp=10)
mean((Col.test[, "Apps"] - pls.pred)^2)
## [1] 994703.4

Solution: The test error obtained is 994703.

  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?

Summary of results obtained: Linear Model: 984743 Ridge Regression: 984731 Lasso: 998042 PCR: 1682909 PLS: 994703

We can see that the results for our linear model, ridge regression, lasso, and PLS are not much different. PCR, compared to our other models, has a worse performance.

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.

First, we select the best subset

attach(Boston)
predict.regsubsets <- function(object, newdata, id, ...) {
  formula <- as.formula(object$call[[2]])
  model_matrix <- model.matrix(formula, newdata)
  coefi <- coef(object, id = id)
  model_matrix[, names(coefi)] %*% coefi
}
k <- 10
p <- ncol(Boston)-1
folds <- sample(rep(1:k, length = nrow(Boston)))
cv_errors <- matrix(NA, k, p)

for (i in 1:k) {
  best_fit <- regsubsets(crim ~ . , data = Boston[folds!=i,], nvmax = p)
  for (j in 1:p) {
    prediction <- predict(best_fit, Boston[folds==i, ], id = j)
    cv_errors[i,j] <-  mean((Boston$crim[folds==i] - prediction)^2)    
  }
}

cv_rmse <- sqrt(apply(cv_errors, 2, mean))
plot(cv_rmse, type = "b")

summary(best_fit)
## Subset selection object
## Call: regsubsets.formula(crim ~ ., data = Boston[folds != i, ], nvmax = p)
## 12 Variables  (and intercept)
##         Forced in Forced out
## zn          FALSE      FALSE
## indus       FALSE      FALSE
## chas        FALSE      FALSE
## nox         FALSE      FALSE
## rm          FALSE      FALSE
## age         FALSE      FALSE
## dis         FALSE      FALSE
## rad         FALSE      FALSE
## tax         FALSE      FALSE
## ptratio     FALSE      FALSE
## lstat       FALSE      FALSE
## medv        FALSE      FALSE
## 1 subsets of each size up to 12
## Selection Algorithm: exhaustive
##           zn  indus chas nox rm  age dis rad tax ptratio lstat medv
## 1  ( 1 )  " " " "   " "  " " " " " " " " "*" " " " "     " "   " " 
## 2  ( 1 )  " " " "   " "  " " " " " " " " "*" " " " "     "*"   " " 
## 3  ( 1 )  " " " "   " "  " " " " " " " " "*" " " " "     "*"   "*" 
## 4  ( 1 )  "*" " "   " "  " " " " " " "*" "*" " " " "     " "   "*" 
## 5  ( 1 )  "*" " "   " "  "*" " " " " "*" "*" " " " "     " "   "*" 
## 6  ( 1 )  "*" " "   " "  "*" " " " " "*" "*" " " "*"     " "   "*" 
## 7  ( 1 )  "*" " "   " "  "*" " " " " "*" "*" " " "*"     "*"   "*" 
## 8  ( 1 )  "*" "*"   " "  "*" " " " " "*" "*" " " "*"     "*"   "*" 
## 9  ( 1 )  "*" "*"   " "  "*" "*" " " "*" "*" " " "*"     "*"   "*" 
## 10  ( 1 ) "*" " "   "*"  "*" "*" " " "*" "*" "*" "*"     "*"   "*" 
## 11  ( 1 ) "*" "*"   "*"  "*" "*" " " "*" "*" "*" "*"     "*"   "*" 
## 12  ( 1 ) "*" "*"   "*"  "*" "*" "*" "*" "*" "*" "*"     "*"   "*"
which.min(cv_rmse)
## [1] 11
bos.sub.err=(cv_rmse[which.min(cv_rmse)]^2)
bos.sub.err
## [1] 42.99104

Solution: The mean standard error is 42.4191. This includes the following variables: zn, indus, chas, nox, rmc, dis, rad, tax, ptratio, lstat, medv (“age” has been thrown out)

Next, we continue with Lasso

bos_x <- model.matrix(crim ~ . -1, data = Boston)
bos_lasso <- cv.glmnet(bos_x, Boston$crim, type.measure = "mse", alpha = 1)
plot(bos_lasso)

coef(bos_lasso)
## 13 x 1 sparse Matrix of class "dgCMatrix"
##                    s1
## (Intercept) 1.0894283
## zn          .        
## indus       .        
## chas        .        
## nox         .        
## rm          .        
## age         .        
## dis         .        
## rad         0.2643196
## tax         .        
## ptratio     .        
## lstat       .        
## medv        .
bos_lasso_error <- (bos_lasso$cvm[bos_lasso$lambda == bos_lasso$lambda.1se])
bos_lasso_error
## [1] 55.3469

Solution: The mean standard error is 54.88244.

Next, we continue with ridge regression:

bos_x <- model.matrix(crim ~ . -1, data = Boston)
bos_ridge <- cv.glmnet(bos_x, Boston$crim, type.measure = "mse", alpha = 0)
plot(bos_ridge)

coef(bos_ridge)
## 13 x 1 sparse Matrix of class "dgCMatrix"
##                       s1
## (Intercept)  0.713143934
## zn          -0.002994073
## indus        0.028593701
## chas        -0.157594395
## nox          1.825119224
## rm          -0.138551061
## age          0.006029211
## dis         -0.091316045
## rad          0.043640145
## tax          0.001994740
## ptratio      0.068418435
## lstat        0.034406354
## medv        -0.022668857
bos_ridge_error <- (bos_ridge$cvm[bos_ridge$lambda == bos_ridge$lambda.1se])
bos_ridge_error
## [1] 60.03422

Solution: The mean standard error for the ridge regression is 55.14448, which is higher than the Lasso MSE.

We continue with PCR

bos_pcr <- pcr(crim ~ . , data = Boston, scale = TRUE, validation = "CV")
summary(bos_pcr)
## Data:    X dimension: 506 12 
##  Y dimension: 506 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            8.61    7.265    7.265    6.860    6.838    6.785    6.777
## adjCV         8.61    7.262    7.263    6.854    6.835    6.782    6.773
##        7 comps  8 comps  9 comps  10 comps  11 comps  12 comps
## CV       6.645    6.657    6.647     6.630     6.590     6.524
## adjCV    6.640    6.651    6.641     6.624     6.583     6.516
## 
## TRAINING: % variance explained
##       1 comps  2 comps  3 comps  4 comps  5 comps  6 comps  7 comps  8 comps
## X       49.93    63.64    72.94    80.21    86.83    90.26    92.79    94.99
## crim    29.39    29.55    37.39    37.85    38.85    39.23    41.73    41.82
##       9 comps  10 comps  11 comps  12 comps
## X       96.78     98.33     99.48    100.00
## crim    42.12     42.43     43.58     44.93
  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, crossvalidation, or some other reasonable alternative, as opposed to using training error.

Solution: Based on the results, I would propose to use the subset selection model because it has the lowest standard error rate (42.4191) and is therefore the best model to predict the crime rate in Boston.

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

Solution: My chosen model, the subset selection model, does not involve all of the features in the dataset because “age” was thrown out.