Results

Foundational count data models
Poisson regression model
The Poisson regression model serves as the fundamental approach for count data analysis. For a response variable Yᵢ, the Poisson distribution has the probability mass function:

P(Yᵢ = yᵢ | μᵢ) = exp(-μᵢ) μᵢ^(yᵢ) / yᵢ!

where yᵢ = 0, 1, 2, … and μᵢ > 0 is the expected count. The Poisson distribution has the property that E(Yᵢ) = Var(Yᵢ) = μᵢ, which constitutes the equidispersion assumption. The log-linear link function relates the expected count to the linear predictor:

log(μᵢ) = β₀ + β₁xᵢ₁ + β₂xᵢ₂ + … + βₚxᵢₚ

where β are regression coefficients, and x are covariates. Maximum likelihood estimation is used to obtain parameter estimates [17]. The log-likelihood function for independent observations is:

l(β) = Σᵢ [yᵢ log(μᵢ) – μᵢ – log(yᵢ!)]

The model’s computational simplicity and interpretability make it appealing for initial count data analysis; however, the equidispersion assumption often fails in practice, necessitating more flexible approaches [18].

Negative binomial regression model
The Negative Binomial (NB) model addresses overdispersion by introducing an additional dispersion parameter. The NB distribution can be derived as a Poisson-gamma mixture. The probability mass function is:

P(Yᵢ = yᵢ | μᵢ, θ) = Γ(yᵢ + θ) / (Γ(θ) yᵢ!) × (θ / (θ + μᵢ))^θ × (μᵢ / (θ + μᵢ))^yᵢ

where θ > 0 is the dispersion parameter. This mixture yields a variance structure where Var(Yᵢ) = μᵢ + μᵢ²/θ, which exceeds the mean. As θ → ∞, the NB distribution approaches the Poisson distribution. The NB model maintains the same log-linear link function as Poisson regression [19].

Generalized poisson model
The Generalized Poisson (GP) distribution provides additional flexibility by accommodating both overdispersion and underdispersion [20]. The probability mass function is:

P(Yᵢ = yᵢ | μᵢ, λ) = (μᵢ / (1 + λ μᵢ))^yᵢ × (1 + λ yᵢ)^(yᵢ-1) / yᵢ! × exp[-μᵢ(1 + λ yᵢ)/(1 + λ μᵢ)]

where λ is the dispersion parameter. When λ = 0, the GP distribution reduces to the Poisson distribution; when λ > 0, it exhibits overdispersion; and when -1/μ < λ < 0, it exhibits underdispersion. This parameterization allows the GP model to serve as a unified framework for count data with varying dispersion patterns [21].

Offset and exposure adjustment
Count data often arise from observation units with differing exposure sizes, such as person-time, population size, area, or number of trials [22]. When exposure varies across observations, the expected count should be modelled relative to exposure rather than as an absolute count. This is accomplished by including an offset variable in the linear predictor [23]. For the Poisson model, the log-linear link becomes: log(μᵢ) = log(tᵢ) + β₀ + β₁xᵢ₁ + β₂xᵢ₂ + … + βₚxᵢₚ where tᵢ is the known exposure or offset for observation i. Equivalently, the model can be written as: log (μᵢ / tᵢ) = β₀ + β₁xᵢ₁ + β₂xᵢ₂ + … + βₚxᵢₚ The coefficient of log(tᵢ) is fixed at 1 rather than estimated as an ordinary regression coefficient. Offset variables are used in the same way for Negative Binomial, Generalized Poisson, zero-inflated, hurdle, and other count regression models, because the offset adjusts the expected count for differential exposure while the covariate effects are estimated on the rate scale [24]. Failure to include an appropriate offset when exposure varies can bias estimated covariate effects and lead to misleading inference [25].

Models for zero-inflated count data
Zero-inflated Poisson model
The Zero-Inflated Poisson (ZIP) model addresses excess zeros in count data by assuming that the data arise from a mixture of two distributions. The probability mass function is: P(Yᵢ = yᵢ | πᵢ, μᵢ) = πᵢ + (1 – πᵢ) exp(-μᵢ), if yᵢ = 0 P(Yᵢ = yᵢ | πᵢ, μᵢ) = (1 – πᵢ) exp(-μᵢ) μᵢ^(yᵢ) / yᵢ!, if yᵢ > 0

where πᵢ represents the probability of belonging to the zero state (i.e., the observation is structurally zero). The zero-inflation probability can depend on covariates through a logistic model:

logit(πᵢ) = γ₀ + γ₁zᵢ₁ + … + γₘzᵢₘ,

while the count component uses a log-linear model. The ZIP model accommodates overdispersion arising from excess zeros [26].

Zero-inflated negative binomial model
The Zero-Inflated Negative Binomial (ZINB) model extends the ZIP framework by using an NB distribution for the count component. The probability mass function is:

P(Yᵢ = 0 | πᵢ, μᵢ, θ) = πᵢ + (1 – πᵢ) (θ / (θ + μᵢ))^θ

P(Yᵢ = yᵢ | πᵢ, μᵢ, θ) = (1 – πᵢ) Γ(yᵢ + θ) / (Γ(θ) yᵢ!) × (θ / (θ + μᵢ))^θ × (μᵢ / (θ + μᵢ))^yᵢ, for yᵢ > 0

The ZINB model provides greater flexibility than ZIP because it accommodates both excess zeros and overdispersion simultaneously through inclusion of both a zero-inflation parameter and a dispersion parameter. This dual capability makes the ZINB model particularly suitable for datasets where both sources of variability are present [27].

Hurdle models
Hurdle models consist of two parts: a binary process determining whether the count is zero or positive, and a truncated count process for positive observations. The probability mass function is: P(Yᵢ = 0) = φᵢ

P(Yᵢ = yᵢ) = (1 – φᵢ) P(Y = yᵢ | y > 0), for yᵢ > 0

where φᵢ is the hurdle probability. The positive counts are modeled with a count distribution truncated to exclude zeros [28]. For a Poisson hurdle model, the positive part is: P(Y = y | y > 0) = exp(-μ) μ^y / (y! (1 – exp(-μ)))

Hurdle models differ from zero-inflated models in their interpretation: zeros in hurdle models arise from a single process (the hurdle), whereas zeros in zero-inflated models arise from two possible processes. This distinction has important implications for model selection and interpretation, particularly when the mechanisms generating zeros are of substantive interest [11].

Advanced models
Finite mixture models
Finite mixture models capture unobserved heterogeneity by assuming the population consists of multiple latent subpopulations, each with its own count distribution. The overall distribution is:

f(yᵢ | Θ) = Σ_{k=1}^K π_k f_k(yᵢ | θ_k)

where K is the number of components, π_k are mixing proportions (0 < π_k < 1, Σπ_k = 1), and f_k are component distributions from the exponential family [29]. The mixture model provides flexible modeling of complex distributions and can accommodate multiple modes and other non-standard features. This flexibility makes finite mixture models valuable for exploring latent structures in count data, though they require careful specification of the number of components and are computationally intensive [30].

Panel count models with random effects
Panel count models account for correlation across repeated observations within individuals or groups by incorporating random effects into the linear predictor. The mixed-effects framework adds subject-specific terms that capture within-subject correlation:

log(μᵢⱼ) = β₀ + β₁xᵢⱼ₁ + … + βₚxᵢⱼₚ + bᵢ

where bᵢ ~ N(0, σ²) are random effects [31]. The model accommodates both overdispersion and correlation among repeated measures, making it appropriate for longitudinal count data. Estimation is often performed using numerical integration techniques or through Bayesian hierarchical approaches [32].

Estimation methods
Maximum likelihood estimation
Maximum likelihood estimation (MLE) is the primary approach for parameter estimation in count data models. The log-likelihood function is specified based on the assumed distribution. For the Poisson model, the log-likelihood is: l(β) = Σᵢ [yᵢ log(μᵢ) – μᵢ – log(yᵢ!)] Parameter estimates are obtained by maximising the log-likelihood function, typically using iterative numerical optimization algorithms. MLE provides consistent and asymptotically efficient estimates under regularity conditions; furthermore, it provides a framework for constructing confidence intervals and hypothesis tests based on the asymptotic normality of the estimators [33].

Expectation maximization algorithm
The Expectation Maximization (EM) algorithm is particularly useful for mixture models and zero-inflated models where the likelihood function involves latent variables. The algorithm alternates between two steps: Expectation Step (E-step): Compute the expected complete-data log-likelihood conditioning on observed data and current parameter estimates:

Q(θ | θ^(t)) = E[log L(θ | Y, Z) | Y, θ^(t)]

where Z represents latent variables (e.g., component membership indicators in mixture models or zero-state indicators in ZIP models). Maximization Step (M-step): Maximize the expected complete-data log-likelihood to obtain updated parameter estimates:

θ^(t+1) = argmax Q(θ | θ^(t))

The EM algorithm provides stable convergence and naturally handles missing data structures, including latent class memberships and zero-inflation indicators [34]. Although convergence can be slow, the EM algorithm is widely used because each maximization step typically has a closed-form solution or involves simple standard estimation procedures [35].

Model diagnostics and selection  Residual analysis
Residual analysis provides critical information about model fit and helps identify violations of model assumptions. For count data models, Pearson residuals are defined as:

r_i = (yᵢ – μ̂ᵢ) / √Var(Yᵢ)

Deviance residuals are defined as:

d_i = sign(yᵢ – μ̂ᵢ) √{2 [yᵢ log(yᵢ / μ̂ᵢ) – (yᵢ – μ̂ᵢ)]}

Randomized quantile residuals are often preferred for their improved performance in detecting model misspecification. These residuals transform the data to approximate normality, thereby facilitating standard diagnostic procedures [36].

Overdispersion tests
Tests for overdispersion help determine whether the NB model is more appropriate than the Poisson model. Common approaches include:

Likelihood Ratio Test (LRT): Compares the log-likelihood of the Poisson model (constrained, θ = 0) against the NB model (unconstrained). The test statistic is: LR = 2(l_NB – l_Poisson) which follows a mixture of chi-square distributions.
Score test: Assesses whether the dispersion parameter significantly differs from zero without estimating the alternative model [37].
Wald test: Evaluates the significance of the dispersion parameter using: W = θ̂² / Var(θ̂) These tests are important because selecting a model that incorrectly assumes equidispersion can lead to invalid inference, while using an overdispersed model when not needed may result in a loss of efficiency [38].

Zero-inflation tests
Tests for zero-inflation help determine whether zero-inflated models are necessary. The Vuong test is commonly used for comparing non-nested models such as ZIP versus Poisson. The test statistic is:

Vuong = (1/√n) Σᵢ log[f₁(yᵢ | xᵢ) / f₂(yᵢ | xᵢ)] / SD[log(f₁/f₂)]

where f₁ and f₂ are the probability mass functions of the competing models [39]. Positive values favor the zero-inflated model. Other approaches include score tests and information criteria such as AIC and BIC, which can be used to compare nested and non-nested models [40].

Additional count regression approaches
Quasi-Poisson regression
Quasi-Poisson regression addresses overdispersion by introducing a dispersion parameter φ into the variance function: Var(Yᵢ) = φ μᵢ where φ > 0 is estimated from the data. This approach does not specify a full probability distribution; instead, it uses quasi-likelihood estimation. The regression coefficients are estimated identically to Poisson regression; however, standard errors are multiplied by √φ to account for overdispersion. Quasi-Poisson is computationally simple and provides robust inference; nevertheless, it does not provide a likelihood for model comparison [41].

Conway-maxwell-poisson regression
The Conway-Maxwell-Poisson (CMP) distribution generalizes the Poisson distribution by adding a dispersion parameter ν that controls the rate of decay of successive probabilities. The probability mass function is:

P(Y = y) = (λ^y) / ((y!)^ν Z(λ, ν))

where Z(λ, ν) = Σ_{j=0}^∞ λ^j / (j!)^ν is a normalizing constant. When ν = 1, the CMP reduces to the Poisson distribution; when ν < 1, it accommodates overdispersion; when ν > 1, it accommodates underdispersion. This model provides a unified framework for modeling count data with varying dispersion patterns [42].

Generalized estimating equations
Generalized Estimating Equations (GEE) extend generalized linear models to correlated count data. The estimating equations are:

Σᵢ Dᵢᵀ Vᵢ⁻¹ (yᵢ – μᵢ) = 0

where Dᵢ = ∂μᵢ/∂β, Vᵢ = Aᵢ^(1/2) R(α) Aᵢ^(1/2) is the working covariance matrix, Aᵢ is a diagonal matrix of variance functions, and R(α) is a working correlation matrix. GEE provides consistent estimates even when the correlation structure is misspecified; however, inference depends on robust standard errors [43].

Model selection framework
As outlined in Figure 2, the methodological decision framework for count data modeling follows a structured, sequential workflow. First, Step 1 involves assessing the distribution of the outcome variable: if there are no or few zeros (zero proportion < 5%), Poisson or negative binomial (NB) models are considered; if there are excess zeros (zero proportion > 15%), zero-inflated or hurdle models are appropriate; and if all counts are positive (no zeros), zero-truncated models should be considered.

Second, Step 2 requires assessing dispersion by calculating the variance-to-mean ratio, where a variance approximately equal to the mean points to a Poisson model, a variance greater than the mean suggests NB, generalized Poisson (GP), or quasi-Poisson models, and a variance less than the mean indicates GP or Conway-Maxwell-Poisson (CMP) models. Third, Step 3 accounts for data structure, noting that independent observations use standard count models, clustered or repeated measures call for generalized estimating equations (GEE), random effects, or panel models, and longitudinal designs with subject-specific effects require mixed-effects count models.

Finally, Step 4 evaluates the model through diagnostic tests, which involve conducting overdispersion tests (likelihood ratio, score, or Wald tests), running zero-inflation tests (Vuong or score tests), examining residual diagnostics, and comparing information criteria such as AIC and BIC (Figure 2).

Model comparison summary
Table 2 presents a comparison of count data models, including their distributions, handling of overdispersion, underdispersion, zero inflation, and correlation, along with key limitations.

Strengths and limitations of this review
This review offers several strengths that enhance its value for researchers and practitioners. First, the review provides comprehensive and systematic coverage of count data models, ranging from foundational Poisson regression to advanced finite mixture and panel models; therefore, readers gain a complete picture of the available analytical toolkit. Second, the review clearly distinguishes between related but conceptually different approaches, such as zero-inflated models versus hurdle models, which helps researchers avoid common misapplications. Third, the inclusion of both estimation techniques and diagnostic tools offers practical guidance beyond theoretical description; consequently, readers can implement these models effectively. Fourth, the review synthesizes recent developments in the field, including generalized Poisson distributions and Bayesian estimation approaches, thereby ensuring that the content is current and relevant. Fifth, the structured decision-making framework assists researchers in selecting appropriate models based on specific data characteristics.

Despite these strengths, several limitations should be acknowledged. First, this review is methodological in nature and does not include empirical data analysis or simulation studies; therefore, the practical performance of the models under different conditions is not demonstrated directly. Second, the scope is primarily focused on frequentist estimation approaches; although Bayesian methods are mentioned, their detailed implementation and prior specification are not thoroughly covered. Third, emerging machine learning approaches for count data, such as tree-based methods and neural networks, are not included; hence, the review may not fully represent the most recent computational advances. Fourth, this review assumes that readers have a basic statistical background; hence, beginners may find some technical details challenging without additional supporting materials. Fifth, the coverage of spatial and spatiotemporal count models remains limited. Acknowledging these limitations helps position this review as a starting point rather than a definitive guide, and it encourages readers to consult specialised literature for deeper exploration of specific topics.

Funding

The authors did not receive any specific funding for this work.