Sergio Nolazco Sergio Nolazco - 1 month ago 9
R Question

NaNs produced in negative binomial regression when using dnbinom()

I am using

dnbinom()
for writing the log-likelihood function and then estimate parameters using
mle2()
{bbmle} in R.

The problem is that I got 16 warnings for my negative binomial model, all of them NaNs produced like this one:


1: In dnbinom(y, mu = mu, size = k, log = TRUE) : NaNs produced


My code:

# data
x <- c(0.35,0.45,0.90,0.05,1.00,0.50,0.45,0.25,0.15,0.40,0.26,0.37,0.43,0.34,0.00,0.11,0.00,0.00,0.00,0.41,0.14,0.80,0.60,0.23,0.17,0.31,0.30,0.00,0.23,0.33,0.30,0.00,0.00)
y <- c(1,10,0,0,67,0,9,5,0,0,0,82,36,0,32,7,7,132,14,33,0,67,11,39,41,67,9,1,44,62,111,52,0)

# log-likelihood function
negbinglmLL = function(beta,gamma,k) {
mu= exp(beta+gamma*x)
-sum(dnbinom(y,mu=mu, size=k, log=TRUE))
}

# maximum likelihood estimator
model <- mle2(negbinglmLL, start=list(beta=mean(y), gamma= 0, k=mean(y)^2/(var(y)-mean(y))))


What do these warnings mean, and if this is a serious problem how can I avoid it?

Answer

You're not restricting the negative log-likelihood function from trying negative values of k. This probably doesn't mess up your final answer, but it's always best to avoid these kinds of warnings if you can. Two simple strategies:

  • put a lower bound on k (switching to method=L-BFGS-B)
  • fit the k parameter on the log scale, as follows:
negbinglmLL = function(beta,gamma,logk) { 
  mu= exp(beta+gamma*x)
  -sum(dnbinom(y,mu=mu, size=exp(logk), log=TRUE))
}

model <- mle2(negbinglmLL,
              start=list(beta=mean(y),
                         gamma= 0, 
                      logk=log(mean(y)^2/(var(y)-mean(y)))))

By the way, for simple problems like this you can use a formula-based shortcut as follows:

mle2(y~dnbinom(mu=exp(logmu),size=exp(logk)),
     parameters=list(logmu~x),
     start=list(logmu=0,logk=0),
     data=data.frame(x,y))

For this simple case MASS::glm.nb should also work perfectly well (but perhaps this is the simplest version of something that will get more complicated/beyond the scope of glm.nb).

Comments