A Copula-Based Framework for Multivariate Zero-Inflated Mixed Poisson Models

Abstract: Multivariate count data often contain overdispersion, excess zeros, and complex dependence. Existing multivariate zero-inflated count models usually use a single dependence structure to jointly model structural zeros and count outcomes. This makes the two sources of dependence difficult to interpret separately and often leads to slow computation or intractable likelihoods. This paper proposes a general framework for multivariate zero-inflated mixed Poisson models that explicitly separates these two sources of dependence within a unified likelihood-based formulation. The proposed hierarchical model combines zero-inflated mixed Poisson marginals with separate dependence models for the structural-zero and latent-intensity components. Dependence among structural zeros is modeled by standard parametric copulas, while dependence among latent mixing variables is modeled by checkerboard copulas. Likelihood-based inference is carried out using an inference-functions-for-margins procedure, and the proposed likelihood also supports full maximum likelihood estimation when computationally feasible. Simulation studies show accurate parameter estimation and good finite-sample performance under different dependence structures and latent mixing distributions. Applications to benchmark datasets from existing statistical software and a healthcare utilization dataset demonstrate the flexibility, computational efficiency, and practical usefulness of the proposed framework.

• This paper presents a new hierarchical copula-based framework for multivariate zero-inflated mixed Poisson models, enhancing flexibility, scalability, and separability of dependence layers (for 10 - 34000 + demonstrations).
• For empirically provided examples: MEPS datasets in the healthcare sector, terror-related data in actuarial settings, and scRNA data in genomics, the framework's method outshines traditional methods, offering superior log-likelihood results and computational efficiency, making it valuable across insurance, healthcare, and genomics fields.
• The framework utilizes standard parametric copulas for structural-zero indicators and checkerboard copulas for latent-intensity variables, supports both IFM estimation and full maximum likelihood.

Motivation and contribution

Multivariate count data in insurance, healthcare, and genomics frequently exhibit overdispersion, excess zeros, and dependence among the response components simultaneously. Existing multivariate zero-inflated count models typically introduce dependence through a single mechanism—a copula, shared random effect, or latent factor—applied jointly to structural zeros and count outcomes. This conflates two substantively distinct sources of dependence: correlation in whether units enter the count-generating process at all, and correlation in latent intensities among those that do. In actuarial settings these sources carry different implications for risk classification, portfolio segmentation, and premium calculation. A second obstacle is computational: copula models for discrete responses require evaluating high-dimensional rectangle probabilities, which becomes expensive or infeasible as dimension grows.

This paper proposes a hierarchical framework that separates the two dependence layers explicitly while retaining a unified likelihood-based formulation (2608.12732). The framework accommodates a broad class of mixed Poisson marginals (Gamma, Lognormal, Pareto, Weibull, Inverse Gaussian), uses standard parametric copulas for the structural-zero indicators and checkerboard copulas for the latent-intensity variables, and supports both inference-functions-for-margins (IFM) estimation and full maximum likelihood when feasible.

Model structure

For each margin $j=1,\ldots,k$, an observation is generated through three layers: a Bernoulli indicator $B_j$ determining membership in the structural-zero state, a latent intensity $\Lambda_j$ drawn from a positive continuous mixing distribution $F_j(\mu_j,\phi_j)$ with mean $\mu_j$ and variance function $\phi_jV(\mu_j)$, and a conditional Poisson layer with $Y_j = B_jN_j$, where $N_j \mid \Lambda_j \sim \text{Poisson}(\Lambda_j)$. Regression structures are specified via logistic links for the zero-state probabilities $\pi_{ij}$ and log links for the conditional means $\mu_{ij}$. The indicators $B_j$0 and intensities $B_j$1 are assumed mutually independent; dependence within each layer is induced by separate copulas.

The key methodological device concerns the latent-intensity layer. The exact joint probability mass function of the latent counts involves a $B_j$2-dimensional integral of the product of Poisson kernels against an arbitrary copula density—tractable only at prohibitive cost. Replacing the copula density by a checkerboard copula (piecewise constant on an equal-width partition of the unit cube with resolution $B_j$3) converts this integral into a finite weighted sum: the joint pmf becomes $B_j$4, where the cell probabilities $B_j$5 depend only on the marginal mixing distributions and are computable by one-dimensional integration or closed form.

Gaussian factor construction for scalability

Direct checkerboard weight computation requires evaluating $B_j$6 cell probabilities, which is impractical beyond small $B_j$7. The paper parameterizes the checkerboard weights through a $B_j$8-factor Gaussian copula, reducing dependence parameters from $B_j$9 to $\Lambda_j$0. Under this representation, the joint pmf of either the latent counts alone or the full zero-inflated response vector reduces to a $\Lambda_j$1-dimensional Gaussian integral over the factor vectors of the two dependence layers, evaluated by Gauss–Hermite quadrature. Crucially, quadrature dimension depends on the number of factors rather than the number of responses $\Lambda_j$2.

Estimation proceeds by IFM: marginals are estimated first per margin, then checkerboard cell probabilities are fixed and dependence parameters optimized in a second stage. Closed-form marginal probabilities exist for Gamma (via negative binomial) and Inverse Gaussian (via half-integer Bessel expansions); the remaining mixtures use Gauss–Hermite, Gauss–Legendre, or Gauss–Laguerre quadrature. The appendixes provide semi-analytical score functions throughout.

Simulation results

Three experiment blocks support the methodology. First, with $\Lambda_j$3 per scenario, the correctly specified mixing distribution achieved the highest log-likelihood in every scenario across five candidate mixtures, under both mixed Poisson and zero-inflated specifications and two heterogeneity levels ($\Lambda_j$4 and $\Lambda_j$5). Regression coefficients were recovered accurately (e.g., true $\Lambda_j$6 estimated within roughly 0.05 across mixtures). One notable exception: under Pareto mixing with higher heterogeneity, the dispersion estimate was substantially biased ($\Lambda_j$7 versus true 1.500), which the authors attribute to weak identifiability of the Lomax shape rather than estimation failure—an assumption stated but not formally established.

Second, bivariate copula recovery experiments spanning Gaussian, Student-$\Lambda_j$8, Clayton, Gumbel, Frank, and rotated families at $\Lambda_j$9 showed correct identification in nearly all cases, with dependence parameters closely recovered. Under zero inflation, recovery remained accurate except that Gumbel-generated data were occasionally better fit by the survival Clayton copula, with negligible log-likelihood differences—the authors acknowledge these two families produce near-equivalent dependence here. Third, Gaussian factor loadings were recovered accurately for $F_j(\mu_j,\phi_j)$0 margins under both one-factor and two-factor structures.

Empirical performance

Four applications benchmark the framework against established software:

• MEPS (bivariate healthcare utilization): fit comparable to GJRM despite dependence modeled through latent intensities rather than observed responses; the Gamma–Lognormal marginal combination outperformed homogeneous specifications, and the survival Gumbel copula was selected.
• Terror data: the proposed model achieved a higher log-likelihood than bizicount ($F_j(\mu_j,\phi_j)$1 vs.\ $F_j(\mu_j,\phi_j)$2) while separating countermonotone structural-zero dependence from weakly positive latent-intensity dependence—an interpretable decomposition unavailable under single-copula alternatives.
• VHLSS (>34,000 observations): the best specification (IG/LN margins, survival Clayton checkerboard) improved log-likelihood from $F_j(\mu_j,\phi_j)$3 (bizicount Frank) to $F_j(\mu_j,\phi_j)$4, and complete model selection took minutes versus several hours for bizicount. This computational advantage follows directly from replacing discrete rectangle probabilities with continuous copula densities plus univariate mixed-Poisson likelihoods.
• scRNA data ($F_j(\mu_j,\phi_j)$5 and $F_j(\mu_j,\phi_j)$6): the proposed framework dominated PLNmodels on AIC and BIC. At $F_j(\mu_j,\phi_j)$7, the optimal specification used 450 dependence parameters under three factors versus 1,525 in PLN's unrestricted covariance matrix. For ten highly zero-inflated genes, modeling both dependence layers jointly yielded the best fit (log-likelihood $F_j(\mu_j,\phi_j)$8).

A caveat applies to the scRNA comparison: because PLNmodels optimizes an evidence lower bound rather than the observed-data log-likelihood, cross-method log-likelihood comparisons are not strictly like-for-like; the information criteria comparisons partially mitigate this but inherit the same caveat.

Limitations and open questions

Several restrictions bound the applicability of the results. The mutual independence assumption between $F_j(\mu_j,\phi_j)$9 and $\mu_j$0 enables interpretability but precludes models in which zero-propensity and intensity correlate—for example, individuals systematically less likely to claim who would also claim less often if they did. The framework's efficiency rests on equally spaced checkerboard partitions; more general partitions are acknowledged possible but undeveloped. Only the Gaussian factor copula is implemented for weight construction, although the likelihood formulation admits arbitrary factor copulas whose conditional probabilities can be evaluated. Quadrature cost grows with the number of latent factors, so very high dimensions with complex dependence remain challenging. Finally, identifiability of the Pareto dispersion parameter under strong heterogeneity is empirically weak, and no theoretical identifiability analysis of the checkerboard weights themselves is provided.

Conclusion

The paper delivers a modular, likelihood-based framework in which excess zeros, overdispersion, and multivariate dependence are addressed by separately specified components: flexible mixing distributions for marginals, parametric copulas for structural-zero dependence, and checkerboard copulas—factor-parameterized for scalability—for latent-intensity dependence. Simulations demonstrate consistent recovery of mixing distributions, regression parameters, copula families, and factor loadings; applications show competitive or superior fit relative to GJRM, bizicount, and PLNmodels with substantial computational savings in large datasets. The principal open questions concern relaxing the independence assumption between the two latent layers, developing non-Gaussian or semiparametric factor constructions within the same likelihood machinery, and establishing formal identifiability conditions for the checkerboard representation.