A Comparison of Bivariate Zero-Inflated Poisson Inverse Gaussian Regression Models with and without Exposure Variables

: In this paper, we focus on the comparison of the bivariate zero-inﬂated Poisson inverse Gaussian regression (BZIPIGR) type II model in two cases: with and without exposure variables. The BZIPIGR type II model is applied to analyze the occurrence of maternal and early neonatal mortality in South Sulawesi Province, Indonesia using 2019 data, which contain many zero values and have the issue of overdispersion in the response variable. Furthermore, to analyze the number of deaths in various areas, the exposure variable is considered. The maximum likelihood estimation (MLE) is used in parameter estimation, which involves numerical iteration and application of the Berndt–Hall–Hall–Hausman (BHHH) algorithm. Sum square error (SSE) serves as the criterion of model selection when exposure variables are included. The existence of exposure variables strongly affects the model’s accuracy, especially using the BZIPIGR type II model. According to the SSE and RMSE values, the BZIPIGR type II model with exposure variables performs better than the model without exposure variables in estimating parameter values. All predictors with exposure variables in this study had a signiﬁcant inﬂuence on the number of maternal and early neonatal mortalities.


Introduction
Count data that follow a Poisson distribution can be modeled by Poisson regression. Poisson regression requires the assumption that the mean and variance of the response variables must be equal (equidispersion), i.e., one parameter gives both the mean and variance of the distribution. In reality, this assumption may not hold, either with a variance larger than the mean (overdispersion) or otherwise (underdispersion) [1][2][3][4][5]. Such a violation can result in errors in decision making in hypothesis testing due to the occurrence of underestimates [6][7][8]. A model was developed based on mixed Poisson distribution to overcome the overdispersion problem, constituting a combination of Poisson distribution and other distributions, both discrete and continuous, including Poisson Inverse Gaussian (PIG) distribution [9].
Poisson regression can be used to model count data, but between events within a certain period, the data must be independent. For example, maternal mortality at a certain time does not affect maternal mortality at that time or another time. However, the number of mothers who die within a certain period in a given area may be influenced by other factors; for example, more births results in a higher occurrence of maternal mortality in an area. Thus, when comparing the number of mortalities between one area and another, it is necessary to consider the existence of other variables that can affect the number of occurrences, referred to as the exposure variable [10]. Frome and Morris (1989) [11] conducted research involving exposure variables with Poisson regression modeling. Violation of the equidispersion assumption, especially overdispersion within the Poisson distribution, can be caused by the number of zero values (excess zero) in the response variable [11][12][13].
Phang and Loh recommend a zero-inflated Poisson inverse Gaussian regression (ZIP-IGR) model for data with excess zeros and overdispersion [14]. Zero-inflated count data contains many zeros in the response variable. Hilbe conducted simulations with data containing excess zeros and overdispersion. Several models were included in this simulation, such as zero-inflated negative binomial regression (ZINBR), zero-inflated hurdle negative binomial regression (ZIHNBR) and zero-inflated generalized regression (ZIGPR). The simulation obtained the lowest AIC and BIC values with the ZIPIGR model, so the conclusion is that the ZIPIGR model is preferable to the other models [5]. Several other researchers have also used ZIPIGR to model with count data containing excess zeros and overdispersion [15,16].
The ZIPIGR model is limited to univariate cases; therefore, if the relevant application has more than one response variable, excess zeros or overdispersion, this model cannot be used. In this study, we used two correlated response variables: the number of maternal and early neonatal mortalities. Maternal mortality markedly affects the survival of the child, especially in cases of babies younger than one month old (early neonatal), where the basic needs of the baby are dependent on the mother, such as breast feeding. More than 50% of response variables in this study contained zero values and exhibited overdispersion. Therefore, a bivariate zero-inflated Poisson inverse Gaussian regression (BZIPIGR) model was developed.
The BZIPIGR type I model consists of two probability functions that accommodate the state of the response-variable pairs: all response variables are zero, and all response variables are not zero. However, not all cases consist of two response conditions. In the case of the number of maternal and early neonatal mortalities in South Sulawesi Province in 2018, there were four combinations of zero values in the response variable. Therefore, to identify the factors that significantly affected the number of maternal and early neonatal mortalities, a BZIPIGR type II model was developed to accommodate all possible combinations of zero values in the response variable. This is similar to the formation of the probability function in the BZIGPR model [17].
Several researchers have used Newton-Raphson iterations to maximize the likelihood function using the MLE method when the derivative of the likelihood function for the model parameters was in non-closed form. This algorithm requires a second derivative to form a Hessian matrix. Therefore, for a more complex probability distribution function, it is difficult to obtain a second derivative. According to Nugraha, the first-derivative approach can be used to inform the Hessian matrix using a Berndt-Hall-Hall-Hausman (BHHH) iteration [18]. The complexity of the probability function in the BZIPIGR type II model makes it difficult to obtain the second derivative. Therefore, we used the BHHH iteration method as an alternative to overcome this problem.

Bivariate Zero-Inflated Poisson Inverse Gaussian Regression Type II (BZIPIGR)
The BZIPIGR type II model is derived from the ZIPIGR model. The ZIPIGR model is a mixed of zero-inflated (ZI) and Poisson inverse Gaussian (PIG) distributions [16]. The zero value in the ZIPIG distribution can be divided into two sources, namely structural zero (zero state) with probability p 0 and sampling zero (mixed Poisson distribution state or PIG state) with probability (1 − p 0 ). If there are two random variables, Y 1 and Y 2 , with BZIPIG distribution with parameters (η h , γ h , τ) where h = 1, 2 or (Y 1 , Y 2 ) ∼ BZIPIG(η 1 , η 2 , τ, γ 1 , γ 2 ), and where η 1 and η 2 are the parameters of the Poisson distribution; τ is the dispersion parameter, and γ 1 , γ 2 is the parameter of the zero-inflated distribution, the probability function of Y 1 and Y 2 can be obtained with: Hence, the joint probability function of Y 1 and Y 2 can be written as where is the modified Bessel function of the third kind [19]. Therefore, Equation (1) becomes where K y 2 − 1 2 1 τ 1 + 2τη 2 is the modified Bessel function of the third kind [20]; where K y 1 − 1 2 1 τ 1 + 2τ(η 1 + η 2 ) is the modified Bessel function of the third kind [20] and where K y 1 +y 2 − 1 2 1 τ 1 + 2τ(η 1 + η 2 ) is the modified Bessel function of the third kind [20].

BZIPIGR Type II Model with Exposure Variables
The regression model for the BZIPIG distribution is divided into two equations: the regression model with exposure variable, as in Equation (7) [10]: where q hi is the exposure variable for the h-th response variable and the i-th observation; x is the vector of the predictor variable, and λ is the model parameters. The model for the ZI distribution can be written as Equation (9): where α is the vector of the zero-inflated parameters.
Based on the model in Equations (7) and (8), the probability function of the BZIPIGR type II model with exposures variable can be written as follows:

BZIPIGR Type II Model without Exposure Variable
The regression model for the BZIPIGR type II model without exposure variable, as in Equation (10) is [10] where x is the vector of the predictor variable for the i-th observation, and λ is the model parameters for the h-th response variable. Based on the model in Equations (9) and (10), the probability function of the BZIPIGR type II model without exposure variables can be written as follows:

Parameter Estimation of BZIPIGR Type II Model
Parameter estimation using the MLE method begins with taking n random samples, namely (Y 1i , Y 2i , X 1i , X 2i , . . . , X ki ), as there are four independent models in the BZIPIGR type II model, namely (y 1 = 0, y 2 = 0), (y 1 = 0, y 2 > 0), (y 1 > 0, y 2 = 0) and (y 1 > 0, y 2 > 0), where the number of observations for each model is n 1 , n 2 , n 3 and n 4 , respectively. The likelihood equation can be written as follows: Then, the log likelihood function is determined from Equation (11), and the following equation is obtained: T is a parameter vector of the BZIPIGR type II model, then, to maximize the log likelihood function in Equation (12), these equations are derived for each parameter and equalized to zero. The first derivative obtained for each parameter is in non-closed form (Appendix A), so the maximization process continues with numerical iteration, which uses the Berndt-Hall-Hall-Hausman (BHHH) algorithm according to the following steps [21]: Step 1. Determine the initial value of each parameter of the BZIPIGR type II model The initial values of model parameters are obtained from the regression model parameter values using OLS, whereas the dispersion parameter is obtained from the variance of the PIG distribution.
Step 2. Define gradient vector Step 3. Determine the Hessian matrix (H*) for the BHHH algorithm where g i θ (m) is the gradient vector at the i-th observation.
Step 4. Start the BHHH iteration using the following formula: The iteration of the BHHH algorithm starts at m = 0 and stops if θ (m+1) −θ (m) ≤ ε where ε is a small positive close to 0.

Hypothesis Testing of the BZIPIGR Type II Model
Simultaneous hypothesis testing for the BZIPIGR model type II parameters can be performed using the maximum likelihood ratio test (MLRT) method. The hypothesis to be tested is as follows: Before determining the test statistic to be used, the parameter set under the population (Ω) and H 0 (ω) are defined as follows: This simultaneous test uses test statistics where δ is the level of significance; a is the number of parameters under the population, and b is the number of parameters under . Partial hypothesis testing to determine the significance of the effect of each predictor in the model is divided into two tests, namely a partial test of parameters λ hj and α hj . The hypothesis statistical test for each parameter is as follows: a.
Partial test of parameter λ hj Partial test of parameter α hj Partial test of parameter τ H 0 : τ = 0;

Application
The BZIPIGR Type II model was applied to analyze the number of maternal and early neonatal mortalities in South Sulawesi Province, Indonesia, in 2019. Modeling in this study was carried out for two cases: with and without exposure variables. The response variables used in this study were two correlated variables: the number of maternal mortalities (Y 1 ) and the number of early neonatal mortalities (Y 2 ). There were six predictor variables: the percentage of visits by pregnant women (X 1 ), the percentage of pregnant women who received Fe3 tablets (X 2 ), percentage of births assisted by health workers (X 3 ), the percentage of mother who attended at least three postpartum maternal visits (X 4 ), the percentage of active integrated service posts (X 5 ) and the percentage of obstetric complications (X 6 ). The exposure variable used in this study was the number of live births in each sub-district (q). Table 1 shows the descriptive statistics of the response variables and predictor variables. The average number of maternal and early neonatal mortalities was 0.587 and 0.533, respectively. This indicates that the maternal mortality rate was higher than the early neonatal mortality rate in six districts in South Sulawesi Province, Indonesia. Based on the coefficient of variance (CoV) value, maternal mortality was more homogeneous than early neonatal mortality because the CoV value of variable Y 1 was lower than that of variable Y 2 . Before applying the BZIPIGIR type II model, it is necessary to test the overdispersion assumption and calculate the percentage of zero values in the response variable. Overdispersion occurs when the variance of the response variable is larger than the mean. This hypothesis can be tested with the following formula [22]: Overdispersion occurs if the deviance value divided by the degrees of freedom is larger than one. Table 2 shows that all response variables are overdispersed because the value of D h /df is greater than one, and there are excess zeros, with more than 61.33% of all observations on Y 1 and 68% on Y 2 equaling zero. Therefore, the BIZIPIGR type II model should be used to model these data. Based on the results of modeling using the BZIPIGR type II model, the statistical value, G 2 = 1801.85, was greater than χ 2 0.95;24 = 36.415, so the conclusion was that at least one of the predictor variables influenced the number of maternal and early neonatal mortalities in South Sulawesi, Indonesia. The model and partial test statistics for each predictor are presented in Table 3. Based on the results of the partial test in Table 3, all predictor variables were found to have a p-value < 0.05, so the conclusion was that all predictor variables had a significant effect, both for the variable Y 1 and Y 2 , with a significance level of 0.05.

Modeling the Number of Maternal and Early Neonatal Mortalities Using the BZIPIGR Type II Model without Exposure Variables
Based on the results of modeling using the BZIPIGR type II model without exposure variables, the statistical value, G 2 = 776.754, was larger than χ 2 0.95;24 = 36.415, so the conclusion was that at least one of the predictor variables had an influence on the number of maternal and early neonatal mortalities in South Sulawesi, Indonesia. The model and partial test statistics for each predictor are presented seen in Table 4.
Based on the results of the partial test in Table 4, several predictor variables have a p-value > 0.05, namely the X 1 , X 2 , X 3 and X 4 variables for the variable Y 2 , so the conclusion was that the effect of these four variables on the number of early neonatal mortalities (Y 2 ) was not significant. The effect of all variables on the number of maternal mortalities (Y 1 ) was significant.

Comparison of the BZIPIGR Type II Model with and without Exposure Variables
The comparison of the BZIPIGR type II model with and without the exposure variables was performed based on the sum square error (SSE) value between the two models. The SSE values of these two models are presented in Table 5. The above real data analysis obtained similar results to the simulation data analysis. Data were generated as much as n = 75 and n = 100, with 50% of response variables containing zero values. Data were analyzed using the BZIPIGR type II model under two conditions: with and without exposure variables. The simulation SSE and MSE values are presented in Table 6. As shown in Tables 5 and 6, the SSE and RMSE values of the BZIPIGR type II model with exposure variables are less than those of the model without exposure variables. Therefore, the model with exposure variables performed better than the model without exposure variables. Tables 5 and 6, the model with exposure variables performed better than the model without exposure variables. This discussion focuses on the BZIPIGR type II model with exposure variables. There are two BZIPIGR model types II, the log model and the logit model. The log model states that the probability of Y i from the PIG state is significantly affected by a variable, whereas the logit model states that the probability of Y i from the zero state is significantly influenced by a variable. The two models for each response variable are as follows: a.

Based on a comparison of the SSE and RMSE values in
Regression model for the number of maternal mortalities: Regression model for the number of early neonatal mortalities: (17) Equations (13) and (14) show the model for the number of maternal mortalities. In this equation, if other variables are held constant, then every 1% increase in variable X 1 results in an increase in the number of maternal mortalities by 0.005 times, with the same explanation for variables X 3 and X 6 . For variable X 2 , if the other variables are held constant, then every 1% increase in variable X 2 results in a decrease in the number of maternal mortalities by 0.006 times, with the same explanation for variables X 4 and X 5 . However, for the logit model in Equation (14), the coefficient values of X 1 , X 3 and X 6 are positive. The means that, if other predictor variables are held constant, with every increase in the percentage of visits by pregnant women, the probability that maternal mortality will not occur increases, with the same explanation for X 3 (percentage of births assisted by health workers) and X 6 (percentage of obstetric complications).
Based on the log and logit models in Equations (15) and (16), the predictor variable that affects the PIG state also affects the zero state (ZI). Based on the coefficient values of the predictor variables, the influence of the predictor variable on the response variable can be interpreted. For example, X 2 (percentage of pregnant women who received Fe3 tablets) in Equation (15) has a positive value, which means that the number of early neonatal mortalities increases with the percentage of pregnant women who received Fe3 tablets. Variables X 4 (percentage of mother who attended at least three postpartum maternal visits), X 5 (percentage of active integrated service posts) and X 6 (percentage of obstetric complications) increase with the number of early neonatal mortalities. That means that if the other predictor variables are held constant, then every 1% increase in the X 2 variable (percentage of births assisted by health workers) increases the number of maternal mortalities by 0.066 times the number of early neonatal mortalities, with the same explanation for variables X 4 , X 5 and X 6 . For the logit model in Equation (16), where the coefficient value of the X 1 is negative, every increase in the percentage of visits by pregnant women reduces the likelihood of early neonatal mortalities not occurring if other predictor variables are held constant, with the same explanations for X 2 , X 4 and X 6 .

Conclusions
The BZIPIGR model is derived from the Poisson regression model and can be used when the equidispersion assumption is violated and excess zeros are present. A violation of this assumption occurs when asymmetry exists between the mean and variance of the response variable, especially when the variance is larger than the mean (overdispersion). The BZIPIGR type I model only accommodates the possibility that all response variables are zero or that all response variables are non-zero. Therefore, the BZIPIGR type II model was developed to overcome other possibilities of the zero value of the response variable, namely that one of the response variables is zero. The BZIPIGR type II model accommodates all possible combinations of zero values in the response variable. Thus, the BZIPIGR type II model accounts for bivariate cases that experience overdispersion and excess zeros, such as maternal and early neonatal mortalities in South Sulawesi Province, Indonesia, in 2019.
According to the SSE and RMSE values of real data in Table 5, the BZIPIGR type II model with exposure variables performed better than the model without exposure variables in estimating the parameter values. This conclusion can be generally accepted based on the results of the analysis using the simulation data in Table 6, where the same results are obtained, namely lower SSE and RMSE values in the BZIPIGR type II model with exposure variables. The existence of the exposure variables in the BZIPIGR type II model strongly affects the model accuracy. The empirical results of modeling the number of maternal and early neonatal mortalities at the sub-district level in the province of South Sulawesi, Indonesia with exposure variables using the BZIPIGR type II model showed that all predictor variables used in this study, i.e., percentage of visits by pregnant women, percentage of pregnant women who received Fe3 tablets, percentage of births assisted by health workers, percentage of mothers who attended at least three postpartum maternal visits, percentage of active integrated service posts and percentage of obstetric complications, had a significant effect on the response variable.
The limited availability of data at the sub-district level means it is possible that there are other predictor variables that were not included in the study. Therefore, for further research, the BZIPIGR type II model can be applied to the same case by replacing or adding other variables or other cases that meet the assumptions of this model. The existence of differences in the characteristics of each region can cause differences in the incidence in each region; therefore, further research using the BZIPIGR type II model can add spatial effects to the model. Author Contributions: Conceptualization, E.E., P.P. and S.P.R.; methodology, E.E., P.P. and S.P.R.; software, E.E. and P.P.; validation, E.E., P.P. and S.P.R.; formal analysis, E.E., P.P. and S.P.R.; investigation, E.E., P.P. and S.P.R.; data curation, E.E.; writing-original draft preparation, E.E.; writing-review and editing, P.P. and S.P.R.; supervision, P.P. and S.P.R.; visualization, E.E. and S.P.R.; project administration, P.P. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The joint probability function of Y 1 and Y 2 can be written as in the following equation: The First Derivative of Likelihood functions for parameters of the BZIPIGR type II model. The likelihood function of the BZIPIGR type II is T is the vector parameters of the BZIPIGR model. Therefore, the log-likelihood function in Equation (10) is derived from each parameter and equated with zero, as follows: The first derivative of the first log likelihood (log L(•)) function to parameter λ 1 is ; B 32 = 2 πτ 1 2 y −1 1i !; 1i exp x T i λ 1 y 1i + τ −1 ; The first derivative of the log L(•) function to parameter λ 2 is obtained as The first derivative of the log-likelihood function to parameter α 1 is obtained as The first derivative of the log-likelihood function to parameter α 2 is The first derivative of the first log-likelihood function (log L(•)) to parameter τ is obtained as 1i exp x T i λ 1 y 1i exp 1 τ τ −2 ; The first derivative of the likelihood for the model without exposure variables is like the BZIPIGR type II model with exposure variables, but q i is not involved in the equation above.