Assignment 5
Ashley Torres
2026-07-05
Question 2:
The lasso, relative to least squares is less flexible and hence will give improved prediction accuracy when its increase in bias is less than its decrease in variance.
The ridge regression relative to least squares is less flexible and hence will give improved prediction accuracy when its increase in bias is less than its decrease in variance.
Non-linear methods relative to least squares are more flexible and hence will give improved prediction accuracy when its increase in variance is less than its decrease in bias.
Question: 9
a. split the data set into a training set and a test set.
#load libraries
#install packages if not already done so
library(ISLR2)
library(glmnet)
library(pls)
library(ggplot2)
library(DT)
library(leaps)
#load data
data("College")
summary(College)## Private Apps Accept Enroll Top10perc
## No :212 Min. : 81 Min. : 72 Min. : 35 Min. : 1.00
## Yes:565 1st Qu.: 776 1st Qu.: 604 1st Qu.: 242 1st Qu.:15.00
## Median : 1558 Median : 1110 Median : 434 Median :23.00
## Mean : 3002 Mean : 2019 Mean : 780 Mean :27.56
## 3rd Qu.: 3624 3rd Qu.: 2424 3rd Qu.: 902 3rd Qu.:35.00
## Max. :48094 Max. :26330 Max. :6392 Max. :96.00
## Top25perc F.Undergrad P.Undergrad Outstate
## Min. : 9.0 Min. : 139 Min. : 1.0 Min. : 2340
## 1st Qu.: 41.0 1st Qu.: 992 1st Qu.: 95.0 1st Qu.: 7320
## Median : 54.0 Median : 1707 Median : 353.0 Median : 9990
## Mean : 55.8 Mean : 3700 Mean : 855.3 Mean :10441
## 3rd Qu.: 69.0 3rd Qu.: 4005 3rd Qu.: 967.0 3rd Qu.:12925
## Max. :100.0 Max. :31643 Max. :21836.0 Max. :21700
## Room.Board Books Personal PhD
## Min. :1780 Min. : 96.0 Min. : 250 Min. : 8.00
## 1st Qu.:3597 1st Qu.: 470.0 1st Qu.: 850 1st Qu.: 62.00
## Median :4200 Median : 500.0 Median :1200 Median : 75.00
## Mean :4358 Mean : 549.4 Mean :1341 Mean : 72.66
## 3rd Qu.:5050 3rd Qu.: 600.0 3rd Qu.:1700 3rd Qu.: 85.00
## Max. :8124 Max. :2340.0 Max. :6800 Max. :103.00
## Terminal S.F.Ratio perc.alumni Expend
## Min. : 24.0 Min. : 2.50 Min. : 0.00 Min. : 3186
## 1st Qu.: 71.0 1st Qu.:11.50 1st Qu.:13.00 1st Qu.: 6751
## Median : 82.0 Median :13.60 Median :21.00 Median : 8377
## Mean : 79.7 Mean :14.09 Mean :22.74 Mean : 9660
## 3rd Qu.: 92.0 3rd Qu.:16.50 3rd Qu.:31.00 3rd Qu.:10830
## Max. :100.0 Max. :39.80 Max. :64.00 Max. :56233
## Grad.Rate
## Min. : 10.00
## 1st Qu.: 53.00
## Median : 65.00
## Mean : 65.46
## 3rd Qu.: 78.00
## Max. :118.00#set seed
set.seed(1)
#Split into training and test
train <- sample(1:nrow(College), nrow(College)/2)
college.train <- College[train, ]
college.test <- College[-train, ]
y.test <- college.test$Apps
y.train <- college.train$Appsb. Fit a linear model using least squares on the training set, and report the test error obtained.
lm.fit <- lm(Apps~.,data = college.train)
lm.pred <- predict(lm.fit, newdata = college.test)
lm.testError <- mean((lm.pred - y.test)^2)
cat("Error:", lm.testError, "\n")## Error: 1135758c. Fit a ridge regression model on the training set with λ chosen by cross-validation. Report the test error obtained.
x.train <- model.matrix(Apps ~ ., college.train)[,-1]
x.test <- model.matrix(Apps ~ ., college.test)[,-1]
set.seed(1)
cv.ridge <- cv.glmnet(x.train,
y.train,
alpha = 0)
plot(cv.ridge)
ridge.bestlam <- cv.ridge$lambda.min
cat("Chosen lambda:", ridge.bestlam, "\n")
ridge.pred <- predict(cv.ridge,
s = ridge.bestlam,
newx = x.test)
ridge.MSE <- mean((y.test - ridge.pred)^2)
cat("Error:", ridge.MSE, "\n")## Chosen lambda: 405.8404
## Error: 976261.5
d. Fit a lasso model on the training set, with λ chosen by cross validation. Report the test error obtained, along with the number on non-zero coefficient estimates.
set.seed(1)
#fit using cv
cv.out <- cv.glmnet(x.train,
y.train,
alpha = 1)
plot(cv.out) #plots cv error
#choosing lambda
bestlam.lasso <- cv.out$lambda.min
#MSE
lasso.pred <- predict(cv.out,
s = bestlam.lasso,
newx = x.test)
lasso.mse <- mean((lasso.pred - y.test)^2) #test MSE
cat("Lasso MSE:", lasso.mse, "\n") #prints the MSE
#finding the non-zero coefficient estimates
lasso.coef <- predict(cv.out,
type = "coefficients",
s = bestlam.lasso)
lasso.coef
lasso.numcoef <- sum(lasso.coef != 0)-1 #non zero coef excluding intercept
cat("Number of Non-Zero Coefficients:", lasso.numcoef, "\n")## Lasso MSE: 1115901
## 18 x 1 sparse Matrix of class "dgCMatrix"
## s=1.97344
## (Intercept) -7.688896e+02
## PrivateYes -3.127034e+02
## Accept 1.762718e+00
## Enroll -1.318195e+00
## Top10perc 6.482356e+01
## Top25perc -2.081406e+01
## F.Undergrad 7.119149e-02
## P.Undergrad 1.246161e-02
## Outstate -1.049091e-01
## Room.Board 2.088305e-01
## Books 2.926466e-01
## Personal 3.955068e-03
## PhD -1.455463e+01
## Terminal 5.395858e+00
## S.F.Ratio 2.171398e+01
## perc.alumni 5.088260e-01
## Expend 4.824455e-02
## Grad.Rate 7.036148e+00
## Number of Non-Zero Coefficients: 17
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.
set.seed(1)
pcr.fit <- pcr(Apps ~ .,
data = college.train,
scale = TRUE,
validation = "CV")
validationplot(pcr.fit, val.type = "MSEP")
summary(pcr.fit) #summary table shows 17 comps is best## 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 4006 2373 2372 2069 1961 1919
## adjCV 4288 4007 2368 2369 1999 1948 1911
## 7 comps 8 comps 9 comps 10 comps 11 comps 12 comps 13 comps
## CV 1919 1921 1876 1832 1832 1836 1837
## adjCV 1912 1915 1868 1821 1823 1827 1827
## 14 comps 15 comps 16 comps 17 comps
## CV 1853 1759 1341 1270
## adjCV 1850 1733 1326 1257
##
## 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
pcr.pred <- predict(pcr.fit,
college.test,
ncomp = 17)
pcr.mse <- mean((y.test - pcr.pred)^2)
cat("PCR MSE:", pcr.mse, "\n")## PCR MSE: 1135758f. 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.train,
scale = TRUE,
validation = "CV")
validationplot(pls.fit, val.type = "MSEP")
summary(pls.fit) #M = 14## 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
pls.pred <- predict(pls.fit,
college.test,
ncomp = 14)
pls.mse <- mean((y.test - pls.pred)^2)
cat("Partial Least Squares MSE:", pls.mse, "\n")## Partial Least Squares MSE: 1137701g. 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 approached?
The five models showed similar results, however ridge regression produced the lowest mean squared error. Similar results produced from the models indicate a prediction of college applications with reasonable accuracy. Since the difference in results isn’t drastically different, this indicates that no mode out preforms the others, yet ridge regression has a slight advantage in reducing the test error.
Question: 11
We will now try to predict per capita 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.
data set up
summary(Boston)
y <- Boston$crim
x <- model.matrix(crim ~., Boston)[, -1]
set.seed(1)
#splitting the data
train <- sample(1:nrow(Boston), nrow(Boston)/2)
test <- (-train)
boston.train <- Boston[train,]
boston.test <- Boston[-train,]
y.test <- boston.test$crim
y.train <- boston.train$crim## crim zn indus chas
## Min. : 0.00632 Min. : 0.00 Min. : 0.46 Min. :0.00000
## 1st Qu.: 0.08205 1st Qu.: 0.00 1st Qu.: 5.19 1st Qu.:0.00000
## Median : 0.25651 Median : 0.00 Median : 9.69 Median :0.00000
## Mean : 3.61352 Mean : 11.36 Mean :11.14 Mean :0.06917
## 3rd Qu.: 3.67708 3rd Qu.: 12.50 3rd Qu.:18.10 3rd Qu.:0.00000
## Max. :88.97620 Max. :100.00 Max. :27.74 Max. :1.00000
## nox rm age dis
## Min. :0.3850 Min. :3.561 Min. : 2.90 Min. : 1.130
## 1st Qu.:0.4490 1st Qu.:5.886 1st Qu.: 45.02 1st Qu.: 2.100
## Median :0.5380 Median :6.208 Median : 77.50 Median : 3.207
## Mean :0.5547 Mean :6.285 Mean : 68.57 Mean : 3.795
## 3rd Qu.:0.6240 3rd Qu.:6.623 3rd Qu.: 94.08 3rd Qu.: 5.188
## Max. :0.8710 Max. :8.780 Max. :100.00 Max. :12.127
## rad tax ptratio lstat
## Min. : 1.000 Min. :187.0 Min. :12.60 Min. : 1.73
## 1st Qu.: 4.000 1st Qu.:279.0 1st Qu.:17.40 1st Qu.: 6.95
## Median : 5.000 Median :330.0 Median :19.05 Median :11.36
## Mean : 9.549 Mean :408.2 Mean :18.46 Mean :12.65
## 3rd Qu.:24.000 3rd Qu.:666.0 3rd Qu.:20.20 3rd Qu.:16.95
## Max. :24.000 Max. :711.0 Max. :22.00 Max. :37.97
## medv
## Min. : 5.00
## 1st Qu.:17.02
## Median :21.20
## Mean :22.53
## 3rd Qu.:25.00
## Max. :50.00Best Subset
p <- ncol(Boston)-1
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
}
regfit.best <- regsubsets(crim ~.,
data = boston.train,
nvmax = p)
test.errors <- rep(NA, p)
for (i in 1:p){
pred <- predict.regsubsets(regfit.best,
newdata = boston.test,
id = i)
test.errors[i] <- mean((y.test - pred)^2)
}
test.errors
best.size <- which.min(test.errors)
best.size
bestsub.mse <- min(test.errors)
bestsub.mse
coef(regfit.best, best.size)
cat("Best Subset MSE:", bestsub.mse, "\n")## [1] 40.14557 41.87706 42.02985 41.65722 41.84624 40.78108 40.84996 41.13602
## [9] 41.08239 41.08846 41.08306 41.19923
## [1] 1
## [1] 40.14557
## (Intercept) rad
## -2.5242729 0.6702535
## Best Subset MSE: 40.14557Ridge Regression:
x.train <- model.matrix(crim ~ ., boston.train)[,-1]
x.test <- model.matrix(crim ~., boston.test)[,-1]
set.seed(1)
cv.ridge <- cv.glmnet(x.train,
y.train,
alpha = 0)
plot(cv.ridge)
ridge.bestlam <- cv.ridge$lambda.min
cat("Chosen lambda:", ridge.bestlam, "\n")
ridge.pred <- predict(cv.ridge,
s = ridge.bestlam,
newx = x.test)
ridge.mse <- mean((y.test - ridge.pred)^2)
cat("MSE:", ridge.mse, "\n")## Chosen lambda: 0.5919159
## MSE: 40.17065
Lasso Regression:
set.seed(1)
cv.out <- cv.glmnet(x.train,
y.train,
alpha = 1)
plot(cv.out)
#choose lambda
bestlam <- cv.out$lambda.min
#MSE
lasso.pred <- predict(cv.out,
s = bestlam,
newx = x.test)
lasso.mse <- mean((lasso.pred - y.test)^2) #test mse
cat("Lasso MSE:", lasso.mse, "\n")## Lasso MSE: 40.89958
PCR
set.seed(1)
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 6 comps
## CV 9.275 7.681 7.682 7.408 7.144 7.146 7.114
## adjCV 9.275 7.675 7.677 7.405 7.136 7.138 7.106
## 7 comps 8 comps 9 comps 10 comps 11 comps 12 comps
## CV 6.942 6.949 6.896 6.854 6.807 6.728
## adjCV 6.932 6.942 6.892 6.845 6.796 6.717
##
## TRAINING: % variance explained
## 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps 7 comps 8 comps
## X 50.61 63.48 72.71 80.29 86.43 90.37 92.94 95.04
## crim 32.91 33.06 37.89 42.35 42.62 43.12 45.76 45.79
## 9 comps 10 comps 11 comps 12 comps
## X 96.80 98.34 99.50 100.00
## crim 47.16 47.87 48.83 50.21pcr.pred <- predict(pcr.fit,
boston.test,
ncomp = 12)
pcr.mse <- mean((y.test - pcr.pred)^2)
cat("PCR MSE:", pcr.mse, "\n")## PCR MSE: 41.19923b. 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 models demonstrated, the two chosen would be Best subset selection or Ridge Regression. While Best subset had the lowest mean squared error, there was only a slight difference when compared to ridge regression, indicating both models perform well on this data.
–
c. Does your chosen model involve all of the features in the data set? Why or why not?
With Best Subset Selection, not all of the features are included from the data set. The model that was selected based on the lowest MSE score uses Rad as the only predictor. Adding more than this predictor did not prove vital when building the model, since there was no significant improvement to the MSE.
–