23 Which count data model to use

Published

October 10, 2017

23.1 A comparison of various count data models with extra zeros

In empirical studies, data sets with a lot of zeros are often hard to model. There are various models to deal with it: zero-inflated Poisson model, Negative Binomial (NB)model, hurdle model, etc.

Here we are following a zero-inflated model’s thinking: model the data with two processes. One is a Bernoulli process, the other one is a count data process (Poisson or NB).

We’d like to see, in this simulation exercise, how different models perform with changes of sample size and percentage of zeros (we expect the less zero, the better a plain Poisson model would perform). Therefore we vary sample size $n$ and an indicator of how much percentage of zeros in the data $\theta$.

For the count data process ($y_c$): \[ y_c = \lfloor \exp(2 x + u) \rfloor, \qquad u \sim N(0, \sigma^2) \tag{23.1}\]

Note what that is and is not. It is a floored lognormal variable, not a draw from a Poisson distribution, and not normal either. The flooring matters: it makes the conditional mean only approximately $\exp(2x + \sigma^2/2)$, so the log-linear index is not exactly $2x$ and “bias relative to 2” is not a clean common target across the models compared below.

For the Bernoulli process ($y_b$):

\[ z_1 = 4 z + \theta \tag{23.2}\] \[ logit(y_b) = z_1 \tag{23.3}\] \[ p_y = \frac{e^{z_1}}{1+e^{z_1}} \tag{23.4}\]

where $y_b \sim Bernoulli(p_y)$ is the realized 0/1 draw. Combining these two processes:

\[ y = y_c \cdot y_b \tag{23.5}\] so $y = y_c$ when $y_b=1$ and $y = 0$ when $y_b=0$.

23.1.1 Zero-inflated Poisson models

A zero-inflated Poisson needs specifying both the binary process and the count process correctly. Often than not, we don’t have a model for the binary process. Many people simply use the same explanatory variables for both processes. We simulate both situations. Case 1: suppose we observe $z$, and case 2: suppose we don’t observe $z$. In the graph below, they are labeled zip1 and zip2.

23.1.2 Poisson model

A plain Poisson model returns a consistent estimator for the coefficients, with or without Poisson-distributed data. We expect Poisson model’s performance improve with sample size. Note that the standard errors from a Poisson model needs adjustment, which we do not discuss in this post.

23.1.3 NB model

NB model is used widely to handle “overdispersion” problem. That is, the variance far exceeds the mean, therefore the Poisson model is considered inappropriate. NB model addresses that by allowing an extra parameter. However, many people also use it to model “extra zero” situation, we’ll see in our simulation it may not be better than a plain Poisson model.

23.1.4 Log-linear model

What about an OLS model with $log(y+1)$?

23.1.5 hurdle model

A hurdle model models the zero’s and other values separately; that is, the zero’s are from a binomial process only, the other positive values are from a truncated count data process. We assume here, in the simulation, we don’t observe $z$. Therefore, $x$ is determining both binary and count processes. In the graph below, it’s labeled hurdle.

Code

library(MASS)
library(pscl)
library(parallel)
set.seed(666)

gen.sim <- function(df) {
  nobs <- df["nobs"]; th <- df["th"]
  z    <- rnorm(nobs, 0, 1)
  x    <- rnorm(nobs, 0, 1)
  u    <- rnorm(nobs, 0, 1)
  y.count  <- floor(exp(2 * x + u))
  z1       <- 4 * z + th
  prob     <- plogis(z1)
  y.logit  <- rbinom(nobs, size = 1, prob = prob)
  y        <- ifelse(y.logit == 1, y.count, y.logit)

  m1 <- zeroinfl(y ~ x | z)
  m4 <- zeroinfl(y ~ x | x)
  m2 <- glm(y ~ x, family = "poisson")
  m3 <- lm(log(y + 1) ~ x)
  m5 <- glm.nb(y ~ x)
  m6 <- hurdle(y ~ x)

  c(zip1       = coef(m1, "count")["x"] - 2,
    poisson     = coef(m2)["x"] - 2,
    log.linear  = coef(m3)["x"] - 2,
    zip2        = coef(m4, "count")["x"] - 2,
    nb          = coef(m5)["x"] - 2,
    hurdle      = coef(m6, "count")["x"] - 2)
}

data.grid <- expand.grid(
  nobs = ceiling(exp(seq(4, 9, 1)) / 100) * 100,
  nsim = seq(1, 100, 1),
  th   = seq(-4, 4, 2)
)

cl      <- makeCluster(detectCores() - 1)
clusterEvalQ(cl, { library(MASS); library(pscl) })
clusterExport(cl, "gen.sim")
results <- parApply(cl, data.grid, 1, gen.sim)
stopCluster(cl)

forshiny <- cbind(data.grid, t(results))
write.csv(forshiny, "results.csv", row.names = FALSE)

Count data models can be used even if data is not “counts”; for example, some non-negative non-integer numbers. In fact, Poisson model is consistent even if data is not Poisson-distributed, if the model specification is correct on modeling the log of expected counts. We simulate both scenarios: Case 1, data is generated from a Poisson process. Case 2, data is a floored lognormal – $\lfloor \exp(2x+u) \rfloor$ with $u$ normal – so it is non-negative and integer-valued without being a count from any count distribution, and we use count data models on it anyway. (Earlier versions of this text called case 2 “Normal”; the variable being modelled is not normal, only the $u$ inside the exponential is.) The above code is for case 2.

We simulate 100 times with $\theta$ ranging from -4 to 4, lower number means higher percentage of zeros; number of observations from $e^4$ to $e^9$ (roughly from 50 to 8000).

Since there are many simulations, we use base R’s parallel library to speed things up.

For raw code, please visit case1: poisson and case2: normal. The code above updates the original snowfall-based version to use base R parallel.

Status of these simulation findings. The grid above is retained for its code rather than as current evidence, and the ranking it once carried should not be read as verified. Two things undercut it. The metric for the log-linear row was wrong in the original run – see the note below – and the case-2 data-generating process is a floored lognormal whose log-linear index is not exactly $2x$, so subtracting 2 from each model’s slope does not compare like with like. The chunk is marked eval: false, so nothing here regenerates it.

What the run suggested, and what remains plausible on grounds independent of it: plain Poisson QMLE holds up well in large samples even with unmodelled excess zeros when the conditional mean is correctly specified; zero-inflated Poisson does better when the zero process is known; and log-linear regression on $\log(y+1)$ targets a different quantity altogether. For evidence rather than recollection, the worked example in the next section is the part of this chapter that actually runs.

(Note: the log.linear bias line above previously read exp(coef(m3)["x"]) - 2 instead of coef(m3)["x"] - 2, which — since the log-linear coefficient is already on the log scale — manufactured a spurious ≈5.4 apparent bias regardless of the model’s actual performance. That has been corrected here to match the other rows. Log-linear regression on log(y+1) is still expected to be biased relative to Poisson-family models in a zero-inflated count setting, for real reasons (the log transform’s nonlinearity near zero, and the fact that $E[\log(y+1)] \neq \log(E[y]+1)$), but the magnitude of “consistently the worst” above reflects the original (buggy) run and has not been re-verified against a corrected simulation.)

23.2 Quick worked example

All five models on a single small dataset (n = 1000, ~40% zeros).

The DGP is a zero-inflated Poisson built in two independent pieces. We draw $x \sim N(0,1)$ and $z \sim N(0,1)$, independent of each other. The count part is $y_c \sim \text{Poisson}(\exp(1.5 + 0.8x))$ and the zero part is $y_b \sim \text{Bernoulli}(\Lambda(1 + 2z))$, with the observed outcome $y = y_c \cdot y_b$. So the true slope on $x$ is 0.8, and — this is the feature that drives the results — the zero process depends on $z$ alone, never on $x$.

Code

library(MASS)
library(pscl)
set.seed(42)

n  <- 1000
x  <- rnorm(n)
z  <- rnorm(n)                          # zero-inflation driver
mu <- exp(1.5 + 0.8 * x)               # count mean
y_c <- rpois(n, lambda = mu)
y_b <- rbinom(n, 1, plogis(1 + 2 * z)) # 1 = non-zero
y   <- y_c * y_b                        # zero-inflated count

cat("Zero fraction:", round(mean(y == 0), 2), "\n")

Zero fraction: 0.39

Code

cat("Mean (non-zero):", round(mean(y[y > 0]), 2), "\n\n")

Mean (non-zero): 6.79

Code

# 1. Poisson
m_pois <- glm(y ~ x, family = poisson)

# 2. Negative binomial
m_nb <- glm.nb(y ~ x)

# 3. Zero-inflated Poisson (ZIP) — correct zero model
m_zip <- zeroinfl(y ~ x | z)

# 4. ZIP — wrong zero model (only x)
m_zip2 <- zeroinfl(y ~ x | x)

# 5. Hurdle
m_hurdle <- hurdle(y ~ x)

# 6. Log-linear OLS on log(y + 1)
m_loglin <- lm(log(y + 1) ~ x)

# Compare x coefficient (true = 0.8)
coefs <- c(
  Poisson    = coef(m_pois)["x"],
  NB         = coef(m_nb)["x"],
  `ZIP (correct z)` = coef(m_zip, "count")["x"],
  `ZIP (wrong)` = coef(m_zip2, "count")["x"],
  Hurdle     = coef(m_hurdle, "count")["x"],
  `Log-linear log(y+1)` = coef(m_loglin)["x"]
)

knitr::kable(
  data.frame(Model = names(coefs), `x coefficient` = round(coefs, 3),
             Bias = round(coefs - 0.8, 3)),
  caption = "True x coefficient = 0.8"
)

True x coefficient = 0.8
	Model	x.coefficient	Bias
Poisson.x	Poisson.x	0.850	0.050
NB.x	NB.x	0.824	0.024
ZIP (correct z).x	ZIP (correct z).x	0.803	0.003
ZIP (wrong).x	ZIP (wrong).x	0.802	0.002
Hurdle.x	Hurdle.x	0.803	0.003
Log-linear log(y+1).x	Log-linear log(y+1).x	0.458	-0.342

The sample comes out 39% zeros with a non-zero mean of 6.79. The estimated slopes are 0.850 for Poisson, 0.824 for negative binomial, 0.803 for ZIP with the correct zero model, 0.802 for ZIP with the wrong one, 0.803 for the hurdle model, and 0.458 for log-linear OLS.

The first five models recover the true x slope of 0.8 reasonably well; the log-linear row is included to make the earlier claim checkable, and needs reading with care — see the note after the table. This is the key point of this DGP: the zero-inflation process depends on z, which is generated independently of x, so omitting the zero model distorts the zero mass and the intercept, not the slope on x. The hurdle model and even the ZIP with a misspecified (x-only) zero model still recover the slope closely; the plain Poisson shows a small upward bias in this particular sample (n = 1000), but that is finite-sample noise, not a systematic property of the model. Since the zero-inflation driver z is independent of x, the conditional mean here is correctly specified for Poisson QMLE — increasing n from 1,000 to 1,000,000 shrinks the Poisson x-coefficient bias from about +0.04 to effectively 0, confirming consistency. The unmodeled excess zeros create overdispersion, but overdispersion under a correctly-specified conditional mean affects only the standard errors (which should be corrected with a robust/sandwich estimator), not the consistency of the slope.

About the log-linear row: its coefficient is not estimating the same thing as the others, which is the real reason to avoid it rather than any particular bias number. The count models target $\log E[y \mid x]$, whereas OLS on $\log(y+1)$ targets $E[\log(y+1) \mid x]$, and $E[\log(y+1)] \neq \log(E[y]+1)$ by Jensen’s inequality. The $+1$ offset makes the discrepancy worst exactly where the mass of zeros is. So comparing its coefficient to 0.8 is an apples-to-oranges comparison, and the honest statement is that log-linear answers a different question — not that it answers this one badly by some measured amount. A misspecified zero process biases the count slope only when the zero process is correlated with the count covariates — which is not the case here. (Change z to depend on x and the omitted-zero-process slope bias appears.)

23.3 Mixed count models with glmmTMB

When data has both overdispersion and random effects (clustered counts, panel data, repeated measures), glmmTMB handles the full range in one package:

Code

library(glmmTMB)

# Zero-inflated negative binomial with random intercept per subject
m_glmmtmb <- glmmTMB(
  count ~ treatment + time + (1 | subject),   # count model
  ziformula = ~ treatment,                     # zero-inflation model
  family = nbinom2,
  data = your_data
)

summary(m_glmmtmb)

# Conditional vs marginal predictions
predict(m_glmmtmb, type = "conditional")  # excluding zero-inflation
predict(m_glmmtmb, type = "response")     # including zero-inflation

glmmTMB supports: Poisson, NB1, NB2, zero-inflated variants of all of the above, hurdle models, beta-binomial, tweedie, and more. The ziformula argument takes a one-sided formula for the zero-inflation component — use ~ 1 for a constant zero-inflation rate, or covariates to model what drives the structural zeros.

--- title: "Which count data model to use" date: "2017-10-10" --- ## A comparison of various count data models with extra zeros In empirical studies, data sets with a lot of zeros are often hard to model. There are various models to deal with it: zero-inflated Poisson model, Negative Binomial (NB)model, hurdle model, etc. Here we are following a zero-inflated model's thinking: model the data with two processes. One is a Bernoulli process, the other one is a count data process (Poisson or NB). We'd like to see, in this simulation exercise, how different models perform with changes of sample size and percentage of zeros (we expect the less zero, the better a plain Poisson model would perform). Therefore we vary sample size $n$ and an indicator of how much percentage of zeros in the data $\theta$. For the count data process ($y_c$): $$ y_c = \lfloor \exp(2 x + u) \rfloor, \qquad u \sim N(0, \sigma^2) $$ {#eq-count-data-1} Note what that is and is not. It is a floored lognormal variable, not a draw from a Poisson distribution, and not normal either. The flooring matters: it makes the conditional mean only approximately $\exp(2x + \sigma^2/2)$, so the log-linear index is not exactly $2x$ and "bias relative to 2" is not a clean common target across the models compared below. For the Bernoulli process ($y_b$): $$ z_1 = 4 z + \theta $$ {#eq-count-data-2} $$ logit(y_b) = z_1 $$ {#eq-count-data-3} $$ p_y = \frac{e^{z_1}}{1+e^{z_1}} $$ {#eq-count-data-4} where $y_b \sim Bernoulli(p_y)$ is the realized 0/1 draw. Combining these two processes: $$ y = y_c \cdot y_b$$ {#eq-count-data-5} so $y = y_c$ when $y_b=1$ and $y = 0$ when $y_b=0$. ### Zero-inflated Poisson models A zero-inflated Poisson needs specifying both the binary process and the count process correctly. Often than not, we don't have a model for the binary process. Many people simply use the same explanatory variables for both processes. We simulate both situations. Case 1: suppose we observe $z$, and case 2: suppose we don't observe $z$. In the graph below, they are labeled zip1 and zip2. ### Poisson model A plain Poisson model returns a consistent estimator for the coefficients, with or without Poisson-distributed data. We expect Poisson model's performance improve with sample size. Note that the standard errors from a Poisson model needs adjustment, which we do not discuss in this post. ### NB model NB model is used widely to handle "overdispersion" problem. That is, the variance far exceeds the mean, therefore the Poisson model is considered inappropriate. NB model addresses that by allowing an extra parameter. However, many people also use it to model "extra zero" situation, we'll see in our simulation it may not be better than a plain Poisson model. ### Log-linear model What about an OLS model with $log(y+1)$? ### hurdle model A hurdle model models the zero's and other values separately; that is, the zero's are from a binomial process only, the other positive values are from a truncated count data process. We assume here, in the simulation, we don't observe $z$. Therefore, $x$ is determining both binary and count processes. In the graph below, it's labeled hurdle. ```{r} #| label: count-sim #| eval: false #| echo: true library(MASS) library(pscl) library(parallel) set.seed(666) gen.sim <- function(df) { nobs <- df["nobs"]; th <- df["th"] z <- rnorm(nobs, 0, 1) x <- rnorm(nobs, 0, 1) u <- rnorm(nobs, 0, 1) y.count <- floor(exp(2 * x + u)) z1 <- 4 * z + th prob <- plogis(z1) y.logit <- rbinom(nobs, size = 1, prob = prob) y <- ifelse(y.logit == 1, y.count, y.logit) m1 <- zeroinfl(y ~ x | z) m4 <- zeroinfl(y ~ x | x) m2 <- glm(y ~ x, family = "poisson") m3 <- lm(log(y + 1) ~ x) m5 <- glm.nb(y ~ x) m6 <- hurdle(y ~ x) c(zip1 = coef(m1, "count")["x"] - 2, poisson = coef(m2)["x"] - 2, log.linear = coef(m3)["x"] - 2, zip2 = coef(m4, "count")["x"] - 2, nb = coef(m5)["x"] - 2, hurdle = coef(m6, "count")["x"] - 2) } data.grid <- expand.grid( nobs = ceiling(exp(seq(4, 9, 1)) / 100) * 100, nsim = seq(1, 100, 1), th = seq(-4, 4, 2) ) cl <- makeCluster(detectCores() - 1) clusterEvalQ(cl, { library(MASS); library(pscl) }) clusterExport(cl, "gen.sim") results <- parApply(cl, data.grid, 1, gen.sim) stopCluster(cl) forshiny <- cbind(data.grid, t(results)) write.csv(forshiny, "results.csv", row.names = FALSE) ``` Count data models can be used even if data is not "counts"; for example, some non-negative non-integer numbers. In fact, Poisson model is consistent even if data is not Poisson-distributed, if the model specification is correct on modeling the log of expected counts. We simulate both scenarios: Case 1, data is generated from a Poisson process. Case 2, data is a floored lognormal -- $\lfloor \exp(2x+u) \rfloor$ with $u$ normal -- so it is non-negative and integer-valued without being a count from any count distribution, and we use count data models on it anyway. (Earlier versions of this text called case 2 "Normal"; the variable being modelled is not normal, only the $u$ inside the exponential is.) The above code is for case 2. We simulate 100 times with $\theta$ ranging from -4 to 4, lower number means higher percentage of zeros; number of observations from $e^4$ to $e^9$ (roughly from 50 to 8000). Since there are many simulations, we use base R's `parallel` library to speed things up. For raw code, please visit [case1: poisson](https://github.com/xiangao/poisson) and [case2: normal](https://github.com/xiangao/poisson2). The code above updates the original `snowfall`-based version to use base R `parallel`. **Status of these simulation findings.** The grid above is retained for its code rather than as current evidence, and the ranking it once carried should not be read as verified. Two things undercut it. The metric for the log-linear row was wrong in the original run -- see the note below -- and the case-2 data-generating process is a floored lognormal whose log-linear index is not exactly $2x$, so subtracting 2 from each model's slope does not compare like with like. The chunk is marked `eval: false`, so nothing here regenerates it. What the run suggested, and what remains plausible on grounds independent of it: plain Poisson QMLE holds up well in large samples even with unmodelled excess zeros *when the conditional mean is correctly specified*; zero-inflated Poisson does better when the zero process is known; and log-linear regression on $\log(y+1)$ targets a different quantity altogether. For evidence rather than recollection, the worked example in the next section is the part of this chapter that actually runs. (Note: the `log.linear` bias line above previously read `exp(coef(m3)["x"]) - 2` instead of `coef(m3)["x"] - 2`, which — since the log-linear coefficient is already on the log scale — manufactured a spurious ≈5.4 apparent bias regardless of the model's actual performance. That has been corrected here to match the other rows. Log-linear regression on `log(y+1)` is still expected to be biased relative to Poisson-family models in a zero-inflated count setting, for real reasons (the log transform's nonlinearity near zero, and the fact that $E[\log(y+1)] \neq \log(E[y]+1)$), but the *magnitude* of "consistently the worst" above reflects the original (buggy) run and has not been re-verified against a corrected simulation.) ## Quick worked example All five models on a single small dataset (n = 1000, ~40% zeros). The DGP is a zero-inflated Poisson built in two independent pieces. We draw $x \sim N(0,1)$ and $z \sim N(0,1)$, independent of each other. The count part is $y_c \sim \text{Poisson}(\exp(1.5 + 0.8x))$ and the zero part is $y_b \sim \text{Bernoulli}(\Lambda(1 + 2z))$, with the observed outcome $y = y_c \cdot y_b$. So the true slope on $x$ is 0.8, and — this is the feature that drives the results — the zero process depends on $z$ alone, never on $x$. ```{r} #| label: quick-example #| cache: true #| warning: false #| message: false library(MASS) library(pscl) set.seed(42) n <- 1000 x <- rnorm(n) z <- rnorm(n) # zero-inflation driver mu <- exp(1.5 + 0.8 * x) # count mean y_c <- rpois(n, lambda = mu) y_b <- rbinom(n, 1, plogis(1 + 2 * z)) # 1 = non-zero y <- y_c * y_b # zero-inflated count cat("Zero fraction:", round(mean(y == 0), 2), "\n") cat("Mean (non-zero):", round(mean(y[y > 0]), 2), "\n\n") # 1. Poisson m_pois <- glm(y ~ x, family = poisson) # 2. Negative binomial m_nb <- glm.nb(y ~ x) # 3. Zero-inflated Poisson (ZIP) — correct zero model m_zip <- zeroinfl(y ~ x | z) # 4. ZIP — wrong zero model (only x) m_zip2 <- zeroinfl(y ~ x | x) # 5. Hurdle m_hurdle <- hurdle(y ~ x) # 6. Log-linear OLS on log(y + 1) m_loglin <- lm(log(y + 1) ~ x) # Compare x coefficient (true = 0.8) coefs <- c( Poisson = coef(m_pois)["x"], NB = coef(m_nb)["x"], `ZIP (correct z)` = coef(m_zip, "count")["x"], `ZIP (wrong)` = coef(m_zip2, "count")["x"], Hurdle = coef(m_hurdle, "count")["x"], `Log-linear log(y+1)` = coef(m_loglin)["x"] ) knitr::kable( data.frame(Model = names(coefs), `x coefficient` = round(coefs, 3), Bias = round(coefs - 0.8, 3)), caption = "True x coefficient = 0.8" ) ``` The sample comes out 39% zeros with a non-zero mean of 6.79. The estimated slopes are 0.850 for Poisson, 0.824 for negative binomial, 0.803 for ZIP with the correct zero model, 0.802 for ZIP with the wrong one, 0.803 for the hurdle model, and 0.458 for log-linear OLS. The first five models recover the true `x` slope of 0.8 reasonably well; the log-linear row is included to make the earlier claim checkable, and needs reading with care --- see the note after the table. This is the key point of *this* DGP: the zero-inflation process depends on `z`, which is generated **independently of `x`**, so omitting the zero model distorts the zero mass and the intercept, not the slope on `x`. The hurdle model and even the ZIP with a *misspecified* (x-only) zero model still recover the slope closely; the plain Poisson shows a small upward bias *in this particular sample* (n = 1000), but that is finite-sample noise, not a systematic property of the model. Since the zero-inflation driver `z` is independent of `x`, the conditional mean here is correctly specified for Poisson QMLE — increasing `n` from 1,000 to 1,000,000 shrinks the Poisson `x`-coefficient bias from about +0.04 to effectively 0, confirming consistency. The unmodeled excess zeros create overdispersion, but overdispersion under a correctly-specified conditional mean affects only the standard errors (which should be corrected with a robust/sandwich estimator), not the consistency of the slope. About the log-linear row: its coefficient is **not** estimating the same thing as the others, which is the real reason to avoid it rather than any particular bias number. The count models target $\log E[y \mid x]$, whereas OLS on $\log(y+1)$ targets $E[\log(y+1) \mid x]$, and $E[\log(y+1)] \neq \log(E[y]+1)$ by Jensen's inequality. The $+1$ offset makes the discrepancy worst exactly where the mass of zeros is. So comparing its coefficient to 0.8 is an apples-to-oranges comparison, and the honest statement is that log-linear answers a different question --- not that it answers this one badly by some measured amount. A misspecified zero process biases the count *slope* only when the zero process is correlated with the count covariates — which is not the case here. (Change `z` to depend on `x` and the omitted-zero-process slope bias appears.) ## Mixed count models with glmmTMB When data has both overdispersion *and* random effects (clustered counts, panel data, repeated measures), `glmmTMB` handles the full range in one package: ```{r} #| eval: false #| echo: true library(glmmTMB) # Zero-inflated negative binomial with random intercept per subject m_glmmtmb <- glmmTMB( count ~ treatment + time + (1 | subject), # count model ziformula = ~ treatment, # zero-inflation model family = nbinom2, data = your_data ) summary(m_glmmtmb) # Conditional vs marginal predictions predict(m_glmmtmb, type = "conditional") # excluding zero-inflation predict(m_glmmtmb, type = "response") # including zero-inflation ``` `glmmTMB` supports: Poisson, NB1, NB2, zero-inflated variants of all of the above, hurdle models, beta-binomial, tweedie, and more. The `ziformula` argument takes a one-sided formula for the zero-inflation component — use `~ 1` for a constant zero-inflation rate, or covariates to model what drives the structural zeros.