1
$\begingroup$

I have a problem where glm.nb (R version 3.1.0, MASS version 7.3.33) converges on some data, but adding only one 0 it does not converge any more. This is the data

x <- c(3908,2729,10,803,1893,27,1312,1457,4534,3420,3,1608,903,1702,
       3041,1267,1381,3983,203,2202,1021,1550,1293,2572,1868,877,2317,
       2442, 1174,2450,1183,349)

glm.nb(x~1) converges fine, but when I run glm.nb(c(x,0)~1) it does not converge. zeroinfl(c(x,0)~1, dist="negbin") converges estimating the zero probability at 0.029 (this is roughly 1/(length(x)+1)). It seems that the problem is with theta.ml, which glm.nb uses. More precisely theta.ml(c(x,0), 1681) (1681 is the poisson estimate of mu) does not converge and this fails glm.nb.

To me adding one 0 seems like a benign thing to do (in this case), for such a dramatic change in behavior. My problem is bigger than the indication above, because I have many other pieces of data where glm.nb/theta.ml does not converge (most have more than one 0) and I am not sure what to do. I am trying to compare the negative binomial fit with its zero-inflated version (zeroinf) and am getting foiled because of this. Is the failure of glm.nb an indication that negative binomial is not appropriate? This might be the case for the examples with more 0's, but the above example with only one 0 is confusing me, because it makes me think that the problem is with the theta.ml code.

Any comments/suggestions? theta.ml seems to employ a simple iterative procedure and perhaps someone who understands it better can comment on its convergence properties.

Improve this question
$\endgroup$

1 Answer 1

Reset to default
2
$\begingroup$

Numerical solvers are very difficult to visualize. The negative binomial GLM is an exact numerical solver for the somewhat complicated negative binomial likelihood. It is easier to visualize how the likelihood "appears" for regular plain vanilla GLMs like Poisson or quasipoisson models. It's quadratic. For NB likelihoods, any two different beta-coefficients for the "rate" parameter will lead to different optimal gamma "coefficients" for the shape parameter.

The estimation process, like mixed effects models, is done with expectation-maximization. What is the big implication for NB likelihoods? The existence and uniqueness of solutions to the likelihood is not guaranteed. The 0 that you are adding is so highly inconsistent with the likelihood for the existing data that the solutions are either on the boundary of the parameter space ($\pm \infty$) or there are multiple solutions. Adding one spurious 0 is far from benign. Even for "sane" estimation procedures, adding random values is never really benign, since if they are influence points, they will drastically influence results.

Adding a zero-inflated component to the likelihood simplifies things greatly... since the original NB likelihood sans 0 converged, you are making a heirarchical model that "explains" the 0 via another parameter. Effectively, you have removed the 0 observation from the NB model with that process.

Improve this answer
$\endgroup$
8
  • $\begingroup$ Thanks, AdamO. This makes sense: I should be careful calling things benign. As far as the 0: the actual data (that came from an experiment) is with the 0 and glm.nb diverged. I just noticed that the removing it makes glm.nb converge and phrased my question as adding it.. $\endgroup$
    – user77084
    Commented May 13, 2015 at 21:05
  • $\begingroup$ In practical terms: what choice do I have? I wanted to compare the NB fit with the ZINB fit. Is it safe for me to say that because of the zero, the data is not NB distributed? I was hoping to get a more quantitative answer, but perhaps failure to converge is the best I can get. Its just that divergence seems to me an absence of evidence kind of thing. What if a different solver converges? $\endgroup$
    – user77084
    Commented May 13, 2015 at 21:36
  • $\begingroup$ I just tried fitdistr from MASS (it just uses optim) and it seems to converge. It produced some warnings (nan's produced), but the final values seem fine... These warnings appear even if I run it with an actual generated NB data, so it seems to be an issue with this solver. $\endgroup$
    – user77084
    Commented May 13, 2015 at 21:40
  • $\begingroup$ @user77084 the issue of creating a probability model that "explains the data" is a very non-scientific way of going about things. Turn your attention to the "why" of this 0 valued observation rather than the "what". Is there an unobserved process that determines whether someone experiences this event of interest? For instance, you may be measuring mRNA gene expression for a trait that some subjects in the sample do not have. With a 0-inf model, you are estimating a rate only among individuals with the trait of interest. If you are interested in fitting a negbinom model that (ctned) $\endgroup$
    – AdamO
    Commented May 13, 2015 at 21:48
  • $\begingroup$ averages over all observations, the fitdistr should theoretically alleviate some of that since it uses a more "shotgun" approach to maximizing the likelihood rather than the profile likelihood approach used in the EM algorithm of glm.nb. This one, however, can be quite imprecise for log-linear models in several parameters, but it looks like you are only estimating a rate here, so I imagine that it is an adequate approach to the problem. $\endgroup$
    – AdamO
    Commented May 13, 2015 at 21:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.