A more interpretable regression model for count data with excess of zeros

Abstract

Count data are common in medical research. When these data have more zeros than expected by the most used count distributions, it is common to employ a zero-inflated regression model. However, the interpretability of these models is much lower than the most used count regression models. We present a more interpretable regression model that estimates the mean event rate and models covariate-dependent dispersion directly. Additionally, the dispersion parameter can be interpreted as an index of clumping, a control parameter for overdispersion. We discuss inferential and diagnostic tools and perform a Monte Carlo simulation study to evaluate the performance of the maximum likelihood estimator. Finally, the usefulness of the proposed regression model is illustrated through an application on antenatal care visits.

Keywords

Count data quantile residual zero-inflated data zero-inflated Poisson regression model

1. Introduction

Count data are common in several areas of medicine such as immunology,¹ cardiology,² urology,³ pediatrics,⁴ psychiatry⁵ and neurology.⁶ A possible way to study the relationship between these variables and a set of explanatory variables is assuming that the response has a certain count probability distribution such as Poisson, negative binomial, Poisson inverse Gaussian, among many others. However, it is not unusual that the response variable has more zeros than expected by the best known count distributions.

When there is an excess of zeros, the response variable is usually modeled using a zero-inflated regression model. These models assume that the response variable is a mixture of a degenerated distribution at zero and a known count probability distribution. Historically, several zero-inflated regression models were proposed. For example, the zero-inflated Poisson (ZIP) regression model was proposed by Lambert.⁷ Subsequently, the zero-inflated negative binomial (ZINB) regression model was proposed by Ridout et al.⁸ and the zero-inflated logarithmic (ZIL) regression model was introduced by Rigby,⁹ among others.

In the zero-inflated regression models, there are two means related to the response variable. The first is the mean of the response variable. The other is the mean of the count probability distribution that composes the distribution of the response variable. We usually want to model the former mean as a function of covariates. However, in almost all of the zero-inflated regression models, the latter mean is modeled. The reason is that the vector of parameters of the distribution of the response variable in these models includes the mean of the count probability distribution and not the mean of the response variable.

To solve this issue, Long et al.¹⁰ proposed a reparameterization of the ZIP distribution that they called marginalized ZIP (MZIP). In the MZIP parameterization, one of the parameters is the mean of response variable (marginalized mean) and so it is more convenient to use in a regression model that they introduced in the same work. After this work, marginalized zero-inflated models based on negative binomial distribution (MZINB)¹¹ and generalized Poisson distribution (MZIGP)¹² were proposed. Todem et al.¹³ and Martin and Hall¹⁴ introduced more general marginalized forms of models for zero-inflated count data. Extensions of MZIP and MZINB to models with random effects were proposed by Long et al.¹⁵ and Burgette et al.,¹⁶ respectively, while an extension to MZINB to spatial data was proposed by Mutiso et al.¹⁷ From a software standpoint,¹⁸ present new commands in Stata for modeling count data using marginalized Zero-inflated distributions. As an applied illustration,¹⁹ investigate the association between sugar sweetened beverage intake and dental caries among underserved Black adolescents utilizing a marginalized ZIP (MZIP) regression analysis. Finally, in Mun et al.²⁰ the MZIP and MZINB models are used in a two-step individual participant data (IPD) meta-analysis.

To the best of our knowledge, there is only one zero-inflated regression model based on the Poisson distribution in which the mean of the response variable is directly modeled. As mentioned before, this model was proposed by Long et al.¹⁰ It can be fitted in the gamlss²¹ and in the mcount package²² of the software R and was considered, for example, by Seidel et al.,²³ Laurence et al.,¹⁹ Mun et al.,²⁰ and Sims et al.,²⁴ However, this parameterization does not change the other parameters of the distribution. As a result, it is not parameterized in terms of a readily interpretable dispersion parameter. The same issue is observed in the other marginalized zero-inflated regression models.

This work proposes a zero-inflated regression model based on a novel reparameterization of the zero-inflated Poisson distribution. This parameterization is more useful than the existing ones because one of the parameters of the distribution is the mean, and the other is a dispersion parameter. As a result, the proposed model is more interpretable than other zero-inflated regression models. We assume that both parameters of the distribution of the response are functions of covariates, and so the proposed model has a structure similar to a double generalized linear model.²⁵ As we assume that the mean and a dispersion parameter of the response are functions of covariates, our model can also be seen as a link between marginalized models and dispersion models.²⁶ Some contributions to dispersion models on count data include Bonat et al.²⁷ and Bar-Lev and Ridder.²⁸

In the same way as Long et al.,¹⁰ our model formulation offers meaningful statements about an exposure effect on an entire population in contrast to the traditional ZIP model whose regression parameters have interpretations for unobservable latent classes.¹⁰ However, our model not only estimates the mean intensity of events (as in the traditional Poisson regression) but also captures how the variability around this mean depends on observable factors.

The remainder of this paper is organized as follows. Section 2 proposes a reparameterization of the zero-inflated Poisson distribution. Section 3 introduces the regression model associated with this novel reparameterization. Diagnostic tools for this regression model are discussed in Section 4. In Section 5, Monte Carlo simulation studies are performed to evaluate the performance of the maximum likelihood estimators of the parameters of the proposed regression model. The usefulness of our model is illustrated through an application on antenatal care visits presented in Section 6. Concluding remarks are provided in Section 7.

2. Parameterizations of the ZIP

The first version of the ZIP distribution was introduced by Lambert.⁷ In this setup, the probability mass function is given by:

Pr (Y = y | λ, p) = {\begin{cases} p + (1 - p) e^{- λ}, & if y = 0, \\ (1 - p) e^{- λ} λ^{y} / y!, & if y = 1, 2, 3, \dots, \end{cases}

(1)

where

λ > 0

0 < p < 1

. In this parameterization,

E (Y) = (1 - p) λ

and

Var (Y) = λ (1 - p) (1 + λ p)

. This probability mass function is derived based on the assumption that data are composed of two (unobservable) subpopulations. The response variable is zero for the first subpopulation and is Poisson distributed for the other. The parameter

p

is the probability of belonging to the first subpopulation and

λ

is the mean of the response variable for the second subpopulation. However, since the two subpopulations are not observable, the parameters

p

and

λ

are of less interest in practice.

A regression model is more easily interpretable when one of the parameters of the distribution of the response variable is the mean or the median. Long et al.¹⁰ proposed a new parameterization of the ZIP distribution, in which one of the parameters is the mean. From the ZIP parameterization, they considered that $μ = (1 - p) λ$ , that is, $λ = μ / (1 - p)$ . Therefore they obtain from (1) the following probability mass function

Pr (Y = y | μ, p) = {\begin{cases} p + (1 - p) e^{- (μ / (1 - p))}, & if y = 0, \\ \frac{(μ)^{y}}{y! (1 - p)^{y - 1}} e^{- (μ / (1 - p))}, & if y = 1, 2, 3, \dots, \end{cases}

where

μ > 0

0 < p < 1

. In this parameterization,

E (Y) = μ

and

Var (Y) = μ [1 + μ p / (1 - p)]

. We will refer to this parameterization as ZIP2, because this is the name used in the gamlss package. Note that from the ZIP to the ZIP2 parameterization, the parameter

p

was not changed. As a result,

p

is not a readily interpretable dispersion parameter.

Here we propose a novel parameterization of the ZIP distribution, in which both parameters are of direct interest. Moreover, in this parameterization, the variance of the ZIP distribution is a simple function of the parameters. The ZIP distribution in our proposed parameterization is indexed by the mean and a dispersion parameter. From the ZIP parameterization, we consider that $μ = (1 - p) λ$ and $ϕ = p λ$ , that is, $p = ϕ / (μ + ϕ)$ and $λ = μ + ϕ$ . Therefore, we obtain from (1) the following probability mass function

Pr (Y = y ∣ μ, ϕ) = {\begin{cases} \frac{ϕ + μ e^{- (μ + ϕ)}}{μ + ϕ}, & if y = 0, \\ \frac{μ (μ + ϕ)^{y - 1} e^{- (μ + ϕ)}}{y!}, & if y = 1, 2, 3, \dots, \end{cases}

(2)

where

μ > 0

ϕ > 0

. In this parameterization,

E (Y) = μ

and

Var (Y) = μ (1 + ϕ)

and therefore the variance is a linear function of

ϕ

. We will refer to this parameterization as ZIP3. In the ZIP3 parameterization, the mean and the variance are similar to those of the parameterization of the negative binomial known as NBII in the gamlss package.⁹ However, the probability functions of ZIP3 and NBII are considerably distinct and the fit of their respective regression model to specific data can be very different as seen in Section 6.

Note that in the ZIP3 parameterization, $ϕ$ is the index of clumping (IOC). This index was proposed by David and Moore²⁹ and is defined as $IOC (Y) = Var (Y) / E (Y) - 1$ . Using the expressions for $E (Y)$ and $Var (Y)$ presented above, we easily obtain that $ϕ = IOC (Y)$ in the ZIP3 parameterization. Observe that the IOC corresponds to the relative increase in the variance compared to the variance of Poisson distribution. Therefore, in the ZIP3 parameterization, if $ϕ = l$ , then the variance is $(l \times 100) %$ higher than the mean. The IOC has been widely used in ecology (see for example Hasan et al.³⁰ and Zhang et al.³¹) and also in works in the area of statistics.^32,33

Note that ZIP, ZIP2 and ZIP3 are different parameterizations of the same probability distribution (zero-inflated Poisson distribution). As a result, any parameter of the three parameterizations can be linked with simple formula to the parameters of the other two. While the fits of the regression models based on these three parameterizations are often numerically similar in practice, the models (and their fits) are distinct except in special cases, that is, intercept only models and saturated models based solely on categorical explanatory variables. The primary appeal of the ZIP3 parameterization is in the interpretation of the regression parameters.

3. A regression model

Considering the interpretability advantages of ZIP3 parameterization, it is very convenient to use when the response variable has a high proportion of zeros. We define in this section a regression model based on the ZIP3 parameterization and use it to fit real data in Section 6. However, before introducing our regression model, we obtain some results that enable us to use the model in an useful computational framework.

The expression (2) can be rewritten as follows:

Pr (Y = y ∣ μ, ϕ) = {[\frac{ϕ + μ e^{- (μ + ϕ)}}{μ + ϕ}]}^{I (y = 0)} {[\frac{μ (μ + ϕ)^{y - 1} e^{- (μ + ϕ)}}{y!}]}^{(1 - I (y = 0))}, y \in Z_{0}^{+},

(3)

where

I (\cdot)

is an indicator function, that is,

I (y = 0) = 1

y = 0

and 0 otherwise. Equation (3) enables us to obtain the logarithm of the probability function. The expression can be stated as follows:

\begin{aligned} ℓ (μ, ϕ ∣ y) & = I (y = 0) [\log (ϕ + μ e^{- (μ + ϕ)}) - \log (μ + ϕ)] + (1 - I (y = 0)) [\log (μ) + (y - 1) \times \\ \log (μ + ϕ) - μ - ϕ - \log (y!)] . \end{aligned}

(4)

Additionally, the partial derivatives of first and second order of Equation (4) with respect to the parameters indexing the ZIP3 distribution are given by:

\begin{aligned} \partial_{μ} ℓ (μ, ϕ ∣ y) & = I (y = 0) [\frac{(1 - μ)}{(μ + ϕ e^{μ + ϕ})} - \frac{1}{(μ + ϕ)}] + (1 - I (y = 0)) [\frac{1}{μ} + \frac{(y - 1)}{(μ + ϕ)} - 1], \\ \partial_{μ μ}^{2} ℓ (μ, ϕ ∣ y) & = I (y = 0) [\frac{ϕ (μ - 2) e^{μ + ϕ} - 1}{(μ + ϕ e^{μ + ϕ})^{2}} + \frac{1}{(μ + ϕ)^{2}}] + (1 - I (y = 0)) [- \frac{1}{μ^{2}} - \frac{(y - 1)}{(μ + ϕ)^{2}}], \\ \partial_{ϕ} ℓ (μ, ϕ ∣ y) & = I (y = 0) [\frac{e^{μ + ϕ} - μ}{(μ + ϕ e^{μ + ϕ})} - \frac{1}{(μ + ϕ)}] + (1 - I (y = 0)) [\frac{(y - 1)}{(μ + ϕ)} - 1], \\ \partial_{ϕ ϕ}^{2} ℓ (μ, ϕ ∣ y) & = I (y = 0) [\frac{e^{μ + ϕ} [μ (ϕ + 2) - e^{μ + ϕ}]}{(μ + ϕ e^{μ + ϕ})^{2}} + \frac{1}{(μ + ϕ)^{2}}] + (1 - I (y = 0)) [- \frac{(y - 1)}{(μ + ϕ)^{2}}], \\ \partial_{μ ϕ}^{2} ℓ (μ, ϕ ∣ y) & = I (y = 0) [\frac{(μ - 1) (ϕ + 1) e^{μ + ϕ}}{(μ + ϕ e^{μ + ϕ})^{2}} + \frac{1}{(μ + ϕ)^{2}}] + (1 - I (y = 0)) [- \frac{(y - 1)}{(μ + ϕ)^{2}}] . \end{aligned}

(5)

3.1. The ZIP3 regression model

In ZIP3 regression, the response variables $Y = (Y_{1}, \dots, Y_{n})^{⊤}$ are independent and follow a ZIP distribution with parameters $μ_{i}$ and $ϕ_{i}$ as defined in the ZIP3 parameterization presented in (2). Moreover, the parameters $μ = (μ_{1}, \dots, μ_{i}, \dots, μ_{n})^{⊤}$ and $ϕ = (ϕ_{1}, \dots, ϕ_{i}, \dots, ϕ_{n})^{⊤}$ satisfy

g_{1} (μ_{i}) = x_{i}^{⊤} β = η_{i}, and g_{2} (ϕ_{i}) = z_{i}^{⊤} γ = ς_{i}, (i = 1, \dots, n),

(6)

for (vectors of) covariates

x_{i}^{⊤} = (1, x_{i 2}, \dots, x_{i q_{1}})

, and

z_{i}^{⊤} = (1, z_{i 2}, \dots, z_{i q_{2}})

with

β = (β_{1}, \dots, β_{q_{1}})^{⊤}

and

γ = (γ_{1}, \dots, γ_{q_{2}})^{⊤}

being the parameter vectors associated with

x_{i}

and

z_{i}

, respectively. Additionally,

g_{k} (\cdot), (k = 1, 2)

specifies the link between the random and the systematic components and is a strictly monotonic and twice differentiable function. It follows that the ZIP3 regression log-likelihood is:

\begin{aligned} ℓ_{1} (μ, ϕ ∣ y) & = \sum_{i = 1}^{n} ℓ (μ_{i}, ϕ_{i} ∣ y_{i}) \\ = \sum_{i = 1}^{n} ϱ_{i} [\log (ϕ_{i} + μ_{i} e^{- (μ_{i} + ϕ_{i})}) - \log (μ_{i} + ϕ_{i})] + (1 - ϱ_{i}) [\log (μ_{i}) + (y_{i} - 1) \times \\ \log (μ_{i} + ϕ_{i}) - μ_{i} - ϕ_{i} - \log (y_{i}!)] . \end{aligned}

(7)

where

ϱ_{i} = I (y_{i} = 0)

μ_{i} = g_{1}^{- 1} (x_{i}^{⊤} β)

, and

ϕ_{i} = g_{2}^{- 1} (z_{i}^{⊤} γ)

Let $θ = (β^{⊤}, γ^{⊤})^{⊤}$ be the unknown $s$ -dimensional $(s := q_{1} + q_{2})$ parameter in models (6). The maximum likelihood (ML) estimator $\hat{θ} = ({\hat{β}}^{⊤}, {\hat{γ}}^{⊤})^{⊤}$ of $θ$ is the solution of the $s$ -dimensional score equation

U (\hat{θ}) = 0,

(8)

where

U (θ) = (U_{β} (θ)^{⊤}, U_{γ} (θ)^{⊤})^{⊤}

is the score vector. The score function is given by taking the first derivative of the log-likelihood function, (7), with respect to each element of

θ

. By the chain rule, it follows that

U_{β_{j}} (θ) = \sum_{i = 1}^{n} \frac{\partial ℓ (μ_{i}, ϕ_{i} ∣ y_{i})}{\partial μ_{i}} \frac{\partial μ_{i}}{\partial η_{i}} \frac{\partial η_{i}}{\partial β_{j}} = \sum_{i = 1}^{n} d_{μ_{i}} l_{μ_{i}} x_{i j}, j = 1, \dots, q_{1},

and

U_{γ_{t}} (θ) = \sum_{i = 1}^{n} \frac{\partial ℓ (μ_{i}, ϕ_{i} ∣ y_{i})}{\partial ϕ_{i}} \frac{\partial ϕ_{i}}{\partial ς_{i}} \frac{\partial ς_{i}}{\partial γ_{t}} = \sum_{i = 1}^{n} d_{ϕ_{i}} l_{ϕ_{i}} z_{i t}, t = 1, \dots, q_{2},

where

\begin{aligned} d_{μ_{i}} & = ϱ_{i} [\frac{(1 - μ_{i})}{(μ_{i} + ϕ_{i} e^{μ_{i} + ϕ_{i}})} - \frac{1}{(μ_{i} + ϕ_{i})}] + (1 - ϱ_{i}) [\frac{1}{μ_{i}} + \frac{(y_{i} - 1)}{(μ_{i} + ϕ_{i})} - 1]; \\ l_{μ_{i}} & = \frac{1}{g_{1}^{'} (μ_{i})}; \\ d_{ϕ_{i}} & = ϱ_{i} [\frac{e^{μ_{i} + ϕ_{i}} - μ_{i}}{(μ_{i} + ϕ_{i} e^{μ_{i} + ϕ_{i}})} - \frac{1}{(μ_{i} + ϕ_{i})}] + (1 - ϱ_{i}) [\frac{(y_{i} - 1)}{(μ_{i} + ϕ_{i})} - 1]; \\ l_{ϕ_{i}} & = \frac{1}{g_{2}^{'} (ϕ_{i})}, \end{aligned}

and the score vector can be written compactly as

U_{β} (θ) = X^{⊤} L_{μ} [(y * - μ *) + ϱ ⊙ c_{1}] and U_{γ} (θ) = Z^{⊤} L_{ϕ} [(y * - 1_{n}) + ϱ ⊙ c_{2}],

where

⊙

represents the Hadamard product,³⁴

X

is a

n \times q_{1}

matrix with the

i

th row given by

x_{i}^{⊤}

Z

is a

n \times q_{2}

matrix with the

i

th row given by

z_{i}^{⊤}

L_{μ} = diag (l_{μ_{1}}, \dots, l_{μ_{n}})

L_{ϕ} = diag (l_{ϕ_{1}}, \dots, l_{ϕ_{n}})

{y *}^{⊤} = (\frac{y_{1} - 1}{μ_{1} + ϕ_{1}}, \dots, \frac{y_{n} - 1}{μ_{n} + ϕ_{n}})

1_{n}^{⊤} = (1, \dots, 1)

μ * = (1 - \frac{1}{μ_{1}}, \dots, 1 - \frac{1}{μ_{n}})

c_{1}^{⊤} = (\frac{1 - μ_{1}}{μ_{1} + ϕ_{1} e^{μ_{1} + ϕ_{1}}} + \frac{μ_{1} - 1}{μ_{1}} - \frac{y_{1}}{μ_{1} + ϕ_{1}}), \dots, \frac{1 - μ_{n}}{μ_{n} + ϕ_{n} e^{μ_{n} + ϕ_{n}}} + \frac{μ_{n} - 1}{μ_{n}} - \frac{y_{n}}{μ_{n} + ϕ_{n}})

c_{2}^{⊤} = (\frac{e^{μ_{1} + ϕ_{1}} - μ_{1}}{μ_{1} + ϕ_{1} e^{μ_{1} + ϕ_{1}}} + \frac{μ_{1} + ϕ_{1} - y_{1}}{μ_{1} + ϕ_{1}}, \dots, \frac{e^{μ_{n} + ϕ_{n}} - μ_{n}}{μ_{n} + ϕ_{n} e^{μ_{n} + ϕ_{n}}} + \frac{μ_{n} + ϕ_{n} - y_{n}}{μ_{n} + ϕ_{n}})

, and

ϱ^{⊤} = (ϱ_{1}, \dots, ϱ_{n})

With the results presented in Equation (5), we can readily integrate the ZIP3 distribution into the distribution family framework of the R package gamlss.³⁵ To make this integration, we developed a suite of functions, encompassing the pseudo-random number generator, the probability density function, the quantile function, and the cumulative distribution function. All codes are available on Github through the link https://github.com/statlab-oficial/ZIP3. With this integration, one can fit the ZIP3 regression model taking advantage of all inferential and diagnostic tools of the gamlss package. We intend to include the ZIP3 regression model in the gamlss package after this work is published.

4. Diagnostic analysis

We propose using a residual and a global influence measure for model diagnostics for ZIP3 regression model.

4.1. Residual analysis

When the response of a regression model is discrete, Pearson and deviance residuals are also discrete. As a result, these residuals have a considerable probability of not detecting lack of fit.³⁶ For this reason, when the response is discrete, it is better to use the randomized quantile residual³⁷ to evaluate the goodness of fit of the regression model. For the ZIP3 regression model, the randomized quantile residual is given by

q_{i} = Φ^{- 1} (u_{i}),

(9)

where

u_{i}

is a uniform random variable on the interval

(F (y_{i} - 1; {\hat{μ}}_{i}, {\hat{ϕ}}_{i}), F (y_{i}; {\hat{μ}}_{i}, {\hat{ϕ}}_{i}))

Φ (\cdot)

and

F (\cdot)

are the cumulative distribution function of the standard normal distribution and of the ZIP3 distribution, respectively, and

{\hat{μ}}_{i}

and

{\hat{ϕ}}_{i}

are the ML estimates of the parameters

μ_{i}

and

ϕ_{i}

, respectively. The residual

q_{i}

is asymptotically standard normally distributed under the correct model

4.2. Global influence

According to Cook et al.,³⁸ the likelihood displacement^{39, p. 182} is the most useful measure for identifying influential observations. It has a similar expression for all parametric regression models and has been widely used in recent works.^40–42 For the ZIP3 regression model, the likelihood displacement is given by

{LD}_{i} = 2 [ℓ_{1} (\hat{μ}, \hat{ϕ} ∣ y) - ℓ_{1} ({\hat{μ}}_{(i)}, {\hat{ϕ}}_{(i)} ∣ y)],

(10)

where

{\hat{μ}}_{(i)}

and

{\hat{ϕ}}_{(i)}

denote the ML estimators of

μ

and

ϕ

, respectively, obtained by dropping the

i

th observation. The calculation of

{LD}_{i}

for the

n

observations requires the estimation of

(n + 1)

ZIP3 regression models. However, this is not an issue, since the fit of a ZIP3 regression model is fast when using the code developed in this work.

5. Simulation studies

We conducted Monte Carlo (MC) simulation studies to evaluate the performance of the ML estimators of the parameters of the ZIP3 regression model using small and moderate sample sizes. The scenario under analysis considers the following: sample sizes $n \in {50, 100, 200, 500}$ and values for the parameters as presented in (11),

\log (μ_{i}) = - 1.0 + 1.0 x_{i 1} + 0.5 x_{i 2} and \log (ϕ_{i}) = 1.0 + 0.5 z_{i 1},

(11)

where the covariates

x_{i 1}

and

z_{i 1}

, for

i = 1, \dots, n

, were generated from the standard uniform distribution, and

x_{i 2}

was generated from the uniform distribution on an interval

(1, 3)

. The number of MC replications was 5,000 and all simulations were performed using the R programming language.

For each value of the parameter and sample size, we report the bias (B) and mean squared error (MSE) of the ML estimators in Table 1. Note that, as the sample size increases, the bias and mean squared error of the ML estimators decrease, as expected. The biases are small, except for ${\hat{β}}_{0}$ , for which the bias is small only for $n \geq 100$ . The mean squared errors are moderate when $n = 50$ , but small for all ML estimators when $n = 500$ .

Table 1.

Bias and MSE of the ML estimator in considered scenario.

	Mean parameter						Precision parameter
$n$	$B ({\hat{β}}_{0})$	$B ({\hat{β}}_{1})$	$B ({\hat{β}}_{2})$	$MSE ({\hat{β}}_{0})$	$MSE ({\hat{β}}_{1})$	$MSE ({\hat{β}}_{2})$	$B ({\hat{γ}}_{0})$	$B ({\hat{γ}}_{1})$	$MSE ({\hat{γ}}_{0})$	$MSE ({\hat{γ}}_{1})$
50	−0.0918	0.0092	0.0050	1.3104	0.6522	0.2002	−0.0372	−0.0046	0.1338	0.3250
100	−0.0499	−0.0006	0.0047	0.5922	0.3544	0.0887	−0.0209	0.0011	0.0590	0.1674
200	−0.0289	0.0009	0.0033	0.2623	0.1631	0.0391	−0.0059	−0.0032	0.0286	0.0763
500	−0.0137	0.0022	0.0012	0.0981	0.0601	0.0153	−0.0017	−0.0026	0.0110	0.0296

To assess the calibration of the empirical distribution of the fitted model, we compare the overall observed mean with the predicted mean for the location component, and verify that the empirical variances align with the model-implied variances for the dispersion component. In addition, we evaluated the average proportion of zeros in the generated sample and compared it with the average proportion of zeros obtained from the fitted model, where we can observe that the values are very close. These results are presented in Table 2, where we can observe that the fitted model is well calibrated.

Table 2.

Distributional calibration checks and observed and predicted zero proportion in considered scenario.

	Distributional calibration				Proportion of zeros
$n$	$\bar{Y}$	$\hat{E (Y)}$	$S_{Y}^{2}$	$\hat{Var (Y)}$	$\bar{p}$	$\hat{p}$
50	1.8931	1.8947	9.4995	9.7338	0.6797	0.6859
100	1.8147	1.8136	8.8941	8.9655	0.6825	0.6850
200	1.7921	1.7913	8.8763	8.9152	0.6880	0.6895
500	1.7861	1.7855	8.8889	8.8967	0.6859	0.6895

In Table 3 we observe that the sample standard deviations (SD) of the estimates are close to the estimated asymptotic standard errors (SE) and they are closer as $n$ increases.

Table 3.

Standard errors and standard deviation estimates.

	Mean parameter						Precision parameter
$n$	$SD ({\hat{β}}_{0})$	$SD ({\hat{β}}_{1})$	$SD ({\hat{β}}_{2})$	$SE ({\hat{β}}_{0})$	$SE ({\hat{β}}_{1})$	$SE ({\hat{β}}_{2})$	$SD ({\hat{γ}}_{0})$	$SD ({\hat{γ}}_{1})$	$SE ({\hat{γ}}_{0})$	$SE ({\hat{γ}}_{1})$
50	1.1410	0.8075	0.4473	1.0598	0.7568	0.4162	0.3638	0.5700	0.3539	0.5596
100	0.7678	0.5953	0.2978	0.7404	0.5780	0.2871	0.2422	0.4091	0.2425	0.4085
200	0.5113	0.4039	0.1977	0.4971	0.3959	0.1939	0.1691	0.2762	0.1675	0.2752
500	0.3129	0.2450	0.1238	0.3059	0.2471	0.1214	0.1048	0.1719	0.1043	0.1709

To estimate the parameters by maximum likelihood, we used the RS algorithm.³⁵ Table 4 shows that the RS algorithm exhibits highly stable performance, requiring around 3.6 to 3.8 iterations across all sample sizes. Besides that, according to the results presented in Table 4, it can be observed that, among the 5000 Monte Carlo replications, non-convergence occurred in only two cases when the sample size is equal to 50.

Table 4.

Iterations and non-convergences.

$n$	Average number of iterations	Non-convergences
50	3.81	2
100	3.65	0
200	3.63	0
500	3.67	0

Finally, for $n \geq 100$ convergence is attained in all replications, demonstrating numerical stability. Therefore, the RS algorithm is efficient and reliable.

Table 5 presents the empirical coverage probabilities of the asymptotic confidence intervals of level $100 (1 - α) %$ , denoted by $CI (\cdot; 100 (1 - α) %)$ . Note that, for the mean parameter, the asymptotic confidence intervals have an empirical coverage probability slightly lower than the nominal levels in small samples. However, the empirical coverages are very close to the nominal levels when $n = 500$ . For the precision parameter, the empirical coverages are close to the nominal levels, even in small samples.

Table 5.

Coverage probabilities of asymptotic confidence intervals.

	Mean parameter						Precision parameter
$n$	$CI (β_{0}; 99 %)$	$CI (β_{0}; 95 %)$	$CI (β_{1}; 99 %)$	$CI (β_{1}; 95 %)$	$CI (β_{2}; 99 %)$	$CI (β_{2}; 95 %)$	$CI (γ_{0}; 99 %)$	$CI (γ_{0}; 95 %)$	$CI (γ_{1}; 99 %)$	$CI (γ_{1}; 95 %)$
50	0.983	0.938	0.988	0.941	0.987	0.945	0.990	0.955	0.990	0.949
100	0.989	0.942	0.988	0.949	0.988	0.945	0.990	0.946	0.991	0.951
200	0.989	0.945	0.991	0.949	0.990	0.947	0.993	0.956	0.992	0.952
500	0.989	0.945	0.992	0.953	0.989	0.947	0.993	0.949	0.993	0.952

We also considered two other scenarios, considering the following vectors for $β$ and $γ$ : Case 2: $β = (0.5, 0.5, 1.0)^{⊤}$ and $γ = (1.0, 0.5)^{⊤}$ ; Case 3: $β = (- 1, 0.5, 1.0)^{⊤}$ and $γ = (- 0.5, 0.5)^{⊤}$ . However, due to space limitations, these results are not included. The findings for these additional scenarios are nevertheless consistent with those presented in the text.

6. Application to data on number of antenatal care visits in Oromia—Ethiopia

Antenatal care (ANC) is the service delivered by skilled healthcare providers to pregnant women and adolescent girls to ensure the best health conditions for both mother and baby during pregnancy.⁴³ ANC reduces maternal and newborn mortalities and makes it easier to detect infections early and prevent them from progressing.⁴⁴ The World Health Organization (WHO) has recommended eight ANC visits during the course of pregnancy since 2016.⁴⁵

Women in sub-Saharan Africa have a significantly higher risk of dying from causes related to pregnancy and childbirth than in any other region, representing for 70% of all maternal deaths in 2020 globally.⁴⁶ To increase the ANC visits and to reduce this mortality, it is valuable to identify covariates that are related to the number of ANC visits. Here, we consider data on the number of antenatal care visits in the region of Oromia, Ethiopia. Data refers to the year of 2019 and was obtained from the Demographic and Health Surveys (DHS) Program.⁴⁷ Data have information on 491 mothers who had a baby in the previous 5 years and the response variable is the number of antenatal care visits in the last pregnancy of the women. Note that, for this response variable, it is desirable that the mean of the response is high (8 or even greater). However, dispersion is also important. Dispersion should be low to ensure that few mothers have low ANC visits when the mean is high. The following covariates are available: age, type of residence, education level, birth order number, wealth index and number of household members,

Figure 1 presents a bar plot of the response variable. The number of antenatal care visits by mother in Oromia has an asymmetric distribution and the sample has many zeros (29.5%). Unfortunately, few mothers in the sample had the recommended 8 ANC visits. The high proportion of zeros, coupled with the fact that the number of mothers who had 3 or 4 ANC visits is greater than the number who had 1 or 2 ANC visits, suggests that a zero-inflated regression model may be needed to adequately fit the response variable. In addition, the ZIP3 regression model may be a good choice to model the ANC visits in Oromia, since dispersion is important for this study.

Figure 1.

Bar plot for the number of antenatal care visits in Oromia.

The ZIP3 regression model was fitted considering a logarithmic link function for $μ$ and for $ϕ$ . We selected covariates for the marginal mean submodel and for $ϕ$ especially based on the results of likelihood ratio tests. Table 6 presents the parameter estimates, standard errors and $p$ -values of the likelihood ratio tests for the final ZIP3 regression model. Note that the estimates of the parameters associated with educational level and wealth index for the submodel for $μ$ are positive. Therefore, as expected, the mean of antenatal care visits is higher for women who have a higher education and are from a family with a higher income. In addition, the estimates of the parameters associated with wealth index for the submodel for $ϕ$ are negative. Hence, the index of clumping of the antenatal care visits is lower for richer women.

Table 6.

The final ZIP3 model for the antenatal care visits in Oromia—Ethiopia.

Submodel	Covariate	Category	Estimate	Std. Error	$p$ -value
$μ$	Intercept		0.521	0.123	$< 0.0001$
	Educational	No education (ref)
	level	Primary	0.208	0.080
		Secondary	0.243	0.117	0.0473
		Higher	0.250	0.232
	Wealth	Poorest (ref)
	index	Poorer	0.180	0.147
		Middle	0.419	0.139
		Richer	0.432	0.144	0.0073
		Richest	0.649	0.146
$ϕ$	Intercept		0.398	0.149	0.0080
	Wealth	Poorest (ref)
	index	Poorer	−0.233	0.211
		Middle	−0.667	0.260	0.0007
		Richer	−0.763	0.292
		Richest	−1.462	0.507

For a better interpretation of the results of the final ZIP3 regression model, the exponential of the parameter estimates were calculated (ratio effects). We also calculated the $95 %$ confidence intervals for these ratio effects and the results are presented in Table 7. For example, it is estimated that the mean of antenatal care visits by mother is $91.3 % (CI (95 %) = (43.7, 154.5))$ higher for mothers from the richest families than for those from the poorest families considering mother with the same education level. It is also estimated that the IOC antenatal care visits is $76.8 %$ lower for mothers from the richest families than for those from the poorest families. These findings suggest that access to antenatal care visits in Oromia, Ethiopia, is much lower for poor women than for richer ones.

Table 7.

Confidence intervals (95%) for the ratio effects (exponential of the parameters) for the final model for the antenatal care visits in Oromia—Ethiopia.

			Ratio	Confidence interval
Submodel	Covariate	Category	effect	Lower limit	Upper limit
$μ$	Intercept		1.683	1.323	2.141
	Educational	No education (ref)
	level	Primary	1.232	1.053	1.441
		Secondary	1.275	1.013	1.605
		Higher	1.284	0.815	2.021
	Wealth	Poorest (ref)
	index	Poorer	1.198	0.899	1.596
		Middle	1.520	1.157	1.998
		Richer	1.540	1.162	2.041
		Richest	1.913	1.437	2.545
$ϕ$	Intercept		1.489	1.111	1.995
	Wealth	Poorest (ref)
	index	Poorer	0.793	0.524	1.199
		Middle	0.513	0.308	0.855
		Richer	0.466	0.263	0.827
		Richest	0.232	0.086	0.626

We used the tools discussed in Section 4 to conduct the diagnostic analysis in the final ZIP3 regression model. The left plot of Figure 2 presents a normal probability plot with a simulated envelope⁴⁸ using the randomized quantile residual. The plot does not suggest model misspecification. We also obtained the likelihood displacement for the 491 observations (right plot of Figure 2). Observation {94} has considerably higher value of the likelihood displacement. This woman was in the higher educational level and had 6 antenatal care visits that is more than expected, considering that she was not in the richest families. To study the impact on model inference after removing this case identified as potentially influential, we fitted the model without this observation.

Figure 2.

Normal probability plot with simulated envelope and index plot of the likelihood displacement for the Oromia—Ethiopia data. (a) Envelope and (b) Likelihood displacement

Table 8 presents the relative changes (RC) in the parameter estimates and their corresponding changes in the estimated standard errors (RCSE), based on antenatal care visit data. These changes are calculated from

RC ({\hat{θ}}_{j})_{(i)} = | \frac{{\hat{θ}}_{j} - {\hat{θ}}_{j (i)}}{{\hat{θ}}_{j}} | \times 100 % and RCSE ({\hat{θ}}_{j})_{(i)} = | \frac{SE ({\hat{θ}}_{j}) - SE ({\hat{θ}}_{j})_{(i)}}{SE ({\hat{θ}}_{j})} | \times 100 %,

where

{\hat{θ}}_{j (i)}

and

SE ({\hat{θ}}_{j})_{(i)}

represent the ML estimates of the

j

th parameter of the model and the estimates of the standard error of the corresponding estimator, respectively, obtained after removing the

i

th observation. Note that all RC and RCSE in the fitted model, without the influential observation, are lower than

3 %

, except the RC and RCSE associated with the higher education level class. This estimate is more affected because there are only 6 women in the higher educational level in the dataset. However, the change in the corresponding p-value is small from slight below 5% to slight above 5%. This small change is not important in practical terms. To investigate if other simple count regression model provides a better fit to the number of antenatal care visits by mother in Oromia than the ZIP3, we fitted four other regression models considering the same covariates of the final model presented in Table 6. The considered regression models assume the following distributions for the response variable: ZIP2, Poisson (PO), negative binomial (NB) and zero-inflated negative binomial (ZINB). Table 9 presents the AIC and BIC of the five regression models. Note that the ZIP3 regression model has the lowest AIC and BIC in all cases, suggesting that this model provides a better fit to the number of antenatal care visits by mother in Oromia than its competitors. As expected, ZIP2 and ZIP3 have a similar AIC and BIC, suggesting that they provide a comparable fit.

Table 8.

RCs (in %) in ML estimates and in the corresponding estimated standard errors for the indicated removed case(s), and respective $p$ -values using data on antenatal care visits in Oromia—Ethiopia.

Remove cases	Submodel	Covariate	Category	$RC (\hat{θ})$	$RCSE (\hat{θ})$	$p$ -value
None	$μ$	Intercept		$\times$	$\times$	$< 0.0001$
		Education level	No education (ref)	$\times$	$\times$	0.0473
			Primary	$\times$	$\times$
			Secondary	$\times$	$\times$
			Higher	$\times$	$\times$
		Wealth	Poorest (ref)	$\times$	$\times$	0.0073
			Poorer	$\times$	$\times$
			Middle	$\times$	$\times$
			Richer	$\times$	$\times$
			Richest	$\times$	$\times$
	$ϕ$	Intercept		$\times$	$\times$	0.0080
		Wealth	Poorest (ref)	$\times$	$\times$	0.0007
			Poorer	$\times$	$\times$
			Middle	$\times$	$\times$
			Richer	$\times$	$\times$
			Richest	$\times$	$\times$
{94}	$μ$	Intercept		0.13	0.01	$< 0.0001$
		Education level	No education (ref)	$\times$	$\times$	0.0516
			Primary	0.54	0.04
			Secondary	0.97	0.07
			Higher	43.66	10.44
		Wealth	Poorest (ref)	$\times$	$\times$	$< 0.0001$
			Poorer	0.01	0.01
			Middle	0.07	0.01
			Richer	2.25	0.33
			Richest	1.61	0.01
	$ϕ$	Intercept		0.01	0.00	0.0080
		Wealth	Poorest (ref)	$\times$	$\times$	0.0007
			Poorer	0.00	0.00
			Middle	0.01	0.01
			Richer	0.36	0.12
			Richest	0.29	0.56

Table 9.

AIC and BIC for the final model for the antenatal care visits in Oromia—Ethiopia.

Model	AIC	BIC
ZIP	1809.7	1864.2
ZIP2	1809.7	1864.2
ZIP3	1803.7	1858.3
PO	2074.9	2108.5
NB	1984.2	2038.8
ZINB	1820.3	1895.9

7. Concluding remarks

In this work, we introduced a more interpretable regression model for count data with excess of zeros based on a reparameterization of the zero-inflated Poisson distribution. Inferential and diagnostic tools for this novel model were discussed. An application on number of antenatal care visits in Oromia, Ethiopia illustrated the usefulness of the proposed regression model.

The existing zero-inflated regression models are less interpretable than the other common models for count data. The parameters of our regression model are easily interpreted, especially when using the logarithmic link function as it was done in Section 6. Therefore, the proposed ZIP3 regression model will be very useful in medical research and also in other areas, when the response is a count variable with a high proportion of zeros.

When the response variable has more zeros than expected by the best known count distributions and a zero-inflated Poisson regression model is under consideration, it is of interest to know when to use which of the possible models. The ZIP regression model is the most suitable when the excess of zeros results from the existence of two subpopulations and the main interest is to study the mean response of the subpopulation that can assume values greater than zero. When there are two subpopulations, but the interest is the mean of the response variable instead of the mean of the response in one of the subpopulations, the ZIP2 regression model is more suitable. However, in practice, the excess of zeros is not usually the result of the existence of two subpopulations. When there is a single population, the ZIP3 regression model is the most appropriate when not only the mean but also the dispersion of the response variable is important. The application presented in Section 6 exemplifies a situation where the mean and dispersion of the response variable are important. Finally, when there is a single population and we are interested only in the mean of the response, then both ZIP3 and ZIP2 regression models can be used.

Apart from the ZIP regression model, there are several other regression models for count data with excess of zeros, such as the ZINB regression model, the ZIGP regression model,⁴⁹ the ZICMP regression model,⁵⁰ the ZIDP regression model,⁵¹ and the ZIBELL regression model.⁵² These models may outperform the ZIP3 regression model in some applications. Therefore, future research can propose reparameterizations for the distribution of the response variable of these models in order to make all parameters of the distribution of response variable interpretable.

We also acknowledge an anonymous reviewer for suggesting a shared-parameter⁵³ extension of the ZIP3 model, in which the linear predictor of the marginal mean is incorporated as the sole covariate in the dispersion component. Such a formulation offers several appealing features, including a reduction in the number of parameters, improved numerical stability in maximum likelihood estimation, and closer alignment with modeling strategies that avoid extensive model selection procedures. Moreover, under a log-link specification for the dispersion component, this approach may recover the NBII variance structure as a special case, providing a convenient framework for testing nested models. Taken together, these aspects suggest that a shared-parameter ZIP3 formulation may represent a parsimonious and competitive alternative to existing models such as ZINB2, constituting a promising avenue for future research, particularly in applications where interpretability, parameter parsimony, and model stability are of primary concern.

Footnotes

Acknowledgments

The authors thank the editors and reviewers for their constructive comments on an earlier version of this manuscript that resulted in significant improvement of the paper.

ORCID iDs

Gustavo H. A. Pereira

Jeremias Leao

Manoel Santos-Neto

Jianwen Cai

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was partially supported by São Paulo Research Foundation (FAPESP), grant number 2020/16334-9. The research was also partially supported by CNPq and CAPES grants from the Brazilian federal government, by FAPEAM grants from the government of the State of Amazonas. In addition, Jeremias Le\∼ao is supported by the Brazilian agencies CNPq (grant number 300393/2025-3).

Declaration of conflicting interest

The authors have no conflicting interests to declare with respect to the research, authorship, and/or publication of this article.

References

Brazel

Grant

Cabal

, et al. Baseline immunoglobulin G and immune function in non-Hodgkin lymphoma: a retrospective analysis. Front Immunol 2024; 15: 1334899.

Doan

Howell

Ball

, et al. Identifying areas of australia with high out-of-hospital cardiac arrest incidence and low bystander cardiopulmonary resuscitation rates: a retrospective, observational study. PLoS ONE 2024; 19: e0301176.

Hutchison

Jones

Ghosal

, et al. Comparison of in-person versus online comprehensive pelvic floor rehabilitation program following prostatectomy. Urology 2024; 190: 90–96.

Størdal

Tapia

Lund-Blix

, et al. Genotypes predisposing for celiac disease and autoimmune diabetes and risk of infections in early childhood. J Pediatr Gastroenterol Nutr 2024; 78: 295–303.

Byhoff

Dinh

Lucas

, et al. Mental health care use by ethnicity and preferred language in a national cohort of community health center patients. Psychiat Serv 2024; 75: 363–368.

Joundi

Hill

Stang

, et al. Association between time to treatment with endovascular thrombectomy and home-time after acute ischemic stroke. Neurology 2024; 102: e209454.

Lambert

. Zero-inflated poisson regression, with an application to defects in manufacturing. Technometrics 1992; 34: 1–14.

Ridout

Hinde

Demétrio

. A score test for testing a zero-inflated Poisson regression model against zero-inflated negative binomial alternatives. Biometrics 2001; 57: 219–223.

Rigby

Stasinopoulos

Heller

, et al. Distributions for modeling location, scale, and shape: using GAMLSS in R. CRC press, 2019.

10.

Long

Preisser

Herring

, et al. A marginalized zero-inflated Poisson regression model with overall exposure effects. Stat Med 2014; 33: 5151–5165.

11.

Preisser

Das

Long

, et al. Marginalized zero-inflated negative binomial regression with application to dental caries. Stat Med 2016a; 35: 1722–1735.

12.

Famoye

Preisser

. Marginalized zero-inflated generalized Poisson regression. J Appl Stat 2018; 45: 1247–1259.

13.

Todem

Kim

Hsu

. Marginal mean models for zero-inflated count data. Biometrics 2016; 72: 986–994. doi:https://doi.org/10.1111/biom.12492.

14.

Martin

Hall

. Marginal zero-inflated regression models for count data. J Appl Stat 2017; 44: 1807–1826.

15.

Long

Preisser

Herring

, et al. A marginalized zero-inflated Poisson regression model with random effects. J R Stat Soc Ser C: Appl Stat 2015; 64: 815–830.

16.

Burgette

Preisser

Weinberger

, et al. Impact of early head start in North Carolina on dental use among children younger than 3 years. Am J Public Health 2017; 107: 614–620.

17.

Mutiso

Pearce

Benjamin-Neelon

, et al. A marginalized zero-inflated negative binomial model for spatial data: modeling COVID-19 deaths in Georgia. Biomet J 2024; 66: e202300182.

18.

Cummings

Hardin

. Modeling count data with marginalized zero-inflated distributions. Stata J 2019; 19: 499–509.

19.

Laurence

Farmer-Dixon

Southwell

, et al. Sugar-sweetened beverage consumption and caries prevalence in underserved black adolescents. Pediatr Dent 2021; 43: 363–370.

20.

Mun

Zhou

Huh

, et al. Brief alcohol interventions are effective through 6 months: findings from marginalized zero-inflated Poisson and negative binomial models in a two-step IPD meta-analysis. Prev Sci 2023; 24: 1608–1621.

21.

Stasinopoulos

Rigby

. Generalized additive models for location scale and shape (gamlss) in r. J Stat Softw 2008; 23: 1–46.

22.

Zhou

Huh

, et al. mcount: Marginalized count regression models. 2022. DOI: 10.32614/CRAN.package.mcount. https://CRAN.R-project.org/package=mcount. R package.

23.

Seidel

Pazini

Tomazella

, et al. Predicting rice stem stink bug population dynamics based on GAMLSS models. Environ Entomol 2020; 49: 1145–1154.

24.

Sims

Tiwari

Levitan

, et al. Application of marginalized zero-inflated models when mediators have excess zeroes. Stat Methods Med Res 2024; 33: 148–161.

25.

Smyth

. Generalized linear models with varying dispersion. J R Stat Soc: Ser B (Methodol) 1989; 51: 47–60.

26.

Jorgensen

. The theory of dispersion models. CRC Press, 1997.

27.

Bonat

Jørgensen

Kokonendji

, et al. Extended Poisson–Tweedie: properties and regression models for count data. Stat Modell 2018; 18: 24–49.

28.

Bar-Lev

Ridder

. Exponential dispersion models for overdispersed zero-inflated count data. Commun Stat-Simul Comput 2023; 52: 3286–3304.

29.

David

Moore

. Notes on contagious distributions in plant populations. Ann Bot 1954; 18: 47–53.

30.

Hasan

Athanassiou

Hossain

. Estimating long-term spatial distribution of plodia interpunctella in various food facilities at Rajshahi municipality, Bangladesh, through pheromone-baited traps. Sci Rep 2022; 12: 15986.

31.

Zhang

Cai

Jin

, et al. Activity patterns, population dynamics, and spatial distribution of the stick tea thrips, dendrothrips minowai, in tea plantations. Insects 2023; 14: 152.

32.

Arnold

Manjunath

. Statistical inference for distributions with one Poisson conditional. J Appl Stat 2021; 48: 2306–2325.

33.

Vigas

Ortega

Cordeiro

, et al. The new neyman type a generalized odd log-logistic-g-family with cure fraction. J Appl Stat 2022; 49: 2805–2824.

34.

Johnson

. Hadamard products of matrices. Line Multil Algeb 1974; 1: 295–307.

35.

Rigby

Stasinopoulos

. Generalized additive models for location, scale and shape (with discussion). Appl Stat 2005; 54: 507–554.

36.

Feng

Sadeghpour

. A comparison of residual diagnosis tools for diagnosing regression models for count data. BMC Med Res Methodol 2020; 20: 1–21.

37.

Dunn

Smyth

. Randomized quantile residuals. J Comput Graph Stat 1996; 5: 236–244.

38.

Cook

Peña

Weisberg

. The likelihood displacement: a unifying principle for influence measures. Commun Stat-Theory Methods 1988; 17: 623–640.

39.

Cook

Peña

Weisberg

. Residuals and influence in regression. Chapman-Hall, 1982.

40.

Cortés

de Castro

Gallardo

. A new family of quantile regression models applied to nutritional data. J Appl Stat 2023; 51: 1378–1398.

41.

Fabio

Villegas

Carrasco

, et al. Diagnostic tools for a multivariate negative binomial model for fitting correlated data with overdispersion. Commun Stat-Theory Methods 2023; 52: 1833–1853.

42.

Ibacache-Pulgar

Villegas

López-Gonzales

, et al. Influence measures in nonparametric regression model with symmetric random errors. Stat Methods Appt 2023; 32: 1–25.

43.

Edessa

Dida

Teferi

. Early initiation of antenatal care and its associated factors among antenatal care followers at public health facilities in ambo town administration, Central Ethiopia. J Family Med Prim Care 2023; 12: 67–75.

44.

Woldeyohannes

Yohannes

Likassa

, et al. Count regression models analysis of factors affecting antenatal care utilization in Ethiopia. Ann Med Surgery 2023; 85: 4882–4886.

45.

Dinagde

Feyisa

Afework

, et al. Level of optimal antenatal care utilization and its associated factors among pregnant women in Arba Minch Town, Southern Ethiopia: new WHO-recommended ANC 8+ model. Front Global Women’s Health 2024; 5: 1259637.

46.

Bobo

Asante

Woldie

, et al. Evaluating equity across the continuum of care for maternal health services: analysis of national health surveys from 25 sub-Saharan African countries. Int J Equity Health 2023; 22: 239.

47.

The DHS Program. Demographic and health surveys (DHS) program. https://dhsprogram.com/ (2025, accessed 02 October 2025).

48.

Atkinson

. Two graphical displays for outlying and influential observations in regression. Biometrika 1981; 68: 13–20.

49.

Famoye

Singh

. Zero-inflated generalized poisson regression model with an application to domestic violence data. J Data Sci 2006; 4: 117–130.

50.

Sellers

Raim

. A flexible zero-inflated model to address data dispersion. Comput Stat Data Anal 2016; 99: 68–80.

51.

Xie

Lin

Wei

. Score tests for zero-inflated double Poisson regression models. Acta Mathemat Appl Sin English Ser 2017; 33: 851–864.

52.

Lemonte

Moreno-Arenas

Castellares

. Zero-inflated bell regression models for count data. J Appl Stat 2020; 47: 265–286.

53.

Preisser

Das

Long

, et al. Marginalized zero-inflated negative binomial regression with application to dental caries. Statistics in MedicineStat Med 2016b; 35: 1722–1735.