Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

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 stackoverflow.com Mar 26 at 16:59

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

5  
This really should be migrated back to stats.stackexchange since it's actually a statistics question –  hadley Mar 26 at 12:52
1  
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
2  
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

6 Answers 6

up vote 24 down vote accepted

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.

share|improve this answer
2  
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".

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...

share|improve this answer
1  
As it stands this answer is a little brief. Could you expand a little further about how the Bayesian Lasso works and perhaps give a reference which discusses it? –  Glen_b Apr 7 at 7:22
    
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. –  Sandipan Karmakar Apr 8 at 7:44
    
Then one can use Gibbs Sampling for simulating the chain. See <stat.ufl.edu/archived/casella/Papers/Lasso.pdf>;. –  Sandipan Karmakar Apr 8 at 7:45
    
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}$... and 3. The hierarchical structure of the problem at hand cannot be encoded using frequentist model, which is quite easy in Bayesian framework. If you need any help write to me at my mail id... –  Sandipan Karmakar Apr 8 at 7:46
    
This is the sort of information that is useful (thank you), but please note that your actual answer is still as brief as before. Could you edit the information into your answer please, rather than discuss it in comments? The aim is to have a collection of useful, reasonably complete answers. I previously upvoted your answer in anticipation of such an edit. –  Glen_b Apr 8 at 8:47

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
1  
Bootstrapping can both estimate and correct for bias, though samples need to be fairly large. –  Glen_b Apr 7 at 7:20

Bayesian LASSO needs prior distributions to be assigned to the parameters of the model. In LASSO model, we have the objective function $||y−X\beta||_2^2+\lambda||\beta||_1$ with $\lambda$ as the regularization parameter. Here as we have $\ell_1$-norm for $\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 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 $\beta$ and

  3. The hierarchical structure of the problem at hand cannot be encoded using frequentist model, which is quite easy in Bayesian framework. If you need any help write to me at my mail id

share|improve this answer
    
I upvoted both your answers, so you can edit this answer into your first one which I think would be better for the Wiki feeling of the site. –  Momo Apr 9 at 11:05

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

Your Answer

 
discard

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.