# 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?

-

## migrated from stackoverflow.comMar 26 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 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 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 at 23:17
You can do this easily in a Bayesian framework using package monomvn, see my answer below. –  fabians Apr 9 at 15:37

## 5 Answers

This paper suggests 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 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".

-

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

-

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 http://www.stat.ufl.edu/archived/casella/Papers/Lasso.pdf.

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.

-

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 Apr 7 at 7:20