Sparse model: using spike and slab prior

Sparse model for expression and genotype

\[Y = \tilde{X}\beta + G\theta + \epsilon\] \(Y\): quantitative traits, \(N\) x 1 vector, \(N\) is number of individuals, the columns are centered.

\(\tilde{X}\): the genetic component of gene expression, \(N\) x \(J\) matrix, \(J\) is number of genes, the columns are centered.

\(G\): genotype (standardized). \(N\) x \(M\) matrix, \(M\) is number of SNPs.

The error term \(\epsilon\) has a normal distribution,

\[\epsilon|\tau \sim N(0, \tau^{-1})\]

We use the spike and slab prior for gene expression effect size \(\beta_j\) (effect size for gene \(j\)). \[\beta_j|\pi_\beta,\tau \sim (1- \pi_\beta)\delta_0 + \pi_\beta N(0, \sigma_\beta^2/\tau)\] Here \(\delta_0\) is point mass at 0 (\(\delta_0(\beta_j)=1\) if \(\beta_j=0\), otherwise \(\delta_0(\beta_j)=0\) ).

We use the spike and slab prior for SNP effect size \(\theta_m\) (effect size for SNP \(m\)). \[ \theta_m | \tau, \pi_\theta\sim (1- \pi_\theta)\delta_0 + \pi_\theta N(0, \sigma_\theta^2/\tau)\]

We use the following prior distributions for the hyperparameters \(\omega = (\pi_\beta,\pi_\theta, \tau, \sigma_\beta^2,\sigma_\theta^2)\): \[\tau \sim Gamma(\kappa_1, \kappa_2)\]

\[log(\pi_\theta) \sim U(log(1/M), log(1))\] \[log(\pi_\beta) \sim U(log(1/J), log(1))\] We define \(h_{SNP}\) and \(h_{expr}\) as follows: \[ h_{SNP} := \frac{\pi_\theta M\sigma_\theta^2}{\pi_\theta M\sigma_\theta^2 + \pi_\beta Jvar(\tilde{X})\sigma_\beta^2+ 1}\] \[ h_{expr} := \frac{\pi_\beta Jvar(\tilde{X})\sigma_\beta^2}{\pi_\theta M\sigma_\theta^2 + \pi_\beta Jvar(\tilde{X})\sigma_\beta^2+ 1}\]

Instead of assigning priors for \(\sigma_\beta^2\) and \(\sigma_\theta^2\), we assign uniform priors for \(h_{SNP}\) and \(h_{expr}\):

\[ h_{SNP} \sim U(0,1), h_{expr} \sim U(0,1) \]

Inference

We are interested in the hyperparameters \(\omega = (\pi_\beta, \pi_\theta, \tau, \sigma_\beta^2,\sigma_\theta^2)\). We introduce two vectors of binary indicators (\(\gamma_\beta\) = \(\{\gamma_{\beta1}, \gamma_{\beta2}, ..., \gamma_{\beta J}\} \in \{0,1\}^J\) and \(\gamma_\theta\) = \(\{\gamma_{\theta1}, \gamma_{\theta2}, ..., \gamma_{\theta M}\} \in \{0,1\}^M\)) that indicates whether the corresponding \(\beta_j\) or \(\theta_m\) is non-zero.

\[\gamma_{\beta j} \sim Bernoulli(\pi_\beta)\] \[\gamma_{\theta m} \sim Bernoulli(\pi_\theta)\] \[ \beta_j|\gamma_{\beta j}=1 \sim N(0, \sigma_\beta^2/\tau), \beta_j|\gamma_{\beta j}=0 \sim \delta_0\] \[ \theta_m|\gamma_{\theta m}=1 \sim N(0, \sigma_\theta^2/\tau), \theta_m|\gamma_{\theta m}=0 \sim \delta_0\] The posterior distrition is

\[P(\pi_\beta, \pi_\theta, \tau, h_{expr}, h_{SNP}, \gamma_\beta, \gamma_\theta|y) \propto P(y|\pi_\beta, \pi_\theta, \tau, h_{expr}, h_{SNP}, \gamma_\beta, \gamma_\theta) P(h_{expr})P(h_{SNP})P(\gamma_\beta|\pi_\beta)P( \gamma_\theta|\pi_\theta)P(\pi_\theta)P(\pi_\beta)\]

Session information

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Scientific Linux 7.4 (Nitrogen)

Matrix products: default
BLAS/LAPACK: /software/openblas-0.2.19-el7-x86_64/lib/libopenblas_haswellp-r0.2.19.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
 [1] workflowr_1.6.0 Rcpp_1.0.4.6    digest_0.6.18   later_0.7.5    
 [5] rprojroot_1.3-2 R6_2.3.0        backports_1.1.2 git2r_0.26.1   
 [9] magrittr_1.5    evaluate_0.12   stringi_1.3.1   fs_1.3.1       
[13] promises_1.0.1  whisker_0.3-2   rmarkdown_1.10  tools_3.5.1    
[17] stringr_1.4.0   glue_1.3.0      httpuv_1.4.5    yaml_2.2.0     
[21] compiler_3.5.1  htmltools_0.3.6 knitr_1.20