Sparse model

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{M\sigma_\theta^2}{M\sigma_\theta^2 + Jvar(\tilde{X})\sigma_\beta^2+ 1}$$ $$h_{expr} := \frac{Jvar(\tilde{X})\sigma_\beta^2}{M\sigma_\theta^2 + 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.0      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