15  Zero Inflation

In a prior chapter, we used a negative binomial regression to model the number of enforcement operations in Holland’s (2015) data.

Holland, Alisha C. 2015. “The Distributive Politics of Enforcement.” American Journal of Political Science 59 (2): 357–71. https://doi.org/10.1111/ajps.12125.

# load data
sant <- crdata::holland2015 |>
  filter(city == "santiago")

# formula
f <- operations ~ lower + vendors + budget + population

# nb regression
nb_fit <- MASS::glm.nb(
  formula = f,
  data = sant
)

# simulate from predictive distribution
nb_sims <- simulate(nb_fit, nsim = 5)

# plot
bind_cols(sant, nb_sims) |>
  pivot_longer(cols = c(operations, starts_with("sim_"))) |>
  separate(name, into = c("type", "sim_id"), sep = "_", remove = FALSE) |>
  ggplot(aes(x = value)) + 
  facet_wrap(vars(name)) +
  geom_histogram(center = 0, binwidth = 1)

Notice that the negative binomial model has too many large values. I have a subtle intuition that the negative binomial is struggling to simultaneously capture the (1) mean, (2) the large spike at zero, and (3) the relative absence of large values.

To address this, we might add a component to the model to inflate the number of zeros relative to the usual negative binomial distribution.

It is fairly common for count data to have lots of zeros. We can imagine this substantively by thinking of two different processes. One process determines whether the observation is at risk of any events at all. Another process determines the number of events, which might still be zero.¹

¹ This “two process” motivation is simply a story to motivate the model—don’t take it too literally, especially when interpreting parameter estimates.

We can model an excess of zeros using a zero-inflation model. A zero-inflation model has two steps.

1. In the first step, we have a logit model that determines whether the observation will equal zero or be a draw from some distribution.
2. In the second step, we take the draw from some distribution if needed.

15.1 Negative Binomial case

Let’s first consider zero-inflation in the negative binomial. How could we model an excess of zeros for an NB baseline distribution?

15.1.1 pmf

Let \(f(y_i; \mu_i, \theta)\) denote the pmf of a negative binomial distribution with mean \(\mu_i\) and dispersion \(\theta\).² We assume that observation \(i\) comes from

² This is the nbinom2 parameterization with overdispersion parameter \(\theta\), where \(\operatorname{Var}(y_i) = \mu_i + \mu_i^2/\theta\).

\[ y_i \sim \begin{cases} 0 & \text{with probability } \pi_i, \\ \text{NegBin}(\mu_i, \theta) & \text{with probability } 1 - \pi_i, \end{cases} \]

where \(\pi_i \in (0,1)\) is the probability of a structural zero. The remaining values follow the usual negative binomial distribution.³

³ Important note: even if an observation is not a structural zero, it still might equal zero, because the negative binomial has positive probability on zero as well. So there are two ways to get a zero.

Equivalently, the pmf of \(y_i\) can be written as

\[ f(y_i; \pi_i, \mu_i, \theta) = \begin{cases} \pi_i + (1 - \pi_i) f(0; \mu_i, \theta) & \text{if } y_i = 0 \\[6pt] (1 - \pi_i) f(y_i; \mu_i, \theta) & \text{if } y_i \neq 0 \end{cases}. \]

15.1.2 Likelihood

Create an indicator variable \(d_i = \mathbf{1}(y_i = 0)\) that equals one if \(y_i = 0\) and zero otherwise. Then for independent \(y_1, \dots, y_N\), the likelihood is

\[ L(\pi, \mu, \theta) = \prod_{i=1}^N \Bigg\{ \big[ \pi_i + (1 - \pi_i) f(0; \mu_i, \theta) \big]^{d_i} \cdot \big[ (1 - \pi_i) f(y_i; \mu_i, \theta) \big]^{1 - d_i} \Bigg\}. \]

15.1.3 Log-Likelihood

Taking logs, the log-likelihood is

\[ \ell(\pi, \mu, \theta) = \sum_{i=1}^N \Bigg\{ d_i \cdot \log\big( \pi_i + (1 - \pi_i) f(0; \mu_i, \theta) \big) + (1 - d_i) \cdot \big[ \log(1 - \pi_i) + \log f(y_i; \mu_i, \theta) \big] \Bigg\}. \]

We can imagine that \(\pi_i = \operatorname{logit}^{-1}(Z_i\gamma)\) and \(\mu_i = \exp(X_i\beta)\), with \(\theta\) modeled as in a standard negative binomial regression.

15.1.4 In R

We can fit the zero-inflated negative binomial using the glmmTMB() function in R. Create the design matrix X with formula (for \(\mu_i = \exp (X_i\beta)\)) and the design matrix Z with ziformula (for \(\pi_i = \operatorname{logit}^{-1} (Z_i\gamma)\)). The family nbinom2 specifies the negative binomial with variance \(\mu + \mu^2/\theta\) (which matches MASS::glm.nb()).

• glmmTMB() has a simulate() method to simulate from the predictive distribution.
• glmmTMB() works well with {marginaleffects}.

15.1.4.1 Constant \(\pi\)

# formula
f <- operations ~ lower + vendors + budget + population

# zinb regression
library(glmmTMB)
zinb_fit <- glmmTMB(
  formula = f,  # usual nb part
  ziformula = ~ 1, # zero-inflation (logit) w/ intercept only
  family = nbinom2,  # "theta parameterization"
  data = sant
)

Then we can use simulate() to simulate from the predictive distribution and plot it in the usual way.

# simulate from predictive distribution
zinb_sims <- simulate(zinb_fit, nsim = 5)

# plot
bind_cols(sant, zinb_sims) |>
  pivot_longer(cols = c(operations, starts_with("sim_"))) |>
  separate(name, into = c("type", "sim_id"), sep = "_", remove = FALSE) |>
  ggplot(aes(x = value)) + 
  facet_wrap(vars(name)) +
  geom_histogram(center = 0, binwidth = 1)

Notice that this looks just like the negative binomial model. The estimated probability of a structural zero is essentially zero, so the model collapses to the negative binomial.

# estimated probability of a structural zero
plogis(fixef(zinb_fit)$zi)

 (Intercept) 
2.080873e-09 

15.1.4.2 \(\pi_i = \operatorname{logit}^{-1}(Z_i\gamma)\)

# formulas
f <- operations ~ lower + vendors + budget + population
zi_f <-  ~ lower + vendors + budget + population

# zinb regression
zinb2_fit <- glmmTMB(
  formula = f,  
  ziformula = zi_f, # 🆕 model zi using covariates
  family = nbinom2,  
  data = sant
)

Then we can use simulate() to simulate from the predictive distribution and plot it in the usual way.

# simulate from predictive distribution
zinb2_sims <- simulate(zinb2_fit, nsim = 5)

# plot
bind_cols(sant, zinb2_sims) |>
  pivot_longer(cols = c(operations, starts_with("sim_"))) |>
  separate(name, into = c("type", "sim_id"), sep = "_", remove = FALSE) |>
  ggplot(aes(x = value)) + 
  facet_wrap(vars(name)) +
  geom_histogram(center = 0, binwidth = 1)

The predictive distribution for the ZINB model with covariates looks much better. Because the zero-inflated part of the model handles the spike at zero, the negative binomial part of the model can better match the remainder of the distribution.

15.1.5 IC

The BIC likes the ZINB model with covariates best, and the model without zero-inflation is a close second.

BIC(nb_fit, zinb_fit, zinb2_fit) |>
  mutate(diff_min = BIC - min(BIC),
         post_prob = exp(-0.5*diff_min)/sum(exp(-0.5*diff_min)))

          df      BIC  diff_min  post_prob
nb_fit     6 142.8430 0.5379471 0.40320644
zinb_fit   7 146.3694 4.0643077 0.06914933
zinb2_fit 11 142.3051 0.0000000 0.52764423

The AIC likes the ZINB model with covariates best and the other two models aren’t particularly close.

AIC(nb_fit, zinb_fit, zinb2_fit) |>
  mutate(diff_min = AIC - min(AIC),
         akaike_weights = exp(-0.5*diff_min)/sum(exp(-0.5*diff_min)))

          df      AIC diff_min akaike_weights
nb_fit     6 133.6848  8.16975    0.016446724
zinb_fit   7 135.6848 10.16975    0.006050412
zinb2_fit 11 125.5151  0.00000    0.977502864

15.1.6 {marginaleffects}

We can use the avg_comparisons() function to compute the average first difference as we move lower from its min to its max, averaging across the observed values of all the other covariates. For the zero-inflated models, {marginaleffects} works with the expected value of \(y_i\) itself, \(E(y_i) = (1 - \pi_i)\,\mu_i\), which combines both parts of the model.

# compute avg. fd
avg_comparisons(nb_fit, variables = list(lower = "minmax"))

 Estimate Std. Error      z Pr(>|z|)   S 2.5 % 97.5 %
      -13       15.4 -0.845    0.398 1.3 -43.1   17.1

Term: lower
Type: response
Comparison: Max - Min

avg_comparisons(zinb_fit, variables = list(lower = "minmax"))

 Estimate Std. Error      z Pr(>|z|)   S 2.5 % 97.5 %
      -13         15 -0.867    0.386 1.4 -42.3   16.4

Term: lower
Type: response
Comparison: Max - Min

avg_comparisons(zinb2_fit, variables = list(lower = "minmax"))

 Estimate Std. Error      z Pr(>|z|)   S 2.5 % 97.5 %
    -5.72       6.32 -0.905    0.365 1.5 -18.1   6.67

Term: lower
Type: response
Comparison: Max - Min

# compute avg. fd
nb_fd <- avg_comparisons(nb_fit, variables = list(lower = "minmax"))
zinb_fd <- avg_comparisons(zinb_fit, variables = list(lower = "minmax"))
zinb2_fd <- avg_comparisons(zinb2_fit, variables = list(lower = "minmax"))

# plot
list("NB" = nb_fd, 
     "ZINB; no covariates" = zinb_fd,
     "ZINB; w/ covariates" = zinb2_fd) |>
  bind_rows(.id = "model") |>
  ggplot(aes(x = estimate, xmin = conf.low, xmax = conf.high,
             y = model)) + 
  geom_errorbar(orientation = "y", width = 0.2) + 
  geom_point()

Notice that all three intervals include zero. The negative binomial and the ZINB with constant \(\pi\) give the same estimate (about \(-13\) operations). The ZINB with covariates gives an estimate less than half as large (about \(-5.7\)) and a much narrower interval.

Warning

It’s worth pointing out that the sant dataset only has 34 observations and the zero-inflated negative binomial model with covariates has 11 parameters. This is small enough that we might worry about asymptotic approximations.

15.2 Poisson case

We could similarly use a zero-inflated Poisson. It follows the same logic. The zero-inflated negative binomial is often a safer default because it also allows overdispersion in the count part, so I’ll just mention the possibility.

15.3 Hurdle models

A close alternative is a hurdle model. It also uses a logit model for zero versus positive, but it draws the positive counts from a zero-truncated distribution, so every zero comes from the logit step. glmmTMB() fits these too (e.g., family = truncated_nbinom2 with a ziformula).

15.4 Other cases

We can easily extend this logic to other distributions as well. For example, we might imagine that some variables, like trade (in dollars) between two countries, might have a continuous distribution with a spike at zero. In this case, we could imagine a zero-inflated normal, for example.