## 2 Specifications of Geographically Weighted Multivariate Poisson Inverse Gaussian Regression
The joint probability mass function of \(Y_{i1}\), \(Y_{i2}\),..., \(Y_{im}\) with \(j=1,2,\ldots,m\), \(i=1,2\),..., \(n\  
\[=\prod\limits_{i=1}^{n}\left(e^{\frac{1}{\tau(u_{i^{*}},v_{i^{*}})}} \left(\frac{2}{\pi\tau(u_{i^{*}},v_{i^{*}})}\right)^{\frac{1}{2}}K_{s_{i}}(z_{i} (u_{i^{*}},v_{i^{*}}))\right.\] \[\left(1+2\tau(u_{i^{*}},v_{i^{*}})\sum\limits_{j=1}^{m}q_{i}e^{ \mathbf{x}_{i}^{T}\mathbf{\beta}_{j}(u_{i^{*}},v_{i^{*}})}\right)^{-\frac{(2 \sum\limits_{j=1}^{n}v_{i^{*}})}{4}}\prod\limits_{j=1}^{m}\frac{\left(q_{i}e^{ \mathbf{x}_{i}^{T}\mathbf{\beta}_{j}(u_{i^{*}},v_{i^{*}})}\right)^{v_{ij}}}{y_{ ij}!}\left.\right)\]  
The log-likelihood function of Equation (6) is as follows:  
\[\ell^{*}(\Theta_{i^{*}}) =\log\!\left\{\prod\limits_{i=1}^{n}\left(P(y_{i1},\ldots,y_{imi} ;\mathbf{\beta}_{1}(u_{i^{*}},v_{i^{*}}),...,\mathbf{\beta}_{m}(u_{i^{*}},v_{i^ {*}}),\tau(u_{i^{*}},v_{i^{*}}))\right)^{w_{ii*}}\right\} \tag{7}\] \[=\sum\limits_{i=1}^{n}w_{ii^{*}}\log(P(y_{i1},\ldots,y_{imi}; \mathbf{\beta}_{1}(u_{i^{*}},v_{i^{*}}),...,\mathbf{\beta}_{m}(u_{i^{*}},v_{i ^{*}}),\tau(u_{i^{*}},v_{i^{*}})))\] \[=\sum\limits_{i=1}^{n}\left(\frac{w_{ii^{*}}}{\tau(u_{i^{*}},v_{ i^{*}})}+\frac{1}{2}\log\left(\frac{2}{\pi\tau(u_{i^{*}},v_{i^{*}})}\right)w_{ii^{ *}}+\log K_{s_{i}}(z(u_{i^{*}},v_{i^{*}}))w_{ii^{*}}-\right.\] \[\left.\left(\frac{2\sum\limits_{j=1}^{n}y_{ij}-1}{4}\right)\log \!\left(1+2\tau(u_{i^{*}},v_{i^{*}})\sum\limits_{j=1}^{m}q_{i}e^{\mathbf{x}_{ i}^{T}\mathbf{\beta}_{j}(u_{i^{*}},v_{i^{*}})}\right)w_{ii^{*}}+\sum\limits_{j=1}^{m}y_ {ij}\log(q_{i})w_{ii^{*}}+\right.\] \[\left.\sum\limits_{j=1}^{m}y_{ij}\mathbf{x}_{i}^{T}\mathbf{\beta }_{j}(u_{i^{*}},v_{i^{*}})w_{ii^{*}}-\sum\limits_{j=1}^{m}\log(y_{ij}!)w_{ii^{ *}}\right)\]  
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*}(u_{i},v_{i})=\exp\!\left(-\frac{1}{2}\!\left(\frac{d_{ii*}}{h}\right)^ {2}\right) \tag{8}\]  
Meanwhile, the fixed bisquare kernel function is  
\[w_{ii*}(u_{i},v_{i})=\left\{ \begin{array}{ll}\left(1-\left(d_{ii*}/h\right)^ {2}\right)^{2},&\text{untuk }\;d_{ii*}\leq h\\ 0,&\text{untuk }\;d_{ii*}>h\end{array} \right. \tag{9}\]  
with \(d_{ii*}=\sqrt{\left(u_{i}-u_{i*}\right)^{2}+\left(v_{i}-v_{i*}\right)^{2}}\) is the euclidean distance between location \((u_{i},v_{i})\) and location \((u_{i*},v_{i*})\) 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.