2
$\begingroup$

I am trying to test effects of 3 predictors on overdispersed count data with many zeros, and a Vuong test suggested that a zero-inflated neg. binomial model would fit better than a negative binomial model; however, I'm at a loss to interpret the outcome, also as the zeroinfl-model apparently produces opposed results to those of the glm.nb model (which again seem to be more intuitive, looking at the raw data only). So, for instance below, from the zeroinfl model it seems to me that predictor 'ha' had a significant negative effect on the counts, but in the glm.nb model, this is positive. I looked up the interpretation of zeroinfl in several sources, but ended up utterly confused. Any help with this would be greatly appreciated!

zeroinfl(formula = count ~ ha + re + an, data = data, dist = "negbin")

Pearson residuals:
    Min      1Q  Median      3Q     Max 
-0.8027 -0.4015 -0.3759 -0.1546  4.5059 

Count model coefficients (negbin with log link):
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   4.9183     0.2530  19.443  < 2e-16 ***
ha       -0.2135     0.2927  -0.729  0.46584    
re        0.6506     0.2787   2.335  0.01955 *  
an       -0.2620     0.3114  -0.842  0.40002    
Log(theta)    0.6663     0.2542   2.621  0.00877 ** 

Zero-inflation model coefficients (binomial with logit link):
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  1.49631    0.48655   3.075  0.00210 **
ha      -1.67048    0.51640  -3.235  0.00122 **
re      -0.17539    0.51943  -0.338  0.73561   
an      -0.03723    0.60535  -0.061  0.95096   

versus

glm.nb(formula = count ~ ha + re + an, data = data,  
    init.theta = 1264914.211, link = log)
    Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-15.927   -9.712   -6.934   -2.476   48.338  
    Coefficients:
Estimate Std. Error z value Pr(>|z|)     
(Intercept)  3.17986    0.03335  95.362  < 2e-16 *** 
ha       0.98945    0.03035  32.598  < 2e-16 *** 
re       0.67368    0.03212  20.977  < 2e-16 *** 
an      -0.15229    0.03675  -4.144 3.42e-05 ***
$\endgroup$
1
$\begingroup$

When you fit MASS:negbin, and there are too many zeros and it cannot estimate the dispersion parameter well and these throws all your estimate off. You can see that your standard errors for the coefficients are too small.

In negative binomial, the variance is parameterized as :

var = mean + (mean^2) / theta

If you look at the output of theta, = 1264914 you get for the negative binomial, this makes it almost a poisson, which doesn't quite makes sense, especially when you have a lot of zeros, indicating extra dispersion.

I can demonstrate a very simple example below:

set.seed(999)
data = data.frame(ha=runif(300,min=1,max=2),re=runif(300,min=1,max=2)) 
data$counts = rnbinom(300,mu=exp(3 + 1*data$ha + 0.6*data$re),size=1) 
summary(MASS::glm.nb(counts ~ ha+re,data=data))

    MASS::glm.nb(formula = counts ~ ha + re, data = data, init.theta = 1.049918539, 
    link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.2162  -1.0459  -0.3770   0.4384   2.4797  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   2.9833     0.4222   7.066 1.59e-12 ***
ha            0.9325     0.1967   4.740 2.14e-06 ***
re            0.6971     0.1972   3.536 0.000407 ***

We get back the coefficients, and theta we expect.

Now we assume some of the data will be missing at random and they are recorded as zero, and we fit it on this data:

data$counts[sample(nrow(data),120)] = 0
summary(MASS::glm.nb(counts ~ ha+re,data=data))

MASS::glm.nb(formula = counts ~ ha + re, data = data, init.theta = 518053.5004, 
    link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-25.691  -15.071   -7.891    4.156   63.963  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.89495    0.03942   48.07   <2e-16 ***
ha           1.02324    0.01737   58.90   <2e-16 ***
re           0.97123    0.01756   55.31   <2e-16 ***

You can see I get a very high theta, which in turns throws the std.error estimate off. If we fit a zeroinflated negbin, it will more or less return you the correct estimate:

library(pscl)
summary(zeroinfl(counts ~ ha + re | 1, data = data, dist = "negbin"))

Count model coefficients (negbin with log link):
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   2.3368     0.5407   4.322 1.55e-05 ***
ha            1.0493     0.2499   4.200 2.67e-05 ***
re            0.9854     0.2402   4.102 4.09e-05 ***
Log(theta)    0.1280     0.1011   1.266    0.205    

Zero-inflation model coefficients (binomial with logit link):
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -0.4128     0.1184  -3.486 0.000491 ***

The theta estimate is exp(0.1280) = 1.13 which is what we started with. Coefficients are a bit off, but I guess this is expected giving the amount of non zero data available. This is a simplified model, but it shows that the negative binomial is dependent on you estimating the theta correctly. So you should be careful about interpreting the negative binomial.

In your case, I think most likely a increase in ha increases your counts and decreases the probability of getting zero, giving you the result with zero-inflated. So based on the nature of your data, if you can say most like it has nothing to do with it, I would suggest fitting it like my above example, assuming a flat logit for the zero-inflated part. Then you will see the theta

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you so much for your elaborate answer to my problem! It really helped me a lot! $\endgroup$ – Don Joe Apr 18 at 6:35
  • $\begingroup$ you're welcome ! glad it was useful for you $\endgroup$ – StupidWolf Apr 18 at 16:25

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.