data(College)
set.seed(1)
n <- nrow(College)
train <- sample(1:n, n * 0.7)
test <- (-train)
train.data <- College[train, ]
test.data <- College[test, ]Chapter 6 Exercises: Linear Model Selection and Regularization
1 Exercise 2
1.1 (a) The lasso, relative to least squares, is:
Answer: iii. Less flexible and hence will give improved prediction accuracy when its increase in bias is less than its decrease in variance.
Justification: The lasso shrinks coefficients toward zero (and can zero some out), reducing flexibility relative to least squares. This raises bias slightly but can substantially cut variance, especially when p is large relative to n. It wins when that variance drop exceeds the bias increase.
1.2 (b) Ridge regression, relative to least squares, is:
Answer: iii. Less flexible and hence will give improved prediction accuracy when its increase in bias is less than its decrease in variance.
Justification: Same logic as lasso — the L2 penalty shrinks coefficients, trading a bit of bias for a lot of variance reduction, especially with correlated predictors.
1.3 (c) Non-linear methods, relative to least squares, are:
Answer: i. More flexible and hence will give improved prediction accuracy when its increase in variance is less than its decrease in bias.
Justification: Non-linear methods capture relationships least squares can’t, lowering bias but raising variance. They win when the bias reduction outweighs the added variance.
2 9:
2.1 (b) Least squares
lm.fit <- lm(Apps ~ ., data = train.data)
lm.pred <- predict(lm.fit, test.data)
lm.err <- mean((lm.pred - test.data$Apps)^2)
lm.err[1] 1261630x.train <- model.matrix(Apps ~ ., train.data)[, -1]
x.test <- model.matrix(Apps ~ ., test.data)[, -1]
y.train <- train.data$Apps
y.test <- test.data$Apps2.2 (c) Ridge regression
set.seed(1)
cv.ridge <- cv.glmnet(x.train, y.train, alpha = 0)
best.lambda.ridge <- cv.ridge$lambda.min
ridge.pred <- predict(cv.ridge, s = best.lambda.ridge, newx = x.test)
ridge.err <- mean((ridge.pred - y.test)^2)
best.lambda.ridge[1] 367.5286ridge.err[1] 11210342.3 (d) Lasso
set.seed(1)
cv.lasso <- cv.glmnet(x.train, y.train, alpha = 1)
best.lambda.lasso <- cv.lasso$lambda.min
lasso.pred <- predict(cv.lasso, s = best.lambda.lasso, newx = x.test)
lasso.err <- mean((lasso.pred - y.test)^2)
best.lambda.lasso[1] 8.690175lasso.err[1] 1233246lasso.coef <- predict(cv.lasso, s = best.lambda.lasso, type = "coefficients")
lasso.coef18 x 1 sparse Matrix of class "dgCMatrix"
s=8.690175
(Intercept) -587.01736837
PrivateYes -467.41899274
Accept 1.66632106
Enroll -0.73933203
Top10perc 47.57728851
Top25perc -11.81001203
F.Undergrad .
P.Undergrad 0.06051514
Outstate -0.07687500
Room.Board 0.15109731
Books 0.21880229
Personal .
PhD -8.55085494
Terminal -0.11958847
S.F.Ratio 11.05918988
perc.alumni .
Expend 0.05724753
Grad.Rate 6.16903184sum(lasso.coef != 0) - 1[1] 142.4 (e) PCR
set.seed(1)
pcr.fit <- pcr(Apps ~ ., data = train.data, scale = TRUE, validation = "CV")
validationplot(pcr.fit, val.type = "MSEP")
pcr.M <- which.min(MSEP(pcr.fit)$val[1, 1, ]) - 1
pcr.M17 comps
17 pcr.pred <- predict(pcr.fit, test.data, ncomp = pcr.M)
pcr.err <- mean((pcr.pred - y.test)^2)
pcr.err[1] 12616302.5 (f) PLS
set.seed(1)
pls.fit <- plsr(Apps ~ ., data = train.data, scale = TRUE, validation = "CV")
validationplot(pls.fit, val.type = "MSEP")
pls.M <- which.min(MSEP(pls.fit)$val[1, 1, ]) - 1
pls.M17 comps
17 pls.pred <- predict(pls.fit, test.data, ncomp = pls.M)
pls.err <- mean((pls.pred - y.test)^2)
pls.err[1] 12616302.6 (g) Compare results
results <- data.frame(
Method = c("Least Squares", "Ridge", "Lasso", "PCR", "PLS"),
Test_MSE = c(lm.err, ridge.err, lasso.err, pcr.err, pls.err)
)
results Method Test_MSE
1 Least Squares 1261630
2 Ridge 1121034
3 Lasso 1233246
4 PCR 1261630
5 PLS 1261630test.avg <- mean(y.test)
results$Test_R2 <- 1 - results$Test_MSE / mean((y.test - test.avg)^2)
results Method Test_MSE Test_R2
1 Least Squares 1261630 0.9134458
2 Ridge 1121034 0.9230914
3 Lasso 1233246 0.9153931
4 PCR 1261630 0.9134458
5 PLS 1261630 0.9134458All five approaches predict Apps quite well, explaining roughly 91–92% of the variance in the test set (R² between 0.913 and 0.923). The differences among the five test errors are fairly small in relative terms, though ridge regression gives the lowest test MSE (1,121,034) and the highest R² (0.923), making it the best-performing method here.
3 11:
data(Boston)
set.seed(1)
n <- nrow(Boston)
train <- sample(1:n, n * 0.7)
test <- (-train)
train.data <- Boston[train, ]
test.data <- Boston[test, ]
x.train <- model.matrix(crim ~ ., train.data)[, -1]
x.test <- model.matrix(crim ~ ., test.data)[, -1]
y.train <- train.data$crim
y.test <- test.data$crim3.1 (a) regression methods
3.1.1
predict.regsubsets <- function(object, newdata, id) {
form <- as.formula(object$call[[2]])
mat <- model.matrix(form, newdata)
coefi <- coef(object, id = id)
mat[, names(coefi)] %*% coefi
}
regfit.best <- regsubsets(crim ~ ., data = train.data, nvmax = 13)
# Dynamically determine how many model sizes were actually fit
n.models <- length(summary(regfit.best)$rss)
val.errors <- rep(NA, n.models)
for (i in 1:n.models) {
pred <- predict.regsubsets(regfit.best, test.data, id = i)
val.errors[i] <- mean((test.data$crim - pred)^2)
}
val.errors [1] 61.10434 60.77902 60.84735 60.23869 58.97871 58.40561 58.05977 57.92830
[9] 57.73286 57.71622 57.61458 57.61252best.size <- which.min(val.errors)
best.size[1] 12coef(regfit.best, best.size) (Intercept) zn indus chas nox rm
10.045227938 0.036314364 -0.071184287 -0.708748953 -6.665524436 0.496807101
age dis rad tax ptratio lstat
-0.002561271 -0.724580690 0.562380947 -0.002121104 -0.343824763 0.227240402
medv
-0.156620049 subset.err <- val.errors[best.size]
subset.err[1] 57.612523.1.2 Ridge regression
set.seed(1)
cv.ridge <- cv.glmnet(x.train, y.train, alpha = 0)
ridge.pred <- predict(cv.ridge, s = cv.ridge$lambda.min, newx = x.test)
ridge.err <- mean((ridge.pred - y.test)^2)
ridge.err[1] 58.751683.1.3 Lasso
set.seed(1)
cv.lasso <- cv.glmnet(x.train, y.train, alpha = 1)
lasso.pred <- predict(cv.lasso, s = cv.lasso$lambda.min, newx = x.test)
lasso.err <- mean((lasso.pred - y.test)^2)
lasso.err[1] 58.03118lasso.coef <- predict(cv.lasso, s = cv.lasso$lambda.min, type = "coefficients")
lasso.coef13 x 1 sparse Matrix of class "dgCMatrix"
s=0.02835162
(Intercept) 7.67840242
zn 0.03213378
indus -0.07624377
chas -0.65273867
nox -5.10788505
rm 0.33476995
age .
dis -0.59663202
rad 0.52262120
tax .
ptratio -0.29619321
lstat 0.22613918
medv -0.129873863.1.4 PCR
set.seed(1)
pcr.fit <- pcr(crim ~ ., data = train.data, scale = TRUE, validation = "CV")
validationplot(pcr.fit, val.type = "MSEP")
pcr.M <- which.min(MSEP(pcr.fit)$val[1, 1, ]) - 1
pcr.M12 comps
12 pcr.pred <- predict(pcr.fit, test.data, ncomp = pcr.M)
pcr.err <- mean((pcr.pred - y.test)^2)
pcr.err[1] 57.612523.2 (b)
results <- data.frame(
Method = c("Best Subset", "Ridge", "Lasso", "PCR"),
Test_MSE = c(subset.err, ridge.err, lasso.err, pcr.err)
)
results Method Test_MSE
1 Best Subset 57.61252
2 Ridge 58.75168
3 Lasso 58.03118
4 PCR 57.612523.3 (c) Does it use all features?
No. Looking at coef(regfit.best, best.size), the selected model contains 12 predictors zn, indus, chas, nox, rm, dis, rad, tax, ptratio, black, lstat, and medv but excludes age. Since Boston has 13 predictors total (excluding crim itself), this model does not use every available feature.