## 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, [PERSON]. [20] developed the Geographically Weighted Bivariate Poisson Inverse Gaussian Regression (GWBPIGR) model. To accommodate a combination of global and local parameter estimates, [PERSON].