Simulation Design

In the simulation, \[c_{j2...R} = c_{j1} + \delta_{j}\] \[\hat{c}_{j}|c_{j} \sim N_{R}(c_{j}, \frac{1}{2}I)\] Let L be the contrast matrix. Therefore, \[\hat{\delta}_{j}|\delta_{j} = L\hat{c}_{j}|c_{j} \sim N_{R-1}(\delta_{j}, \frac{1}{2}LL')\]

We first generate the data:

set.seed(2018)
data = sim.contrast2(nsamp = 21000, ncond = 8)

This simulation routine creates a dataset with 8 conditions, and four different types of deviations \(\delta_{j}\): null, independent among conditions, condition-specific in condition 2, and shared (equal effects in all conditions). The data contains 10% non-null effects \[\delta_{j} \sim \frac{9}{10}N(0,0)+\frac{1}{30}N(0,I)+\frac{1}{30}N(0,11^{T})+\frac{1}{30}N(0,e_{1}e_{1}^{T})\]

Set up the contrast matrix and the mash contrast data object

R = 8
L = diag(R-1)
L = cbind(rep(-1, R-1), L)
row.names(L) = c('2-1','3-1','4-1','5-1','6-1','7-1','8-1')

mash_data = mash_set_data(Bhat=data$Chat, Shat=data$Shat)
mash_data_L = mash_set_data_contrast(mash_data, L)

Set up the covariance matrices:

# canonical
U.c = cov_canonical(mash_data_L)
# data driven
m.1by1 = mash_1by1(mash_data_L, alpha=1)
strong = get_significant_results(m.1by1,0.05)
# only 6 strong samples

# center Z
# mash_data_L.center = mash_data_L
# mash_data_L.center$Bhat = mash_data_L$Bhat/mash_data_L$BShat # obtain z
# mash_data_L.center$BShat = matrix(1, nrow(mash_data_L$Bhat),ncol(mash_data_L$Bhat))
# mash_data_L.center$Bhat = apply(mash_data_L.center$Bhat, 2, function(x) x - mean(x))
# U.pca = cov_pca(mash_data_L.center,5,strong)
# U.ed = cov_ed(mash_data_L.center, U.pca, strong)

Fit mashcontrast model

mashcontrast.model = mash(mash_data_L, U.c, algorithm.version = 'R')

 - Computing 21000 x 181 likelihood matrix.
 - Likelihood calculations took 1.66 seconds.
 - Fitting model with 181 mixture components.
 - Model fitting took 4.76 seconds.
 - Computing posterior matrices.
 - Computation allocated took 0.18 seconds.

print(get_loglik(mashcontrast.model),digits = 10)

[1] -184139.2001

Use get_significant_results to find the indices of effects that are “significant”:

length(get_significant_results(mashcontrast.model))

[1] 111

The number of false positive is 0.

barplot(get_estimated_pi(mashcontrast.model),las = 2,cex.names = 0.7)

Comparing discovery result with \(\hat{\delta}_{j} \sim N_{7}(\delta_{j},I)\):

Indep.data = mash_set_data(Bhat = data$Chat%*%t(L))
Indep.model = mash(Indep.data, U.c, algorithm.version = 'R')

 - Computing 21000 x 181 likelihood matrix.
 - Likelihood calculations took 1.53 seconds.
 - Fitting model with 181 mixture components.
 - Model fitting took 5.09 seconds.
 - Computing posterior matrices.
 - Computation allocated took 0.15 seconds.

print(get_loglik(Indep.model),digits = 10)

[1] -194668.353

length(get_significant_results(Indep.model))

[1] 6349

There are more discoveries here, but lots of them are false discoveries. The number of false positive is 5549.

barplot(get_estimated_pi(Indep.model),las = 2, cex.names = 0.7)

The RRMSE plot:

We check the False Positive Rate and True Positive Rate. \[FPR = \frac{|N\cap S|}{|N|} \quad TPR = \frac{|CS\cap S|}{|T|} \]

Session information

sessionInfo()

R version 3.4.3 (2017-11-30)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.3

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/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] mashr_0.2-4 ashr_2.2-4 

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.15      knitr_1.17        magrittr_1.5     
 [4] REBayes_1.2       MASS_7.3-47       doParallel_1.0.11
 [7] pscl_1.5.2        SQUAREM_2017.10-1 lattice_0.20-35  
[10] foreach_1.4.4     plyr_1.8.4        stringr_1.2.0    
[13] tools_3.4.3       parallel_3.4.3    grid_3.4.3       
[16] rmeta_2.16        htmltools_0.3.6   iterators_1.0.9  
[19] assertthat_0.2.0  yaml_2.1.16       rprojroot_1.2    
[22] digest_0.6.13     Matrix_1.2-12     codetools_0.2-15 
[25] evaluate_0.10.1   rmarkdown_1.8     stringi_1.1.6    
[28] compiler_3.4.3    Rmosek_8.0.69     backports_1.1.2  
[31] mvtnorm_1.0-7     truncnorm_1.0-7