Last updated: 2018-10-24
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.
✔ Repository version: 9f8adcd
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
Ignored: docs/figure/Test.Rmd/
Untracked files:
Untracked: analysis/mashEMlk.Rmd
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
File | Version | Author | Date | Message |
---|---|---|---|---|
Rmd | 9f8adcd | zouyuxin | 2018-10-24 | wflow_publish(“analysis/mashrEM.Rmd”) |
html | d02d30d | zouyuxin | 2018-10-23 | Build site. |
Rmd | 3f88eab | zouyuxin | 2018-10-23 | wflow_publish(“analysis/mashrEM.Rmd”) |
html | e5aab50 | zouyuxin | 2018-10-23 | Build site. |
Rmd | 9875360 | zouyuxin | 2018-10-23 | wflow_publish(“analysis/mashrEM.Rmd”) |
The mashr model is \[ \hat{\mathbf{b}}_{j}|\mathbf{b}_{j} \sim N_{R}(\hat{\mathbf{b}}_{j}; \mathbf{b}_{j}, \mathbf{S}_{j}\mathbf{V}\mathbf{S}_{j}) \\ \mathbf{b}_{j} \sim \sum_{l,k} \pi_{l,k} N_{R}(\mathbf{b}_{j}; \mathbf{0}, \omega_{l}\mathbf{U}_{k}) = \sum_{p=1}^{P} \pi_{p} N_{R}(\mathbf{b}_{j}; \mathbf{0}, \Sigma_{p}) \]
\[ p(\hat{\mathbf{B}}, \mathbf{B}, \mathbf{z}) h(\boldsymbol{\pi}) = \prod_{j=1}^{J} \prod_{p = 1}^{P}\left[\pi_{p} N_{R}(\hat{\mathbf{b}}_{j}; \mathbf{b}_{j}, \mathbf{S}_{j}\mathbf{V}\mathbf{S}_{j})N_{R}(\mathbf{b}_{j}; \mathbf{0}, \Sigma_{p})\right]^{\mathbb{I}(z_{j} = p)} \prod_{p=1}^{P} \pi_{p}^{\lambda_{p}-1} \]
\[ \begin{align*} P(z_{j}=p, \mathbf{b}_{j}|\hat{\mathbf{b}}_{j}) &= \frac{P(z_{j}=p, \mathbf{b}_{j},\hat{\mathbf{b}}_{j})}{P(\hat{\mathbf{b}}_{j})} = \frac{P(\hat{\mathbf{b}}_{j}|\mathbf{b}_{j})P(\mathbf{b}_{j}|z_{j}=p) P(z_{j}=p)}{P(\hat{\mathbf{b}}_{j})} \\ &= \frac{\pi_{p} N_{R}(\hat{\mathbf{b}}_{j}; \mathbf{b}_{j}, \mathbf{S}_{j}\mathbf{V}\mathbf{S}_{j})N_{R}(\mathbf{b}_{j}; \mathbf{0}, \Sigma_{p})}{\sum_{p'}\pi_{p'} N_{R}(\hat{\mathbf{b}}_{j}; \mathbf{0}, \mathbf{S}_{j}\mathbf{V}\mathbf{S}_{j} + \Sigma_{p'})} \\ &= \frac{\pi_{p} N_{R}(\hat{\mathbf{b}}_{j}; \mathbf{0}, \mathbf{S}_{j}\mathbf{V}\mathbf{S}_{j} + \Sigma_{p})}{\sum_{p'}\pi_{p'} N_{R}(\hat{\mathbf{b}}_{j}; \mathbf{0}, \mathbf{S}_{j}\mathbf{V}\mathbf{S}_{j} + \Sigma_{p'})} \frac{N_{R}(\hat{\mathbf{b}}_{j}; \mathbf{b}_{j}, \mathbf{S}_{j}\mathbf{V}\mathbf{S}_{j})N_{R}(\mathbf{b}_{j}; \mathbf{0}, \Sigma_{p})}{N_{R}(\hat{\mathbf{b}}_{j}; \mathbf{0}, \mathbf{S}_{j}\mathbf{V}\mathbf{S}_{j} + \Sigma_{p})} \\ &= \gamma_{jp} P(\mathbf{b}_{j}|z_{j}=p, \hat{\mathbf{b}}_{j}) \\ &= P(z_{j}=p|\hat{\mathbf{b}}_{j}) P(\mathbf{b}_{j}|z_{j}=p, \hat{\mathbf{b}}_{j}) \end{align*} \]
E step: \[ \begin{align*} \mathbb{E}_{\mathbf{z}, \mathbf{B}|\hat{\mathbf{B}}}\log p(\hat{\mathbf{B}}, \mathbf{B}, \mathbf{z}) h(\boldsymbol{\pi}) &= \mathbb{E}_{\mathbf{z}, \mathbf{B}|\hat{\mathbf{B}}} \{ \sum_{j=1}^{J}\sum_{p = 1}^{P} \mathbb{I}(z_{j} = p)\left[\log \pi_{p} + \log N_{R}(\hat{\mathbf{b}}_{j}; \mathbf{b}_{j}, \mathbf{S}_{j}\mathbf{V}\mathbf{S}_{j}) + \log N_{R}(\mathbf{b}_{j}; \mathbf{0}, \Sigma_{p})\right] + \sum_{p=1}^{P} (\lambda_{p}-1) \log \pi_{p} \} \\ &= \sum_{j=1}^{J} \sum_{p=1}^{P} \gamma_{jp} \left[\log \pi_{p} - \frac{1}{2}\log |\mathbf{V}| - \log |\mathbf{S}_{j}| - \frac{1}{2}\mathbb{E}_{\mathbf{b}_{j}|\hat{\mathbf{b}}_{j}, z_{j}=p}\left((\hat{\mathbf{b}}_{j}-\mathbf{b}_{j})^{T}\mathbf{S}_{j}^{-1}\mathbf{V}^{-1}\mathbf{S}_{j}^{-1}(\hat{\mathbf{b}}_{j}-\mathbf{b}_{j})\right) - \frac{1}{2}\log |\Sigma_{p}| - \frac{1}{2}\mathbb{E}_{\mathbf{b}_{j}|\hat{\mathbf{b}}_{j}, z_{j}=p}\left(\mathbf{b}_{j}^{T}\Sigma_{p}^{-1}\mathbf{b}_{j} \right) \right] + \sum_{p=1}^{P} (\lambda_{p}-1)\log \pi_{p} \end{align*} \] M step:
Update \(\mathbf{\pi}\): \[ \hat{\pi}_{p} = \frac{(\sum_{j=1}^{J} \gamma_{jp}) + \lambda_{p} - 1}{J - P + \sum_{p=1}^{P}\lambda_{p} } \] Update \(\Sigma_{p}\): \[ \begin{align*} f(\Sigma_{p})'&= \sum_{j=1}^{J}\gamma_{jp} \left[-\frac{1}{2} \Sigma_{p}^{-1} + \frac{1}{2} \Sigma_{p}^{-1}\mathbb{E}_{\mathbf{b}_{j}|\hat{\mathbf{b}}_{j}, z_{j}=p}\left[\mathbf{b}_{j}\mathbf{b}_{j}^{T} \right]\Sigma_{p}^{-1} \right] = 0 \\ \sum_{j=1}^{J}\gamma_{jp}\Sigma_{p} &= \sum_{j=1}^{J} \gamma_{jp} \mathbb{E}_{\mathbf{b}_{j}|\hat{\mathbf{b}}_{j}, z_{j}=p}\left[\mathbf{b}_{j}\mathbf{b}_{j}^{T} \right] \\ \hat{\Sigma}_{p} &= \frac{\sum_{j=1}^{J} \gamma_{jp} \mathbb{E}\left[\mathbf{b}_{j}\mathbf{b}_{j}^{T} | \hat{\mathbf{b}}_{j}, z_{j}=p\right]}{\sum_{j=1}^{J}\gamma_{jp}} \end{align*} \]
Update \(\mathbf{V}\):
If we require the diagonal of \(\mathbf{V}\) be 1, we can decompose \(V = DCD\), C is the covariance matrix, D = \(diag(1/sqrt(C_{jj}))\). \[ \begin{align*} f(\mathbf{C}) &= \sum_{j=1}^{J} \sum_{p=1}^{P} \gamma_{jp} \left[- \frac{1}{2}\log |\mathbf{D}\mathbf{C}\mathbf{D}| - \frac{1}{2}\mathbb{E}_{\mathbf{b}_{j}|\hat{\mathbf{b}}_{j}, z_{j}=p}\left((\hat{\mathbf{b}}_{j}-\mathbf{b}_{j})^{T}\mathbf{S}_{j}^{-1}\mathbf{D}^{-1}\mathbf{C}^{-1}\mathbf{D}^{-1}\mathbf{S}_{j}^{-1}(\hat{\mathbf{b}}_{j}-\mathbf{b}_{j})\right) \right] \\ &= \sum_{j=1}^{J} \sum_{p=1}^{P} \gamma_{jp} \left[- \frac{1}{2}\log |\mathbf{C}| - \log |\mathbf{D}|- \frac{1}{2}\mathbb{E}_{\mathbf{b}_{j}|\hat{\mathbf{b}}_{j}, z_{j}=p}\left((\hat{\mathbf{b}}_{j}-\mathbf{b}_{j})^{T}\mathbf{S}_{j}^{-1}\mathbf{D}^{-1}\mathbf{C}^{-1}\mathbf{D}^{-1}\mathbf{S}_{j}^{-1}(\hat{\mathbf{b}}_{j}-\mathbf{b}_{j})\right) \right] \\ f(\mathbf{C})' &= \sum_{j=1}^{J} \sum_{p=1}^{P} \gamma_{jp}\left[ -\frac{1}{2} \mathbf{C}^{-1} + \frac{1}{2} \mathbf{C}^{-1} \mathbf{D}^{-1}\mathbf{S}_{j}^{-1}\mathbb{E}\left((\hat{\mathbf{b}}_{j}-\mathbf{b}_{j}) (\hat{\mathbf{b}}_{j}-\mathbf{b}_{j})^{T}|\hat{\mathbf{b}}_{j}, z_{j} = p \right)\mathbf{S}_{j}^{-1}\mathbf{D}^{-1} \mathbf{C}^{-1} \right] = 0 \\ \mathbf{C} &= \frac{1}{J} \sum_{j=1}^{J} \sum_{p=1}^{P} \gamma_{jp}\mathbf{D}^{-1}\mathbf{S}_{j}^{-1}\mathbb{E}\left((\hat{\mathbf{b}}_{j}-\mathbf{b}_{j}) (\hat{\mathbf{b}}_{j}-\mathbf{b}_{j})^{T}|\hat{\mathbf{b}}_{j}, z_{j} = p \right)\mathbf{S}_{j}^{-1}\mathbf{D}^{-1} \\ &= \frac{1}{J} \mathbf{D}^{-1}\sum_{j=1}^{J} \mathbf{S}_{j}^{-1}\mathbb{E}\left((\hat{\mathbf{b}}_{j}-\mathbf{b}_{j}) (\hat{\mathbf{b}}_{j}-\mathbf{b}_{j})^{T}|\hat{\mathbf{b}}_{j}\right)\mathbf{S}_{j}^{-1}\mathbf{D}^{-1} \end{align*} \]
We can update \(\mathbf{C}\) and \(\mathbf{V}\) as \[ \hat{\mathbf{C}}_{(t+1)} = \hat{\mathbf{D}}^{-1}_{(t)}\frac{1}{J} \left[\sum_{j=1}^{J} \mathbf{S}_{j}^{-1}\mathbb{E}\left[ (\hat{\mathbf{b}}_{j} - \mathbf{b}_{j})(\hat{\mathbf{b}}_{j} - \mathbf{b}_{j})^{T} | \hat{\mathbf{b}}_{j}\right]\mathbf{S}_{j}^{-1} \right] \hat{\mathbf{D}}^{-1}_{(t)} \\ \hat{\mathbf{D}}_{(t+1)} = diag(1/\sqrt{\hat{\mathbf{C}}_{(t+1)jj}}) \\ \hat{\mathbf{V}}_{(t+1)} = \hat{\mathbf{D}}_{(t+1)}\hat{\mathbf{C}}_{(t+1)}\hat{\mathbf{D}}_{(t+1)} \]
The resulting \(\hat{\mathbf{V}}_{(t+1)}\) is equivalent as \[ \hat{\mathbf{C}}_{(t+1)} =\frac{1}{J} \left[\sum_{j=1}^{J} \mathbf{S}_{j}^{-1}\mathbb{E}\left[ (\hat{\mathbf{b}}_{j} - \mathbf{b}_{j})(\hat{\mathbf{b}}_{j} - \mathbf{b}_{j})^{T} | \hat{\mathbf{b}}_{j}\right]\mathbf{S}_{j}^{-1} \right] \\ \hat{\mathbf{D}}_{(t+1)} = diag(1/\sqrt{\hat{\mathbf{C}}_{(t+1)jj}}) \\ \hat{\mathbf{V}}_{(t+1)} = \hat{\mathbf{D}}_{(t+1)}\hat{\mathbf{C}}_{(t+1)}\hat{\mathbf{D}}_{(t+1)} \]
This reproducible R Markdown analysis was created with workflowr 1.1.1