38
$\begingroup$

I want to use Lasso or ridge regression for a model with more than 50,000 variables. I want do so using software package in R. How can I estimate the shrinkage parameter ($\lambda$)?

Edits:

Here is the point I got up to:

set.seed (123)
Y <- runif (1000)
Xv <- sample(c(1,0), size= 1000*1000,  replace = T)
X <- matrix(Xv, nrow = 1000, ncol = 1000)

mydf <- data.frame(Y, X)

require(MASS)
lm.ridge(Y ~ ., mydf)

plot(lm.ridge(Y ~ ., mydf,
              lambda = seq(0,0.1,0.001)))

enter image description here

My question is: How do I know which $\lambda$ is best for my model?

Improve this question
$\endgroup$
7
  • 3
    $\begingroup$ The smoothing parameter is not statistically estimable but is using chosen to maximize out-of-sample fit using, for example, cross validation. I think the standard packages for LASSO and Ridge Regression in R have built in functionality to do this for you - have you looked into that? $\endgroup$
    – Macro
    Apr 16, 2012 at 12:22
  • 4
    $\begingroup$ I disagree - you can estimate the smoothing parameter, using a mixed model approach. Reml methods exist as do heirarchical bayes methods. You do not need expensive cross validation. $\endgroup$
    – probabilityislogic
    Apr 16, 2012 at 12:35
  • 1
    $\begingroup$ @probabilityislogic thank you for the information. It would be great if there is few detail on script how we can do this using reml $\endgroup$
    – John
    Apr 16, 2012 at 12:44
  • 3
    $\begingroup$ For ridge regression leave-one-out cross-validation is essentially free (Allen's PRESS statistic) and I have found it to be a reasonably good method. However with so may features pretty much any method you use will be unstable and a fully Bayesian approach marginalising over both the parameters and the regularisation parameters is likely to be a more reliable solution (as I think probabilityislogic was suggesting). If you don't like Bayesian methods, then use bagging and re-estimate lambda each time. $\endgroup$
    – Dikran Marsupial
    Apr 16, 2012 at 19:27
  • 2
    $\begingroup$ @Macro - (nothing like a response 18 months later). In mixed model approach there are two additional terms which depend only on $\lambda$ but not $\beta$. These are $-k\log(\lambda)$ and $\log|X^TX+\lambda I|$ where $k$ is the number of betas and X the predictor matrix. The first term comes from $\beta\sim N(0,\sigma^2\lambda^{-1})$ where $\sigma^2$ is error variance. The second term is a REML-correction to account for the uncertainty in plugging in $\beta=\hat{\beta}$. $\endgroup$
    – probabilityislogic
    Sep 8, 2013 at 12:36

3 Answers 3

Reset to default
25
$\begingroup$

The function cv.glmnet from the R package glmnet does automatic cross-validation on a grid of $\lambda$ values used for $\ell_1$-penalized regression problems. In particular, for the lasso. The glmnet package also supports the more general elastic net penalty, which is a combination of $\ell_1$ and $\ell_2$ penalization. As of version 1.7.3. of the package taking the $\alpha$ parameter equal to 0 gives ridge regression (at least, this functionality was not documented until recently).

Cross-validation is an estimate of the expected generalization error for each $\lambda$ and $\lambda$ can sensibly be chosen as the minimizer of this estimate. The cv.glmnet function returns two values of $\lambda$. The minimizer, lambda.min, and the always larger lambda.1se, which is a heuristic choice of $\lambda$ producing a less complex model, for which the performance in terms of estimated expected generalization error is within one standard error of the minimum. Different choices of loss functions for measuring the generalization error are possible in the glmnet package. The argument type.measure specifies the loss function.

Alternatively, the R package mgcv contains extensive possibilities for estimation with quadratic penalization including automatic selection of the penalty parameters. Methods implemented include generalized cross-validation and REML, as mentioned in a comment. More details can be found in the package authors book: Wood, S.N. (2006) Generalized Additive Models: an introduction with R, CRC.

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ You might probably want to add that cv.glmnet returns two values for $\lambda$, lambda.min and lambda.1se (“one-standard-error” rule). $\endgroup$
    – chl
    Apr 17, 2012 at 10:34
  • $\begingroup$ @chl, thanks for the suggestion. I should have added that. $\endgroup$
    – NRH
    Apr 17, 2012 at 12:00
16
$\begingroup$

This answer is MATLAB specific, however, the basic concepts should be quite similar to what you're used to with R...

In the case of MATLAB, you have the option to run lasso with cross validation enabled.

If you do so, the lasso function will report two critical parameter values

  1. The lambda value that minimizes the cross validated mean squared error
  2. The lambda value with the greatest amount of shrinkage whose CVMSE is within one standard error of the minimum.

You also get a nice little chart that you can use to inspect the relationship between lambda and CVMSE

enter image description here

In general, you'll chose a value of lambda that falls between the blue line and the green line.

The following blog posting includes some demo code based on some examples in

Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. J. Royal. Statist. Soc B., Vol. 58, No. 1, pages 267-288).

http://blogs.mathworks.com/loren/2011/11/29/subset-selection-and-regularization-part-2/

Improve this answer
$\endgroup$
10
$\begingroup$

I have had good success using effective AIC, that is using AIC with the effective degrees of freedom - see Gray JASA 87:942 1992 for effective d.f. This is implemented for $L_{2}$ penalty in the R rms package for linear and logistic models, and the rms pentrace function can be used to solve for the shrinkage coefficient that optimizes effective AIC. A case study that shows how to do differential shrinkage (e.g. more shrinkage for interactions) is Harrell et al Stat in Med 17:909, 1998.

Improve this answer
$\endgroup$
2

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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