Simulation without signal

\[c_{j} = \mu_{j} 1\] \[\hat{c}_{j} \sim N_{R}(c_{j}, \frac{1}{2}I)\] Let L be the contrast matrix, which comparing each condition with the mean. \[L = \left(\begin{array}{c c c c} \frac{R-1}{R} & -\frac{1}{R} & \cdots & -\frac{1}{R} \\ -\frac{1}{R} & \frac{R-1}{R} & \cdots & -\frac{1}{R} \\ \vdots & \ddots & \ddots & \vdots \\ -\frac{1}{R} & \cdots & \frac{R-1}{R} & -\frac{1}{R} \end{array} \right)_{R-1 \times R}\] There are only R-1 rows, intead of R. Since \(c_{j,R}-\bar{c_{j}} = -\sum_{r=1}^{R-1} (c_{j,r}-\bar{c_{j}})\).

Therefore, \[\hat{\delta}_{j} = L\hat{c}_{j} \sim N_{R-1}(0, \frac{1}{2}LL')\]

We first generate the data:

set.seed(1)
data = sim_contrast1(nsamp = 10000, ncond = 8)

Mash contrast model

Set up the contrast matrix and the mash contrast data object

L = rbind(c(7/8,-1/8,-1/8,-1/8,-1/8,-1/8,-1/8,-1/8),
          c(-1/8,7/8,-1/8,-1/8,-1/8,-1/8,-1/8,-1/8),
          c(-1/8,-1/8,7/8,-1/8,-1/8,-1/8,-1/8,-1/8),
          c(-1/8,-1/8,-1/8,7/8,-1/8,-1/8,-1/8,-1/8),
          c(-1/8,-1/8,-1/8,-1/8,7/8,-1/8,-1/8,-1/8),
          c(-1/8,-1/8,-1/8,-1/8,-1/8,7/8,-1/8,-1/8),
          c(-1/8,-1/8,-1/8,-1/8,-1/8,-1/8,7/8,-1/8))
row.names(L) = seq(1,7)
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:

U.c = cov_canonical(mash_data_L)

Fit mashcontrast model

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

 - Computing 10000 x 169 likelihood matrix.
 - Likelihood calculations took 1.05 seconds.
 - Fitting model with 169 mixture components.
 - Model fitting took 1.94 seconds.
 - Computing posterior matrices.
 - Computation allocated took 0.05 seconds.

length(get_significant_results(mashcontrast.model))

[1] 0

There is no discovery, which is as we expected. The true deviations from mean are zero for all samples.

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

Subtract mean directly

If we subtract the mean from the data directly, ignoring the correlation structure. \[Var(\hat{c}_{j,r}-\bar{\hat{c}_{j}}) = \frac{1}{2} - \frac{1}{2R}\]

Indep.data = mash_set_data(Bhat = (data$Chat - apply(data$Chat,1, mean))[,1:7],
                           Shat = matrix(sqrt(0.5-1/(8*2)), nrow(data$Chat), ncol(data$Chat)-1))
U.c = cov_canonical(mash_data_L)
Indep.model = mash(Indep.data, U.c, algorithm.version = 'R')

 - Computing 10000 x 169 likelihood matrix.
 - Likelihood calculations took 0.86 seconds.
 - Fitting model with 169 mixture components.
 - Model fitting took 2.05 seconds.
 - Computing posterior matrices.
 - Computation allocated took 0.08 seconds.

length(get_significant_results(Indep.model))

[1] 0

There are no false positives.

barplot(get_estimated_pi(Indep.model),las = 2)

When there are no signal, subtracting mean directly from the data performs as good as the mashcommonbaseline model.

Session information

sessionInfo()

R version 3.4.4 (2018-03-15)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.4

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-8 ashr_2.2-7 

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.16      knitr_1.20        magrittr_1.5     
 [4] REBayes_1.3       MASS_7.3-50       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.3.0    
[13] tools_3.4.4       parallel_3.4.4    grid_3.4.4       
[16] rmeta_3.0         git2r_0.21.0      htmltools_0.3.6  
[19] iterators_1.0.9   assertthat_0.2.0  yaml_2.1.19      
[22] rprojroot_1.3-2   digest_0.6.15     Matrix_1.2-14    
[25] codetools_0.2-15  evaluate_0.10.1   rmarkdown_1.9    
[28] stringi_1.2.2     compiler_3.4.4    Rmosek_8.0.69    
[31] backports_1.1.2   mvtnorm_1.0-7     truncnorm_1.0-8