Bayesian variable selection for mixed effects model with shrinkage prior

Abstract

Recently, many shrinkage priors have been proposed and studied in linear models to address massive regression problems. However, shrinkage priors are rarely used in mixed effects models. In this article, we address the problem of joint selection of both fixed effects and random effects with the use of several shrinkage priors in linear mixed models. The idea is to shrink small coefficients to zero while minimally shrink large coefficients due to the heavy tails. The shrinkage priors can be obtained via a scale mixture of normal distributions to facilitate computation. We use a stochastic search Gibbs sampler to implement a fully Bayesian approach for variable selection. The approach is illustrated using simulated data and a real example.

Appendix

We outline the detailed Gibbs sampler for posterior sampling as follows.

• 1. Given the data and the current values of $g, \varvec{\varLambda },\varvec{\varGamma }, \varvec{\xi }, \sigma ^2, J$, the posterior of $\varvec{\beta }_J$ is sampled from a multivariate normal $N(\varvec{\mu }_{\varvec{\beta }_J}, \varvec{\varSigma }_{\varvec{\beta }_J})$:

  $$\begin{aligned} \varvec{\varSigma }_{\varvec{\beta }_J}^{-1}= & {} {\varvec{X}}^{{\varvec{J}}'} {\varvec{X}}^{{\varvec{J}}} /\sigma ^2 + g/\sigma ^2 {\varvec{\zeta }}_{\varvec{J}}^{{\varvec{-1}}}, \quad {\varvec{\zeta }}_{\varvec{J}}=\text{ Diag }(\kappa _{k:J_k=1}, k=1, \ldots , l ),\\ \varvec{\mu }_{\varvec{\beta }_J}= & {} \varvec{\varSigma }_{\varvec{\beta }_J} \left\{ \sum _i \sum _j {\varvec{Z}}_{{\varvec{ij}}}(y_{ij}- {\varvec{X}}_{{\varvec{ij}}}' \varvec{\varLambda }\varvec{\varGamma }\varvec{\xi }_i )/\sigma ^2 \right\} . \end{aligned}$$

• 2. Given the data and the current values of $\varvec{\varLambda },\varvec{\varGamma }, \varvec{\xi }, g, \sigma ^2, \varvec{\beta }$, the full conditional posterior of $J_k$ for $k=1, \ldots , l$ is calculated by integrating out $\varvec{\beta }$ and $\sigma ^2$:

  $$\begin{aligned} p(J_k=1|\text{ else })= & {} \frac{1.0}{1.0 + \varpi _k}, \quad \varpi _k= \frac{S_1(J_k=0)}{S_1(J_k=1)} \frac{S_2(J_k=0)}{S_2(J_k=1)}\frac{1-p_{k0}}{p_{k0}} \sqrt{\frac{\kappa _k}{g}},\\ S_1(J)= & {} \left\{ \sum _i \sum _j (y_{ij}- {{\varvec{X}}_{{\varvec{ij}}}'} \varvec{\varLambda }\varvec{\varGamma }\varvec{\xi }_i )^2 - \sigma ^2 \varvec{\mu }_{\varvec{\beta }_J}' \varvec{\varSigma }_{\varvec{\beta }_J}^{-1}\varvec{\mu }_{\varvec{\beta }_J} \right\} ^{-\sum _i n_i/2}, \\ S_2(J)= & {} \left| \sum _i\sum _j {\varvec{X}}_{{\varvec{ij}}} {\varvec{X}}_{{\varvec{ij}}} + g {\varvec{\zeta }}_{\varvec{J}}^{{\varvec{-1}}} \right| ^{-1/2}. \end{aligned}$$

• 3. Given the data and updated values of $J, \varvec{\beta }, \sigma ^2$, the posterior of g is sampled from $G(a^*, b^*)$, where

  $$\begin{aligned} a^*= & {} \left( 1.0 + \sum _k J_k \right) /2.0, \\ b^*= & {} N/2 + \varvec{\beta }^{J'} {\varvec{\zeta }}_{\varvec{J}}^{{\varvec{-1}}}\varvec{\beta }^{J}/(2\sigma ^2). \end{aligned}$$

• 4. Given the data and updated values of $J, \varvec{\beta }, \sigma ^2, \varvec{\varLambda }, \varvec{\varGamma }, \varvec{\xi }$, the posterior of $\sigma ^2$ is sampled from $IG(h_1, h_2)$, where

  $$\begin{aligned} h_1= & {} \frac{ N+\sum _k J_k}{2}, \\ h_2= & {} \sum _i\sum \tau _{ij}^2 + \varvec{\beta }^{J'} g\zeta _J^{-1}\varvec{\beta }^{J}/2,\\ \tau _{ij}= & {} y_{ij}- {\varvec{X}}_{{\varvec{ij}}}' \varvec{\beta }- {\varvec{Z}}_{{\varvec{ij}}}' \varvec{\varLambda }\varvec{\varGamma }\varvec{\xi }_i. \end{aligned}$$

• 5. Reparameterize model (2) as $y_{ij}^*=y_{ij}-\sum _{k=1}^q z_{ijk}\lambda _k \varvec{\xi }_{ik}$ and $y_{ij}^*= {\varvec{X}}_{{\varvec{ij}}}^{{\varvec{J}}'} \varvec{\beta }^J + {\varvec{o}}_{{\varvec{ij}}}' {\varvec{r}}$, with the kth element ${\varvec{r}}_k=\varvec{\varGamma }(fm), f=2, \ldots , q, m=1, \ldots f-1,$ and ${\varvec{o}}_{{\varvec{ij}}}'= z_{ijf}\lambda _f\xi _{im}$. With prior ${\varvec{r}} \sim N(\varvec{\mu }_{r0}, \varvec{\varSigma }_{r0})$, the posterior is also multivariate normal $N(\varvec{\mu }_{r}^*, \varvec{\varSigma }_{r}^*)1( {\varvec{r}} \in R_r)$ where ${\varvec{r}} \in R_r$ constrains elements of ${\varvec{r}}$ to be zero when the corresponding random effects are zeroes.

  $$\begin{aligned} \varvec{\varSigma }_{r}^{*-1}= & {} \varvec{\varSigma }_{r0}^{-1} + \sum _i \sum _j {\varvec{o}}_{{\varvec{ij}}} {\varvec{o}}_{{\varvec{ij}}}'/\sigma ^2, \\ \varvec{\mu }_{r}^*= & {} \varvec{\varSigma }_r^* \left\{ \varvec{\varSigma }_{r0}^{-1} \varvec{\mu }_{r0} + \sum _i \sum _j {\varvec{o}}_{{\varvec{ij}}}' ( y_{ij}^* -{\varvec{X}}_{{\varvec{ij}}}^{{\varvec{J}}'} \varvec{\beta }^J)/\sigma ^2 \right\} . \end{aligned}$$

• 6. Reparameterize $T_{ij}=y_{ij} - {\varvec{X}}_{{\varvec{ij}}}^{{\varvec{J}}'}\varvec{\beta }^{J}-\sum _{k \ne l}t_{ijk}\lambda _k$, where $t_{ijk}=\sum _{m=1}^k z_{ijk}\gamma _{km}\varvec{\xi }_{im}$. For prior $\lambda _k \sim ZIN_+(p_{k0}, \mu _0, \sigma _0^2)$, the full conditional posterior is $ZIN_+(p_k^*, \mu ^*,\sigma ^{*2})$:

  $$\begin{aligned} \sigma ^{*-2}= & {} \sum _i \sum _j t_{ijl}^2/\sigma ^2 + 1.0,\\ \mu ^*= & {} \sigma ^{*2}\left( \sum _i \sum _j t_{ijl}T_{ij}/\sigma ^2\right) ,\\ p_k^*= & {} \frac{p_{k0}}{p_{k0}+(1-p_{k0}) \frac{N(0|\mu _0, \sigma _0^2)}{N(0|\mu ^*, \sigma ^{*2})} \frac{1-\varPhi (0, \mu ^*, \sigma ^{*2}) }{1-\varPhi (0, \mu _0, \sigma _0^2)} }. \end{aligned}$$

  where $\varPhi (.)$ is the CDF of standard normal distribution, and $N(a|\mu , \sigma ^2)$ is the pdf of normal distribution with mean $\mu$, variance $\sigma ^2$ at a.

• 7. Sample the posterior $\varvec{\xi }_i$ from a multivariate normal distribution with parameters

  $$\begin{aligned} \varvec{\varSigma }^{*-1}= & {} \varvec{\varSigma }_{\varvec{\xi }}^{-1} + \sum _{i=1}^{n_i} \varvec{\varGamma }'\varvec{\varLambda }' {\varvec{Z}}_{{\varvec{ij}}} {\varvec{Z}}_{{\varvec{ij}}}'\varvec{\varLambda }\varvec{\varGamma }\sigma ^2 , \quad \varvec{\varSigma }_{\varvec{\xi }}=\text{ Diag }( \tau _k), \ k=1, \ldots , q\\ \varvec{\mu }^*= & {} \varvec{\varSigma }^* \left\{ \sum _{j=1}^{n_i} \varvec{\varGamma }'\varvec{\varLambda }' {\varvec{Z}}_{{\varvec{ij}}}(y_{ij} - {\varvec{X}}_{{\varvec{ij}}}^{J'} \varvec{\beta }^{J} )/\sigma ^2 \right\} \end{aligned}$$

• 8. Given the updated parameter $\lambda _k$, update $\phi _k^2$ from the inverse gamma distribution $IG(h_3, h_4)$, where

  $$\begin{aligned} h_3= & {} 1.0\\ h_4= & {} (1.0 + \lambda _k)/2 \end{aligned}$$

• 9. Given the updated parameters $\varphi$, $\sigma$, $\beta$ and g, update $\kappa$:

  $$\begin{aligned} \kappa _k|\text{ else } \sim \text{ Inv-Gau } \left( \left| \frac{\varphi _k \sigma }{\beta _k \sqrt{g}} \right| , \varphi _k^2 \right) , \end{aligned}$$

  where $\text{ Inv-Gau } \text{(a, } \text{ b) }$ is the inverse Gaussian distribution with location and scale parameters a and b.

• 10. Given the updated parameter $\sigma$, g and $\beta$, update $\varphi$ from

  $$\begin{aligned} \varphi _k \sim Gamma(2.0, |\beta _k \sqrt{g}/\sigma | + 1.0 ) \end{aligned}$$

• 11. Propose the candidate density $\exp (1)$ for $\tau _k$, the updated value $\tau _k^*$ is accepted with probability $\min { \left\{ 1, \frac{f(\tau _k^*)\exp ( -\tau _k^0)}{f(\tau _k^0)\exp ( -\tau _k^*)}\right\} }$, where $f(\tau _k)=\tau _k^{-n/2} \exp { (-\sum _i \varvec{\xi }_{ik}^2/(2\tau _k)} -m_k^2\tau _k/2 )$.

• 12. Propose the candidate density $\exp (1)$ for $m_k$, the updated value $m_k^*$ is accepted with probability $\min { \left\{ 1, \frac{g(m_k^*)\exp {( -m_k^0)}}{g(m_k^0)\exp ( -m_k^*)}\right\} }$, where $g(m_k)=m_k^2 \exp { (-m_k^2\tau _k/2 -m_k ) }$.