Suppose that y=(y1,y2,…,yn) ′ is a response variable for subjects i = 1, 2, …, n, and that x=(x1,x2,…,xp) ′ is covariate. For 0 < τ < 1, let Q_yi(τ|x_i) denote the τ-th quantile regression function of y_i with associated p dimensional vector of covariates x_i. The quantile regression function is expressed in the form of Qyi(τ|xi)=xiTβ, for i = 1, 2, …, n, where β is a p × 1 vector of coefficients for indicator variables that depend on τ. Then, a linear quantile regression model can be expressed as

Here, e_i is the error term whose distribution is restricted so much that τ-th quantile is equal to zero. Then, quantile regression estimation for β is obtained by minimizing

where ρ_τ(u) is the check function defined by

Here, I(.) is an indicator function that takes the value of unity when I(.) is true and zero otherwise and here u=yi-xiTβ. However, this indicator function is not differentiable at zero, and explicit solutions to minimization problems are unobtainable [26, 27]. In quantile regression methods, linear programming is often implemented for parameter estimation.

Yu and Moyeed [27] found that minimizing the expression (2) is equivalent to maximizing a likelihood function formed by combining the independently distributed asymmetric Laplace error distribution. The asymmetric Laplace distribution is employed in likelihood distribution in order to make Bayesian estimation more natural and effective [27, 28], since this distribution is a possible parametric link between the minimization problem in Equation (2) and maximum likelihood theory. Therefore, a random variable ε_i is said to be distributed as a skewed Laplace distribution with density [29].

where δ is a scale parameter. It is also known that the mean and variance of the asymmetric Laplace distribution are given respectively by

One property of the asymmetric Laplace distribution is this distribution can be represented with various mixture representations. The Gibbs sampling algorithm is then utilized for Bayesian analysis of the quantile regression model based on a skewed Laplace distribution's theoretical derivation. The mixture of the exponential and normal distribution of the skewed Laplace distribution allows efficient Gibbs sampling [29–31]. More specifically let

where

Equation (1) and (6) lead to

The random variables z_i and ξ_i are mutually independent. Variable z_i is exponentially distributed with a mean of 1/δ. Variable ξ_i is a standard normal distribution. Thus, the conditional distribution of y_i given z_i is normally distributed with a mean xiTβ+θzi and variance ϕ2zi. Then the quantile regression model here is represented as a normal regression model. This representation provides an easy way to construct a Gibbs sampler and save time in sampling the regression coefficients.

One crucial problem in building a quantile regression model is the selection of predictors. The prediction accuracy is often improved by choosing an appropriate subset of predictors. It is often desired to identify a smaller subset of predictors from a large set of predictors to obtain better interpretation in practice. There have been active studies on the sparse representation of linear regression. Li et al. [19] showed that the least absolute shrinkage and selection operator (Lasso) technique could simultaneously perform variable selection and parameter estimation. The Lasso estimate, which is also known as an L1-regularized least squares estimate, involves the sum of the coefficients' absolute values as the penalty. This L1-regularized has the advantage of simultaneously controlling the fitted coefficients' variance and performing the automatic variable selection. The Lasso estimates are defined as Alhamzawi and Yu [9, 20] and Zou [32].

where λ > 0 is a regularized parameter that controls the degree of penalization. The second term in the expression (9) is an L₁-regularized term, which could be interpreted as a Bayesian posterior mode estimate under independent Laplace priors for the coefficients. As a nonnegative regularization parameter λ increases, the Lasso estimates continuously the shrinks quantile regression coefficient toward zero. The appropriate prior distribution for β_k ∈ β(k = 1, ..., p) is Laplace distribution, defined as follows

Here, it is assumed that the residuals ε_i have a skewed Laplace distribution as represented in Equation (4). Then, we employed the Bayesian Adaptive LASSO quantile regression (BALQR) to estimate the unknown parameter model. The penalty function in (9) can be made “adaptive” by choosing different shrinkages for different coefficients:

where λ_k > 0 is the tuning parameter for the k^th coefficient. Here, the new proposed Laplace prior for β_k is formed by

This equation can be represented as a scale mixture of normal with an exponential mixing density [29, 33, 34].

Let ak=δ1/2λk. The prior distribution for β_k is expressed in the form

or in the form

Then, we consider the class of inverse gamma priors on λk2 of the form

By combining Equations (13) and (14), we obtain the posterior density function of λk2 is inverse gamma with shape parameter 1 + σ and rate parameter δsk2+ρ. The value of two hyperparameters σ and ρ affect the amount of shrinkage in the prior Equation (14), e.g., larger σ and smaller ρ lead to bigger penalization. This BALQR uses a Laplace prior for β_k such that each β_k has a Lasso type of penalization parameter δ1/2λk. More detailed explanations in formulating posterior density function for all parameter models are written in Alhamzawi et al. [29] and Xu and Tang [34].

To estimate credible intervals, it is not automatically valid by constructing the posterior. Yang et al. [35] argued to employ the Wald method based on the asymptotic approximation to the variance-covariance matrix of the posterior sequences to estimate the Bayes credible interval, as also reported in Li et al. [19] and Yue and Hang [36]. In this present study, we implemented the Wald method based on the asymptotic approximation to the variance-covariance matrix of the posterior sequences to estimate the Bayes credible interval.

In this section, we demonstrate the application of the Bayesian, Bayesian Lasso, and Bayesian Adaptive Lasso in quantile regression to several different generating processes. The goal of this simulation study here is to reveal the performance of the proposed methods and their associated algorithm in recovering the true parameters.

For the proposed model, the MCMC simulations were implemented in R version 3.6.1 [37]. In this simulation study, the response variable y_i is generated from the following regression model.

where each covariate x_ik is simulated from a standard normal distribution and β = (5, 1.5, 3). Three different distributions for e_i were considered: (i) heteroscedastic normal, (1 + x₁)N(0, 1), (ii) autocorrelated error, sin(seq(0.1π, 18.3π, 0.1π) + Z) with Z (0, 0.1), (iii) nonnormal error, the mixture of two Student's t distribution, 0.1t₍₁₎ + 0.9t₍₃₎. For each choice of the error distribution, we employed a Bootstrap resampling method with 100 simulations were carried out. In each simulation, 150 observations were generated. Four different values of the given quantile τ = 0.25, 0.50, 0.75, and 0.95 were considered. To assess the sampling efficiency of the proposed algorithm, the Monte Carlo standard errors for each β_k, k = 1, 2, 3 were calculated by running the Gibbs sampler for 5,000 iterations with an initial burn-in of 1,000 iterations to lessen the effect of initial simulations. The process resulted in 4,000 final posterior samples for each regression parameter. Then the width of 95% Bayes credible interval was estimated for each selected quantile for each proposed method. The absolute bias is also estimated when the estimates of β^k were compared with the true value of β_k. The results are summarized in Tables 1–3.

Table 1 presents the results of all three proposed methods at several selected quantiles (i.e., τ = 0.25, τ = 0.50, τ = 0.75, and τ = 0.95). The table also presents the absolute bias and the width of 95% Bayes credible interval. The table shows that all proposed methods yielded very similar results. The table shows that, in general, Bayesian Adaptive Lasso quantile regression performs well compared to two other methods, BLQR and BQR.

Table 2 shows the results for autocorrelated errors. The table informs us that all proposed methods yielded almost similar values for corresponding quantiles. In general, BALQR performed best among the three methods because of the number of the smallest values of the entries.

Table 3 presents the absolute bias and the width of 95% Bayes credible interval for the three-parameter models of a nonnormal error condition. The results show that all values of the absolute Bias for corresponding quantiles are almost similar. Note that BALQR yielded the smallest values among the other three methods except for β^3. The Bayes credible intervals for all proposed methods at the same quantiles are almost similar. In general, the credible interval for BALQR is narrower than the other two methods.

The results concerning three different error distributions inform us of two things. First, although the performance of all the three methods proposed in this study is very close in general, the Bayes Adaptive Lasso quantile regression method performs better than the Bayes Lasso quantile regression and Bayes quantile method. Second, the BALQR method is robust to the error distribution assumptions, such as the normality assumption, the homogeneous assumption, or the non-correlated assumption. We assumed that the BALQR tends to produce more suitable values for the parameter estimated than other methods.

For the next analysis, we have to evaluate our algorithm's convergence used in the BALQR method. We use convergence diagnostics such as the posterior plots (trace plot and density plot) and autocorrelation analysis. Figures 1–3 present the trace plot, density plot, and autocorrelation plot of β^1 at quantile τ = 0.25 for all three conditions of error, i.e., heteroscedastic normal, autocorrelated error, and non-normal error, respectively. The author saves other plots for limited space.

As shown in the trace plots in Figures 1A, 2A, and 3A that all generated samples lie within two parallel horizontal lines, centered at respective values, and no trends are detected. The histograms of marginal posteriors in Figures 1B, 2B, and 3B above inform us that the conditional posterior distributions are the desired stationary univariate normal. All posterior distributions shrink at the true parameter value (trace plot and density plot). Furthermore, Figures 1C, 2C, and 3C inform that the decrease in the empirical autocorrelation of posterior samples proves that the underlying chains are stationary.

The results obtained from these convergence diagnostics indicate that our algorithm used in the BALQR approach could produce adequate and acceptable values of the estimated parameter.