Consider a hurdle model predicting count data y
from a normal predictor x
:
set.seed(1839)
# simulate poisson with many zeros
x <- rnorm(100)
e <- rnorm(100)
y <- rpois(100, exp(-1.5 + x + e))
# how many zeroes?
table(y == 0)
FALSE TRUE
31 69
In this case, I have count data with 69 zeros and 31 positive counts. Nevermind for the moment that this is, by definition of the data-generation procedure, a Poisson process, because my question is about hurdle models.
Let's say I want to handle these excess zeros by a hurdle model. From my reading about them, it seemed like hurdle models aren't actual models per se—they are just doing two different analyses sequentially. First, a logistic regression predicting whether or not the value is positive versus zero. Second, a zero-truncated Poisson regression with only including the non-zero cases. This second step felt wrong to me because it is (a) throwing away perfectly good data, which (b) could lead to power issues since much of the data are zeros, and (c) basically not a "model" in and of itself, but just sequentially running two different models.
So I tried a "hurdle model" versus just running the logistic and zero-truncated Poisson regression separately. They gave me identical answers (I'm abbreviating the output, for brevity's sake):
> # hurdle output
> summary(pscl::hurdle(y ~ x))
Count model coefficients (truncated poisson with log link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.5182 0.3597 -1.441 0.1497
x 0.7180 0.2834 2.533 0.0113 *
Zero hurdle model coefficients (binomial with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.7772 0.2400 -3.238 0.001204 **
x 1.1173 0.2945 3.794 0.000148 ***
> # separate models output
> summary(VGAM::vglm(y[y > 0] ~ x[y > 0], family = pospoisson()))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.5182 0.3597 -1.441 0.1497
x[y > 0] 0.7180 0.2834 2.533 0.0113 *
> summary(glm(I(y == 0) ~ x, family = binomial))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.7772 0.2400 3.238 0.001204 **
x -1.1173 0.2945 -3.794 0.000148 ***
---
This seems off to me since many different mathematical representations of the model include the probability that an observation is non-zero in the estimation of positive count cases, but the models I ran above completely ignore one another. For example, this is from Chapter 5, page 128 of Smithson & Merkle's Generalized Linear Models for Categorical and Continuous Limited Dependent Variables:
...Second, the probability that $y$ assumes any value (zero and the positive integers) must equal one. This is not guaranteed in Equation (5.33). To deal with this issue, we multiply the Poisson probability by the Bernoulli success probability $\pi$.
These issues require us to express the above hurdle model as
$$ P(Y=y|\boldsymbol{x,z,\beta,\gamma}) = \begin{cases} 1-\hat\pi &\text{for } y=0 \\ \hat\pi\times\frac{\exp(-\hat\lambda)\hat\lambda^y/y!}{1-\exp(-\hat\lambda)} &\text{for } y=1,2,\ldots \end{cases} \tag{5.34} $$ where $\hat\lambda=\exp(\boldsymbol{x\beta})$, $\hat\pi = {\rm logit}^{-1}(\boldsymbol{z\gamma})$, $\boldsymbol x$ are the covariates for the Poisson model, $\boldsymbol z$ are the covariates for the logistic regression model, and $\hat{\boldsymbol{\beta}}$ and $\hat{\boldsymbol{\gamma}}$ are the respective regression coefficients....
By doing the two models completely separate from one another—which seems to be what hurdle models do—I don't see how $\hat{\pi}$ is incorporated into the prediction of positive count cases. But based on how I was able to replicate the hurdle
function by just running two different models, I don't see how $\text{logit}^{-1}(z\hat{\gamma})$ plays a role in the truncated Poisson regression at all.
Am I understanding hurdle models correctly? They seem two be just running two sequential models: First, a logistic; Second, a Poisson, completely ignoring cases where $y = 0$. I would appreciate if someone could clear-up my confusion with the $\hat{\pi}$ business.
If I am correct that that is what hurdle models are, what is the definition of a "hurdle" model, more generally? Imagine two different scenarios:
Imagine modeling competitiveness of electoral races by looking at competitiveness scores (1 - (winner's proportion of vote - runner up's proportion of vote)). This is [0, 1), because there are no ties (e.g., 1). A hurdle model makes sense here, because there is one process (a) was the election uncontested? and (b) if it wasn't, what predicted competitiveness? So we first do a logistic regression to analyze 0 vs. (0, 1). Then we do beta regression to analyze the (0, 1) cases.
Imagine a typical psychological study. Responses are [1, 7], like a traditional Likert scale, with a huge ceiling effect at 7. One could do a hurdle model that's logistic regression of [1, 7) vs. 7, and then a Tobit regression for all cases where observed responses are < 7.
Would it be safe to call both of these situations "hurdle" models, even if I estimate them with two sequential models (logistic and then beta in the first case, logistic and then Tobit in the second)?
pscl::hurdle
, but it looks the same in Equation 5 here: cran.r-project.org/web/packages/pscl/vignettes/countreg.pdf Or perhaps I'm still missing something basic that would make it click for me? $\endgroup$hurdle()
. In our paired/ vignette, we try to emphasize the more general building blocks, though. $\endgroup$