[SOLVED] STAT340 Discussion 09: the LASSO

STAT340 Discussion 09: the LASSO

STAT340 Discussion 09: the LASSO

Copyright By Assignmentchef assignmentchef

Keith Levin and Wu

link to source

XKCD comic

Refresher: the LASSO

In lecture, we discussed the LASSO, which is a technique for regularizing linear regression while simultaneously doing variable selection.

Recall that under the LASSO, we try to choose our coefficients (beta_0,beta_1,beta_2,dots,beta_p) to minimize [
sum_{i=1}^n left( Y_i beta^T X_i right)^2 + lambda sum_{j=1}^p |beta_j|,
] where (lambda ge 0) controls the amount of regularization (i.e., how much we care about the second term) and we are assuming that the vectors of predictors (X_1,X_2,dots,X_n in mathbb{R}^{p+1}) include an entry equal to 1 to account for the intercept term: [
X_i = left( 1, X_{i,1}, X_{i,2}, dots, X_{i,p} right)^T in mathbb{R}^{p+1}.

The LASSO objective (an objective is just a quantity that we want to optimize in this case, we are trying to minimize the LASSO cost) is very similar to ridge regression, which we discussed at length in lecture, and tries to minimize: [
operatorname{RSS} + lambda sum_{j=1}^p beta_j^2.

The difference is that ridge regression penalizes the sum of squares of the coefficients, while the LASSO penalizes the sum of their absolute values. It turns out that this small change has a surprisingly large effect on the solution that we find.

Specifically, as discussed in lecture, the LASSO penalty actually encourages us to set unhelpful coefficients to zero. This is in contrast to ridge regression, which encourages coefficients that do not help much with prediction to be close to zero, but if you examine our worked example in lecture, youll see that the coefficients in our estimated models, even for large values of (lambda), are not equal to zero.

In lecture, we waved our hands at the LASSO, and we mentioned that glmnet is a popular package for fitting the LASSO, but we didnt dive into implementation like we did with ridge regression. So lets get a bit of practice with glmnet now.

mtcars yet again

Weve been using the mtcars data set as a running example in lecture for our regression problems, so why stop now? Lets load the mtcars data set and just remind ourselves of the variables involved.

data(mtcars);
head(mtcars);

##mpg cyl disphp dratwtqsec vs am gear carb
## Mazda RX4 21.0 6160 110 3.90 2.620 16.460144
## Mazda RX4 Wag 21.0 6160 110 3.90 2.875 17.020144
## Datsun 71022.8 410893 3.85 2.320 18.611141
## Hornet 4 Drive21.4 6258 110 3.08 3.215 19.441031
## Hornet Sportabout 18.7 8360 175 3.15 3.440 17.020032
## Valiant 18.1 6225 105 2.76 3.460 20.221031

In the event that this output doesnt fit on one screen, heres the full list of variables:

names(mtcars);

##[1] mpgcyldisp hp drat wt qsec vs am gear
## [11] carb

Our goal is still to predict the fuel efficiency (mpg) from all of the other variables, this time using the LASSO.

Installation and overview of glmnet

To start, we need to install glmnet. So lets do that.

# Uncomment this line and run it if you need to install glmnet.
#install.packages(glmnet);

If you run into weird errors about clang or anything like that, you might need to try including the dependencies=TRUE argument in install.packages. Also, if your version of R is too old (something like older than 3.4, if I remember correctly), you may need to update R.

Once you have glmnet installed, lets have a look at the documentation.

library(glmnet);

## Warning: package glmnet was built under R version 3.6.2

## Loading required package: Matrix

## Loaded glmnet 4.1-3

Read over the help documentation. Note that the documentation refers to generalized linear model via penalized maximum likelihood. This is a fancy way of saying that glmnet covers a much wider array of regression problems than just linear regression with the LASSO. Well return to that below.

Skipping to the Arguments section of the documentation, we see that we need to first specify the input data matrix x (i.e., the matrix of predictors) and the responses y.

In our data, then the y argument will be the column of car mpg values, and the matrix x will be everything else.

We also need to specify a family. This tells glmnet what kind of regression we want to do (e.g., linear vs logistic).

Skipping down a few arguments, we see that there are some arguments related to specifying the (lambda) parameter that controls the amount of regularization. Well just set (lambda) by hand, using the lambda parameter, but you can see that there are other options in the documentation.

Reading through more of the docs, youll see that are plenty more arguments available to modify the behavior of our fitting procedure. Well leave most of those for another day when you know more about the nuts and bolts of regression and optimization.

However, there is one very important one: if you read the documentation for the alpha parameter, this is described as The elasticnet mixing parameter, with (0 le alpha le 1). We need to set (alpha = 1) to get the LASSO.

Fitting a model

So lets get down to it and fit a model, already!

From reading the documentation, we see that we need to split the mtcars data set into a column of responses and everything else. Recall that the mpg column is the 1st column,

head(mtcars)

##mpg cyl disphp dratwtqsec vs am gear carb
## Mazda RX4 21.0 6160 110 3.90 2.620 16.460144
## Mazda RX4 Wag 21.0 6160 110 3.90 2.875 17.020144
## Datsun 71022.8 410893 3.85 2.320 18.611141
## Hornet 4 Drive21.4 6258 110 3.08 3.215 19.441031
## Hornet Sportabout 18.7 8360 175 3.15 3.440 17.020032
## Valiant 18.1 6225 105 2.76 3.460 20.221031

So lets split it out.

y_mtc <- mtcars[,1]; # Grab just the first column# … and the predictors are everything else.x_mtc <- mtcars[, -c(1)];# Just to verify:head(x_mtc);## cyl disphp dratwtqsec vs am gear carb## Mazda RX4 6160 110 3.90 2.620 16.460144## Mazda RX4 Wag 6160 110 3.90 2.875 17.020144## Datsun 710410893 3.85 2.320 18.611141## Hornet 4 Drive6258 110 3.08 3.215 19.441031## Hornet Sportabout 8360 175 3.15 3.440 17.020032## Valiant 6225 105 2.76 3.460 20.221031Lets start with a sanity check. Well fit the LASSO with (lambda=0). This is just linear regression, so the predictions should be the same (up to rounding errors) if we use the built-in linear regression in lm.# Remember, alpha=1 for the LASSO.mtc_lasso_lambda0 <- glmnet(x_mtc, y_mtc, alpha = 1, lambda=0);Note: if you get an error along the lines of missing functions in Rcpp, try running install.packages(‘Rcpp’); library(Rcpp). This will reinstall the Rcpp package and should correct the issue. This was a common problem for students earlier in the semester, so you are unlikely to encounter it now, but just in case.Lets extract the coefficients just to see what theyre like.coef( mtc_lasso_lambda0 )## 11 x 1 sparse Matrix of class “dgCMatrix”##s0## (Intercept) 12.19850081## cyl -0.09882217## disp 0.01307841## hp-0.02142912## drat 0.79812453## wt-3.68926778## qsec 0.81769993## vs 0.32109677## am 2.51824708## gear 0.66755681## carb-0.21040602Now, lets fit linear regression and compare the coefficients.#TODO: code goes here.# Fit linear regression to the mtcars data set and extract the coefficients.mtc_vanilla_lm <- lm( mpg ~ ., mtcars )coef(mtc_vanilla_lm)## (Intercept) cyldisphpdratwt ## 12.30337416 -0.111440480.01333524 -0.021482120.78711097 -3.71530393 ##qsecvsamgearcarb ##0.821040750.317762812.520226890.65541302 -0.19941925The coefficients are largely the same, though they are often only in agreement to one or two decimal points. This has a lot to do with the fact that the scalings of the variables in the mtcars data set are very different. Lets compare the model predictions. Note that this is still just prediction on the train set. We are doing this not to evaluate model fit but just to verify that the two different models we fit are reasonably similar (as they should be, because they are both actually the same linear regression model with no regularization).# Note that the lasso model object requires its input prediction examples# to be in a matrix, so we are obliging it with as.matrix().predict( mtc_lasso_lambda0, newx=as.matrix(x_mtc) )## s0## Mazda RX4 22.59390## Mazda RX4 Wag 22.11105## Datsun 71026.25301## Hornet 4 Drive21.22927## Hornet Sportabout 17.68818## Valiant 20.38336## Duster 36014.37072## Merc 240D 22.49403## Merc 23024.41218## Merc 28018.71290## Merc 280C 19.20352## Merc 450SE14.19206## Merc 450SL15.60995## Merc 450SLC 15.75257## Cadillac Fleetwood12.02110## 10.93297## Chrysler Imperial 10.49791## Fiat 12827.77790## Honda Civic 29.88808## Toyota Corolla29.50987## Toyota Corona 23.62986## Dodge Challenger16.94555## AMC Javelin 17.73892##2813.30725## Pontiac Firebird16.68582## Fiat X1-9 28.29339## Porsche 914-2 26.15796## Lotus Europa27.62799## Ford Pantera L18.88532## Ferrari Dino19.68972## Maserati Bora 13.93165## Volvo 142E24.37203And compare withpredict( mtc_vanilla_lm, mtcars)## Mazda RX4 Mazda RX4 WagDatsun 710Hornet 4 Drive ##22.5995122.1118926.2506421.23740 ## Hornet Sportabout 360 Merc 240D ##17.6934320.3830414.3862622.49601 ##Merc 230Merc 280 Merc 280CMerc 450SE ##24.4190918.6990319.1916514.17216 ##Merc 450SL Merc 450SLCCadillac Fleetwood ##15.5995715.7422212.0340110.93644 ## Chrysler ImperialFiat 128 Honda CivicToyota Corolla ##10.4936327.7729129.8967429.51237 ## Toyota CoronaDodge Challenger AMC Javelin 28 ##23.6431016.9430517.7321813.30602 ##Pontiac Firebird Fiat X1-9 Porsche 914-2Lotus Europa ##16.6916828.2934726.1529527.63627 ##Ford Pantera LFerrari 142E ##18.8700419.6938313.9411224.36827Okay, its a little annoying to scroll back and forth, but youll see that the predictions agree quite well.Adding some regularizationAll right, now, modify the above code to run LASSO with (lambda=1) and extract the coefficients. Compare the coefficients with those returned by linear regression. Any observations? You may find it useful to plot the different sets of coefficients in a plot.#TODO: code goes here.TODO: discussion/observations go here.Sending coefficients to zeroOkay, now try running the same LASSO code with increasing values of (lambda). What happens as you increase (lambda)?Try a few different choices of (lambda) and compare the resulting coefficients. Note that you can pass a vector of values to the lambda argument in glmnet to use multiple values for (lambda) at once.Question 1: What is (approximately) the smallest value of (lambda) for which all of the coefficients get set to zero?#TODO: code goes here.TODO: discussion/observations go here.Question 2: When (lambda = 0), we saw that all of the coefficients were non-zero, and we know that as (lambda) increases, the LASSO penalizes non-zero coefficients more and more. So there should exist a point, as (lambda) increases from (0), at which the first coefficient gets set to zero. What is (approximately) this value of (lambda)?#TODO: code goes here.TODO: discussion/observations go here. CS: assignmentchef QQ: 1823890830 Email: [email protected]