glm.nb fails to converge when adding one zero 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.
 A: 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.
