To try out some simulations that don’t match the canonical covariance matrices and illustrate how the data driven matrices help.
Here the function
simple_sims_2 simulates data in five
conditions with just two types of effect:
shared effects only in the first two conditions; and
shared effects only in the last three conditions.
library(ashr) library(mashr) set.seed(1) = simple_sims2(1000,1) simdata = cbind(c(1,1,0,0,0),c(1,1,0,0,0),rep(0,5),rep(0,5),rep(0,5)) true.U1 = cbind(rep(0,5),rep(0,5),c(0,0,1,1,1),c(0,0,1,1,1),c(0,0,1,1,1)) true.U2 = list(true.U1 = true.U1, true.U2 = true.U2)U.true
Run 1-by-1 to add the strong signals and ED covariances.
= mash_set_data(simdata$Bhat, simdata$Shat) data .1by1 = mash_1by1(data) m= get_significant_results(m.1by1) strong = cov_canonical(data) U.c = cov_pca(data,5,strong) U.pca = cov_ed(data,U.pca,strong) U.ed # Computes covariance matrices based on extreme deconvolution, # initialized from PCA. = mash(data, U.c) m.c = mash(data, U.ed) m.ed = mash(data, c(U.c,U.ed)) m.c.ed = mash(data, U.true) m.true print(get_loglik(m.c),digits = 10) print(get_loglik(m.ed),digits = 10) print(get_loglik(m.c.ed),digits = 10) print(get_loglik(m.true),digits = 10)
The log-likelihood is much better from data-driven than canonical covariances. This is good! Indeed, here the data-driven fit is very slightly better fit than the true matrices, but only very slightly.