The Geographically Weighted Multivariate Poisson Inverse Gaussian Regression Model and Its Applications

Abstract

This study aims to develop a method for multivariate spatial overdispersion count data with mixed Poisson distribution, namely the Geographically Weighted Multivariate Poisson Inverse Gaussian Regression (GWMPIGR) model. The parameters of the GWMPIGR model are estimated locally using the maximum likelihood estimation (MLE) method by considering spatial effects. Therefore, the significance of the regression parameter differs for each location. In this study, four GWMPIGR models are evaluated based on the exposure variable and the spatial weighting function. We compare the performance of those four models in real-world application using data on the number of infant, under-5 and maternal deaths in East Java in 2019 using five predictor variables. In this study, the GWMPIGR model uses one exposure variable and three exposure variables. Compared to the fixed kernel Gaussian weighting function, the GWMPIGR model with the fixed kernel bisquare weighting function and one exposure variable has a better fit based on the AICc value. Furthermore, according to the best GWMPIGR model, there are several regional groups formed based on predictors that significantly affected each event in East Java in 2019.

1. Introduction

This study aims to develop a method for multivariate spatial overdispersion count data. In general, count data can be analyzed using Poisson regression models. However, the assumption of equidispersion in Poisson regression is difficult to fulfill. On the other hand, variance that exceeds the mean (overdispersion) is frequently found when analyzing real data [1,2]. Therefore, alternative methods are needed to model data with overdispersion conditions. Mixed Poisson models are often used as an alternative to Poisson regression models, such as the negative binomial (NBR) regression model [3,4,5] and Poisson-inverse Gaussian regression (PIGR) model [6,7,8]. In this study, the PIGR model was chosen because this model performs better when modeling data with high overdispersion. In this study, the PIGR model was developed into a multivariate model with two or more response variables [9,10].

The PIGR model is a global regression model that assumes each observation location is influenced by the same predictor. However, in some cases, effects of location cannot be ignored because each location has different characteristics, such as geography and culture, among others. Other spatial regression models, specifically a point-based spatial regression model for overdispersed count data, have been studied. The Geographically Weighted Bivariate Generalized Poisson Regression (GWBGPR) model was used by [11] to analyze the factors affecting the amount of infant and maternal mortality in East Java, Indonesia. Parameter estimation and hypothesis testing of the Geographically Weighted Multivariate Generalized Poisson Regression (GWMGPR) model has been studied by [12], theoretically. Both studies use the maximum likelihood estimation (MLE) method and Newton-Raphson (NR) iterative algorithm to estimate the model parameters. The NR methods require the second derivatives of the log-likelihood function with respect to each parameter in the model.

Generalized Poisson regression deals with problem related to over- or underdispersion. Meanwhile, the negative binomial regression (NBR) and PIGR models only deal with overdispersion. The Geographically Weighted Negative Binomial Regression (GWNBR) model was studied by [13], who applied it to simulation data and real data to show the superiority of the model compared to global regression models, such as the Poisson and negative binomial regression models. The iteratively reweighted least squares (IRLS) and NR methods are often used in parameter estimation. The GWNBR model was also used by [14] to model the confirmed COVID-19 cases in East Java. The model uses the adaptive bisquare Kernel weighting function and the parameter estimation is performed using a combination of the IRLS and NR methods. The results showed that COVID-19 spread quickly in locations with a high population density. The authors of [15] evaluated the Geographically Weighted Multivariate Negative Binomial Regression (GWMNBR) in a multivariate way, and the parameters were estimated by combination of the MLE and NR methods. The results showed that the GWMNBR method performed better than the global method.

As mentioned before, the PIGR method only deals with overdispersion, and according to [1] and other previous empirical studies [7,8,16], the PIGR model has larger variance than the NBR model. Therefore, the PIGR model can accommodate greater overdispersion than the NBR model can. The PIGR model has been developed into a multivariate model with two or more response variables [10,17]. Globally, the bivariate Poisson inverse Gaussian regression (BPIGR) model was developed by [9], while the multivariate Poisson inverse Gaussian regression (MPIGR) was developed by [10,18]. Locally, the Geographically Weighted Poisson Inverse Gaussian Regression (GWPIGR) model was created by [19]. For bivariate cases, Ref. [20] developed the Geographically Weighted Bivariate Poisson Inverse Gaussian Regression (GWBPIGR) model. To accommodate a combination of global and local parameter estimates, Ref. [21] developed an alternative GWBPIGR model, namely the Mixed Geographically Weighted Bivariate Poisson Inverse Gaussian Regression (MGWBPIGR) model. In all those studies, parameter estimation was carried out using the MLE and NR methods and the spatial kernel weighting function.

The purpose of this study was to develop a spatial based-point model for multivariate count data with PIG distribution, namely the Geographically Weighted Multivariate Poisson Inverse Gaussian Regression (GWMPIGR) model. The proposed model was developed based on the MPIGR model, which was developed by [10]. There are two main differences in the model that we developed in this study compared to those that have been developed in previous studies. First, by considering the probability of each occurrence in each research unit, several exposures were used for the case study in this study. We have previously studied univariance, via case study, using the PIGR model with exposure, which performed better than the PIGR model without exposure [9]. Second, we used the iterative Berndt-Hall-Hall-Hausman (BHHH) methods for parameter estimation. Compared to the NR method, the BHHH method does not need the second derivative of the log-likelihood function with respect to each parameter in the model. This second derivative has the potential to cause problems such as the Hessian matrix and may not be positively definite. It also typically results in a complex nonlinear function that needs to be solved, meaning that the expected value remains unknown. Therefore, the second derivative was eliminated, and a simpler computation with guaranteed convergence was obtained, the BHHH method [22,23].

For the real-world application of the proposed model, we analyzed one of the health problems in Indonesia, namely the high number of infant, under-5 and maternal deaths. Numerous studies have been conducted to produce effective policies to reduce the mortality rate of the three groups of individuals. This problem is also one of the targets in the third goal of the SDGs for 2030, living a healthy and prosperous life, which can partially be achieved by reducing the neonatal mortality rate to 12 per 1000 live births and the under-5 mortality rate by 25 per 1000 children under five in the middle of the same year. Furthermore, the maternal mortality rate should also be reduced to less than 70 per 100,000 live births.

Indonesia, the fourth most populous country in the world, is struggling to achieve the 2030 SDG targets for maternal and under-5 mortality. Based on the Indonesian Demographic and Health Survey (IDHS) in 2017, the Maternal Mortality Rate (MMR) was 305 per 100,000 live births. In 2020, the Infant Mortality Rate (IMR) in Indonesia was 17.6 per 1000 live births and the under-5 Mortality Rate (U5MR) was 23 per 1000 live births [24,25]. East Java is one of the cities in Indonesia with the highest number of infant, under-5 and maternal deaths compared to other regions in Indonesia. Therefore, it is necessary to study the causal factors to reduce infant, under-5 and maternal mortality in East Java.

In terms of the factors that affect infant, under-5 and maternal mortality, various aspects of life such as socioeconomic status, health facilities, environment, education and culture play equally important roles. Several factors are often the closest factor or the direct cause of the death. These factors include maternal factors, nutrition and disease control in infant and children under five, while bleeding and hypertension are the main causes of death in mothers [26,27]. In 2018, Ref. [28] conducted a study on the factors of neonatal mortality in Indonesia and found that low birth weight (LBW), giving birth to children too close together, medical personnel and postnatal/postpartum visits have a significant effect on neonatal mortality in Indonesia.

In this study, we use four GWMPIGR models, the exposure variable and the spatial weighting function. This study aims to compare on the performance of each model in modeling the number of infant, under-5 and maternal deaths and determine the factors that cause those events based on the best model. District/city level data with different characteristics (spatial heterogeneity) are used. Thus, the form of functions and parameters between locations are divergent. As a result, the parameter significance is also different for each location. Through the GWMPIGR model, districts/cities in Java can be grouped based on predictor variables that significantly affect the number of infant, under-5 and maternal deaths. This grouping is expected to help the relevant agencies make policies to solve this issue.

The specifications of the GWMPIGR model are described in Section 2, the real-world application of the GWMPIGR model is described in Section 3 and a discussion of the factors affecting the number of infant, under-5 and maternal deaths in East Java in 2019 based on the selected model will be discussed in Section 4.

2. Specifications of Geographically Weighted Multivariate Poisson Inverse Gaussian Regression

The MPIGR model is a model that was expanded from the PIGR model, which is a GLMs model with a logarithmic link function. Let (Y₁, Y₂, …, Y_m) ~ MPIG (μ₁, μ₂, …, μ_m), and the joint probability mass function is as follows:

$$P\left(y_1,y_2,\dots ,y_m\right)={\left(2z{\pi }^{-1}\right)}^{\frac{1}{2}}{e}^{\frac{1}{\tau }}{K}_s(z){\left(z\tau \right)}^{{}^{-{\displaystyle \sum _{j=1}^m{y}_j}}}{\displaystyle \prod _{j=1}^m\frac{{\left({\mu }_j\right)}^{y_j}}{y_j!}}$$

(1)

with y_j = 0, 1, 2, … and $\tau$ is the dispersion parameter $\left(\tau >0\right)$, $s={\displaystyle \sum _{j=1}^m{y}_j}-\frac{1}{2}$, $z=\frac{1}{\tau }\sqrt{1+2\tau {\displaystyle \sum _{j=1}^m{\mu }_j}}$, and $K_s\left(z\right)=K_{\left({\displaystyle \sum _{j=1}^m{y}_j}\right)-\frac{1}{2}}\frac{1}{\tau }\sqrt{1+2\tau {\displaystyle \sum _{j=1}^m{\mu }_j}}$ is the third modification of the Bessel function [8]. The properties of the PIGD can be found in [29,30]. The MPIGR model can be described as

$$E\left(Y_j\right)={\mu }_j=q_j\exp \left(x^T{\beta }_j\right)$$

(2)

where j = 1, 2, …, m, ${\mu }_j$ is the predicted mean of the jth response variable, $q_j$ is the exposure for the jth response variable, and τ is the overdispersion parameter. ${\beta }_j$ is the regression parameter that corresponds to the predictor variable k, for k = 1, 2, …, p [10,18].

The GWMPIGR model is based on MPIGR model when it incorporates spatial or location-based aspects into the model. The GWMPIGR model is a point-based spatial model where the regression parameters depend on the geographic location. Therefore, each location has different functions and parameters. The GWMPIGR model for the jth response variable can be written as follows:

$${\mu }_j\left(u_i,v_i\right)=q_j{e}^{x_{}^T{\beta }_j\left(u_i,v_i\right)}$$

(3)

where $\left(u_i,v_i\right)$ refer to the longitude and latitude coordinates of the ith location, ${\beta }_j\left(u_i,v_i\right)={\left[\begin{array}{ccc}{\beta }_{j0}\left(u_i,v_i\right)& {\beta }_{j1}\left(u_i,v_i\right)& \begin{array}{cc}...& {\beta }_{jp}\left(u_i,v_i\right)\end{array}\end{array}\right]}^T$ is a regression parameter vector of the jth response variable for location i with the following dimension (p + 1) × 1.

For the GWMPIGR model, parameter estimation is performed using the maximum likelihood estimation (MLE) method and the Berndt-Hall-Hall-Hausman (BHHH) iteration method. At each observation location (point), parameter estimation is carried out using a spatial weighted matrix (W). The elements of matrix W show how much influence the observations at a particular location have on the surrounding observations (w_ii*).

The joint probability mass function of Y_i1, Y_i2, …, Y_im with j = 1, 2, …, m, i = 1, 2, …, n and $\left(u_i,v_i\right)$ coordinates is as follows:

$$\begin{array}{l}P\left(y_{i1},y_{i2},\dots ,y_{im}\right)\\ ={\left(\frac{2}{\pi \tau \left(u_i,v_i\right)}\right)}^{\frac{1}{2}}{e}^{\frac{1}{{\tau }_i}}{K}_{s_i}\left(z\left(u_i,v_i\right)\right){\left(1+2\tau \left(u_i,v_i\right){\displaystyle \sum _{j=1}^m{\mu }_j\left(u_i,v_i\right)}\right)}^{-\frac{\left(2{\displaystyle \sum _{j=1}^m{y}_{ij}}-1\right)}{4}}{\displaystyle \prod _{j=1}^m\frac{{\mu }_j{\left(u_i,v_i\right)}^{y_{ij}}}{y_{ij}!}}\end{array}$$

(4)

with $s_i={\displaystyle \sum _{j=1}^m{y}_{ij}}-\frac{1}{2}$, $z\left(u_i,v_i\right)=\frac{1}{\tau \left(u_i,v_i\right)}\sqrt{1+2\tau \left(u_i,v_i\right){\displaystyle \sum _{j=1}^m{q}_i{e}^{x_i^T{\beta }_j\left(u_i,v_i\right)}}}$, $K_{s_i}\left(z\left(u_i,v_i\right)\right)=K_{{\displaystyle \sum _{j=1}^m{y}_{ij}-\frac{1}{2}}}\left(\frac{1}{\tau \left(u_i,v_i\right)}\sqrt{1+2\tau \left(u_i,v_i\right){\displaystyle \sum _{j=1}^m{q}_i{e}^{x_i^T{\beta }_j\left(u_i,v_i\right)}}}\right)$, j = 1, 2, …, m, $y_{ij}>0$, and $\tau >0$.

Let ${\mathsf{\theta }}_i={\left[\begin{array}{ccccc}{\beta }_1^T\left(u_i,v_i\right)& {\beta }_2^T\left(u_i,v_i\right)& \dots & {\beta }_m^T\left(u_i,v_i\right)& \tau \left(u_i,v_i\right)\end{array}\right]}^T$ is the vector of parameter the GWMPIGR model. Then, the likelihood function for the GWMPIGR model is

$$L\left({\mathsf{\theta }}_i,i=1,2,...,n\right)={\displaystyle \prod _{i=1}^{n}P\left(y_{i1},\dots ,y_{im};{\beta }_1\left(u_i,v_i\right),...,{\beta }_m\left(u_i,v_i\right),\tau \left(u_i,v_i\right)\right)}$$

$$\begin{array}{l}={\displaystyle \prod _{i=1}^n({\left(\frac{2}{\pi }\right)}^{\frac{1}{2}}{e}^{\frac{1}{\tau \left(u_i,v_i\right)}}{\left(\frac{1}{\tau \left(u_i,v_i\right)}\right)}^{\frac{1}{2}}{K}_{s_i}\left(z_i\left(u_i,v_i\right)\right)}\times \\ {\left(1+2\tau \left(u_i,v_i\right){\displaystyle \sum _{j=1}^m{q}_i{e}^{x_i^T{\beta }_j\left(u_i,v_i\right)}}\right)}^{-\frac{\left(2{\displaystyle \sum _{j=1}^m{y}_{ij}}-1\right)}{4}}{\displaystyle \prod _{j=1}^m\frac{{\left(q_i{e}^{x_i^T{\beta }_j\left(u_i,v_i\right)}\right)}^{y_{ij}}}{y_{ij}!}})\end{array}$$

(5)

Let ${\mathsf{\theta }}_{i^{*}}={\left[\begin{array}{ccccc}{\beta }_1^T\left(u_{i^{*}},v_{i^{*}}\right)& {\beta }_2^T\left(u_{i^{*}},v_{i^{*}}\right)& \dots & {\beta }_m^T\left(u_{i^{*}},v_{i^{*}}\right)& \tau \left(u_{i^{*}},v_{i^{*}}\right)\end{array}\right]}^T$ is the vector of parameter for location i*. Then, the likelihood function to estimate the parameter at location i* with the spatial weighted matrix w_ii* is as follows:

$$L\left({\mathsf{\theta }}_{i^{\ast }}\right)={\displaystyle \prod _{i=1}^n{\left(P\left(y_{i1},y_{i2},\dots ,y_{im};{\beta }_1\left(u_{i^{\ast }},v_{i^{\ast }}\right),...,{\beta }_m\left(u_{i^{\ast }},v_{i^{\ast }}\right),\tau \left(u_{i^{\ast }},v_{i^{\ast }}\right)\right)\right)}^{w_{ii*}}}$$

$$\begin{array}{l}={\displaystyle \prod _{i=1}^n(e^{\frac{1}{\tau \left(u_{i^{\ast }},v_{i^{\ast }}\right)}}{\left(\frac{2}{\pi \tau \left(u_{i^{\ast }},v_{i^{\ast }}\right)}\right)}^{\frac{1}{2}}{K}_{s_i}\left(z_i\left(u_{i^{\ast }},v_{i^{\ast }}\right)\right)}\\ {\left(1+2\tau \left(u_{i^{\ast }},v_{i^{\ast }}\right){\displaystyle \sum _{j=1}^m{q}_i{e}^{x_i^T{\beta }_j\left(u_{i^{\ast }},v_{i^{\ast }}\right)}}\right)}^{-\frac{\left(2{\displaystyle \sum _{j=1}^m{y}_{ij}}-1\right)}{4}}{{\displaystyle \prod _{j=1}^m\frac{{\left(q_i{e}^{x_i^T{\beta }_j\left(u_{i^{\ast }},v_{i^{\ast }}\right)}\right)}^{y_{ij}}}{y_{ij}!}})}^{w_{ii*}}\end{array}$$

(6)

The log-likelihood function of Equation (6) is as follows:

$$\begin{array}{ll}{\ell }^{\ast }\left({\mathsf{\theta }}_{i^{*}}\right)& =\log \left\{{\displaystyle \prod _{i=1}^n{\left(P\left(y_{i1},\dots ,y_{im};{\beta }_1\left(u_{i^{*}},v_{i^{*}}\right),...,{\beta }_m\left(u_{i^{*}},v_{i^{*}}\right),\tau \left(u_{i^{*}},v_{i^{*}}\right)\right)\right)}^{w_{ii*}}}\right\}\\ & ={\displaystyle \sum _{i=1}^n{w}_{i{i}^{\ast }}}\log \left(P\left(y_{i1},\dots ,y_{im};{\beta }_1\left(u_{i^{*}},v_{i^{*}}\right),...,{\beta }_m\left(u_{i^{*}},v_{i^{*}}\right),\tau \left(u_{i^{*}},v_{i^{*}}\right)\right)\right)\end{array}$$

$$\begin{array}{l}={\displaystyle \sum _{i=1}^n(\frac{w_{i{i}^{\ast }}}{\tau \left(u_{i^{\ast }},v_{i^{\ast }}\right)}+\frac{1}{2}\log \left(\frac{2}{\pi \tau \left(u_{i^{\ast }},v_{i^{\ast }}\right)}\right)w_{i{i}^{\ast }}+\log K_{s_i}\left(z\left(u_{i^{\ast }},v_{i^{\ast }}\right)\right)w_{i{i}^{\ast }}-}\\ \left(\frac{2{\displaystyle \sum _{j=1}^m{y}_{ij}}-1}{4}\right)\log \left(1+2\tau \left(u_{i^{\ast }},v_{i^{\ast }}\right){\displaystyle \sum _{j=1}^m{q}_i{e}^{x_i{}^T{\beta }_j\left(u_{i^{\ast }},v_{i^{\ast }}\right)}}\right)w_{i{i}^{\ast }}+{\displaystyle \sum _{j=1}^m{y}_{ij}\log \left(q_i\right)}{w}_{i{i}^{\ast }}+\\ {\displaystyle \sum _{j=1}^m{y}_{ij}{x}_i{}^T{\beta }_j\left(u_{i^{\ast }},v_{i^{\ast }}\right)}{w}_{i{i}^{\ast }}-{\displaystyle \sum _{j=1}^m\log \left(y_{ij}!\right)w_{i{i}^{\ast }}})\end{array}$$

(7)

Several methods can be used to determine the spatial weights for each different location in the GWMPIGR model, including the kernel function. There are two types of the kernel functions: the fixed kernel function and the adaptive kernel function. The difference between those two functions lies in the bandwidth value (h). The bandwidth value for the fixed kernel function is the same for every location (h). Meanwhile, the bandwidth value for the adaptive kernel function is different for every location (h_i). In this study, we used fixed kernel functions, specifically the fixed Gaussian kernel functions and fixed bisquare kernel functions. The fixed Gaussian kernel function is formulated as follows:

$$w_{ii*}\left(u_i,v_i\right)=\exp \left(-\frac{1}{2}{\left(\frac{d_{ii*}}{h}\right)}^2\right)$$

(8)

Meanwhile, the fixed bisquare kernel function is

$$w_{ii*}\left(u_i,v_i\right)=\{\begin{array}{l}{\left(1-{\left(d_{ii*}/h\right)}^2\right)}^2,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{untuk}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{d}_{ii*}\le h\\ 0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em} \mathrm{untuk}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{d}_{ii*}>h\end{array}$$

(9)

with $d_{ii*}=\sqrt{{\left(u_i-u_{i*}\right)}^2+{\left(v_i-v_{i*}\right)}^2}$ is the eucliden distance between location $\left(u_i,v_i\right)$ and location $\left(u_{i*},v_{i*}\right)$ and h is a smoothing parameter (bandwidth).

The bandwidth value is related to the accuracy of the model. A very small bandwidth value causes the variance to become larger. On the other hand, a large bandwidth value can cause a larger bias. Therefore, it is very important to choose the appropriate bandwidth value. According to [31], the selection of the optimal bandwidth value can be achieved using the cross validation (CV) criteria, which can be formulated mathematically as follows:

$$CV\left(h\right)=\min \left\{{\displaystyle \sum _{i=1}^n{\left[y_i-{\widehat{y}}_{\left(\ne i\right)}(h)\right]}^T\left[y_i-{\widehat{y}}_{\left(\ne i\right)}(h)\right]}\right\}$$

(10)

The ML estimator that maximizes the log-likelihood function of the GWMPIGR model in Equation (7) is obtained by solving the system of equations for all of the first partial derivatives of the log-likelihood function for each parameter and then equating it with zero as follows:

$$\frac{\partial {\ell }^{\ast }\left({\mathsf{\theta }}_{i^{*}}\right)}{\partial {\mathsf{\theta }}_{i^{*}}}=0$$

(11)

The first-order partial derivatives of the log-likelihood function are as follows:

$$\frac{\partial \ell \left({\mathsf{\theta }}_{i^{\ast }}\right)}{\partial {\beta }_1\left(u_{i^{\ast }},v_{i^{\ast }}\right)}={\displaystyle \sum _{i=1}^n\left[y_{i1}-M\left(y_{i1},y_{i2},\dots ,y_{im}\right)\exp \left(x_i^T{\beta }_1\left(u_{i^{\ast }},v_{i^{\ast }}\right)\right)\right]x_i{}^T{w}_{i{i}^{\ast }}}$$

(12)

$$\frac{\partial \ell \left({\mathsf{\theta }}_{i^{\ast }}\right)}{\partial {\beta }_2\left(u_{i^{\ast }},v_{i^{\ast }}\right)}={\displaystyle \sum _{i=1}^n\left[y_{i2}-M\left(y_{i1},y_{i2},\dots ,y_{im}\right)\exp \left(x_i^T{\beta }_2\left(u_{i^{\ast }},v_{i^{\ast }}\right)\right)\right]x_i{}^T{w}_{i{i}^{\ast }}}$$

(13)

$$\frac{\partial \ell \left({\mathsf{\theta }}_{i^{*}}\right)}{\partial {\beta }_m\left(u_{i^{\ast }},v_{i^{\ast }}\right)}={\displaystyle \sum _{i=1}^n\left[y_{im}-M\left(y_{i1},y_{i2},\dots ,y_{im}\right)\exp \left(x_i^T{\beta }_m\left(u_{i^{\ast }},v_{i^{\ast }}\right)\right)\right]x_i{}^T{w}_{i{i}^{\ast }}}$$

(14)

$$\begin{array}{ll}\frac{\partial \ell \left({\mathsf{\theta }}_{i^{\ast }}\right)}{\partial \tau \left(u_{i^{\ast }},v_{i^{\ast }}\right)}& =\frac{1}{{\tau }^2\left(u_{i^{\ast }},v_{i^{\ast }}\right)}{\displaystyle \sum _{i=1}^n(M\left(y_{i1},y_{i2},\dots ,y_{im}\right)}\times \\ & \left(1+\tau \left(u_{i^{\ast }},v_{i^{\ast }}\right){\displaystyle \sum _{j=1}^m\exp \left(x_i^T{\beta }_j\left(u_{i^{\ast }},v_{i^{\ast }}\right)\right)}\right)-1-\frac{1}{\tau \left(u_{i^{\ast }},v_{i^{\ast }}\right)}{\displaystyle \sum _{j=1}^m{y}_{ij}})w_{ii*}\end{array}$$

(15)

where

$$M\left(y_{i1},y_{i2},\dots ,y_{im}\right)=\frac{1}{\sqrt{1+2\tau \left(u_{i^{\ast }},v_{i^{\ast }}\right){\displaystyle \sum _{j=1}^m{q}_i\exp \left(x_i^T{\beta }_j\left(u_{i^{\ast }},v_{i^{\ast }}\right)\right)}}}\frac{K_{{\displaystyle \sum _{j=1}^m{y}_{ij}}+\frac{1}{2}}\left(z\left(u_{i^{\ast }},v_{i^{\ast }}\right)\right)}{K_{{\displaystyle \sum _{j=1}^m{y}_{ij}}-\frac{1}{2}}\left(z\left(u_{i^{\ast }},v_{i^{\ast }}\right)\right)}$$

(16)

The first derivative of the log-likelihood function of the GWMPIGR model parameters (Equations (12)–(15)) is non-closed-form equation. Thus, the estimator value is obtained through the Berndt-Hall-Hall-Hausman (BHHH) iteration method using the following algorithm:

Step-1

Determine the initial value of each parameter using the results of the MPIGR model parameter estimator.

Step-2

Define the gradient vector $g\left({\mathsf{\theta }}_i\right)$ with the first derivative (Equation (11)) as the element.

Step-3

Define the Hessian matrix $H\left({\mathsf{\theta }}_i\right)$,

$$H\left({\mathsf{\theta }}_i\right)=-\left[{\displaystyle \sum _{i=1}^n{g}_i\left(\mathsf{\theta }\right){\left(g_i\left(\mathsf{\theta }\right)\right)}^T}\right]$$

where $g_i\left(\mathsf{\theta }\right)$ is the individual gradient vector.

Step-4

Define the tolerance limits (ε = 10⁻³) and maximum iteration (t^* = 1000).

Step-5

Start the BHHH iteration using the following equation:

$${\widehat{\mathsf{\theta }}}_i^{\left(t+1\right)}={\widehat{\mathsf{\theta }}}_i^{\left(t\right)}-H^{-1}\left({\widehat{\mathsf{\theta }}}_i^{\left(t\right)}\right)g\left({\widehat{\mathsf{\theta }}}_i^{\left(t\right)}\right),$$

(17)

Step-6

The iteration will stop at the t* iteration or at the value of $\Vert {\widehat{\mathsf{\theta }}}_i^{(t+1)}-{\widehat{\mathsf{\theta }}}_i^{(t)}\Vert \le \epsilon .$

Repeat this algorithm for each location i (i = 1, 2, …, n).

The simultaneous hypothesis testing of the GWMPIGR model parameters is carried out by testing H₀: ${\beta }_{j1}\left(u_i,v_i\right)={\beta }_{j2}\left(u_i,v_i\right)=...={\beta }_{jk}\left(u_i,v_i\right)=\dots ={\beta }_{jp}\left(u_i,v_i\right)=0$ versus H₁: at least one ${\beta }_{jk}\left(u_i,v_i\right)\ne 0$, where j = 1, 2, …, m; k = 1, 2, …, p; and i = 1, 2, …, n. The following statistical test used is:

$$G^2=2\left(\ell \left(\widehat{\Omega }\right)-\ell \left(\widehat{\omega }\right)\right)$$

(18)

where $\ell \left(\widehat{\Omega }\right)$ is the log-likelihood function for a parameter set under population $\left(\Omega =\left\{{\beta }_1^{}\left(u_i,v_i\right),{\beta }_2^{}\left(u_i,v_i\right),\dots ,{\beta }_m^{}\left(u_i,v_i\right),\tau \left(u_i,v_i\right);i=1,2,\dots ,n\right\}\right)$ and $\ell \left(\widehat{\omega }\right)$ is the log-likelihood function for a parameter set under the null hypothesis $\omega =\left\{{\beta }_{01\omega }\left(u_i,v_i\right),{\beta }_{02\omega }\left(u_i,v_i\right),\dots ,{\beta }_{0m\omega }\left(u_i,v_i\right),{\tau }_{\omega }\left(u_i,v_i\right);i=1,2,...,n\right\}$.

The log-likelihood function of GWMPIGR model for parameter set under population is:

$$\begin{array}{ll}\ell \left(\widehat{\Omega }\right)& =\frac{n}{2}\log \left(\frac{2}{\pi }\right)+{\displaystyle \sum _{i=1}^n\frac{1}{\widehat{\tau }\left(u_i,v_i\right)}}-\frac{1}{2}{\displaystyle \sum _{i=1}^n\log \widehat{\tau }\left(u_i,v_i\right)}+{\displaystyle \sum _{i=1}^n\log \left(K_{s_i}\left(z_i\left(u_i,v_i\right)\right)\right)}\\ & -{\displaystyle \sum _{i=1}^n\left(\frac{2{\displaystyle \sum _{j=1}^m{y}_{ij}}-1}{4}\right)}\log \left(1+2\widehat{\tau }\left(u_i,v_i\right){\displaystyle \sum _{j=1}^m{q}_i{e}^{x_i{}^T{\widehat{\beta }}_j\left(u_i,v_i\right)}}\right)+\\ & {\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^m{y}_{ij}\log \left(q_i\right)}}+{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^m{y}_{ij}{x}_i{}^T{\widehat{\beta }}_j\left(u_i,v_i\right)}}-\left({\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^m\log \left(y_{ij}!\right)}}\right)\end{array}$$

(19)

The log-likelihood function of GWMPIGR model for parameter set under H₀ is:

$$\begin{array}{ll}\ell \left(\widehat{\omega }\right)& =\frac{n}{{\widehat{\tau }}_{\widehat{\omega }}\left(u_i,v_i\right)}+\frac{n}{2}\log \left(\frac{2}{\pi {\widehat{\tau }}_{\widehat{\omega }}\left(u_i,v_i\right)}\right)+{\displaystyle \sum _{i=1}^n\log \left(K_{s_i}\left(z_{\widehat{\omega }}\left(u_i,v_i\right)\right)\right)}-{\displaystyle \sum _{i=1}^n\frac{\left(2{\displaystyle \sum _{j=1}^m{y}_{ij}}-1\right)}{4}\times }\\ & \log \left(1+2{\tau }_{\widehat{\omega }}\left(u_i,v_i\right){\displaystyle \sum _{j=1}^m{q}_i{e}^{{\widehat{\beta }}_{j0\widehat{\omega }}\left(u_i,v_i\right)}}\right)+{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^m{y}_{ij}\log \left(q_i\right)}}+{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^m{y}_{ij}{\widehat{\beta }}_{j0\widehat{\omega }}\left(u_i,v_i\right)}}-\\ & {\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^m\log \left(y_{ij}!\right)}}\end{array}$$

(20)

The estimator value for the parameter set under H₀ are obtained using the same steps as those for estimating the parameters of the GWMPIGR model.

The G² test statistic (Equation (18)) has a ${\chi }_{trace\left(R\right)}^2$, where trace(R) is the number of effective parameters for the GWMPIGR model, in which the elements of the matrix R are formulated as follows:

$$r_i=X_i{\left(X^{T}W\left(u_{i^{*}},v_{i^{*}}\right)A_{\left(t\right)}\left(u_{i^{*}},v_{i^{*}}\right)X\right)}^{-1}{X}^{T}W\left(u_{i^{*}},v_{i^{*}}\right)A_{\left(t\right)}\left(u_{i^{*}},v_{i^{*}}\right)$$

The critical area for testing the regression parameter hypothesis can be determined for the GWMPIGR model simultaneously by rejecting H₀ if the value of $G_{}^2>{\chi }_{\alpha ,trace\left(R\right)}^2$ [32,33].

This study will compare several models based on the exposure variable and the spatial weighting function used. The exposure variable in the GWMPIGR model in Equation (3) can consist of the same or different exposure levels for all of the response variables. Therefore, the best model will be selected using the corrected AIC (AICc) value:

$$AI{C}_c=-2L\left({\mathsf{\theta }}_{MPIGR}\right)+2k^{*}+\frac{2k^{*}\left(k^{*}+1\right)}{n-k^{*}-1}$$

(21)

where k* is the number of effective parameters in the GWMPIGR model, and $L\left({\mathsf{\theta }}_{MPIGR}\right)$ is the likelihood function of the GWMPIGR model (Equation (5)).

3. Results

This section shows the fit of the GWMPIGR model on the number of infant, under-5 and maternal deaths in East Java in 2019. The data were sourced from the Public Health Office and Statistics Indonesia. The Health Office of Indonesia provides data on health indicators from the provincial and district/city levels, and these data are published regularly every year. Meanwhile, Statistics Indonesia provides data for all indicators (including health), and these data are published annually. District/city level data are provided by the Public Health Office and Statistics Indonesia for every province in Indonesia.

In this study, the data consisted of 38 districts/cities in the province of East Java Province, and there were three response variables: the number of infant deaths (Y₁), the number of under-5 deaths (Y₂) and the number of maternal deaths (Y₃); five predictor variables: the percentage of active integrated service post (X₁), the percentage of active family planning participants (X₂), the percentage of the population with BPJS health insurance (X₃), education index (X₄) and the percentage of household that has improved sanitation (X₅); and three exposure variables sourced from the Public Health Office and Statistics Indonesia. Four models will be used based on the exposure and the spatial weighting function, as presented in

Table 1. Model Specifications.

Model	Number of Exposure	Exposure	Spatial Weighting Function
1	1	q: The number of live births	Gaussian Kernel fixed function
2			Bisquare Kernel fixed function
3	3	q1: The number of live births q2: The number of people aged 1-4 years old q3: The number of pregnant women	Gaussian Kernel fixed function
4			Bisquare Kernel fixed function

.

The first two GWMPIGR models only had one exposure variable, the number of live births, which was used by considering the definitions of the mortality rate for each individual group. Meanwhile, the third and fourth models use three exposure variables by considering the number of people who are at risk in each individual group. The best of the four models will be selected based on the AICc value. Then, the factors affecting the number of infant, under-5 and maternal deaths in East Java in 2019 are discussed based on the best model.

Figure 1, Figure 2 and Figure 3 are thematic maps that provide descriptions of the response variables. The number of deaths in each population is divided into four levels from the lowest to the highest number of deaths in the district/city.

Figure 1 and Figure 3 show that most of the districts/cities in East Java have a high number of infant and maternal deaths. Meanwhile, Figure 2 shows that the Probolinggo District and Surabaya City have the highest number of under-5 deaths. Furthermore, the summary statistics of the predictor variables and exposure variables can be seen in

Table 2. Summary statistics of predictor variables and exposure variables.

Variable	Mean	Standard Deviation	Minimum	Maximum	VIF
The percentage of active integrated service post (X1)	80.95	14.91	31.00	98.00	1.871
The percentage of active family planning participants (X2)	75.40	3.39	67.03	82.97	1.186
The percentage of the population with BPJS health insurance (X3)	50.47	10.70	32.83	74.87	1.998
Education index (X4)	0.628	0.075	0.490	0.770	2.999
The percentage of household that has improved sanitation (X5)	72.20	18.87	25.51	97.43	3.104
The number of live birth (q1)	15,240	10,370	2012	44,378
The number of under-5 population (q2)	14,787	10,008	2079	41,646
The number of pregnant women (q3)	59,252	39,641	8013	168,060

.

According to Table 2, the mean of the active integrated service post in East Java in 2019 was 80.95%, with a minimum percentage of 31%. Furthermore, the mean of active family planning participants is 75.4%, with a minimum percentage of 67.03%. This means around 33% of people do not participate in family planning. The mean of the population with BPJS health insurance is 50.47%, which means that the other half of population uses other types of health insurances. The mean of the education index is 0.628. The education index is calculated based on the expected years of schooling and the mean years of schooling. The higher the value of the education index, the better the quality of education and society, which will consequently improve the quality of life. As for the environment variables, there was a large gap in the percentage of households with improved sanitation, with the minimum percentage being 25.25%. This means that there are still many households with inadequate sanitation.

Furthermore, a multicollinearity test, which assumes that there is no relationship among predictor variables (mutually independent), was carried out as a requirement for the regression analysis. In this study, the variance inflation factor (VIF) criteria was used to identify multicollinearity among predictor variables. The last column in Table 2 shows that the VIF value of all of the predictor variables is lower than ten (VIF < 10). Thus, there are no multicollinearity issues among the predictor variables.

An initial examination of the relationship between the response and the predictor variables is important because it relates to the initial description of the relationships between variables found in the modeling. Based on Figure 4, the relationship between the predictor variable and the response variable seems non-linear. Log(Y₁/q₁) has the strongest correlation with X₃ (−0.23). Meanwhile, log(Y₂/q₂) has the strongest correlation with X₁ (−0.36) and log(Y₃/q₃) has the strongest correlation with X₂ (−0.29).

Several assumptions must be met before carrying out GWMPIGR modelling. Specifically, the response variable is not Poisson distributed, and there is spatial heterogeneity. Testing the distribution of the response variables was carried out using the Crockett test, and it was found that the response variables did not have a trivariate Poisson distribution (because there were three responses). This result is also supported by the issue of overdispersion on the response variables. The methods used to test the overdispersion are deviance per degree of freedom (deviance/df) and the Lagrange multiplier (LM) test, the test results of which are shown in Table 2.

Overdispersion exists if the value of deviance/df is greater than one and the LM method produces a value of ${\chi }^2>Z_{\alpha }^{}$ or a p-value of <α. Based on

Table 3. Overdispersion test.

Variable	Deviance/df	LM
		χ2	p-Value
Y1	38.78	211,028.5	0.000 *
Y2	8.13	1336.02	0.000 *
Y3	5.44	3532.65	0.000 *

* Significant at α = 5%.

, there are cases of overdispersion on the response variable. Therefore, the GWMPIGR model can be used to model data on the number of infant, under-5 and maternal deaths in East Java in 2019. Furthermore, the spatial heterogeneity test was carried out by using the Glejser test and it was found that there is spatial heterogeneity when the test statistics value G = 210.5993 is larger than ${\chi }_{0.05;15}^2=24.996$.

The parameter of the GWMPIGR model were estimated locally using spatial weighting so that each observation location had a different parameter estimate value. Spatial weighting represents the location between observations. The closer an observation location is to another observation location, the greater the weighting value and the greater the influence on these observations. In this study, spatial weighting was determined using a fixed Gaussian and bisquare kernel function. The optimum bandwidth, minimum CV and AICc values for each model are presented in

Table 4. Bandwidth optimum, CV minimum and AICc for each model.

Model	Bandwidth Optimum	CV Minimum	AICc
1	0.975566	505,601.5	64,309.95
2	2.043404	496,850.6	64,445.95
3	0.880499	529,872.7	74,930.46
4	1.834768	524,106.5	74,962.95

: The best model.

.

Based on Table 4, model 2 has the lowest CV minimum value. However, this does not mean that model 2 is the best model. The CV value is used to determine the optimum bandwidth for the spatial weight matrix. Therefore, the CV values in Table 4 refer to the minimum CV value of the optimum bandwidth for each model. On the other hand, the AICc can be used to evaluate how well a model fits the data and to determine the best model from multiple models for the same dataset. AICc uses the MLE of the model (log-likelihood) as a measure of fit. The model that fits the data the best has the maximum likelihood. Therefore, the model with high log-likelihood has a low AICc value. Based on Table 4, model 1 has the smallest AICc value. Thus, in this study, the GWMPIGR model with the fixed Gaussian kernel spatial weighting function and the number of live births as the exposure is the best model to determine the number of infant, under-5 and maternal deaths in East Java in 2019. Therefore, only model 1 will be discussed further.

The value of the parameter estimator of the GWMPIGR model is different for each observation location, as is the significance of the parameters. The GWMPIGR model generated 722 coefficient estimates for 38 districts/cities in the province of East Java. Simultaneous testing of the GWMPIGR model parameters shows that at least one parameter has a significant effect on the response variable (G² = 384.17 > ${\chi }_{0.05;12}^2=21.026$). Meanwhile, partial hypothesis testing was carried out on the parameters of the GWMPIGR model locally for each observation’s location. For example, the parameter estimates for the GWMPGR model with the fixed Gaussian kernel spatial weighted function for Banyuwangi District (i = 10) and Surabaya City (i = 37) can be seen in

Table 5. Parameter estimate of the GWMPGR model with the fixed Gaussian kernel spatial weighted function for Banyuwangi District (i = 10) and Surabaya City (i = 37).

Parameter	Banyuwangi District (i = 10)			Surabaya City (i = 37)
	Estimate	Standard Error	Z Value	Estimate	Standard Error	Z Value
${\beta }_{01} \left(u_i,v_i\right)$	−4.828	0.007138	−676.437 *	−4.729	0.074794	−63.238 *
${\beta }_{11} \left(u_i,v_i\right)$	0.007	0.002722	2.689 *	0.010	0.001115	8.738 *
${\beta }_{21} \left(u_i,v_i\right)$	0.015	0.004518	3.397 *	0.009	0.001725	5.119 *
${\beta }_{31} \left(u_i,v_i\right)$	0.008	0.007047	1.110	0.015	0.001777	8.597 *
${\beta }_{41} \left(u_i,v_i\right)$	−3.686	0.008304	−443.894 *	−3.680	0.061363	−59.976 *
${\beta }_{51} \left(u_i,v_i\right)$	0.0006	0.003648	0.175	−0.003	0.001202	−2.139 *
${\beta }_{02} \left(u_i,v_i\right)$	1.242	0.00081	1532.867 *	1.249	0.005383	232.158 *
${\beta }_{12} \left(u_i,v_i\right)$	−0.025	0.015424	−1.595	−0.042	0.003436	−12.198 *
${\beta }_{22} \left(u_i,v_i\right)$	−0.159	0.037767	−4.225 *	−0.114	0.007923	−14.398 *
${\beta }_{32} \left(u_i,v_i\right)$	0.048	0.062046	0.769	−0.002	0.010041	−0.176
${\beta }_{42} \left(u_i,v_i\right)$	7.193	0.000808	8902.213 *	7.195	0.005621	1280.048 *
${\beta }_{52} \left(u_i,v_i\right)$	−0.020	0.018594	−1.083	−0.011	0.003638	−2.918 *
${\beta }_{03} \left(u_i,v_i\right)$	0.046	0.00162	28.508 *	0.059	0.009337	6.349 *
${\beta }_{13} \left(u_i,v_i\right)$	0.014	0.014956	0.962	0.009	0.006989	1.405
${\beta }_{23} \left(u_i,v_i\right)$	−0.142	0.020068	−7.092 *	−0.123	0.006156	−19.954 *
${\beta }_{33} \left(u_i,v_i\right)$	0.023	0.028089	0.826	−0.014	0.01074	−1.298
${\beta }_{43} \left(u_i,v_i\right)$	3.652	0.001021	3574.817 *	3.653	0.007974	458.116 *
${\beta }_{53} \left(u_i,v_i\right)$	−0.010	0.018	−0.569	0.0003	0.00712	0.004
$\tau \left(u_i,v_i\right)$	41.709	0.000228	182,747 *	41.706	0.002148	19,415.9 *

* Significant at α = 5%.

.

Table 5 shows that there are differences in the parameter coefficient values and in the parameter significance for Banyuwangi District and Surabaya City. Some of the predictors that were significant in Surabaya City were not significant in Banyuwangi District. Variable X₃ and X₅ did not significantly affect Y₁ in the Banyuwangi District, while all of the predictors significantly affected Y₁ in Surabaya City. Variables X₁, X₃ and X₅ did not significantly affect Y₂ in the Banyuwangi District, while only variable X₃ did not significantly affect Y₂ in Surabaya City. Meanwhile, similar parameter significance was found for the variable Y₃ More specifically, variable Y₃ is equally influenced by variables X₂ and X₄ in the Banyuwangi District and Surabaya City. The significance of the parameters of the GWMPIGR model in other regencies/cities in the province of East Java that are different is that they can then be formed into several regional groups according to the predictors that significantly affect the number of infant, under-5 and maternal deaths. This grouping will be discussed further in the Discussion section.

To provide an example, the GWMPIGR model for Surabaya City based on Table 5 is presented as follows:

$$\begin{array}{l}{\widehat{\mu }}_1\left(u_{37},v_{37}\right)=q_{37}^{}\exp \left(-4.729+0.010X_{1,37}+0.009X_{2,37}+0.015X_{3,37}-3.680X_{4,37}-0.003X_{5,37}\right)\\ {\widehat{\mu }}_2\left(u_{37},v_{37}\right)=q_{37}^{}\exp \left(1.249-0.042X_{1,37}-0.114X_{2,37}-0.002X_{3,37}+7.195X_{4,37}-0.011X_{5,37}\right)\\ {\widehat{\mu }}_3\left(u_{37},v_{37}\right)=q_{37}^{}\exp \left(0.059+0.009X_{1,37}-0.123X_{2,37}-0.014X_{3,37}+3.653X_{4,37}+0.0003X_{5,37}\right)\end{array}$$

The value of the regression parameter coefficients in the model above shows the magnitude of the change in the average number of infant, under-5 and maternal deaths in Surabaya due to the influence of each predictor variable. The interpretation of the model will be discussed further in the discussion section.

4. Discussion

Model 1’s performance for Surabaya City was interpreted based on each predictor variable. Increasing the percentage of active integrated service posts (X₁) will increase the average number of infant deaths (Y₁), reduce the average number of under-5 deaths (Y₂) and increase the average number of maternal deaths (Y₃). However, X₁ does not significantly affect the number of maternal deaths (Y₃).

Increasing the percentage of active family planning participants (X₂) will increase the average number of infant deaths (Y₁), reduce the average number of under-5 deaths (Y₂) and reduce the average number of maternal deaths (Y₃). Increasing the population with BPJS health insurance (X₃) will increase the average number of infant deaths (Y₁), reduce the average number of under-5 deaths (Y₂) and reduce the average number of maternal deaths (Y₃). However, X₃ does not significantly affect the number of under-5 deaths (Y₂) or the number of maternal deaths (Y₃).

Increasing in education index (X₄) will reduce the average number of infant deaths (Y₁), increase the average number of under-5 deaths (Y₂) and increase the average number of maternal deaths (Y₃). Moreover, increasing the percentage of households with improved sanitation (X₅) will reduce the average number of infant deaths (Y₁), reduce the average number of under-5 deaths (Y₂) and increase the average number of maternal deaths (Y₃). However, X₅ does not significantly affect the number of maternal deaths (Y₃).

Table 6. The comparison of signs between correlation and regression parameters coefficient of the GWMPIGR model for Surabaya City.

Response Variable	Predictors	Cor (Y/q, X)	GWMPIGR Modeling
			X1	X2	X3	X4	X5
Y1	X1	+ *	+ *
	X2	+		+ *
	X3	+			+ *
	X4	- *				- *
	X5	-					- *
Y2	X1	- *	- *
	X2	-		- *
	X3	-			-
	X4	-				+ *
	X5	+					- *
Y3	X1	+	+ *
	X2	-		- *
	X3	+			-
	X4	+				+ *
	X5	+ *					+

* Significant at α = 5%; + Positive relationship; - Negative relationship.

shows the differences between correlation and GWMPIGR modelling for Surabaya City. Some predictors have a consistent relationship with correlation, and some do not. As we can see that the correlation between the number of under-5 deaths (Y₂) and education index (X₄) or the percentage of household that has improved sanitation (X₅) is not the same as that observed in the GWMPIGR modelling. Additionally, the correlation between the number of under-5 deaths (Y₃) and the percentage of the population with BPJS health insurance (X₃) is not consistent with the GWMPIGR modelling. However, X₃ shows that it is not significantly correlated with Y₃ and that it does not significantly affect Y₃. However, the pattern of the relationship between the response variable and the predictor variable is different for each location.

Furthermore, the results of the parameter testing for the GWMPIGR model 1 produced several regional groups based on significant predictors, which are presented in

Table 7. Regional groups based on significant predictors by the GWMPIGR model 1.

Response Variable	Number of Groups	District/City Order Number	Total Number District/City	Significant Predictor
Y1	3	6–8, 14–18, 25, 30, 32–35, 38	16	X1, X2, X3, X4, X5
		1–5, 9, 11–13, 19–24, 26–29, 31, 36	21	X1, X2, X3, X4
		10	1	X1, X2, X4
Y2	3	1–3, 9, 11–12, 19–21, 29, 36	11	X1, X2, X4, X5
		4–8, 13–18, 22–28, 30–35, 37–38	26	X1, X2, X4
		10	1	X2, X4
Y3	2	1–27, 30–38	36	X1, X2, X4
		28–29	2	X2, X4

.

There are three regional groups according to the significant predictors for the number of infant deaths and the number of under-5 deaths, while there are two regional groups based on significant predictors for the number of maternal deaths. For the number of infant deaths (Y₁), 16 districts/cities in East Java are affected by all of the predictor variables, 21 districts/cities are influenced by the percentage of active integrated service post (X₁), the percentage of active family planning participants (X₂), the percentage of the population with BPJS health insurance (X₃), education index (X₄) and only one district/city is affected by the percentage of active integrated service post (X₁), the percentage of active family planning participants (X₂) and education index (X₄).

For the number of under-5 deaths (Y₂), 11 districts/cities in East Java are affected by the percentage of active integrated service post (X₁), the percentage of active family planning participants (X₂), education index (X₄) and the percentage of households that have improved sanitation (X₅), 26 districts/cities are influenced by variables the percentage of active integrated service post (X₁), the percentage of active family planning participants (X₂) and education index (X₄) and only one district/city is affected by the percentage of active family planning participants (X₂) and education index (X₄). Moreover, for the number of maternal deaths (Y₃), 36 districts/cities in East Java are affected by the percentage of active integrated service post (X₁), the percentage of active family planning participants (X₂) and education index (X₄) and two districts/cities are affected by the variables the percentage of active family planning participants (X₂) and education index (X₄). Visually, the groupings are presented in the thematic maps shown in Figure 5, Figure 6 and Figure 7.

Based on Figure 5, Figure 6 and Figure 7, the GWMPIGR modeling forms several regional groups. Each group is assumed to have the same characteristics, so the predictors that significantly affect each event are also the same. Thus, this will make it easier for the government to make policies to decrease the number of infant, under-5 and maternal deaths. For example, based on Figure 7, the percentage of active family planning participants and the education index affect the number of maternal deaths in the Pamekasan and Sumenep districts. Therefore, those two districts should focus on increasing the number of people using family planning services and on improving the quality of education. The same applies to other regional groups.

5. Conclusions

The GWMPIGR model is a point-based spatial regression model where the parameter estimator is influenced by location effects. In this study, four GWMPIGR models were constructed based on exposure and spatial weighting functions: the fixed Gaussian kernel function and the Bisquare fixed kernel function. Based on the AICc value, the GWMPIGR model with the fixed Gaussian kernel weighting function and the number of live births as the exposure was better at modeling data on the number of infant, under-5 and maternal deaths in East Java in 2019. Several regional groups in East Java were formed based on predictors that significantly affect each event. However, further research, such as simulation studies, is needed to evaluate the proposed method.