Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Multiple inflated negative binomial regression for correlated multivariate count data

Multiple inflated negative binomial regression for correlated multivariate count data 1IntroductionGeneralized linear models (GLMs) have become popular for modeling count data in which the underlying distribution of the response variables is assumed Poisson or negative binomial distribution. In cases where the equidispersion assumption of the Poisson model is violated, the negative binomial distribution with size r(>0)r\left(\gt 0)and mean μ(>0)\mu \left(\gt 0)given by (1.1)fNB(y∣r,μ)=Γ(y+r)Γ(r)y!rr+μrμr+μy,y=0,1,…{f}_{{\rm{NB}}}(y| r,\mu )=\frac{\Gamma (y+r)}{\Gamma \left(r)y\hspace{0.1em}\text{&#x0021;}\hspace{0.1em}}{\left(\frac{r}{r+\mu }\right)}^{r}{\left(\frac{\mu }{r+\mu }\right)}^{y},\hspace{1em}y=0,1,\ldots is often preferred. From [10], the variance of YYis given by (1.2)Var(Y)=μ+μ2r.\hspace{0.1em}\text{Var}\hspace{0.1em}\left(Y)=\mu +\frac{{\mu }^{2}}{r}.See [10] for a full derivation of (1.1) using a Poisson-gamma mixture. Notice that as r→∞r\to \infty , the distribution reduces to the Poisson in which E(Y)=Var(Y)=μE\left(Y)=\hspace{0.1em}\text{Var}\hspace{0.1em}\left(Y)=\mu . Thus, the dispersion parameter α=1/r\alpha =1\hspace{0.1em}\text{/}\hspace{0.1em}rallows one to directly model overdispersed data, unlike the traditional Poisson regression model. [3] showed that (1.1) is often referred to as the negative binomial-II (NB-II) model due to the quadratic term in (1.2). In general, negative binomial-P (NB-P) models have been given in the literature, such as in [8]. However, they are not explored in this article. See [10] for more on traditional negative binomial regression models.In certain cases, count data contain excessive zero-valued counts and violate both the Poisson and negative binomial distribution assumptions. Both hurdle and zero-inflation Poisson and negative binomial distribution models such as in [9,14,18] have been proposed to solve this problem. Essentially, both models assume zero-valued counts are generated from a second process. However, unlike hurdle models, zero-inflation models assume that zeroes are generated from both a second process and the underlying data generating process. For example, consider the work-week hours for males/females generated from a distribution with nonnegative integer-valued support (discussed in Section 7.2). Then there are two types of zero-valued counts: those unemployed and those unemployed due to a disability. That is, one type cannot work, while the other can, though both counts are zero-valued. Zero-inflation models allow one to model these situations through a binomial-count mixture model. More recently, doubly inflated models have been proposed by [23,24] for the Poisson distribution for cases when there are an excessive amount of zero-valued and kk-valued counts, where kkis a positive integer, such that the distributional assumptions of zero-inflation models also are not met. Mathews et al. [16] provided a method of finding doubly inflated multivariate negative binomial distributions using Gaussian copula. However, they did not consider the possibility of having multiple inflation points in the data. Here, we provide a method of finding multivariate distribution functions allowing arbitrary inflation points in the data using a copula. In this case, the standard zero-inflated negative binomial distribution is generalized to a multinomial-negative binomial mixture distribution. We also propose a regression method for correlated multiple inflated multivariate count data.Care must be taken in choosing an appropriate form of the multivariate negative binomial distribution. Numerous bivariate forms of the negative binomial distribution have been proposed in the literature. These are often constructed using a product of Poisson mass functions paired with a mixing distribution, such as the gamma distribution or log-normal model [1,13]. For example, [15] derived the bivariate negative binomial distribution from the probability mass functions of Y1∼Pois(αθ){Y}_{1}\hspace{0.33em} \sim \hspace{0.33em}{\rm{Pois}}\left(\alpha \theta ), Y2∼Pois(βθ){Y}_{2}\hspace{0.33em} \sim \hspace{0.33em}{\rm{Pois}}\left(\beta \theta ), where α>0\alpha \gt 0, β>0\beta \gt 0, θ>0\theta \gt 0, and a gamma mixing distribution with shape parameter rrand scale parameter λ=1\lambda =1: (1.3)∫0∞(αθ)y1e−αθy1!(βθ)y2e−βθy2!θr−1e−θΓ(r)dθ.\underset{0}{\overset{\infty }{\int }}\frac{{\left(\alpha \theta )}^{{y}_{1}}{e}^{-\alpha \theta }}{{y}_{1}\hspace{0.1em}\text{&#x0021;}\hspace{0.1em}}\frac{{\left(\beta \theta )}^{{y}_{2}}{e}^{-\beta \theta }}{{y}_{2}\hspace{0.1em}\text{&#x0021;}\hspace{0.1em}}\frac{{\theta }^{r-1}{e}^{-\theta }}{\Gamma \left(r)}{\rm{d}}\theta .Unfortunately, the aforementioned model only allows for nonnegative correlation values. Efforts have been made by [5,17] and others to relax this restriction. One possible solution is found by forming the multivariate negative binomial distribution through a copula. Indeed, copula theory has been growing in popularity for modeling count data, for example, [19] and [27]. A copula is a multivariate probability distribution with uniform marginal distributions and dependence parameter Ω{\boldsymbol{\Omega }}. A fundamental result in copula theory is Sklar’s theorem, which allows one to form a multivariate probability distribution through a copula with only knowledge of the form of the marginal distributions. In this article, we form a multivariate negative binomial distribution using a copula, in which Ω{\boldsymbol{\Omega }}represents the copula parameters.This article is organized as follows: in Section 2, we review the copula function, Sklar’s theorem, and the Gaussian copula; in Section 3, we review multivariate discrete distributions derived using copula functions and the marginal negative binomial distribution; in Section 4, we derive the multivariate multipleinflated negative binomial (MMINB) model using the copula function; in Section 5, we extend the results of Section 4 to a regression setting; in Section 6, we provide an iterative procedure for estimating the MMINB regression model parameters using maximum likelihood estimation (MLE); in Section 7.1, we consider simulated data and fit the data with proposed MMINB model and the usual multivariate negative binomial regression model using Gaussian copula; in Section 7.2, we apply the model to the cps91 dataset from the Wooldridge package in R, see [31] and [25]. We end with concluding remarks in Section 8.2Copulas and Sklar’s theoremThis section describes how the copula function may be used to form a multivariate probability mass function. We refer to [12] and [29] for more information on copulas than presented here. A copula is a multivariate probability distribution with uniform marginal distributions on the interval [0,1]\left[0,1].Definition 1A function C:[0,1]d→[0,1]C:{\left[0,1]}^{d}\to \left[0,1]is a dd-dimensional copula if the following are true: (1)C(1,…,yj,…,1)=yj,∀j=1,2,…,dC\left(1,\ldots ,{y}_{j},\ldots ,1)={y}_{j},\hspace{1em}\forall j=1,2,\ldots ,dwith (1,…,yj,…,1)∈[0,1]d\left(1,\ldots ,{y}_{j},\ldots ,1)\in {\left[0,1]}^{d},(2)C(y1,y2,…,yd)=0C({y}_{1},{y}_{2},\ldots ,{y}_{d})=0for (y1,y2,…,yd)∈[0,1]d({y}_{1},{y}_{2},\ldots ,{y}_{d})\in {\left[0,1]}^{d}with at least one of yj=0{y}_{j}=0for j=1,2,…,dj=1,2,\ldots ,d,(3)CCis dd-nondecreasing, i.e., for any u1=(u1(1),u2(1),…,ud(1)){{\bf{u}}}_{1}=\left({u}_{1}^{\left(1)},{u}_{2}^{\left(1)},\ldots ,{u}_{d}^{\left(1)}), u2=(u1(2),u2(2),…,ud(2)){{\bf{u}}}_{2}=\left({u}_{1}^{\left(2)},{u}_{2}^{\left(2)},\ldots ,{u}_{d}^{\left(2)})belong to [0,1]d{\left[0,1]}^{d}with uj(1)≤uj(2){u}_{j}^{\left(1)}\le {u}_{j}^{\left(2)}, for all j=1,2,…,dj=1,2,\ldots ,d:(2.1)∑a1=12∑a2=12…∑ad=12(−1)a1+a2+…+adC(u1(a1),u2(a2),…,ud(ad))≥0.\mathop{\sum }\limits_{{a}_{1}=1}^{2}\mathop{\sum }\limits_{{a}_{2}=1}^{2}\ldots \mathop{\sum }\limits_{{a}_{d}=1}^{2}{\left(-1)}^{{a}_{1}+{a}_{2}+\ldots +{a}_{d}}C\left({u}_{1}^{\left({a}_{1})},{u}_{2}^{\left({a}_{2})},\ldots ,{u}_{d}^{\left({a}_{d})})\ge 0.Sklar’s theorem, given in [26], is a fundamental result in copula theory as it allows one to form a joint distribution using known marginal distributions.Theorem 1Let Y1,Y2,…,Yd{Y}_{1},{Y}_{2},\ldots ,{Y}_{d}be random variables with marginal distribution functions F1,F2,…,Fd{F}_{1},{F}_{2},\ldots ,{F}_{d}and joint cumulative distribution function F, then the following are true: (1)There exists a d-dimensional copula C such that for all y1,y2,…,yd∈R{y}_{1},{y}_{2},\ldots ,{y}_{d}\in {\mathbb{R}}(2.2)F(y1,y2,…,yd)=C(F1(y1),F2(y2),…,Fd(yd)).F({y}_{1},{y}_{2},\ldots ,{y}_{d})=C\left({F}_{1}({y}_{1}),{F}_{2}({y}_{2}),\ldots ,{F}_{d}({y}_{d})).(2)If Y1,Y2,…,Yd{Y}_{1},{Y}_{2},\ldots ,{Y}_{d}are continuous then the copula C is unique. Otherwise, C can be uniquely determined on a d-dimensional rectangle Range(F1)×Range(F2)×…×Range(Fd){\rm{Range}}\left({F}_{1})\times {\rm{Range}}\left({F}_{2})\times \ldots \times {\rm{Range}}\left({F}_{d}).Copulas are popular because of their ability to model dependence among Y1,Y2,…,Yd{Y}_{1},{Y}_{2},\ldots ,{Y}_{d}. One popular copula is the Gaussian copula since this dependence is modeled through a correlation matrix RR.Definition 2The Gaussian copula is given by the function: (2.3)C(u1,u2,…,ud∣R)=ΦR(Φ−1(u1),Φ−1(u2),…,Φ−1(ud)),C\left({u}_{1},{u}_{2},\ldots ,{u}_{d}| R)={{\boldsymbol{\Phi }}}_{R}\left({\Phi }^{-1}\left({u}_{1}),{\Phi }^{-1}\left({u}_{2}),\ldots ,{\Phi }^{-1}\left({u}_{d})),where Φ−1{\Phi }^{-1}is the inverse cumulative distribution function of the standard Gaussian distribution and ΦR{{\boldsymbol{\Phi }}}_{R}is the joint cumulative distribution function of a standard multivariate Gaussian distribution with correlation matrix RR.Definition 3The Gaussian copula density [2] is given by (2.4)c(u1,u2,…,ud∣R)=1∣R∣exp−12UT×(R−1−Id)×U,c\left({u}_{1},{u}_{2},\ldots ,{u}_{d}| R)=\frac{1}{\sqrt{| R| }}\exp \left\{-\frac{1}{2}{{\bf{U}}}^{T}\times ({R}^{-1}-{I}_{d})\times {\bf{U}}\right\},where U=[Φ−1(u1),Φ−1(u2),…,Φ−1(ud)]T{\bf{U}}={{[}{\Phi }^{-1}\left({u}_{1}),{\Phi }^{-1}\left({u}_{2}),\ldots ,{\Phi }^{-1}\left({u}_{d})]}^{T}, Φ\Phi is the univariate standard Gaussian distribution, RRis the correlation matrix of the standard multivariate Gaussian distribution, and Id{I}_{d}be the d×dd\times didentity matrix.For high-dimensional problems, the Gaussian copula model quickly becomes intractable. One reason for this intractability is the large correlation matrix that must be estimated from the data. To overcome this issue, we impose a correlation structure of two types: (1)Equi-correlation structure: Under this structure, we assume RRas a function of ρ\rho given by(2.5)R(ρ)=ρ11T−(1−ρ)Id,R\left(\rho )=\rho {\bf{1}}{{\bf{1}}}^{T}-\left(1-\rho ){I}_{d},where Id{I}_{d}is a dd-dimensional identity matrix, ρ∈−1d−1,1\rho \in \left(-\frac{1}{d-1},1\right), and 1{\bf{1}}is a dd-dimensional column vector of ones. From [20], we have(2.6)R−1(ρ)=11−ρId−ρ(1−ρ){1+(d−1)ρ}11t,{R}^{-1}\left(\rho )=\frac{1}{1-\rho }{I}_{d}-\frac{\rho }{\left(1-\rho )\left\{1+\left(d-1)\rho \right\}}{{\bf{11}}}^{t},where M2=diag(0,1,…,1,0){M}_{2}={\rm{diag}}\left(0,1,\ldots ,1,0)and M1{M}_{1}is a tri-diagonal matrix with 0 on the main diagonal and 1 on the upper and lower diagonals.(2)AR-1 structure: Under this structure, the (i,ji,j)th element of RRis given by ρ∣i−j∣{\rho }^{| i-j| }, with ρ∈(−1,1)\rho \in \left(-1,1). By [4], the inverse of RRas a function of ρ\rho is given by(2.7)R−1(ρ)=11−ρ2(Id−ρ2M2−ρM1),{R}^{-1}\left(\rho )=\frac{1}{1-{\rho }^{2}}\left({I}_{d}-{\rho }^{2}{M}_{2}-\rho {M}_{1}),where M2=diag(0,1,…,1,0){M}_{2}={\rm{diag}}\left(0,1,\ldots ,1,0)and M1{M}_{1}is tri-diagonal matrix with 0 on the main diagonal and 1 on the upper and lower diagonals. Gaussian copula with different correlation matrices by taking different values of ρ\rho are given in Figure 1.Figure 1Gaussian copula density functions with different values of ρ\rho : (a) ρ=0.9\rho =0.9and (b) ρ=−0.25\rho =-0.25.Several researchers [7,30,32] discussed the identifiability issues of the model parameters for multivariate discrete outcomes using copulas in the literature. Since the copula is not unique for multivariate discrete responses, [7] discussed that the model parameters are not identifiable uniquely. [30] diminished the identification problem of copulas using extensive simulation studies by considering several possibilities of generating data in the presence of continuous regressors. [32] argued that the continuous covariates widen the range of distribution function of the responses, ensuring the uniqueness of the copula.3Multivariate discrete probability mass function using a copulaHere, we discuss using a copula to construct a form of multivariate discrete probability mass function of a correlated dd-dimensional random vector from its marginals. Let Y=[Y1,Y2,…,Yd]T{\bf{Y}}={\left[{Y}_{1},{Y}_{2},\ldots ,{Y}_{d}]}^{T}be discrete random vectors with respective marginal cumulative distribution functions F1(y1∣θ1),F2(y2∣θ2),…,Fd(yd∣θd){F}_{1}({y}_{1}| {{\boldsymbol{\theta }}}_{1}),{F}_{2}({y}_{2}| {{\boldsymbol{\theta }}}_{2}),\ldots ,{F}_{d}({y}_{d}| {{\boldsymbol{\theta }}}_{d}). Then their joint probability mass function P(Y=y)P\left({\bf{Y}}={\bf{y}})formed using a copula CCis given by (3.1)fY∣C(y∣Ψ)=P(Y1=y1,Y2=y2,…,Yd=yd∣Ψ)=∑a1=12∑a2=12…∑ad=12(−1)a1+a2+…+adC(u1(a1),u2(a2),…,ud(ad)∣Ψ),\begin{array}{rcl}{f}_{{\bf{Y}}| C}({\bf{y}}| {\boldsymbol{\Psi }})& =& P\left({Y}_{1}={y}_{1},{Y}_{2}={y}_{2},\ldots ,{Y}_{d}={y}_{d}| {\boldsymbol{\Psi }})\\ & =& \mathop{\displaystyle \sum }\limits_{{a}_{1}=1}^{2}\mathop{\displaystyle \sum }\limits_{{a}_{2}=1}^{2}\ldots \mathop{\displaystyle \sum }\limits_{{a}_{d}=1}^{2}{\left(-1)}^{{a}_{1}+{a}_{2}+\ldots +{a}_{d}}C\left({u}_{1}^{\left({a}_{1})},{u}_{2}^{\left({a}_{2})},\ldots ,{u}_{d}^{\left({a}_{d})}| {\boldsymbol{\Psi }}),\end{array}where y∈D={0,1,2,…}d{\bf{y}}\in D={\left\{0,1,2,\ldots \right\}}^{d}, Ψ=[θ1T,θ2T,…,θdT,ΩT]T{\boldsymbol{\Psi }}={\left[{{\boldsymbol{\theta }}}_{1}^{T},{{\boldsymbol{\theta }}}_{2}^{T},\ldots ,{{\boldsymbol{\theta }}}_{d}^{T},{{\boldsymbol{\Omega }}}^{T}]}^{T}, Ω{\boldsymbol{\Omega }}is the copula parameters, uj(1)=Fj(yj∣θj){u}_{j}^{\left(1)}={F}_{j}({y}_{j}| {{\boldsymbol{\theta }}}_{j}), and uj(2)=Fj(yj−1∣θj){u}_{j}^{\left(2)}={F}_{j}({y}_{j}-1| {{\boldsymbol{\theta }}}_{j}), where j=1,…,dj=1,\ldots ,d. Here, we restrict ourselves to the case where F1(y1∣θ1),F2(y2∣θ2),…,Fd(yd∣θd){F}_{1}({y}_{1}| {{\boldsymbol{\theta }}}_{1}),{F}_{2}({y}_{2}| {{\boldsymbol{\theta }}}_{2}),\ldots ,{F}_{d}({y}_{d}| {{\boldsymbol{\theta }}}_{d})are negative binomial distribution functions with size parameter rj{r}_{j}, mean parameter μj{\mu }_{j}, and θj=[rj,μj]T{{\boldsymbol{\theta }}}_{j}={\left[{r}_{j},{\mu }_{j}]}^{T}. Recall that the probability mass function in Section 1 is given by (3.2)fNB(yj∣rj,μj)=Γ(yj+rj)Γ(rj)yj!rjrj+μjrjμjrj+μjyj,yj=0,1,….{f}_{{\rm{NB}}}({y}_{j}| {r}_{j},{\mu }_{j})=\frac{\Gamma ({y}_{j}+{r}_{j})}{\Gamma \left({r}_{j}){y}_{j}\hspace{0.1em}\text{&#x0021;}\hspace{0.1em}}{\left(\frac{{r}_{j}}{{r}_{j}+{\mu }_{j}}\right)}^{{r}_{j}}{\left(\frac{{\mu }_{j}}{{r}_{j}+{\mu }_{j}}\right)}^{{y}_{j}},\hspace{1em}{y}_{j}=0,1,\ldots .The corresponding distribution function is given by (3.3)FNB(y)=P(Yj≤y)=∑{yj:yj≤y}fNB(yj∣rj,μj),y∈R,{F}_{{\rm{NB}}}(y)=P\left({Y}_{j}\le y)=\sum _{\{{y}_{j}:{y}_{j}\le y\}}{f}_{{\rm{NB}}}({y}_{j}| {r}_{j},{\mu }_{j}),\hspace{1em}y\in {\mathbb{R}},where rj,μj>0{r}_{j},{\mu }_{j}\gt 0and yj∈{0,1,…}{y}_{j}\in \left\{0,1,\ldots \right\}. Bivariate negative binomial distributions with different correlation values and parameters are given in Figure 2.Figure 2Bivariate negative binomial distribution using a Gaussian copula with (a) r1=1.5,μ1=20,ρ=−0.30{r}_{1}=1.5,{\mu }_{1}=20,\rho =-0.30and (b) r2=1.5,μ2=20,ρ=0.80{r}_{2}=1.5,{\mu }_{2}=20,\rho =0.80.4Multivariate multiple-inflated negative binomial distribution functionHere, we describe the method of constructing a multivariate multiple-inflated negative binomial distribution function with the help of a multinomial distributed latent variable ZZ.Let Y{\bf{Y}}be a dd-dimensional random vector that follows a multivariate multiple-inflated negative binomial distribution with mm(a nonnegative integer) inflation points. We find a form of the distribution function FY(y){F}_{{\bf{Y}}}({\bf{y}})by considering a latent variable ZZwith the probability mass function given by (4.1)fZ(z)=pmifz=mpm−1ifz=m−1⋮p1ifz=11−qifz=0,{f}_{Z}\left(z)=\left\{\begin{array}{ll}{p}_{m}\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}z=m\\ {p}_{m-1}\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}z=m-1\\ \hspace{1em}\vdots \hspace{1.0em}\\ {p}_{1}\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}z=1\\ 1-q\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}z=0,\end{array}\right.where q=pm+pm−1+…+p1<1q={p}_{m}+{p}_{m-1}+\ldots +{p}_{1}\lt 1with pm,pm−1,…,p1∈(0,1){p}_{m},{p}_{m-1},\ldots ,{p}_{1}\in \left(0,1), and m is a nonnegative integer. Then the conditional probability mass function is given by (4.2)fY∣Z,Ψ(y∣z,Ψ)=1ify=kmandz=m1ify=km−1andz=m−1⋮1ify=k1andz=1fY∣C(y∣Ψ)ify∈Dandz=0,{f}_{{\bf{Y}}| Z,{\boldsymbol{\Psi }}}({\bf{y}}| z,{\boldsymbol{\Psi }})=\left\{\begin{array}{ll}1\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}{\bf{y}}={{\bf{k}}}_{m}\hspace{1em}\hspace{0.1em}\text{and}\hspace{0.1em}\hspace{1em}z=m\\ 1\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}{\bf{y}}={{\bf{k}}}_{m-1}\hspace{1em}\hspace{0.1em}\text{and}\hspace{0.1em}\hspace{1em}z=m-1\\ \vdots \hspace{1.0em}\\ 1\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}{\bf{y}}={{\bf{k}}}_{1}\hspace{1em}\hspace{0.1em}\text{and}\hspace{0.1em}\hspace{1em}z=1\\ {f}_{{\bf{Y}}| C}({\bf{y}}| {\boldsymbol{\Psi }})\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}{\bf{y}}\in D\hspace{1em}\hspace{0.1em}\text{and}\hspace{0.1em}\hspace{1em}z=0,\end{array}\right.where D={0,1,2,…}dD={\left\{0,1,2,\ldots \right\}}^{d}, ddis the dimension of Y{\bf{Y}}, fY∣C(⋅∣Ψ){f}_{{\bf{Y}}| C}\left(\cdot | {\boldsymbol{\Psi }})is a multivariate negative binomial probability mass function formed using a copula, Ψ=[r1,…,rd,μ1,…,μd,ΩT]T{\boldsymbol{\Psi }}={\left[{r}_{1},\ldots ,{r}_{d},{\mu }_{1},\ldots ,{\mu }_{d},{{\boldsymbol{\Omega }}}^{T}]}^{T}is a parameter list, Ω{\boldsymbol{\Omega }}is the copula parameters, and kl=[kl1,kl2,…,kld]T{{\bf{k}}}_{l}={\left[{k}_{l1},{k}_{l2},\ldots ,{k}_{ld}]}^{T}is a dd-dimensional inflation point for l=1,2,…,ml=1,2,\ldots ,m. By (4.1) and (4.2), we have the joint probability mass function: (4.3)fY,Z∣Ψ(y,z∣Ψ)=pmifz=mandy=kmpm−1ifz=m−1andy=km−1⋮p1ifz=1andy=k1(1−q)fY∣C(y∣Ψ)ifz=0andy∈D.{f}_{{\bf{Y}},Z| {\boldsymbol{\Psi }}}({\bf{y}},z| {\boldsymbol{\Psi }})=\left\{\begin{array}{ll}{p}_{m}\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}z=m\hspace{1em}\hspace{0.1em}\text{and}\hspace{0.1em}\hspace{1em}{\bf{y}}={{\bf{k}}}_{m}\\ {p}_{m-1}\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}z=m-1\hspace{1em}\hspace{0.1em}\text{and}\hspace{0.1em}\hspace{1em}{\bf{y}}={{\bf{k}}}_{m-1}\\ \hspace{1em}\vdots \hspace{1.0em}\\ {p}_{1}\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}z=1\hspace{1em}\hspace{0.1em}\text{and}\hspace{0.1em}\hspace{1em}{\bf{y}}={{\bf{k}}}_{1}\\ \left(1-q){f}_{{\bf{Y}}| C}({\bf{y}}| {\boldsymbol{\Psi }})\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}z=0\hspace{1em}\hspace{0.1em}\text{and}\hspace{0.1em}\hspace{1em}{\bf{y}}\in D.\end{array}\right.Finally, by (4.3), we have the multivariate multiple-inflated negative binomial probability mass function given by (4.4)fY∣Ψ(y∣Ψ)=pm+(1−q)fY∣C(y∣Ψ)ify=kmpm−1+(1−q)fY∣C(y∣Ψ)ify=km−1⋮p1+(1−q)fY∣C(y∣Ψ)ify=k1(1−q)fY∣C(y∣Ψ)otherwise,{f}_{{\bf{Y}}| {\boldsymbol{\Psi }}}({\bf{y}}| {\boldsymbol{\Psi }})=\left\{\begin{array}{ll}{p}_{m}+\left(1-q){f}_{{\bf{Y}}| C}({\bf{y}}| {\boldsymbol{\Psi }})\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}{\bf{y}}={{\bf{k}}}_{m}\\ {p}_{m-1}+\left(1-q){f}_{{\bf{Y}}| C}({\bf{y}}| {\boldsymbol{\Psi }})\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}{\bf{y}}={{\bf{k}}}_{m-1}\\ \hspace{1em}\vdots \hspace{1.0em}\\ {p}_{1}+\left(1-q){f}_{{\bf{Y}}| C}({\bf{y}}| {\boldsymbol{\Psi }})\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}{\bf{y}}={{\bf{k}}}_{1}\\ \left(1-q){f}_{{\bf{Y}}| C}({\bf{y}}| {\boldsymbol{\Psi }})\hspace{1.0em}& \hspace{0.1em}\text{otherwise}\hspace{0.1em},\end{array}\right.where y∈D{\bf{y}}\in D.Multivariate multiple-inflated negative binomial probability mass functions with different correlation values and parameters are given in Figure 3. The marginal probability mass function of the dd-dimensional random vector Y{\bf{Y}}can be found using (4.4): (4.5)fYj(yj)=∑y1…∑yj−1∑yj+1…∑ydfY(y∣Ψ),{f}_{{Y}_{j}}({y}_{j})=\sum _{{y}_{1}}\ldots \sum _{{y}_{j-1}}\sum _{{y}_{j+1}}\ldots \sum _{{y}_{d}}{f}_{{\bf{Y}}}({\bf{y}}| {\boldsymbol{\Psi }}),where yj{y}_{j}is the jjth component of y=[y1,y2,…,yd]T∈D{\bf{y}}={[{y}_{1},{y}_{2},\ldots ,{y}_{d}]}^{T}\in Dand Yj{Y}_{j}is the jjth components of the random vector Y{\bf{Y}}for j=1,2,…,dj=1,2,\ldots ,d. Therefore, the model simplifies to (4.6)fYj∣μj,rj(yj∣μj,rj)=pm+(1−q)fNB(yj∣μj,rj)ifyj=kmjpm−1+(1−q)fNB(yj∣μj,rj)ifyj=k(m−1)j⋮p1+(1−q)fNB(yj∣μj,rj)ifyi=k1j(1−q)fNB(yj∣μj,rj)otherwise,{f}_{{Y}_{j}| {\mu }_{j},{r}_{j}}({y}_{j}| {\mu }_{j},{r}_{j})=\left\{\begin{array}{ll}{p}_{m}+\left(1-q){f}_{{\rm{NB}}}({y}_{j}| {\mu }_{j},{r}_{j})\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}{y}_{j}={k}_{mj}\\ {p}_{m-1}+\left(1-q){f}_{{\rm{NB}}}({y}_{j}| {\mu }_{j},{r}_{j})\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}{y}_{j}={k}_{\left(m-1)j}\\ \hspace{1em}\vdots \hspace{1.0em}\\ {p}_{1}+\left(1-q){f}_{{\rm{NB}}}({y}_{j}| {\mu }_{j},{r}_{j})\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}{y}_{i}={k}_{1j}\\ \left(1-q){f}_{{\rm{NB}}}({y}_{j}| {\mu }_{j},{r}_{j})\hspace{1.0em}& \hspace{0.1em}\text{otherwise}\hspace{0.1em},\end{array}\right.where yj∈{0,1,2,…}{y}_{j}\in \left\{0,1,2,\ldots \right\}, kl=[kl1,kl2,…,kld]T{{\bf{k}}}_{l}={\left[{k}_{l1},{k}_{l2},\ldots ,{k}_{ld}]}^{T}is a dd-dimensional inflation point for l=1,2,…,ml=1,2,\ldots ,m, and fNB(yj∣μj,rj){f}_{{\rm{NB}}}({y}_{j}| {\mu }_{j},{r}_{j})is a negative binomial probability mass function with size rj{r}_{j}and mean μj{\mu }_{j}. From (4.6), the marginal expected value of the component Yj{Y}_{j}is given by (4.7)E(Yj)=∑yj∞yjfYj∣μj,rj(yj∣μj,rj)=0fYj∣μj,rj(0∣μj,rj)+k1jfYj∣μj,rj(k1j∣μj,rj)+…+kmjfYj∣μj,rj(kmj∣μj,rj)+∑yjyj≠0,k1j,…,kmjyjfYj(yj∣μj,rj)=k1jp1+…+kmjpj+(1−q)∑yj∞yjfNB(yj∣μj,rj)=k1jp1+…+kmjpj+(1−q)μj=λj,[say].\begin{array}{rcl}E\left({Y}_{j})& =& \mathop{\displaystyle \sum }\limits_{{y}_{j}}^{\infty }{y}_{j}{f}_{{Y}_{j}| {\mu }_{j},{r}_{j}}({y}_{j}| {\mu }_{j},{r}_{j})\\ & =& 0{f}_{{Y}_{j}| {\mu }_{j},{r}_{j}}\left(0| {\mu }_{j},{r}_{j})+{k}_{1j}{f}_{{Y}_{j}| {\mu }_{j},{r}_{j}}\left({k}_{1j}| {\mu }_{j},{r}_{j})+\ldots +{k}_{mj}{f}_{{Y}_{j}| {\mu }_{j},{r}_{j}}\left({k}_{mj}| {\mu }_{j},{r}_{j})+\displaystyle \sum _{\begin{array}{c}{y}_{j}\\ {y}_{j}\ne 0,{k}_{1j},\ldots ,{k}_{mj}\end{array}}{y}_{j}{f}_{{Y}_{j}}({y}_{j}| {\mu }_{j},{r}_{j})\\ & =& {k}_{1j}{p}_{1}+\ldots +{k}_{mj}{p}_{j}+\left(1-q)\mathop{\displaystyle \sum }\limits_{{y}_{j}}^{\infty }{y}_{j}{f}_{{\rm{NB}}}({y}_{j}| {\mu }_{j},{r}_{j})\\ & =& {k}_{1j}{p}_{1}+\ldots +{k}_{mj}{p}_{j}+\left(1-q){\mu }_{j}\\ & =& {\lambda }_{j},\hspace{0.33em}\left[{\rm{say}}].\end{array}As mentioned earlier, we have μj=λj−k1jp1−k2jp2−…−kmjpj1−q.{\mu }_{j}=\frac{{\lambda }_{j}-{k}_{1j}{p}_{1}-{k}_{2j}{p}_{2}-\ldots -{k}_{mj}{p}_{j}}{1-q}.In the next section, we provide a method of fitting a multivariate multiple inflated negative binomial regression model using data.Figure 3Bivariate multiple-inflated negative binomial distribution using a Gaussian copula with parameters (a) r1=1.5,μ1=20{r}_{1}=1.5,{\mu }_{1}=20and (b) r2=1.5,μ2=20{r}_{2}=1.5,{\mu }_{2}=20. (a) MMINB distribution with k1=(0,0){{\bf{k}}}_{1}=\left(0,0)and k2=(25,20){{\bf{k}}}_{2}=\left(25,20). (b) MMINB distribution with k1=(0,0){{\bf{k}}}_{1}=\left(0,0), k2=(0,20){{\bf{k}}}_{2}=\left(0,20), k3=(20,0){{\bf{k}}}_{3}=\left(20,0), and k4=(20,20){{\bf{k}}}_{4}=\left(20,20).5Multivariate multiple-inflated negative binomial regression modelIn this section, we propose a regression method for multivariate multiple-inflated negative binomial dd-dimensional responses y∈Rd{\bf{y}}\in {{\mathbb{R}}}^{d}with ν\nu -dimensional covariates vector x∈Rν{\bf{x}}\in {{\mathbb{R}}}^{\nu }. Let yi=[yi1,yi2,…,yid]T{{\bf{y}}}_{i}={[{y}_{i1},{y}_{i2},\ldots ,{y}_{id}]}^{T}be a vector of response variables conditional on the iith observation xi=[xi1,xi2,…,xiν]T,{{\bf{x}}}_{i}={\left[{x}_{i1},{x}_{i2},\ldots ,{x}_{i\nu }]}^{T},i=1,2,…,ni=1,2,\ldots ,n. We assume that yi∼MMINB(r1,…,rd,μi1,…μid,pi1,…,pim){{\bf{y}}}_{i}\hspace{0.33em} \sim \hspace{0.33em}{\rm{MMINB}}\left({r}_{1},\ldots ,{r}_{d},{\mu }_{i1},\ldots {\mu }_{id},{p}_{i1},\ldots ,{p}_{im}), where μij{\mu }_{ij}and rj{r}_{j}are the marginal mean and size, respectively, of fNB(yij){f}_{{\rm{NB}}}({y}_{ij})and pil{p}_{il}is the probability of observing an inflated value at kl{{\bf{k}}}_{l}, l=1,2,…,ml=1,2,\ldots ,m.From equation (4.7), we obtain, E(Yij∣xi)=k1jpi1+k2jpi2+⋯+kmjpim+(1−qi)μij=λij.E\left({Y}_{ij}| {{\bf{x}}}_{i})={k}_{1j}{p}_{i1}+{k}_{2j}{p}_{i2}+\cdots +{k}_{mj}{p}_{im}+\left(1-{q}_{i}){\mu }_{ij}={\lambda }_{ij}.We note that only μij{\mu }_{ij}and pil{p}_{il}are conditional on xi{{\bf{x}}}_{i}with i=1,…,ni=1,\ldots ,n, j=1,…,dj=1,\ldots ,d, and l=1,…,ml=1,\ldots ,m. We use the log-link function for modeling (4.7) : (5.1)log(λi1)=η1(xi)=f1(xi)β1log(λi2)=η2(xi)=f2(xi)β2⋮log(λid)=ηd(xi)=fd(xi)βd,\begin{array}{rcl}\log \left({\lambda }_{i1})& =& {\eta }_{1}\left({{\bf{x}}}_{i})={f}_{1}\left({{\bf{x}}}_{i}){{\boldsymbol{\beta }}}_{1}\\ \log \left({\lambda }_{i2})& =& {\eta }_{2}\left({{\bf{x}}}_{i})={f}_{2}\left({{\bf{x}}}_{i}){{\boldsymbol{\beta }}}_{2}\\ \vdots \\ \log \left({\lambda }_{id})& =& {\eta }_{d}\left({{\bf{x}}}_{i})={f}_{d}\left({{\bf{x}}}_{i}){{\boldsymbol{\beta }}}_{d},\end{array}where ηj(xi){\eta }_{j}\left({{\bf{x}}}_{i})is the linear predictor, i.e., a polynomial function of xi{{\bf{x}}}_{i}with coefficient βj{{\boldsymbol{\beta }}}_{{\bf{j}}}, j=1,2,…,dj=1,2,\ldots ,d, fj(xi){f}_{j}\left({{\bf{x}}}_{i})is a vector-valued function and βj{{\boldsymbol{\beta }}}_{j}is an τj×1{\tau }_{j}\times 1dimensional vector. Note that τj{\tau }_{j}is a nonnegative integer representing the dimension of the vector βj{{\boldsymbol{\beta }}}_{j}. It follows that the marginal mean of fNB(yij){f}_{{\rm{NB}}}({y}_{ij})is modeled as follows: (5.2)μij=exp{ηj(xi)}−k1jpi1−⋯−kmjpim1−qi.{\mu }_{ij}=\frac{\exp \left\{{\eta }_{j}\left({{\bf{x}}}_{i})\right\}-{k}_{1j}{p}_{i1}-\cdots -{k}_{mj}{p}_{im}}{1-{q}_{i}}.In a similar fashion, we model the probability of inflation points using the multinomial logit model.We model the logarithmic odds of each pil{p}_{il}using a set of linear predictors: (5.3)lnpimpi0=lnP(Z=m)P(Z=0)=ξm(xi)=gm(xi)γmlnpi(m−1)pi0=lnP(Z=m−1)P(Z=0)=ξm−1(xi)=gm−1(xi)γm−1⋮lnpi1pi0=lnP(Z=1)P(Z=0)=ξ1(xi)=g1(xi)γ1,\begin{array}{rcl}\mathrm{ln}\left\{\phantom{\rule[-1.25em]{}{0ex}},\frac{{p}_{im}}{{p}_{i0}}\right\}& =& \mathrm{ln}\left\{\phantom{\rule[-1.25em]{}{0ex}},\frac{P\left(Z=m)}{P\left(Z=0)}\right\}={\xi }_{m}\left({{\bf{x}}}_{i})={g}_{m}\left({{\bf{x}}}_{i}){{\boldsymbol{\gamma }}}_{m}\\ \mathrm{ln}\left\{\phantom{\rule[-1.25em]{}{0ex}},\frac{{p}_{i\left(m-1)}}{{p}_{i0}}\right\}& =& \mathrm{ln}\left\{\phantom{\rule[-1.25em]{}{0ex}},\frac{P\left(Z=m-1)}{P\left(Z=0)}\right\}={\xi }_{m-1}\left({{\bf{x}}}_{i})={g}_{m-1}\left({{\bf{x}}}_{i}){{\boldsymbol{\gamma }}}_{m-1}\\ & \vdots & \\ \mathrm{ln}\left\{\phantom{\rule[-1.25em]{}{0ex}},\frac{{p}_{i1}}{{p}_{i0}}\right\}& =& \mathrm{ln}\left\{\phantom{\rule[-1.25em]{}{0ex}},\frac{P\left(Z=1)}{P\left(Z=0)}\right\}={\xi }_{1}\left({{\bf{x}}}_{i})={g}_{1}\left({{\bf{x}}}_{i}){{\boldsymbol{\gamma }}}_{1},\end{array}where ξl(xi){\xi }_{l}\left({{\bf{x}}}_{i})is the linear predictor, i.e., a polynomial function of xi{{\bf{x}}}_{i}with coefficient γl{{\boldsymbol{\gamma }}}_{{\bf{l}}}, l=1,2,…,ml=1,2,\ldots ,m, gl(xi){g}_{l}\left({{\bf{x}}}_{i})is a vector-valued function and γl{{\boldsymbol{\gamma }}}_{l}is an ψl×1{\psi }_{l}\times 1dimensional vector of covariates l=1,2,…,ml=1,2,\ldots ,mand pi0=pi1+pi2+⋯+pim{p}_{i0}={p}_{i1}+{p}_{i2}+\cdots +{p}_{im}. From (5.3), we obtain (5.4)pim=exp{gm(xi)γm}1+∑l=1mexp{gl(xi)γl}pi(m−1)=exp{gm−1(xi)γm−1}1+∑l=1mexp{gl(xi)γl}⋮pi1=exp{g1(xi)γ1}1+∑l=1mexp{gl(xi)γl}.\begin{array}{rcl}{p}_{im}& =& \frac{\exp \left\{{g}_{m}\left({{\bf{x}}}_{i}){{\boldsymbol{\gamma }}}_{m}\right\}}{1+{\displaystyle \sum }_{l=1}^{m}\exp \left\{{g}_{l}\left({{\bf{x}}}_{i}){{\boldsymbol{\gamma }}}_{l}\right\}}\\ {p}_{i\left(m-1)}& =& \frac{\exp \left\{{g}_{m-1}\left({{\bf{x}}}_{i}){{\boldsymbol{\gamma }}}_{m-1}\right\}}{1+{\displaystyle \sum }_{l=1}^{m}\exp \left\{{g}_{l}\left({{\bf{x}}}_{i}){{\boldsymbol{\gamma }}}_{l}\right\}}\\ & \vdots & \\ {p}_{i1}& =& \frac{\exp \left\{{g}_{1}\left({{\bf{x}}}_{i}){{\boldsymbol{\gamma }}}_{1}\right\}}{1+{\displaystyle \sum }_{l=1}^{m}\exp \left\{{g}_{l}\left({{\bf{x}}}_{i}){{\boldsymbol{\gamma }}}_{l}\right\}}.\end{array}Next, we provide a method of estimating unknown parameters of the proposed regression model.6Estimation of parametersWe estimate model parameters Θ=[r1,…,rd,ℬT,ΓT,ΩT]T{\boldsymbol{\Theta }}={\left[{r}_{1},\ldots ,{r}_{d},{{\mathcal{ {\mathcal B} }}}^{T},{\Gamma }^{T},{{\boldsymbol{\Omega }}}^{T}]}^{T}using the method of maximum likelihood estimation (MLE), where ℬ=[β1,β2,…,βd]T{\mathcal{ {\mathcal B} }}={\left[{{\boldsymbol{\beta }}}_{1},{{\boldsymbol{\beta }}}_{2},\ldots ,{{\boldsymbol{\beta }}}_{d}]}^{T}, Γ=[γ1,γ2,…,γm]T\Gamma ={\left[{{\boldsymbol{\gamma }}}_{1},{{\boldsymbol{\gamma }}}_{2},\ldots ,{{\boldsymbol{\gamma }}}_{m}]}^{T}and Ω{\boldsymbol{\Omega }}is the copula parameters. The likelihood function is given by (6.1)l(Θ∣X,y)=∏{yi∣yi=km}log{pim+(1−qi)fCΦ∣R(yi∣Θ)}×∏{yi∣yi=km−1}log{pi(m−1)+(1−qi)fCΦ∣R(yi∣Θ)}×…×∏{yi∣yi≠{km,km−1,…,k1}}log{(1−qi)fCΦ∣R(yi∣Θ)},\begin{array}{rcl}l\left({\boldsymbol{\Theta }}| {\bf{X}},{\bf{y}})& =& \displaystyle \prod _{\{{{\bf{y}}}_{i}| {{\bf{y}}}_{i}={{\bf{k}}}_{m}\}}\log \left\{{p}_{im}+\left(1-{q}_{i}){f}_{{C}_{\Phi }| R}({{\bf{y}}}_{i}| {\boldsymbol{\Theta }})\right\}\\ & & \times \displaystyle \prod _{\{{{\bf{y}}}_{i}| {{\bf{y}}}_{i}={{\bf{k}}}_{m-1}\}}\log \left\{{p}_{i\left(m-1)}+\left(1-{q}_{i}){f}_{{C}_{\Phi }| R}({{\bf{y}}}_{i}| {\boldsymbol{\Theta }})\right\}\\ & & \times \ldots \times \displaystyle \prod _{\{{{\bf{y}}}_{i}| {{\bf{y}}}_{i}\ne \left\{{{\bf{k}}}_{m},{{\bf{k}}}_{m-1},\ldots ,{{\bf{k}}}_{1}\right\}\}}\log \left\{\left(1-{q}_{i}){f}_{{C}_{\Phi }| R}({{\bf{y}}}_{i}| {\boldsymbol{\Theta }})\right\},\end{array}where X{\bf{X}}is a design matrix. It follows that the log-likelihood function is given by (6.2)l(Θ∣X,y)=∑{yi∣yi=km}log{pim+(1−qi)fCΦ∣R(yi∣Θ)}+∑{yi∣yi=km−1}log{pi(m−1)+(1−qi)fCΦ∣R(yi∣Θ)}+…+∑{yi∣yi≠{km,km−1,…,k1}}log{(1−qi)fCΦ∣R(yi∣Θ)}\begin{array}{rcl}l\left({\boldsymbol{\Theta }}| {\bf{X}},{\bf{y}})& =& \displaystyle \sum _{\{{{\bf{y}}}_{i}| {{\bf{y}}}_{i}={{\bf{k}}}_{m}\}}\log \left\{{p}_{im}+\left(1-{q}_{i}){f}_{{C}_{\Phi }| R}({{\bf{y}}}_{i}| {\boldsymbol{\Theta }})\right\}\\ & & +\displaystyle \sum _{\{{{\bf{y}}}_{i}| {{\bf{y}}}_{i}={{\bf{k}}}_{m-1}\}}\log \left\{{p}_{i\left(m-1)}+\left(1-{q}_{i}){f}_{{C}_{\Phi }| R}({{\bf{y}}}_{i}| {\boldsymbol{\Theta }})\right\}+\ldots \\ & & +\displaystyle \sum _{\{{{\bf{y}}}_{i}| {{\bf{y}}}_{i}\ne \left\{{{\bf{k}}}_{m},{{\bf{k}}}_{m-1},\ldots ,{{\bf{k}}}_{1}\right\}\}}\log \left\{\left(1-{q}_{i}){f}_{{C}_{\Phi }| R}({{\bf{y}}}_{i}| {\boldsymbol{\Theta }})\right\}\end{array}The maximum likelihood estimates are found by maximizing (6.2). A complicated expression such as the one mentioned earlier does not have a closed form solution for parameter estimates. Moreover, the number of estimated parameters grows quickly for high-dimensional cases with many covariates and inflation points. In addition, care must be taken to ensure that μij{\mu }_{ij}given in (5.2) remains positive. Therefore, we maximize l(Θ∣X,y)l\left({\boldsymbol{\Theta }}| {\bf{X}},{\bf{y}})subject to the constraint (6.3)hj(Θ∣xi)=exp{ηj(xi)}−k1jpi1−…−kmjpim≥0∀i,j.{h}_{j}\left({\boldsymbol{\Theta }}| {{\bf{x}}}_{i})=\exp \left\{{\eta }_{j}\left({{\bf{x}}}_{i})\right\}-{k}_{1j}{p}_{i1}\hspace{0.33em}-\ldots -{k}_{mj}{p}_{im}\ge 0\hspace{1em}\forall i,j.This optimization problem can be solved using most nonlinear programming optimization algorithms. Here, we optimize ℬ{\mathcal{ {\mathcal B} }}and Γ\Gamma both subject to the constraint given in (6.3) and r1,…,rd{r}_{1},\ldots ,{r}_{d}subject to the constraint that r1,…,rd>0{r}_{1},\ldots ,{r}_{d}\gt 0using augmented Lagrangian methods and constrained optimization by linear approximation (COBYLA) methods. This type of optimization can be done in R programming using the auglag function in the nloptr package, (see [33]). For more information regarding these two methods, see [6,21]. Due to a large number of parameters, we use an iterative maximization procedure similar to [11] and [28] for obtaining estimates: (1)Obtain initial size estimates r10,…,rd0{r}_{1}^{0},\ldots ,{r}_{d}^{0}by maximizing the marginal log-likelihood functions with respect to r1,…,rd{r}_{1},\ldots ,{r}_{d}.(2)Obtain initial estimates ℬ0=[β10,…,βd0]T{{\mathcal{ {\mathcal B} }}}^{0}={\left[{{\boldsymbol{\beta }}}_{{\bf{1}}}^{0},\ldots ,{{\boldsymbol{\beta }}}_{d}^{0}]}^{T}by fitting a negative binomial regression model to the marginal distributions.(3)Obtain initial estimates Γ0=[γ10,…,γm0]T{\Gamma }^{0}={\left[{{\boldsymbol{\gamma }}}_{{\bf{1}}}^{0},\ldots ,{{\boldsymbol{\gamma }}}_{m}^{0}]}^{T}by fitting a multinomial logistic model to Z given in equation (4.3).(4)Obtain Ω0{{\boldsymbol{\Omega }}}^{0}by maximizing l(Ω∣X,y,r10,…,rd0,ℬ0,Γ0)l\left({\boldsymbol{\Omega }}| {\bf{X}},{\bf{y}},{r}_{1}^{0},\ldots ,{r}_{d}^{0},{{\mathcal{ {\mathcal B} }}}^{0},{\Gamma }^{0}).(5)Having obtained initial estimates, obtain r1i,…,rdi{r}_{1}^{i},\ldots ,{r}_{d}^{i}by maximizing l(r1,…,rd∣X,y,Ωi−1,ℬi−1,Γi−1)l\left({r}_{1},\ldots ,{r}_{d}| {\bf{X}},{\bf{y}},{{\boldsymbol{\Omega }}}^{i-1},{{\mathcal{ {\mathcal B} }}}^{i-1},{\Gamma }^{i-1})(6)Similarly, obtain ℬi{{\mathcal{ {\mathcal B} }}}^{i}, Γi{\Gamma }^{i}, Ωi{{\boldsymbol{\Omega }}}^{i}by maximizing l(ℬ∣X,y,r1i−1,…,rdi−1,Ωi−1,Γi−1)l\left({\mathcal{ {\mathcal B} }}| {\bf{X}},{\bf{y}},{r}_{1}^{i-1},\ldots ,{r}_{d}^{i-1},{{\boldsymbol{\Omega }}}^{i-1},{\Gamma }^{i-1}), l(Γ∣X,y,r1i−1,…,rdi−1,Ωi−1,ℬi−1)l\left(\Gamma | {\bf{X}},{\bf{y}},{r}_{1}^{i-1},\ldots ,{r}_{d}^{i-1},{{\boldsymbol{\Omega }}}^{i-1},{{\mathcal{ {\mathcal B} }}}^{i-1}), and l(Ω∣X,y,r1i−1,…,rdi−1,ℬi−1,Γi−1)l\left({\boldsymbol{\Omega }}| {\bf{X}},{\bf{y}},{r}_{1}^{i-1},\ldots ,{r}_{d}^{i-1},{{\mathcal{ {\mathcal B} }}}^{i-1},{\Gamma }^{i-1}), respectively, for i=1,2,3,…i=1,2,3,\ldots .(7)Repeat steps 5 and 6 until some convergence criteria are met.Here, convergence is obtained if the maximum of the differences between the current and previous parameter estimates is less than some threshold ε\varepsilon . However, we note that convergence in Akaike Information Criteria (AIC) or Bayesian Information Criteria may be used as well. We may also use existing packages such as “optim” and “maxLik” packages in the R environment to find MLEs by maximizing the aforementioned log-likelihood function.7Numerical exampleHere, we provide numerical examples to illustrate the proposed methodologies using simulated data and apply the method to real data.7.1Simulation studiesHere, we consider 500 equidistance points of independent variables xxfrom 0 to 1 and simulate responses yi=[yi1,yi2]′∼MMINB(r1,r2,μi1,μi2,pi1,pi2,ρ){{\bf{y}}}_{i}=[{y}_{i1},{y}_{i2}]^{\prime} \hspace{0.33em} \sim \hspace{0.33em}{\rm{MMINB}}\left({r}_{1},{r}_{2},{\mu }_{i1},{\mu }_{i2},{p}_{i1},{p}_{i2},\rho )considering two inflation points [0,0]′\left[0,0]^{\prime} and [1,1]′\left[1,1]^{\prime} (i.e., m=2m=2) with probabilities pi1{p}_{i1}and pi2{p}_{i2}, respectively, and we consider the Gaussian copula with parameter Ω={ρ}{\boldsymbol{\Omega }}=\left\{\rho \right\}, where ρ\rho is the correlation coefficient of the marginal distributions . The size parameters are considered as r1=5,r2=10{r}_{1}=5,{r}_{2}=10, and the correlation parameter for Gaussian copula is considered as ρ=0.7\rho =0.7. The linear predictors for λik{\lambda }_{ik}are assumed as follows: η1(x)=1.4−2.0xη2(x)=−0.5+1.5x,\begin{array}{rcl}{\eta }_{1}\left(x)& =& 1.4-2.0x\\ {\eta }_{2}\left(x)& =& -0.5+1.5x,\end{array}and the link function is assumed as the log link function.Again, the linear predictor, for odds of pil{p}_{il}are as follows: ξ1(x)=−1.7+3.0xξ2(x)=−1.5+2.5x,\begin{array}{rcl}{\xi }_{1}\left(x)& =& -1.7+3.0x\\ {\xi }_{2}\left(x)& =& -1.5+2.5x,\end{array}and the link function is assumed as a log link function for generating the data. We simulate n=500n=500responses y{\bf{y}}at design points xx.Researchers usually fit the multivariate count responses assuming multivariate negative binomial distribution without considering the inflation points in the data. We may obtain better performance from the fitted model if we consider the presence of inflation points in the data. So, here, we consider two models, M1{M}_{1}and M2{M}_{2}, to fit the data, where the model M1{M}_{1}is the usual multivariate negative binomial regression model and the model M2{M}_{2}is the proposed multivariate multiple inflated negative binomial regression model, where in both models, we use Gaussian copula for the dependence structure. Models M1{M}_{1}and M2{M}_{2}are described as follows:Model M1{M}_{1}: Here, we fit the simulated data with a multivariate negative binomial regression model without considering the presence of inflation points in the data. The form of the multivariate negative binomial distribution is found using the Gaussian copula. The linear predictors are assumed as follows: η1(x)=β10+β11xη2(x)=β20+β21x,\begin{array}{rcl}{\eta }_{1}\left(x)& =& {\beta }_{10}+{\beta }_{11}x\\ {\eta }_{2}\left(x)& =& {\beta }_{20}+{\beta }_{21}x,\end{array}and the link function is considered a log link function for fitting the data. The unknown parameters of the model are estimated using the MLE. The estimated values of the parameters, along with their standard errors and pp-values, are reported in Table 1. The AIC value of the fitted model M1{M}_{1}is 2795.562.Table 1Parameter estimates, standard errors, pp-values, and the AIC value of the fitted model M1{M}_{1}using 500 data pointsParametersEstimatesSEpp-valueβ10=1.4{\beta }_{10}=1.41.52950.1049<0.0001\lt 0.0001β11=−2{\beta }_{11}=-2−2.2583-2.25830.2000<0.0001\lt 0.0001β20=−0.5{\beta }_{20}=-0.5−0.1078-0.10780.14080.0443β21=1.5{\beta }_{21}=1.50.59770.22070.0067r1=5{r}_{1}=50.76870.0783<0.0001\lt 0.0001r2=10{r}_{2}=100.57910.0603<0.0001\lt 0.0001ρ=0.7\rho =0.70.83350.0150<0.0001\lt 0.0001Log-likelihood−1390.781-1390.781AIC2795.562Model M2{M}_{2}: Here, we fit the simulated data with the proposed model, assuming the presence of inflation points in the data. Since the data contain two inflation points (0,0) and (1,0), we fit a multivariate doubly inflated negative binomial regression model using the data. The forms of the linear predictors are assumed as follows: η1(x)=β10+β11x,η2(x)=β20+β21x,\begin{array}{rcl}{\eta }_{1}\left(x)& =& {\beta }_{10}+{\beta }_{11}x,\\ {\eta }_{2}\left(x)& =& {\beta }_{20}+{\beta }_{21}x,\end{array}and ξ1(x)=γ10+γ11x,ξ2(x)=γ20+γ21x,\begin{array}{rcl}{\xi }_{1}\left(x)& =& {\gamma }_{10}+{\gamma }_{11}x,\\ {\xi }_{2}\left(x)& =& {\gamma }_{20}+{\gamma }_{21}x,\end{array}and the corresponding link functions are assumed as log link function and multivariate logistic link function, respectively, as described in Section 5. The estimated values of the parameters using the MLE and the AIC value of the fitted model M2{M}_{2}are given in Table 2.Table 2Parameter estimates, standard errors, pp-values, and the AIC value of the fitted model M2{M}_{2}using 500 simulated data pointsParametersEstimatesSEpp-valueβ10=1.4{\beta }_{10}=1.41.48600.0902<0.0001\lt 0.0001β11=−2{\beta }_{11}=-2−2.2178-2.21780.1838<0.0001\lt 0.0001β20=−0.5{\beta }_{20}=-0.5−0.2385-0.23850.10380.0215β21=1.5{\beta }_{21}=1.50.69880.23310.0027γ10=−1.7{\gamma }_{10}=-1.7−1.9139-1.91390.2509<0.0001\lt 0.0001γ11=3.0{\gamma }_{11}=3.04.16350.4403<0.0001\lt 0.0001γ20=−1.5{\gamma }_{20}=-1.5−1.2426-1.24260.2113<0.0001\lt 0.0001γ21=2.5{\gamma }_{21}=2.53.00980.3594<0.0001\lt 0.0001r1=5{r}_{1}=55.79181.56300.0002r2=10{r}_{2}=108.88435.71330.0119ρ=0.7\rho =0.70.68970.0446<0.0001\lt 0.0001Log-likelihood−1059.785-1059.785AIC2141.569From the simulation study, we observe that the AIC value of model M2{M}_{2}is smaller than the AIC value of model M1{M}_{1}. The difference of the deviances of models M1{M}_{1}and M2{M}_{2}is 661.992 with 4 degrees of freedom (pp-value <0.0001\lt 0.0001). This shows that we obtain a significantly better fit using model M2{M}_{2}compared to model M1{M}_{1}by considering the existence of the multiple inflated points in the data.7.2Application to work week hours for couplesHere, we apply the proposed method of fitting a multivariate multiple inflated negative binomial regression model to real data. We model the work week hours for husbands and wives from the cps91 data set found in the Wooldridge package in R. The data set comes from the May 1991 Current Population Survey with a sample size of n=5,634n=\hspace{0.1em}\text{5,634}\hspace{0.1em}and contains factors such as education, age, ethnicity, and income for husbands and wives. We consider Yi1=“hushrs”{Y}_{i1}={\mbox{''}hushrs\mbox{''}}(husband’s weekly hours) and Yi2=“hours”{Y}_{i2}=\hspace{0.1em}\text{&#x201C;}\hspace{0.1em}{hours}\hspace{0.0em}\text{&#x201D;}\hspace{0.1em}(wife’s weekly hours) given in the data and assume Yi=[Yi1,Yi2]T∼MMINB(r1,r2,μi1,μi2,pi1,…,pim,ρ){{\bf{Y}}}_{i}={\left[{Y}_{i1},{Y}_{i2}]}^{T}\hspace{0.33em} \sim \hspace{0.33em}{\rm{MMINB}}\left({r}_{1},{r}_{2},{\mu }_{i1},{\mu }_{i2},{p}_{i1},\ldots ,{p}_{im},\rho ). We model λi1{\lambda }_{i1}, λi2{\lambda }_{i2}and pi1,pi2,…,pim{p}_{i1},{p}_{i2},\ldots ,{p}_{im}using the covariates, x1=“huseduc”{x}_{1}=\hspace{0.1em}\text{&#x201C;}\hspace{0.1em}{huseduc}\hspace{0.0em}\text{&#x201D;}\hspace{0.1em}, x2=“educ”{x}_{2}=\hspace{0.1em}\text{&#x201C;}\hspace{0.1em}{educ}\hspace{0.0em}\text{&#x201D;}\hspace{0.1em}, and x3=“expersq”{x}_{3}=\hspace{0.1em}\text{&#x201C;}\hspace{0.1em}{expersq}\hspace{0.0em}\text{&#x201D;}\hspace{0.1em}, which are corresponding to the husband’s years of education, wife’s years of education, and wife’s years of experience squared, respectively.We consider five different models to fit the data corresponding to five different number of inflation points m=0,1,…,4m=0,1,\ldots ,4. The inflation points are chosen according to their frequency of occurring given in Table 3. For all models, we use the Gaussian copula to find the form of the multivariate distribution function. The models are described as follows:Table 3Ten points with highest countsPointCount(40,0)\left(40,0)851(40,40)\left(40,40)781(0,0)\left(0,0)468(0,40)\left(0,40)214(50,0)\left(50,0)178(60,0)\left(60,0)156(50,40)\left(50,40)140(60,40)\left(60,40)93(45,0)\left(45,0)84(40,35)\left(40,35)79Model M(0){M}^{\left(0)}: Here, we consider the proposed model with m=0m=0, i.e., the usual multivariate negative binomial regression model without considering the presence of inflation points in the data. The linear predictors are assumed as follows: η1(x)=β10+β11x1+β12x2+β13x3η2(x)=β20+β21x1+β22x2+β23x3,\begin{array}{rcl}{\eta }_{1}\left({\bf{x}})& =& {\beta }_{10}+{\beta }_{11}{x}_{1}+{\beta }_{12}{x}_{2}+{\beta }_{13}{x}_{3}\\ {\eta }_{2}\left({\bf{x}})& =& {\beta }_{20}+{\beta }_{21}{x}_{1}+{\beta }_{22}{x}_{2}+{\beta }_{23}{x}_{3},\end{array}and the link function is considered a log link function for fitting the data. We estimate the parameters using MLE, and the estimated parameter values are presented in Table 4. The log-likelihood and AIC values are −46554.65-46554.65and 93131.3, respectively.Table 4Parameter estimates for m=0m=0ParameterEstimateStandard errorpp-valuer1{r}_{1}1.11340.0242<0.0001\lt 0.0001r2{r}_{2}0.26540.0058<0.0001\lt 0.0001ρ\rho 0.08410.0129<0.0001\lt 0.0001β10{\beta }_{10}3.61830.0128<0.0001\lt 0.0001β11{\beta }_{11}−0.0160-0.01600.01660.3380β12{\beta }_{12}0.10200.0162<0.0001\lt 0.0001β13{\beta }_{13}−0.1124-0.11240.0139<0.0001\lt 0.0001β20{\beta }_{20}3.00770.0260<0.0001\lt 0.0001β21{\beta }_{21}0.17860.0332<0.0001\lt 0.0001β22{\beta }_{22}−0.0432-0.04320.03350.1965β23{\beta }_{23}−0.0715-0.07150.02870.0127Log-likelihood−46,554.65AIC93,131.3Model M(1){M}^{\left(1)}: We fit the data with Model M1 assuming the presence of one inflation point, i.e., m=1m=1in the data and the inflation point is (40,0)\left(40,0)having the maximum number of occurrence among all observations. The linear predictors of the model are considered as follows: η1(x)=β10+β11x1+β12x2+β13x3,η2(x)=β20+β21x1+β22x2+β23x3,\begin{array}{rcl}{\eta }_{1}\left({\bf{x}})& =& {\beta }_{10}+{\beta }_{11}{x}_{1}+{\beta }_{12}{x}_{2}+{\beta }_{13}{x}_{3},\\ {\eta }_{2}\left({\bf{x}})& =& {\beta }_{20}+{\beta }_{21}{x}_{1}+{\beta }_{22}{x}_{2}+{\beta }_{23}{x}_{3},\end{array}and ξ1(x)=γ10+γ11x1+γ12x2+γ13x3,{\xi }_{1}\left({\bf{x}})={\gamma }_{10}+{\gamma }_{11}{x}_{1}+{\gamma }_{12}{x}_{2}+{\gamma }_{13}{x}_{3},and the corresponding link functions are assumed as log link function and multivariate logistic link function, respectively, as described in Section 5. The estimated parameter values using the MLE are given in Table 5. We observe that the log-likelihood value is −43644.01-43644.01, and the AIC value is 87318.02 for the fitted model.Table 5Parameter estimates for m=1m=1ParameterEstimateStandard errorpp-valuer1{r}_{1}0.87830.0204<0.0001\lt 0.0001r2{r}_{2}0.41260.0096<0.0001\lt 0.0001ρ\rho 0.18830.0135<0.0001\lt 0.0001β10{\beta }_{10}3.61560.0130<0.0001\lt 0.0001β11{\beta }_{11}−0.0482-0.04820.01530.0016β12{\beta }_{12}0.11620.0156<0.0001\lt 0.0001β13{\beta }_{13}−0.1026-0.10260.0137<0.0001\lt 0.0001β20{\beta }_{20}3.00790.0234<0.0001\lt 0.0001β21{\beta }_{21}0.18150.0305<0.0001\lt 0.0001β22{\beta }_{22}−0.0439-0.04390.03020.1464β23{\beta }_{23}0.00600.02620.8191γ10{\gamma }_{10}−1.7672-1.76720.0384<0.0001\lt 0.0001γ11{\gamma }_{11}−0.3118-0.31180.0433<0.0001\lt 0.0001γ12{\gamma }_{12}−0.0491-0.04910.04430.2681γ13{\gamma }_{13}−0.00003-0.000030.03050.9991Log-likelihood−43,644.01AIC87,318.02Model M(2){M}^{\left(2)}: This model fits the data considering presence of two inflation points (40,0)\left(40,0)and (40,40)\left(40,40)(i.e., m=2m=2) with linear predictors η1(x)=β10+β11x1+β12x2+β13x3,η2(x)=β20+β21x1+β22x2+β23x3,\begin{array}{rcl}{\eta }_{1}\left({\bf{x}})& =& {\beta }_{10}+{\beta }_{11}{x}_{1}+{\beta }_{12}{x}_{2}+{\beta }_{13}{x}_{3},\\ {\eta }_{2}\left({\bf{x}})& =& {\beta }_{20}+{\beta }_{21}{x}_{1}+{\beta }_{22}{x}_{2}+{\beta }_{23}{x}_{3},\end{array}and ξ1(x)=γ10+γ11x1+γ12x2+γ13x3,ξ2(x)=γ20+γ21x1+γ22x2+γ23x3,\begin{array}{rcl}{\xi }_{1}\left({\bf{x}})& =& {\gamma }_{10}+{\gamma }_{11}{x}_{1}+{\gamma }_{12}{x}_{2}+{\gamma }_{13}{x}_{3},\\ {\xi }_{2}\left({\bf{x}})& =& {\gamma }_{20}+{\gamma }_{21}{x}_{1}+{\gamma }_{22}{x}_{2}+{\gamma }_{23}{x}_{3},\end{array}and the corresponding link functions as log link function and multivariate logistic link function, respectively, as given in Section 5. The MLEs of the model parameters are given in Table 6. The log-likelihood and AIC values for the fitted model are reported as −37803.86-37803.86and 75645.72, respectively.Table 6Parameter estimates for m=2m=2ParameterEstimateStandard errorpp-valuer1{r}_{1}0.68230.0171<0.0001\lt 0.0001r2{r}_{2}0.31370.0081<0.0001\lt 0.0001ρ\rho 0.14660.0153<0.0001\lt 0.0001β10{\beta }_{10}3.61850.0134<0.0001\lt 0.0001β11{\beta }_{11}−0.0446-0.04460.01510.0033β12{\beta }_{12}0.10860.0151<0.0001\lt 0.0001β13{\beta }_{13}−0.1057-0.10570.0139<0.0001\lt 0.0001β20{\beta }_{20}3.01300.0221<0.0001\lt 0.0001β21{\beta }_{21}0.18570.0288<0.0001\lt 0.0001β22{\beta }_{22}−0.0436-0.04360.02880.1298β23{\beta }_{23}0.00980.02540.6989γ10{\gamma }_{10}−1.5908-1.59080.0390<0.0001\lt 0.0001γ11{\gamma }_{11}−0.2921-0.29210.0447<0.0001\lt 0.0001γ12{\gamma }_{12}0.02880.04460.5182γ13{\gamma }_{13}−0.0362-0.03620.03830.3454γ20{\gamma }_{20}−1.6585-1.65850.0399<0.0001\lt 0.0001γ21{\gamma }_{21}0.15580.05130.0024γ22{\gamma }_{22}−0.1196-0.11960.04870.0140γ23{\gamma }_{23}−0.13794-0.137940.04480.0021Log-likelihood−37803.86-37803.86AIC75645.72Model M(3){M}^{\left(3)}: Here, we consider the presence of three inflation points (40,0)\left(40,0), (40,40)\left(40,40), and (0,0)\left(0,0)(i.e., m=3m=3) in the data. The link functions are the same as model M(2){M}^{\left(2)}, and the linear predictors are as follows: η1(x)=β10+β11x1+β12x2+β13x3,η2(x)=β20+β21x1+β22x2+β23x3,\begin{array}{rcl}{\eta }_{1}\left({\bf{x}})& =& {\beta }_{10}+{\beta }_{11}{x}_{1}+{\beta }_{12}{x}_{2}+{\beta }_{13}{x}_{3},\\ {\eta }_{2}\left({\bf{x}})& =& {\beta }_{20}+{\beta }_{21}{x}_{1}+{\beta }_{22}{x}_{2}+{\beta }_{23}{x}_{3},\end{array}and ξ1(x)=γ10+γ11x1+γ12x2+γ13x3,ξ2(x)=γ20+γ21x1+γ22x2+γ23x3,ξ3(x)=γ30+γ31x1+γ32x2+γ33x3.\begin{array}{rcl}{\xi }_{1}\left({\bf{x}})& =& {\gamma }_{10}+{\gamma }_{11}{x}_{1}+{\gamma }_{12}{x}_{2}+{\gamma }_{13}{x}_{3},\\ {\xi }_{2}\left({\bf{x}})& =& {\gamma }_{20}+{\gamma }_{21}{x}_{1}+{\gamma }_{22}{x}_{2}+{\gamma }_{23}{x}_{3},\\ {\xi }_{3}\left({\bf{x}})& =& {\gamma }_{30}+{\gamma }_{31}{x}_{1}+{\gamma }_{32}{x}_{2}+{\gamma }_{33}{x}_{3}.\end{array}The estimated model parameters, AIC, and log-likelihood values are reported in Table 7. We see that the log-likelihood is −36672.42-36672.42, and the AIC value is 73390.84 for the fitted model.Table 7Parameter estimates for m=3m=3ParameterEstimateStandard errorpp-valuer1{r}_{1}1.57660.0436<0.0001\lt 0.0001r2{r}_{2}0.46010.0123<0.0001\lt 0.0001ρ\rho −0.2087-0.20870.0146<0.0001\lt 0.0001β10{\beta }_{10}3.63040.0103<0.0001\lt 0.0001β11{\beta }_{11}0.00010.01320.9936β12{\beta }_{12}0.08550.0131<0.0001\lt 0.0001β13{\beta }_{13}−0.1128-0.11280.0114<0.0001\lt 0.0001β20{\beta }_{20}3.02410.0205<0.0001\lt 0.0001β21{\beta }_{21}0.20460.0276<0.0001\lt 0.0001β22{\beta }_{22}−0.0180-0.01800.02790.5174β23{\beta }_{23}0.00410.02410.8658γ10{\gamma }_{10}−2.2440-2.24400.0571<0.0001\lt 0.0001γ11{\gamma }_{11}−0.1402-0.14020.05500.0108γ12{\gamma }_{12}−0.3246-0.32460.0519<0.0001\lt 0.0001γ13{\gamma }_{13}0.44120.0437<0.0001\lt 0.0001γ20{\gamma }_{20}−1.4689-1.46890.0396<0.0001\lt 0.0001γ21{\gamma }_{21}−0.3672-0.36720.0521<0.0001\lt 0.0001γ22{\gamma }_{22}−0.0001-0.00010.06350.9985γ23{\gamma }_{23}0.00090.03990.9817γ30{\gamma }_{30}−1.5335-1.53350.0404<0.0001\lt 0.0001γ31{\gamma }_{31}0.14850.05370.0057γ32{\gamma }_{32}−0.1768-0.17680.05220.0007γ33{\gamma }_{33}−0.0596-0.05960.04580.1939Log-likelihood−36672.42-36672.42AIC73,390.84Model M(4){M}^{\left(4)}: For this model, we consider the presence of four inflation points (40,0)\left(40,0), (40,40)\left(40,40), (0,0)\left(0,0), and (0,40)\left(0,40)(i.e., m=4m=4) in the data. The linear predictors are assumed as follows: η1(x)=β10+β11x1+β12x2+β13x3,η2(x)=β20+β21x1+β22x2+β23x3,\begin{array}{rcl}{\eta }_{1}\left({\bf{x}})& =& {\beta }_{10}+{\beta }_{11}{x}_{1}+{\beta }_{12}{x}_{2}+{\beta }_{13}{x}_{3},\\ {\eta }_{2}\left({\bf{x}})& =& {\beta }_{20}+{\beta }_{21}{x}_{1}+{\beta }_{22}{x}_{2}+{\beta }_{23}{x}_{3},\end{array}and ξ1(x)=γ10+γ11x1+γ12x2+γ13x3,ξ2(x)=γ20+γ21x1+γ22x2+γ23x3,ξ3(x)=γ30+γ31x1+γ32x2+γ33x3,ξ4(x)=γ40+γ41x1+γ42x2+γ43x3,\begin{array}{rcl}{\xi }_{1}\left({\bf{x}})& =& {\gamma }_{10}+{\gamma }_{11}{x}_{1}+{\gamma }_{12}{x}_{2}+{\gamma }_{13}{x}_{3},\\ {\xi }_{2}\left({\bf{x}})& =& {\gamma }_{20}+{\gamma }_{21}{x}_{1}+{\gamma }_{22}{x}_{2}+{\gamma }_{23}{x}_{3},\\ {\xi }_{3}\left({\bf{x}})& =& {\gamma }_{30}+{\gamma }_{31}{x}_{1}+{\gamma }_{32}{x}_{2}+{\gamma }_{33}{x}_{3},\\ {\xi }_{4}\left({\bf{x}})& =& {\gamma }_{40}+{\gamma }_{41}{x}_{1}+{\gamma }_{42}{x}_{2}+{\gamma }_{43}{x}_{3},\end{array}and the corresponding link functions are assumed as log link function and multivariate logistic link function, respectively, as described in Section 5. The values of the log-likelihood and AIC are −34957.68-34957.68and 69969.36, respectively, and the estimated model parameters are reported in Table 8.Table 8Parameter estimates for m=4m=4ParameterEstimateStandard errorpp-valuer1{r}_{1}3.48840.1045<0.0001\lt 0.0001r2{r}_{2}0.41510.0116<0.0001\lt 0.0001ρ\rho −0.1276-0.12760.0153<0.0001\lt 0.0001β10{\beta }_{10}3.63410.0083<0.0001\lt 0.0001β11{\beta }_{11}0.00070.01030.9434β12{\beta }_{12}0.07790.0103<0.0001\lt 0.0001β13{\beta }_{13}−0.0922-0.09220.0093<0.0001\lt 0.0001β20{\beta }_{20}3.02560.0204<0.0001\lt 0.0001β21{\beta }_{21}0.17800.0267<0.0001\lt 0.0001β22{\beta }_{22}−0.0082-0.00820.02650.7567β23{\beta }_{23}0.02210.02430.3618γ10{\gamma }_{10}−2.1693-2.16930.0572<0.0001\lt 0.0001γ11{\gamma }_{11}−0.0052-0.00520.05620.9267γ12{\gamma }_{12}−0.4043-0.40430.0508<0.0001\lt 0.0001γ13{\gamma }_{13}0.46680.0434<0.0001\lt 0.0001γ20{\gamma }_{20}−1.4054-1.40540.0400<0.0001\lt 0.0001γ21{\gamma }_{21}−0.2712-0.27120.0445<0.0001\lt 0.0001γ22{\gamma }_{22}−0.1092-0.10920.04680.0197γ23{\gamma }_{23}0.07330.03990.0659γ30{\gamma }_{30}−1.4692-1.46920.0408<0.0001\lt 0.0001γ31{\gamma }_{31}0.19170.05340.0003γ32{\gamma }_{32}−0.2320-0.23200.0525<0.0001\lt 0.0001γ33{\gamma }_{33}−0.0072-0.00720.04660.8766γ40{\gamma }_{40}−2.8408-2.84080.0764<0.0001\lt 0.0001γ41{\gamma }_{41}0.14460.09150.1138γ42{\gamma }_{42}−0.3998-0.39980.0792<0.0001\lt 0.0001γ43{\gamma }_{43}0.39150.0656<0.0001\lt 0.0001Log-likelihood−34957.68-34957.68AIC69,969.36Note that the same set of covariates is used for the marginal means and inflation probabilities since we are primarily interested in the change in the log-likelihood value and AIC for different values of mm.Histograms of the bivariate distribution and marginal distributions are given in Figures 4 and 5, respectively. We observe that the model with m=4m=4provides lowest AIC value among all the models with m∈{0,1,2,3,4}m\in \left\{0,1,2,3,4\right\}. Also, according to the chi-square test, the model considering four inflation points in this data provides a significantly best fit than the other models with m∈{0,1,2,3}m\in \left\{0,1,2,3\right\}.Figure 4Histogram of hushrs and hours variables from the cps91 data set.Figure 5Histogram of marginal distributions with estimated density. (a) Hushrs variable. (b) Hours variable.8ConclusionWe have presented a regression model in which the response vector y=[y1,y2,…]{\bf{y}}=[{y}_{1},{y}_{2},\ldots ]is assumed to follow a multiple-inflated negative binomial distribution. Copula methods were introduced to model the dependence structure of y{\bf{y}}through a correlation matrix RR. Due to the complexity of the likelihood function, an iterative procedure has been proposed for estimating model parameters. First, we consider simulated data and fit the proposed multivariate multiple inflated negative binomial model and the usual multivariate negative binomial regression model using Gaussian copula without considering inflation points in the data. We observe that the proposed model considering the presence of inflation points in the data, provides a significantly better fit than the usual model without considering the same. Next, five models accounting for mminflation points were fit to a current population survey data set, with m=0,…4m=0,\ldots 4. Accounting for the multiple-inflated structure generally reduced the AIC and log-likelihood value.One problem of the model is the difficulty of estimating parameters, despite the proposed iterative procedure. In particular, derivative-based optimization methods often fail at some point during the stepwise estimation process, and hence, the use of (COBYLA) methods mentioned in Section 6. With this in mind, as well as the constraints imposed in Section 6, there is no reason to expect parameter estimates to be global maxima.By using the proposed methodologies, one can choose any copula to find the multivariate distribution function from marginals. We used two-dimensional copulas for our numerical examples and found expected results by estimating parameters using the proposed algorithms. However, the estimation of model parameters becomes cumbersome as the response dimension increases. To address this problem, we may use vine copulas and other existing concepts in the literature to have computationally stable estimates from the model. One of the future directions may be fitting the model using different data sets and improving the estimation process. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Dependence Modeling de Gruyter

Multiple inflated negative binomial regression for correlated multivariate count data

Loading next page...
 
/lp/de-gruyter/multiple-inflated-negative-binomial-regression-for-correlated-2bpD4BZ52t
Publisher
de Gruyter
Copyright
© 2022 Joseph Mathews et al., published by De Gruyter
ISSN
2300-2298
eISSN
2300-2298
DOI
10.1515/demo-2022-0149
Publisher site
See Article on Publisher Site

Abstract

1IntroductionGeneralized linear models (GLMs) have become popular for modeling count data in which the underlying distribution of the response variables is assumed Poisson or negative binomial distribution. In cases where the equidispersion assumption of the Poisson model is violated, the negative binomial distribution with size r(>0)r\left(\gt 0)and mean μ(>0)\mu \left(\gt 0)given by (1.1)fNB(y∣r,μ)=Γ(y+r)Γ(r)y!rr+μrμr+μy,y=0,1,…{f}_{{\rm{NB}}}(y| r,\mu )=\frac{\Gamma (y+r)}{\Gamma \left(r)y\hspace{0.1em}\text{&#x0021;}\hspace{0.1em}}{\left(\frac{r}{r+\mu }\right)}^{r}{\left(\frac{\mu }{r+\mu }\right)}^{y},\hspace{1em}y=0,1,\ldots is often preferred. From [10], the variance of YYis given by (1.2)Var(Y)=μ+μ2r.\hspace{0.1em}\text{Var}\hspace{0.1em}\left(Y)=\mu +\frac{{\mu }^{2}}{r}.See [10] for a full derivation of (1.1) using a Poisson-gamma mixture. Notice that as r→∞r\to \infty , the distribution reduces to the Poisson in which E(Y)=Var(Y)=μE\left(Y)=\hspace{0.1em}\text{Var}\hspace{0.1em}\left(Y)=\mu . Thus, the dispersion parameter α=1/r\alpha =1\hspace{0.1em}\text{/}\hspace{0.1em}rallows one to directly model overdispersed data, unlike the traditional Poisson regression model. [3] showed that (1.1) is often referred to as the negative binomial-II (NB-II) model due to the quadratic term in (1.2). In general, negative binomial-P (NB-P) models have been given in the literature, such as in [8]. However, they are not explored in this article. See [10] for more on traditional negative binomial regression models.In certain cases, count data contain excessive zero-valued counts and violate both the Poisson and negative binomial distribution assumptions. Both hurdle and zero-inflation Poisson and negative binomial distribution models such as in [9,14,18] have been proposed to solve this problem. Essentially, both models assume zero-valued counts are generated from a second process. However, unlike hurdle models, zero-inflation models assume that zeroes are generated from both a second process and the underlying data generating process. For example, consider the work-week hours for males/females generated from a distribution with nonnegative integer-valued support (discussed in Section 7.2). Then there are two types of zero-valued counts: those unemployed and those unemployed due to a disability. That is, one type cannot work, while the other can, though both counts are zero-valued. Zero-inflation models allow one to model these situations through a binomial-count mixture model. More recently, doubly inflated models have been proposed by [23,24] for the Poisson distribution for cases when there are an excessive amount of zero-valued and kk-valued counts, where kkis a positive integer, such that the distributional assumptions of zero-inflation models also are not met. Mathews et al. [16] provided a method of finding doubly inflated multivariate negative binomial distributions using Gaussian copula. However, they did not consider the possibility of having multiple inflation points in the data. Here, we provide a method of finding multivariate distribution functions allowing arbitrary inflation points in the data using a copula. In this case, the standard zero-inflated negative binomial distribution is generalized to a multinomial-negative binomial mixture distribution. We also propose a regression method for correlated multiple inflated multivariate count data.Care must be taken in choosing an appropriate form of the multivariate negative binomial distribution. Numerous bivariate forms of the negative binomial distribution have been proposed in the literature. These are often constructed using a product of Poisson mass functions paired with a mixing distribution, such as the gamma distribution or log-normal model [1,13]. For example, [15] derived the bivariate negative binomial distribution from the probability mass functions of Y1∼Pois(αθ){Y}_{1}\hspace{0.33em} \sim \hspace{0.33em}{\rm{Pois}}\left(\alpha \theta ), Y2∼Pois(βθ){Y}_{2}\hspace{0.33em} \sim \hspace{0.33em}{\rm{Pois}}\left(\beta \theta ), where α>0\alpha \gt 0, β>0\beta \gt 0, θ>0\theta \gt 0, and a gamma mixing distribution with shape parameter rrand scale parameter λ=1\lambda =1: (1.3)∫0∞(αθ)y1e−αθy1!(βθ)y2e−βθy2!θr−1e−θΓ(r)dθ.\underset{0}{\overset{\infty }{\int }}\frac{{\left(\alpha \theta )}^{{y}_{1}}{e}^{-\alpha \theta }}{{y}_{1}\hspace{0.1em}\text{&#x0021;}\hspace{0.1em}}\frac{{\left(\beta \theta )}^{{y}_{2}}{e}^{-\beta \theta }}{{y}_{2}\hspace{0.1em}\text{&#x0021;}\hspace{0.1em}}\frac{{\theta }^{r-1}{e}^{-\theta }}{\Gamma \left(r)}{\rm{d}}\theta .Unfortunately, the aforementioned model only allows for nonnegative correlation values. Efforts have been made by [5,17] and others to relax this restriction. One possible solution is found by forming the multivariate negative binomial distribution through a copula. Indeed, copula theory has been growing in popularity for modeling count data, for example, [19] and [27]. A copula is a multivariate probability distribution with uniform marginal distributions and dependence parameter Ω{\boldsymbol{\Omega }}. A fundamental result in copula theory is Sklar’s theorem, which allows one to form a multivariate probability distribution through a copula with only knowledge of the form of the marginal distributions. In this article, we form a multivariate negative binomial distribution using a copula, in which Ω{\boldsymbol{\Omega }}represents the copula parameters.This article is organized as follows: in Section 2, we review the copula function, Sklar’s theorem, and the Gaussian copula; in Section 3, we review multivariate discrete distributions derived using copula functions and the marginal negative binomial distribution; in Section 4, we derive the multivariate multipleinflated negative binomial (MMINB) model using the copula function; in Section 5, we extend the results of Section 4 to a regression setting; in Section 6, we provide an iterative procedure for estimating the MMINB regression model parameters using maximum likelihood estimation (MLE); in Section 7.1, we consider simulated data and fit the data with proposed MMINB model and the usual multivariate negative binomial regression model using Gaussian copula; in Section 7.2, we apply the model to the cps91 dataset from the Wooldridge package in R, see [31] and [25]. We end with concluding remarks in Section 8.2Copulas and Sklar’s theoremThis section describes how the copula function may be used to form a multivariate probability mass function. We refer to [12] and [29] for more information on copulas than presented here. A copula is a multivariate probability distribution with uniform marginal distributions on the interval [0,1]\left[0,1].Definition 1A function C:[0,1]d→[0,1]C:{\left[0,1]}^{d}\to \left[0,1]is a dd-dimensional copula if the following are true: (1)C(1,…,yj,…,1)=yj,∀j=1,2,…,dC\left(1,\ldots ,{y}_{j},\ldots ,1)={y}_{j},\hspace{1em}\forall j=1,2,\ldots ,dwith (1,…,yj,…,1)∈[0,1]d\left(1,\ldots ,{y}_{j},\ldots ,1)\in {\left[0,1]}^{d},(2)C(y1,y2,…,yd)=0C({y}_{1},{y}_{2},\ldots ,{y}_{d})=0for (y1,y2,…,yd)∈[0,1]d({y}_{1},{y}_{2},\ldots ,{y}_{d})\in {\left[0,1]}^{d}with at least one of yj=0{y}_{j}=0for j=1,2,…,dj=1,2,\ldots ,d,(3)CCis dd-nondecreasing, i.e., for any u1=(u1(1),u2(1),…,ud(1)){{\bf{u}}}_{1}=\left({u}_{1}^{\left(1)},{u}_{2}^{\left(1)},\ldots ,{u}_{d}^{\left(1)}), u2=(u1(2),u2(2),…,ud(2)){{\bf{u}}}_{2}=\left({u}_{1}^{\left(2)},{u}_{2}^{\left(2)},\ldots ,{u}_{d}^{\left(2)})belong to [0,1]d{\left[0,1]}^{d}with uj(1)≤uj(2){u}_{j}^{\left(1)}\le {u}_{j}^{\left(2)}, for all j=1,2,…,dj=1,2,\ldots ,d:(2.1)∑a1=12∑a2=12…∑ad=12(−1)a1+a2+…+adC(u1(a1),u2(a2),…,ud(ad))≥0.\mathop{\sum }\limits_{{a}_{1}=1}^{2}\mathop{\sum }\limits_{{a}_{2}=1}^{2}\ldots \mathop{\sum }\limits_{{a}_{d}=1}^{2}{\left(-1)}^{{a}_{1}+{a}_{2}+\ldots +{a}_{d}}C\left({u}_{1}^{\left({a}_{1})},{u}_{2}^{\left({a}_{2})},\ldots ,{u}_{d}^{\left({a}_{d})})\ge 0.Sklar’s theorem, given in [26], is a fundamental result in copula theory as it allows one to form a joint distribution using known marginal distributions.Theorem 1Let Y1,Y2,…,Yd{Y}_{1},{Y}_{2},\ldots ,{Y}_{d}be random variables with marginal distribution functions F1,F2,…,Fd{F}_{1},{F}_{2},\ldots ,{F}_{d}and joint cumulative distribution function F, then the following are true: (1)There exists a d-dimensional copula C such that for all y1,y2,…,yd∈R{y}_{1},{y}_{2},\ldots ,{y}_{d}\in {\mathbb{R}}(2.2)F(y1,y2,…,yd)=C(F1(y1),F2(y2),…,Fd(yd)).F({y}_{1},{y}_{2},\ldots ,{y}_{d})=C\left({F}_{1}({y}_{1}),{F}_{2}({y}_{2}),\ldots ,{F}_{d}({y}_{d})).(2)If Y1,Y2,…,Yd{Y}_{1},{Y}_{2},\ldots ,{Y}_{d}are continuous then the copula C is unique. Otherwise, C can be uniquely determined on a d-dimensional rectangle Range(F1)×Range(F2)×…×Range(Fd){\rm{Range}}\left({F}_{1})\times {\rm{Range}}\left({F}_{2})\times \ldots \times {\rm{Range}}\left({F}_{d}).Copulas are popular because of their ability to model dependence among Y1,Y2,…,Yd{Y}_{1},{Y}_{2},\ldots ,{Y}_{d}. One popular copula is the Gaussian copula since this dependence is modeled through a correlation matrix RR.Definition 2The Gaussian copula is given by the function: (2.3)C(u1,u2,…,ud∣R)=ΦR(Φ−1(u1),Φ−1(u2),…,Φ−1(ud)),C\left({u}_{1},{u}_{2},\ldots ,{u}_{d}| R)={{\boldsymbol{\Phi }}}_{R}\left({\Phi }^{-1}\left({u}_{1}),{\Phi }^{-1}\left({u}_{2}),\ldots ,{\Phi }^{-1}\left({u}_{d})),where Φ−1{\Phi }^{-1}is the inverse cumulative distribution function of the standard Gaussian distribution and ΦR{{\boldsymbol{\Phi }}}_{R}is the joint cumulative distribution function of a standard multivariate Gaussian distribution with correlation matrix RR.Definition 3The Gaussian copula density [2] is given by (2.4)c(u1,u2,…,ud∣R)=1∣R∣exp−12UT×(R−1−Id)×U,c\left({u}_{1},{u}_{2},\ldots ,{u}_{d}| R)=\frac{1}{\sqrt{| R| }}\exp \left\{-\frac{1}{2}{{\bf{U}}}^{T}\times ({R}^{-1}-{I}_{d})\times {\bf{U}}\right\},where U=[Φ−1(u1),Φ−1(u2),…,Φ−1(ud)]T{\bf{U}}={{[}{\Phi }^{-1}\left({u}_{1}),{\Phi }^{-1}\left({u}_{2}),\ldots ,{\Phi }^{-1}\left({u}_{d})]}^{T}, Φ\Phi is the univariate standard Gaussian distribution, RRis the correlation matrix of the standard multivariate Gaussian distribution, and Id{I}_{d}be the d×dd\times didentity matrix.For high-dimensional problems, the Gaussian copula model quickly becomes intractable. One reason for this intractability is the large correlation matrix that must be estimated from the data. To overcome this issue, we impose a correlation structure of two types: (1)Equi-correlation structure: Under this structure, we assume RRas a function of ρ\rho given by(2.5)R(ρ)=ρ11T−(1−ρ)Id,R\left(\rho )=\rho {\bf{1}}{{\bf{1}}}^{T}-\left(1-\rho ){I}_{d},where Id{I}_{d}is a dd-dimensional identity matrix, ρ∈−1d−1,1\rho \in \left(-\frac{1}{d-1},1\right), and 1{\bf{1}}is a dd-dimensional column vector of ones. From [20], we have(2.6)R−1(ρ)=11−ρId−ρ(1−ρ){1+(d−1)ρ}11t,{R}^{-1}\left(\rho )=\frac{1}{1-\rho }{I}_{d}-\frac{\rho }{\left(1-\rho )\left\{1+\left(d-1)\rho \right\}}{{\bf{11}}}^{t},where M2=diag(0,1,…,1,0){M}_{2}={\rm{diag}}\left(0,1,\ldots ,1,0)and M1{M}_{1}is a tri-diagonal matrix with 0 on the main diagonal and 1 on the upper and lower diagonals.(2)AR-1 structure: Under this structure, the (i,ji,j)th element of RRis given by ρ∣i−j∣{\rho }^{| i-j| }, with ρ∈(−1,1)\rho \in \left(-1,1). By [4], the inverse of RRas a function of ρ\rho is given by(2.7)R−1(ρ)=11−ρ2(Id−ρ2M2−ρM1),{R}^{-1}\left(\rho )=\frac{1}{1-{\rho }^{2}}\left({I}_{d}-{\rho }^{2}{M}_{2}-\rho {M}_{1}),where M2=diag(0,1,…,1,0){M}_{2}={\rm{diag}}\left(0,1,\ldots ,1,0)and M1{M}_{1}is tri-diagonal matrix with 0 on the main diagonal and 1 on the upper and lower diagonals. Gaussian copula with different correlation matrices by taking different values of ρ\rho are given in Figure 1.Figure 1Gaussian copula density functions with different values of ρ\rho : (a) ρ=0.9\rho =0.9and (b) ρ=−0.25\rho =-0.25.Several researchers [7,30,32] discussed the identifiability issues of the model parameters for multivariate discrete outcomes using copulas in the literature. Since the copula is not unique for multivariate discrete responses, [7] discussed that the model parameters are not identifiable uniquely. [30] diminished the identification problem of copulas using extensive simulation studies by considering several possibilities of generating data in the presence of continuous regressors. [32] argued that the continuous covariates widen the range of distribution function of the responses, ensuring the uniqueness of the copula.3Multivariate discrete probability mass function using a copulaHere, we discuss using a copula to construct a form of multivariate discrete probability mass function of a correlated dd-dimensional random vector from its marginals. Let Y=[Y1,Y2,…,Yd]T{\bf{Y}}={\left[{Y}_{1},{Y}_{2},\ldots ,{Y}_{d}]}^{T}be discrete random vectors with respective marginal cumulative distribution functions F1(y1∣θ1),F2(y2∣θ2),…,Fd(yd∣θd){F}_{1}({y}_{1}| {{\boldsymbol{\theta }}}_{1}),{F}_{2}({y}_{2}| {{\boldsymbol{\theta }}}_{2}),\ldots ,{F}_{d}({y}_{d}| {{\boldsymbol{\theta }}}_{d}). Then their joint probability mass function P(Y=y)P\left({\bf{Y}}={\bf{y}})formed using a copula CCis given by (3.1)fY∣C(y∣Ψ)=P(Y1=y1,Y2=y2,…,Yd=yd∣Ψ)=∑a1=12∑a2=12…∑ad=12(−1)a1+a2+…+adC(u1(a1),u2(a2),…,ud(ad)∣Ψ),\begin{array}{rcl}{f}_{{\bf{Y}}| C}({\bf{y}}| {\boldsymbol{\Psi }})& =& P\left({Y}_{1}={y}_{1},{Y}_{2}={y}_{2},\ldots ,{Y}_{d}={y}_{d}| {\boldsymbol{\Psi }})\\ & =& \mathop{\displaystyle \sum }\limits_{{a}_{1}=1}^{2}\mathop{\displaystyle \sum }\limits_{{a}_{2}=1}^{2}\ldots \mathop{\displaystyle \sum }\limits_{{a}_{d}=1}^{2}{\left(-1)}^{{a}_{1}+{a}_{2}+\ldots +{a}_{d}}C\left({u}_{1}^{\left({a}_{1})},{u}_{2}^{\left({a}_{2})},\ldots ,{u}_{d}^{\left({a}_{d})}| {\boldsymbol{\Psi }}),\end{array}where y∈D={0,1,2,…}d{\bf{y}}\in D={\left\{0,1,2,\ldots \right\}}^{d}, Ψ=[θ1T,θ2T,…,θdT,ΩT]T{\boldsymbol{\Psi }}={\left[{{\boldsymbol{\theta }}}_{1}^{T},{{\boldsymbol{\theta }}}_{2}^{T},\ldots ,{{\boldsymbol{\theta }}}_{d}^{T},{{\boldsymbol{\Omega }}}^{T}]}^{T}, Ω{\boldsymbol{\Omega }}is the copula parameters, uj(1)=Fj(yj∣θj){u}_{j}^{\left(1)}={F}_{j}({y}_{j}| {{\boldsymbol{\theta }}}_{j}), and uj(2)=Fj(yj−1∣θj){u}_{j}^{\left(2)}={F}_{j}({y}_{j}-1| {{\boldsymbol{\theta }}}_{j}), where j=1,…,dj=1,\ldots ,d. Here, we restrict ourselves to the case where F1(y1∣θ1),F2(y2∣θ2),…,Fd(yd∣θd){F}_{1}({y}_{1}| {{\boldsymbol{\theta }}}_{1}),{F}_{2}({y}_{2}| {{\boldsymbol{\theta }}}_{2}),\ldots ,{F}_{d}({y}_{d}| {{\boldsymbol{\theta }}}_{d})are negative binomial distribution functions with size parameter rj{r}_{j}, mean parameter μj{\mu }_{j}, and θj=[rj,μj]T{{\boldsymbol{\theta }}}_{j}={\left[{r}_{j},{\mu }_{j}]}^{T}. Recall that the probability mass function in Section 1 is given by (3.2)fNB(yj∣rj,μj)=Γ(yj+rj)Γ(rj)yj!rjrj+μjrjμjrj+μjyj,yj=0,1,….{f}_{{\rm{NB}}}({y}_{j}| {r}_{j},{\mu }_{j})=\frac{\Gamma ({y}_{j}+{r}_{j})}{\Gamma \left({r}_{j}){y}_{j}\hspace{0.1em}\text{&#x0021;}\hspace{0.1em}}{\left(\frac{{r}_{j}}{{r}_{j}+{\mu }_{j}}\right)}^{{r}_{j}}{\left(\frac{{\mu }_{j}}{{r}_{j}+{\mu }_{j}}\right)}^{{y}_{j}},\hspace{1em}{y}_{j}=0,1,\ldots .The corresponding distribution function is given by (3.3)FNB(y)=P(Yj≤y)=∑{yj:yj≤y}fNB(yj∣rj,μj),y∈R,{F}_{{\rm{NB}}}(y)=P\left({Y}_{j}\le y)=\sum _{\{{y}_{j}:{y}_{j}\le y\}}{f}_{{\rm{NB}}}({y}_{j}| {r}_{j},{\mu }_{j}),\hspace{1em}y\in {\mathbb{R}},where rj,μj>0{r}_{j},{\mu }_{j}\gt 0and yj∈{0,1,…}{y}_{j}\in \left\{0,1,\ldots \right\}. Bivariate negative binomial distributions with different correlation values and parameters are given in Figure 2.Figure 2Bivariate negative binomial distribution using a Gaussian copula with (a) r1=1.5,μ1=20,ρ=−0.30{r}_{1}=1.5,{\mu }_{1}=20,\rho =-0.30and (b) r2=1.5,μ2=20,ρ=0.80{r}_{2}=1.5,{\mu }_{2}=20,\rho =0.80.4Multivariate multiple-inflated negative binomial distribution functionHere, we describe the method of constructing a multivariate multiple-inflated negative binomial distribution function with the help of a multinomial distributed latent variable ZZ.Let Y{\bf{Y}}be a dd-dimensional random vector that follows a multivariate multiple-inflated negative binomial distribution with mm(a nonnegative integer) inflation points. We find a form of the distribution function FY(y){F}_{{\bf{Y}}}({\bf{y}})by considering a latent variable ZZwith the probability mass function given by (4.1)fZ(z)=pmifz=mpm−1ifz=m−1⋮p1ifz=11−qifz=0,{f}_{Z}\left(z)=\left\{\begin{array}{ll}{p}_{m}\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}z=m\\ {p}_{m-1}\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}z=m-1\\ \hspace{1em}\vdots \hspace{1.0em}\\ {p}_{1}\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}z=1\\ 1-q\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}z=0,\end{array}\right.where q=pm+pm−1+…+p1<1q={p}_{m}+{p}_{m-1}+\ldots +{p}_{1}\lt 1with pm,pm−1,…,p1∈(0,1){p}_{m},{p}_{m-1},\ldots ,{p}_{1}\in \left(0,1), and m is a nonnegative integer. Then the conditional probability mass function is given by (4.2)fY∣Z,Ψ(y∣z,Ψ)=1ify=kmandz=m1ify=km−1andz=m−1⋮1ify=k1andz=1fY∣C(y∣Ψ)ify∈Dandz=0,{f}_{{\bf{Y}}| Z,{\boldsymbol{\Psi }}}({\bf{y}}| z,{\boldsymbol{\Psi }})=\left\{\begin{array}{ll}1\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}{\bf{y}}={{\bf{k}}}_{m}\hspace{1em}\hspace{0.1em}\text{and}\hspace{0.1em}\hspace{1em}z=m\\ 1\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}{\bf{y}}={{\bf{k}}}_{m-1}\hspace{1em}\hspace{0.1em}\text{and}\hspace{0.1em}\hspace{1em}z=m-1\\ \vdots \hspace{1.0em}\\ 1\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}{\bf{y}}={{\bf{k}}}_{1}\hspace{1em}\hspace{0.1em}\text{and}\hspace{0.1em}\hspace{1em}z=1\\ {f}_{{\bf{Y}}| C}({\bf{y}}| {\boldsymbol{\Psi }})\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}{\bf{y}}\in D\hspace{1em}\hspace{0.1em}\text{and}\hspace{0.1em}\hspace{1em}z=0,\end{array}\right.where D={0,1,2,…}dD={\left\{0,1,2,\ldots \right\}}^{d}, ddis the dimension of Y{\bf{Y}}, fY∣C(⋅∣Ψ){f}_{{\bf{Y}}| C}\left(\cdot | {\boldsymbol{\Psi }})is a multivariate negative binomial probability mass function formed using a copula, Ψ=[r1,…,rd,μ1,…,μd,ΩT]T{\boldsymbol{\Psi }}={\left[{r}_{1},\ldots ,{r}_{d},{\mu }_{1},\ldots ,{\mu }_{d},{{\boldsymbol{\Omega }}}^{T}]}^{T}is a parameter list, Ω{\boldsymbol{\Omega }}is the copula parameters, and kl=[kl1,kl2,…,kld]T{{\bf{k}}}_{l}={\left[{k}_{l1},{k}_{l2},\ldots ,{k}_{ld}]}^{T}is a dd-dimensional inflation point for l=1,2,…,ml=1,2,\ldots ,m. By (4.1) and (4.2), we have the joint probability mass function: (4.3)fY,Z∣Ψ(y,z∣Ψ)=pmifz=mandy=kmpm−1ifz=m−1andy=km−1⋮p1ifz=1andy=k1(1−q)fY∣C(y∣Ψ)ifz=0andy∈D.{f}_{{\bf{Y}},Z| {\boldsymbol{\Psi }}}({\bf{y}},z| {\boldsymbol{\Psi }})=\left\{\begin{array}{ll}{p}_{m}\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}z=m\hspace{1em}\hspace{0.1em}\text{and}\hspace{0.1em}\hspace{1em}{\bf{y}}={{\bf{k}}}_{m}\\ {p}_{m-1}\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}z=m-1\hspace{1em}\hspace{0.1em}\text{and}\hspace{0.1em}\hspace{1em}{\bf{y}}={{\bf{k}}}_{m-1}\\ \hspace{1em}\vdots \hspace{1.0em}\\ {p}_{1}\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}z=1\hspace{1em}\hspace{0.1em}\text{and}\hspace{0.1em}\hspace{1em}{\bf{y}}={{\bf{k}}}_{1}\\ \left(1-q){f}_{{\bf{Y}}| C}({\bf{y}}| {\boldsymbol{\Psi }})\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}z=0\hspace{1em}\hspace{0.1em}\text{and}\hspace{0.1em}\hspace{1em}{\bf{y}}\in D.\end{array}\right.Finally, by (4.3), we have the multivariate multiple-inflated negative binomial probability mass function given by (4.4)fY∣Ψ(y∣Ψ)=pm+(1−q)fY∣C(y∣Ψ)ify=kmpm−1+(1−q)fY∣C(y∣Ψ)ify=km−1⋮p1+(1−q)fY∣C(y∣Ψ)ify=k1(1−q)fY∣C(y∣Ψ)otherwise,{f}_{{\bf{Y}}| {\boldsymbol{\Psi }}}({\bf{y}}| {\boldsymbol{\Psi }})=\left\{\begin{array}{ll}{p}_{m}+\left(1-q){f}_{{\bf{Y}}| C}({\bf{y}}| {\boldsymbol{\Psi }})\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}{\bf{y}}={{\bf{k}}}_{m}\\ {p}_{m-1}+\left(1-q){f}_{{\bf{Y}}| C}({\bf{y}}| {\boldsymbol{\Psi }})\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}{\bf{y}}={{\bf{k}}}_{m-1}\\ \hspace{1em}\vdots \hspace{1.0em}\\ {p}_{1}+\left(1-q){f}_{{\bf{Y}}| C}({\bf{y}}| {\boldsymbol{\Psi }})\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}{\bf{y}}={{\bf{k}}}_{1}\\ \left(1-q){f}_{{\bf{Y}}| C}({\bf{y}}| {\boldsymbol{\Psi }})\hspace{1.0em}& \hspace{0.1em}\text{otherwise}\hspace{0.1em},\end{array}\right.where y∈D{\bf{y}}\in D.Multivariate multiple-inflated negative binomial probability mass functions with different correlation values and parameters are given in Figure 3. The marginal probability mass function of the dd-dimensional random vector Y{\bf{Y}}can be found using (4.4): (4.5)fYj(yj)=∑y1…∑yj−1∑yj+1…∑ydfY(y∣Ψ),{f}_{{Y}_{j}}({y}_{j})=\sum _{{y}_{1}}\ldots \sum _{{y}_{j-1}}\sum _{{y}_{j+1}}\ldots \sum _{{y}_{d}}{f}_{{\bf{Y}}}({\bf{y}}| {\boldsymbol{\Psi }}),where yj{y}_{j}is the jjth component of y=[y1,y2,…,yd]T∈D{\bf{y}}={[{y}_{1},{y}_{2},\ldots ,{y}_{d}]}^{T}\in Dand Yj{Y}_{j}is the jjth components of the random vector Y{\bf{Y}}for j=1,2,…,dj=1,2,\ldots ,d. Therefore, the model simplifies to (4.6)fYj∣μj,rj(yj∣μj,rj)=pm+(1−q)fNB(yj∣μj,rj)ifyj=kmjpm−1+(1−q)fNB(yj∣μj,rj)ifyj=k(m−1)j⋮p1+(1−q)fNB(yj∣μj,rj)ifyi=k1j(1−q)fNB(yj∣μj,rj)otherwise,{f}_{{Y}_{j}| {\mu }_{j},{r}_{j}}({y}_{j}| {\mu }_{j},{r}_{j})=\left\{\begin{array}{ll}{p}_{m}+\left(1-q){f}_{{\rm{NB}}}({y}_{j}| {\mu }_{j},{r}_{j})\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}{y}_{j}={k}_{mj}\\ {p}_{m-1}+\left(1-q){f}_{{\rm{NB}}}({y}_{j}| {\mu }_{j},{r}_{j})\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}{y}_{j}={k}_{\left(m-1)j}\\ \hspace{1em}\vdots \hspace{1.0em}\\ {p}_{1}+\left(1-q){f}_{{\rm{NB}}}({y}_{j}| {\mu }_{j},{r}_{j})\hspace{1.0em}& \hspace{0.1em}\text{if}\hspace{0.1em}\hspace{0.33em}{y}_{i}={k}_{1j}\\ \left(1-q){f}_{{\rm{NB}}}({y}_{j}| {\mu }_{j},{r}_{j})\hspace{1.0em}& \hspace{0.1em}\text{otherwise}\hspace{0.1em},\end{array}\right.where yj∈{0,1,2,…}{y}_{j}\in \left\{0,1,2,\ldots \right\}, kl=[kl1,kl2,…,kld]T{{\bf{k}}}_{l}={\left[{k}_{l1},{k}_{l2},\ldots ,{k}_{ld}]}^{T}is a dd-dimensional inflation point for l=1,2,…,ml=1,2,\ldots ,m, and fNB(yj∣μj,rj){f}_{{\rm{NB}}}({y}_{j}| {\mu }_{j},{r}_{j})is a negative binomial probability mass function with size rj{r}_{j}and mean μj{\mu }_{j}. From (4.6), the marginal expected value of the component Yj{Y}_{j}is given by (4.7)E(Yj)=∑yj∞yjfYj∣μj,rj(yj∣μj,rj)=0fYj∣μj,rj(0∣μj,rj)+k1jfYj∣μj,rj(k1j∣μj,rj)+…+kmjfYj∣μj,rj(kmj∣μj,rj)+∑yjyj≠0,k1j,…,kmjyjfYj(yj∣μj,rj)=k1jp1+…+kmjpj+(1−q)∑yj∞yjfNB(yj∣μj,rj)=k1jp1+…+kmjpj+(1−q)μj=λj,[say].\begin{array}{rcl}E\left({Y}_{j})& =& \mathop{\displaystyle \sum }\limits_{{y}_{j}}^{\infty }{y}_{j}{f}_{{Y}_{j}| {\mu }_{j},{r}_{j}}({y}_{j}| {\mu }_{j},{r}_{j})\\ & =& 0{f}_{{Y}_{j}| {\mu }_{j},{r}_{j}}\left(0| {\mu }_{j},{r}_{j})+{k}_{1j}{f}_{{Y}_{j}| {\mu }_{j},{r}_{j}}\left({k}_{1j}| {\mu }_{j},{r}_{j})+\ldots +{k}_{mj}{f}_{{Y}_{j}| {\mu }_{j},{r}_{j}}\left({k}_{mj}| {\mu }_{j},{r}_{j})+\displaystyle \sum _{\begin{array}{c}{y}_{j}\\ {y}_{j}\ne 0,{k}_{1j},\ldots ,{k}_{mj}\end{array}}{y}_{j}{f}_{{Y}_{j}}({y}_{j}| {\mu }_{j},{r}_{j})\\ & =& {k}_{1j}{p}_{1}+\ldots +{k}_{mj}{p}_{j}+\left(1-q)\mathop{\displaystyle \sum }\limits_{{y}_{j}}^{\infty }{y}_{j}{f}_{{\rm{NB}}}({y}_{j}| {\mu }_{j},{r}_{j})\\ & =& {k}_{1j}{p}_{1}+\ldots +{k}_{mj}{p}_{j}+\left(1-q){\mu }_{j}\\ & =& {\lambda }_{j},\hspace{0.33em}\left[{\rm{say}}].\end{array}As mentioned earlier, we have μj=λj−k1jp1−k2jp2−…−kmjpj1−q.{\mu }_{j}=\frac{{\lambda }_{j}-{k}_{1j}{p}_{1}-{k}_{2j}{p}_{2}-\ldots -{k}_{mj}{p}_{j}}{1-q}.In the next section, we provide a method of fitting a multivariate multiple inflated negative binomial regression model using data.Figure 3Bivariate multiple-inflated negative binomial distribution using a Gaussian copula with parameters (a) r1=1.5,μ1=20{r}_{1}=1.5,{\mu }_{1}=20and (b) r2=1.5,μ2=20{r}_{2}=1.5,{\mu }_{2}=20. (a) MMINB distribution with k1=(0,0){{\bf{k}}}_{1}=\left(0,0)and k2=(25,20){{\bf{k}}}_{2}=\left(25,20). (b) MMINB distribution with k1=(0,0){{\bf{k}}}_{1}=\left(0,0), k2=(0,20){{\bf{k}}}_{2}=\left(0,20), k3=(20,0){{\bf{k}}}_{3}=\left(20,0), and k4=(20,20){{\bf{k}}}_{4}=\left(20,20).5Multivariate multiple-inflated negative binomial regression modelIn this section, we propose a regression method for multivariate multiple-inflated negative binomial dd-dimensional responses y∈Rd{\bf{y}}\in {{\mathbb{R}}}^{d}with ν\nu -dimensional covariates vector x∈Rν{\bf{x}}\in {{\mathbb{R}}}^{\nu }. Let yi=[yi1,yi2,…,yid]T{{\bf{y}}}_{i}={[{y}_{i1},{y}_{i2},\ldots ,{y}_{id}]}^{T}be a vector of response variables conditional on the iith observation xi=[xi1,xi2,…,xiν]T,{{\bf{x}}}_{i}={\left[{x}_{i1},{x}_{i2},\ldots ,{x}_{i\nu }]}^{T},i=1,2,…,ni=1,2,\ldots ,n. We assume that yi∼MMINB(r1,…,rd,μi1,…μid,pi1,…,pim){{\bf{y}}}_{i}\hspace{0.33em} \sim \hspace{0.33em}{\rm{MMINB}}\left({r}_{1},\ldots ,{r}_{d},{\mu }_{i1},\ldots {\mu }_{id},{p}_{i1},\ldots ,{p}_{im}), where μij{\mu }_{ij}and rj{r}_{j}are the marginal mean and size, respectively, of fNB(yij){f}_{{\rm{NB}}}({y}_{ij})and pil{p}_{il}is the probability of observing an inflated value at kl{{\bf{k}}}_{l}, l=1,2,…,ml=1,2,\ldots ,m.From equation (4.7), we obtain, E(Yij∣xi)=k1jpi1+k2jpi2+⋯+kmjpim+(1−qi)μij=λij.E\left({Y}_{ij}| {{\bf{x}}}_{i})={k}_{1j}{p}_{i1}+{k}_{2j}{p}_{i2}+\cdots +{k}_{mj}{p}_{im}+\left(1-{q}_{i}){\mu }_{ij}={\lambda }_{ij}.We note that only μij{\mu }_{ij}and pil{p}_{il}are conditional on xi{{\bf{x}}}_{i}with i=1,…,ni=1,\ldots ,n, j=1,…,dj=1,\ldots ,d, and l=1,…,ml=1,\ldots ,m. We use the log-link function for modeling (4.7) : (5.1)log(λi1)=η1(xi)=f1(xi)β1log(λi2)=η2(xi)=f2(xi)β2⋮log(λid)=ηd(xi)=fd(xi)βd,\begin{array}{rcl}\log \left({\lambda }_{i1})& =& {\eta }_{1}\left({{\bf{x}}}_{i})={f}_{1}\left({{\bf{x}}}_{i}){{\boldsymbol{\beta }}}_{1}\\ \log \left({\lambda }_{i2})& =& {\eta }_{2}\left({{\bf{x}}}_{i})={f}_{2}\left({{\bf{x}}}_{i}){{\boldsymbol{\beta }}}_{2}\\ \vdots \\ \log \left({\lambda }_{id})& =& {\eta }_{d}\left({{\bf{x}}}_{i})={f}_{d}\left({{\bf{x}}}_{i}){{\boldsymbol{\beta }}}_{d},\end{array}where ηj(xi){\eta }_{j}\left({{\bf{x}}}_{i})is the linear predictor, i.e., a polynomial function of xi{{\bf{x}}}_{i}with coefficient βj{{\boldsymbol{\beta }}}_{{\bf{j}}}, j=1,2,…,dj=1,2,\ldots ,d, fj(xi){f}_{j}\left({{\bf{x}}}_{i})is a vector-valued function and βj{{\boldsymbol{\beta }}}_{j}is an τj×1{\tau }_{j}\times 1dimensional vector. Note that τj{\tau }_{j}is a nonnegative integer representing the dimension of the vector βj{{\boldsymbol{\beta }}}_{j}. It follows that the marginal mean of fNB(yij){f}_{{\rm{NB}}}({y}_{ij})is modeled as follows: (5.2)μij=exp{ηj(xi)}−k1jpi1−⋯−kmjpim1−qi.{\mu }_{ij}=\frac{\exp \left\{{\eta }_{j}\left({{\bf{x}}}_{i})\right\}-{k}_{1j}{p}_{i1}-\cdots -{k}_{mj}{p}_{im}}{1-{q}_{i}}.In a similar fashion, we model the probability of inflation points using the multinomial logit model.We model the logarithmic odds of each pil{p}_{il}using a set of linear predictors: (5.3)lnpimpi0=lnP(Z=m)P(Z=0)=ξm(xi)=gm(xi)γmlnpi(m−1)pi0=lnP(Z=m−1)P(Z=0)=ξm−1(xi)=gm−1(xi)γm−1⋮lnpi1pi0=lnP(Z=1)P(Z=0)=ξ1(xi)=g1(xi)γ1,\begin{array}{rcl}\mathrm{ln}\left\{\phantom{\rule[-1.25em]{}{0ex}},\frac{{p}_{im}}{{p}_{i0}}\right\}& =& \mathrm{ln}\left\{\phantom{\rule[-1.25em]{}{0ex}},\frac{P\left(Z=m)}{P\left(Z=0)}\right\}={\xi }_{m}\left({{\bf{x}}}_{i})={g}_{m}\left({{\bf{x}}}_{i}){{\boldsymbol{\gamma }}}_{m}\\ \mathrm{ln}\left\{\phantom{\rule[-1.25em]{}{0ex}},\frac{{p}_{i\left(m-1)}}{{p}_{i0}}\right\}& =& \mathrm{ln}\left\{\phantom{\rule[-1.25em]{}{0ex}},\frac{P\left(Z=m-1)}{P\left(Z=0)}\right\}={\xi }_{m-1}\left({{\bf{x}}}_{i})={g}_{m-1}\left({{\bf{x}}}_{i}){{\boldsymbol{\gamma }}}_{m-1}\\ & \vdots & \\ \mathrm{ln}\left\{\phantom{\rule[-1.25em]{}{0ex}},\frac{{p}_{i1}}{{p}_{i0}}\right\}& =& \mathrm{ln}\left\{\phantom{\rule[-1.25em]{}{0ex}},\frac{P\left(Z=1)}{P\left(Z=0)}\right\}={\xi }_{1}\left({{\bf{x}}}_{i})={g}_{1}\left({{\bf{x}}}_{i}){{\boldsymbol{\gamma }}}_{1},\end{array}where ξl(xi){\xi }_{l}\left({{\bf{x}}}_{i})is the linear predictor, i.e., a polynomial function of xi{{\bf{x}}}_{i}with coefficient γl{{\boldsymbol{\gamma }}}_{{\bf{l}}}, l=1,2,…,ml=1,2,\ldots ,m, gl(xi){g}_{l}\left({{\bf{x}}}_{i})is a vector-valued function and γl{{\boldsymbol{\gamma }}}_{l}is an ψl×1{\psi }_{l}\times 1dimensional vector of covariates l=1,2,…,ml=1,2,\ldots ,mand pi0=pi1+pi2+⋯+pim{p}_{i0}={p}_{i1}+{p}_{i2}+\cdots +{p}_{im}. From (5.3), we obtain (5.4)pim=exp{gm(xi)γm}1+∑l=1mexp{gl(xi)γl}pi(m−1)=exp{gm−1(xi)γm−1}1+∑l=1mexp{gl(xi)γl}⋮pi1=exp{g1(xi)γ1}1+∑l=1mexp{gl(xi)γl}.\begin{array}{rcl}{p}_{im}& =& \frac{\exp \left\{{g}_{m}\left({{\bf{x}}}_{i}){{\boldsymbol{\gamma }}}_{m}\right\}}{1+{\displaystyle \sum }_{l=1}^{m}\exp \left\{{g}_{l}\left({{\bf{x}}}_{i}){{\boldsymbol{\gamma }}}_{l}\right\}}\\ {p}_{i\left(m-1)}& =& \frac{\exp \left\{{g}_{m-1}\left({{\bf{x}}}_{i}){{\boldsymbol{\gamma }}}_{m-1}\right\}}{1+{\displaystyle \sum }_{l=1}^{m}\exp \left\{{g}_{l}\left({{\bf{x}}}_{i}){{\boldsymbol{\gamma }}}_{l}\right\}}\\ & \vdots & \\ {p}_{i1}& =& \frac{\exp \left\{{g}_{1}\left({{\bf{x}}}_{i}){{\boldsymbol{\gamma }}}_{1}\right\}}{1+{\displaystyle \sum }_{l=1}^{m}\exp \left\{{g}_{l}\left({{\bf{x}}}_{i}){{\boldsymbol{\gamma }}}_{l}\right\}}.\end{array}Next, we provide a method of estimating unknown parameters of the proposed regression model.6Estimation of parametersWe estimate model parameters Θ=[r1,…,rd,ℬT,ΓT,ΩT]T{\boldsymbol{\Theta }}={\left[{r}_{1},\ldots ,{r}_{d},{{\mathcal{ {\mathcal B} }}}^{T},{\Gamma }^{T},{{\boldsymbol{\Omega }}}^{T}]}^{T}using the method of maximum likelihood estimation (MLE), where ℬ=[β1,β2,…,βd]T{\mathcal{ {\mathcal B} }}={\left[{{\boldsymbol{\beta }}}_{1},{{\boldsymbol{\beta }}}_{2},\ldots ,{{\boldsymbol{\beta }}}_{d}]}^{T}, Γ=[γ1,γ2,…,γm]T\Gamma ={\left[{{\boldsymbol{\gamma }}}_{1},{{\boldsymbol{\gamma }}}_{2},\ldots ,{{\boldsymbol{\gamma }}}_{m}]}^{T}and Ω{\boldsymbol{\Omega }}is the copula parameters. The likelihood function is given by (6.1)l(Θ∣X,y)=∏{yi∣yi=km}log{pim+(1−qi)fCΦ∣R(yi∣Θ)}×∏{yi∣yi=km−1}log{pi(m−1)+(1−qi)fCΦ∣R(yi∣Θ)}×…×∏{yi∣yi≠{km,km−1,…,k1}}log{(1−qi)fCΦ∣R(yi∣Θ)},\begin{array}{rcl}l\left({\boldsymbol{\Theta }}| {\bf{X}},{\bf{y}})& =& \displaystyle \prod _{\{{{\bf{y}}}_{i}| {{\bf{y}}}_{i}={{\bf{k}}}_{m}\}}\log \left\{{p}_{im}+\left(1-{q}_{i}){f}_{{C}_{\Phi }| R}({{\bf{y}}}_{i}| {\boldsymbol{\Theta }})\right\}\\ & & \times \displaystyle \prod _{\{{{\bf{y}}}_{i}| {{\bf{y}}}_{i}={{\bf{k}}}_{m-1}\}}\log \left\{{p}_{i\left(m-1)}+\left(1-{q}_{i}){f}_{{C}_{\Phi }| R}({{\bf{y}}}_{i}| {\boldsymbol{\Theta }})\right\}\\ & & \times \ldots \times \displaystyle \prod _{\{{{\bf{y}}}_{i}| {{\bf{y}}}_{i}\ne \left\{{{\bf{k}}}_{m},{{\bf{k}}}_{m-1},\ldots ,{{\bf{k}}}_{1}\right\}\}}\log \left\{\left(1-{q}_{i}){f}_{{C}_{\Phi }| R}({{\bf{y}}}_{i}| {\boldsymbol{\Theta }})\right\},\end{array}where X{\bf{X}}is a design matrix. It follows that the log-likelihood function is given by (6.2)l(Θ∣X,y)=∑{yi∣yi=km}log{pim+(1−qi)fCΦ∣R(yi∣Θ)}+∑{yi∣yi=km−1}log{pi(m−1)+(1−qi)fCΦ∣R(yi∣Θ)}+…+∑{yi∣yi≠{km,km−1,…,k1}}log{(1−qi)fCΦ∣R(yi∣Θ)}\begin{array}{rcl}l\left({\boldsymbol{\Theta }}| {\bf{X}},{\bf{y}})& =& \displaystyle \sum _{\{{{\bf{y}}}_{i}| {{\bf{y}}}_{i}={{\bf{k}}}_{m}\}}\log \left\{{p}_{im}+\left(1-{q}_{i}){f}_{{C}_{\Phi }| R}({{\bf{y}}}_{i}| {\boldsymbol{\Theta }})\right\}\\ & & +\displaystyle \sum _{\{{{\bf{y}}}_{i}| {{\bf{y}}}_{i}={{\bf{k}}}_{m-1}\}}\log \left\{{p}_{i\left(m-1)}+\left(1-{q}_{i}){f}_{{C}_{\Phi }| R}({{\bf{y}}}_{i}| {\boldsymbol{\Theta }})\right\}+\ldots \\ & & +\displaystyle \sum _{\{{{\bf{y}}}_{i}| {{\bf{y}}}_{i}\ne \left\{{{\bf{k}}}_{m},{{\bf{k}}}_{m-1},\ldots ,{{\bf{k}}}_{1}\right\}\}}\log \left\{\left(1-{q}_{i}){f}_{{C}_{\Phi }| R}({{\bf{y}}}_{i}| {\boldsymbol{\Theta }})\right\}\end{array}The maximum likelihood estimates are found by maximizing (6.2). A complicated expression such as the one mentioned earlier does not have a closed form solution for parameter estimates. Moreover, the number of estimated parameters grows quickly for high-dimensional cases with many covariates and inflation points. In addition, care must be taken to ensure that μij{\mu }_{ij}given in (5.2) remains positive. Therefore, we maximize l(Θ∣X,y)l\left({\boldsymbol{\Theta }}| {\bf{X}},{\bf{y}})subject to the constraint (6.3)hj(Θ∣xi)=exp{ηj(xi)}−k1jpi1−…−kmjpim≥0∀i,j.{h}_{j}\left({\boldsymbol{\Theta }}| {{\bf{x}}}_{i})=\exp \left\{{\eta }_{j}\left({{\bf{x}}}_{i})\right\}-{k}_{1j}{p}_{i1}\hspace{0.33em}-\ldots -{k}_{mj}{p}_{im}\ge 0\hspace{1em}\forall i,j.This optimization problem can be solved using most nonlinear programming optimization algorithms. Here, we optimize ℬ{\mathcal{ {\mathcal B} }}and Γ\Gamma both subject to the constraint given in (6.3) and r1,…,rd{r}_{1},\ldots ,{r}_{d}subject to the constraint that r1,…,rd>0{r}_{1},\ldots ,{r}_{d}\gt 0using augmented Lagrangian methods and constrained optimization by linear approximation (COBYLA) methods. This type of optimization can be done in R programming using the auglag function in the nloptr package, (see [33]). For more information regarding these two methods, see [6,21]. Due to a large number of parameters, we use an iterative maximization procedure similar to [11] and [28] for obtaining estimates: (1)Obtain initial size estimates r10,…,rd0{r}_{1}^{0},\ldots ,{r}_{d}^{0}by maximizing the marginal log-likelihood functions with respect to r1,…,rd{r}_{1},\ldots ,{r}_{d}.(2)Obtain initial estimates ℬ0=[β10,…,βd0]T{{\mathcal{ {\mathcal B} }}}^{0}={\left[{{\boldsymbol{\beta }}}_{{\bf{1}}}^{0},\ldots ,{{\boldsymbol{\beta }}}_{d}^{0}]}^{T}by fitting a negative binomial regression model to the marginal distributions.(3)Obtain initial estimates Γ0=[γ10,…,γm0]T{\Gamma }^{0}={\left[{{\boldsymbol{\gamma }}}_{{\bf{1}}}^{0},\ldots ,{{\boldsymbol{\gamma }}}_{m}^{0}]}^{T}by fitting a multinomial logistic model to Z given in equation (4.3).(4)Obtain Ω0{{\boldsymbol{\Omega }}}^{0}by maximizing l(Ω∣X,y,r10,…,rd0,ℬ0,Γ0)l\left({\boldsymbol{\Omega }}| {\bf{X}},{\bf{y}},{r}_{1}^{0},\ldots ,{r}_{d}^{0},{{\mathcal{ {\mathcal B} }}}^{0},{\Gamma }^{0}).(5)Having obtained initial estimates, obtain r1i,…,rdi{r}_{1}^{i},\ldots ,{r}_{d}^{i}by maximizing l(r1,…,rd∣X,y,Ωi−1,ℬi−1,Γi−1)l\left({r}_{1},\ldots ,{r}_{d}| {\bf{X}},{\bf{y}},{{\boldsymbol{\Omega }}}^{i-1},{{\mathcal{ {\mathcal B} }}}^{i-1},{\Gamma }^{i-1})(6)Similarly, obtain ℬi{{\mathcal{ {\mathcal B} }}}^{i}, Γi{\Gamma }^{i}, Ωi{{\boldsymbol{\Omega }}}^{i}by maximizing l(ℬ∣X,y,r1i−1,…,rdi−1,Ωi−1,Γi−1)l\left({\mathcal{ {\mathcal B} }}| {\bf{X}},{\bf{y}},{r}_{1}^{i-1},\ldots ,{r}_{d}^{i-1},{{\boldsymbol{\Omega }}}^{i-1},{\Gamma }^{i-1}), l(Γ∣X,y,r1i−1,…,rdi−1,Ωi−1,ℬi−1)l\left(\Gamma | {\bf{X}},{\bf{y}},{r}_{1}^{i-1},\ldots ,{r}_{d}^{i-1},{{\boldsymbol{\Omega }}}^{i-1},{{\mathcal{ {\mathcal B} }}}^{i-1}), and l(Ω∣X,y,r1i−1,…,rdi−1,ℬi−1,Γi−1)l\left({\boldsymbol{\Omega }}| {\bf{X}},{\bf{y}},{r}_{1}^{i-1},\ldots ,{r}_{d}^{i-1},{{\mathcal{ {\mathcal B} }}}^{i-1},{\Gamma }^{i-1}), respectively, for i=1,2,3,…i=1,2,3,\ldots .(7)Repeat steps 5 and 6 until some convergence criteria are met.Here, convergence is obtained if the maximum of the differences between the current and previous parameter estimates is less than some threshold ε\varepsilon . However, we note that convergence in Akaike Information Criteria (AIC) or Bayesian Information Criteria may be used as well. We may also use existing packages such as “optim” and “maxLik” packages in the R environment to find MLEs by maximizing the aforementioned log-likelihood function.7Numerical exampleHere, we provide numerical examples to illustrate the proposed methodologies using simulated data and apply the method to real data.7.1Simulation studiesHere, we consider 500 equidistance points of independent variables xxfrom 0 to 1 and simulate responses yi=[yi1,yi2]′∼MMINB(r1,r2,μi1,μi2,pi1,pi2,ρ){{\bf{y}}}_{i}=[{y}_{i1},{y}_{i2}]^{\prime} \hspace{0.33em} \sim \hspace{0.33em}{\rm{MMINB}}\left({r}_{1},{r}_{2},{\mu }_{i1},{\mu }_{i2},{p}_{i1},{p}_{i2},\rho )considering two inflation points [0,0]′\left[0,0]^{\prime} and [1,1]′\left[1,1]^{\prime} (i.e., m=2m=2) with probabilities pi1{p}_{i1}and pi2{p}_{i2}, respectively, and we consider the Gaussian copula with parameter Ω={ρ}{\boldsymbol{\Omega }}=\left\{\rho \right\}, where ρ\rho is the correlation coefficient of the marginal distributions . The size parameters are considered as r1=5,r2=10{r}_{1}=5,{r}_{2}=10, and the correlation parameter for Gaussian copula is considered as ρ=0.7\rho =0.7. The linear predictors for λik{\lambda }_{ik}are assumed as follows: η1(x)=1.4−2.0xη2(x)=−0.5+1.5x,\begin{array}{rcl}{\eta }_{1}\left(x)& =& 1.4-2.0x\\ {\eta }_{2}\left(x)& =& -0.5+1.5x,\end{array}and the link function is assumed as the log link function.Again, the linear predictor, for odds of pil{p}_{il}are as follows: ξ1(x)=−1.7+3.0xξ2(x)=−1.5+2.5x,\begin{array}{rcl}{\xi }_{1}\left(x)& =& -1.7+3.0x\\ {\xi }_{2}\left(x)& =& -1.5+2.5x,\end{array}and the link function is assumed as a log link function for generating the data. We simulate n=500n=500responses y{\bf{y}}at design points xx.Researchers usually fit the multivariate count responses assuming multivariate negative binomial distribution without considering the inflation points in the data. We may obtain better performance from the fitted model if we consider the presence of inflation points in the data. So, here, we consider two models, M1{M}_{1}and M2{M}_{2}, to fit the data, where the model M1{M}_{1}is the usual multivariate negative binomial regression model and the model M2{M}_{2}is the proposed multivariate multiple inflated negative binomial regression model, where in both models, we use Gaussian copula for the dependence structure. Models M1{M}_{1}and M2{M}_{2}are described as follows:Model M1{M}_{1}: Here, we fit the simulated data with a multivariate negative binomial regression model without considering the presence of inflation points in the data. The form of the multivariate negative binomial distribution is found using the Gaussian copula. The linear predictors are assumed as follows: η1(x)=β10+β11xη2(x)=β20+β21x,\begin{array}{rcl}{\eta }_{1}\left(x)& =& {\beta }_{10}+{\beta }_{11}x\\ {\eta }_{2}\left(x)& =& {\beta }_{20}+{\beta }_{21}x,\end{array}and the link function is considered a log link function for fitting the data. The unknown parameters of the model are estimated using the MLE. The estimated values of the parameters, along with their standard errors and pp-values, are reported in Table 1. The AIC value of the fitted model M1{M}_{1}is 2795.562.Table 1Parameter estimates, standard errors, pp-values, and the AIC value of the fitted model M1{M}_{1}using 500 data pointsParametersEstimatesSEpp-valueβ10=1.4{\beta }_{10}=1.41.52950.1049<0.0001\lt 0.0001β11=−2{\beta }_{11}=-2−2.2583-2.25830.2000<0.0001\lt 0.0001β20=−0.5{\beta }_{20}=-0.5−0.1078-0.10780.14080.0443β21=1.5{\beta }_{21}=1.50.59770.22070.0067r1=5{r}_{1}=50.76870.0783<0.0001\lt 0.0001r2=10{r}_{2}=100.57910.0603<0.0001\lt 0.0001ρ=0.7\rho =0.70.83350.0150<0.0001\lt 0.0001Log-likelihood−1390.781-1390.781AIC2795.562Model M2{M}_{2}: Here, we fit the simulated data with the proposed model, assuming the presence of inflation points in the data. Since the data contain two inflation points (0,0) and (1,0), we fit a multivariate doubly inflated negative binomial regression model using the data. The forms of the linear predictors are assumed as follows: η1(x)=β10+β11x,η2(x)=β20+β21x,\begin{array}{rcl}{\eta }_{1}\left(x)& =& {\beta }_{10}+{\beta }_{11}x,\\ {\eta }_{2}\left(x)& =& {\beta }_{20}+{\beta }_{21}x,\end{array}and ξ1(x)=γ10+γ11x,ξ2(x)=γ20+γ21x,\begin{array}{rcl}{\xi }_{1}\left(x)& =& {\gamma }_{10}+{\gamma }_{11}x,\\ {\xi }_{2}\left(x)& =& {\gamma }_{20}+{\gamma }_{21}x,\end{array}and the corresponding link functions are assumed as log link function and multivariate logistic link function, respectively, as described in Section 5. The estimated values of the parameters using the MLE and the AIC value of the fitted model M2{M}_{2}are given in Table 2.Table 2Parameter estimates, standard errors, pp-values, and the AIC value of the fitted model M2{M}_{2}using 500 simulated data pointsParametersEstimatesSEpp-valueβ10=1.4{\beta }_{10}=1.41.48600.0902<0.0001\lt 0.0001β11=−2{\beta }_{11}=-2−2.2178-2.21780.1838<0.0001\lt 0.0001β20=−0.5{\beta }_{20}=-0.5−0.2385-0.23850.10380.0215β21=1.5{\beta }_{21}=1.50.69880.23310.0027γ10=−1.7{\gamma }_{10}=-1.7−1.9139-1.91390.2509<0.0001\lt 0.0001γ11=3.0{\gamma }_{11}=3.04.16350.4403<0.0001\lt 0.0001γ20=−1.5{\gamma }_{20}=-1.5−1.2426-1.24260.2113<0.0001\lt 0.0001γ21=2.5{\gamma }_{21}=2.53.00980.3594<0.0001\lt 0.0001r1=5{r}_{1}=55.79181.56300.0002r2=10{r}_{2}=108.88435.71330.0119ρ=0.7\rho =0.70.68970.0446<0.0001\lt 0.0001Log-likelihood−1059.785-1059.785AIC2141.569From the simulation study, we observe that the AIC value of model M2{M}_{2}is smaller than the AIC value of model M1{M}_{1}. The difference of the deviances of models M1{M}_{1}and M2{M}_{2}is 661.992 with 4 degrees of freedom (pp-value <0.0001\lt 0.0001). This shows that we obtain a significantly better fit using model M2{M}_{2}compared to model M1{M}_{1}by considering the existence of the multiple inflated points in the data.7.2Application to work week hours for couplesHere, we apply the proposed method of fitting a multivariate multiple inflated negative binomial regression model to real data. We model the work week hours for husbands and wives from the cps91 data set found in the Wooldridge package in R. The data set comes from the May 1991 Current Population Survey with a sample size of n=5,634n=\hspace{0.1em}\text{5,634}\hspace{0.1em}and contains factors such as education, age, ethnicity, and income for husbands and wives. We consider Yi1=“hushrs”{Y}_{i1}={\mbox{''}hushrs\mbox{''}}(husband’s weekly hours) and Yi2=“hours”{Y}_{i2}=\hspace{0.1em}\text{&#x201C;}\hspace{0.1em}{hours}\hspace{0.0em}\text{&#x201D;}\hspace{0.1em}(wife’s weekly hours) given in the data and assume Yi=[Yi1,Yi2]T∼MMINB(r1,r2,μi1,μi2,pi1,…,pim,ρ){{\bf{Y}}}_{i}={\left[{Y}_{i1},{Y}_{i2}]}^{T}\hspace{0.33em} \sim \hspace{0.33em}{\rm{MMINB}}\left({r}_{1},{r}_{2},{\mu }_{i1},{\mu }_{i2},{p}_{i1},\ldots ,{p}_{im},\rho ). We model λi1{\lambda }_{i1}, λi2{\lambda }_{i2}and pi1,pi2,…,pim{p}_{i1},{p}_{i2},\ldots ,{p}_{im}using the covariates, x1=“huseduc”{x}_{1}=\hspace{0.1em}\text{&#x201C;}\hspace{0.1em}{huseduc}\hspace{0.0em}\text{&#x201D;}\hspace{0.1em}, x2=“educ”{x}_{2}=\hspace{0.1em}\text{&#x201C;}\hspace{0.1em}{educ}\hspace{0.0em}\text{&#x201D;}\hspace{0.1em}, and x3=“expersq”{x}_{3}=\hspace{0.1em}\text{&#x201C;}\hspace{0.1em}{expersq}\hspace{0.0em}\text{&#x201D;}\hspace{0.1em}, which are corresponding to the husband’s years of education, wife’s years of education, and wife’s years of experience squared, respectively.We consider five different models to fit the data corresponding to five different number of inflation points m=0,1,…,4m=0,1,\ldots ,4. The inflation points are chosen according to their frequency of occurring given in Table 3. For all models, we use the Gaussian copula to find the form of the multivariate distribution function. The models are described as follows:Table 3Ten points with highest countsPointCount(40,0)\left(40,0)851(40,40)\left(40,40)781(0,0)\left(0,0)468(0,40)\left(0,40)214(50,0)\left(50,0)178(60,0)\left(60,0)156(50,40)\left(50,40)140(60,40)\left(60,40)93(45,0)\left(45,0)84(40,35)\left(40,35)79Model M(0){M}^{\left(0)}: Here, we consider the proposed model with m=0m=0, i.e., the usual multivariate negative binomial regression model without considering the presence of inflation points in the data. The linear predictors are assumed as follows: η1(x)=β10+β11x1+β12x2+β13x3η2(x)=β20+β21x1+β22x2+β23x3,\begin{array}{rcl}{\eta }_{1}\left({\bf{x}})& =& {\beta }_{10}+{\beta }_{11}{x}_{1}+{\beta }_{12}{x}_{2}+{\beta }_{13}{x}_{3}\\ {\eta }_{2}\left({\bf{x}})& =& {\beta }_{20}+{\beta }_{21}{x}_{1}+{\beta }_{22}{x}_{2}+{\beta }_{23}{x}_{3},\end{array}and the link function is considered a log link function for fitting the data. We estimate the parameters using MLE, and the estimated parameter values are presented in Table 4. The log-likelihood and AIC values are −46554.65-46554.65and 93131.3, respectively.Table 4Parameter estimates for m=0m=0ParameterEstimateStandard errorpp-valuer1{r}_{1}1.11340.0242<0.0001\lt 0.0001r2{r}_{2}0.26540.0058<0.0001\lt 0.0001ρ\rho 0.08410.0129<0.0001\lt 0.0001β10{\beta }_{10}3.61830.0128<0.0001\lt 0.0001β11{\beta }_{11}−0.0160-0.01600.01660.3380β12{\beta }_{12}0.10200.0162<0.0001\lt 0.0001β13{\beta }_{13}−0.1124-0.11240.0139<0.0001\lt 0.0001β20{\beta }_{20}3.00770.0260<0.0001\lt 0.0001β21{\beta }_{21}0.17860.0332<0.0001\lt 0.0001β22{\beta }_{22}−0.0432-0.04320.03350.1965β23{\beta }_{23}−0.0715-0.07150.02870.0127Log-likelihood−46,554.65AIC93,131.3Model M(1){M}^{\left(1)}: We fit the data with Model M1 assuming the presence of one inflation point, i.e., m=1m=1in the data and the inflation point is (40,0)\left(40,0)having the maximum number of occurrence among all observations. The linear predictors of the model are considered as follows: η1(x)=β10+β11x1+β12x2+β13x3,η2(x)=β20+β21x1+β22x2+β23x3,\begin{array}{rcl}{\eta }_{1}\left({\bf{x}})& =& {\beta }_{10}+{\beta }_{11}{x}_{1}+{\beta }_{12}{x}_{2}+{\beta }_{13}{x}_{3},\\ {\eta }_{2}\left({\bf{x}})& =& {\beta }_{20}+{\beta }_{21}{x}_{1}+{\beta }_{22}{x}_{2}+{\beta }_{23}{x}_{3},\end{array}and ξ1(x)=γ10+γ11x1+γ12x2+γ13x3,{\xi }_{1}\left({\bf{x}})={\gamma }_{10}+{\gamma }_{11}{x}_{1}+{\gamma }_{12}{x}_{2}+{\gamma }_{13}{x}_{3},and the corresponding link functions are assumed as log link function and multivariate logistic link function, respectively, as described in Section 5. The estimated parameter values using the MLE are given in Table 5. We observe that the log-likelihood value is −43644.01-43644.01, and the AIC value is 87318.02 for the fitted model.Table 5Parameter estimates for m=1m=1ParameterEstimateStandard errorpp-valuer1{r}_{1}0.87830.0204<0.0001\lt 0.0001r2{r}_{2}0.41260.0096<0.0001\lt 0.0001ρ\rho 0.18830.0135<0.0001\lt 0.0001β10{\beta }_{10}3.61560.0130<0.0001\lt 0.0001β11{\beta }_{11}−0.0482-0.04820.01530.0016β12{\beta }_{12}0.11620.0156<0.0001\lt 0.0001β13{\beta }_{13}−0.1026-0.10260.0137<0.0001\lt 0.0001β20{\beta }_{20}3.00790.0234<0.0001\lt 0.0001β21{\beta }_{21}0.18150.0305<0.0001\lt 0.0001β22{\beta }_{22}−0.0439-0.04390.03020.1464β23{\beta }_{23}0.00600.02620.8191γ10{\gamma }_{10}−1.7672-1.76720.0384<0.0001\lt 0.0001γ11{\gamma }_{11}−0.3118-0.31180.0433<0.0001\lt 0.0001γ12{\gamma }_{12}−0.0491-0.04910.04430.2681γ13{\gamma }_{13}−0.00003-0.000030.03050.9991Log-likelihood−43,644.01AIC87,318.02Model M(2){M}^{\left(2)}: This model fits the data considering presence of two inflation points (40,0)\left(40,0)and (40,40)\left(40,40)(i.e., m=2m=2) with linear predictors η1(x)=β10+β11x1+β12x2+β13x3,η2(x)=β20+β21x1+β22x2+β23x3,\begin{array}{rcl}{\eta }_{1}\left({\bf{x}})& =& {\beta }_{10}+{\beta }_{11}{x}_{1}+{\beta }_{12}{x}_{2}+{\beta }_{13}{x}_{3},\\ {\eta }_{2}\left({\bf{x}})& =& {\beta }_{20}+{\beta }_{21}{x}_{1}+{\beta }_{22}{x}_{2}+{\beta }_{23}{x}_{3},\end{array}and ξ1(x)=γ10+γ11x1+γ12x2+γ13x3,ξ2(x)=γ20+γ21x1+γ22x2+γ23x3,\begin{array}{rcl}{\xi }_{1}\left({\bf{x}})& =& {\gamma }_{10}+{\gamma }_{11}{x}_{1}+{\gamma }_{12}{x}_{2}+{\gamma }_{13}{x}_{3},\\ {\xi }_{2}\left({\bf{x}})& =& {\gamma }_{20}+{\gamma }_{21}{x}_{1}+{\gamma }_{22}{x}_{2}+{\gamma }_{23}{x}_{3},\end{array}and the corresponding link functions as log link function and multivariate logistic link function, respectively, as given in Section 5. The MLEs of the model parameters are given in Table 6. The log-likelihood and AIC values for the fitted model are reported as −37803.86-37803.86and 75645.72, respectively.Table 6Parameter estimates for m=2m=2ParameterEstimateStandard errorpp-valuer1{r}_{1}0.68230.0171<0.0001\lt 0.0001r2{r}_{2}0.31370.0081<0.0001\lt 0.0001ρ\rho 0.14660.0153<0.0001\lt 0.0001β10{\beta }_{10}3.61850.0134<0.0001\lt 0.0001β11{\beta }_{11}−0.0446-0.04460.01510.0033β12{\beta }_{12}0.10860.0151<0.0001\lt 0.0001β13{\beta }_{13}−0.1057-0.10570.0139<0.0001\lt 0.0001β20{\beta }_{20}3.01300.0221<0.0001\lt 0.0001β21{\beta }_{21}0.18570.0288<0.0001\lt 0.0001β22{\beta }_{22}−0.0436-0.04360.02880.1298β23{\beta }_{23}0.00980.02540.6989γ10{\gamma }_{10}−1.5908-1.59080.0390<0.0001\lt 0.0001γ11{\gamma }_{11}−0.2921-0.29210.0447<0.0001\lt 0.0001γ12{\gamma }_{12}0.02880.04460.5182γ13{\gamma }_{13}−0.0362-0.03620.03830.3454γ20{\gamma }_{20}−1.6585-1.65850.0399<0.0001\lt 0.0001γ21{\gamma }_{21}0.15580.05130.0024γ22{\gamma }_{22}−0.1196-0.11960.04870.0140γ23{\gamma }_{23}−0.13794-0.137940.04480.0021Log-likelihood−37803.86-37803.86AIC75645.72Model M(3){M}^{\left(3)}: Here, we consider the presence of three inflation points (40,0)\left(40,0), (40,40)\left(40,40), and (0,0)\left(0,0)(i.e., m=3m=3) in the data. The link functions are the same as model M(2){M}^{\left(2)}, and the linear predictors are as follows: η1(x)=β10+β11x1+β12x2+β13x3,η2(x)=β20+β21x1+β22x2+β23x3,\begin{array}{rcl}{\eta }_{1}\left({\bf{x}})& =& {\beta }_{10}+{\beta }_{11}{x}_{1}+{\beta }_{12}{x}_{2}+{\beta }_{13}{x}_{3},\\ {\eta }_{2}\left({\bf{x}})& =& {\beta }_{20}+{\beta }_{21}{x}_{1}+{\beta }_{22}{x}_{2}+{\beta }_{23}{x}_{3},\end{array}and ξ1(x)=γ10+γ11x1+γ12x2+γ13x3,ξ2(x)=γ20+γ21x1+γ22x2+γ23x3,ξ3(x)=γ30+γ31x1+γ32x2+γ33x3.\begin{array}{rcl}{\xi }_{1}\left({\bf{x}})& =& {\gamma }_{10}+{\gamma }_{11}{x}_{1}+{\gamma }_{12}{x}_{2}+{\gamma }_{13}{x}_{3},\\ {\xi }_{2}\left({\bf{x}})& =& {\gamma }_{20}+{\gamma }_{21}{x}_{1}+{\gamma }_{22}{x}_{2}+{\gamma }_{23}{x}_{3},\\ {\xi }_{3}\left({\bf{x}})& =& {\gamma }_{30}+{\gamma }_{31}{x}_{1}+{\gamma }_{32}{x}_{2}+{\gamma }_{33}{x}_{3}.\end{array}The estimated model parameters, AIC, and log-likelihood values are reported in Table 7. We see that the log-likelihood is −36672.42-36672.42, and the AIC value is 73390.84 for the fitted model.Table 7Parameter estimates for m=3m=3ParameterEstimateStandard errorpp-valuer1{r}_{1}1.57660.0436<0.0001\lt 0.0001r2{r}_{2}0.46010.0123<0.0001\lt 0.0001ρ\rho −0.2087-0.20870.0146<0.0001\lt 0.0001β10{\beta }_{10}3.63040.0103<0.0001\lt 0.0001β11{\beta }_{11}0.00010.01320.9936β12{\beta }_{12}0.08550.0131<0.0001\lt 0.0001β13{\beta }_{13}−0.1128-0.11280.0114<0.0001\lt 0.0001β20{\beta }_{20}3.02410.0205<0.0001\lt 0.0001β21{\beta }_{21}0.20460.0276<0.0001\lt 0.0001β22{\beta }_{22}−0.0180-0.01800.02790.5174β23{\beta }_{23}0.00410.02410.8658γ10{\gamma }_{10}−2.2440-2.24400.0571<0.0001\lt 0.0001γ11{\gamma }_{11}−0.1402-0.14020.05500.0108γ12{\gamma }_{12}−0.3246-0.32460.0519<0.0001\lt 0.0001γ13{\gamma }_{13}0.44120.0437<0.0001\lt 0.0001γ20{\gamma }_{20}−1.4689-1.46890.0396<0.0001\lt 0.0001γ21{\gamma }_{21}−0.3672-0.36720.0521<0.0001\lt 0.0001γ22{\gamma }_{22}−0.0001-0.00010.06350.9985γ23{\gamma }_{23}0.00090.03990.9817γ30{\gamma }_{30}−1.5335-1.53350.0404<0.0001\lt 0.0001γ31{\gamma }_{31}0.14850.05370.0057γ32{\gamma }_{32}−0.1768-0.17680.05220.0007γ33{\gamma }_{33}−0.0596-0.05960.04580.1939Log-likelihood−36672.42-36672.42AIC73,390.84Model M(4){M}^{\left(4)}: For this model, we consider the presence of four inflation points (40,0)\left(40,0), (40,40)\left(40,40), (0,0)\left(0,0), and (0,40)\left(0,40)(i.e., m=4m=4) in the data. The linear predictors are assumed as follows: η1(x)=β10+β11x1+β12x2+β13x3,η2(x)=β20+β21x1+β22x2+β23x3,\begin{array}{rcl}{\eta }_{1}\left({\bf{x}})& =& {\beta }_{10}+{\beta }_{11}{x}_{1}+{\beta }_{12}{x}_{2}+{\beta }_{13}{x}_{3},\\ {\eta }_{2}\left({\bf{x}})& =& {\beta }_{20}+{\beta }_{21}{x}_{1}+{\beta }_{22}{x}_{2}+{\beta }_{23}{x}_{3},\end{array}and ξ1(x)=γ10+γ11x1+γ12x2+γ13x3,ξ2(x)=γ20+γ21x1+γ22x2+γ23x3,ξ3(x)=γ30+γ31x1+γ32x2+γ33x3,ξ4(x)=γ40+γ41x1+γ42x2+γ43x3,\begin{array}{rcl}{\xi }_{1}\left({\bf{x}})& =& {\gamma }_{10}+{\gamma }_{11}{x}_{1}+{\gamma }_{12}{x}_{2}+{\gamma }_{13}{x}_{3},\\ {\xi }_{2}\left({\bf{x}})& =& {\gamma }_{20}+{\gamma }_{21}{x}_{1}+{\gamma }_{22}{x}_{2}+{\gamma }_{23}{x}_{3},\\ {\xi }_{3}\left({\bf{x}})& =& {\gamma }_{30}+{\gamma }_{31}{x}_{1}+{\gamma }_{32}{x}_{2}+{\gamma }_{33}{x}_{3},\\ {\xi }_{4}\left({\bf{x}})& =& {\gamma }_{40}+{\gamma }_{41}{x}_{1}+{\gamma }_{42}{x}_{2}+{\gamma }_{43}{x}_{3},\end{array}and the corresponding link functions are assumed as log link function and multivariate logistic link function, respectively, as described in Section 5. The values of the log-likelihood and AIC are −34957.68-34957.68and 69969.36, respectively, and the estimated model parameters are reported in Table 8.Table 8Parameter estimates for m=4m=4ParameterEstimateStandard errorpp-valuer1{r}_{1}3.48840.1045<0.0001\lt 0.0001r2{r}_{2}0.41510.0116<0.0001\lt 0.0001ρ\rho −0.1276-0.12760.0153<0.0001\lt 0.0001β10{\beta }_{10}3.63410.0083<0.0001\lt 0.0001β11{\beta }_{11}0.00070.01030.9434β12{\beta }_{12}0.07790.0103<0.0001\lt 0.0001β13{\beta }_{13}−0.0922-0.09220.0093<0.0001\lt 0.0001β20{\beta }_{20}3.02560.0204<0.0001\lt 0.0001β21{\beta }_{21}0.17800.0267<0.0001\lt 0.0001β22{\beta }_{22}−0.0082-0.00820.02650.7567β23{\beta }_{23}0.02210.02430.3618γ10{\gamma }_{10}−2.1693-2.16930.0572<0.0001\lt 0.0001γ11{\gamma }_{11}−0.0052-0.00520.05620.9267γ12{\gamma }_{12}−0.4043-0.40430.0508<0.0001\lt 0.0001γ13{\gamma }_{13}0.46680.0434<0.0001\lt 0.0001γ20{\gamma }_{20}−1.4054-1.40540.0400<0.0001\lt 0.0001γ21{\gamma }_{21}−0.2712-0.27120.0445<0.0001\lt 0.0001γ22{\gamma }_{22}−0.1092-0.10920.04680.0197γ23{\gamma }_{23}0.07330.03990.0659γ30{\gamma }_{30}−1.4692-1.46920.0408<0.0001\lt 0.0001γ31{\gamma }_{31}0.19170.05340.0003γ32{\gamma }_{32}−0.2320-0.23200.0525<0.0001\lt 0.0001γ33{\gamma }_{33}−0.0072-0.00720.04660.8766γ40{\gamma }_{40}−2.8408-2.84080.0764<0.0001\lt 0.0001γ41{\gamma }_{41}0.14460.09150.1138γ42{\gamma }_{42}−0.3998-0.39980.0792<0.0001\lt 0.0001γ43{\gamma }_{43}0.39150.0656<0.0001\lt 0.0001Log-likelihood−34957.68-34957.68AIC69,969.36Note that the same set of covariates is used for the marginal means and inflation probabilities since we are primarily interested in the change in the log-likelihood value and AIC for different values of mm.Histograms of the bivariate distribution and marginal distributions are given in Figures 4 and 5, respectively. We observe that the model with m=4m=4provides lowest AIC value among all the models with m∈{0,1,2,3,4}m\in \left\{0,1,2,3,4\right\}. Also, according to the chi-square test, the model considering four inflation points in this data provides a significantly best fit than the other models with m∈{0,1,2,3}m\in \left\{0,1,2,3\right\}.Figure 4Histogram of hushrs and hours variables from the cps91 data set.Figure 5Histogram of marginal distributions with estimated density. (a) Hushrs variable. (b) Hours variable.8ConclusionWe have presented a regression model in which the response vector y=[y1,y2,…]{\bf{y}}=[{y}_{1},{y}_{2},\ldots ]is assumed to follow a multiple-inflated negative binomial distribution. Copula methods were introduced to model the dependence structure of y{\bf{y}}through a correlation matrix RR. Due to the complexity of the likelihood function, an iterative procedure has been proposed for estimating model parameters. First, we consider simulated data and fit the proposed multivariate multiple inflated negative binomial model and the usual multivariate negative binomial regression model using Gaussian copula without considering inflation points in the data. We observe that the proposed model considering the presence of inflation points in the data, provides a significantly better fit than the usual model without considering the same. Next, five models accounting for mminflation points were fit to a current population survey data set, with m=0,…4m=0,\ldots 4. Accounting for the multiple-inflated structure generally reduced the AIC and log-likelihood value.One problem of the model is the difficulty of estimating parameters, despite the proposed iterative procedure. In particular, derivative-based optimization methods often fail at some point during the stepwise estimation process, and hence, the use of (COBYLA) methods mentioned in Section 6. With this in mind, as well as the constraints imposed in Section 6, there is no reason to expect parameter estimates to be global maxima.By using the proposed methodologies, one can choose any copula to find the multivariate distribution function from marginals. We used two-dimensional copulas for our numerical examples and found expected results by estimating parameters using the proposed algorithms. However, the estimation of model parameters becomes cumbersome as the response dimension increases. To address this problem, we may use vine copulas and other existing concepts in the literature to have computationally stable estimates from the model. One of the future directions may be fitting the model using different data sets and improving the estimation process.

Journal

Dependence Modelingde Gruyter

Published: Jan 1, 2022

Keywords: Copulas; inflated distribution; multivariate count data; negative binomial distributions; regression models; 62H10

There are no references for this article.