Cross validating lasso regression in R

Question

The R function cv.glm (library: boot) calculates the estimated K-fold cross-validation prediction error for generalized linear models and returns delta. Does it make sense to use this function for a lasso regression (library: glmnet) and if so, how can it be carried out? The glmnet library uses a cross-validation to get the best turning parameter, but I did not find any example that cross-validates the final glmnet equation.

Answer 1

An example on how to do vanilla plain cross-validation for lasso in glmnet on mtcars data set.

1. Load data set.

2. Prepare features (independent variables). They should be of matrix class. The easiest way to convert df containing categorical variables into matrix is via model.matrix. Mind you, by default glmnet fits intercept, so you'd better strip intercept from model matrix.

3. Prepare response (dependent variable). Let's code cars with above average mpg as efficient ('1') and the rest as inefficient ('0'). Convert this variable to factor.

4. Run cross-validation via cv.glmnet. It will pickup alpha=1 from default glmnet parameters, which is what you asked for: lasso regression.

5. By examining the output of cross-validation you may be interested in at least 2 pieces of information:

   • lambda, that minimizes cross-validated error. glmnet actually provides 2 lambdas: lambda.min and lambda.1se. It's your judgement call as a practicing statistician which to use.

   • resulting regularized coefficients.

Please see the R code per the above instructions:

# Load data set
data("mtcars")

# Prepare data set 
x   <- model.matrix(~.-1, data= mtcars[,-1])
mpg <- ifelse( mtcars$mpg < mean(mtcars$mpg), 0, 1)
y   <- factor(mpg, labels = c('notEfficient', 'efficient'))

library(glmnet)

# Run cross-validation
mod_cv <- cv.glmnet(x=x, y=y, family='binomial')

mod_cv$lambda.1se
[1] 0.108442

coef(mod_cv, mod_cv$lambda.1se)
                     1
(Intercept)  5.6971598
cyl         -0.9822704
disp         .        
hp           .        
drat         .        
wt           .        
qsec         .        
vs           .        
am           .        
gear         .        
carb         .  

mod_cv$lambda.min
[1] 0.01537137

coef(mod_cv, mod_cv$lambda.min)
                      1
(Intercept)  6.04249733
cyl         -0.95867199
disp         .         
hp          -0.01962924
drat         0.83578090
wt           .         
qsec         .         
vs           .         
am           2.65798203
gear         .         
carb        -0.67974620

Final comments:

• note, the model's output says nothing about statistical significance of the coefficients, only values.

• l1 penalizer (lasso), which you asked for, is notorious for instability as evidenced in this blog post and this stackexchange question. A better way could be to cross-validate on alpha too, which would let you decide on proper mix of l1 and l2 penalizers.

• an alternative way to do cross-validation could be to turn to caret's train( ... method='glmnet')

• and finally, the best way to learn more about cv.glmnet and it's defaults coming from glmnet is of course ?glmnet in R's console )))