glmnet - compute maximal lambda value

Question

I would like to be able to compute the glmnet "lambda.max" value, in a logistic regression model. The lambda_max value stands for the smallest value for which all coefficients are zero.

According to the glmnet package vignette, "lambda.max is not given, but easily computed from the input x and y". I have unfortunately no idea on the way to compute this parameter.

Answer 1

The smallest value of lambda for which no parameters are selected may be computed by

$\max_j \frac{1}{\alpha n} \sum_{i=1}^n [Y_i - \bar Y (1- \bar Y)] X_{ij}$

See my example:

n <- 500
p <- 3
b <- c(-5,3,2,0)

X <- cbind(rep(1,n),scale(matrix(rnorm(p*n),nrow=n)))
Y <- rbinom(n,1,prob = exp(X%*%b)/(1 + exp(X%*%b)))

alpha <- .5

max( abs(t(Y - mean(Y)*(1-mean(Y))) %*% X ) )/ ( alpha * n) # largest lambda value
glmnet(x=X,y=Y,family="binomial",alpha = alpha,standardize=FALSE)$lambda[1] # largest lambda value

This comes out of the coordinate descent algorithm from this paper: Friedman, Jerome, Trevor Hastie, and Rob Tibshirani. "Regularization paths for generalized linear models via coordinate descent." Journal of statistical software 33.1 (2010): 1.

Answer 2

Here is an example:

library(glmnet)

n <- 500L

x1 <- rnorm(n, 2.0, 0.5)
x2 <- rnorm(n, -1.0, 2)
y <- factor(rbinom(n, 1L, plogis(-0.6 + 1.0 * x1 - 0.8 * x2)))

X <- matrix(c(x1, x2), ncol = 2)

mod <- glmnet(X, y, "binomial")

Now you can see the degrees of freedom and corresponding lambda by simply:

> print(mod)

Call:  glmnet(x = X, y = y, family = "binomial") 

      Df       %Dev   Lambda
 [1,]  0 -1.026e-14 0.149300
 [2,]  1  3.314e-02 0.136000
 [3,]  1  6.073e-02 0.123900
...

So in this case the coefficients are all 0 when lambda is 0.149300 (or higher). To confirm:

> coef(mod, s = 0.149300)
3 x 1 sparse Matrix of class "dgCMatrix"
                   1
(Intercept) 1.688296
V1          .       
V2          .     

> coef(mod, s = 0.136000)
3 x 1 sparse Matrix of class "dgCMatrix"
                      1
(Intercept)  1.64233062
V1           .         
V2          -0.04946726

Note that the vectors for the degrees of freedom and lambda can also be accessed via mod$df and mod$lambda, and you can also change the values of lambda which glmnet tries (if say you wanted to home in on "lambda_max") using for example:

mod <- glmnet(X, y, "binomial", lambda = seq(0.149, 0.151, by = 0.0001))