# Standard errors for lasso prediction using R

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?

• 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 - Reinstate Monica 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 meta.stats.stackexchange.com/q/2007/159 – 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

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.

• That would explain why the packages don't implement standard errors. – Rob Hyndman Mar 26 '14 at 5:44
• Has there been any progress in estimating lasso and elastic net standard errors and confidence intervals since the time of this question? If not, what are the most common statistical procedures for assessing the accuracy of lasso and elastic net point estimates? – develarist Nov 29 at 18:01

On a related note, which may be helpful, Tibshirani and colleagues have proposed a significance test for the lasso. The paper is available, and titled "A significance test for the lasso". A free version of the paper can be found here

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.

> 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)
[1]  52.772443 -78.690610  24.234753   9.717777 -23.360369 -45.477199
>
> ## sd of y:
> sd_y <- apply(post_pred_y, 1, sd)
[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)))
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


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?

• Bootstrapping can both estimate and correct for bias, though samples need to be fairly large. – Glen_b -Reinstate Monica Apr 7 '14 at 7:20

There is the selectiveInference package in R, https://cran.r-project.org/web/packages/selectiveInference/index.html, that provides confidence intervals and p values for your coefficients fitted by the LASSO, based on the following paper:

Stephen Reid, Jerome Friedman, and Rob Tibshirani (2014). A study of error variance estimation in lasso regression. arXiv:1311.5274

PS: just realise that this produces error estimates for your parameters, not sure for the error on your final prediction, if that's what you're after... I suppose you could use "population prediction intervals" for that if you like (by resampling parameters according to the fit following a multivariate normal distribution).