# A Bayesian-Weighted Inverse Gaussian Regression Model with Application to Seismological Data

Authors: Ehsan Mesdaghi; Afshin Fallah; Rahman Farnoosh; Gholam Hossein Yari
Journal: Advances in Mathematical Physics
 (2022)

---

## Abstract

Developing an efficient model for analyzing right-skewed positive observations has a long history, and many authors have attempt in this direction. This is because the common analytic modeling procedures such as linear regression are often inappropriate for such data and leads to inadequate results. In this article, we proposed a new model for regression analysis of the right-skewed data by assuming the weighted inverse Gaussian, as a great flexible distribution, for response observations. In the proposed model, the complementary reciprocal of the location parameter of response variable is considered to be a linear function of the explanatory variables. We developed a fully Bayesian framework to infer about the model parameters based on a general noninformative prior structure and employed a Gibbs sampler to derive the posterior inferences by using the Markov chain Monte Carlo methods. A comparative simulation study is worked out to assess and compare the proposed model with other usual competitor models, and it is observed that efficiency is quite satisfactory. A real seismological data set is also analyzed to explain the applicability of the proposed Bayesian model and to access its performance. The results indicate to the more accuracy of proposed regression model in estimation of model parameters and prediction of future observations in comparison to its usual competitors in literature. Particularly, the relative prediction efficiency of the proposed regression model to the inverse Gaussian and log-normal regression models has been obtained to be 1.16 and 64, respectively, for the real-world example discussed in this paper.


---

## Body

## 1. Introduction

Modeling right-skewed positive observations has been considered repeatedly by many authors during the last decades; hence, there is a vast literature on different probability distributions utilized in order to take into account the uncertainty behind this type of data. Many distributions have been candidated and studied for this purpose from the exponential, gamma and Weibull to log-normal (LN), inverse Gaussian (IG), etc. The IG distribution, as a member of exponential distributions family, benefits from considerable flexibility along with good theoretical and computational properties. A look at the papers was considered that the IG distribution indicates to the desirability of this distribution in statistical literature (see, for example, [1–5]). The IG distribution has been recommended for analyzing right-skewed positive data as a serious competitor for some of the most well-known distributions such as log-normal and gamma. It is employed in various applications such as reliability (e.g., [2, 6]), social sciences (e.g., [7, 8]), marketing (e.g., [9]), engineering (e.g., [10]), industrial management (e.g., [11]), cardiology [12], and civil engineering (e.g., [13]). In addition, some extensions to the IG distribution have been proposed in order to provide more flexible and interpretable distributions. In this direction, [14] proposed a three-parameter extension of the IG distribution, named weighted inverse Gaussian (WIG) distribution with a great ability for modeling positive right-skewed data specially in the field of reliability and lifetime analysis. As pointed out by [15], WIG is a versatile lifetime distribution and a serious competitor for other lifetime distributions. The most important inferential aspects of WIG distribution have been considered by Gupta and Kundu. They proposed an expectation maximization (EM) algorithm to estimate the parameters of distribution.On the other hand, some authors have been discussed the regression analysis by assuming the IG as distribution of response observations. For example, [16] considered a simple regression model with a single explanatory variable and zero intercept. [17] suggested a log-normal approximation in exponential regression. [6, 18] proposed a model which offers physical interpretation in survival analysis. This model further has been used by [6] in analysis of variance and survey sampling, respectively. Up to our best knowledge, although the IG distribution has been considered vastly in regression analysis, there is no published research considering the WIG distribution for modeling purposes. In this paper, we developed a new model for regression analysis of positive right-skewed response observations under the WIG distribution. A fully Bayesian framework has been utilized to estimate the model parameters and predict the future observations. The proposed methodology employs a set of noninformative vague prior distributions from a well-known family of flexible distributions and proceeds to the posterior inferences by using a Gibbs sampler. The main advantages of the proposed model are as follows: (i) due to using the WIG as distribution of response variable, it has great flexibility for taking into account the variability of positively skewed response observations. This leads to more accurate estimations and predictions for the proposed model in comparison to its competitors in the literature. (ii) Despite serious difficulties related to the frequentist approach for estimating the parameters of the WIG distribution (see [15] for more details), the proposed model provides a straightforward Bayesian approach with a tractable posterior distribution and interpretable posterior inferences about the parameters. (iii) Due to the Bayesian nature of the proposed model, it is possible to incorporate personal prior beliefs about uncertainty behind the parameters in order to construct a more reliable model.The paper organized as follows: Section2 provides some background about the WIG distribution containing its most important distributional and inferential properties. In Section 3, a WIG regression model is proposed. The essential theory for the Bayesian analysis of the proposed model is developed in Section 4. In Sections 5.1 and 5.2, we provide empirical evidence for assessing the proposed methodology via analyzing simulated and real-world data, respectively. The paper is closed by Conclusion.


## 2. A Look at Weighted Inverse Gaussian Distribution

A random variable,Y, is said to have a WIG distribution with location parameter μ, scale parameter λ, and shape parameter δ, Y∼WIGμ,λ,δ, if its density function is given by
(1)fWIGy∣μ,λ,δ=wy;δ,μfIGy∣μ,λ,wherewy;δ,μ=δ+y/δ+μ and
(2)fIGy∣μ,λ=λ2πy3exp−λy−μ22yμ2show the density function of IG distribution. The shape of WIG distribution is always unimodal, and the mean and variance of distribution are given, respectively, by
(3)EY=μ+μ3λμ+δ,(4)VarY=μ3λ+μ5λμ+δ3−μμ+δ.As it can be seen from Equation (4), despite the Gaussian distribution, the mean and variance of WIG distribution are correlated. The cumulative distribution function (CDF) of WIG distribution can be written as
(5)FWIGy∣μ,δ,λ=Φηy+1−2μμ+δexp2λμΦξy,where Φ· denoted the CDF of standard normal distribution and ηy=λ1/2y−μ/μy1/2 and ξy=−λ1/2y−μ/μy1/2. The key point in derivation of these results is that the CDF of IG μ,σ2 distribution can be expressed as a linear combination of standard normal CDF of the form Fx=Φαx+exp2/σ2μΦβx, where αx=x−μ/σμx1/2 and βx=x+μ/σμx1/2are solutions of Euler differential equation 4x2y"+4xy′−y=0 which are of the form yx=c1x1/2+c2x−1/2, for constants c1 and c2. If we consider that a pair of linearity-independent solutions αx and βx correspond to IG distribution, then a function of the form Gx=AΦαx+BΦβx with chosen constants A and B is CDF of the WIG distribution. See [14] for more details.The WIG family becomes the IG family asδ⟶∞ and the reciprocal inverse Gaussian (RIG) family as δ⟶0, where the RIG, like the IG, is a special case of the generalized inverse Gaussian (GIG) distribution (see, for example, [19] for more details about GIG distribution).The problem of estimating the WIG distribution parameters has been discussed in the literature. [14] showed that the maximum likelihood (ML) estimators of the WIG distribution parameters can be obtained as responses of a nonlinear challenging three-dimensional optimization problem. [15] formulated the problem in a missing data structure and developed an EM algorithm for computing the ML estimates. They used asymptotic optimality properties of the ML estimators to construct asymptotic confidence intervals for the unknown parameters. We refer the readers to [14, 15, 20, 21] and other above-mentioned references for more details about the WIG distribution.Before concluding this short review section, it is useful to clarify two points. First, despite the IG distribution, the WIG distribution does not belong to the exponential family of distributions. Of course in some special cases, for example, for knownδ, it has a near relation to the exponential family of distributions.Second, the WIG distribution is a reparameterized version of the mixed inverse Gaussian (MIG) distribution. The MIG distribution defined as a two-component finite mixture of IG and length-biased inverse Gaussian (LBIG) distribution as(6)fMIGy∣μ,λ,δ=1−αfIGy∣μ,λ+αfLBIGy∣μ,λ,δ,where α∈0,1 is the mixing coefficient and
(7)fLBIGy∣μ,λ=y/μfIGy∣μ,λdenotes the pdf of LBIG distribution. From this point of view, the WIG density function obtained by adopting a reparameterization of the form α=μ/μ+δ in the MIG density function is given in Equation (6).


## 3. The WIG Regression Model

The traditional theory of regression modeling is based on normality assumption for the response observations. While in many applications, these assumptions may not be valid due to positive right-skewed structure of data. An elementary solution, in this situation, is to use the log-transformed data that provide the possibility of using standard inferences based on the normal distribution and leads to a log-normal regression model. But applying any transformation to data increases the complexity of the model building procedure and makes meanings of the parameters less clear. In these situations, the WIG family of distributions, due to its flexibility and desired distributional properties, can be a suitable candidate for data analysis.Consider a set of paired observations given byyi,xi,i=1,⋯,n, where for a given i, i=xi1,⋯,xip′ denotes the vector of explanatory variables and yi indicates the corresponding random response variable. Assuming WIGμ,λ,δ as distribution of the response observations, the WIG regression model is given by
(8)yi∣xi∼WIGμxi,λ,δ,μi=μxi=1xi′β,i=1,⋯,n,where β=β1,⋯,βp′ denotes the vector of regression coefficients. The likelihood function associated with the model (8) is given by
(9)Lμ,λ,δ∣y=∏i=1nδ+yiδ+μifIGyi∣μI,λ,∝λnexp−λ2∑i=1nyi−μi2yiμi2∏i=1nδ+yiδ+μi.LettingY=diagy1,⋯,yn and X=x1′,⋯,xn′′=xijn×p, it is easy to show that the following identities hold:
(10)∑i=1n∑j=1pyijxijβj′xijβj=βX′YyXβ,∑i=1n∑j=1pxijβj=1n′Xβ,∑i=1nμi−1=1n′Xβ,where 1n is a n-dimensional vector with all elements equal to one. Therefore, the summation which appeared in the exponential function in the right-hand side of Equation (9) can be written as
(11)∑i=1nyi−μi2yiμi2=∑i=1nyiμi−2−2μi−1+yi−1=∑i=1nxi′βyixi′β−2xi′β+yi−1=∑i=1n∑j=1pxijβjyi∑j=1pxijβj−2∑j=1pxijβj+yi−1=Qβ,where Qβ=YXβ−1n′Y−1YXβ−1n denotes a quadratic form in terms of parameter β. Hence, the likelihood and log-likelihood functions of the model are given by
(12)Lβ,λ,δ∣y∝λnexp−λ2Qβ∏i=1nδ+yiδ+xi′β−1,ℓβ,λ,δ∣y∝nlogλ−λ2Qβ+∑i=1nlogδ+yiδ+xi′β−1,respectively. Equating the partial derivations of the log-likelihood function with respect to the model parameters to zero leads to the following system of equations:
(13)∂∂δℓβ,δ,λ∣y=∑i=1nyi+xi′β−1δ+yiδ+xi′β−1=0,∂∂βℓβ,δ,λ∣y=∑i=1n−∂/∂βxi′β−1δ+yiδ+xi′β−12−12∂∂βQβ=0,∂∂λℓβ,δ,λ∣y=nλ−12Qβ=0,where Qβ is a quadratic form in terms of parameter β introduced previously. As it can be seen, we have two nonlinear equations in terms of parameters β and δ and a linear equation in terms of parameter λ. Therefore, the ML estimator of parameter λ is obtained in the closed form to be λ^=2n/Qβ, and there is no closed form expression for the ML estimators of parameters β and δ. Hence, a usual iterative optimization procedure like the Newton-Raphson could be employed to calculate the ML estimates of the model parameters. It should be pointed out that it is possible to use the other frequentist point estimation approaches such as the method of moments (MM) to estimate the model parameters. See [15] for more details about difficulties related to the estimation of parameters in WIG distribution by using the frequentist approaches. In what follows, we developed a Bayesian approach for this purpose.Before going further, we tend to give some explanation about the link function of the proposed model. As it can be seen, here we have adopted an inverse form link function to specify the relation of the explanatory variables and the location parameter of the response observations. Although, as it is mentioned previously, there is no background on using WIG distribution in modeling context, the link function given in (8) has been previously used by some authors for modeling purposes using the IG distribution. For example, [6, 22, 23] used the inverse form link function for regression analysis and analysis of covariance (ANCOVA) under the IG distribution, respectively. Of course, there are other possible choices for link function of this regression model. For example, [24] used a linear link function as μxi=xi′β, and [25] utilized an exponential form link function, μxi=expxi′β, in modeling of censored observations by IG distribution in order to guarantee the positivity of response variable mean. On the other hand since IG, despite WIG, is a member of exponential family with canonical parameter μ−2, one can use a link function of the form μxi=xi′β−1/2 for modeling purposes based on the IG distribution.Regardless of differences between the IG and WIG distributions, each of the above-mentioned link functions can be considered, probably with some modifications, for modeling process with the WIG distribution, provided that the estimated value of vector of regression coefficients should not result in a negative value for the mean of response observations. However, from the theoretical point of view, there are two reasons supporting this choice for link function. First, the canonical parameter in the inverse Gaussian family of distribution appears in the form of a function of inverse mean. The second reason is related to the mathematical computations. This choice enables us to provide posterior distribution in a well-formed mathematical expression in terms of a quadratic form of regression coefficients and gives us closed form expressions for some of posterior full-conditional distributions.


## 4. Bayesian Analysis of the Model

In order to provide a Bayesian framework for analysis of the proposed model, it is necessary to set suitable prior distributions for the model parameters. These prior distributions should be able to explain the prior beliefs of data analyst about the model parameters.


### 4.1. Prior Setting

In regression analysis, the main interested parameter that should be estimated is the vector of regression coefficient,β=β1,⋯,βp. Theoretically, any multivariate continuous distribution on ℝp could be employed as a prior distribution for this parameter. For example, in a vague or a noninformative Bayesian analysis, a flat multivariate normal or a multivariate uniform distribution with large variance is a usual choice for this purpose. But, in the proposed model, considering such prior distributions for β may lead to the nonpositive values for μi=1/xi′β, i=1,⋯,n, which is inconsistent with the positivity constrain for the location parameter of WIG distribution. Taking these considerations into account, we consider the following joint prior distribution of the model parameters as
(14)τi∣λ∼Gamma2,λ,λ∼Gammac0,d0,δ∼Gammae0,f0,where τi=μi+δ and zero indices indicate the hyperparameters of the model that should be determined based on prior beliefs. As it can be seen, in this prior setting, all of the prior knowledge about the model parameters is expressed by using three gamma distributions. The family of gamma distributions seems to be sufficiently flexible that one can expect that a member of this family could describe the prior opinion of the data analyst about each of the interested parameters. In addition, it is always a potential choice to describe various positive populations from exponential to normal shape (see, for example, [26]). Under the above prior distributions setting, the prior distribution of parameter τi−δ appears as distribution of difference of two nonindependent gamma random variables. While deriving the distribution of sum of two gamma random variables is an elementary work, this is not a trivial work for the case of difference of random variables especially when they are nonindependent. Generally, there is a vast literature on the sum of independent and nonindependent gamma random variables, and many authors [27–29] discussed this problem under various assumptions such as independent and identically or nonidentically distributed random variables. However, the difference of the nonindependent gamma random variables has been less considered in the literature. Holm and Alouini [30] showed that the general formulas for the probability density function of the difference of two not necessarily identically distributed gamma random variables can be expressed in the form of the McKay Bessel function distribution of type II [31]. The random variable V follows the type II MacKay distribution with parameters a>−1/2, b>0, and c<1, denoted by McKayII (a, b, and c), when the density function of V is given by
(15)fVv=1−c2a+1/2vaπ2aba+1Γa+1/2e−vc/bKavb,v≠0,1−c2a+1/2π2bΓaΓa+1/2,v=0,where Ka· is the modified Bessel function of the second kind and of order a. He provided the moments of this distribution in terms of the Gauss hypergeometric function. In addition, there are relatively simple expressions for the moments generating function, cumulants, and the moments of the McKay distributions. It should be noted that several famous distributions such as normal; distributions of Pearson system, e.g., type III Pearson (see [32]); and distribution of first product moment coefficient in a sample from normal population are special cases of McKay distribution of type II. The support of the McKay distribution is the line of real numbers, and its density function has a skew-symmetric shape. This implies that the prior distributions given in Equation (14) are in consistent with our knowledge about the parameters of the model. Specifically, the prior distribution of parameter β holds the necessary constrain μxi>0,i=1,⋯,n.


### 4.2. Posterior Inference

Considering the likelihood function (9) and the prior distributions (14), the joint posterior distribution of the model parameters in terms of (μ,λ,δ) is obtained to be
(16)πμ,λ,δ∣y∝Lμ,λ,δ∣yπμ,λ,δ=Lμ,λ,δ∣yπμ∣λπλπδ,∝∏i=1nδ+yiδ+μi−1λ1/2fIGyi,μi,λ×∏i=1nλ2δ+μi×exp−λδ+μi×λc0−1exp−λd0×δe0−1exp−f0δ∝∏i=1nδ+yiλ1/2fIGyi∣μI,λ×λ2nexp−λ∑i=1nδ+μI×λc0−1exp−λd0×δe0−1exp−f0δ.Due to the equalityμi−1=xi′β, the posterior distribution (16) can be written in terms of interested parameters β,λ,δ as
(17)πβ,λ,δ∣y∝λ5/2n+c0−1exp−λQβ2+∑i=1n1xi′β+nδ+d0×δe0−1exp−f0δ∏i=1nδ+yi,where Qβ is a quadratic form in terms of the vector of regression coefficients, β, given in Equation (11). Clearly, the complexity of joint posterior distribution (17) precludes deriving analytical posterior inferences. In what follows, we employed the Markov chain Monte Carlo (MCMC) methods to sample form posterior distribution in order to deriving numerical estimates for interested posterior quantities.


### 4.3. Gibbs Sampler

To construct a Gibbs sampler, it is necessary to determine the full conditional posterior distributions. The kernel of full conditional distribution for a typical parameter is easily obtained by vanishing all unrelated quantities to this parameter from the joint posterior distributions. According to (17), the full conditional posterior distribution of parameter β in proportional form is given by
(18)πβ∣δ,λ∝exp−λQβ2+∑i=1n1xi′β,which cannot be written down in a closed form density function. Also, the full conditional distribution of parameter λ is
(19)πλ∣β,δ,y∝λ5/2n+c0−1exp−λQβ2+∑i=1n1xi′β+nδ+d0.Therefore, the full conditional distribution for this parameter, in closed form, is obtained to be
(20)λ∣β,δ∼Gamma52n+c0,Qβ2+∑i=1n1i′β+nδ+d0.For parameterδ, one can write
(21)πδ∣β,λ,y∝∏i=1nδ+yiδe0−1exp−δnλ+f0=∏i=1nδ1/ne0−1δ+yiexp−δnλ+f0=∏i=1nδ1/ne0−1+1exp−δnλ+f0+∏i=1nyiδ1/ne0−1exp−δnλ+f0=δe0+n−1exp−δnλ+f0+∏i=1nyi∏i=1nδ1/ne0−1exp−δnλ+f0∝Γe0+nnλ+f0e0+nfGδ∣e0+n,nλ+f0+∏i=1nyiΓe0−1nλ+f0e0−1×fGδ∣e0−1,nλ+f0,where fG·∣α,β denotes the density function of a gamma distribution with mean αβ. As it can be seen, the full conditional posterior distribution of parameter δ is proportional to a linear combination of two gamma density functions with equal shape parameters with positive coefficients. It is easy to show that the full conditional posterior distribution of δ is a two-component gamma mixture distribution given by
(22)πδ∣β,λ,η=ηfGδ∣e0+n,nλ+f0+1−ηfGδ∣e0−1,nλ+f0,where η=c1/c1+c2 with
(23)c1=Γe0+nΓe0+n+∏i=1nyiΓe0−1nλ+f0n+1,c2=∏i=1nyiΓe0−1Γe0+n/nλ+f0n+1+∏i=1nyiΓe0−1.Knowing the full conditional distributions of the model parameters, it is possible to sample from the joint posterior distribution of the model parameter. Since there is no closed form expression for the full conditional distribution of parameterβ, a nested Metropolis-Hastings algorithm should be used for sampling from this parameter within the iterations of Gibbs algorithm.Also, consideringy=y1,⋯,yn as old observations, the posterior predictive distribution for a new observation yn+1 which corresponds to vector of explanatory variablesyn+1 is given by
(24)pyn+1∣y=∫∫∫fWIGyn+1∣y,β,λ,δπβ,λ,δ∣ydβdλdδyn+1>0,where fWIGβ,λ,δ and πβ,λ,δ∣y are the density function of WIG distribution and joint posterior distribution of the model parameters given in Equations (1) and (17), respectively. The MCMC estimate of the posterior predictive distribution is obtained to be
(25)p^yn+1∣y=1S∑s=1SfWIGyn+1∣βs,λs,δsyn+1>0,where βs,λs,δss=1S denotes a MCMC sample from the joint posterior distribution of parameters.


## 4.1. Prior Setting

In regression analysis, the main interested parameter that should be estimated is the vector of regression coefficient,β=β1,⋯,βp. Theoretically, any multivariate continuous distribution on ℝp could be employed as a prior distribution for this parameter. For example, in a vague or a noninformative Bayesian analysis, a flat multivariate normal or a multivariate uniform distribution with large variance is a usual choice for this purpose. But, in the proposed model, considering such prior distributions for β may lead to the nonpositive values for μi=1/xi′β, i=1,⋯,n, which is inconsistent with the positivity constrain for the location parameter of WIG distribution. Taking these considerations into account, we consider the following joint prior distribution of the model parameters as
(14)τi∣λ∼Gamma2,λ,λ∼Gammac0,d0,δ∼Gammae0,f0,where τi=μi+δ and zero indices indicate the hyperparameters of the model that should be determined based on prior beliefs. As it can be seen, in this prior setting, all of the prior knowledge about the model parameters is expressed by using three gamma distributions. The family of gamma distributions seems to be sufficiently flexible that one can expect that a member of this family could describe the prior opinion of the data analyst about each of the interested parameters. In addition, it is always a potential choice to describe various positive populations from exponential to normal shape (see, for example, [26]). Under the above prior distributions setting, the prior distribution of parameter τi−δ appears as distribution of difference of two nonindependent gamma random variables. While deriving the distribution of sum of two gamma random variables is an elementary work, this is not a trivial work for the case of difference of random variables especially when they are nonindependent. Generally, there is a vast literature on the sum of independent and nonindependent gamma random variables, and many authors [27–29] discussed this problem under various assumptions such as independent and identically or nonidentically distributed random variables. However, the difference of the nonindependent gamma random variables has been less considered in the literature. Holm and Alouini [30] showed that the general formulas for the probability density function of the difference of two not necessarily identically distributed gamma random variables can be expressed in the form of the McKay Bessel function distribution of type II [31]. The random variable V follows the type II MacKay distribution with parameters a>−1/2, b>0, and c<1, denoted by McKayII (a, b, and c), when the density function of V is given by
(15)fVv=1−c2a+1/2vaπ2aba+1Γa+1/2e−vc/bKavb,v≠0,1−c2a+1/2π2bΓaΓa+1/2,v=0,where Ka· is the modified Bessel function of the second kind and of order a. He provided the moments of this distribution in terms of the Gauss hypergeometric function. In addition, there are relatively simple expressions for the moments generating function, cumulants, and the moments of the McKay distributions. It should be noted that several famous distributions such as normal; distributions of Pearson system, e.g., type III Pearson (see [32]); and distribution of first product moment coefficient in a sample from normal population are special cases of McKay distribution of type II. The support of the McKay distribution is the line of real numbers, and its density function has a skew-symmetric shape. This implies that the prior distributions given in Equation (14) are in consistent with our knowledge about the parameters of the model. Specifically, the prior distribution of parameter β holds the necessary constrain μxi>0,i=1,⋯,n.


## 4.2. Posterior Inference

Considering the likelihood function (9) and the prior distributions (14), the joint posterior distribution of the model parameters in terms of (μ,λ,δ) is obtained to be
(16)πμ,λ,δ∣y∝Lμ,λ,δ∣yπμ,λ,δ=Lμ,λ,δ∣yπμ∣λπλπδ,∝∏i=1nδ+yiδ+μi−1λ1/2fIGyi,μi,λ×∏i=1nλ2δ+μi×exp−λδ+μi×λc0−1exp−λd0×δe0−1exp−f0δ∝∏i=1nδ+yiλ1/2fIGyi∣μI,λ×λ2nexp−λ∑i=1nδ+μI×λc0−1exp−λd0×δe0−1exp−f0δ.Due to the equalityμi−1=xi′β, the posterior distribution (16) can be written in terms of interested parameters β,λ,δ as
(17)πβ,λ,δ∣y∝λ5/2n+c0−1exp−λQβ2+∑i=1n1xi′β+nδ+d0×δe0−1exp−f0δ∏i=1nδ+yi,where Qβ is a quadratic form in terms of the vector of regression coefficients, β, given in Equation (11). Clearly, the complexity of joint posterior distribution (17) precludes deriving analytical posterior inferences. In what follows, we employed the Markov chain Monte Carlo (MCMC) methods to sample form posterior distribution in order to deriving numerical estimates for interested posterior quantities.


## 4.3. Gibbs Sampler

To construct a Gibbs sampler, it is necessary to determine the full conditional posterior distributions. The kernel of full conditional distribution for a typical parameter is easily obtained by vanishing all unrelated quantities to this parameter from the joint posterior distributions. According to (17), the full conditional posterior distribution of parameter β in proportional form is given by
(18)πβ∣δ,λ∝exp−λQβ2+∑i=1n1xi′β,which cannot be written down in a closed form density function. Also, the full conditional distribution of parameter λ is
(19)πλ∣β,δ,y∝λ5/2n+c0−1exp−λQβ2+∑i=1n1xi′β+nδ+d0.Therefore, the full conditional distribution for this parameter, in closed form, is obtained to be
(20)λ∣β,δ∼Gamma52n+c0,Qβ2+∑i=1n1i′β+nδ+d0.For parameterδ, one can write
(21)πδ∣β,λ,y∝∏i=1nδ+yiδe0−1exp−δnλ+f0=∏i=1nδ1/ne0−1δ+yiexp−δnλ+f0=∏i=1nδ1/ne0−1+1exp−δnλ+f0+∏i=1nyiδ1/ne0−1exp−δnλ+f0=δe0+n−1exp−δnλ+f0+∏i=1nyi∏i=1nδ1/ne0−1exp−δnλ+f0∝Γe0+nnλ+f0e0+nfGδ∣e0+n,nλ+f0+∏i=1nyiΓe0−1nλ+f0e0−1×fGδ∣e0−1,nλ+f0,where fG·∣α,β denotes the density function of a gamma distribution with mean αβ. As it can be seen, the full conditional posterior distribution of parameter δ is proportional to a linear combination of two gamma density functions with equal shape parameters with positive coefficients. It is easy to show that the full conditional posterior distribution of δ is a two-component gamma mixture distribution given by
(22)πδ∣β,λ,η=ηfGδ∣e0+n,nλ+f0+1−ηfGδ∣e0−1,nλ+f0,where η=c1/c1+c2 with
(23)c1=Γe0+nΓe0+n+∏i=1nyiΓe0−1nλ+f0n+1,c2=∏i=1nyiΓe0−1Γe0+n/nλ+f0n+1+∏i=1nyiΓe0−1.Knowing the full conditional distributions of the model parameters, it is possible to sample from the joint posterior distribution of the model parameter. Since there is no closed form expression for the full conditional distribution of parameterβ, a nested Metropolis-Hastings algorithm should be used for sampling from this parameter within the iterations of Gibbs algorithm.Also, consideringy=y1,⋯,yn as old observations, the posterior predictive distribution for a new observation yn+1 which corresponds to vector of explanatory variablesyn+1 is given by
(24)pyn+1∣y=∫∫∫fWIGyn+1∣y,β,λ,δπβ,λ,δ∣ydβdλdδyn+1>0,where fWIGβ,λ,δ and πβ,λ,δ∣y are the density function of WIG distribution and joint posterior distribution of the model parameters given in Equations (1) and (17), respectively. The MCMC estimate of the posterior predictive distribution is obtained to be
(25)p^yn+1∣y=1S∑s=1SfWIGyn+1∣βs,λs,δsyn+1>0,where βs,λs,δss=1S denotes a MCMC sample from the joint posterior distribution of parameters.


## 5. Data Illustration

In this section, we preformed a simulation study and analyzed a real data set to access the performance of the proposed model and to explain its applicability. To provide a situation for an impartial judgment, we considered the IG and log-normal regression models as two famous competitors for the proposed model in computations.(i)
Bayesian IG Regression Model. The IG regression model has been studied repeatedly in the literature and can be written as yi∣xi∼IGμxi,λ;i=1,⋯,n, with the same link function as the proposed WIG regression model, i.e., μxi−1=xi′β. For the Bayesian analysis of this model, we used a vague structure for the prior distributions by considering β∣λ∼Npu0,λ−1V0 and λ∼Gammaa0/2,b0/2, where hyperparameters were set to be a0=b0=0.1,u0=0,⋯,0′∈Rp, and V0=k0Ip with k0=102. Under these prior settings, it can be shown that the posterior distributions of the parameters λ and βare members of gamma and noncentral t-student family of distributions, and their corresponding Bayes estimators are given by λ^=a0+na0+n/b0+∑i=1nyi−1−n1n′Xβb0+∑i=1nyi−1−n1n′Xβ and β=nk0−1Ip+X′YX−1X′1n, respectively. We refer the readers to [6, 22] for more details about the modeling by IG distribution(ii)
Bayesian Log-Normal Regression Model. As it can be found in any Bayesian textbook (see, for example, [33]), to construct a Bayesian regression model for a log-normally distributed response variable, Y, it is suffice to construct a normal regression model for the logarithm of the response observations as z∣X,β,σ2∼NnXβ,σ2I, where z=z1,⋯,zn′ with zi=logyi;i=1,⋯,n. Using a semiconjugate structure for the prior distributions of the formβ∼Npβ0,τ0Ip, 1/σ2∼Gammaν0/2,ν0/2σ02, it is straightforward to contract a Gibbs sampler in order to approximate the joint posterior distribution and interested posterior quantities by employing the following full conditional distributions:(26)β∣σ2,z∼Npτ0−1σ2Ip+X′X−1X′y,τ0−1σ2Ip+X′X−1,1σ2∣β,z∼Gammaν0+n2,12ν0σ02+z−Xβ′z−Xβ.It should be noted that the hyperparametersτ0, ν0, and σ02 have been set, in our studies, in such manner that produce a vague prior setting.


### 5.1. Simulation Study

We considered a WIG regression model with two explanatory variables given by(27)yi∣xi∼WIGμxi,λ,δ,μxi=1xi1β1+xi2β2,with xi=xi1,xi2, λ=2, and δ=3. The values of explanatory variables simulated from independent standard normal distributions and the values of regression coefficients are set to be β=β1,β2′=0.1,0.2′. Then, the response observations are simulated from model (27) for different sample sizes 25, 50, 75, and 100. Given simulated data, we used the proposed Gibbs sampler to create a MCMC sample in order to estimate the posterior interested parameters after 5000 burn-in period. The Geweke [34] diagnostic was used to test convergence of the algorithm. The root mean square error (RMSE=1/R∑r=1Rθ∧r−θ2 and absolute value of the relative bias RB=1/R∑r=1Rθ∧r/θ−1 of Bayes estimators of the regression coefficients for the WIG model along with corresponding values for IG and log-normal model were calculated, where θ is the real value of the parameter, θ∧r denotes the value of the parameter θ in r-th repetition, and θ¯ is the estimated posterior mean. In all computations, the number of repetitions, R, is fixed to be 4000 in order to taking into account the uncertainty of the simulated data. The results are presented in Table 1. In Table 2, a general comparison of strengths and limitations of different discussed models has been presented.Table 1
The values of RMSE and RB of the Bayes estimators of the regression coefficients for WIG model along with corresponding values for IG and LN models.
Sample sizeParameterCriterionRegression modelWIGIGLN25β1RB0.23200.29801.9473RMSE1.61031.81961.9476β2RB1.16861.19222.9271RMSE1.51671.59732.927350β1RB0.78691.02801.5233RMSE0.82081.03301.5320β2RB1.80002.25012.5713RMSE1.84022.33912.578175β1RB0.80161.11441.7650RMSE0.84621.17611.7662β2RB1.63321.89602.8198RMSE1.68421.96942.8013100β1RB0.62180.96091.8919RMSE0.74701.09041.8726β2RB1.50201.67892.8914RMSE1.60141.81232.8917Table 2
A general comparison of strengths and limitations of different discussed models.
Regression modelStrengthsLimitationsLN (see, for example, [32])(i) The LN regression model is the most famous traditional candidate for analyzing positive right-skewed response observations(ii) It can be fitted and interpreted easily based on a normal model by using a simple log transformation(i) Moments of the response variable depend overly on exact log-normality assumption(ii) Small sample sequential moments of the response variable oscillate to excess as the sample size increases(iii) The distribution of response variable has too thick right-hand tail to be plausible(iv) LN distribution is also almost symmetrical, and there is no way to control the skewness or asymmetry of the distributionIG (see, for example, [18, 22, 23, 35])(i) The sampling theory of the response observations is tractable(ii) There has been a growing attention to the use of IG distribution, and recent researches have revealed rigid evidence supporting the IG distribution in comparison with the log-normal in most applicationsFitting procedure under IG distribution is more complicated than corresponding procedure for LN modelWIG (see, for example, [18, 22, 23, 35])(i) The sampling theory of the response observations is tractable(ii) There has been a growing attention to the use of IG distribution, and recent researches have revealed rigid evidence supporting the IG distribution in comparison with the log-normal in most applicationsFitting procedure for WIG model is more complicated than corresponding procedure in IG modelAs it can be seen, the absolute value of relative biases and the root mean square errors of the proposed WIG regression model are less than the corresponding values for both IG and log-normal regression models. This implies that the Bayes estimator of regression coefficients for WIG model is more efficient than the counterpart estimators in IG and log-normal models.


### 5.2. Application to the Seismological Data

To investigate the effect of soil shear wave velocity (VS30) in thirty meters at top of the site and the Joyner-Boore distance (RJB) as explanatory variables on the peak ground acceleration (PGA) as response variable, 25 different sites are selected and interested variables are measured. The data, presented in Table 3, are a subset of the Imperial Valley earthquake in 1979. The histogram and box plot of response observations, given in Figure 1, show the right-skewed structure of data. This indicates that the WIG and some of other right-skewed distributions such as IG and log-normal could be suitable candidate distributions for fitting to response observations.Table 3
A subset of the Imperial Valley earthquake data in 1979.
SitePGAVS30RJBSitePGAVS30RJB10.290.860.44140.260.910.9820.690.030.54150.380.630.7930.180.710.45160.410.220.9040.240.780.17170.450.510.7650.100.530.29180.440.080.5660.270.490.52190.530.260.8670.150.570.03200.090.890.0980.220.330.31210.170.690.3590.310.850.07220.050.640.69100.200.250.76230.160.380.33110.370.880.60240.130.050.61120.120.940.94250.140.920.56130.080.670.75Figure 1
The histogram and box plot of peak ground acceleration data.
(a)(b)The result of the Kolmogorov-Smirnov goodness of fit test for the WIG distribution along with corresponding values for the IG distributions is given in Table4. We also provided the corresponding values for the log-normal distribution as a traditional candidate distribution for modeling positive right-skewed data. As it can be seen, the null hypotheses of the WIG and IG distributions for response distributions both are not rejected at 0.05 level of significance. The Bayesian estimates of the model parameters along with their corresponding 95% HPD credible intervals for different models are presented in Table 5. It should be pointed out that the estimates also are obtained as the posterior mean of the MCMC samples after a 5000 burn-in period. The value of hyperparameters in (14) is set to be c0=d0=e0=f0=0.01 in order to provide noninformative vague prior distributions with large variances for parameters of the model. The Geweke diagnostic test [34] was used to assess the convergence of the sampler.Table 4
The values of the Kolmogorov-Smirnov goodness of fit test for WIG, IG, and LN distributions.
DistributionTest statisticP valueWIG0.14270.6366IG0.24620.0805LN0.54220.0024Table 5
The Bayes estimates and 95% HPD credible intervals of regression coefficients for WIG model along with corresponding values for IG and LN models.
Regression modelβ1β2λδσWIG5.01 (1.30, 6.06)1.50 (1.25, 3.40)1.16 (0.57, 1.56)0.88 (1.82, 0.51)—IG5.12 (2.24, 8.00)1.80 (-0.82,4.43)0.51 (0.25, 1.81)——LN5.81 (2.52, 6.01)1.78 (1.56, 4.01)——0.6184 (0.34, 1.11)As it can be seen, from the95% HPD credible intervals, for the WIG and LN models, both of the regression coefficients β1 and β2 are different from zero at a significance level of 0.05, and this is the same for other parameters of these models. But, for the model IG, the coefficient β2 is not significantly different from zero at the same level. It should be noticed that the estimated values of the regression coefficients for LN model are not comparable to corresponding values in the WIG and IG models due to the differences between the functional forms of the link functions of these models. In such situations, one can use some information-based criteria to assess the fit and complexity of models or provide a criterion for evaluating and comparing the predictive power of them. In this direction, the Akaike information criterion (AIC; −2ℓθ+2k), Bayesian information criterion (BIC; −ℓθ+klogn), and Hannan-Quinn information criterion (HQIC; −2klθ+2kloglogn) are employed, where in definition of these criteria, ℓ·, θ, n, and k denote the log-likelihood function of the model, vector of estimated parameters, sample size, and the number of model parameters, respectively. The lesser values of these criteria, presented in Table 6, indicate better parsimonious fits. We also used the cross-validation mean square error of prediction criterion,
(28)CVMSEP=1n∑i=1nyi−y∧−i2,to calculate the predictive performance of different models. Where in Equation (28), y^−i denotes the predicated value of yi based on the data yj,xj,j≠i=1,⋯,n, i.e., all observations expect i-th observation.Table 6
The estimated regression models along with the corresponding values of AIC, BIC, HQIC, andCVMSEP criteria for different models.
ModelEstimated modelAICBICHIQCCVMSEPWIGμ∧i−1=5.01x1+1.50x2193.94196.38194.620.0322IGμ∧i−1=5.12x1+1.89x2229.62232.06230.290.0372LNlogμ^i=5.81x1+1.87x2260.63270.07276.522.0781The values of AIC, BIC, and HQIC show that the data provide more support for the WIG regression model than other models considering the parsimonious principle. Also, the values ofCVMSEP indicate to the higher prediction accuracy, in average, for the WIG model.


## 5.1. Simulation Study

We considered a WIG regression model with two explanatory variables given by(27)yi∣xi∼WIGμxi,λ,δ,μxi=1xi1β1+xi2β2,with xi=xi1,xi2, λ=2, and δ=3. The values of explanatory variables simulated from independent standard normal distributions and the values of regression coefficients are set to be β=β1,β2′=0.1,0.2′. Then, the response observations are simulated from model (27) for different sample sizes 25, 50, 75, and 100. Given simulated data, we used the proposed Gibbs sampler to create a MCMC sample in order to estimate the posterior interested parameters after 5000 burn-in period. The Geweke [34] diagnostic was used to test convergence of the algorithm. The root mean square error (RMSE=1/R∑r=1Rθ∧r−θ2 and absolute value of the relative bias RB=1/R∑r=1Rθ∧r/θ−1 of Bayes estimators of the regression coefficients for the WIG model along with corresponding values for IG and log-normal model were calculated, where θ is the real value of the parameter, θ∧r denotes the value of the parameter θ in r-th repetition, and θ¯ is the estimated posterior mean. In all computations, the number of repetitions, R, is fixed to be 4000 in order to taking into account the uncertainty of the simulated data. The results are presented in Table 1. In Table 2, a general comparison of strengths and limitations of different discussed models has been presented.Table 1
The values of RMSE and RB of the Bayes estimators of the regression coefficients for WIG model along with corresponding values for IG and LN models.
Sample sizeParameterCriterionRegression modelWIGIGLN25β1RB0.23200.29801.9473RMSE1.61031.81961.9476β2RB1.16861.19222.9271RMSE1.51671.59732.927350β1RB0.78691.02801.5233RMSE0.82081.03301.5320β2RB1.80002.25012.5713RMSE1.84022.33912.578175β1RB0.80161.11441.7650RMSE0.84621.17611.7662β2RB1.63321.89602.8198RMSE1.68421.96942.8013100β1RB0.62180.96091.8919RMSE0.74701.09041.8726β2RB1.50201.67892.8914RMSE1.60141.81232.8917Table 2
A general comparison of strengths and limitations of different discussed models.
Regression modelStrengthsLimitationsLN (see, for example, [32])(i) The LN regression model is the most famous traditional candidate for analyzing positive right-skewed response observations(ii) It can be fitted and interpreted easily based on a normal model by using a simple log transformation(i) Moments of the response variable depend overly on exact log-normality assumption(ii) Small sample sequential moments of the response variable oscillate to excess as the sample size increases(iii) The distribution of response variable has too thick right-hand tail to be plausible(iv) LN distribution is also almost symmetrical, and there is no way to control the skewness or asymmetry of the distributionIG (see, for example, [18, 22, 23, 35])(i) The sampling theory of the response observations is tractable(ii) There has been a growing attention to the use of IG distribution, and recent researches have revealed rigid evidence supporting the IG distribution in comparison with the log-normal in most applicationsFitting procedure under IG distribution is more complicated than corresponding procedure for LN modelWIG (see, for example, [18, 22, 23, 35])(i) The sampling theory of the response observations is tractable(ii) There has been a growing attention to the use of IG distribution, and recent researches have revealed rigid evidence supporting the IG distribution in comparison with the log-normal in most applicationsFitting procedure for WIG model is more complicated than corresponding procedure in IG modelAs it can be seen, the absolute value of relative biases and the root mean square errors of the proposed WIG regression model are less than the corresponding values for both IG and log-normal regression models. This implies that the Bayes estimator of regression coefficients for WIG model is more efficient than the counterpart estimators in IG and log-normal models.


## 5.2. Application to the Seismological Data

To investigate the effect of soil shear wave velocity (VS30) in thirty meters at top of the site and the Joyner-Boore distance (RJB) as explanatory variables on the peak ground acceleration (PGA) as response variable, 25 different sites are selected and interested variables are measured. The data, presented in Table 3, are a subset of the Imperial Valley earthquake in 1979. The histogram and box plot of response observations, given in Figure 1, show the right-skewed structure of data. This indicates that the WIG and some of other right-skewed distributions such as IG and log-normal could be suitable candidate distributions for fitting to response observations.Table 3
A subset of the Imperial Valley earthquake data in 1979.
SitePGAVS30RJBSitePGAVS30RJB10.290.860.44140.260.910.9820.690.030.54150.380.630.7930.180.710.45160.410.220.9040.240.780.17170.450.510.7650.100.530.29180.440.080.5660.270.490.52190.530.260.8670.150.570.03200.090.890.0980.220.330.31210.170.690.3590.310.850.07220.050.640.69100.200.250.76230.160.380.33110.370.880.60240.130.050.61120.120.940.94250.140.920.56130.080.670.75Figure 1
The histogram and box plot of peak ground acceleration data.
(a)(b)The result of the Kolmogorov-Smirnov goodness of fit test for the WIG distribution along with corresponding values for the IG distributions is given in Table4. We also provided the corresponding values for the log-normal distribution as a traditional candidate distribution for modeling positive right-skewed data. As it can be seen, the null hypotheses of the WIG and IG distributions for response distributions both are not rejected at 0.05 level of significance. The Bayesian estimates of the model parameters along with their corresponding 95% HPD credible intervals for different models are presented in Table 5. It should be pointed out that the estimates also are obtained as the posterior mean of the MCMC samples after a 5000 burn-in period. The value of hyperparameters in (14) is set to be c0=d0=e0=f0=0.01 in order to provide noninformative vague prior distributions with large variances for parameters of the model. The Geweke diagnostic test [34] was used to assess the convergence of the sampler.Table 4
The values of the Kolmogorov-Smirnov goodness of fit test for WIG, IG, and LN distributions.
DistributionTest statisticP valueWIG0.14270.6366IG0.24620.0805LN0.54220.0024Table 5
The Bayes estimates and 95% HPD credible intervals of regression coefficients for WIG model along with corresponding values for IG and LN models.
Regression modelβ1β2λδσWIG5.01 (1.30, 6.06)1.50 (1.25, 3.40)1.16 (0.57, 1.56)0.88 (1.82, 0.51)—IG5.12 (2.24, 8.00)1.80 (-0.82,4.43)0.51 (0.25, 1.81)——LN5.81 (2.52, 6.01)1.78 (1.56, 4.01)——0.6184 (0.34, 1.11)As it can be seen, from the95% HPD credible intervals, for the WIG and LN models, both of the regression coefficients β1 and β2 are different from zero at a significance level of 0.05, and this is the same for other parameters of these models. But, for the model IG, the coefficient β2 is not significantly different from zero at the same level. It should be noticed that the estimated values of the regression coefficients for LN model are not comparable to corresponding values in the WIG and IG models due to the differences between the functional forms of the link functions of these models. In such situations, one can use some information-based criteria to assess the fit and complexity of models or provide a criterion for evaluating and comparing the predictive power of them. In this direction, the Akaike information criterion (AIC; −2ℓθ+2k), Bayesian information criterion (BIC; −ℓθ+klogn), and Hannan-Quinn information criterion (HQIC; −2klθ+2kloglogn) are employed, where in definition of these criteria, ℓ·, θ, n, and k denote the log-likelihood function of the model, vector of estimated parameters, sample size, and the number of model parameters, respectively. The lesser values of these criteria, presented in Table 6, indicate better parsimonious fits. We also used the cross-validation mean square error of prediction criterion,
(28)CVMSEP=1n∑i=1nyi−y∧−i2,to calculate the predictive performance of different models. Where in Equation (28), y^−i denotes the predicated value of yi based on the data yj,xj,j≠i=1,⋯,n, i.e., all observations expect i-th observation.Table 6
The estimated regression models along with the corresponding values of AIC, BIC, HQIC, andCVMSEP criteria for different models.
ModelEstimated modelAICBICHIQCCVMSEPWIGμ∧i−1=5.01x1+1.50x2193.94196.38194.620.0322IGμ∧i−1=5.12x1+1.89x2229.62232.06230.290.0372LNlogμ^i=5.81x1+1.87x2260.63270.07276.522.0781The values of AIC, BIC, and HQIC show that the data provide more support for the WIG regression model than other models considering the parsimonious principle. Also, the values ofCVMSEP indicate to the higher prediction accuracy, in average, for the WIG model.


## 6. Conclusions

Modeling right-skewed positive data has been at the core of interest in lifetime analysis and many other fields. This paper developed a regression model by considering the weighted inverse Gaussian distribution for response observations as a versatile candidate for analysis of the right-skewed positive data. A Bayesian framework employed to analyze the model. The proposed method uses a set of noninformative flexible priors from the gamma family of distributions and provides a tractable joint posterior distribution for parameters. It is computationally straightforward and proceeds by constructing a Gibbs sampler to derive the posterior inferences. The comparative experimental studies, simulation, and real-world example showed that the proposed weighted inverse Gaussian regression model has the ability to dominate inverse Gaussian and log-normal as two most famous competitor models in the literature. Particularly, considering that the relative prediction efficiency (REP) of two given predictorsz^1 and z^2 in terms of cross-validation mean square error of prediction criterion is defined as REPz^1,z^2=CVMSEPz^2/CVMSEPz^1, it can be seen that the relative prediction efficiency of the proposed WIG regression model to the IG and LN regression models is given by 1.16 and 64, respectively, for the real-world example discussed in this paper.



---

*Source: 3943930-2022-11-17.xml*