# Validity of negative binomial model: algorithm did not converge and alternation limit reached

I have strongly skewed, zero-inflated count time-series. I want to regress the counts against a continuous variable. Here's what they look like when plotted against each other year by year. (X lags behind Y by one year.) The blue print is the average followed by variance. Since variance is much higher than mean, I have selected to model them as a Negative Binomial rather than a Poisson. I am ignoring the time dimension for now. Here is what the estimation looks like:

library(MASS)
summary(glm.nb(hq_ct~ihh_emp, temp))

Call:
glm.nb(formula = hq_ct ~ ihh_emp, data = temp, init.theta = 0.2482159799,

Deviance Residuals:
Min       1Q   Median       3Q      Max
-1.7784  -1.2604  -0.7383  -0.2973   7.2707

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.469081   0.096494  -4.861 1.17e-06 ***
ihh_emp      0.043719   0.001353  32.311  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(0.2482) family taken to be 1)

Null deviance: 7686.4  on 5714  degrees of freedom
Residual deviance: 6104.3  on 5713  degrees of freedom
AIC: 33855

Number of Fisher Scoring iterations: 1

Theta:  0.24822
Std. Err.:  0.00483
Warning while fitting theta: alternation limit reached

2 x log-likelihood:  -33848.72400
Warning messages:
1: glm.fit: algorithm did not converge
2: In glm.nb(hq_ct ~ ihh_emp, temp) : alternation limit reached


I get two warnings. I ran the same model again, this time including a size variable which should explain many of the zero values. (There are two reasons for zero values: 1) Voluntary reporting, which means that more counts are reported as time goes by and that fewer were reported in the past; 2) Individual size, where smaller individuals tend to have zero counts.)

summary(glm.nb(hq_ct~ihh_emp+emp_app, temp))

Call:
glm.nb(formula = hq_ct ~ ihh_emp + emp_app, data = temp, init.theta = 0.4944747188,

Deviance Residuals:
Min       1Q   Median       3Q      Max
-4.3025  -1.3073  -0.4733   0.0794   5.0559

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.873e-01  7.484e-02   3.838 0.000124 ***
ihh_emp     7.555e-03  1.095e-03   6.901 5.15e-12 ***
emp_app     3.083e-06  3.252e-08  94.780  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(0.4945) family taken to be 1)

Null deviance: 13912.7  on 5714  degrees of freedom
Residual deviance:  5979.4  on 5712  degrees of freedom
AIC: 30492

Number of Fisher Scoring iterations: 1

Theta:  0.4945
Std. Err.:  0.0112
Warning while fitting theta: alternation limit reached

2 x log-likelihood:  -30484.0090
Warning messages:
1: glm.fit: algorithm did not converge
2: In glm.nb(hq_ct ~ ihh_emp + emp_app, temp) : alternation limit reached


I get the same two warnings again. How should I deal with them? Apart from those, does the model look good? (Note: this is only a first stage in my analysis and I plan to include panel analysis and more control variables).

Here's my data.

EDIT: I just found out that the number of iterations can be tweaked so I did, which dismissed the warnings (see below). However, the question remains as to the interpretation of the estimation: does it look good? Should I graduate to a zero-inflated model?

summary(glm.nb(hq_ct~ihh_emp, temp, control=glm.control(maxit=50)))

Call:
glm.nb(formula = hq_ct ~ ihh_emp, data = temp, control = glm.control(maxit = 50),
init.theta = 0.2482159657, link = log)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-1.7784  -1.2604  -0.7383  -0.2973   7.2707

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.469089   0.096494  -4.861 1.17e-06 ***
ihh_emp      0.043720   0.001353  32.311  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(0.2482) family taken to be 1)

Null deviance: 7686.4  on 5714  degrees of freedom
Residual deviance: 6104.3  on 5713  degrees of freedom
AIC: 33855

Number of Fisher Scoring iterations: 1

Theta:  0.24822
Std. Err.:  0.00483

2 x log-likelihood:  -33848.72400

• Your original problem seems to indeed only have been a problem of too few iterations. Regarding zero-inflated versus stand negative binomial: How badly off is the number of zeros from what your model predicts? Another perspective: If you have tons of data, often you might just as well fit the more complex model. – Björn Oct 28 '17 at 7:23
• @Björn Thanks for your feedback! Forgive me but I'm not sure I understand your question. How can I know how badly off the number of zeros is? For example, all I think I know is that some zeros are probably due to relatively greater underreporting in early years. – syre Oct 28 '17 at 7:40
• Does your software let you predict data from your fitted model? If for a really large dataset in several predicted datasets you have a small proportion of zeros and in your real data you have a huge proportion, then that's a pretty strong hint. – Björn Oct 28 '17 at 8:24
• @Björn Yes R can probably predict data but I'm not there yet. – syre Oct 28 '17 at 8:50