High-Dimensional Variable Selection for Quantile Regression based on Variational Bayesian Method

1. Introduction

Quantile regression, as introduced by Koenker and Bassett (1978) [1], is an important statistical inference about the relationship between quantiles of the response distribution and available covariates, can offer a practically significant alternative to traditional mean regression because it provides a more comprehensive description of the response distribution than the mean. Moreover, quantile regression can capture the heterogeneous impact of regressors on different parts of the distribution [2], has excellent computational properties [3], exhibits robustness to outliers and has wide applicability [4]. For these reasons, quantile regression has attracted extensive attention in the literature. For example, see [5] for Bayesian quantile regression with asymmetric Laplace distribution to specify the likelihood, [6] for Bayesian Nonparametric approach to inference for quantile regression, [7] for a mechanism for Bayesian inference of quantile regression models, [8] for model selection in quantile regression, among others.

Although there is a growing literature on quantile regression, to the best of our knowledge, little of the existing quantile regression models focused on high-dimensional data with variable number is larger sample size. In practice, a large number of variables may be collected and some of them are insignificant and should be excluded from the final model. In the last two decades, there has been active methodological research in penalized methods for significant variable selection in linear parametric models. For example, see [9] for ridge regression, [10] for the least absolute shrinkage and selection operator (Lasso), [11] for smoothly clipped absolute deviation penalty (SCAD), [12] for elastic net penalty, [13] for adaptive lasso. These methods have been extended in quantile regression, see [14] for the $L_1$-regularization method for quantile regression, [15] for variable selection of quantile regression using SCAD and adaptive-Lasso penalties, however the aforementioned regularization methods are computationally complex and unstable, and they do not consider the prior information of parameters, which may lead to unsatisfactory parametric estimation accuracy. Bayesian approaches for variable selection and parameter estimation have received much attention over the past decades, because they can largely improve the accuracy and efficiency of parametric estimation by imposing various priors on model parameters, consistently select important variables and provide more information for variable selection than the corresponding penalization method with a highly non-conbex optimization problem. For example, see [16] for Bayesian Lasso where $L_1$ peanalty is involved in Laplace prior, [17] for the Bayesian form of adaptive Lasso, [18] for Bayesian Lasso quantile regression (BLQR), and [19] for Bayesian adaptive Lsso quantile regression (BALQR). The above-mentioned literature involves the implementation of the standard Gibbs sampler for posterior computation, which is no so scalable for high-dimensional data where number of variables is large compared with the sample size [20].

To address the issue, the Bayesian variable selection method with spike-and-slab prior [20] is favored by researchers, and it can be applied to high-dimensional data at cost of heavy computation burden. As a computationally efficient alternative to Markov chain Monte Carlo (MCMC) simulation, variational Bayes (VB) methods are increasingly used in machine learning and statistics for approximating posterior distributions in Bayesian inference. High-efficient variational Bayesian Spike-and-Slab Lasso(VBSSL) has been explored for some high-dimensional models. Ray and Szabo (2022) [21] used VBSSL method in high-dimensional linear model with regression coefficient’s prior specified as the mixture of Laplace distribution and Dirac mass. Yi and Tang (2022) [22] used VBSSL technology in the high-dimensional linear mixed models with interesting parameter’s prior as the mixture of two Laplace distributions. However, to the best of our knowledge, there is little work done on VBSSL mothod in quantile regression. Xi et al. (2016) [23] consider Bayesian variable selection in the nonparametric quantile regression with small variable dimension, in which the spike-and-slab prior is chosen as the mixture of point mass at zero and normal distribution. In this paper, in order to alleviate the computational burden and quantify the parametric uncertainty, we propose the quantile regression with Spike-and-Slab Lasso penalty based on variational Bayesian (VBSSLQR), in which the prior is mixture of two Laplace distributions respectively with smaller or larger variance.

The main contributions of this paper are as follows. First, our proposed VBSSLQR can carry out variable selection for high-dimensional quantile regression with quite a low computation cost, without requiring the non-convex optimization and avoiding the curse of dimensionality problem. Second, compared to the mean regressions, our quantile approach provides a more systematic strategy for examining how covariates affect different quantiles of the response distribution. Third, the mean regression errors are often assumed to be sub-Gaussian in ultra-high dimensional regression, which is not required in our setting.

The rest of the paper is organized as follows. In Section 2, for high-dimensional data, we propose an efficient quantile regression with Spike-and-Slab Lasso penalty based on variational Bayes (VBSSLQR). In Section 3, we randomly generate high-dimensional data with n=200 and p=500 (excluding intercept items), and perform 500 simulation experiments to explore the performance of our algorithm, and compare it with other quantile regression methods (Bayesian and non-Bayesian) and non-quantile regression methods. The results show that our method is significantly superior to other approachs in the case of high-dimensional data. We applied VBSSLQR to the real dataset that it contains information about crime in various cities in the United States, and compared it with other quantile regression methods. The results show that our method also has well performance and excellent efficiency in the real data, and the relevant results are shown in Section 4. Some concluding remarks are given in Section 5. Technical details are presented in the Appendix.

2. Models and Methods

2.1. Quantile Regression

Consider a dataset of n independent subjects, For the $i_{th}$ subject, let $y_i$ be the response and ${\mathbf{x}}_i={(1,x_{i1},...,x_{ip})}^{\top }$ is a $(r+1)\times 1$ predictor vector, a simple linear regression model is defined as follows:

$$\begin{array}{c}\hfill y_i={\mathbf{x}}_i^{\top }\mathbf{\beta }+{\epsilon }_i,i=1,2,\cdots ,n,\end{array}$$

(1)

where $\mathbf{\beta }={({\beta }_0,{\beta }_1,\cdots ,{\beta }_r)}^{\top }$ is the regression coefficient vector with ${\beta }_0$ corresponding to the intercept terms, and ${\epsilon }_i$ represents error term with unknown distribution. It is usual to assume that the $\tau$th quantile of the random error term is 0, that is, $Q_{\tau }\left({\epsilon }_i\right)=0$ for $0<\tau <1$. According to this assumption, the form of ${\tau }_{th}$ quantile regression of model (1) is specified as follows:

$$\begin{array}{c}\hfill Q_{y_i}\left(\tau \right|{\mathbf{x}}_i)={\mathbf{x}}_i^{\top }\mathbf{\beta },\end{array}$$

(2)

where $Q_{y_i}\left(\tau \right|{\mathbf{x}}_i)$ is the inverse cumulative distribution function of $y_i$ given ${\mathbf{x}}_i$ evaluated at $\tau$. The estimate of the regression coefficient vector $\mathbf{\beta }$ in equation (2) is

$$\begin{array}{cc}\hfill \widehat{\mathbf{\beta }}& =\underset{\mathbf{\beta }\in {\mathbf{R}}^{r+1}}{\mathrm{argmin}}\left[\sum _{i:y_i\ge {\mathbf{x}}_i\mathbf{\beta }}\tau |y_i-{\mathbf{x}}_i^{\top }\mathbf{\beta }|+\sum _{i:y_i<{\mathbf{x}}_i\beta }(1-\tau )|y_i-{\mathbf{x}}_i^{\top }\mathbf{\beta }|\right]\hfill \\ \hfill & =\underset{\mathbf{\beta }\in {\mathbf{R}}^{r+1}}{\mathrm{argmin}}\sum _{i=1}^n{\rho }_{\tau }(y_i-{\mathbf{x}}_i^{\top }\mathbf{\beta }),\hfill \end{array}$$

(3)

where the loss function ${\rho }_{\tau }\left(\mu \right)=\mu \times (\tau -\mathrm{I}(\mu <0))$ with the indicator function $\mathrm{I}(·)$.

In light of [24] and [5], minimizing equation (3) is equivalent to maximizing the likelihood of n independent individuals with the $i_{th}$ one distributed as asymmetric Laplace distribution (ALD) specified as

$$\begin{array}{c}\hfill p\left(y_i\right|{\mathbf{x}}_i,\mathbf{\beta },\sigma ,\tau )=\frac{\tau (1-\tau )}{\sigma }\exp \left\{-{\rho }_{\tau }\left(\frac{y_i-{\mu }_i}{\sigma }\right)\right\},\end{array}$$

(4)

where the local parameter ${\mu }_i={\mathbf{x}}_i^{\top }\mathbf{\beta }$, the scale parameter $\sigma >0$ and the skewness parameter $\tau$ is between 0 and 1, obviously ALD is Laplace distribution when $\tau =0.5$ and ${\mu }_i=0$. However, it is computationally infeasible to carry out statistical inference based directly on equation (4) involving non-differentiable point ${\mu }_i$. Following [25], equation (4) can be rewritten as the following hierarchical fashion:

$$\left\{\begin{array}{c}{y}_i={\mathbf{x}}_i^{\top }\mathbf{\beta }+k_1{z}_i+\sqrt{k_2\sigma z_i}{\xi }_i,\hfill \\ z_i|\sigma \stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }\mathrm{Exp}\left(\frac{1}{\sigma }\right),\hfill \\ {\xi }_i\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }N(0,1),\hfill \\ z_i\mathrm{is}\mathrm{independent}\mathrm{of}{\xi }_i,\hfill \end{array}\right.$$

(5)

where $k_1=\frac{1-2\tau }{\tau (1-\tau )}$, $k_2=\frac{2}{\tau (1-\tau )}$ and $\mathrm{Exp}\left(\frac{1}{\sigma }\right)$ denotes the exponential distribution with mean $\sigma$, whose specific density function is $p\left(z_i\right|\sigma )=\frac{1}{\sigma }exp\left(-\frac{1}{\sigma }{z}_i\right)\mathrm{I}(z_i>0)$. Equation (5) illustrates that an asymmetric Laplace distribution can also be represented as a mixture of exponential and standard normal distributions, which allows us to express a quantile regression model as a normal regression model in which the response has following conditional distribution

$$y_i|{\mathbf{x}}_i,\mathbf{\beta },z_i,\sigma \stackrel{\mathrm{ind}}{\sim }N\left({\mathbf{x}}_i^{\top }\mathbf{\beta }+k_1{z}_i,k_2\sigma z_i\right).$$

For the above-defined quantile regression model with high-dimensional covariate vector (r is enough large), it is of major interest to estimate parameter vector $\mathbf{\beta }$ and to identify the important covariates. To this end, we consider Bayesian Quantile Regression based on Spike-and-Slab Lasso as follows.

2.2. Bayesian Quantile Regression based on Spike-and-Slab Lasso

As early as 2016, Xi et al. [23] applied Spike-and-Slab prior to Bayesian quantile regression, but their proposed prior was a mixture of zero particle and normal distribution with large variance, and the estimate of posterior density was obtained through Gibbs sampler. In order to provide new theoretical insights for a class of continuous spike-and-slab priors, Rockova (2018) [26] introduced the novel family of spike-and-slab priors, which is a mixture of two density functions respectively with spike or slab probability. In this paper, we adopt spike-and-slab Lasso prior with a mixture of two Laplace distributions respectively with large or small variance [26], which facilitate the variational Bayesian technique to approximate the posterior density of parameters and to improve the efficiency of the algorithm. In light of ref [26], given the indicator ${\gamma }_j=0$ or 1, the prior of $\mathbf{\beta }$ in Bayesian quantile regression model (5) can be written as:

$$\begin{array}{c}\hfill \pi \left(\mathbf{\beta }\right|\mathbf{\gamma })=\prod _{j=0}^r\pi \left({\beta }_j\right|{\gamma }_j)=\prod _{j=0}^r\left[{\gamma }_j{\mathsf{\Psi }}_0\left({\beta }_j\right|{\lambda }_0)+(1-{\gamma }_j){\mathsf{\Psi }}_1\left({\beta }_j\right|{\lambda }_1)\right],\end{array}$$

(6)

where Laplace density ${\mathsf{\Psi }}_0\left({\beta }_j\right|{\lambda }_0)=\frac{{\lambda }_0}{2}exp(-{\lambda }_0|{\beta }_j\left|\right)$ and ${\mathsf{\Psi }}_1\left({\beta }_j\right|{\lambda }_1)=\frac{{\lambda }_1}{2}exp(-{\lambda }_1|{\beta }_j\left|\right)$ with precision parameters ${\lambda }_0$ and ${\lambda }_1$ satisfying that ${\lambda }_0\gg {\lambda }_1$, the indicator variable set $\mathbf{\gamma }=\left\{{\gamma }_j\right|j=0,1,\cdots ,r\}$, the $j_{th}$ variable is active when ${\gamma }_j=0$, and inactive otherwise. Similar to [27], Laplace distribution for the regression coefficient ${\beta }_j$ can be represented as a mixture of a normal distribution and an exponential distribution, specifically the distribution of ${\beta }_j$ can be expressed as a hierarchical structure as follows:

$$\begin{array}{cc}\hfill {\beta }_j|h_{0j}^2,h_{1j}^2,{\gamma }_j& \stackrel{\mathrm{ind}.}{\sim }{\gamma }_{j}N\left(0,h_{0j}^2\right)+\left(1-{\gamma }_j\right)N\left(0,h_{1j}^2\right),\hfill \\ \hfill h_{0j}^2|{\lambda }_0^2& \stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }\mathrm{Exp}\left(\frac{{\lambda }_0^2}{2}\right),\hfill \\ \hfill h_{1j}^2|{\lambda }_1^2& \stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }\mathrm{Exp}\left(\frac{{\lambda }_1^2}{2}\right),\hfill \\ \hfill {\gamma }_j& \stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }B\left({\pi }_{\gamma }\right),\hfill \end{array}$$

(7)

where $B\left({\pi }_{\gamma }\right)$ denotes Bernoulli distribution with ${\pi }_{\gamma }$ being the probability that indicator variable ${\gamma }_j$ equals one for $j=0,1,\cdots ,r$, and specify the prior of ${\pi }_{\gamma }$ as Beta distribution $\mathrm{Be}\left(a_{{\pi }_{\gamma }},b_{{\pi }_{\gamma }}\right)$ with hyperparameters $a_{{\pi }_{\gamma }}$ and $b_{{\pi }_{\gamma }}$, ${\lambda }_0^2$ and ${\lambda }_1^2$ are regularization parameters to identify important variables, for which we consider the following conjugate priors:

$$\begin{array}{ccc}{\lambda }_0^2\hfill & \sim \mathrm{Ga}\left({\nu }_{{\lambda }_0},1\right),{\lambda }_1^2\hfill & \sim \mathrm{Ga}\left({\nu }_{{\lambda }_1},1\right),\hfill \end{array}$$

where $\mathrm{Ga}\left(a,b\right)$ denotes the gamma distribution with shape parameter a and scale parameter b. As mentioned above, ${\lambda }_0$ and ${\lambda }_1$ should satisfy that ${\lambda }_0\gg {\lambda }_1$, to this end, we select hyperparmeters ${\nu }_{{\lambda }_0}$ and ${\nu }_{{\lambda }_1}$ to satisfy that ${\nu }_{{\lambda }_0}\gg {\nu }_{{\lambda }_1}$. The prior of scale parameter $\sigma$ in (5) is inverse gamma distribution $\mathrm{IG}\left(a_{\sigma },b_{\sigma }\right)$ the hyperparmeters $a_{\sigma }=1$ and $b_{\sigma }=0.01$ in the paper leading to almost non-informative prior.

Under Bayesian statistical paradise, based on the above priors and likelihood of quantile regression, it is direct to induce the following posterior distribution $\pi \left(\mathbf{\theta }\right|\mathbf{D})\propto p(\mathbf{\theta },\mathbf{D})$, where $\mathbf{\theta }=\{\mathbf{\beta },\mathbf{z},\sigma ,\mathbf{\gamma },{\mathbf{h}}_0^2,{\mathbf{h}}_1^2,{\lambda }_0^2,{\lambda }_1^2,{\pi }_{\gamma }\}$, the latent variable set $\mathbf{z}=\left\{z_i\right|i=1,\cdots ,n\}$, ${\mathbf{h}}_0^2=\left\{h_{0j}^2\right|j=0,1,\cdots ,r\}$, ${\mathbf{h}}_1^2=\left\{h_{1j}^2\right|j=0,1,\cdots ,r\}$, observing set $\mathbf{D}=\{\mathbf{y},\mathbf{x}\}$ with the response set $\mathbf{y}=\left\{y_i\right|i=1,\cdots ,n\}$ and covariate set $\mathbf{x}=\left\{x_i\right|i=1,\cdots ,n\}$. Based on the hierarchical structure (5) of the quantile regression likelihood $p\left(y_i\right|{\mathbf{x}}_i,\mathbf{\beta })$ and the hierarchical structure (7) of the Spike-and-Slab prior for regression coefficient vector $\mathbf{\beta }$, we derive the joint density

$$\begin{array}{cc}\hfill p(\mathbf{\theta },\mathbf{D})& =\prod _{i=1}^{n}N\left(y_i\right|{\mathbf{x}}_i^{\top }\mathbf{\beta }+k_1{z}_i,k_2\sigma z_i)\mathrm{Exp}\left(z_i\right|{\sigma }^{-1})\mathrm{Ga}\left({\lambda }_0^2\right|{\nu }_{{\lambda }_0},1)\mathrm{Ga}\left({\lambda }_1^2\right|{\nu }_{{\lambda }_1},1)\hfill \\ \hfill & \times \prod _{j=0}^r{\left[N\left({\beta }_j\right|0,h_{0j}^2)\right]}^{{\gamma }_j}{\left[N\left({\beta }_j\right|0,h_{1j}^2)\right]}^{(1-{\gamma }_j)}\mathrm{Exp}\left(h_{0j}^2\right|\frac{{\lambda }_0^2}{2})\mathrm{Exp}\left(h_{1j}^2\right|\frac{{\lambda }_1^2}{2})\hfill \\ \hfill & \times \prod _{j=0}^{r}B\left({\gamma }_j\right|1,{\pi }_{\gamma })\mathrm{Be}\left({\pi }_{\gamma }\right|a_{{\pi }_{\gamma }},b_{{\pi }_{\gamma }})\mathrm{IG}\left(\sigma \right|a_{\sigma },b_{\sigma })\hfill \\ \hfill & \propto \prod _{i=1}^n{\left(k_2\sigma z_i\right)}^{-\frac{1}{2}}\exp \left(-\frac{{(y_i-{\mathbf{x}}_i^{\top }\mathbf{\beta }-k_1{z}_i)}^2}{2k_2\sigma z_i}\right)\frac{1}{\sigma }\exp \left(-\frac{z_i}{\sigma }\right){\left({\lambda }_0^2\right)}^{v_{{\lambda }_0}-1}\exp (-{\lambda }_0^2)\hfill \\ \hfill & \times {\left({\lambda }_1^2\right)}^{(v_{{\lambda }_1}-1)}\exp (-{\lambda }_1^2)\prod _{j=0}^r{\left[\frac{1}{\sqrt{h_{0j}^2}}\exp \left(-\frac{{\beta }_j^2}{2h_{0j}^2}\right)\right]}^{{\gamma }_j}{\left[\frac{1}{\sqrt{h_{1j}^2}}\exp \left(-\frac{{\beta }_j^2}{2h_{1j}^2}\right)\right]}^{(1-{\gamma }_j)}\hfill \\ \hfill & \times \frac{{\lambda }_0^2}{2}\exp \left(-\frac{{\lambda }_0^2}{2}{h}_{0j}^2\right)\frac{{\lambda }_1^2}{2}\exp \left(-\frac{{\lambda }_1^2}{2}{h}_{1j}^2\right)\prod _{j=0}^r{\pi }_{\gamma }^{{\gamma }_j}{(1-{\pi }_{\gamma })}^{1-{\gamma }_j}{\pi }_{\gamma }^{a_{{\pi }_{\gamma }}-1}{(1-{\pi }_{\gamma })}^{b_{{\pi }_{\gamma }}-1}\hfill \\ \hfill & \times \frac{b_{\sigma }^{a_{\sigma }}}{\Gamma \left(a_{\sigma }\right)}{\sigma }^{-a_{\sigma }-1}\exp \left(-\frac{b_{\sigma }}{\sigma }\right).\hfill \end{array}$$

(8)

Although it is easy to sample from the above presented posterior can, it is very time-consuming for a sufficiently large r. To address the issue, we investigate a fast yet efficient approach as follows, i.e., variational Bayesian method.

2.3. Quantile Regression with Spike-and-Slab Lasso Penalty Based on Variational Bayesian

At present, the most commonly variational Bayesian approximation posterior distribution methods is mean field approximation theory [28] which is of the highest efficiency among variational methods, especially for those parameters or parameter block with conjugate priors. In view of the fact that Bayesian quantile regression needs to take into account that the variance of each observation value is different, and each $y_i$ corresponds to a potential variable $z_i$, which will result in the algorithm efficiency of quantile regression being lower than that of general mean regression. Therefore, in this paper, we use the variational Bayesian algorithm of the mean field approximation, which is the most efficient algorithm, to derive the quantile regression model with Spike-and-Slab Lasso penalty.

According to the principle of variational inference, we require specifying a variational family $\mathcal{F}$ of densities for random variables $\mathbf{\theta }$ having the same support $\Theta$ as the posterior density $\pi \left(\mathbf{\theta }\right|\mathbf{D})$. Let $\mathrm{q}\left(\mathbf{\theta }\right)\in \mathcal{F}$ be any variational density approximating the posterior density $\pi \left(\mathbf{\theta }\right|\mathbf{D})$. The main purpose of variational Bayes is to find the best approximation to $\pi \left(\mathbf{\theta }\right|\mathbf{D})$ via the Kullback-Leibler divergence between $\mathrm{q}\left(\mathbf{\theta }\right)$ and $\pi \left(\mathbf{\theta }\right|\mathbf{D})$, which is equivalent to solving the following optimization problem:

$${\mathrm{q}}^{*}\left(\mathbf{\theta }\right)=\underset{\mathrm{q}\left(\mathbf{\theta }\right)\in \mathcal{F}}{\mathrm{argmin}}\mathrm{KL}\left[\mathrm{q}\left(\mathbf{\theta }\right)\Vert \pi \left(\mathbf{\theta }\right|\mathbf{D})\right],$$

where $\mathrm{KL}\left[\mathrm{q}\left(\mathbf{\theta }\right)\Vert \pi \left(\mathbf{\theta }\right|\mathbf{D})\right]={\int }_{\Theta }\mathrm{q}\left(\mathbf{\theta }\right)\log \frac{\mathrm{q}\left(\mathbf{\theta }\right)}{\pi \left(\mathbf{\theta }\right|\mathbf{D})}d\mathbf{\theta }$, which is not less than zero and equal to zero if and only if $\mathrm{q}\left(\mathbf{\theta }\right)\equiv \pi \left(\mathbf{\theta }\right|\mathbf{D})$. The posterior density $\pi \left(\mathbf{\theta }\right|\mathbf{D})=\frac{p(\mathbf{\theta },\mathbf{D})}{p\left(\mathbf{D}\right)}$ with the joint distribution $p(\mathbf{\theta },\mathbf{D})$ of parameter $\mathbf{\theta }$ and data $\mathbf{D}$ and the marginal distribution $p\left(\mathbf{D}\right)$ data $\mathbf{D}$. Since $p\left(\mathbf{D}\right)={\int }_{\Theta }p(\mathbf{\theta },\mathbf{D})d\mathbf{\theta }$, which has not analytic expression for our considered model, it is rather difficult to implement the optimization problem presented above. It is easy to induce that

$$\log p\left(\mathbf{D}\right)=\mathrm{KL}\left[\mathrm{q}\left(\mathbf{\theta }\right)\Vert \pi \left(\mathbf{\theta }\right|\mathbf{D})\right]+\mathcal{L}\left\{\mathrm{q}\left(\mathbf{\theta }\right)\right\}\ge \mathcal{L}\left\{\mathrm{q}\left(\mathbf{\theta }\right)\right\},$$

in which the evidence lower bound (ELOB) $\mathcal{L}\left\{\mathrm{q}\left(\mathbf{\theta }\right)\right\}={\mathrm{E}}_{\mathrm{q}\left(\mathbf{\theta }\right)}\left[\log p(\mathbf{\theta },\mathbf{D})\right]-{\mathrm{E}}_{\mathrm{q}\left(\mathbf{\theta }\right)}\left[\log \mathrm{q}\left(\mathbf{\theta }\right)\right]$ with ${\mathrm{E}}_{\mathrm{q}\left(\mathbf{\theta }\right)}\left[·\right]$ representing the expectation taken with respect to the variational density $\mathrm{q}\left(\mathbf{\theta }\right)$. Thus, minimizing $\mathrm{KL}\left[\mathrm{q}\left(\mathbf{\theta }\right)\Vert \pi \left(\mathbf{\theta }\right|\mathbf{D})\right]$ is equivalent to maximizing $\mathcal{L}\left\{\mathrm{q}\right(\mathbf{\theta }\left)\right\}$ because log$p\left(\mathbf{D}\right)$ does not depend on $\mathbf{\theta }$. That is,

$${\mathrm{q}}^{*}\left(\mathbf{\theta }\right)=\underset{\mathrm{q}\left(\mathbf{\theta }\right)\in \mathcal{F}}{\mathrm{argmin}}\mathrm{KL}\left[\mathrm{q}\left(\mathbf{\theta }\right)\Vert \pi \left(\mathbf{\theta }\right|\mathbf{D})\right]=\underset{\mathrm{q}\left(\mathbf{\theta }\right)\in \mathcal{F}}{\mathrm{argmax}}\mathcal{L}\left\{\mathrm{q}\left(\mathbf{\theta }\right)\right\},$$

which implies that finding the best approximation problem for $\pi \left(\mathbf{\theta }\right|\mathbf{D})$ is transformed into maximizing $\mathcal{L}\left\{\mathrm{q}\right(\mathbf{\theta }\left)\right\}$ over the variational family $\mathcal{F}$. The complexity of the optimization procedure heavily depends on that of the variational family $\mathcal{F}$. Hence, it is quite desirable to specify a relatively simple variational family $\mathcal{F}$ to optimize the objective function $\mathcal{L}\left\{\mathrm{q}\right(\mathbf{\theta }\left)\right\}$ with respect to $q\left(\mathbf{\theta }\right)$ over $\mathcal{F}$.

Following the commonly used approach to specify a simple variational family $\mathcal{F}$ in the variational literature, we consider the well-known mean-field variational family in which blocks of $\mathbf{\theta }$ are mutually independent and each is dominated by a distinct factor in the variational density. That is, the variational density $\mathrm{q}\left(\mathbf{\theta }\right)$ is assumed to be factorized across the components of $\mathbf{\theta }$:

$$\mathrm{q}\left(\mathbf{\theta }\right)=\prod _{j=0}^r\mathrm{q}({\beta }_j,{\gamma }_j)\mathrm{q}\left(h_{0j}^2\right)\mathrm{q}\left(h_{1j}^2\right)\prod _{i=1}^n\mathrm{q}\left(z_i\right)\mathrm{q}\left({\lambda }_0^2\right)\mathrm{q}\left({\lambda }_1^2\right)\mathrm{q}\left({\pi }_{\gamma }\right)\equiv \prod _s^S{\mathrm{q}}_s\left({\mathbf{\theta }}_s\right),$$

where the individual variational densities ${\mathrm{q}}_s\left({\mathbf{\theta }}_s\right)$’s are unspecified, but the above assumed factorization across components is pre-given. Moreover, the optimal solutions of ${\mathrm{q}}_s\left({\mathbf{\theta }}_s\right)$’s are to be determined by maximizing $\mathcal{L}\left\{\mathrm{q}({\mathbf{\theta }}_1,\cdots ,{\mathbf{\theta }}_S)\right\}$ with respect to densities ${\mathrm{q}}_1\left({\mathbf{\theta }}_1\right),\cdots ,{\mathrm{q}}_S\left({\mathbf{\theta }}_S\right)$ via the coordinate ascent method, where $\mathbf{\theta }=\{{\mathbf{\theta }}_1,\cdots ,{\mathbf{\theta }}_S\}$ where ${\mathbf{\theta }}_s$ can be either a scalar or a vector. This means that when the correlation between several unknown parameters or potential variables cannot be ignored, they should be put in the same block and merged into ${\mathbf{\theta }}_s$. Following the idea of the coordinate ascent method given in ref [29], when fixing other variational factors ${\mathrm{q}}_k\left({\mathbf{\theta }}_k\right)$ for $k\ne s$, i.e., the optimal density ${\mathrm{q}}_s^{*}\left({\mathbf{\theta }}_s\right)$, that maximizes $\mathcal{L}\left\{\mathrm{q}\right(\mathbf{\theta }\left)\right\}$ with respect to ${\mathrm{q}}_s\left({\mathbf{\theta }}_s\right)$, is shown to take the form

$$\begin{array}{c}\hfill {\mathrm{q}}_s^{*}\left({\mathbf{\theta }}_s\right)\propto \exp \left\{{\mathrm{E}}_{-{\theta }_s}\left[\log p(\mathbf{\theta },\mathbf{D})\right]\right\},\end{array}$$

(9)

where $\log p(\mathbf{\theta },\mathbf{D})$ is the logarithm of the joint density function and ${\mathrm{E}}_{-s}[·]$ is the expectation taken with respect to the density ${\prod }_{k\ne s}{\mathrm{q}}_k\left({\mathbf{\theta }}_k\right)$ for $s=1,\cdots ,S$.

According to equations (8)-(9), we can derive the variational posterior for each parameter as follows (see Appendix A for the details) :

$$\begin{array}{cc}\hfill {\beta }_j& \stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }{\mu }_{{\gamma }_j}N\left({\mu }_{0j},{\sigma }_{0j}^2\right)+(1-{\mu }_{{\gamma }_j})N\left({\mu }_{1j},{\sigma }_{1j}^2\right),\mathrm{for}j=1,\cdots ,r,\hfill \\ \hfill h_{0j}^2& \stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }{\mu }_{{\gamma }_j}\mathrm{GIG}\left(\frac{1}{2},{\mathrm{E}}_{{\lambda }_0^2}\left({\lambda }_0^2\right),{\sigma }_{0j}^2+{\mu }_{0j}^2\right)+(1-{\mu }_{{\gamma }_j})\mathrm{Exp}\left(\frac{1}{2}{\mathrm{E}}_{{\lambda }_0^2}\left({\lambda }_0^2\right)\right),\hfill \\ \hfill h_{1j}^2& \stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }(1-{\mu }_{{\gamma }_j})\mathrm{GIG}\left(\frac{1}{2},{\mathrm{E}}_{{\lambda }_1^2}\left({\lambda }_1^2\right),{\sigma }_{1j}^2+{\mu }_{1j}^2\right)+{\mu }_{{\gamma }_j}\mathrm{Exp}\left(\frac{1}{2}{\mathrm{E}}_{{\lambda }_1^2}\left({\lambda }_1^2\right)\right),\hfill \\ \hfill {\lambda }_0^2& \sim \mathrm{Ga}\left(r+1+{\nu }_{{\lambda }_0},1+\frac{1}{2}\sum _{j=1}^r{\mathrm{E}}_{h_{0j}^2}\left(h_{0j}^2\right)\right),\hfill \\ \hfill {\lambda }_1^2& \sim \mathrm{Ga}\left(r+1+{\nu }_{{\lambda }_1},1+\frac{1}{2}\sum _{j=1}^r{\mathrm{E}}_{h_{1j}^2}\left(h_{1j}^2\right)\right),\hfill \\ \hfill {\gamma }_j& \stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }B\left(1,{(1+{\mathrm{e}}^{{\zeta }_i})}^{-1}\right),\mathrm{for}j=0,1,\cdots ,r,\hfill \\ \hfill {\pi }_{\gamma }& \sim \mathrm{Be}\left(a+\sum _{j=0}^r{\mu }_{{\gamma }_j},r+1+b-\sum _{j=0}^r{\mu }_{{\gamma }_j}\right),\hfill \\ \hfill z_i& \stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }\mathrm{GIG}\left(\frac{1}{2},a_{z_i},b_{z_i}\right),\mathrm{for}i=1,2,\cdots ,n,\hfill \\ \hfill \sigma & \sim \mathrm{IG}\left(\frac{3}{2}n+a_{\sigma },c_{\sigma }\right),\hfill \end{array}$$

where $\mathrm{GIG}(·,·,·)$ denotes generalized inverse Gaussian distribution,

$$\begin{array}{cc}\hfill & {\mu }_{0j}={\sigma }_{oj}^2\frac{{\mathrm{E}}_{\sigma }\left({\sigma }^{-1}\right)}{k_2}\sum _{i=1}^n\left[{\mathrm{E}}_{z_i}\left(z_i^{-1}\right)\left(y_i-{\mathbf{x}}_{i,(-j)}^{\top }{\mathbf{\mu }}_{{\mathbf{\beta }}_{(-j)}}\right)-k_1\right]x_{ij},\hfill \\ \hfill & {\sigma }_{0j}^2={\left[{\mathrm{E}}_{h_{0j}^2|{\gamma }_j=1}\left(h_{0j}^{-2}\right)+\frac{{\mathrm{E}}_{\sigma }\left({\sigma }^{-1}\right)}{k_2}\sum _{i=1}^n{x}_{ij}^2{\mathrm{E}}_{z_i}\left(z_i^{-1}\right)\right]}^{-1},\hfill \\ \hfill & {\mu }_{1j}={\sigma }_{1j}^2\frac{{\mathrm{E}}_{\sigma }\left({\sigma }^{-1}\right)}{k_2}\sum _{i=1}^n\left[{\mathrm{E}}_{z_i}\left(z_i^{-1}\right)\left(y_i-{\mathbf{x}}_{i,(-j)}^{\top }{\mathbf{\mu }}_{{\mathbf{\beta }}_{(-j)}}\right)-k_1\right]x_{ij},\hfill \\ \hfill & {\sigma }_{1j}^2={\left[{\mathrm{E}}_{h_{1j}^2|{\gamma }_j=0}\left(h_{1j}^{-2}\right)+\frac{{\mathrm{E}}_{\sigma }\left({\sigma }^{-1}\right)}{k_2}\sum _{i=1}^n{x}_{ij}^2{\mathrm{E}}_{z_i}\left(z_i^{-1}\right)\right]}^{-1},\hfill \\ \hfill \end{array}$$

(10)

and

$$\begin{array}{cc}\hfill & {\zeta }_j={\mathrm{E}}_{{\pi }_{\gamma }}\left(\log (1-{\pi }_{\gamma })\right)-{\mathrm{E}}_{{\pi }_{\gamma }}\left(\log {\pi }_{\gamma }\right)+\frac{1}{2}{\mathrm{E}}_{h_{0j}^2|{\gamma }_j=1}\left(\log h_{0j}^2\right)-\frac{1}{2}{\mathrm{E}}_{h_{1j}^2|{\gamma }_j=0}\left(\log h_{1j}^2\right)\hfill \\ \hfill & +\frac{1}{2}\{\log {\sigma }_{1j}^2-\log {\sigma }_{0j}^2\}+\frac{1}{2}(\log {\sigma }_{1j}^2-\log {\sigma }_{0j}^2)+\frac{1}{2}({\mu }_{1j}^2{\left({\sigma }_{1j}^2\right)}^{-1}-{\mu }_{0j}^2{\left({\sigma }_{0j}^2\right)}^{-1}),\hfill \\ \hfill & a_{z_i}=\frac{{\mathrm{E}}_{\sigma }\left({\sigma }^{-1}\right)}{k_2}\left(k_1^2+2k_2\right),\hfill \\ \hfill & b_{z_i}=\frac{{\mathrm{E}}_{\sigma }\left({\sigma }^{-1}\right)}{k_2}[{(y_i-{\mathbf{x}}_i^{\top }{\mathbf{\mu }}_{\mathbf{\beta }})}^2+{\mathbf{x}}_i^{\top }{\Sigma }_{\mathbf{\beta }}{\mathbf{x}}_i],\hfill \\ \hfill & c_{\sigma }=\frac{1}{2k_2}\sum _{i=1}^n\{k_1^2{\mu }_{z_i}+{\mathrm{E}}_{z_i}\left(z_i^{-1}\right)[{(y_i-{\mathbf{x}}_i^{\top }{\mathbf{\mu }}_{\mathbf{\beta }})}^2+{\mathbf{x}}_i^{\top }{\Sigma }_{\mathbf{\beta }}{\mathbf{x}}_i]\hfill \\ \hfill & -2k_1(y_i-{\mathbf{x}}_i^{\top }{\mathbf{\mu }}_{\mathbf{\beta }})\}+\sum _{i=1}^n{\mu }_{z_i}+b_{\sigma }.\hfill \end{array}$$

In the above equation, ${\mu }_{{\lambda }_0^2}={\mathrm{E}}_{{\lambda }_0^2}\left({\lambda }_0^2\right)$, and ${\mu }_{{\lambda }_1^2}$, ${\mu }_{{\beta }_j}$, ${\mu }_{{\gamma }_j}$, ${\mu }_{z_i}$ are similar to ${\mu }_{{\lambda }_0^2}$, ${\mathbf{\mu }}_{\mathbf{\beta }}={({\mu }_{{\beta }_0},{\mu }_{{\beta }_1},\cdots ,{\mu }_{{\beta }_p})}^{\top }$ with ${\mu }_{{\beta }_j}={\mu }_{{\gamma }_j}{\mu }_{0j}+(1-{\mu }_{{\gamma }_j}){\mu }_{1j}$, and ${\mathbf{\mu }}_{{\mathbf{\beta }}_{(-j)}}$ is a $p\times 1$ vector with the $j_{th}$ component of vector ${\mathbf{\mu }}_{\mathbf{\beta }}$ deleted, ${\Sigma }_{\mathbf{\beta }}=\mathrm{diag}\left({\sigma }_{{\beta }_0}^2{\sigma }_{{\beta }_1}^2,\cdots {\sigma }_{{\beta }_p}^2\right)$ with ${\sigma }_{{\beta }_j}^2={\mu }_{{\gamma }_j}{\sigma }_{0j}^2+(1-{\mu }_{{\gamma }_j}){\sigma }_{1j}^2+{\mu }_{{\gamma }_j}(1-{\mu }_{{\gamma }_j}){({\mu }_{0j}-{\mu }_{1j})}^2$.

In the section above, we have derived the variational posterior of each parameter. Using the idea of coordinate axis optimization, we can update each variational distribution iteratively until it converges.

For this reason, we list the variational Bayesian Spike-and-Slab Lasso quantile regression (VBSSLQR) algorithm as shown in the following Algorithm 1:

In Algorithm 1 above, the $\mathsf{\Psi }(·)$ is the digamma fuction and the expectation ${\mathrm{E}}_{h_{0j}^2|{\gamma }_j=1}^{(t+1)}\left(\log h_{0j}^2\right)$ of $\log h_{0j}^2$ with respect to generalized inverse Gaussian distribution. So, we assume that $x\sim \mathrm{GIG}(p,a,b)$, then:

$$\begin{array}{c}\hfill {\mathrm{E}}_x\left(\log x\right)=\frac{\frac{\mathrm{d}{\mathrm{K}}_p\left(\sqrt{ab}\right)}{\mathrm{d}p}}{{\mathrm{K}}_p\left(\sqrt{ab}\right)}-\frac{1}{2}\ln \left(\frac{a}{b}\right),\end{array}$$

(11)

where ${\mathrm{K}}_p(·)$ represents Bessel function of the second kind. Note that there is no analytic solution or function to the differential of the Modified Bessel function. Therefore, We approcimate ${\mathrm{E}}_x\left(\log x\right)$ by the second-order Taylor expansion of $\log x$ . This paper lists the expectation of some parameter functions about variational posteriors involved in algorithm Section 2.3. See Appendix B for details. Based on our proposed VBSSLQR algorithm, in the next section, we randomly generate high-dimensional data, conducting simulation research and comparing the performance of other methods.

3. Simulation Studies

In this section, we use simulated high-dimensional data with sample size $n=200$ and variable number $r=500$, in order to study the performance of VBSSLQR and compare it with existing methods, including linear regression with the Lasso penalty (Lasso), linear regression with the adaptive Lasso penalty (ALasso), quantile regression with the Lasso penalty (QRL), quantile regression with the adaptive Lasso penalty (QRAL), Bayesian regularized quantile regression with the Lasso penalty (BLQR) and Bayesian regularized quantile regression with the adaptive Lasso penalty (BALQR). The data in the simulation studies are generated by equation (1), in which the covariate vector ${\mathbf{x}}_i$ is randomly generated from the multivariate normal distribution $N(\mathbf{0},\Sigma )$ with the ${(k,l)}_{th}$ of $\Sigma$ being $0.5^{|k-l|}$. Among these covariates, we only consider ten important explanatory variables that have significant impact on the dependent variable. We set the 1, 51, 101, 151, 201, 251, 301, 351, 401, 451 predictors to be active and their regression coefficients are -3, -2.5, -2, -1.5, -1, 1, 1.5, 2, 2.5, 3, and the rest are zero. In addition, we discuss the performance of various approachs in the case of two kinds of random errors, namely, independent and identically distributed (i.i.d.) random errors and heterogenous random errors.

3.1. Independent and Identically Distributed Random Errors Random

In this subsection, with reference to [19] and [23], we set the random errors ${\epsilon }_i$’s in equation (1) to be independently and identically distributed and consider following five different distributions with $\tau$ quantile being zero for it:

• The error ${\epsilon }_i\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }N(-\mu ,1)$ with $\mu$ being the $\tau$ quantile of $N(0,1)$, for $i=1,\cdots ,n$;
• The error ${\epsilon }_i\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }\mathrm{Laplace}(-\mu ,1)$ with $\mu$ being the $\tau$ quantile of $\mathrm{Laplace}(0,1)$, where $\mathrm{Laplace}(a,b)$ denotes the Laplace distribution with location parameter a and scale parameter b;
• The error ${\epsilon }_i\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }0.1N(-{\mu }_1,9)+0.9N(-{\mu }_2,1)$ with ${\mu }_1$ and ${\mu }_2$ respectively being the $\tau$ quantile of $N(0,9)$ and $N(0,1)$;
• The error ${\epsilon }_i\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }0.1\mathrm{Laplace}$$(-{\mu }_1,9)+0.9\mathrm{Laplace}(-{\mu }_2,1)$ with ${\mu }_1$ and ${\mu }_2$ respectively being the $\tau$ quantile of $\mathrm{Laplace}(0,9)$ and $\mathrm{Laplace}(0,1)$;
• The error ${\epsilon }_i\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }\mathrm{Cauchy}(-\mu ,0.2)$ with $\mu$ being the $\tau$ quantile of $\mathrm{Cauchy}(0,0.2)$, where $\mathrm{Cauchy}(a,b)$ denotes the Cauchy distribution with location parameter a and scale parameter b;

For above any error distribution with any $\tau \in (0.3,0.5,0.7)$, we run 500 replications for each method, and evaluate the performance by two criterion. The first criterion is the median of mean absolute deviations (MMAD), whcih quantify the general distance between the estimated conditional quantile and the true conditional quantile. Specifically, mean absolute deviation (MAD) in any replication is defined as $\frac{1}{n}{\sum }_{i=1}^n|{\mathbf{x}}_i\mathbf{\beta }-{\mathbf{x}}_i{\widehat{\mathbf{\beta }}}_{\tau }\left)\right|,$ where ${\widehat{\mathbf{\beta }}}_{\tau }$ is the estimate of regression coefficient $\mathbf{\beta }$ given $\tau$. The second criterion is the mean of true positives (TP) and false positives (FP) selected by each method.

Table 1 shows the median and stndard deviation (SD) of MADs estimated by each method for simulations with homogeneous errors. It is clear that our method is either optimal (bold) in all cases, especially in quantile $\tau =0.3,0.7$, our approach is significantly superior to the other six methods. When the error distribution is normal, Laplace and normal mixture, the MMAD of BALQR method is suboptimal. When the error distribution is Cauchy, the MMAD of the QRL approach is suboptimal, but Table 2 shows that it sacrifices the complexity of the model. When the error distribution is Laplace mixture, the MMAD of the Lasso approach is suboptimal, however it selects the over-fitting quantile regression model with about 30 FP variables, as can be seen from the Table 2. It is particularly important to note that in the case of high dimensional data, MMAD for the quantile regression model with Lasso penalty or adaptive lasso penalty is not less than MMAD for general linear models with Lasso penalty. Therefore, it is inappropriate to still use Lasso penalty or adaptive lasso in this case.

In order to show the results of variable selection more intuitively, we introduced TP (True positives) and FP (False positives) to calculate the mean of TP and FP of 500 repeated simulations respectively. The detailed results are shown in Table 2 below:

It is not difficult to conclude from the results in Table 2: Lasso method can basically select all true active variables, but it fits many false active variables, and when the random errors is Cauchy distribution, it cannot select all true active variables; ALasso approach can identify true active variables, but there are also some misjudged behaviors, especially when the random errors distribution is Cauchy distribution; Although QRL approach can identify real active variables, it still contains many false active variables, and the model it identifies has high complexity; QRAL method can not identify all true active variables, and incorrectly identifies some inactive variables; In the case of high-dimensional data, BLQR approach cannot select all the true active variables, but it also does not select the inactive variables incorrectly; The true active variable selected by BALQR method is better than that by BLQR method, but some inactive variables are selected incorrectly, especially when the random errors distribution is Cauchy distribution; Our VBSSLQR method not only has the smallest MMAD, but also can select true active variables and eliminate most false active variables. As a whole, our method is superior to the other six methods in variable selection, especially in quantile $\tau =0.3,0.7$ it is significantly better than BLQR, BALQR and QRL.

3.2. Heterogenous Random Errors

Now we consider the case of heterogenous random errors to demonstrate the performance of our method. In this subsection, the data are generated from the following model:

$$\begin{array}{c}\hfill y_i={\mathbf{x}}_i^{\top }\mathbf{\beta }+(1+u_i){\epsilon }_i,\end{array}$$

(12)

where $u_i\stackrel{i.i.d.}{\sim }U(0,1)$ in which $U(a,b)$ represents the uniform distribution with support set $(a,b)$. The design matrix ${\mathbf{x}}_i$ is generated in the same way as above, and the regression coefficient $\mathbf{\beta }$ is set as before. Furthermore, in the simulation study, we combine ${\mathbf{x}}_i$ and $u_i$, that is, $u_i$ is also a covariate. Finally, random error ${\epsilon }_i$ is also generated from five different distributions defined in subsection (Section 3.1).

We also studied the performance of the quantile $\tau \in (0.3,0.5,0.7)$ under different methods and simulated 500 times to calculate the MMAD and mean of TP/FP. We list the experimental results in the following two tables.

Table 3. The Median and standard deviation of 500 MADs estimated via various methods for simulations with heterogenous random errors.

Table 3. The Median and standard deviation of 500 MADs estimated via various methods for simulations with heterogenous random errors.

In heterogenous random errors, our approach is still the best one similar to i.i.d. random errors. It is noteworthy that the VBSSLQR method is more robust than other methods, and our method has the smallest MMAD change compared to the case i.i.d. random errors. We can see that the MMAD of our method remains basically unchanged, while the MMAD of the other six methods differs even more than 0.5 compared to i.i.d. random errors in some states.

We also investigate the mean of TP and FP in heterogenous random errors under different quantiles $\tau$. The results are listed in Table 4, whichs show that the effect for variable selection of heterogenous random errors is slightly lower than the effect for variable selection of i.i.d. random errors under the same sample size, but our method still provides the best selection result.

We also calculate the mean execution times of various Bayesian quantile regressions under different quantile $\tau$ for different distributions of random errors, and list the results in Table A1 of Appendix C, which illustrate our proposed VBSSLQR approach is a lot more efficient than BLQR and BALQR , for which we sampled MCMC 1000 times and discarded the first 500 times, and made statistical inference based on 500 samples (the experimental study shows that the algorithm has converged after 500 samples).

4. Examples

In this section, we analyze the real dataset that it contains information about crime in various cities in the United States. It is available at the University of Irvine machine learning repository (http://archive.ics.uci.edu/ml/datasets/communities+and+crime). The per capita violent crimes was calculated using population and the sum of crimevariables considered violent crimes in the United States: murder, rape, robbery, and assault. The observed individuals are communities. The dataset has 116 variables where the first four columns are response, name of community, code of county, code of community, the middle features are demographic information about each community such as population, age, race, income, and the final columns are region. According to the source: "the data combines socio-economic data from the 1990 US Census, law enforcement data from the 1990 US LEMAS survey, and crime data from the 1995 FBI UCR". This dataset has been applied to quantile regression [30]. The dataset is available at:

• train set: https://academics.hamilton.edu/mathematics/ckuruwit/Data/Crime/train.csv.
• test set: https://academics.hamilton.edu/mathematics/ckuruwit/Data/Crime/train.csv.

Our response variable of interest is the murder rate of each city (response variable $y_i\in [0,1]$). We choose this variable because it is the most dangerous one among the list of violent crimes. Identifying factors correlated with murder rate will be of immense importance to law enforcement agencies and the public at large.

To adapt the data to our model, we preprocess the data as follows:

• Delete columns from the data set that contain missing data.
• Delete the data when the response variable $y_i$ equals 0 because this is not an issue of interest to us.
• Transform $y_i$: $y_i^{\prime }=\log \frac{y_i}{1-y_i}$ and letting $y_i^{\prime }$ as the new response variable.
• Convert some qualitative variables into quantitative variables.
• Standardized covariates.

After the above data preprocessing, we get 1060 observation objects and 95 covariates in train set, and we get 122 observation objects and 95 covariates in test set. We estimate a quantile regression model between the response $y_i^{\prime }$ and the 96 predictors with intercept.

In this section, we only compare QRL, QRAL, BLQR, BALQR with our method VBSSLQR for real dataset, they are all quantile penalty regression obviously. Under different quantile $\tau \in (0.1,0.3,0.5,0.7,0.9)$, we compare the performance of different approachs. We counted the Root Mean Squared Error (RMSE) of each method under each quantile $\tau$ and the number of selected active variables to evaluate the performance of each approach in test set, where RMSE are evaluated by $\mathrm{RMSE}=\sqrt{\frac{1}{n}{\sum }_{i=1}^n{(y_i-{\widehat{y}}_{\tau i})}^2}$ with ${\widehat{y}}_{\tau i}$ is the fitted value of response $y_i$ under quantile $\tau$. Finally, the results are listed in the Table 5 below.

To visually show the results listed in Table 4, we will use best results are bolded under each quantile $\tau$. Obviously, our method performs better than other methods at all quantile and the active variables selected by our method are suitable, which means that our method is very competitive compared with other methods. BLQR and BALQR approachs can only identify the intercept, and they cannot identify the variables that really affect the response quantile. QRL and QRAL methods are prone to overfitting because they recognize too many active variables. Finally, the efficiency of our method in real data is also significantly higher than that of other approachs. So far, we believe that there are sufficient reasons to show that our method is enough competitive compared with other approachs.

Similar to [30], we list the active variables selected under each quantile in Table 6. So, Table 6 shows the variable selection of our proposed method in real data.

In the above table, common variables represents the five variables that appear most frequently. As evident in the Table 6, our method selected only a small number of predictors at each quantile level. The common 5 variables that the model picked are:

• PctVacantBoarded: percentage of households that are vacant and boarded up to prevent from vandalism.
• NumInShelters: number of shelters in the community.
• FemalePctDiv: percentage of females who are divorced.
• TotalPctDiv: percentage of peoples who are divorced.
• racePctWhite: percentage of person who are white race.

In order to get a better understanding of these five common variables we plot their correlation to the response.

Figure 1 depicts the relationship of five common variables and murder rate. Only the correlation between the NumInShelters and murder rate is not obvious, and the MalePctDivorce, RentLowQ, racePctWhite, PctVacantBoarded significantly affecting murder rate. The racePctWhite and murder rate are negatively correlated, while FemalePctDiv, TotalPctDiv, PctVacantBoarded and murder rate are positively correlated. The result was in accord with the practical situation and it can be seen that our results are basically consistent with that of [30], besides our method can more comprehensively select important variables under the same quantile. Therefore, percentage of females who are divorced, percentage of persons who are divorced, percentage of households that are vacant and boarded up to prevent from vandalism and percentage of person who are white race affect the murder rate.

5. Conclusions

In this paper, we propose variational Bayesian Spike-and-Slab Lasso quantile regression variable selection and estimation. This method applies Spike-and-Slab Lasso prior to each regression coefficient $\mathit{\beta }$, thus punishing the regression coefficient $\mathit{\beta }$. The spike prior distribution we choose is a small variance Laplace distribution, while the Slab prior distribution is a large variance Laplace distribution [26]. It is precisely because it punishes each regression coefficient $\mathit{\beta }$, giving it a powerful variable selection function, but it also brings a problem of low efficiency, especially when introducing Spike-and-Slab prior into quantile regression (Note, the algorithm efficiency of quantile regression is inherently lower than that of general regression approachs). In order to solve the problem of inefficiency (influenced by both quantile regression and Spike-and-Slab Lasso prior) and make the algorithm feasible, we introduce the variational Bayesian method to approximate the posterior distribution of each parameter. Simulation studies and real data analyses show that variational Bayesian Spike-and-Slab Lasso quantile regression performs well and has strong competitiveness compared with other approachs (Bayesian method or non-Bayesian method, quantile regression or non-quantile regression), especially in the case of high-dimensional data.