Last updated: 2017-01-15

Code version: 94a3347615361216b39b8a9333bfa8354794e0ff

Here we check that mixfdr actually does OK at fitting the bimodal distribution when noise is small compared with bimodality.

Simulate bimodal data with well-separated peaks. The histogram shows the multi-modal nature of the underlying g is clear.

source("../code/dsc-shrink/datamakers/datamaker.R")
sim.bimodal =
  rnormmix_datamaker(args = list(g = normalmix(c(0.5,0.5),c(-5,5),c(1,1)),
                                 min_pi0   = 0.4,
                                 max_pi0   = 0.4,
                                 nsamp     = 1000,
                                 betahatsd = 1))
hist(sim.bimodal$input$betahat,nclass = 20,xlim = c(-10,10),prob = TRUE,
     main="histogram of simulated betahat values")

Now run mixfdr on the simulated data. Note that the fit captures the multi-modal nature of g, although mixfdr continues to overestimate pi0.

source("../code/dsc-shrink/methods/mixfdr.wrapper.R")
sim.bimodal.mixfdroutput =
  mixfdr.wrapper(sim.bimodal$input,args = list(theonull = TRUE))
Warning in (function (x, J = 3, P = NA, noiseSD = NA, theonull = FALSE, :
Assuming known noise noiseSD = 1 . If needed rerun with noiseSD = NA to fit
noiseSD.
Fitting preliminary model 
Fitting final model 
Warning in (function (x, J = 3, P = NA, noiseSD = NA, theonull = FALSE, :
Null proportion pi0 is small. Consider increasing penalization and/or using
an empirical null.
Warning in (function (x, J = 3, P = NA, noiseSD = NA, theonull = FALSE, :
Using an empirical null with a fitted noiseSD gives a substantially
different model. Consider rerunning with theonull = FALSE and noiseSD = NA.

Fitted Model: J = 3 groups
----------------------------
null?        TRUE   FALSE   FALSE   

pi =        0.5090  0.2369  0.2541  

mu =         0.000   5.020  -5.148  

sigma =     1.000   1.330   1.401   

noiseSD =   1   
g = mixfdr2fitted.g(sim.bimodal.mixfdroutput)$fitted.g
x = seq(-10,10,length = 100)
par(cex.main = 1,font.main = 1)
plot(ecdf(sim.bimodal$meta$beta),cex = 0.6,
     main = "CDF of true beta: empirical (black) and fitted (red)")
lines(x,mixcdf(g,x),lwd = 2,lty = "dotted",col = "orangered")

Session information

sessionInfo()
R version 3.3.2 (2016-10-31)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: macOS Sierra 10.12.2

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] mixfdr_1.0   ashr_2.0.5   knitr_1.15.1

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.8       lattice_0.20-34   codetools_0.2-15 
 [4] digest_0.6.10     foreach_1.4.3     rprojroot_1.1    
 [7] truncnorm_1.0-7   MASS_7.3-45       grid_3.3.2       
[10] backports_1.0.4   magrittr_1.5      evaluate_0.10    
[13] stringi_1.1.2     pscl_1.4.9        doParallel_1.0.10
[16] rmarkdown_1.3     iterators_1.0.8   tools_3.3.2      
[19] stringr_1.1.0     parallel_3.3.2    yaml_2.1.14      
[22] SQUAREM_2016.8-2  htmltools_0.3.5