# Bayesian hurdle model and predictions

I fitted a Bayesian hurdle model assuming a Poisson distribution:

$P(Y_i=y_i) = \begin{cases} 1-\pi_{i} & \mbox{if }y_i=0\\[0.5em] \pi_{i} \frac{\lambda_{i}^{y_i} e^{-\lambda_i}}{y_i!(1-e^{-\lambda_i})} &\mbox{if }y_i > 0 \end{cases}$

Because the log-likelihood is separable with respect to $\pi$ and $\lambda$, I fitted two separate regressions (i.e., logistic and zero-truncated Poisson regression) using a set of predictors. In other words:

$\lambda = \log(\beta_{1}^T X) \\ \pi = \mbox{logit}(\beta_{2}^T X)$

My issue arises in trying to predict new values. In a Bayesian context, I have posterior distributions for the coefficients from the logistic and zero-truncated Poisson regressions. For a new predicted value, I will have a posterior distribution for $\pi$ and $\lambda$. To get the $E(\pi \cdot \lambda)$ is easy, I can just take the product of the posterior means. But what about getting the 2.5 and 97.5% percentiles for the credible interval around for the new predicted value?

Let's call the posterior predictive distribution $\tilde{Y}$.

\begin{align}\tilde{Y}_{\pi} &\sim\text{Bernoulli}(\pi)\\ \tilde{Y}_{\lambda} &\sim\text{Poisson}(\lambda)\\ \tilde{Y} &\sim\tilde{Y}_{\pi}*\tilde{Y}_{\lambda}\end{align}

You would need to truncate $\tilde{Y}_{\lambda}$ to be non-zero. Since you would then have a full predictive distribution, it would be easy to calculate intervals, expectations, etc.

Computation wise, this could be implemented in a number of ways. One possibility (in pseudo-code):

y_new = 0
y_pi = bernoulli_rng(pi)
if(y_pi==1)
while(y_new==0)
y_new = poisson_rng(lambda)
return y_new

• Nice response @c-r-peterson! I guess I should have clarified this, but my real issue is computational. Both $\pi$ and $\lambda$ will have distributions, so you would need to draw from both to construct the posterior predictive distribution of $\tilde{Y}$. That can be a lot, perhaps just random sampling would work? – user13317 Sep 9 '16 at 13:18
• Assuming you're using an MCMC method to calculate your posterior, then you could just use the samples from that. This would be fairly simple to implement in Stan, for example. – C.R. Peterson Sep 9 '16 at 15:06