Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I'm trying to use a LASSO model for prediction, and I need to estimate standard errors. Surely someone has already written a package to do this. But as far as I can see, none of the packages on CRAN that do predictions using a LASSO will return standard errors for those predictions.

So my question is: Is there a package or some R code available to compute standard errors for LASSO predictions?

share|improve this question

migrated from Mar 26 '14 at 16:59

This question came from our site for professional and enthusiast programmers.

This really should be migrated back to stats.stackexchange since it's actually a statistics question – hadley Mar 26 '14 at 12:52
To clarify the underlying nature of this question (since it has been bouncing back & forth b/t CV & SO), I wonder if we might edit the title, Rob. How about 'Why doesn't there seem to be a package for LASSO standard errors, are they difficult to compute?', or something like that, perhaps coupled w/ some minor edits to the body to make it consistent. I think that would make it more clearly on-topic on CV so that this ambiguity does not arise & we don't have to go back & forth. – gung Mar 26 '14 at 14:27
I could make the question more about the statistical methodology, but that wasn't actually what I wanted to know. There should be a place for questions on CV about what software implements a given method. Further discussion at – Rob Hyndman Mar 26 '14 at 23:17
You can do this easily in a Bayesian framework using package monomvn, see my answer below. – fabians Apr 9 '14 at 15:37
up vote 36 down vote accepted

Kyung et al. (2010), "Penalized regression, standard errors, & Bayesian lassos", Bayesian Analysis , 5, 2, suggest that there might not be a consensus on a statistically valid method of calculating standard errors for the lasso predictions. Tibshirani seems to agree (slide 43) that standard errors are still an unresolved issue.

share|improve this answer
That would explain why the packages don't implement standard errors. – Rob Hyndman Mar 26 '14 at 5:44

On a related note, which may be helpful, Tibshirani and colleagues have proposed a significance test for the lasso. The (not yet published) manuscript is available, and titled "A significance test for the lasso".

† Now published: Lockhart(2014), "A significance test for the lasso", Ann. Statist., 42, 2

share|improve this answer

Sandipan Karmakar answer tells you what to do, this should help you on the "how":

> library(monomvn)
> ## following the lars diabetes example
> data(diabetes)
> str(diabetes)
'data.frame':   442 obs. of  3 variables:
 $ x : AsIs [1:442, 1:10] 0.038075.... -0.00188.... 0.085298.... -0.08906.... 0.005383.... ...
      ..- attr(*, "dimnames")=List of 2
      .. ..$ : NULL
  .. ..$ : chr  "age" "sex" "bmi" "map" ...

 $ y : num  151 75 141 206 135 97 138 63 110 310 ...


> ## Bayesian Lasso regression
> reg_blas <- with(diabetes, blasso(x, y))
t=100, m=8
t=200, m=5
t=300, m=8
t=400, m=8
t=500, m=7
t=600, m=8
t=700, m=8
t=800, m=8
t=900, m=5
> ## posterior mean beta (setting those with >50% mass at zero to exactly zero)
> (beta <- colMeans(reg_blas$beta) * (colMeans(reg_blas$beta != 0)  > 0.5))
      b.1       b.2       b.3       b.4       b.5       b.6       b.7       b.8 
   0.0000 -195.9795  532.7136  309.1673 -101.1288    0.0000 -196.4315    0.0000 
      b.9      b.10 
 505.4726    0.0000 
> ## n x nsims matrix of realizations from the posterior predictive:
> post_pred_y <- with(reg_blas, X %*% t(beta))
> ## predictions:
> y_pred <- rowMeans(post_pred_y)
> head(y_pred)
[1]  52.772443 -78.690610  24.234753   9.717777 -23.360369 -45.477199
> ## sd of y:
> sd_y <- apply(post_pred_y, 1, sd)
> head(sd_y)
[1] 6.331673 6.756569 6.031290 5.236101 5.657265 6.150473
> ## 90% credible intervals
> ci_y <- t(apply(post_pred_y, 1, quantile, probs=c(0.05, 0.95)))
> head(ci_y)
             5%       95%
[1,]  42.842535  62.56743
[2,] -88.877760 -68.47159
[3,]  14.933617  33.85679
[4,]   1.297094  18.01523
[5,] -32.709132 -14.13260
[6,] -55.533807 -35.77809
share|improve this answer

Bayesian LASSO is the only alternative to the problem of calculating standard errors. Standard errors are automatically calculated in Bayesian LASSO...You can implement Bayesian LASSO very easily using Gibbs Sampling scheme...

Bayesian LASSO needs prior distributions to be assigned to the parameters of the model. In LASSO model, we have the objective function $||\mathbf{y}-\mathbf{X}\boldsymbol{\beta}||_2^2 + \lambda||\boldsymbol{\beta}||_1$ with $\lambda$ as the regularization parameter. Here as we have $\ell_1$-norm for $\boldsymbol{\beta}$ so, a special type of prior distribution is needed for this, LAPLACE distribution a scale mixture of normal distribution with exponential distribution as mixing density. Based upon the full conditional posteriors of each of the parameters are to be deduced.

Then one can use Gibbs Sampling for simulating the chain. See Park & Cassella (2008), "The Bayesian Lasso", JASA, 103, 482.

There are three inherent drawbacks of frequentist LASSO:

  1. One has to choose $\lambda$ by cross validation or other means.

  2. Standard errors are difficult to calculate as the LARS and other algorithms produce point estimates for $\boldsymbol{\beta}$.

  3. The hierarchical structure of the problem at hand cannot be encoded using frequentist model, which is quite easy in Bayesian framework.

share|improve this answer

To add to the answers above, the issue appears to be that even a bootstrap is likely insufficient as the estimate from the penalized model is biased and bootstrapping will only speak to the variance - ignoring the bias of the estimate. This is nicely summarized in the vignette for the penalized package on Page 18.

If being used for prediction however, why is a standard error from the model required? Can you not cross validate or bootstrap appropriately and produce a standard error around a metric related to prediction such as MSE?

share|improve this answer
Bootstrapping can both estimate and correct for bias, though samples need to be fairly large. – Glen_b Apr 7 '14 at 7:20

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.