Skip to content
BY 4.0 license Open Access Published by De Gruyter Open Access November 2, 2022

Multiple inflated negative binomial regression for correlated multivariate count data

  • Joseph Mathews , Sumangal Bhattacharya , Sumen Sen and Ishapathik Das EMAIL logo
From the journal Dependence Modeling

Abstract

This article aims to provide a method of regression for multivariate multiple inflated count responses assuming the responses follow a negative binomial distribution. Negative binomial regression models are common in the literature for modeling univariate and multivariate count data. However, two problems commonly arise in modeling such data: choice of the multivariate form of the underlying distribution and modeling the zero-inflated structure of the data. Copula functions have become a popular solution to the former problem because they can model the response variables’ dependence structure. The latter problem is often solved by modeling an assumed latent variable Z generating excess zero-valued counts in addition to the underlying distribution. However, despite their flexibility, zero-inflation models do not account for the case of m additional inflated values at k 1 , k 2 , , k m . We propose a multivariate multiple-inflated negative binomial regression model for modeling such cases. Furthermore, we present an iterative procedure for estimating model parameters using maximum likelihood estimation. The multivariate distribution functions considering the dependence structure of the response vectors are found using copula methods. The proposed method is illustrated using simulated data and applied to real data.

MSC 2010: 62H10

1 Introduction

Generalized 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 ) and mean μ ( > 0 ) given by

(1.1) f NB ( y r , μ ) = Γ ( y + r ) Γ ( r ) y ! r r + μ r μ r + μ y , y = 0 , 1 ,

is often preferred. From [10], the variance of Y is given by

(1.2) Var ( Y ) = μ + μ 2 r .

See [10] for a full derivation of (1.1) using a Poisson-gamma mixture. Notice that as r , the distribution reduces to the Poisson in which E ( Y ) = Var ( Y ) = μ . Thus, the dispersion parameter α = 1 / r allows 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 k -valued counts, where k is 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 Y 1 Pois ( α θ ) , Y 2 Pois ( β θ ) , where α > 0 , β > 0 , θ > 0 , and a gamma mixing distribution with shape parameter r and scale parameter λ = 1 :

(1.3) 0 ( α θ ) y 1 e α θ y 1 ! ( β θ ) y 2 e β θ y 2 ! θ r 1 e θ Γ ( r ) d θ .

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 Ω . 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 Ω 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.

2 Copulas and Sklar’s theorem

This 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 ] .

Definition 1

A function C : [ 0 , 1 ] d [ 0 , 1 ] is a d -dimensional copula if the following are true:

  1. C ( 1 , , y j , , 1 ) = y j , j = 1 , 2 , , d with ( 1 , , y j , , 1 ) [ 0 , 1 ] d ,

  2. C ( y 1 , y 2 , , y d ) = 0 for ( y 1 , y 2 , , y d ) [ 0 , 1 ] d with at least one of y j = 0 for j = 1 , 2 , , d ,

  3. C is d -nondecreasing, i.e., for any u 1 = ( u 1 ( 1 ) , u 2 ( 1 ) , , u d ( 1 ) ) , u 2 = ( u 1 ( 2 ) , u 2 ( 2 ) , , u d ( 2 ) ) belong to [ 0 , 1 ] d with u j ( 1 ) u j ( 2 ) , for all j = 1 , 2 , , d :

    (2.1) a 1 = 1 2 a 2 = 1 2 a d = 1 2 ( 1 ) a 1 + a 2 + + a d C ( u 1 ( a 1 ) , u 2 ( a 2 ) , , u d ( a d ) ) 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 1

Let Y 1 , Y 2 , , Y d be random variables with marginal distribution functions F 1 , F 2 , , 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 y 1 , y 2 , , y d R

    (2.2) F ( y 1 , y 2 , , y d ) = C ( F 1 ( y 1 ) , F 2 ( y 2 ) , , F d ( y d ) ) .

  2. If Y 1 , Y 2 , , Y d are continuous then the copula C is unique. Otherwise, C can be uniquely determined on a d-dimensional rectangle Range ( F 1 ) × Range ( F 2 ) × × Range ( F d ) .

Copulas are popular because of their ability to model dependence among Y 1 , Y 2 , , Y d . One popular copula is the Gaussian copula since this dependence is modeled through a correlation matrix R .

Definition 2

The Gaussian copula is given by the function:

(2.3) C ( u 1 , u 2 , , u d R ) = Φ R ( Φ 1 ( u 1 ) , Φ 1 ( u 2 ) , , Φ 1 ( u d ) ) ,

where Φ 1 is the inverse cumulative distribution function of the standard Gaussian distribution and Φ R is the joint cumulative distribution function of a standard multivariate Gaussian distribution with correlation matrix R .

Definition 3

The Gaussian copula density [2] is given by

(2.4) c ( u 1 , u 2 , , u d R ) = 1 R exp 1 2 U T × ( R 1 I d ) × U ,

where U = [ Φ 1 ( u 1 ) , Φ 1 ( u 2 ) , , Φ 1 ( u d ) ] T , Φ is the univariate standard Gaussian distribution, R is the correlation matrix of the standard multivariate Gaussian distribution, and I d be the d × d identity 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 R as a function of ρ given by

    (2.5) R ( ρ ) = ρ 1 1 T ( 1 ρ ) I d ,

    where I d is a d -dimensional identity matrix, ρ 1 d 1 , 1 , and 1 is a d -dimensional column vector of ones. From [20], we have

    (2.6) R 1 ( ρ ) = 1 1 ρ I d ρ ( 1 ρ ) { 1 + ( d 1 ) ρ } 11 t ,

    where M 2 = diag ( 0 , 1 , , 1 , 0 ) and 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 , j )th element of R is given by ρ i j , with ρ ( 1 , 1 ) . By [4], the inverse of R as a function of ρ is given by

    (2.7) R 1 ( ρ ) = 1 1 ρ 2 ( I d ρ 2 M 2 ρ M 1 ) ,

    where M 2 = diag ( 0 , 1 , , 1 , 0 ) and 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 ρ are given in Figure 1.

Figure 1 
               Gaussian copula density functions with different values of 
                     
                        
                        
                           ρ
                        
                        \rho 
                     
                  : (a) 
                     
                        
                        
                           ρ
                           =
                           0.9
                        
                        \rho =0.9
                     
                   and (b) 
                     
                        
                        
                           ρ
                           =
                           −
                           0.25
                        
                        \rho =-0.25
                     
                  .
Figure 1

Gaussian copula density functions with different values of ρ : (a) ρ = 0.9 and (b) ρ = 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.

3 Multivariate discrete probability mass function using a copula

Here, we discuss using a copula to construct a form of multivariate discrete probability mass function of a correlated d -dimensional random vector from its marginals. Let Y = [ Y 1 , Y 2 , , Y d ] T be discrete random vectors with respective marginal cumulative distribution functions F 1 ( y 1 θ 1 ) , F 2 ( y 2 θ 2 ) , , F d ( y d θ d ) . Then their joint probability mass function P ( Y = y ) formed using a copula C is given by

(3.1) f Y C ( y Ψ ) = P ( Y 1 = y 1 , Y 2 = y 2 , , Y d = y d Ψ ) = a 1 = 1 2 a 2 = 1 2 a d = 1 2 ( 1 ) a 1 + a 2 + + a d C ( u 1 ( a 1 ) , u 2 ( a 2 ) , , u d ( a d ) Ψ ) ,

where y D = { 0 , 1 , 2 , } d , Ψ = [ θ 1 T , θ 2 T , , θ d T , Ω T ] T , Ω is the copula parameters, u j ( 1 ) = F j ( y j θ j ) , and u j ( 2 ) = F j ( y j 1 θ j ) , where j = 1 , , d . Here, we restrict ourselves to the case where F 1 ( y 1 θ 1 ) , F 2 ( y 2 θ 2 ) , , F d ( y d θ d ) are negative binomial distribution functions with size parameter r j , mean parameter μ j , and θ j = [ r j , μ j ] T . Recall that the probability mass function in Section 1 is given by

(3.2) f NB ( y j r j , μ j ) = Γ ( y j + r j ) Γ ( r j ) y j ! r j r j + μ j r j μ j r j + μ j y j , y j = 0 , 1 , .

The corresponding distribution function is given by

(3.3) F NB ( y ) = P ( Y j y ) = { y j : y j y } f NB ( y j r j , μ j ) , y R ,

where r j , μ j > 0 and y j { 0 , 1 , } . Bivariate negative binomial distributions with different correlation values and parameters are given in Figure 2.

Figure 2 
               Bivariate negative binomial distribution using a Gaussian copula with (a) 
                     
                        
                        
                           
                              
                                 r
                              
                              
                                 1
                              
                           
                           =
                           1.5
                           ,
                           
                              
                                 μ
                              
                              
                                 1
                              
                           
                           =
                           20
                           ,
                           ρ
                           =
                           −
                           0.30
                        
                        {r}_{1}=1.5,{\mu }_{1}=20,\rho =-0.30
                     
                   and (b) 
                     
                        
                        
                           
                              
                                 r
                              
                              
                                 2
                              
                           
                           =
                           1.5
                           ,
                           
                              
                                 μ
                              
                              
                                 2
                              
                           
                           =
                           20
                           ,
                           ρ
                           =
                           0.80
                        
                        {r}_{2}=1.5,{\mu }_{2}=20,\rho =0.80
                     
                  .
Figure 2

Bivariate negative binomial distribution using a Gaussian copula with (a) r 1 = 1.5 , μ 1 = 20 , ρ = 0.30 and (b) r 2 = 1.5 , μ 2 = 20 , ρ = 0.80 .

4 Multivariate multiple-inflated negative binomial distribution function

Here, we describe the method of constructing a multivariate multiple-inflated negative binomial distribution function with the help of a multinomial distributed latent variable Z .

Let Y be a d -dimensional random vector that follows a multivariate multiple-inflated negative binomial distribution with m (a nonnegative integer) inflation points. We find a form of the distribution function F Y ( y ) by considering a latent variable Z with the probability mass function given by

(4.1) f Z ( z ) = p m if z = m p m 1 if z = m 1 p 1 if z = 1 1 q if z = 0 ,

where q = p m + p m 1 + + p 1 < 1 with p m , p m 1 , , p 1 ( 0 , 1 ) , and m is a nonnegative integer. Then the conditional probability mass function is given by

(4.2) f Y Z , Ψ ( y z , Ψ ) = 1 if y = k m and z = m 1 if y = k m 1 and z = m 1 1 if y = k 1 and z = 1 f Y C ( y Ψ ) if y D and z = 0 ,

where D = { 0 , 1 , 2 , } d , d is the dimension of Y , f Y C ( Ψ ) is a multivariate negative binomial probability mass function formed using a copula, Ψ = [ r 1 , , r d , μ 1 , , μ d , Ω T ] T is a parameter list, Ω is the copula parameters, and k l = [ k l 1 , k l 2 , , k l d ] T is a d -dimensional inflation point for l = 1 , 2 , , m . By (4.1) and (4.2), we have the joint probability mass function:

(4.3) f Y , Z Ψ ( y , z Ψ ) = p m if z = m and y = k m p m 1 if z = m 1 and y = k m 1 p 1 if z = 1 and y = k 1 ( 1 q ) f Y C ( y Ψ ) if z = 0 and y D .

Finally, by (4.3), we have the multivariate multiple-inflated negative binomial probability mass function given by

(4.4) f Y Ψ ( y Ψ ) = p m + ( 1 q ) f Y C ( y Ψ ) if y = k m p m 1 + ( 1 q ) f Y C ( y Ψ ) if y = k m 1 p 1 + ( 1 q ) f Y C ( y Ψ ) if y = k 1 ( 1 q ) f Y C ( y Ψ ) otherwise ,

where y 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 d -dimensional random vector Y can be found using (4.4):

(4.5) f Y j ( y j ) = y 1 y j 1 y j + 1 y d f Y ( y Ψ ) ,

where y j is the j th component of y = [ y 1 , y 2 , , y d ] T D and Y j is the j th components of the random vector Y for j = 1 , 2 , , d . Therefore, the model simplifies to

(4.6) f Y j μ j , r j ( y j μ j , r j ) = p m + ( 1 q ) f NB ( y j μ j , r j ) if y j = k m j p m 1 + ( 1 q ) f NB ( y j μ j , r j ) if y j = k ( m 1 ) j p 1 + ( 1 q ) f NB ( y j μ j , r j ) if y i = k 1 j ( 1 q ) f NB ( y j μ j , r j ) otherwise ,

where y j { 0 , 1 , 2 , } , k l = [ k l 1 , k l 2 , , k l d ] T is a d -dimensional inflation point for l = 1 , 2 , , m , and f NB ( y j μ j , r j ) is a negative binomial probability mass function with size r j and mean μ j . From (4.6), the marginal expected value of the component Y j is given by

(4.7) E ( Y j ) = y j y j f Y j μ j , r j ( y j μ j , r j ) = 0 f Y j μ j , r j ( 0 μ j , r j ) + k 1 j f Y j μ j , r j ( k 1 j μ j , r j ) + + k m j f Y j μ j , r j ( k m j μ j , r j ) + y j y j 0 , k 1 j , , k m j y j f Y j ( y j μ j , r j ) = k 1 j p 1 + + k m j p j + ( 1 q ) y j y j f NB ( y j μ j , r j ) = k 1 j p 1 + + k m j p j + ( 1 q ) μ j = λ j , [ say ] .

As mentioned earlier, we have

μ j = λ j k 1 j p 1 k 2 j p 2 k m j 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 3 
               Bivariate multiple-inflated negative binomial distribution using a Gaussian copula with parameters (a) 
                     
                        
                        
                           
                              
                                 r
                              
                              
                                 1
                              
                           
                           =
                           1.5
                           ,
                           
                              
                                 μ
                              
                              
                                 1
                              
                           
                           =
                           20
                        
                        {r}_{1}=1.5,{\mu }_{1}=20
                     
                   and (b) 
                     
                        
                        
                           
                              
                                 r
                              
                              
                                 2
                              
                           
                           =
                           1.5
                           ,
                           
                              
                                 μ
                              
                              
                                 2
                              
                           
                           =
                           20
                        
                        {r}_{2}=1.5,{\mu }_{2}=20
                     
                  . (a) MMINB distribution with 
                     
                        
                        
                           
                              
                                 k
                              
                              
                                 1
                              
                           
                           =
                           
                              (
                              
                                 0
                                 ,
                                 0
                              
                              )
                           
                        
                        {{\bf{k}}}_{1}=\left(0,0)
                     
                   and 
                     
                        
                        
                           
                              
                                 k
                              
                              
                                 2
                              
                           
                           =
                           
                              (
                              
                                 25
                                 ,
                                 20
                              
                              )
                           
                        
                        {{\bf{k}}}_{2}=\left(25,20)
                     
                  . (b) MMINB distribution with 
                     
                        
                        
                           
                              
                                 k
                              
                              
                                 1
                              
                           
                           =
                           
                              (
                              
                                 0
                                 ,
                                 0
                              
                              )
                           
                        
                        {{\bf{k}}}_{1}=\left(0,0)
                     
                  , 
                     
                        
                        
                           
                              
                                 k
                              
                              
                                 2
                              
                           
                           =
                           
                              (
                              
                                 0
                                 ,
                                 20
                              
                              )
                           
                        
                        {{\bf{k}}}_{2}=\left(0,20)
                     
                  , 
                     
                        
                        
                           
                              
                                 k
                              
                              
                                 3
                              
                           
                           =
                           
                              (
                              
                                 20
                                 ,
                                 0
                              
                              )
                           
                        
                        {{\bf{k}}}_{3}=\left(20,0)
                     
                  , and 
                     
                        
                        
                           
                              
                                 k
                              
                              
                                 4
                              
                           
                           =
                           
                              (
                              
                                 20
                                 ,
                                 20
                              
                              )
                           
                        
                        {{\bf{k}}}_{4}=\left(20,20)
                     
                  .
Figure 3

Bivariate multiple-inflated negative binomial distribution using a Gaussian copula with parameters (a) r 1 = 1.5 , μ 1 = 20 and (b) r 2 = 1.5 , μ 2 = 20 . (a) MMINB distribution with k 1 = ( 0 , 0 ) and k 2 = ( 25 , 20 ) . (b) MMINB distribution with k 1 = ( 0 , 0 ) , k 2 = ( 0 , 20 ) , k 3 = ( 20 , 0 ) , and k 4 = ( 20 , 20 ) .

5 Multivariate multiple-inflated negative binomial regression model

In this section, we propose a regression method for multivariate multiple-inflated negative binomial d -dimensional responses y R d with ν -dimensional covariates vector x R ν . Let y i = [ y i 1 , y i 2 , , y i d ] T be a vector of response variables conditional on the i th observation x i = [ x i 1 , x i 2 , , x i ν ] T , i = 1 , 2 , , n . We assume that y i MMINB ( r 1 , , r d , μ i 1 , μ i d , p i 1 , , p i m ) , where μ i j and r j are the marginal mean and size, respectively, of f NB ( y i j ) and p i l is the probability of observing an inflated value at k l , l = 1 , 2 , , m .

From equation (4.7), we obtain,

E ( Y i j x i ) = k 1 j p i 1 + k 2 j p i 2 + + k m j p i m + ( 1 q i ) μ i j = λ i j .

We note that only μ i j and p i l are conditional on x i with i = 1 , , n , j = 1 , , d , and l = 1 , , m . We use the log-link function for modeling (4.7) :

(5.1) log ( λ i 1 ) = η 1 ( x i ) = f 1 ( x i ) β 1 log ( λ i 2 ) = η 2 ( x i ) = f 2 ( x i ) β 2 log ( λ i d ) = η d ( x i ) = f d ( x i ) β d ,

where η j ( x i ) is the linear predictor, i.e., a polynomial function of x i with coefficient β j , j = 1 , 2 , , d , f j ( x i ) is a vector-valued function and β j is an τ j × 1 dimensional vector. Note that τ j is a nonnegative integer representing the dimension of the vector β j . It follows that the marginal mean of f NB ( y i j ) is modeled as follows:

(5.2) μ i j = exp { η j ( x i ) } k 1 j p i 1 k m j p i m 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 p i l using a set of linear predictors:

(5.3) ln p i m p i 0 = ln P ( Z = m ) P ( Z = 0 ) = ξ m ( x i ) = g m ( x i ) γ m ln p i ( m 1 ) p i 0 = ln P ( Z = m 1 ) P ( Z = 0 ) = ξ m 1 ( x i ) = g m 1 ( x i ) γ m 1 ln p i 1 p i 0 = ln P ( Z = 1 ) P ( Z = 0 ) = ξ 1 ( x i ) = g 1 ( x i ) γ 1 ,

where ξ l ( x i ) is the linear predictor, i.e., a polynomial function of x i with coefficient γ l , l = 1 , 2 , , m , g l ( x i ) is a vector-valued function and γ l is an ψ l × 1 dimensional vector of covariates l = 1 , 2 , , m and p i 0 = p i 1 + p i 2 + + p i m . From (5.3), we obtain

(5.4) p i m = exp { g m ( x i ) γ m } 1 + l = 1 m exp { g l ( x i ) γ l } p i ( m 1 ) = exp { g m 1 ( x i ) γ m 1 } 1 + l = 1 m exp { g l ( x i ) γ l } p i 1 = exp { g 1 ( x i ) γ 1 } 1 + l = 1 m exp { g l ( x i ) γ l } .

Next, we provide a method of estimating unknown parameters of the proposed regression model.

6 Estimation of parameters

We estimate model parameters Θ = [ r 1 , , r d , T , Γ T , Ω T ] T using the method of maximum likelihood estimation (MLE), where = [ β 1 , β 2 , , β d ] T , Γ = [ γ 1 , γ 2 , , γ m ] T and Ω is the copula parameters. The likelihood function is given by

(6.1) l ( Θ X , y ) = { y i y i = k m } log { p i m + ( 1 q i ) f C Φ R ( y i Θ ) } × { y i y i = k m 1 } log { p i ( m 1 ) + ( 1 q i ) f C Φ R ( y i Θ ) } × × { y i y i { k m , k m 1 , , k 1 } } log { ( 1 q i ) f C Φ R ( y i Θ ) } ,

where X is a design matrix. It follows that the log-likelihood function is given by

(6.2) l ( Θ X , y ) = { y i y i = k m } log { p i m + ( 1 q i ) f C Φ R ( y i Θ ) } + { y i y i = k m 1 } log { p i ( m 1 ) + ( 1 q i ) f C Φ R ( y i Θ ) } + + { y i y i { k m , k m 1 , , k 1 } } log { ( 1 q i ) f C Φ R ( y i Θ ) }

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 μ i j given in (5.2) remains positive. Therefore, we maximize l ( Θ X , y ) subject to the constraint

(6.3) h j ( Θ x i ) = exp { η j ( x i ) } k 1 j p i 1 k m j p i m 0 i , j .

This optimization problem can be solved using most nonlinear programming optimization algorithms. Here, we optimize and Γ both subject to the constraint given in (6.3) and r 1 , , r d subject to the constraint that r 1 , , r d > 0 using 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 r 1 0 , , r d 0 by maximizing the marginal log-likelihood functions with respect to r 1 , , r d .

  2. Obtain initial estimates 0 = [ β 1 0 , , β d 0 ] T by fitting a negative binomial regression model to the marginal distributions.

  3. Obtain initial estimates Γ 0 = [ γ 1 0 , , γ m 0 ] T by fitting a multinomial logistic model to Z given in equation (4.3).

  4. Obtain Ω 0 by maximizing l ( Ω X , y , r 1 0 , , r d 0 , 0 , Γ 0 ) .

  5. Having obtained initial estimates, obtain r 1 i , , r d i by maximizing l ( r 1 , , r d X , y , Ω i 1 , i 1 , Γ i 1 )

  6. Similarly, obtain i , Γ i , Ω i by maximizing l ( X , y , r 1 i 1 , , r d i 1 , Ω i 1 , Γ i 1 ) , l ( Γ X , y , r 1 i 1 , , r d i 1 , Ω i 1 , i 1 ) , and l ( Ω X , y , r 1 i 1 , , r d i 1 , i 1 , Γ i 1 ) , respectively, for i = 1 , 2 , 3 , .

  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 ε . 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.

7 Numerical example

Here, we provide numerical examples to illustrate the proposed methodologies using simulated data and apply the method to real data.

7.1 Simulation studies

Here, we consider 500 equidistance points of independent variables x from 0 to 1 and simulate responses y i = [ y i 1 , y i 2 ] MMINB ( r 1 , r 2 , μ i 1 , μ i 2 , p i 1 , p i 2 , ρ ) considering two inflation points [ 0 , 0 ] and [ 1 , 1 ] (i.e., m = 2 ) with probabilities p i 1 and p i 2 , respectively, and we consider the Gaussian copula with parameter Ω = { ρ } , where ρ is the correlation coefficient of the marginal distributions . The size parameters are considered as r 1 = 5 , r 2 = 10 , and the correlation parameter for Gaussian copula is considered as ρ = 0.7 . The linear predictors for λ i k are assumed as follows:

η 1 ( x ) = 1.4 2.0 x η 2 ( x ) = 0.5 + 1.5 x ,

and the link function is assumed as the log link function.

Again, the linear predictor, for odds of p i l are as follows:

ξ 1 ( x ) = 1.7 + 3.0 x ξ 2 ( x ) = 1.5 + 2.5 x ,

and the link function is assumed as a log link function for generating the data. We simulate n = 500 responses y at design points x .

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, M 1 and M 2 , to fit the data, where the model M 1 is the usual multivariate negative binomial regression model and the model 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 M 1 and M 2 are described as follows:

Model 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 + β 11 x η 2 ( x ) = β 20 + β 21 x ,

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 p -values, are reported in Table 1. The AIC value of the fitted model M 1 is 2795.562.

Table 1

Parameter estimates, standard errors, p -values, and the AIC value of the fitted model M 1 using 500 data points

Parameters Estimates SE p -value
β 10 = 1.4 1.5295 0.1049 < 0.0001
β 11 = 2 2.2583 0.2000 < 0.0001
β 20 = 0.5 0.1078 0.1408 0.0443
β 21 = 1.5 0.5977 0.2207 0.0067
r 1 = 5 0.7687 0.0783 < 0.0001
r 2 = 10 0.5791 0.0603 < 0.0001
ρ = 0.7 0.8335 0.0150 < 0.0001
Log-likelihood 1390.781
AIC 2795.562

Model 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 + β 11 x , η 2 ( x ) = β 20 + β 21 x ,

and

ξ 1 ( x ) = γ 10 + γ 11 x , ξ 2 ( x ) = γ 20 + γ 21 x ,

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 M 2 are given in Table 2.

Table 2

Parameter estimates, standard errors, p -values, and the AIC value of the fitted model M 2 using 500 simulated data points

Parameters Estimates SE p -value
β 10 = 1.4 1.4860 0.0902 < 0.0001
β 11 = 2 2.2178 0.1838 < 0.0001
β 20 = 0.5 0.2385 0.1038 0.0215
β 21 = 1.5 0.6988 0.2331 0.0027
γ 10 = 1.7 1.9139 0.2509 < 0.0001
γ 11 = 3.0 4.1635 0.4403 < 0.0001
γ 20 = 1.5 1.2426 0.2113 < 0.0001
γ 21 = 2.5 3.0098 0.3594 < 0.0001
r 1 = 5 5.7918 1.5630 0.0002
r 2 = 10 8.8843 5.7133 0.0119
ρ = 0.7 0.6897 0.0446 < 0.0001
Log-likelihood 1059.785
AIC 2141.569

From the simulation study, we observe that the AIC value of model M 2 is smaller than the AIC value of model M 1 . The difference of the deviances of models M 1 and M 2 is 661.992 with 4 degrees of freedom ( p -value < 0.0001 ). This shows that we obtain a significantly better fit using model M 2 compared to model M 1 by considering the existence of the multiple inflated points in the data.

7.2 Application to work week hours for couples

Here, 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,634 and contains factors such as education, age, ethnicity, and income for husbands and wives. We consider Y i 1 = “hushrs” (husband’s weekly hours) and Y i 2 = hours (wife’s weekly hours) given in the data and assume Y i = [ Y i 1 , Y i 2 ] T MMINB ( r 1 , r 2 , μ i 1 , μ i 2 , p i 1 , , p i m , ρ ) . We model λ i 1 , λ i 2 and p i 1 , p i 2 , , p i m using the covariates, x 1 = huseduc , x 2 = educ , and x 3 = expersq , 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 , , 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 3

Ten points with highest counts

Point Count
( 40 , 0 ) 851
( 40 , 40 ) 781
( 0 , 0 ) 468
( 0 , 40 ) 214
( 50 , 0 ) 178
( 60 , 0 ) 156
( 50 , 40 ) 140
( 60 , 40 ) 93
( 45 , 0 ) 84
( 40 , 35 ) 79

Model M ( 0 ) : Here, we consider the proposed model with m = 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 + β 11 x 1 + β 12 x 2 + β 13 x 3 η 2 ( x ) = β 20 + β 21 x 1 + β 22 x 2 + β 23 x 3 ,

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 and 93131.3, respectively.

Table 4

Parameter estimates for m = 0

Parameter Estimate Standard error p -value
r 1 1.1134 0.0242 < 0.0001
r 2 0.2654 0.0058 < 0.0001
ρ 0.0841 0.0129 < 0.0001
β 10 3.6183 0.0128 < 0.0001
β 11 0.0160 0.0166 0.3380
β 12 0.1020 0.0162 < 0.0001
β 13 0.1124 0.0139 < 0.0001
β 20 3.0077 0.0260 < 0.0001
β 21 0.1786 0.0332 < 0.0001
β 22 0.0432 0.0335 0.1965
β 23 0.0715 0.0287 0.0127
Log-likelihood −46,554.65
AIC 93,131.3

Model M ( 1 ) : We fit the data with Model M1 assuming the presence of one inflation point, i.e., m = 1 in the data and the inflation point is ( 40 , 0 ) having the maximum number of occurrence among all observations. The linear predictors of the model are considered as follows:

η 1 ( x ) = β 10 + β 11 x 1 + β 12 x 2 + β 13 x 3 , η 2 ( x ) = β 20 + β 21 x 1 + β 22 x 2 + β 23 x 3 ,

and

ξ 1 ( x ) = γ 10 + γ 11 x 1 + γ 12 x 2 + γ 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 , and the AIC value is 87318.02 for the fitted model.

Table 5

Parameter estimates for m = 1

Parameter Estimate Standard error p -value
r 1 0.8783 0.0204 < 0.0001
r 2 0.4126 0.0096 < 0.0001
ρ 0.1883 0.0135 < 0.0001
β 10 3.6156 0.0130 < 0.0001
β 11 0.0482 0.0153 0.0016
β 12 0.1162 0.0156 < 0.0001
β 13 0.1026 0.0137 < 0.0001
β 20 3.0079 0.0234 < 0.0001
β 21 0.1815 0.0305 < 0.0001
β 22 0.0439 0.0302 0.1464
β 23 0.0060 0.0262 0.8191
γ 10 1.7672 0.0384 < 0.0001
γ 11 0.3118 0.0433 < 0.0001
γ 12 0.0491 0.0443 0.2681
γ 13 0.00003 0.0305 0.9991
Log-likelihood −43,644.01
AIC 87,318.02

Model M ( 2 ) : This model fits the data considering presence of two inflation points ( 40 , 0 ) and ( 40 , 40 ) (i.e., m = 2 ) with linear predictors

η 1 ( x ) = β 10 + β 11 x 1 + β 12 x 2 + β 13 x 3 , η 2 ( x ) = β 20 + β 21 x 1 + β 22 x 2 + β 23 x 3 ,

and

ξ 1 ( x ) = γ 10 + γ 11 x 1 + γ 12 x 2 + γ 13 x 3 , ξ 2 ( x ) = γ 20 + γ 21 x 1 + γ 22 x 2 + γ 23 x 3 ,

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 and 75645.72, respectively.

Table 6

Parameter estimates for m = 2

Parameter Estimate Standard error p -value
r 1 0.6823 0.0171 < 0.0001
r 2 0.3137 0.0081 < 0.0001
ρ 0.1466 0.0153 < 0.0001
β 10 3.6185 0.0134 < 0.0001
β 11 0.0446 0.0151 0.0033
β 12 0.1086 0.0151 < 0.0001
β 13 0.1057 0.0139 < 0.0001
β 20 3.0130 0.0221 < 0.0001
β 21 0.1857 0.0288 < 0.0001
β 22 0.0436 0.0288 0.1298
β 23 0.0098 0.0254 0.6989
γ 10 1.5908 0.0390 < 0.0001
γ 11 0.2921 0.0447 < 0.0001
γ 12 0.0288 0.0446 0.5182
γ 13 0.0362 0.0383 0.3454
γ 20 1.6585 0.0399 < 0.0001
γ 21 0.1558 0.0513 0.0024
γ 22 0.1196 0.0487 0.0140
γ 23 0.13794 0.0448 0.0021
Log-likelihood 37803.86
AIC 75645.72

Model M ( 3 ) : Here, we consider the presence of three inflation points ( 40 , 0 ) , ( 40 , 40 ) , and ( 0 , 0 ) (i.e., m = 3 ) in the data. The link functions are the same as model M ( 2 ) , and the linear predictors are as follows:

η 1 ( x ) = β 10 + β 11 x 1 + β 12 x 2 + β 13 x 3 , η 2 ( x ) = β 20 + β 21 x 1 + β 22 x 2 + β 23 x 3 ,

and

ξ 1 ( x ) = γ 10 + γ 11 x 1 + γ 12 x 2 + γ 13 x 3 , ξ 2 ( x ) = γ 20 + γ 21 x 1 + γ 22 x 2 + γ 23 x 3 , ξ 3 ( x ) = γ 30 + γ 31 x 1 + γ 32 x 2 + γ 33 x 3 .

The estimated model parameters, AIC, and log-likelihood values are reported in Table 7. We see that the log-likelihood is 36672.42 , and the AIC value is 73390.84 for the fitted model.

Table 7

Parameter estimates for m = 3

Parameter Estimate Standard error p -value
r 1 1.5766 0.0436 < 0.0001
r 2 0.4601 0.0123 < 0.0001
ρ 0.2087 0.0146 < 0.0001
β 10 3.6304 0.0103 < 0.0001
β 11 0.0001 0.0132 0.9936
β 12 0.0855 0.0131 < 0.0001
β 13 0.1128 0.0114 < 0.0001
β 20 3.0241 0.0205 < 0.0001
β 21 0.2046 0.0276 < 0.0001
β 22 0.0180 0.0279 0.5174
β 23 0.0041 0.0241 0.8658
γ 10 2.2440 0.0571 < 0.0001
γ 11 0.1402 0.0550 0.0108
γ 12 0.3246 0.0519 < 0.0001
γ 13 0.4412 0.0437 < 0.0001
γ 20 1.4689 0.0396 < 0.0001
γ 21 0.3672 0.0521 < 0.0001
γ 22 0.0001 0.0635 0.9985
γ 23 0.0009 0.0399 0.9817
γ 30 1.5335 0.0404 < 0.0001
γ 31 0.1485 0.0537 0.0057
γ 32 0.1768 0.0522 0.0007
γ 33 0.0596 0.0458 0.1939
Log-likelihood 36672.42
AIC 73,390.84

Model M ( 4 ) : For this model, we consider the presence of four inflation points ( 40 , 0 ) , ( 40 , 40 ) , ( 0 , 0 ) , and ( 0 , 40 ) (i.e., m = 4 ) in the data. The linear predictors are assumed as follows:

η 1 ( x ) = β 10 + β 11 x 1 + β 12 x 2 + β 13 x 3 , η 2 ( x ) = β 20 + β 21 x 1 + β 22 x 2 + β 23 x 3 ,

and

ξ 1 ( x ) = γ 10 + γ 11 x 1 + γ 12 x 2 + γ 13 x 3 , ξ 2 ( x ) = γ 20 + γ 21 x 1 + γ 22 x 2 + γ 23 x 3 , ξ 3 ( x ) = γ 30 + γ 31 x 1 + γ 32 x 2 + γ 33 x 3 , ξ 4 ( x ) = γ 40 + γ 41 x 1 + γ 42 x 2 + γ 43 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 values of the log-likelihood and AIC are 34957.68 and 69969.36, respectively, and the estimated model parameters are reported in Table 8.

Table 8

Parameter estimates for m = 4

Parameter Estimate Standard error p -value
r 1 3.4884 0.1045 < 0.0001
r 2 0.4151 0.0116 < 0.0001
ρ 0.1276 0.0153 < 0.0001
β