## 2 Methodology  
### Poisson-Inverse-Gaussian Regression Model (PIGM)  
[PERSON] [25] proposed the PIG distribution as a mixture of the Poisson and inverse-Gaussian distributions. The resulting probability mass function is defined as follows:  
\[\mathrm{p}(\mathrm{Y}=\mathrm{y}|\mu,\omega)=\sqrt{\frac{2}{\pi\omega}\frac{e^ {(\mu)^{-1}}}{y!}}\Bigg{(}\sqrt{2\omega\bigg{(}1+\frac{1}{2\omega\mu^{2}} \bigg{)}}\Bigg{)}^{-(y-0.5)}K_{y-0.5}\Bigg{(}\sqrt{\frac{2}{\omega}\bigg{(}1+ \frac{1}{2\omega\mu^{2}}\bigg{)}}\Bigg{)}, \tag{1}\]  
for \(\mathrm{y}=0,\,1,2,\,\ldots,\,\mu>0\) is the mean parameter, \(\omega>0\) is the dispersion parameter, \(K_{y-0.5}(. )\) represents the modified Bessel function of the second kind [26]. The mean and variance are \(E(y)=\mu\) and \(V(y)=\mu+\frac{\mu^{3}}{\omega}\), respectively. The linear predictor is defined as  
\[\mu_{i}=\exp\Big{(}x_{i}^{\prime}\beta\Big{)}\hskip 14.226378 pti=1,2,\ldots,n, \tag{2}\]  
where \(x_{ij}\) is the matrix of the explanatory variables and \(\beta\) is the vector of regression parameters. The maximum likelihood estimator (MLE) method is widely employed for estimating regression parameters due to its desirable statistical properties, such as consistency, efficiency, and asymptotic normality. This method determines the parameter values that maximize the log-likelihood function, identifying the values that best explain the observed data. For the regression model specified in Equation (1), the corresponding log-likelihood function is formulated as follows:  
\[\mathcal{L}(\beta,\omega\mid y_{i})=\sum_{i=1}^{n}\Bigg{[}\frac{1}{2}\mathrm{ ln}\bigg{(}\frac{2}{\pi\omega}\bigg{)}+(\omega\mu)^{-1}-\ln(y_{i}!)-(y-0.5) \Bigg{(}\sqrt{2\omega\bigg{(}1+\frac{1}{2\omega\mu^{2}}\bigg{)}}\Bigg{)}+K_{y- 0.5}\Bigg{(}\sqrt{\frac{2}{\omega}\bigg{(}1+\frac{1}{2\omega\mu^{2}}\bigg{)}} \Bigg{)}\Bigg{]} \tag{3}\]The partial derivative of the log-likelihood function is given as:  
\[S(\beta)=\frac{\partial\epsilon(\beta,\omega\mid y_{i})}{\partial\beta_{j}}=0 \tag{4}\]  
The solution to Equation (4) requires numerical methods such as the Newton-Raphson method and Fisher score method [26; 27]. The Fisher scoring algorithm is commonly employed for estimating the parameters \(\beta_{j},j=1,2,\ldots,q\). Let \(\hat{\beta}^{l}\) denotes the maximum likelihood estimator after \(l\) iterations, then \(\hat{\beta}^{l}\) is given as follows:  
\[\hat{\beta}^{l+1}=\hat{\beta}^{l}+\Big{(}I\Big{(}\beta^{l}\Big{)}\Big{)}^{-1}S \Big{(}\beta^{l}\Big{)} \tag{5}\]  
where \(l\) is the number of iterations until convergence is reached and \(I\Big{(}\beta^{l}\Big{)}=-E\Big{(}\frac{\hat{\beta}^{l}l(\beta)}{\hat{\beta} \hat{\beta}_{l}}\Big{)}\) is the Fisher information matrix. The iterations terminate when \(\hat{\beta}^{l}-\hat{\beta}^{l+1}<\epsilon\) where \(\epsilon\) is a specified small value and usually equals \(10^{-6}\). The final estimate algorithm form is then combined with iterative reweighted least squares as follows, with a few modifications:  
\[\hat{\beta}^{PIG-MLE}=\big{(}X^{\prime}\hat{W}X\big{)}^{-1}\big{(}X^{\prime} \hat{W}\phi\big{)} \tag{6}\]  
where \(\hat{W}\) is a \((n\times n)\) diagonal matrix with its diagonal elements \(\hat{w}_{i}=\hat{\mu}_{i}+\frac{\hat{\mu}_{i}^{3}}{\omega}\), \(i=1,2,\ldots,n\) and \(\phi=\ln\hat{\mu}_{i}+\frac{\hat{w}_{i}-\hat{\mu}_{i}}{\hat{\mu}_{i}+\frac{\hat {w}_{i}^{3}}{\omega}}\). The dispersion parameter \(\omega\) is defined as follows:  
\[\hat{\omega}=\frac{1}{n-p}{\sum_{i=1}^{n}\frac{(y_{i}-\hat{\mu}_{i})^{2}}{Var (\hat{\mu}_{i})}}\]  
In this study, the matrix mean squared error (MMSE) and scalar mean squared error (SMSE) are primarily used to evaluate the performance of the estimator \(\widetilde{\tau}\).