Multiple inflated negative binomial regression for correlated multivariate count data

Joseph Mathews; Sumangal Bhattacharya; Sumen Sen; Ishapathik Das

doi:10.1515/demo-2022-0149

Article Open Access

Multiple inflated negative binomial regression for correlated multivariate count data

Joseph Mathews , Sumangal Bhattacharya , Sumen Sen and Ishapathik Das

Published/Copyright: November 2, 2022

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information

From the journal Dependence Modeling Volume 10 Issue 1

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.

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

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:

C ( 1 , … , y j , … , 1 ) = y j , ∀ j = 1 , 2 , … , d with ( 1 , … , y j , … , 1 ) ∈ [ 0 , 1 ] d ,
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 ,
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:

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 ) ) .
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:

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

Obtain initial size estimates r 1 0 , … , r d 0 by maximizing the marginal log-likelihood functions with respect to r 1 , … , r d .
Obtain initial estimates ℬ 0 = [ β 1 0 , … , β d 0 ] T by fitting a negative binomial regression model to the marginal distributions.
Obtain initial estimates Γ 0 = [ γ 1 0 , … , γ m 0 ] T by fitting a multinomial logistic model to Z given in equation (4.3).
Obtain Ω 0 by maximizing l ( Ω ∣ X , y , r 1 0 , … , r d 0 , ℬ 0 , Γ 0 ) .
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 )
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 , … .
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
β 10	3.6341	0.0083	< 0.0001
β 11	0.0007	0.0103	0.9434
β 12	0.0779	0.0103	< 0.0001
β 13	− 0.0922	0.0093	< 0.0001
β 20	3.0256	0.0204	< 0.0001
β 21	0.1780	0.0267	< 0.0001
β 22	− 0.0082	0.0265	0.7567
β 23	0.0221	0.0243	0.3618
γ 10	− 2.1693	0.0572	< 0.0001
γ 11	− 0.0052	0.0562	0.9267
γ 12	− 0.4043	0.0508	< 0.0001
γ 13	0.4668	0.0434	< 0.0001
γ 20	− 1.4054	0.0400	< 0.0001
γ 21	− 0.2712	0.0445	< 0.0001
γ 22	− 0.1092	0.0468	0.0197
γ 23	0.0733	0.0399	0.0659
γ 30	− 1.4692	0.0408	< 0.0001
γ 31	0.1917	0.0534	0.0003
γ 32	− 0.2320	0.0525	< 0.0001
γ 33	− 0.0072	0.0466	0.8766
γ 40	− 2.8408	0.0764	< 0.0001
γ 41	0.1446	0.0915	0.1138
γ 42	− 0.3998	0.0792	< 0.0001
γ 43	0.3915	0.0656	< 0.0001
Log-likelihood	− 34957.68
AIC	69,969.36

Note 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 m .

Histograms of the bivariate distribution and marginal distributions are given in Figures 4 and 5, respectively. We observe that the model with m = 4 provides lowest AIC value among all the models with m ∈ { 0 , 1 , 2 , 3 , 4 } . 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 } .

Figure 4

Histogram of hushrs and hours variables from the cps91 data set.

Figure 5

Histogram of marginal distributions with estimated density. (a) Hushrs variable. (b) Hours variable.

8 Conclusion

We have presented a regression model in which the response vector y = [ y 1 , y 2 , … ] is assumed to follow a multiple-inflated negative binomial distribution. Copula methods were introduced to model the dependence structure of y through a correlation matrix R . 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 m inflation points were fit to a current population survey data set, with m = 0 , … 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.

Funding information: This work was supported by the Science and Engineering Research Board, Government of India (SRG/2019/001226).
Computer programs used for computations: All the computations in the numerical examples to illustrate the proposed methodologies are performed using a language and environment for statistical computing called R [22].
Conflict of interest: On behalf of all authors, the corresponding author states that there is no conflict of interest.

References

[1] Aitchison, J., & Ho, C. (1989). The multivariate poisson-log normal distribution. Biometrika, 76(4), 643–653. 10.1093/biomet/76.4.643Search in Google Scholar

[2] Arbenz, P. (2013). Bayesian copulae distributions, with application to operational risk management–some comments. Methodology and Computing in Applied Probability, 15(1), 105–108. 10.1007/s11009-011-9224-0Search in Google Scholar

[3] Cameron, A. C., & Trivedi, P. K. (1986). Econometric models based on count data. comparisons and applications of some estimators and tests. Journal of Applied Econometrics, 1(1), 29–53. 10.1002/jae.3950010104Search in Google Scholar

[4] Chaganty, N. R. (1997). An alternative approach to the analysis of longitudinal data via generalized estimating equations. Journal of Statistical Planning and Inference, 63(1), 39–54. 10.1016/S0378-3758(96)00203-0Search in Google Scholar

[5] Famoye, F. (2010). On the bivariate negative binomial regression model. Journal of Applied Statistics, 37(6), 969–981. 10.1080/02664760902984618Search in Google Scholar

[6] Fortin, M., & Glowinski, R. (2000). Augmented Lagrangian methods: Applications to the numerical solution of boundary-value problems. Netherlands: Elsevier. Search in Google Scholar

[7] Genest, C., & Nešlehová, J. (2007). A primer on copulas for count data. ASTIN Bulletin: The Journal of the IAA, 37(2), 475–515. 10.2143/AST.37.2.2024077Search in Google Scholar

[8] Greene, W. (2008). Functional forms for the negative binomial model for count data. Economics Letters, 99(3), 585–590. 10.1016/j.econlet.2007.10.015Search in Google Scholar

[9] Greene, W. H. (1994). Accounting for excess zeros and sample selection in poisson and negative binomial regression models, 94–10.Search in Google Scholar

[10] Hilbe, J. M. (2011). Negative binomial regression. United Kingdom: Cambridge University Press. 10.1017/CBO9780511973420Search in Google Scholar

[11] Joe, H. (2005). Asymptotic efficiency of the two-stage estimation method for copula-based models. Journal of Multivariate Analysis, 94(2), 401–419. 10.1016/j.jmva.2004.06.003Search in Google Scholar

[12] Joe, H. (2014). Dependence modeling with copulas. Florida: CRC Press. 10.1201/b17116Search in Google Scholar

[13] Karlis, D., & Xekalaki, E. (2005). Mixed poisson distributions. International Statistical Review, 73, 35–38. 10.1111/j.1751-5823.2005.tb00250.xSearch in Google Scholar

[14] Lambert, D. (1992). Zero-inflated poisson regression, with an application to defects in manufacturing. Technometrics, 34(2), 1–14. 10.2307/1269547Search in Google Scholar

[15] Marshall, A. W., & Olkin, I. (1990). Multivariate distributions generated from mixtures of convolution and product families. Lecture Notes-Monograph Series, 16, 371–393. 10.1214/lnms/1215457574Search in Google Scholar

[16] Mathews, J., Sen, S., & Das, I. (2019). Multivariate doubly-inflated negative binomial distribution using Gaussian copula. In: Modern statistical methods for spatial and multivariate data, (pp. 147–161). Cham: Springer. 10.1007/978-3-030-11431-2_8Search in Google Scholar

[17] Mitchell, C., & Paulson, A. (1981). A new bivariate negative binomial distribution. Naval Research Logistics Quarterly, 28(3), 359–374. 10.1002/nav.3800280302Search in Google Scholar

[18] Mullahy, J. (1986). Specification and testing of some modified count data models. Journal of Econometrics, 33(3), 341–365. 10.1016/0304-4076(86)90002-3Search in Google Scholar

[19] Nikoloulopoulos, A. K., & Karlis, D. (2009). Modeling multivariate count data using copulas. Communications in Statistics-Simulation and Computation, 39(1), 172–187. 10.1080/03610910903391262Search in Google Scholar

[20] Olkin, I., & Pratt, J. W. (1958). Unbiased estimation of certain correlation coefficients. The Annals of Mathematical Statistics, 29(1), 201–211. 10.1214/aoms/1177706717Search in Google Scholar

[21] Powell, M. J. (1994). A direct search optimization method that models the objective and constraint functions by linear interpolation. In: Advances in optimization and numerical analysis, (pp. 51–67). New York: Springer. 10.1007/978-94-015-8330-5_4Search in Google Scholar

[22] R Core Team (2020). R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. Search in Google Scholar

[23] Sen, S., Sengupta, P., & Diawara, N. (2018). Doubly inflated poisson model using gaussian copula. Communications in Statistics-Theory and Methods, 47(12), 2848–2858. 10.1080/03610926.2017.1342831Search in Google Scholar

[24] Sengupta, P., Chaganty, N. R., & Sabo, R. T. (2016). Bivariate doubly inflated poisson models with applications. Journal of Statistical Theory and Practice, 10(1), 202–215. 10.1080/15598608.2015.1105166Search in Google Scholar

[25] Shea, J. M., & Brown, K. (2018). Wooldridge: 111 data sets from “introductory econometrics: A modern approach, 6e” by jeffrey m. Wooldridge. R package version, 1(1), 5. https://CRAN.R-project.org/package=wooldridge. Search in Google Scholar

[26] Sklar, M. (1959). Fonctions de repartition an dimensions et Leurs marges. Publications of the Institute of Statistics of the University of Paris, 8, 229–231. Search in Google Scholar

[27] So, S., Lee, D.-H., & Jung, B. C. (2011). An alternative bivariate zero-inflated negative binomial regression model using a copula. Economics Letters, 113(2), 183–185. 10.1016/j.econlet.2011.07.017Search in Google Scholar

[28] Song, P. X.-K., Fan, Y., & Kalbfleisch, J. D. (2005). Maximization by parts in likelihood inference. Journal of the American Statistical Association, 100(472), 1145–1158. 10.1198/016214505000000204Search in Google Scholar

[29] Song, X.-K., & Song, P. X. -K. (2007). Correlated data analysis: Modeling, analytics, and applications. New York: Springer Science & Business Media. Search in Google Scholar

[30] Trivedi, P., & Zimmer, D. (2017). A note on identification of bivariate copulas for discrete count data. Econometrics, 5(1), 10. 10.3390/econometrics5010010Search in Google Scholar

[31] Wooldridge, J. (2015). Introductory econometrics: A modern approach. Cengage Learning. Search in Google Scholar

[32] Yang, L., Frees, E. W., & Zhang, Z. (2020). Nonparametric estimation of copula regression models with discrete outcomes. Journal of the American Statistical Association, 115(530), 707–720. 10.1080/01621459.2018.1546586Search in Google Scholar

[33] Ypma, J. (2014). nloptr: R interface to nlopt. R package version 1.04. Search in Google Scholar

Received: 2022-01-09

Revised: 2022-09-07

Accepted: 2022-09-09

Published Online: 2022-11-02

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/demo-2022-0149

Keywords for this article

Copulas; inflated distribution; multivariate count data; negative binomial distributions; regression models

Creative Commons

BY 4.0