Why does glm.nb throw a "missing value" error only on very specific inputs

user_712603 0
Viewed 6293

glm.nb throws an unusual error on certain inputs. While there are a variety of values that cause this error, changing the input even very slightly can prevent the error.

A reproducible example:

set.seed(11)
pop <- rnbinom(n=1000,size=1,mu=0.05)
glm.nb(pop~1,maxit=1000)

Running this code throws the error:

Error in while ((it <- it + 1) < limit && abs(del) > eps) { : 
  missing value where TRUE/FALSE needed

At first I assumed that this had something to do with the algorithm not converging. However, I was surprised to find that changing the input even very slightly can prevent the error. For example:

pop[1000] <- pop[1000] + 1
glm.nb(pop~1,maxit=1000)

I've found that it throws this error on 19.4% of the seeds between 1 and 500:

fit.with.seed = function(s) {
    set.seed(s)
    pop <- rnbinom(n=1000, size=1, mu=0.05)
    m = glm.nb(pop~1, maxit=1000)
}

errors = sapply(1:500, function(s) {
    is.null(tryCatch(fit.with.seed(s), error=function(e) NULL))
})

mean(errors)

I've found only one mention of this error anywhere, on a thread with no responses.

What could be causing this error, and how can it be fixed (other than randomly permuting the inputs every time glm.nb throws an error?)

ETA: Setting control=glm.control(maxit=200,trace = 3) finds that the theta.ml algorithm breaks by getting very large, then becoming -Inf, then becoming NaN:

theta.ml: iter67 theta =5.77203e+15
theta.ml: iter68 theta =5.28327e+15
theta.ml: iter69 theta =1.41103e+16
theta.ml: iter70 theta =-Inf
theta.ml: iter71 theta =NaN
statistics r
2 Answers

Edit: The code and answer has been simplified to one sample, like in the question.

Yes, theta can approach Inf in small samples and sparse data (many zeroes, small mean and large skew). I have found that fitting glm.nb fails when the data are all zeroes and returns:

Error in while ((it <- it + 1) < limit && abs(del) > eps) { : missing value where TRUE/FALSE needed

The following code simulates small samples with a small mean and theta. To prevent the loop from crashing, glm.nb is not fitted when the data are all zeroes.

en1 <- 10
mu1 <- 0.5
size1 <- 0.5
temp <- matrix(nrow=10000, ncol=2)
# theta == Inf is rare so use a large number of reps
for (r in 1:10000){
  dat1 <- rnbinom(n=en1, size=size1, mu=mu1)
  temp[r, 1:2] <- c(mean(dat1), ifelse(max(dat1)!=0, glm.nb(dat1~1)$theta, NA))
}
temp <- as.data.frame(temp)
  names(temp) <- c("mean1","theta1")
  temp[which(is.na(temp$theta1)),]
  # note that it's rare to get all zeroes in the sample
  sum(is.na(temp$theta1))/dim(temp)[1]
# a log scale helps see what's happening
with(temp, plot(mean1, log10(theta1)))
  # estimated thetas should equal size1 = 0.5
  abline(h=log10(0.5), col="red")
  text(2.5, 5, "n1 = n2 = 10", col="red", cex=2, adj=1)
  text(1, 4, "extreme thetas", col="red", cex=2)

See that estimated thetas can be extremely large when the sample size is small (in the first plot below): Lesson learnt: don't expect high quality results from glm.nb for small samples and sparse data; get larger samples (e.g. in the second plot below). Example 1 Example 2

user_195530090
Related