Last updated: 2021-03-08
Checks: 7 0
Knit directory: mmbr-rss-dsc/
This reproducible R Markdown analysis was created with workflowr (version 1.6.2). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
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.
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.
The command set.seed(20200227) 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.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version b7ae692. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use 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: .DS_Store
Ignored: .Rhistory
Ignored: .Rproj.user/
Ignored: data/.DS_Store
Ignored: output/.DS_Store
Untracked files:
Untracked: data/ENSG00000140265.12.Multi_Tissues.rds
Untracked: data/FastQTLSumStats.mash.FL_PC3.rds
Untracked: data/mmbr_ukb_simulation/
Untracked: data/mmbr_ukb_simulation_problem/
Untracked: data/ukb_rss_naive_lfsr_problem.rds
Untracked: data/ukb_rss_naive_pip_problem.rds
Untracked: output/GTExprofile_res.rds
Untracked: output/GTExprofile_resL1.rds
Untracked: output/GTExprofile_resL1_elbo.rds
Untracked: output/GTExprofile_resL3.rds
Untracked: output/GTExprofile_resL3_elbo.rds
Untracked: output/GTExprofile_res_elbo.rds
Untracked: output/GTExprofile_resapprox.rds
Untracked: output/GTExprofile_resapproxL1.rds
Untracked: output/GTExprofile_resapproxL1_elbo.rds
Untracked: output/GTExprofile_resapproxL3.rds
Untracked: output/GTExprofile_resapproxL3_elbo.rds
Untracked: output/GTExprofile_resapprox_elbo.rds
Untracked: output/GTExprofile_resapproxdiag.rds
Untracked: output/GTExprofile_resapproxdiagL1.rds
Untracked: output/GTExprofile_resapproxdiagL1_elbo.rds
Untracked: output/GTExprofile_resapproxdiagL3.rds
Untracked: output/GTExprofile_resapproxdiagL3_elbo.rds
Untracked: output/GTExprofile_resapproxdiag_elbo.rds
Untracked: output/GTExprofile_resdiag.rds
Untracked: output/mmbr_gtex_res.Rprof
Untracked: output/mmbr_gtex_res_approx.Rprof
Untracked: output/mmbr_gtex_res_approx_diag.Rprof
Untracked: output/mmbr_gtex_res_diag.Rprof
Untracked: output/mnm_missing_output.20200527.rds
Untracked: output/test
Untracked: output/tiny_data_211_cond2L2.gif
Untracked: output/tiny_data_211_cond2L2.pdf
Untracked: output/tiny_data_211_cond2L3.gif
Untracked: output/tiny_data_211_cond2L3.pdf
Untracked: output/tiny_data_211_cond2initL3.gif
Untracked: output/tiny_data_211_cond2initL3.pdf
Untracked: output/ukb_rss_20200107/
Unstaged changes:
Modified: analysis/mmbr_missing_rss_problem1.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.
These are the previous versions of the repository in which changes were made to the R Markdown (analysis/mmbr_rss_ukb_simulation_20210107_naiveprior_lfsr_problem.Rmd) and HTML (docs/mmbr_rss_ukb_simulation_20210107_naiveprior_lfsr_problem.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| Rmd | b7ae692 | zouyuxin | 2021-03-08 | wflow_publish(“analysis/mmbr_rss_ukb_simulation_20210107_naiveprior_lfsr_problem.Rmd”) |
| html | 0262806 | zouyuxin | 2021-03-08 | Build site. |
| Rmd | 73c71e4 | zouyuxin | 2021-03-08 | wflow_publish(“analysis/mmbr_rss_ukb_simulation_20210107_naiveprior_lfsr_problem.Rmd”) |
We found some inflation in condition specific pip using canonical prior.
library(mmbr)Loading required package: mashrLoading required package: ashrLoading required package: susieRdat = readRDS('data/ukb_rss_naive_lfsr_problem.rds')
idx = which(rowSums(dat$true_coef != 0)>0)
priorU = dat$priors$naiveThere is one causal SNP with PVE 0.0005. The causal SNP has effect in condition 10, 11. This pattern is not included in canonical priors.
round(dat$true_coef[idx,], 4) [1] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0224
[11] 0.0224 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000All models captures the true signal. The global PIP makes sense. The problem is the lfsr in null conditions at the signal.
With L = 1,
m_init = mmbr::create_mash_prior(mixture_prior = list(matrices=priorU$xUlist, weights=priorU$pi),
null_weight=priorU$null_weight, max_mixture_len=-1)
resultsuff1 = mmbr::msusie_suff_stat(XtX = dat$XtX, XtY = dat$XtY, YtY = dat$YtY, N = dat$N, L=1,
prior_variance=m_init, residual_variance=dat$resid,
compute_objective=TRUE, estimate_residual_variance=F,
estimate_prior_variance=T, estimate_prior_method='EM',
precompute_covariances=T, n_thread=1, max_iter=100,
track_fit = T, verbosity = 0) ## prior unchange, estimate priorThe lfsr at the causal SNP is
round(resultsuff1$lfsr[1507,],3) [1] 0.343 0.071 0.437 0.226 0.301 0.403 0.240 0.425 0.044 0.000 0.000 0.230
[13] 0.113 0.388 0.337 0.434 0.038 0.120 0.304 0.462The lfsr at the causal SNP in condition 10, 11 are very small. In other conditions, the lfsr are around 0.3, which turns to condition specific pip ~0.6-0.7.
At the causal SNP, the posterior mixture weights are mainly in identity matrix and shared matrix (off diagonal 0.25).
names(priorU$xUlist)[which(resultsuff1$mixture_weights[1,1507,]>0.2)-1][1] "shared_1" "shared_2"The estimated prior scalar is
resultsuff1$V[1] 5.294432e-05With L = 2,
resultsuff2 = mmbr::msusie_suff_stat(XtX = dat$XtX, XtY = dat$XtY,
YtY = dat$YtY, N = dat$N, L=2,
prior_variance=m_init, residual_variance=dat$resid,
compute_objective=TRUE, estimate_residual_variance=F,
estimate_prior_variance=T, estimate_prior_method='EM',
precompute_covariances=T, n_thread=1, max_iter=10,
track_fit = T, verbosity = 0) ## prior unchange, estimate prior --> workThe lfsr at the causal SNP is
round(resultsuff2$lfsr[1507,],3) [1] 1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1The condition specific pip for those null conditions are 0.
The first single effect is the causal SNP with posterior mixture weight 1 at singleton 11. The second single effect is the causal SNP with posterior mixture weight 1 at singleton 10.
The estimated prior is
resultsuff2$V[1] 0.0004569579 0.0004530954If we scale the prior variances by 1/sample size (1/248980),
Usmall = lapply(priorU$xUlist, function(x) x/dat$N)
m_initsmall = mmbr::create_mash_prior(mixture_prior = list(matrices=Usmall, weights=priorU$pi),
null_weight=priorU$null_weight, max_mixture_len=-1)
resultsuff2small = mmbr::msusie_suff_stat(XtX = dat$XtX, XtY = dat$XtY,
YtY = dat$YtY, N = dat$N, L=2,
prior_variance=m_initsmall, residual_variance=dat$resid,
compute_objective=TRUE, estimate_residual_variance=F,
estimate_prior_variance=T, estimate_prior_method='EM',
precompute_covariances=T, n_thread=1, max_iter=100,
track_fit = T, verbosity = 0) ## prior/N, estimate priorThe lfsr becomes
round(resultsuff2small$lfsr[1507,],3) [1] 0.343 0.071 0.438 0.226 0.301 0.403 0.240 0.425 0.044 0.000 0.000 0.230
[13] 0.113 0.388 0.337 0.435 0.038 0.120 0.304 0.462The posterior mixture weights are mainly in identity matrix and shared matrix (off diagonal 0.25).
Some conclusion: if the prior doesn’t include the signal pattern and the prior scale is smaller than the signal, the posteior weights will focus on shared patterns. This cause the inflation in condition specific pip.
We try cannonical + a small diagonal (how to choose the small diagonal scalar?). I use 0.01 here.
# U+sI
Us = lapply(priorU$xUlist, function(x) x + 0.01*diag(nrow(x)))
m_inits = mmbr::create_mash_prior(mixture_prior = list(matrices=Us, weights=priorU$pi),
null_weight=priorU$null_weight, max_mixture_len=-1)
resultsuff2s = mmbr::msusie_suff_stat(XtX = dat$XtX, XtY = dat$XtY,
YtY = dat$YtY, N = dat$N, L=2,
prior_variance=m_inits, residual_variance=dat$resid,
compute_objective=TRUE, estimate_residual_variance=F,
estimate_prior_variance=T, estimate_prior_method='EM',
precompute_covariances=T, n_thread=1, max_iter=20,
track_fit = T, verbosity = 0)The lfsr is
round(resultsuff2s$lfsr[1507,],3) [1] 0.331 0.071 0.458 0.221 0.324 0.424 0.262 0.445 0.045 0.000 0.000 0.224
[13] 0.112 0.374 0.360 0.455 0.048 0.138 0.326 0.445The posteior weights for the first single effect is a mixture of singleton 10 and 11.
We fit RSS model with L = 2,
ldeigen = eigen(cov2cor(dat$XtX), symmetric = T)
m_init = mmbr::create_mash_prior(mixture_prior = list(matrices=priorU$xUlist, weights=priorU$pi),
null_weight=priorU$null_weight, max_mixture_len=-1)
resultrss2 = mmbr::msusie_rss(Z = dat$Z, eigenR = ldeigen, L=2,
prior_variance=m_init, residual_variance=dat$resid,
compute_objective=TRUE, estimate_residual_variance=F,
estimate_prior_variance=T, estimate_prior_method='EM',
precompute_covariances=T, n_thread=1, max_iter=100,
track_fit = T, verbosity = 0) ## estimate priorThe lfsr is
round(resultrss2$lfsr[1507,],3) [1] 0.343 0.071 0.438 0.226 0.301 0.403 0.240 0.425 0.044 0.000 0.000 0.230
[13] 0.113 0.388 0.337 0.435 0.038 0.120 0.304 0.462The mixture posterior weights for the causal SNP is on identity matrix and shared matrix (off diagonal 0.25).
If we add the true pattern in the prior,
Ut = priorU$xUlist
Ut$true = matrix(0,20,20)
Ut$true[10:11,10:11] = 1
m_initt = mmbr::create_mash_prior(mixture_prior = list(matrices=Ut, weights=priorU$pi),
null_weight=priorU$null_weight, max_mixture_len=-1)
resultrss2t = mmbr::msusie_rss(Z = dat$Z, eigenR = ldeigen, L=2,
prior_variance=m_initt, residual_variance=dat$resid,
compute_objective=TRUE, estimate_residual_variance=F,
estimate_prior_variance=T, estimate_prior_method='EM',
precompute_covariances=T, n_thread=1, max_iter=100,
track_fit = T, verbosity = 0) The lfsr is
round(resultrss2t$lfsr[1507,],3) [1] 1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1The mixture posterior weights for the causal SNP is all on the true pattern.
We scale the prior using magnitud of Z scores (max(abs(z))^2).
# scale z
sz = max(abs(dat$Z))^2
Uz = lapply(priorU$xUlist, function(x) x*sz)
m_initz = mmbr::create_mash_prior(mixture_prior = list(matrices=Uz, weights=priorU$pi),
null_weight=priorU$null_weight, max_mixture_len=-1)
resultrss2z = mmbr::msusie_rss(Z = dat$Z, eigenR = ldeigen, L=2,
prior_variance=m_initz, residual_variance=dat$resid,
compute_objective=TRUE, estimate_residual_variance=F,
estimate_prior_variance=T, estimate_prior_method='EM',
precompute_covariances=T, n_thread=1, max_iter=100,
track_fit = T, verbosity = 0)The lfsr is
round(resultrss2z$lfsr[1507,],3) [1] 0.343 0.071 0.438 0.226 0.301 0.403 0.240 0.425 0.044 0.000 0.000 0.230
[13] 0.113 0.388 0.337 0.435 0.038 0.120 0.304 0.462We scale the prior using 5000.
# scale 5000
U5000 = lapply(priorU$xUlist, function(x) x*5000)
m_initz5000 = mmbr::create_mash_prior(mixture_prior = list(matrices=U5000, weights=priorU$pi),
null_weight=priorU$null_weight, max_mixture_len=-1)
resultrss25000 = mmbr::msusie_rss(Z = dat$Z, eigenR = ldeigen, L=2,
prior_variance=m_initz5000, residual_variance=dat$resid,
compute_objective=TRUE, estimate_residual_variance=F,
estimate_prior_variance=T, estimate_prior_method='EM',
precompute_covariances=T, n_thread=1, max_iter=100,
track_fit = T, verbosity = 0)The lfsr is
round(resultrss25000$lfsr[1507,],3) [1] 1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1The lfsr looks right!
sessionInfo()R version 4.0.3 (2020-10-10)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur 10.16
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/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
other attached packages:
[1] mmbr_0.0.2.0429 susieR_0.10.0 mashr_0.2.41 ashr_2.2-51
[5] workflowr_1.6.2
loaded via a namespace (and not attached):
[1] progress_1.2.2 softImpute_1.4 tidyselect_1.1.0 xfun_0.19
[5] purrr_0.3.4 reshape2_1.4.4 lattice_0.20-41 colorspace_2.0-0
[9] vctrs_0.3.6 generics_0.1.0 htmltools_0.5.0 yaml_2.2.1
[13] rlang_0.4.10 mixsqp_0.3-46 later_1.1.0.1 pillar_1.4.7
[17] glue_1.4.2 matrixStats_0.58.0 lifecycle_1.0.0 plyr_1.8.6
[21] stringr_1.4.0 munsell_0.5.0 gtable_0.3.0 mvtnorm_1.1-1
[25] evaluate_0.14 knitr_1.30 httpuv_1.5.4 invgamma_1.1
[29] irlba_2.3.3 Rcpp_1.0.6 promises_1.1.1 scales_1.1.1
[33] rmeta_3.0 truncnorm_1.0-8 abind_1.4-5 fs_1.5.0
[37] hms_1.0.0 flashr_0.6-7 ggplot2_3.3.3 digest_0.6.27
[41] stringi_1.5.3 dplyr_1.0.2 grid_4.0.3 rprojroot_2.0.2
[45] tools_4.0.3 magrittr_2.0.1 tibble_3.0.6 crayon_1.4.1
[49] whisker_0.4 pkgconfig_2.0.3 ellipsis_0.3.1 Matrix_1.2-18
[53] prettyunits_1.1.1 SQUAREM_2021.1 reshape_0.8.8 assertthat_0.2.1
[57] rmarkdown_2.5 R6_2.5.0 git2r_0.27.1 compiler_4.0.3