Why doesn't the zero model in a hurdle model exactly match logit result?

I am relatively new to R and I suspect that there is user error here, but I cannot figure out why the output from the logit in the hurdle model does not match the prediction of the "zero" function in the hurdle model. I offer replicable example code below.

library(AER)            # load package
library(pscl)

options(digits=15)
data("NMES1988",package = "AER")        # load the data
nmes <- NMES1988[, c(1, 6,15)]  # select certain columns; Col 1 is number of visits # :8, 13, 15, 18
mod.hurdle <- hurdle(visits ~ ., data = nmes, dist = "poisson"
, zero.dist = "binomial", link="logit")
summary(mod.hurdle)
intercept_zero<-coef(mod.hurdle, model = "zero")
hosp_zero<-coef(mod.hurdle, model = "zero")
school_zero<-coef(mod.hurdle, model = "zero")
intercept_count<-coef(mod.hurdle, model = "count")
hosp_count<-coef(mod.hurdle, model = "count")
school_count<-coef(mod.hurdle, model = "count")

predict(mod.hurdle, type = "count")

# Example entry (entry 2)
nmes[2,1:3]
entr_hosp = nmes[2,2]
entr_schools = nmes[2,3]

# Logit result
1/(1+exp(-(intercept_zero+hosp_zero*entr_hosp+school_zero*entr_schools)))

# Compare to "zero" model (does not match quite)
predict(mod.hurdle, type = "zero")

# Poisson result
exp(intercept_count+hosp_count*entr_hosp+school_count*entr_schools)

# Compare to "count" model (matches exactly)
predict(mod.hurdle, type = "count")

The results are below. The key comparison is that predict(mod.hurdle, type = "zero")2 yields 0.831938525975115, but building the logit out directly yields 0.830278531959031. It's a subtle difference in this example, but I've seen larger differences in my actual use case and more generally, I'm bothered that it's not exact. The fact that the poisson matches exactly tells me it's not just floating point error of some kind.

> mod.hurdle <- hurdle(visits ~ ., data = nmes, dist = "poisson"
+                      , zero.dist = "binomial", link="logit")
> summary(mod.hurdle)

Call:
hurdle(formula = visits ~ ., data = nmes, dist = "poisson", zero.dist = "binomial",

Pearson residuals:
Min              1Q          Median              3Q             Max
-4.198128520664 -1.211622788262 -0.484708375098  0.639181798920 18.669058090565

Count model coefficients (truncated poisson with log link):
Estimate       Std. Error  z value   Pr(>|z|)
(Intercept) 1.71087249806003 0.01968983705141 86.89115 < 2.22e-16 ***
hospital    0.21575500990929 0.00566421605875 38.09089 < 2.22e-16 ***
school      0.01164060586493 0.00172765382124  6.73781 1.6079e-11 ***
Zero hurdle model coefficients (binomial with logit link):
Estimate      Std. Error z value   Pr(>|z|)
(Intercept) 0.8347520292492 0.1145832261207 7.28512 3.2139e-13 ***
hospital    0.5102920065030 0.0894287984171 5.70613 1.1558e-08 ***
school      0.0752850525842 0.0109281136231 6.88912 5.6140e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Number of iterations in BFGS optimization: 8
Log-likelihood: -16843.9648977 on 6 Df
> intercept_zero<-coef(mod.hurdle, model = "zero")
> hosp_zero<-coef(mod.hurdle, model = "zero")
> school_zero<-coef(mod.hurdle, model = "zero")
> intercept_count<-coef(mod.hurdle, model = "count")
> hosp_count<-coef(mod.hurdle, model = "count")
> school_count<-coef(mod.hurdle, model = "count")
>
> predict(mod.hurdle, type = "count")
2
6.21694455344865
>
> # Example entry (entry 2)
> nmes[2,1:3]
visits hospital school
2      1        0     10
> entr_hosp = nmes[2,2]
> entr_schools = nmes[2,3]
>
> # Logit result
> 1/(1+exp(-(intercept_zero+hosp_zero*entr_hosp+school_zero*entr_schools)))
(Intercept)
0.830278531959031
>
> # Compare to "zero" model (does not match quite)
> predict(mod.hurdle, type = "zero")
2
0.831938525975115
>
> # Poisson result
> exp(intercept_count+hosp_count*entr_hosp+school_count*entr_schools)
(Intercept)
6.21694455344865
>
> # Compare to "count" model (matches exactly)
> predict(mod.hurdle, type = "count")
2
6.21694455344865

I think it might have to do with the zero function being a ratio of two functions, as described here, but I haven't been able to figure out how to replicate that ratio myself.

More generally, I need to be able to extract that actual functional forms for both the zero-model and the conditional-mean model, so I would like to be able to puzzle out how those are formed, so I'm starting with this logit discrepancy first.

Per a below comment about convergence, I wanted to add the following addendum: the gap between the two outcomes gets much larger when I limit to a sub-sample that includes proportionally more zeroes. I accomplish this by simply inserting nmes <- nmes[ which(nmes\$visits<=5),] after importing my data, which keeps all the zeroes, but cuts the total sample almost in half, so the zeroes become more important. The rest of the code is the same. Results:

> # Logit result
> 1/(1+exp(-(intercept_zero+hosp_zero*entr_hosp+school_zero*entr_schools)))
(Intercept)
0.74737890643342
>
> # Compare to "zero" model (does not match quite)
> predict(mod.hurdle, type = "zero")
2
0.807471618908057

Now the gap is 5 percentage points, which is very large -- too big to be convergence (especially in a 3 variable set-up, with <5000 data points), I think.

My current thinking is that the hurdle zero-function is, in fact, a ratio of two things: the probability from the logit that the observation is not zero and the probability from the Poisson that the observation is not zero (since zero should never occur conditional on exceeding the hurdle). Here is a passage in that regard from this source: But, I still can't replicate the result from the function itself:

> # Probability of 0 in Poisson
> lambda <- intercept_count+hosp_count*entr_hosp+school_count*entr_schools
> p_zero_poisson <- lambda^0*exp(-lambda)/factorial(0)
>
> # Ratio
> (1/(1+exp(-(intercept_zero+hosp_zero*entr_hosp+school_zero*entr_schools))))/(1-p_zero_poisson)
(Intercept)
1.21506968611375
>
> # Reminder of what function is finding
> predict(mod.hurdle, type = "zero")
2
0.807471618908057

As you can see, the answers are wildly different. So I'm not sure what I'm doing wrong, though I am sure it's something.

• Have you tried playing around with starting values and convergence criteria to see how much difference that makes within one or the other estimation routine? – Glen_b Jan 17 at 8:52
• Thank you for your comment @Glen_b! Unfortunately, I don't think that's it. I've added an edit to the above post to show that when I restrict to a sub-sample that has proportionally more zeroes, the gap gets much bigger. Since the number of variables is small (so estimation should be easy) and the gap is very large, I don't think that convergence issues could possibly explain it. – TGKN Jan 17 at 13:52
• I think that the issue is still most likely that the "zero" component produced by hurdle is actually a ratio between the logit and the zero component of the Poisson. But I can't recreate that components of that ratio, which is what is giving me the issue. – TGKN Jan 17 at 14:00

I have figured out the answer to my question. The concept in terms of the looking at the ratio was correct the entire time, but I was making a mistake in the math.

The ratio is meant to be: $$\begin{equation} ratio = \frac{1-f_{zero}(0|z,\gamma)}{1-f_{count}(0|x,\beta)} \end{equation}$$

So, to put it in words, the numerator is meant to be the probability of exceeding 0 coming from the zero/not zero part of the model (here, a logit). And the denominator is from the poisson component (the count model). The key reason why this knowledge did not solve my earlier confusion is that I was calculating the poisson distribution element incorrectly. The Poisson probability mass function should be from here. I.e. it should be:

$$\begin{equation} p(y|x,\beta) = \frac{e^{y*x\beta}*e^{-e^{x\beta}}}{y!} \end{equation}$$

So to get probability of 0, it is:

$$\begin{equation} p(0|x,\beta) = \frac{e^{0}*e^{-e^{x\beta}}}{0!} = e^{-e^{x\beta}} \end{equation}$$

Before, I was mistakenly thinking that $$\lambda$$ in the Poisson distribution was simply $$x\beta$$ instead of $$e^{x\beta}$$. Once I make this correction, the confusion evaporates, as you can see below.

> # Probability of 0 in logit
> p_zero_logit <- 1-1/(1+exp(-(intercept_zero+hosp_zero*entr_hosp+school_zero*entr_schools)))
>
> # Probability of 0 in Poisson
> lambda <- intercept_count+hosp_count*entr_hosp+school_count*entr_schools
> p_zero_poisson <- exp(0*lambda)*exp(-exp(lambda))/factorial(0)
>
> # Ratio
> (1-p_zero_logit)/(1-p_zero_poisson)
(Intercept)
0.807471618908057
>
> # Compare to model
> predict(mod.hurdle, type = "zero")
2
0.807471618908057

Now, the two match exactly. Thus, the confusion is resolved -- and the average of the response is simply this ratio times the average of the uncensored Poisson, as seen below:

> # Recover uncensored Poisson
> ratio <- (1-p_zero_logit)/(1-p_zero_poisson)
> uncensored_poisson_average <- exp(lambda)
>
> # Product to get full average
> uncensored_poisson_average*ratio
(Intercept)
2.09782687442114
>
> # Compare to prediction from hurdle
> predict(mod.hurdle, type = "response")
2
2.09782687442114