Last updated: 2017-01-15

Code version: b05decc05714ddcb20aa04e56b557d5123d178fb

First, we load the necessary libraries.

library(ashr)
library(reshape2)
library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(dscr)

Compile the tables which we will use to examine the results of the simulation experiments.

load("../output/dsc-shrink-files/res.RData")

neglong = 
  res$negprob %>% 
    select(-user.self,-sys.self,-elapsed,-user.child,-sys.child) %>%
    melt(id.vars = c("method","scenario","seed",".id"),value.name = "negprob") %>%
    filter(negprob > 0.95)
 
poslong = 
  res$posprob %>% 
    select(-user.self,-sys.self,-elapsed,-user.child,-sys.child) %>%
    melt(id.vars = c("method","scenario","seed",".id"),value.name = "posprob") %>%
    filter(posprob > 0.95)

reslong = 
  res$cdf_score %>% 
    select(-user.self,-sys.self,-elapsed,-user.child,-sys.child) %>%    
    melt(id.vars = c("method","scenario","seed",".id"))

reslong.pos = inner_join(reslong,poslong)
Joining, by = c("method", "scenario", "seed", ".id", "variable")
reslong.neg = inner_join(reslong,neglong)
Joining, by = c("method", "scenario", "seed", ".id", "variable")
ash.hu.ft = reslong.pos %>% filter(method == "ash.hu.s" & 
                                   scenario == "flat-top")

Extract an example to illustrate coverage of ash.hu.

eg = load_example(dsc_shrink,42,"flat-top",methodname = "ash.hu.s",
                  homedir = "../code/dsc-shrink")
out        <- eg$output$fitted_g
class(out) <- "list"
out        <- as.data.frame(out)
print(subset(out,pi > 0.01))
        pi      a       b
10 0.02742 -1.795 0.00000
15 0.42292  0.000 0.07932
16 0.50655  0.000 0.11218
24 0.04312  0.000 1.79483

Notice how almost all the inferred weight is on a small positive component. As a result, false sign rate will be small, and there will be a strong tendency to overestimate zero effects. This leads to coverage problems observed.

Now let’s look at an example with u and spiky which seems to be somewhat badly calibrated for the negative discoveries.

ash.u.s.spiky = reslong.neg %>% filter(method == "ash.u.s" & 
                                       scenario == "spiky")
ash.n.s.spiky = reslong.neg %>% filter(method == "ash.n.s" & 
                                         scenario == "spiky")
hist(ash.u.s.spiky$value,nclass = 100,
     main = paste("histogram of quantile where observation falls in its CI"))

hist(ash.n.s.spiky$value,nclass = 100,
     main = paste("histogram of quantile where observation falls in its CI"))

So what seems to be happening here is that the uniform tail is too short; when observation falls outside of this tail it gets a zero quantile of posterior interval.

Can we find an example illustrating this trend? (Note: the next chunk below did not work for me, so I se* eval=FALSE, include=FALSE. Perhaps this can be fixed at some point. -Peter)

For comparison, here are the positive discoveries; here they are not too bad.

ash.u.s.spiky = reslong.pos %>% filter(method == "ash.u.s" & scenario == "spiky")
ash.n.s.spiky = reslong.pos %>% filter(method == "ash.n.s" & scenario == "spiky")
hist(ash.n.s.spiky$value)

hist(ash.u.s.spiky$value)

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] knitr_1.15.1   dscr_0.1.1     dplyr_0.5.0    ggplot2_2.2.0 
[5] reshape2_1.4.2 ashr_2.0.5    

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.8       plyr_1.8.4        iterators_1.0.8  
 [4] tools_3.3.2       digest_0.6.10     evaluate_0.10    
 [7] tibble_1.2        gtable_0.2.0      lattice_0.20-34  
[10] foreach_1.4.3     shiny_0.14.2      DBI_0.5-1        
[13] yaml_2.1.14       parallel_3.3.2    stringr_1.1.0    
[16] rprojroot_1.1     grid_3.3.2        R6_2.2.0         
[19] rmarkdown_1.3     magrittr_1.5      backports_1.0.4  
[22] scales_0.4.1      codetools_0.2-15  htmltools_0.3.5  
[25] MASS_7.3-45       assertthat_0.1    mime_0.5         
[28] colorspace_1.2-6  xtable_1.8-2      httpuv_1.3.3     
[31] stringi_1.1.2     lazyeval_0.2.0    pscl_1.4.9       
[34] doParallel_1.0.10 munsell_0.4.3     truncnorm_1.0-7  
[37] SQUAREM_2016.8-2