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.
Less flexible and hence will give improved prediction accuracy when its increase in bias is less than its decrease in variance.
(b) Repeat (a) for ridge regression relative to least squares. Less flexible and hence will give improved prediction accuracy when its increase in bias is less than its decrease in variance.
(c) Repeat (a) for non-linear methods relative to least squares. More flexible and hence will give improved prediction accuracy when its increase in variance is less than its decrease in bias.
In this exercise, we will predict the number of applications received using the other variables in the College data set.
library(ISLR)
set.seed(1)
(a) Split the data set into a training set and a test set.
train.dim = dim(College)[1]/2
train = sample(1:dim(College)[1], train.dim)
test = -train
CollegeTrain = College[train,]
CollegeTest = 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 = CollegeTrain)
lm.pred = predict(lm.fit, CollegeTest)
mean((CollegeTest[,'Apps']-lm.pred)^2)
## [1] 1135758
(c) Fit a ridge regression model on the training set, with λ chosen by cross-validation. Report the test error obtained.
library(glmnet)
## Warning: package 'glmnet' was built under R version 4.0.4
## Loading required package: Matrix
## Loaded glmnet 4.1-1
train.m = model.matrix(Apps~., data=CollegeTrain)
test.m = model.matrix(Apps~., data=CollegeTest)
grid = 10 ^ seq(4,-2, length=100)
model.ridge = cv.glmnet(train.m, CollegeTrain[,"Apps"], alpha=0, lambda=grid, thresh=1e-12)
lambda.best = model.ridge$lambda.min
lambda.best
## [1] 0.01
ridge.pred = predict(model.ridge, newx = test.m, s=lambda.best)
mean((CollegeTest[,"Apps"]-ridge.pred)^2)
## [1] 1135714
(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.
lasso.model = cv.glmnet(train.m, CollegeTrain[,"Apps"], alpha=1, lambda = grid, thresh=1e-12)
lambda.best = lasso.model$lambda.min
lambda.best
## [1] 0.01
lasso.pred = predict(lasso.model, newx = test.m, s=lambda.best)
mean((CollegeTest[,"Apps"] - lasso.pred)^2)
## [1] 1135660
lasso.model = glmnet(model.matrix(Apps~., data=College), College[,'Apps'], alpha=1)
predict(lasso.model, s=lambda.best, type="coefficients")
## 19 x 1 sparse Matrix of class "dgCMatrix"
## 1
## (Intercept) -471.39372069
## (Intercept) .
## PrivateYes -491.04485135
## Accept 1.57033288
## Enroll -0.75961467
## Top10perc 48.14698891
## Top25perc -12.84690694
## 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
(e) 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.
library(pls)
## Warning: package 'pls' was built under R version 4.0.4
##
## Attaching package: 'pls'
## The following object is masked from 'package:stats':
##
## loadings
fit.pcr = pcr(Apps~., data = CollegeTrain, scale=T, validation = "CV")
validationplot(fit.pcr, val.type = "MSEP")
pred.pcr = predict(fit.pcr, CollegeTest, ncomp = 10)
mean((CollegeTest$Apps-pred.pcr)^2)
## [1] 1723100
(f) 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.
fit.pls = plsr(Apps~., data=CollegeTrain, scale=TRUE, validation = "CV")
validationplot(fit.pls, val.type = "MSEP")
pred.pls = predict(fit.pls, CollegeTest, ncomp=10)
mean((CollegeTest$Apps - pred.pls)^2)
## [1] 1131661
(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? For the most part, it can be seen from all results that the majority of the models have very similar performance in predicting the number of college applications. So there is not much of a difference from the five approaches.
We will now try to predict per capita crime rate in the Boston data set.
library(MASS)
library(leaps)
data(Boston)
set.seed(1)
(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.
Best Subset Selection
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
}
k = 10
folds = sample(1:k, nrow(Boston), replace = TRUE)
cv.errors = matrix(NA, k, 13, dimnames = list(NULL, paste(1:13)))
for (j in 1:k) {
best.fit = regsubsets(crim ~ ., data = Boston[folds != j, ], nvmax = 13)
for (i in 1:13) {
pred = predict(best.fit, Boston[folds == j, ], id = i)
cv.errors[j, i] = mean((Boston$crim[folds == j] - pred)^2)
}
}
mean.cv.errors = apply(cv.errors, 2, mean)
plot(mean.cv.errors, type = "b", xlab = "Number of Variables", ylab = "CV error")
min(mean.cv.errors)
## [1] 42.46014
regfit.best = regsubsets(crim~., data=Boston, nvmax=13)
coef(regfit.best,12)
## (Intercept) zn indus chas nox
## 16.985713928 0.044673247 -0.063848469 -0.744367726 -10.202169211
## rm dis rad tax ptratio
## 0.439588002 -0.993556631 0.587660185 -0.003767546 -0.269948860
## black lstat medv
## -0.007518904 0.128120290 -0.198877768
Lasso
x = model.matrix(crim~., Boston)[,-1]
y = Boston$crim
cv.out = cv.glmnet(x,y, alpha = 1, type.measure = 'mse')
plot(cv.out)
cv.out
##
## Call: cv.glmnet(x = x, y = y, type.measure = "mse", alpha = 1)
##
## Measure: Mean-Squared Error
##
## Lambda Index Measure SE Nonzero
## min 0.051 51 43.11 14.16 11
## 1se 3.376 6 56.89 17.29 1
Ridge
cv.out = cv.glmnet(x, y, alpha = 0, type.measure = 'mse')
cv.out
##
## Call: cv.glmnet(x = x, y = y, type.measure = "mse", alpha = 0)
##
## Measure: Mean-Squared Error
##
## Lambda Index Measure SE Nonzero
## min 0.54 100 43.48 14.33 13
## 1se 74.44 47 57.69 16.76 13
PCR
fit.pcr = pcr(crim~., data=Boston, scale=TRUE, validation="CV")
validationplot(fit.pcr, val.type = 'MSEP')
(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, cross validation, or some other reasonable alternative, as opposed to using training error. The model that comes out on top is the Best Subset Selection because it returns the best results.
(c) Does your chosen model involve all of the features in the data set? Why or why not? The chosen model includes 12 out of the 13 features in the data set.