Error when running glmnet in multinomial

Question

The problem mentioned in this question is fixed in version 1.7.3 of the R package glmnet.

I'm having some problems running glmnet with family=multinomial, and was wondering has encountered something similar or might be able to tell me what I'm doing wrong.

When I put my own dummy data in, the error "Error in apply(nz, 1, median) : dim(X) must have a positive length" gets reported when I run cv.glmnet, which apart from saying "it didn't work" wasn't hugely informative to me.

y=rep(1:3,20) #=> 60 element vector
set.seed(1011)
x=matrix(y+rnorm(20*3*10,sd=0.4),nrow=60) # 60*10 element matrix
glm = glmnet(x,y,family="multinomial")   #=> returns without error
crossval = cv.glmnet(x,y,family="multinomial")   #=> Error in apply(nz, 1, median) : dim(X) must have a positive length
crossval = cv.glmnet(x,y,family="multinomial",type.measure="class")   #=> Error in apply(nz, 1, median) : dim(X) must have a positive length
crossval = cv.glmnet(x,y,family="multinomial",type.measure="mae")   #=> Error in apply(nz, 1, median) : dim(X) must have a positive length
cvglm = cv.glmnet(x,y,family="multinomial",lambda=2)   #=> Error in apply(nz, 1, median) : dim(X) must have a positive length

Here's a visual description of the problem I was trying to get glmnet to solve, if that helps:

my_colours = c('red','green','blue')
plot(x[,1],x[,2],col=my_colours[y])

I'm able to run the example code from the package docs, which makes me suspcious that I'm either misunderstanding something or that there is a bug in glmnet.

library(glmnet)
set.seed(10101)
n=1000;p=30
x=matrix(rnorm(n*p),n,p) #=> 1000*30 element matrix
beta3=matrix(rnorm(30),10,3)
beta3=rbind(beta3,matrix(0,p-10,3))
f3=x%*% beta3
p3=exp(f3)
p3=p3/apply(p3,1,sum)
g3=rmult(p3) #=> 1000 element vector
set.seed(10101)
cvfit=cv.glmnet(x,g3,family="multinomial")

This is using R version 2.13.1 (2011-07-08) and glmnet 1.7.1, though I can generate the same problem on R 2.14.1. Any ideas people?

Answer 1

There is a subtle bug.

What is happening is the following: In your artificial data set, the three group means are on a line, and with the relatively small standard deviation used, the three groups become linearly separable in your 10-dimensional space. As a consequence, all parameters related to the second group are estimated to 0 for all $\lambda$. Check

coef(glm)

Internally in cv.glmnet there is a call to predict to determine for each $\lambda$ the number of non-zero coefficients. Try

predict(glm, type = "nonzero")

The structure is, from reading the cv.glmnet code, supposed to be a list of lists, but the second entry in the list is NULL, and not a list! This causes the error. It happens in this block of code from cv.glmnet

if (inherits(glmnet.object, "multnet")) {
    nz = predict(glmnet.object, type = "nonzero")
    nz = sapply(nz, function(x) sapply(x, length))
    nz = ceiling(apply(nz, 1, median))
}

The result returned from the two nested sapply calls is not a matrix as expected in the last call of apply. This generates the error.

It might be very unlikely to run into the error in practice, but the code should of course be robust to extreme cases. You should report the problem to the maintainer, Trevor Hastie (his email is listed at the link).