Bayesian variable selection and estimation in quantile regression using a quantile-specific prior

Abstract

Asymmetric Laplace (AL) specification has become one of the ideal statistical models for Bayesian quantile regression. In addition to fast convergence of Markov Chain Monte Carlo, AL specification guarantees posterior consistency under model misspecification. However, variable selection under such a specification is a daunting task because, realistically, prior specification of regression parameters should take the quantile levels into consideration. Quantile-specific g-prior has recently been developed for Bayesian variable selection in quantile regression, whereas it comes at a high price of the computational burden due to the intractability of the posterior distributions. In this paper, we develop a novel three-stage computational scheme for the foregoing quantile-specific g-prior, which starts with an expectation-maximization algorithm, followed by Gibbs sampler and ends with an importance re-weighting step that improves the accuracy of approximation. The performance of the proposed procedure is illustrated with simulations and a real-data application. Numerical results suggest that our procedure compares favorably with the Metropolis–Hastings algorithm.

Appendix: The EM algorithm

Consider the linear regression model in (1.1) with the ALD errors, i.e., $\epsilon \sim \mathrm {ALD}(0, \sigma , p)$. To implement an EM algorithm to estimate the maximum likelihood estimates (MLEs) of $\varvec{\beta }$ and $\sigma$, we may treat $\varvec{\nu }$ as the missing data. Based on the non-informative prior $\pi (\varvec{\beta }, \sigma ) \propto 1/\sigma$ and the likelihood function of the complete data $(\mathbf{{y}}, \varvec{\nu })$ in (2.1), the complete posterior distribution of $(\varvec{\beta }, \sigma )$ is given by

$$\begin{aligned}&f(\varvec{\beta }, \sigma \mid \varvec{\nu }, \mathbf{{y}}) \propto \sigma ^{-(3n + 2)/2} \Bigg (\prod _{i = 1}^{n} \nu _i^{-1/2}\Bigg ) \\&\quad \times \exp \Bigg \{-\frac{1}{4 \sigma } \bigg [\sum _{i = 1}^{n} \frac{(y_i - \mathbf{{x}}'_i \varvec{\beta }- \xi \nu _i)^2}{\nu _i} + 4p (1 - p) \sum _{i = 1}^{n} \nu _i \bigg ]\Bigg \}. \end{aligned}$$

The log-likelihood of the complete data can be written as

$$\begin{aligned}&\mathrm {L} (\varvec{\beta }, \sigma \mid \varvec{\nu }, \mathbf{{y}}) \propto -\frac{3n + 2}{2} \mathrm {log} \sigma - \frac{1}{2} \sum _{i = 1}^{n} \mathrm {log} \nu _i \\&\quad - \frac{1}{4 \sigma } \sum _{i = 1}^{n} \frac{(y_i - \mathbf{{x}}'_i \varvec{\beta }- \xi \nu _i)^2}{\nu _i} - \zeta \sum _{i = 1}^{n} \nu _i. \end{aligned}$$

By following the work of Tian et al. (2014), we let $\varvec{\beta }^{(0)}$ and $\sigma ^{(0)}$ be the initial values taken from the OLS estimates. We then obtain the MLEs of the model parameters based on the EM algorithm via the following two steps:

• 1. E step: Given the $(t - 1)$-th iteration values $\varvec{\beta }^{(t - 1)}$ and $\sigma ^{(t - 1)}$, we obtain the mathematical expectation of the complete data, also known as the Q function of the (t)-th iteration, given by

  $$\begin{aligned} Q(\varvec{\beta }, \sigma \mid \varvec{\nu }, \mathbf{{y}})&= \mathrm {E} \Big (\mathrm {L} (\varvec{\beta }, \sigma \mid \mathbf{{y}}, \varvec{\nu }) \mid \mathbf{{y}}, \varvec{\beta }^{(t - 1)}, \sigma ^{(t - 1)}\Big ) \\&\propto -\frac{3n + 2}{2} \mathrm {log} \sigma - \frac{1}{2} \sum _{i = 1}^{n} \mathrm {E} \Big (\mathrm {log} \nu _i \mid \mathbf{{y}}, \varvec{\beta }^{(t - 1)}, \sigma ^{(t - 1)}\Big ) \\&\quad - \frac{1}{4 \sigma } \sum _{i = 1}^{n} (y_i - \mathbf{{x}}'_i \varvec{\beta })^2 \mathrm {E} \Big (\nu _i^{-1} \mid \mathbf{{y}}, \varvec{\beta }^{(t - 1)}, \sigma ^{(t - 1)}\Big ) \\&\quad + \frac{\xi }{2 \sigma } \sum _{i = 1}^{n} (y_i - \mathbf{{x}}'_i \varvec{\beta }) \\&\quad - \bigg (\frac{\xi ^2}{4 \sigma } + \zeta \bigg ) \sum _{i = 1}^{n} \mathrm {E} \Big (\nu _i \mid \mathbf{{y}}, \varvec{\beta }^{(t - 1)}, \sigma ^{(t - 1)}\Big )\\&\propto -\frac{3n + 2}{2} \mathrm {log} \sigma + \frac{\xi }{2 \sigma } \sum _{i = 1}^{n} (y_i - \mathbf{{x}}'_i \varvec{\beta }) - \frac{1}{2} \sum _{i = 1}^{n} \Delta 1_i \\&\quad - \frac{1}{4 \sigma } \sum _{i = 1}^{n} \Delta 2_i \\&\quad - \frac{1}{4 \sigma } \sum _{i = 1}^{n} (y_i - \mathbf{{x}}'_i \varvec{\beta })^2 \Delta 3_i, \end{aligned}$$

  where

  $$\begin{aligned} \Delta 1_i&= \mathrm {E} \Big (\mathrm {log} \nu _i \mid \varvec{\beta }^{(t - 1)}, \sigma ^{(t - 1)}, \mathbf{{y}}\Big ), \\ \Delta 2_i&= \mathrm {E} \Big (\nu _i \mid \varvec{\beta }^{(t - 1)}, \sigma ^{(t - 1)}, \mathbf{{y}}\Big ), \\ \Delta 3_i&= \mathrm {E} \Big (\nu _i^{-1} \mid \varvec{\beta }^{(t - 1)}, \sigma ^{(t - 1)}, \mathbf{{y}}\Big ). \end{aligned}$$

  To compute $\Delta 1_i$, $\Delta 2_i$, and $\Delta 3_i$ for $i = 1, \cdots , n$, we first obtain the conditional pdf $f(\varvec{\nu }\mid \varvec{\beta }, \sigma , \mathbf{{y}})$. We observe that

  $$\begin{aligned} f(\varvec{\nu }\mid \mathbf{{y}})&\propto \Bigg (\prod _{i = 1}^{n} \nu _i^{-1/2}\Bigg ) \exp \Bigg \{-\frac{1}{4 \sigma } \bigg [\sum _{i = 1}^{n} \frac{(y_i - \mathbf{{x}}'_i \varvec{\beta }- \xi \nu _i)^2}{\nu _i} \\&\quad + 4p (1 - p) \sum _{i = 1}^{n} \nu _i \bigg ]\Bigg \} \\&\propto \Bigg (\prod _{i = 1}^{n} \nu _i^{-1/2}\Bigg ) \exp \Bigg \{-\frac{1}{2} \sum _{i = 1}^{n} \bigg [\frac{(y_i - \mathbf{{x}}'_i \varvec{\beta })^2}{2 \sigma } \nu _i^{-1} + \frac{1}{2 \sigma } \nu _i \bigg ]\Bigg \}. \end{aligned}$$

  Since $\nu _i$’s are independent of each other for $i = 1, \cdots , n$, the marginal posterior distribution of $\nu _i$ given $y_i$ is given by

  $$\begin{aligned} f(\nu _i \mid y_i) \propto \nu _i^{-1/2} \exp \Bigg \{-\frac{1}{2} \bigg [\frac{(y_i - \mathbf{{x}}'_i \varvec{\beta })^2}{2 \sigma } \nu _i^{-1} + \frac{1}{2 \sigma } \nu _i \bigg ]\Bigg \}, \end{aligned}$$

  which follows a GIG distribution, denoted by $\mathrm {GIG} \Bigl (1/2, (y_i - \mathbf{{x}}'_i \varvec{\beta })^2/(2 \sigma ), 1/(2 \sigma ) \Bigr )$. Here, if $x \sim \mathrm {GIG}(\lambda , \chi , \psi )$, then the pdf of x is given by

  $$\begin{aligned} f(x) = \frac{\chi ^{-\lambda }(\sqrt{\chi \psi })^\lambda }{2 K_\lambda (\sqrt{\chi \psi })} x^{\lambda - 1} \exp \Bigg \{-\frac{1}{2} \Big (\chi x^{-1} + \psi x\Big ) \Bigg \}, \quad x > 0, \end{aligned}$$

  where $K_\lambda (\cdot )$ denotes a modified Bessel function of the third kind with an index $\lambda$ and the parameters such that $\chi > 0$, and $\psi \ge 0$ if $\lambda < 0$. When $\chi > 0$, $\psi > 0$, and $\alpha \in R$, we have

  $$\begin{aligned} \mathrm {E}(X^\alpha ) = \bigg (\frac{\chi }{\psi }\bigg )^{\alpha /2} \frac{K_{\lambda + \alpha }(\sqrt{\chi \psi })}{K_{\lambda }(\sqrt{\chi \psi })} \quad \mathrm {and} \quad \mathrm {E}(\mathrm {log} X) = \frac{d\mathrm {E}(X^\alpha )}{d\alpha } \Bigr |_{\alpha = 0}. \end{aligned}$$

  Therefore, we obtain the following results given by

  $$\begin{aligned} \Delta 1_i&= \mathrm {E} \Big (\mathrm {log} \nu _i \mid \varvec{\beta }^{(t - 1)}, \sigma ^{(t - 1)} , \mathbf{{y}}\Big ) = \frac{d\mathrm {E} \big (\nu _i^\alpha \mid \mathbf{{y}}, \varvec{\beta }^{(t - 1)}, \sigma ^{(t - 1)} \big )}{d\alpha } \Bigr |_{\alpha = 0}, \\ \Delta 2_i&= \mathrm {E} \Big (\nu _i \mid \varvec{\beta }^{(t - 1)}, \sigma ^{(t - 1)}, \mathbf{{y}}\Big ) = \big \vert y_i - \mathbf{{x}}'_i \varvec{\beta }^{(t - 1)} \big \vert + 2 \sigma ^{(t - 1)}, \\ \Delta 3_i&= \mathrm {E} \Big (\nu _i^{-1} \mid \varvec{\beta }^{(t - 1)}, \sigma ^{(t - 1)}, \mathbf{{y}}\Big ) = \frac{1}{\big \vert y_i - \mathbf{{x}}'_i \varvec{\beta }^{(t - 1)} \big \vert }. \end{aligned}$$

• 2. M step: We update the $\varvec{\beta }^{(t)}$ and $\sigma ^{(t)}$ values by maximizing the Q function $Q(\varvec{\beta }, \sigma \mid \varvec{\beta }^{(t - 1)}, \sigma ^{(t - 1)})$. We first take the derivative with respect to $\varvec{\beta }$ to obtain the MLE of $\varvec{\beta }$ given by

  $$\begin{aligned} \frac{\partial Q}{\partial \varvec{\beta }} = -\frac{\xi }{2 \sigma } \sum _{i = 1}^{n} \mathbf{{x}}_i + \frac{1}{2 \sigma } \sum _{i = 1}^{n} \Delta 3_i (y_i - \mathbf{{x}}'_i \varvec{\beta }) \mathbf{{x}}_i {\mathop {=}\limits ^{set}} 0, \end{aligned}$$

  which indicates that

  $$\begin{aligned} \sum _{i = 1}^{n} \Delta 3_i \mathbf{{x}}_i \mathbf{{x}}'_i \varvec{\beta }= \sum _{i = 1}^{n} \Delta 3_i \mathbf{{x}}_i y_i - \xi \sum _{i = 1}^{n} \mathbf{{x}}_i. \end{aligned}$$

  Letting $\varvec{\nabla 3} = (1/\Delta 3_1, \cdots , 1/\Delta 3_n)',$ $\mathbf{{W}}^{(t - 1)} = \mathrm {diag} (\Delta 3_1, \cdots , \Delta 3_n)$, and $\varvec{1}_n = (1, \cdots , 1)'$, we can rewrite (2) in a matrix form as

  $$\begin{aligned} \mathbf{{X}}' \mathbf{{W}}^{(t - 1)} \mathbf{{X}}\varvec{\beta }= \mathbf{{X}}' \mathbf{{W}}^{(t - 1)} \mathbf{{y}}- \xi \mathbf{{X}}' \varvec{1}_n. \end{aligned}$$

  Thus, we obtain an estimate of $\varvec{\beta }$ given by

  $$\begin{aligned} {\hat{\varvec{\beta }}}^{(t)} = \Big ( \mathbf{{X}}' \mathbf{{W}}^{(t - 1)} \mathbf{{X}}\Big )^{-1} \mathbf{{X}}' \mathbf{{W}}^{(t - 1)} (\mathbf{{y}}- \xi \varvec{\nabla 3}). \end{aligned}$$

  (6.1)

  By taking the derivative with respect to $\sigma$ to obtain the MLE expression of $\sigma$, we have

  $$\begin{aligned} \frac{\partial Q}{\partial \sigma }&= -\frac{3n + 2}{2 \sigma } - \frac{\xi }{2 \sigma ^2} \sum _{i = 1}^{n} (y_i - \mathbf{{x}}'_i \varvec{\beta }) + \frac{1}{4 \sigma ^2} \sum _{i = 1}^{n} \Delta 2_i \\&\quad + \frac{1}{4 \sigma ^2} \sum _{i = 1}^{n} \Delta 3_i (y_i - \mathbf{{x}}'_i \varvec{\beta })^2 {\mathop {=}\limits ^{set}} 0. \end{aligned}$$

  Thus, we obtain an estimate of $\sigma$ give by

  $$\begin{aligned} {\hat{\sigma }}^{(t)}= & {} \frac{1}{2(3n + 2)}\nonumber \\&\left[ \sum _{i = 1}^{n} \Delta 2_i + \sum _{i = 1}^{n} \Delta 3_i (y_i - \mathbf{{x}}'_i \varvec{\beta }^{(t)})^2 - 2 \xi \sum _{i = 1}^{n} (y_i - \mathbf{{x}}'_i \varvec{\beta }^{(t)})\right] . \end{aligned}$$

  (6.2)

We repeatedly run the EM algorithm until a suitable convergence rule has been satisfied; i.e. $\Vert \varvec{\beta }^{(t)} - \varvec{\beta }^{(t-1)}\Vert$ and $\Vert \sigma ^{(t)} - \sigma ^{(t-1)}\Vert$ are sufficiently small, for example. Then we set the final iterative values of the EM algorithm as the posterior modes, denoted by ${\tilde{\varvec{\beta }}}$ and ${{\tilde{\sigma }}}$, for $\varvec{\beta }$ and $\sigma$, respectively.