1
$\begingroup$

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.

enter image description here

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, 
    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.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, 
    link = log)

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
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Björn Oct 28 '17 at 7:23
  • $\begingroup$ @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. $\endgroup$ – syre Oct 28 '17 at 7:40
  • $\begingroup$ 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. $\endgroup$ – Björn Oct 28 '17 at 8:24
  • $\begingroup$ @Björn Yes R can probably predict data but I'm not there yet. $\endgroup$ – syre Oct 28 '17 at 8:50
2
$\begingroup$

To answer your main question, whenever you get a warning like this, you should definitely increase the number of iterations or your estimates may be off.

You should also fit a 0-inflated negative binomial, since your 0s are not random (i.e., they are explained by time). If you have a relatively small number of 0s, there shouldn't be much of a difference, but if the relative number of 0s is large, the difference could be significative. I'd fit a 0-inflated negative binomial then do a simple likelihood ratio test to see which of the two has a better fit.

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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