CODEWITHSUNDEEP
rbigdataglmpoissonglmnet
Eli Holmes
Eli Holmes

Reputation: 676

Lack of convergence of glmnet when lambda=0 for family="poisson"

While getting a handle on glmnet versus glm, I ran into convergence problems for lambda=0 and family="poisson". My understanding is that with lambda=0 (and alpha=1, the default), the answers should be essentially the same.

Below is code changed slightly from the poisson example on the glmnet help page (?glmnet). The only change is that nzc = p so that all variables are in the true model

N=1000; p=50
nzc=p
x=matrix(rnorm(N*p),N,p)
beta=rnorm(nzc)
f = x[,seq(nzc)]%*%beta
mu=exp(f)
y=rpois(N,mu)

#With lambda=0 glmnet throws the convergence error shown below
fit=glmnet(x,y,family="poisson",lambda=0)

#It works with default lambda passed in
# but estimates are quite different from glm.
fit=glmnet(x,y,family="poisson") #use default lambdas
fit2=glm(y~x,family="poisson")
plot(coef(fit2)[2:(p+1)], 
     coef(fit,s=min(fit$lambda))[2:(p+1)],
     xlab="glm",ylab="glmnet")
abline(0,1)

#works fine with gaussian response and lambda=0 or default lambda
#glm and glmnet identical
mu = f
y=rnorm(N,mu)
fit=glmnet(x,y,family="gaussian",lambda=0)
fit2=glm(y~x)
plot(coef(fit2)[2:(p+1)], coef(fit)[2:(p+1)])
abline(0,1)

Here's the error message

Warning messages:
1: from glmnet Fortran code (error code -1); Convergence for 1th lambda value not reached after maxit=100000 iterations; solutions for larger lambdas returned 
2: In getcoef(fit, nvars, nx, vnames) :an empty model has been returned; probably a convergence issue

Updated: The problem seems to be with the intercept being estimated by glmnet when family="poisson" and not related to the setting of lambda per se.

fit=glmnet(x,y,family="poisson")
#intercept should be close to 0
coef(fit)[1,]
#but it is huge
#passing in intercept=FALSE however generates the convergence error again
fit=glmnet(x,y,family="poisson", intercept=FALSE)

Upvotes: 3

Views: 4246

Answers (2)

Frank
Frank

Reputation: 1

You get the error because you try to pass lambda = 0 to glmnet. If you want to select the coefficients from glmnet for lambda = 0, you could use:

coef(fit, s=0)

This automatically selects the last (smallest) value of lambda. I guess you've basically done that already though, with s = min(fit$lambda). If you want to go even smaller than that you might have to manually put in a lambda sequence, but this is a little bit tricky (glmnet seems a little bit stubborn about its lambda's).

Also keep in mind that there might be some bias in glmnet, so it could be slightly different from the results of glm.

Upvotes: 0

IRTFM
IRTFM

Reputation: 263362

I think you are confused about lambda and alpha. alpha is the penalization factor which is set to 0 will give you ridge regression. Typically it is set to something between 0.1 and 1. lambda is typically not set, and there is a warning on the help page NOT to set it to a single value:

WARNING: use with care. Do not supply a single value for lambda

I don't know why you think a lasso penalty should be the same as an unpenalized Poisson model. The whole point of a penalized model is to be less subject to the biases and constraints of an ordinary regression model.

Upvotes: 1

Related Questions