Last updated: 2018-08-20
workflowr checks: (Click a bullet for more information) ✔ R Markdown file: up-to-date
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
✔ Environment: empty
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
✔ Seed:
set.seed(20180529)
The command set.seed(20180529) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.
✔ Session information: recorded
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
✔ Repository version: 0ea299b
wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
Ignored files:
Ignored: .Rhistory
Ignored: .Rproj.user/
Ignored: analysis/.Rhistory
Ignored: docs/.DS_Store
Untracked files:
Untracked: analysis/susie_summarySTAT_Mix.Rmd
Untracked: docs/figure/
Unstaged changes:
Modified: analysis/BayesianLinearRegression.Rmd
Modified: analysis/_site.yml
Modified: analysis/susie_summarySTAT.Rmd
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| Rmd | 0ea299b | zouyuxin | 2018-08-20 | wflow_publish(“analysis/susie.Rmd”) |
| html | b2d26ee | zouyuxin | 2018-07-17 | Build site. |
| Rmd | ed14aba | zouyuxin | 2018-07-17 | wflow_publish(c(“analysis/index.Rmd”, “analysis/susie.Rmd”, |
| html | 295dc11 | zouyuxin | 2018-06-29 | Build site. |
| Rmd | 742ff45 | zouyuxin | 2018-06-29 | wflow_publish(“analysis/susie.Rmd”) |
| html | 0dcada6 | zouyuxin | 2018-06-12 | Build site. |
| Rmd | c515e90 | zouyuxin | 2018-06-12 | wflow_publish(“analysis/susie.Rmd”) |
| html | db9f1ce | zouyuxin | 2018-06-07 | Build site. |
| Rmd | 6ec2be3 | zouyuxin | 2018-06-07 | wflow_publish(“analysis/susie.Rmd”) |
| html | a2803fc | zouyuxin | 2018-06-01 | Build site. |
| Rmd | 2c16cd3 | zouyuxin | 2018-06-01 | wflow_publish(“analysis/susie.Rmd”) |
| html | 44284bd | zouyuxin | 2018-05-29 | Build site. |
| Rmd | 751426d | zouyuxin | 2018-05-29 | wflow_publish(c(“analysis/BayesianLinearRegression.Rmd”, |
\[ y_{i} = \sum_{l=1}^{L} \sum_{j=1}^{p} X_{ij} \gamma_{lj} \beta_{lj} + \epsilon_{i} \quad \epsilon_{i} \sim N(0, \sigma^2) \]
\[ (\gamma_{l1}, \cdots, \gamma_{lp}) \sim Multinomial(1, \pi) \quad \pi = (\frac{1}{p}, \cdots, \frac{1}{p}) \] \[ \beta_{lj} \sim N(0,\sigma^2\sigma_{\beta}^2) \]
\(B = \left[\mathbf{\beta}_{1} \ \cdots \ \mathbf{\beta}_{L} \right]^{T}\); \(\Gamma = \left[\mathbf{\gamma}_{1} \ \cdots \ \mathbf{\gamma}_{L} \right]^{T}\)
There are L non-zero effects. For each of \(l = 1, \cdots, L\) we assume that exactly 1 of the p variables has an effect, as indicated by \(\gamma_{lj}\). \(\pi\) is a p vector of prior probabilities summing to 1.
We want to seek a variational approximation based on \[ q(B, \Gamma) = \prod_{l} q_{l}(\mathbf{\beta}_{l}, \mathbf{\gamma}_{l}) = \prod_{l} q_{l}(\mathbf{\beta}_{l}|\mathbf{\gamma}_{l})q_{l}(\mathbf{\gamma}_{l}) = \prod_{l=1}^{L} q(\mathbf{\gamma}_{l})\prod_{j=1}^{p}q(\beta_{lj}|\mathbf{\gamma}_{l}) \]
\[ q_{l}(\mathbf{\gamma}_{l}) = Multinomial(1,\mathbf{\alpha}_{l}) \quad \mathbf{\alpha}_{l} = (\alpha_{l1}, \cdots,\alpha_{lp}) \quad \sum_{j=1}^{p} \alpha_{lj} = 1 \]
\[ p(\mathbf{y},B, \Gamma|X, \sigma^2, \sigma_{\beta}^2, \pi) = \left( \frac{1}{\sqrt{2\pi \sigma^2}} \right)^{n} \exp \left( -\frac{\|\mathbf{y} - X\mathbf{\beta}_{1} - \cdots - X\mathbf{\beta}_{L}\|^2}{2\sigma^2} \right) \prod_{l=1}^{L}\frac{1}{p}N(\beta_{lj}; 0, \sigma^2\sigma_{\beta}^2) \]
The result \(q\) is conditional distribution for \(B, \Gamma\) (Details in VB)
\[ \beta_{lj}|\gamma_{lj}=1 \sim_{q_{l}} N(\mu_{lj}, s_{lj}^{2}) \quad q(\gamma_{l}) = \alpha_{lj} \]
\[ s_{lj}^{2} = \frac{\sigma^2\sigma_{\beta}^2}{(X^{T}X)_{jj}\sigma_{\beta}^2 + 1} = \left( \frac{(X^{T}X)_{jj}}{\sigma^2} + \frac{1}{\sigma^2\sigma_{\beta}^2} \right)^{-1} \] \[ \mu_{lj} = \frac{(X^{T}\mathbf{y}_{res,-l})_{j}}{(X^{T}X)_{jj} + 1/\sigma_{\beta}^2} = \frac{s_{lj}^{2}}{\sigma^2} (X^{T}\mathbf{y}_{res,-l})_{j} = \frac{s_{lj}^{2}}{\sigma^2} [(X^{T}\mathbf{y})_{j} - \sum_{k\neq j}(X^{T}X)_{kj}\alpha_{lk}\mu_{lk}] \] \[ \alpha_{lj} \propto \frac{1}{p} N(\frac{(X^{T}\mathbf{y}_{res,-l})_{j}}{(X^{T}X)_{jj}}; 0, \frac{\sigma^2}{(X^{T}X)_{jj}} + \sigma^2\sigma_{\beta}^2 ) \] \[ \frac{\alpha_{lj}}{1-\alpha_{lj}} = \frac{1/p}{1-1/p} \frac{s_{lj}}{\sigma_{\beta}\sigma} \exp(\frac{\mu_{lj}^{2}}{2s_{lj}^{2}}) = \frac{1/p}{1-1/p} BF \]
\(\mathbf{r}_{l} = \mathbb{E}(\mathbf{\beta}_{l}) = \alpha_{l.}*\mu_{l.}\) \(\quad\) \(\mathbf{r} = \mathbb{E}(\mathbf{\beta}) = \sum_{l=1}^{L}\mathbf{r}_{l}\)
\[\begin{align*} \mathbb{E}_{q}\left[ \log N_{n}(\mathbf{y} - X\mathbf{\beta}; 0, \sigma^2I) \right] &= -\frac{n}{2}\log 2\pi\sigma^2 - \frac{\mathbb{E}\|\mathbf{y} - X\mathbf{\beta}\|^2}{2\sigma^2} \\ &= -\frac{n}{2}\log 2\pi\sigma^2 - \frac{1}{2\sigma^2}\left[\mathbf{y}^{T}\mathbf{y} - 2 \mathbf{y}^{T}X\mathbf{r} + \|X\mathbf{r}\|^2 - \sum_{l=1}^{L}\|X\mathbf{r}_{l}\|^{2} + \sum_{l}\sum_{j=1}^{p}(X^{T}X)_{jj}\mathbb{E}(\beta_{lj}^2)\right] \end{align*}\]\[ \mathbf{y} = X_{n\times p} \mathbf{\beta} + \epsilon \quad \epsilon \sim N_{n}(0, \sigma^2I) \] For non-zero effect \(\beta_{j}\) \[ \beta_{j}|\gamma_{j} = 1 \sim N(0, \sigma^2\sigma_{\beta}^2) \quad \mathbf{\gamma}\sim Multivariate(1,\pi) \quad \pi = (1/p, \cdots, 1/p) \]
\[ q(\mathbf{\beta}, \mathbf{\gamma}) = q(\mathbf{\beta}| \mathbf{\gamma})q(\mathbf{\gamma}) = q(\mathbf{\gamma}) \prod_{j=1}^{p} q_{j}(\beta_{j}| \mathbf{\gamma}) \]
\[ p(\mathbf{\beta}, \mathbf{\gamma}) = p(\mathbf{\beta}|\mathbf{\gamma})p(\mathbf{\gamma}) = p(\mathbf{\gamma})\prod_{j=1}^{p}p(\beta_{j}|\mathbf{\gamma}) = \frac{1}{p}\prod_{j=1}^{p}N(\beta_{j};0,\sigma^2\sigma_{\beta}^2) \]
The likelihood is \[ p(\mathbf{y},\mathbf{\beta}, \mathbf{\gamma}|X, \sigma^2, \sigma_{\beta}^2, \pi) = \left( \frac{1}{\sqrt{2\pi \sigma^2}} \right)^{n} \exp \left( -\frac{\|\mathbf{y} - X\mathbf{\beta}\|^2}{2\sigma^2} \right) N(\beta_{j}; 0, \sigma^2\sigma_{\beta}^2)\frac{1}{p} \]
\[\begin{align*} p(\mathbf{y}|X, \sigma^2, \sigma_{\beta}^2, \pi) &= \int p(\mathbf{y},\mathbf{\beta}, \mathbf{\gamma}|X, \sigma^2, \sigma_{\beta}^2, \pi) d \mathbf{\beta} d\mathbf{\gamma} \\ &= \sum_{j=1}^{p} \frac{1}{p} \int \left( \frac{1}{\sqrt{2\pi \sigma^2}} \right)^{n} \exp \left( -\frac{\|\mathbf{y} - X\mathbf{\beta}\|^2}{2\sigma^2} \right) N(\beta_{j}; 0, \sigma^2\sigma_{\beta}^2) d\beta_{j} \\ &= \sum_{j=1}^{p} \frac{1}{p} \frac{\int \left( \frac{1}{\sqrt{2\pi \sigma^2}} \right)^{n} \exp \left( -\frac{\|\mathbf{y} - X\mathbf{\beta}\|^2}{2\sigma^2} \right) N(\beta_{j}; 0, \sigma^2\sigma_{\beta}^2) d\beta_{j}}{N(\mathbf{y}; 0, \sigma^2)}N(\mathbf{y}; 0, \sigma^2) \\ &= \sum_{j=1}^{p} \frac{1}{p} \frac{\int N(\hat{\beta}_{j}; \beta_{j}, \frac{\sigma^2}{\mathbf{x}_{j}^{T}\mathbf{x}_{j}}) N(\beta_{j}; 0, \sigma^2\sigma_{\beta}^2) d\beta_{j}}{N(\hat{\beta}_{j}; 0, \frac{\sigma^2}{\mathbf{x}_{j}^{T}\mathbf{x}_{j}})}N(\mathbf{y}; 0, \sigma^2) \\ &= \sum_{j=1}^{p} \frac{1}{p} \frac{N(\hat{\beta}_{j}; 0, \frac{\sigma^2}{\mathbf{x}_{j}^{T}\mathbf{x}_{j}} + \sigma^2\sigma_{\beta}^2)}{N(\hat{\beta}_{j}; 0, \frac{\sigma^2}{\mathbf{x}_{j}^{T}\mathbf{x}_{j}})}N(\mathbf{y}; 0, \sigma^2) \\ &= \sum_{j=1}^{p} \frac{1}{p} BF_{j} N(\mathbf{y}; 0, \sigma^2) \\ \log p(\mathbf{y}|X, \sigma^2, \sigma_{\beta}^2, \pi) &= \mathcal{L}(q) + D_{KL}(q||p(\beta, \gamma|\mathbf{y}, X, \sigma^2, \sigma_{\beta}^2)) \end{align*}\]For the single effect regression, we compute the exact posterior. At the optimum, \(D_{KL}(q||p(\beta, \gamma|\mathbf{y}, X, \sigma^2, \sigma_{\beta}^2)) = 0\), the lower bound is equal to the log-likelihood.
Let’s compute the KL divergence for posterior and the prior
\[ D_{KL}(q||p) = \int q\log q - \int q \log p \]
\[ q(\mathbf{\beta}, \mathbf{\gamma}) = q(\mathbf{\gamma})\prod_{j=1}^{p}q(\beta_{j}|\mathbf{\gamma}) = \prod_{j=1}^{p} \alpha_{j}^{\gamma_{j}} \prod_{j=1}^{p}N(\beta_{j};\mu_{j},s_{j}^{2}) = \alpha_{j} \prod_{j=1}^{p}N(\beta_{j};\mu_{j},s_{j}^{2}) \]
\[\begin{align*} \int \alpha_{j}N(\beta_{j};\mu_{j}, s_{j}^2) \log \alpha_{j}N(\beta_{j};\mu_{j}, s_{j}^2)d\beta_{j}d\mathbf{\gamma} &= \sum_{j=1}^{p}\alpha_{j} \int N(\beta_{j}; \mu_{j}, s_{j}^2) \log \alpha_{j}d\beta_{j} + \sum_{j=1}^{p} \alpha_{j}\int N(\beta_{j};\mu_{j}, s_{j}^2) \log N(\beta_{j};\mu_{j}, s_{j}^2)d\beta_{j} \\ &= \sum_{j=1}^{p} \alpha_{j}\log \alpha_{j} + \alpha_{j}\left[-\frac{1}{2} \left(1 + \log 2\pi s_{j}^2\right)\right] \end{align*}\] \[\begin{align*} \int \alpha_{j}N(\beta_{j};\mu_{j}, s_{j}^2) \log \frac{1}{p} N(\beta_{j}; 0, \sigma^2\sigma_{\beta}^2) d\beta_{j}d\mathbf{\gamma} &= \sum_{j=1}^{p} \alpha_{j} \int N(\beta_{j};\mu_{j}, s_{j}^2) \log \frac{1}{p}d\beta_{j} + \sum_{j=1}^{p}\alpha_{j} \int N(\beta_{j};\mu_{j}, s_{j}^2) \log N(\beta_{j};0, \sigma^2\sigma_{\beta}^2)d\beta_{j} \\ &= \sum_{j=1}^{p} \alpha_{j}\log\frac{1}{p} + \alpha_{j}\left[-\frac{1}{2}\log 2\pi\sigma^2\sigma_{\beta}^2 - \frac{s_{j}^2 + \mu_{j}^2}{2\sigma^2\sigma_{\beta}^2} \right] \end{align*}\] \[\begin{align*} D_{KL}(q||p) &= \sum_{j=1}^{p}\alpha_{j}\log \frac{\alpha_{j}}{1/p} - \alpha_{j} \left[ \frac{1}{2} \left(1 + \log 2\pi s_{j}^2\right)-\frac{1}{2}\log 2\pi\sigma^2\sigma_{\beta}^2 - \frac{s_{j}^2 + \mu_{j}^2}{2\sigma^2\sigma_{\beta}^2} \right] \\ &= \sum_{j=1}^{p} \alpha_{j}\log \frac{\alpha_{j}}{1/p} - \frac{\alpha_{j}}{2} \left[ 1+\log\frac{s_{j}^2}{\sigma^2\sigma_{\beta}^2} - \frac{s_{j}^2 + \mu_{j}^2}{\sigma^2\sigma_{\beta}^2} \right] \end{align*}\]\[ \color{red}{ \mathcal{L}(q) = -\frac{n}{2}\log 2\pi\sigma^2 - \frac{1}{2\sigma^2}\left[\|\mathbf{y} - X\mathbf{r}\|^{2} - \|X\mathbf{r}\|^{2} + \sum_{l}\sum_{j=1}^{p}(X^{T}X)_{jj}\mathbb{E}(\beta_{lj}^2)\right] - \sum_{j=1}^{p}\alpha_{j}\log \frac{\alpha_{j}}{1/p} + \sum_{j=1}^{p} \frac{\alpha_{j}}{2} \left[ 1+\log\frac{s_{j}^2}{\sigma^2\sigma_{\beta}^2} - \frac{s_{j}^2 + \mu_{j}^2}{\sigma^2\sigma_{\beta}^2} \right]} \]
\[ p(B, \Gamma) = \prod_{l=1}^{L} p(\mathbf{\beta}_{l}, \mathbf{\gamma}_{l}) = \prod_{l=1}^{L} p(\mathbf{\beta}_{l}|\mathbf{\gamma}_{l})p(\mathbf{\gamma}_{l}) = \prod_{l=1}^{L} \frac{1}{p}\prod_{j=1}^{p}N(\beta_{lj};0,\sigma^2\sigma_{\beta}^2) \]
\[ q(B, \Gamma) = \prod_{l=1}^{L} q(\mathbf{\beta}_{l}, \mathbf{\gamma}_{l}) = \prod_{l=1}^{L} q(\mathbf{\beta}_{l}|\mathbf{\gamma}_{l})p(\mathbf{\gamma}_{l}) = \prod_{l=1}^{L} \alpha_{lj}\prod_{j=1}^{p}N(\beta_{lj};\mu_{lj},s_{lj}^2) \]
\[ p(\mathbf{y},B, \Gamma|X, \sigma^2, \sigma_{\beta}^2, \pi) = \left( \frac{1}{\sqrt{2\pi \sigma^2}} \right)^{n} \exp \left( -\frac{\|\mathbf{y} - X\mathbf{\beta}_{1}-\cdots-X\mathbf{\beta}_{L}\|^2}{2\sigma^2} \right) \prod_{l=1}^{L}N(\beta_{lj}; 0, \sigma^2\sigma_{\beta}^2)\frac{1}{p} \]
\(\mathbf{\beta} = rowSum(B)\)
\[\begin{align*} \mathcal{L}(q) &= \mathbb{E}_{q}\left[ \log N_{n}(\mathbf{y} - X\mathbf{\beta}; 0, \sigma^2I) \right] + \sum_{l=1}^{L}\mathbb{E}_{q}\left[ \log \frac{p_{l}(\mathbf{\beta}_{l}, \mathbf{\gamma}_{l})}{q_{l}(\mathbf{\beta}_{l}, \mathbf{\gamma}_{l})} \right] \quad p_{l}(\mathbf{\beta}_{l}, \mathbf{\gamma}_{l}) = \text{prior of }\mathbf{\beta}_{l} \ \mathbf{\gamma}_{l} \\ &= -\frac{n}{2}\log 2\pi\sigma^2 - \frac{1}{2\sigma^2}\left[\|\mathbf{y}^{T} - X\mathbf{r}\|^2 - \sum_{l=1}^{L}\|X\mathbf{r}_{l}\|^{2} + \sum_{l}\sum_{j=1}^{p}(X^{T}X)_{jj}\mathbb{E}(\beta_{lj}^2)\right] - \sum_{l=1}^{L} D_{KL}(q_{l}||p_{l}) \end{align*}\]\(\mathbb{E}_{q}(\beta_{lj}) = \alpha_{lj} * \mu_{lj}\) \(\quad\) \(\mathbb{E}_{q}(\beta_{j}) =\mathbb{E}_{q}(\sum_{l=1}^{L} \beta_{lj}) = \sum_{l=1}^{L}\alpha_{lj} * \mu_{lj}\)
\(\mathbb{V}ar_{q}(\beta_{lj}) = \alpha_{lj}(s_{lj}^2 + \mu_{lj}^2) - (\alpha_{lj}\mu_{lj})^2\)
\[ \color{red}{ \mathcal{L}(q) = -\frac{n}{2}\log 2\pi\sigma^2 - \frac{1}{2\sigma^2}\left[\|\mathbf{y} - X\mathbf{r}\|^2 - \sum_{l=1}^{L}\|X\mathbf{r}_{l}\|^{2} + \sum_{l}\sum_{j=1}^{p}(X^{T}X)_{jj}\mathbb{E}(\beta_{lj}^2)\right] - \sum_{l=1}^{L}\sum_{j=1}^{p} \alpha_{lj}\log \frac{\alpha_{lj}}{1/p} + \sum_{l=1}^{L}\sum_{j=1}^{p} \frac{\alpha_{lj}}{2} \left[ 1+\log\frac{s_{lj}^2}{\sigma^2\sigma_{\beta}^2} - \frac{s_{lj}^2 + \mu_{lj}^2}{\sigma^2\sigma_{\beta}^2} \right]} \]
In susieR package, the KL Divergence in the lower bound is computed in another way, but it is equivalent.
The expectation at the right hand side is the log-likelihood for the single effect regression at the optimum.
\[ \color{red}{ \mathcal{L}(q) = -\frac{n}{2}\log 2\pi\sigma^2 - \frac{1}{2\sigma^2}\left[\|\mathbf{y} - X\mathbf{r}\|^2 - \sum_{l=1}^{L}\|X\mathbf{r}_{l}\|^{2} + \sum_{l}\sum_{j=1}^{p}(X^{T}X)_{jj}\mathbb{E}(\beta_{lj}^2)\right] + \sum_{l=1}^{L}\left[ \log p_{l}(\mathbf{y}_{resid,-l}|\sigma^2, \sigma_{\beta}^2) + \frac{n}{2}\log 2\pi\sigma^2 + \frac{1}{2\sigma^2} \mathbb{E}_{q_{l}}\|\mathbf{y}_{resid,-l} - X\mathbf{\beta}_{l}\|^2\right] } \]
sessionInfo()R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.6
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
loaded via a namespace (and not attached):
[1] workflowr_1.1.1 Rcpp_0.12.18 digest_0.6.15
[4] rprojroot_1.3-2 R.methodsS3_1.7.1 backports_1.1.2
[7] git2r_0.23.0 magrittr_1.5 evaluate_0.11
[10] stringi_1.2.4 whisker_0.3-2 R.oo_1.22.0
[13] R.utils_2.6.0 rmarkdown_1.10 tools_3.5.1
[16] stringr_1.3.1 yaml_2.2.0 compiler_3.5.1
[19] htmltools_0.3.6 knitr_1.20 This reproducible R Markdown analysis was created with workflowr 1.1.1