`zeroinfl` from the `pscl` and `countreg` packages give very different results. Why?

Question

I am experimenting with the mdvis dataset from the COUNT package of R for a teaching purpose. I fitted a zero-inflated negative-binomial model using the zeroinfl function from the pscl and countreg packages. However, the results of zeroinfl from the pscl package and from countreg package differ a lot.

The models and the outputs are provided below.

zeroinfl from pscl:

mdvisit_zeroinf<-pscl::zeroinfl(numvisit ~ reform + badh + agegrp + educ + loginc | reform + badh + agegrp + educ + loginc, dist = "negbin", data = mdvis, control = zeroinfl.control(method = "BFGS", EM=F)
summary(mdvisit_zeroinf) 
Call:
pscl::zeroinfl(formula = numvisit ~ reform + badh + agegrp + educ + loginc | 
    reform + badh + agegrp + educ + loginc, data = mdvis, dist = "negbin")
Pearson residuals:
    Min      1Q  Median      3Q     Max 
-0.9667 -0.8037 -0.3597  0.3154  9.4878 

Count model coefficients (negbin with log link):
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.30316    0.56565  -0.536   0.5920    
reform      -0.09916    0.05543  -1.789   0.0736 .  
badh         1.11555    0.07517  14.841   <2e-16 ***
agegrp2      0.11376    0.06990   1.628   0.1036    
agegrp3      0.21126    0.08620   2.451   0.0143 *  
educ2        0.07605    0.07209   1.055   0.2914    
educ3       -0.03979    0.07295  -0.545   0.5854    
loginc       0.13209    0.07499   1.761   0.0782 .  
Log(theta)   0.04747    0.05815   0.816   0.4144    

Zero-inflation model coefficients (binomial with logit link):
            Estimate Std. Error z value Pr(>|z|)
(Intercept) -30.7110   772.2930  -0.040    0.968
reform       14.3274   763.7283   0.019    0.985
badh         -1.9503     3.5392  -0.551    0.582
agegrp2      10.9236   113.6292   0.096    0.923
agegrp3       9.7723   113.6558   0.086    0.931
educ2        -0.6632     1.9878  -0.334    0.739
educ3        -0.8743     1.3501  -0.648    0.517
loginc        0.5437     1.9718   0.276    0.783
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Theta = 1.0486 
Number of iterations in BFGS optimization: 106 
Log-likelihood: -4557 on 17 Df`

zeroinfl from countreg:

mdvisit_zeroinf2<-countreg::zeroinfl(numvisit ~ reform + badh + agegrp + educ + loginc | reform + badh + agegrp + educ + loginc, dist = "negbin", data = mdvis, control = zeroinfl.control(method = "BFGS", EM=F))
summary(mdvisit_zeroinf2)
Call:
countreg::zeroinfl(formula = numvisit ~ reform + badh + agegrp + educ + loginc | 
    reform + badh + agegrp + educ + loginc, data = mdvis, dist = "negbin", 
    control = zeroinfl.control(method = "BFGS", EM = F))
Pearson residuals:
    Min      1Q  Median      3Q     Max 
-0.9523 -0.8057 -0.3615  0.3106  9.6300 

Count model coefficients (negbin with log link):
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.44473    0.53752  -0.827  0.40803    
reform      -0.12962    0.05106  -2.538  0.01113 *  
badh         1.13558    0.07413  15.319  < 2e-16 ***
agegrp2      0.05851    0.06359   0.920  0.35756    
agegrp3      0.18678    0.07176   2.603  0.00925 ** 
educ2        0.06398    0.06621   0.966  0.33385    
educ3       -0.04825    0.07087  -0.681  0.49599    
loginc       0.15338    0.07056   2.174  0.02973 *  
Log(theta)   0.01232    0.04778   0.258  0.79652    

Zero-inflation model coefficients (binomial with logit link):
            Estimate Std. Error z value Pr(>|z|)  
(Intercept) -119.137     95.259  -1.251   0.2111  
reform         4.821      2.646   1.822   0.0684 .
badh         -21.461   4614.989  -0.005   0.9963  
agegrp2        9.261     28.341   0.327   0.7438  
agegrp3        3.502     27.984   0.125   0.9004  
educ2        -19.790    203.802  -0.097   0.9226  
educ3         -7.894      4.758  -1.659   0.0971 .
loginc        13.010     10.463   1.243   0.2137  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
Theta = 1.0124 
Number of iterations in BFGS optimization: 197 
Log-likelihood: -4554 on 17 Df

What could be the reason for such a difference? Results from Stata closely resemble the results of zeroinfl from the pscl package (especially, the results for the count part of the model).

I fited the models exactly similarly by specifying the same options using zeroinfl.control(). I also tried to search online to see if others have reported similar issues previously, and if answers are already available. I searched for answers also on stackoverflow and CV. But I couldn't get any.

Answer 1

Source of the problem:

The problem is that there is no zero inflation but actually less zeros than expected from a negative binomial model - espectially in some subgroups with respect to reform and agegrp.

Symptoms:

Due to the lack of zero inflation, the binary zero inflation part tries as hard as it can to produce predicted zero inflation probabilities that are very close to zero for certain subgroups. Notice especially the intercept in that part of the model that is extremely small with a large standard error.

Overall this looks similar to quasi-complete separation in binary regression models. The likelihood becomes very flat and it depends on the settings of the optimizer where exactly it stops (and these settings differ between pscl and countreg).

As the two components of the zero inflation model (count vs. binary part) cannot be estimated separately, problems in the estimation of one component can lead to problems/differences in the estimation of the other component as well.

Alternative:

The hurdle model has several advantages here: (1) It can not only deal with zero inflation but also with a lack of zeros. (2) The two components of the model can be estimated separately and hence problems cannot spill over.

Illustration:

Let's compare the basic negative binomial model with its zero-inflation and hurdle counterparts:

data("mdvis", package = "COUNT")
mdvis <- transform(mdvis,
  reform = factor(reform),
  badh = factor(badh),
  agegrp = factor(agegrp),
  educ = factor(educ)
)
library("countreg")
f <- numvisit ~ reform + badh + agegrp + educ + loginc
m <- glm.nb(f, data = mdvis)
m2 <- zeroinfl(f, dist = "negbin", data = mdvis)
m3 <- hurdle(f, dist = "negbin", data = mdvis)

In terms of the log-likelihood, the zero inflation model only improves very slightly on the basic model despite needing almost twice as many parameters. The hurdle model improves a bit more but also not much:

c(logLik(m), logLik(m2), logLik(m3))
## [1] -4560.910 -4554.146 -4550.520

In terms of both AIC and BIC the zero inflation model is the worst out of these three models. AIC slightly prefers the hurdle model while BIC prefers the basic model somewhat more clearly.

AIC(m, m2, m3)
##    df      AIC
## m   9 9139.819
## m2 17 9142.293
## m3 17 9135.040
BIC(m, m2, m3)
##    df      BIC
## m   9 9191.195
## m2 17 9239.336
## m3 17 9232.083

Looking at the rootogram as a diagnostic display for the basic model, we see that overall there are actually less observed zeros than expected from a negative binomial model. But by and large the basic model already fits reasonably well.

rootogram(m)

(Note: The version of the rootogram with confidence limits is actually produced by another package under development. The version in countreg does not have these confidence limits.)