Running head: HOMEWORK 8 – STEPWISE REGRESSION, LASSO AND ELASTIC NET 1
Homework 8 – Use Stepwise Regression, Lasso, Elastic net and glmnet
Amitava Chatterjee
OMS Analytics GATECH – Fall 2019HOMEWORK 8 – STEPWISE REGR
...
Running head: HOMEWORK 8 – STEPWISE REGRESSION, LASSO AND ELASTIC NET 1
Homework 8 – Use Stepwise Regression, Lasso, Elastic net and glmnet
Amitava Chatterjee
OMS Analytics GATECH – Fall 2019HOMEWORK 8 – STEPWISE REGRESSION, LASSO AND ELASTIC NET 2
Abstract
Use different Regression models to look at the crime data with scaling to align and assess alpha
across stepwise, lasso, and elastic net regressions.HOMEWORK 8 – STEPWISE REGRESSION, LASSO AND ELASTIC NET 3
Homework 8 – Use Stepwise Regression, Lasso, Elastic net and glmnet
Stepwise Regression
As per Wikipedia - Stepwise regression is a method of fitting regression models in which the choice
of predictive variables is carried out by an automatic procedure. In each step, a variable is considered
for addition to or subtraction from the set of explanatory variables based on some prespecified
criterion.
Lasso
As per Wikipedia - Lasso regression is a type of linear regression that uses shrinkage. Shrinkage is
where data values are shrunk towards a central point, like the mean. The lasso procedure encourages
simple, sparse models (i.e. models with fewer parameters).
Elastic net
The elastic net is a regularized regression method that linearly combines the L1 and L2 penalties of
the lasso and ridge methods.HOMEWORK 8 – STEPWISE REGRESSION, LASSO AND ELASTIC NET 4
Question 11.1
Using the crime data set uscrime.txt from Questions 8.2, 9.1, and 10.1, build a regression
model using:
1. Stepwise regression
2. Lasso
3. Elastic net
For Parts 2 and 3, remember to scale the data first – otherwise, the regression coefficients will
be on different scales and the constraint won’t have the desired effect.
For Parts 2 and 3, use the glmnet function in R.
Notes on R:
• For the elastic net model, what we called λ in the videos, glmnet calls “alpha”; you can get a
range of results by varying alpha from 1 (lasso) to 0 (ridge regression) [and, of course, other values of
alpha in between].
• In a function call like glmnet(x,y,family=”mgaussian”,alpha=1) the predictors x need to be
in R’s matrix format, rather than data frame format. You can convert a data frame to a matrix using
as.matrix – for example, x <- as.matrix(data[,1:n-1])
• Rather than specifying a value of T, glmnet returns models for a variety of values of T.
Answer 11.1
Solution runs stepwise regression, lasso, and elastic net on both the scaled raw data and principal
components found using PCA. The R code for each model covers three things: (1) uses the method to
identify a set of variables to use, (2) builds a regression model using those variables, and (3)
eliminates the insignificant variables in the regression and then builds a regression using theHOMEWORK 8 – STEPWISE REGRESSION, LASSO AND ELASTIC NET 5
remaining variables. After building each model, the code reports the R-squared value on the training
data, and then uses cross-validation to estimate the real R-squared value of the model.
The elastic net models, are used to test 11 different values of alpha, from 0.0 to 1.0 at intervals of 0.1.
One would see slightly different results depending on the random number generator.
# Clear environment
rm(list = ls())
library(caret)
## Loading required package: lattice
## Loading required package: ggplot2
library(glmnet)
## Loading required package: Matrix
## Loading required package: foreach
## Loaded glmnet 2.0-18
# Set the random number generator seed - the results are reproducible
set.seed(1)
# Read in the data
#
ac_data = read.table("uscrime.txt", stringsAsFactors = FALSE, header =
TRUE)
#
# Make sure the data is read correctly
#
head(ac_data)
## M So Ed Po1 Po2 LF M.F Pop NW U1 U2 Wealth Ineq
## 1 15.1 1 9.1 5.8 5.6 0.510 95.0 33 30.1 0.108 4.1 3940 26.1
## 2 14.3 0 11.3 10.3 9.5 0.583 101.2 13 10.2 0.096 3.6 5570 19.4
## 3 14.2 1 8.9 4.5 4.4 0.533 96.9 18 21.9 0.094 3.3 3180 25.0
## 4 13.6 0 12.1 14.9 14.1 0.577 99.4 157 8.0 0.102 3.9 6730 16.7
## 5 14.1 0 12.1 10.9 10.1 0.591 98.5 18 3.0 0.091 2.0 5780 17.4
## 6 12.1 0 11.0 11.8 11.5 0.547 96.4 25 4.4 0.084 2.9 6890 12.6
## Prob Time Crime
## 1 0.084602 26.2011 791
## 2 0.029599 25.2999 1635HOMEWORK 8 – STEPWISE REGRESSION, LASSO AND ELASTIC NET 6
## 3 0.083401 24.3006 578
## 4 0.015801 29.9012 1969
## 5 0.041399 21.2998 1234
## 6 0.034201 20.9995 682
tail(ac_data)
## M So Ed Po1 Po2 LF M.F Pop NW U1 U2 Wealth Ineq
## 42 14.1 0 10.9 5.6 5.4 0.523 96.8 4 0.2 0.107 3.7 4890 17.0
## 43 16.2 1 9.9 7.5 7.0 0.522 99.6 40 20.8 0.073 2.7 4960 22.4
## 44 13.6 0 12.1 9.5 9.6 0.574 101.2 29 3.6 0.111 3.7 6220 16.2
## 45 13.9 1 8.8 4.6 4.1 0.480 96.8 19 4.9 0.135 5.3 4570 24.9
## 46 12.6 0 10.4 10.6 9.7 0.599 98.9 40 2.4 0.078 2.5 5930 17.1
## 47 13.0 0 12.1 9.0 9.1 0.623 104.9 3 2.2 0.113 4.0 5880 16.0
## Prob Time Crime
## 42 0.088904 12.1996 542
## 43 0.054902 31.9989 823
## 44 0.028100 30.0001 1030
## 45 0.056202 32.5996 455
## 46 0.046598 16.6999 508
## 47 0.052802 16.0997 849
# Crime is response, other variables are predictors
—————————- Stepwise Regression ————————————-
Stepwise Regression using original variables and Cross Validation In backward stepwise regression.
Our lower model will have only the intercept and all variables in our full model.
#Scaling the data except the response variable and categorical
scaledData =
as.data.frame(scale(ac_data[,c(1,3,4,5,6,7,8,9,10,11,12,13,14,15)]))
scaledData = cbind(ac_data[,2],scaledData,ac_data[,16]) # Add column 2
back in
colnames(scaledData)[1] = "So"
colnames(scaledData)[16] = "Crime"
# Perform 5 fold CV using the code below
ctrl = trainControl(method = "repeatedcv", number = 5, repeats = 5)
lmFit_Step = train(Crime ~ ., data = scaledData, "lmStepAIC", scope =
list(lower = Crime~1, upper = Crime~.), direction =
"backward",trControl=ctrl, trace =FALSE)
summary(lmFit_Step)HOMEWORK 8 – STEPWISE REGRESSION, LASSO AND ELASTIC NET 7
##
## Call:
## lm(formula = .outcome ~ M + Ed + Po1 + M.F + U1 + U2 + Ineq +
## Prob, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -444.70 -111.07 3.03 122.15 483.30
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 905.09 28.52 31.731 < 2e-16 ***
## M 117.28 42.10 2.786 0.00828 **
## Ed 201.50 59.02 3.414 0.00153 **
## Po1 305.07 46.14 6.613 8.26e-08 ***
## M.F 65.83 40.08 1.642 0.10874
## U1 -109.73 60.20 -1.823 0.07622 .
## U2 158.22 61.22 2.585 0.01371 *
## Ineq 244.70 55.69 4.394 8.63e-05 ***
## Prob -86.31 33.89 -2.547 0.01505 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 195.5 on 38 degrees of freedom
## Multiple R-squared: 0.7888, Adjusted R-squared: 0.7444
## F-statistic: 17.74 on 8 and 38 DF, p-value: 1.159e-10
Fit a new model with 8 variables
mod_Step = lm(Crime ~ M.F+U1+Prob+U2+M+Ed+Ineq+Po1, data = scaledData)
summary(mod_Step)
##
## Call:
## lm(formula = Crime ~ M.F + U1 + Prob + U2 + M + Ed + Ineq + Po1,
## data = scaledData)
##
## Residuals:
## Min 1Q Median 3Q Max
## -444.70 -111.07 3.03 122.15 483.30
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 905.09 28.52 31.731 < 2e-16 ***
## M.F 65.83 40.08 1.642 0.10874
## U1 -109.73 60.20 -1.823 0.07622 .
## Prob -86.31 33.89 -2.547 0.01505 *
## U2 158.22 61.22 2.585 0.01371 *
## M 117.28 42.10 2.786 0.00828 **
## Ed 201.50 59.02 3.414 0.00153 **
## Ineq 244.70 55.69 4.394 8.63e-05 ***HOMEWORK 8 – STEPWISE REGRESSION, LASSO AND ELASTIC NET 8
## Po1 305.07 46.14 6.613 8.26e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 195.5 on 38 degrees of freedom
## Multiple R-squared: 0.7888, Adjusted R-squared: 0.7444
## F-statistic: 17.74 on 8 and 38 DF, p-value: 1.159e-10
We obtain an Adjusted R-SQuared value = 0.7444 using the selected 8 variables using Backward
StepWise regression and Cross Validation
Now let’s use cross-validation to see how good this model really is. Because we only have 47 data
points, let’s use 47-fold cross-validation (equivalently, leave-one-out cross-validation).
SStot = sum((ac_data$Crime - mean(ac_data$Crime))^2)
totsse = 0
for(i in 1:nrow(scaledData)) {
mod_Step_i = lm(Crime ~ M.F+U1+Prob+U2+M+Ed+Ineq+Po1, data =
scaledData[-i,])
pred_i = predict(mod_Step_i,newdata=scaledData[i,])
totsse = totsse + ((pred_i - ac_data[i,16])^2)
}
r2_md = 1 - totsse/SStot
r2_md
## 1
## 0.667621
In the model above, the p-value for M.F is above 0.1.
We might keep it in the model, because it’s close to 0.1 and might be important. That’s what we tested
above.
Or, we might remove it, and re-run the model without it. Let’s see what happens if we do:
mod_Step = lm(Crime ~ U1+Prob+U2+M+Ed+Ineq+Po1, data = scaledData)
summary(mod_Step)
##
## Call:
## lm(formula = Crime ~ U1 + Prob + U2 + M + Ed + Ineq + Po1, data =
scaledData)
##
## Residuals:
## Min 1Q Median 3Q Max
## -520.76 -105.67 9.53 136.28 519.37
##
## Coefficients:HOMEWORK 8 – STEPWISE REGRESSION, LASSO AND ELASTIC NET 9
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 905.09 29.14 31.062 < 2e-16 ***
## U1 -63.86 54.48 -1.172 0.2482
## Prob -84.83 34.61 -2.451 0.0188 *
## U2 134.13 60.71 2.209 0.0331 *
## M 134.20 41.70 3.218 0.0026 **
## Ed 244.38 54.07 4.520 5.62e-05 ***
## Ineq 264.65 55.52 4.767 2.61e-05 ***
## Po1 314.89 46.73 6.738 4.91e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 199.8 on 39 degrees of freedom
## Multiple R-squared: 0.7738, Adjusted R-squared: 0.7332
## F-statistic: 19.06 on 7 and 39 DF, p-value: 8.805e-11
U1 doesn’t look significant… we can take it out too, and re-run the model.
mod_Step = lm(Crime ~ Prob+U2+M+Ed+Ineq+Po1, data = scaledData)
summary(mod_Step)
##
## Call:
## lm(formula = Crime ~ Prob + U2 + M + Ed + Ineq + Po1, data =
scaledData)
##
## Residuals:
## Min 1Q Median 3Q Max
## -470.68 -78.41 -19.68 133.12 556.23
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 905.09 29.27 30.918 < 2e-16 ***
## Prob -86.44 34.74 -2.488 0.01711 *
## U2 75.47 34.55 2.185 0.03483 *
## M 131.98 41.85 3.154 0.00305 **
## Ed 219.79 50.07 4.390 8.07e-05 ***
## Ineq 269.91 55.60 4.855 1.88e-05 ***
## Po1 341.84 40.87 8.363 2.56e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 200.7 on 40 degrees of freedom
## Multiple R-squared: 0.7659, Adjusted R-squared: 0.7307
## F-statistic: 21.81 on 6 and 40 DF, p-value: 3.418e-11
This model looks good, so now let’s see how it cross-validates:
SStot = sum((ac_data$Crime - mean(ac_data$Crime))^2)
totsse = 0
for(i in 1:nrow(scaledData)) {
mod_Step_i = lm(Crime ~ Prob+U2+M+Ed+Ineq+Po1, data = scaledData[-i,])HOMEWORK 8 – STEPWISE REGRESSION, LASSO AND ELASTIC NET 10
pred_i = predict(mod_Step_i,newdata=scaledData[i,])
totsse = totsse + ((pred_i - ac_data[i,16])^2)
}
r2_md = 1 - totsse/SStot
r2_md
## 1
## 0.6661638
So, cross-validation shows that it’s about the same whether we include M.F and U1 (0.668) or not
(0.666). That gives some support to the idea that M.F and U1 really aren’t significant. Since the
quality is about the same, we should probably use the simpler model.
—————————- Lasso Regression ————————————-
#building lasso
XP=data.matrix(scaledData[,-16])
YP=data.matrix(scaledData$Crime)
lasso=cv.glmnet(x=as.matrix(scaledData[,-
16]),y=as.matrix(scaledData$Crime),alpha=1,
nfolds = 5,type.measure="mse",family="gaussian")
#Output the coefficients of the variables selected by lasso
coef(lasso, s=lasso$lambda.min)
## 16 x 1 sparse Matrix of class "dgCMatrix"
## 1
## (Intercept) 889.059998
## So 47.073756
## M 85.028207
## Ed 124.299584
## Po1 308.896175
## Po2 .
## LF 1.112481
## M.F 52.034368
[Show More]
