Last updated: 2018-03-01

Code version: d62d09a

R=5

library(mashr)
Loading required package: ashr
set.seed(2018)
data = simple_sims(500, err_sd = 0.1)

data.miss = data
missing.row = sample(1:nrow(data$Bhat), nrow(data$Bhat)/4)
missing.ind = matrix(0, nrow = nrow(data$B), ncol=ncol(data$B))
missing.ind[missing.row,] = 1

missing = which(missing.ind == 1)
data.miss$Bhat[missing] = NA
data.miss$Shat[missing] = NA
missing1 = is.na(data.miss$Bhat[,1])
# Set the missing value with large standard deviation
data.large.dev = data.miss
data.large.dev$Bhat[is.na(data.miss$Bhat)] = 0
data.large.dev$Shat[is.na(data.miss$Shat)] = 1000

EE model

Warning in REBayes::KWDual(A, rep(1, k), normalize(w), control = control): estimated mixing distribution has some negative values:
               consider reducing rtol
Warning in mixIP(matrix_lik = structure(c(1, 0.204132576568413,
0.736530341806702, : Optimization step yields mixture weights that are
either too small, or negative; weights have been corrected and renormalized
after the optimization.

The log likelihood for model with missing data is -2.111711210^{4}. The weights learned from the data is correct. There are 985 significant samples.

The log likelihood for model without missing vlaues is -1550.3763871. The estimated weights are same. There are 985 significant samples. The missing values with large deviation contain little information.

The estimates for other observed values are similar, as we expected.

sqrt(sum((mash.model.missing$result$PosteriorMean[-missing] - mash.model.missing.na$result$PosteriorMean)^2))
[1] 3.422263e-07

EZ model

FIXME: 'compute_posterior_matrices' in Rcpp does not transfer EZ to EE
Warning in REBayes::KWDual(A, rep(1, k), normalize(w), control = control): estimated mixing distribution has some negative values:
               consider reducing rtol
Warning in mixIP(matrix_lik = structure(c(1, 0.204132576568412,
0.736530341806701, : Optimization step yields mixture weights that are
either too small, or negative; weights have been corrected and renormalized
after the optimization.
FIXME: 'compute_posterior_matrices' in Rcpp does not transfer EZ to EE
par(mfrow=c(1,2))
barplot(get_estimated_pi(mash.model.missing.1), las=2, cex.names = 0.7, main='with missing values')
barplot(get_estimated_pi(mash.model.missing.na.1), las=2, cex.names = 0.7, main='without missing values')

par(mfrow=c(1,1))

The log likelihood for the model with missing values is -2.159352210^{4}. There are 974 significant samples.

The log likelihood for the model without missing values is -1550.3763871. The estimated weights are different. There are 985 significant samples. There are 2 cases for z close to 0. One is the observed effect is close to zero; the other is the value is missing, it has large standard deviation. The EZ model cannot distinguish this two cases. They contain same amount of information.

R=60

EE model

With missing values

The log likelihood is -2.124305510^{5}. The weights learned from the data is correct. There are 1020 significant samples.

Without missing values

Warning in REBayes::KWDual(A, rep(1, k), normalize(w), control = control): estimated mixing distribution has some negative values:
               consider reducing rtol
Warning in mixIP(matrix_lik = structure(c(0.240920767436896,
0.337563330121546, : Optimization step yields mixture weights that are
either too small, or negative; weights have been corrected and renormalized
after the optimization.

The log likelihood is 2.187068210^{4}. The estimated weights are very weird and incorrect. It ignores the null matrix. There are 1014 significant samples.

The estimates for other observed values are different, as we expected.

sqrt(sum((mash.model.missing$result$PosteriorMean[-missing] - mash.model.missing.na$result$PosteriorMean)^2))
[1] 0.7759291

EZ model

With missing values

Warning in REBayes::KWDual(A, rep(1, k), normalize(w), control = control): estimated mixing distribution has some negative values:
               consider reducing rtol
Warning in mixIP(matrix_lik = structure(c(0.240920767436896,
0.33756333012155, : Optimization step yields mixture weights that are
either too small, or negative; weights have been corrected and renormalized
after the optimization.
FIXME: 'compute_posterior_matrices' in Rcpp does not transfer EZ to EE

The log likelihood is -2.13689310^{5}. There are 1010 significant samples.

Without missing values

mash.data.missing.na.1 = mash_set_data(Bhat=data.miss.na$Bhat, Shat= data.miss.na$Shat, alpha = 1)
U.c = cov_canonical(mash.data.missing.na.1)

mash.model.missing.na.1 = mash(mash.data.missing.na.1, U.c, verbose=FALSE)
Warning in REBayes::KWDual(A, rep(1, k), normalize(w), control = control): estimated mixing distribution has some negative values:
               consider reducing rtol
Warning in mixIP(matrix_lik = structure(c(0.240920767436896,
0.33756333012155, : Optimization step yields mixture weights that are
either too small, or negative; weights have been corrected and renormalized
after the optimization.
FIXME: 'compute_posterior_matrices' in Rcpp does not transfer EZ to EE
barplot(get_estimated_pi(mash.model.missing.na.1), las=2, cex.names = 0.7)

The log likelihood is 2.187068210^{4}. The estimated weights are different. There are 1014 significant samples.

Covariance

The weired covariance structure is caused by the high dimension of R.

If we separate conditions into several groups, the covariance structure will be correct.

We show this for the EZ model R=60.

FIXME: 'compute_posterior_matrices' in Rcpp does not transfer EZ to EE
FIXME: 'compute_posterior_matrices' in Rcpp does not transfer EZ to EE
FIXME: 'compute_posterior_matrices' in Rcpp does not transfer EZ to EE

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        git2r_0.20.0      htmltools_0.3.6  
[19] iterators_1.0.9   assertthat_0.2.0  yaml_2.1.16      
[22] rprojroot_1.2     digest_0.6.13     Matrix_1.2-12    
[25] codetools_0.2-15  evaluate_0.10.1   rmarkdown_1.8    
[28] stringi_1.1.6     compiler_3.4.3    Rmosek_8.0.69    
[31] backports_1.1.2   mvtnorm_1.0-7     truncnorm_1.0-7

This R Markdown site was created with workflowr