## 2 Specifications of Geographically Weighted Multivariate Poisson Inverse Gaussian Regression
The following statistical test used is:  
\[G^{2}=2\bigg{(}\ell\bigg{(}\overset{\Lambda}{\mathbf{\Omega}}\bigg{)}-\ell \bigg{(}\overset{\Lambda}{\mathbf{\omega}}\bigg{)}\bigg{)} \tag{18}\]  
where \(\ell\bigg{(}\overset{\Lambda}{\mathbf{\Omega}}\bigg{)}\) is the log-likelihood function for a parameter set under population (\(\mathbf{\Omega}=\{\beta_{1}(u_{i},v_{i}),\beta_{2}(u_{i},v_{i}),\ldots,\beta_{m }(u_{i},v_{i}),\tau(u_{i},v_{i});i=1,2,\ldots,n\}\)) and \(\ell\bigg{(}\overset{\Lambda}{\mathbf{\omega}}\bigg{)}\) is the log-likelihood function for a parameter set under the null hypothesis \(\mathbf{\omega}=\{\beta_{01\omega}(u_{i},v_{i}),\beta_{02\omega}(u_{i},v_{i}), \ldots,\beta_{0m\omega}(u_{i},v_{i}),\tau_{\omega}(u_{i},v_{i});i=1,2,...,n\}\). The log-likelihood function of GWMPIGR model for parameter set under population is:  
\[ \begin{array}{rl}\ell\bigg{(}\overset{\Lambda}{\mathbf{\Omega}}\bigg{)}&= \frac{n}{2}\log\big{(}\frac{2}{n}\big{)}+\sum\limits_{i=1}^{n}\frac{1}{\tau(u_ {i},v_{i})}-\frac{1}{2}\sum\limits_{i=1}^{n}\log\hat{\tau}(u_{i},v_{i})+\sum \limits_{i=1}^{n}\log(K_{s_{i}}(z_{i}(u_{i},v_{i})))\\ &\qquad\qquad-\sum\limits_{i=1}^{n}\left(\frac{2\sum\limits_{j=1}^{n}y_{ij}- 1}{4}\right)\log\left(1+2\hat{\tau}(u_{i},v_{i})\sum\limits_{j=1}^{m}q_{i}e^{ \omega_{i}T^{\Lambda}\beta_{j}(u_{i},v_{i})}\right)+\\ &\qquad\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}y_{ij}\log(q_{i})+\sum \limits_{i=1}^{n}\sum\limits_{j=1}^{m}y_{ij}\mathbf{x}_{i}^{\top}\overset{ \Lambda}{\mathbf{\beta}}_{j}(u_{i},v_{i})-\left(\sum\limits_{i=1}^{n}\sum \limits_{j=1}^{m}\log(y_{ij}!)\right)\end{array}  \tag{19}\]  
The log-likelihood function of GWMPIGR model for parameter set under \(\mathrm{H}_{0}\) is:  
\[ \begin{array}{rl}\ell\bigg{(}\overset{\Lambda}{\mathbf{\omega}}\bigg{)}&= \frac{n}{\tau_{\hat{\omega}}(u_{i},v_{i})}+\frac{n}{2}\log\Big{(}\frac{2}{\tau _{\hat{\omega}}(u_{i},v_{i})}\Big{)}+\sum\limits_{i=1}^{n}\log(K_{s_{i}}(z_{ \hat{\omega}}(u_{i},v_{i})))-\sum\limits_{i=1}^{n}\frac{\left(2\sum\limits_{ j=1}^{m}y_{ij}-1\right)}{4}\times\\ &\log\bigg{(}1+2\tau_{\hat{\omega}}(u_{i},v_{i})\sum\limits_{j=1}^{m}q_{i}e^{ \beta_{j\hat{\omega}}(u_{i},v_{i})}\bigg{)}+\sum\limits_{i=1}^{n}\sum\limits_{ j=1}^{m}y_{ij}\log(q_{i})+\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}y_{ij}\hat{ \beta}_{j0\hat{\omega}}(u_{i},v_{i})-\\ &\qquad\qquad\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\log(y_{ij}!)\end{array}  \tag{20}\]  
The estimator value for the parameter set under \(\mathrm{H}_{0}\) are obtained using the same steps as those for estimating the parameters of the GWMPIGR model. The \(G^{2}\) test statistic (Equation (18)) has a \(\chi^{2}_{trace(\mathbf{R})}\), where \(\mathrm{trace}(\mathbf{R})\) is the number of effective parameters for the GWMPIGR model, in which the elements of the matrix \(\mathbf{R}\) are formulated as follows:  
\[\mathbf{r}_{i}=\mathbf{X}_{i}\Big{(}\mathbf{X}^{T}\mathbf{W}(u_{i^{*}},v_{i^{* }})\mathbf{A}_{(t)}(u_{i^{*}},v_{i^{*}})\mathbf{X}\Big{)}^{-1}\mathbf{X}^{T} \mathbf{W}(u_{i^{*}},v_{i^{*}})\mathbf{A}_{(t)}(u_{i^{*}},v_{i^{*}})\]  
The critical area for testing the regression parameter hypothesis can be determined for the GWMPIGR model simultaneously by rejecting \(\mathrm{H}_{0}\) if the value of \(G^{2}>\chi^{2}_{\alpha,trace(\mathbf{R})}\)[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:  
\[AIC_{c}=-2L(\theta_{MPIGR})+2k^{*}+\frac{2k^{*}(k^{*}+1)}{n-k^{*}-1} \tag{21}\]  
where \(k^{*}\) is the number of effective parameters in the GWMPIGR model, and \(L(\theta_{MPIGR})\) is the likelihood function of the GWMPIGR model (Equation (5)).