Last updated: 2017-01-15
Code version: 92038b69c4d891ab8162d54630d341956dc6e573
Note that rendering this RMarkdown document may take a few minutes because it involves loading and processing tables with millions of rows.
First, we load the necessary libraries.
library(dscr)
library(ashr)
library(reshape2)
library(ggplot2)
library(magrittr)
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, unionlibrary(xtable)Compile the results on the simulated data sets for the summaries below.
load("../output/dsc-shrink-files/res.RData")
coverthresh = 0.05 # threshold at which we look at coverage
findthresh = 0.95 # threshold at we define a discovery significant
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 > findthresh)
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 > findthresh)
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")Overall proportion of negative findings is 0.0776.
Overall proportion of positive findings is 0.0656.
Table of lower tail for all observations:
print(xtabs(lt ~ method + scenario,
reslong %>% group_by(scenario,method) %>%
summarize(lt = mean(value < coverthresh))) %>% round(2)) scenario
method big-normal-nn bimodal bimodal-nn flat-top flat-top-nn
ash.hu 0.06 0.04 0.05 0.06 0.06
ash.hu.nocxx 0.06 0.04 0.05 0.06 0.06
ash.hu.nw2 0.05 0.04 0.05 0.05 0.05
ash.hu.s 0.05 0.07 0.05 0.08 0.06
ash.n 0.05 0.04 0.05 0.05 0.06
ash.n.nw2 0.05 0.04 0.05 0.05 0.05
ash.n.s 0.05 0.04 0.05 0.05 0.05
ash.u 0.06 0.04 0.05 0.06 0.06
ash.u.nw2 0.05 0.04 0.05 0.05 0.05
ash.u.s 0.05 0.04 0.05 0.05 0.05
bayes 0.05 0.04 0.05 0.05 0.05
mixfdr.enull 0.00 0.05 0.05 0.06 0.05
mixfdr.tnull 0.00 0.12 0.20 0.16 0.24
scenario
method near-normal near-normal-nn skew skew-nn spiky spiky-nn
ash.hu 0.07 0.07 0.06 0.06 0.12 0.15
ash.hu.nocxx 0.07 0.07 0.06 0.06 0.12 0.15
ash.hu.nw2 0.05 0.06 0.05 0.06 0.08 0.09
ash.hu.s 0.08 0.06 0.08 0.06 0.12 0.08
ash.n 0.06 0.06 0.06 0.07 0.10 0.11
ash.n.nw2 0.05 0.05 0.06 0.07 0.08 0.08
ash.n.s 0.05 0.05 0.05 0.07 0.05 0.06
ash.u 0.07 0.07 0.07 0.08 0.13 0.15
ash.u.nw2 0.06 0.06 0.06 0.07 0.09 0.10
ash.u.s 0.05 0.06 0.06 0.07 0.06 0.07
bayes 0.04 0.05 0.04 0.05 0.05 0.05
mixfdr.enull 0.06 0.05 0.06 0.06 0.05 0.05
mixfdr.tnull 0.17 0.28 0.14 0.21 0.22 0.40Table of lower tail of positive findings. Because of the unimodal assumption and the favoritism toward the null, this should assess problems with “over shrinkage” toward 0.
print(xtabs(lt ~ method + scenario,
reslong.pos %>% group_by(scenario,method) %>%
summarize(lt = mean(value < coverthresh))) %>% round(2)) scenario
method big-normal-nn bimodal bimodal-nn flat-top flat-top-nn
ash.hu 0.05 0.05 0.04 0.07 0.06
ash.hu.nocxx 0.05 0.05 0.04 0.07 0.06
ash.hu.nw2 0.05 0.05 0.04 0.06 0.05
ash.hu.s 0.05 0.18 0.04 0.48 0.05
ash.n 0.05 0.04 0.03 0.06 0.05
ash.n.nw2 0.05 0.05 0.03 0.05 0.04
ash.n.s 0.05 0.06 0.03 0.08 0.04
ash.u 0.05 0.05 0.04 0.07 0.06
ash.u.nw2 0.05 0.05 0.04 0.06 0.04
ash.u.s 0.05 0.05 0.04 0.08 0.05
bayes 0.05 0.05 0.05 0.05 0.05
scenario
method near-normal near-normal-nn skew skew-nn spiky spiky-nn
ash.hu 0.08 0.07 0.08 0.08 0.08 0.08
ash.hu.nocxx 0.08 0.07 0.08 0.08 0.08 0.08
ash.hu.nw2 0.07 0.06 0.07 0.06 0.07 0.07
ash.hu.s 0.40 0.06 0.46 0.06 0.66 0.13
ash.n 0.06 0.06 0.14 0.12 0.06 0.07
ash.n.nw2 0.05 0.05 0.12 0.10 0.05 0.06
ash.n.s 0.06 0.05 0.12 0.08 0.06 0.05
ash.u 0.07 0.07 0.16 0.14 0.07 0.08
ash.u.nw2 0.06 0.06 0.13 0.11 0.06 0.07
ash.u.s 0.07 0.05 0.12 0.09 0.07 0.06
bayes 0.05 0.05 0.06 0.05 0.04 0.05Table of lower tail of negative findings. This should indicate problems with tail behaviour of \(g\). The uniform methods tend to over-shrink.
print(xtabs(lt ~ method + scenario,reslong.neg %>%
group_by(scenario,method) %>%
summarize(lt = mean(value < coverthresh))) %>% round(2)) scenario
method big-normal-nn bimodal bimodal-nn flat-top flat-top-nn
ash.hu 0.06 0.06 0.04 0.08 0.08
ash.hu.nocxx 0.06 0.06 0.04 0.08 0.08
ash.hu.nw2 0.06 0.06 0.04 0.10 0.09
ash.hu.s 0.06 0.06 0.05 0.09 0.08
ash.n 0.05 0.02 0.02 0.00 0.00
ash.n.nw2 0.05 0.03 0.02 0.01 0.00
ash.n.s 0.05 0.03 0.02 0.02 0.01
ash.u 0.05 0.06 0.04 0.07 0.08
ash.u.nw2 0.06 0.06 0.04 0.12 0.10
ash.u.s 0.06 0.06 0.04 0.10 0.08
bayes 0.05 0.05 0.05 0.05 0.05
scenario
method near-normal near-normal-nn skew skew-nn spiky spiky-nn
ash.hu 0.13 0.12 0.07 0.07 0.13 0.12
ash.hu.nocxx 0.13 0.12 0.07 0.07 0.13 0.12
ash.hu.nw2 0.11 0.10 0.07 0.07 0.13 0.12
ash.hu.s 0.08 0.09 0.06 0.06 0.11 0.08
ash.n 0.06 0.06 0.06 0.07 0.07 0.07
ash.n.nw2 0.06 0.06 0.07 0.07 0.05 0.07
ash.n.s 0.05 0.06 0.07 0.07 0.05 0.07
ash.u 0.10 0.11 0.09 0.09 0.12 0.11
ash.u.nw2 0.09 0.09 0.08 0.08 0.13 0.10
ash.u.s 0.08 0.07 0.08 0.08 0.11 0.09
bayes 0.05 0.05 0.05 0.05 0.05 0.05Compile some more summary tables in Latex format.
save_latex_coverage_table =
function(df, methodnames, filename,
SCENARIONAMES = c("spiky","near-normal","flat-top","skew",
"big-normal","bimodal"),
switch = FALSE) {
df$method <- factor(df$method,levels = methodnames)
df$scenario <- factor(df$scenario,levels = SCENARIONAMES)
mat <- as.matrix(xtabs(lt~method+scenario,df))
if (switch)
mat <- 1 - mat
mat <- xtable(mat,digits = rep(2,ncol(mat) + 1))
write(print(mat,sanitize.text.function = function (x) x,
floating = FALSE,hline.after = NULL,
add.to.row = list(pos = list(-1,0,nrow(mat)),
command = c('\\toprule ',
'\\midrule ',
'\\bottomrule '))),file = filename)
}
save_latex_coverage_table(reslong.neg %>% group_by(scenario,method) %>%
summarize(lt = mean(value < coverthresh)),
c("ash.n","ash.u","ash.hu"),
"table/coverage_neg.tex")% latex table generated in R 3.3.2 by xtable 1.8-2 package
% Sun Jan 15 17:52:28 2017
\begin{tabular}{rrrrrrr}
\toprule & spiky & near-normal & flat-top & skew & big-normal & bimodal \\
\midrule ash.n & 0.07 & 0.06 & 0.00 & 0.06 & 0.00 & 0.02 \\
ash.u & 0.12 & 0.10 & 0.07 & 0.09 & 0.00 & 0.06 \\
ash.hu & 0.13 & 0.13 & 0.08 & 0.07 & 0.00 & 0.06 \\
\bottomrule \end{tabular}save_latex_coverage_table(reslong.pos %>% group_by(scenario,method) %>%
summarize(lt = mean(value < coverthresh)),
c("ash.n","ash.u","ash.hu"),
"table/coverage_pos.tex")% latex table generated in R 3.3.2 by xtable 1.8-2 package
% Sun Jan 15 17:52:29 2017
\begin{tabular}{rrrrrrr}
\toprule & spiky & near-normal & flat-top & skew & big-normal & bimodal \\
\midrule ash.n & 0.06 & 0.06 & 0.06 & 0.14 & 0.00 & 0.04 \\
ash.u & 0.07 & 0.07 & 0.07 & 0.16 & 0.00 & 0.05 \\
ash.hu & 0.08 & 0.08 & 0.07 & 0.08 & 0.00 & 0.05 \\
\bottomrule \end{tabular}save_latex_coverage_table(reslong %>% group_by(scenario,method) %>%
summarize(lt = mean(value < coverthresh)),
c("ash.n","ash.u","ash.hu"),
"table/coverage_all.tex")% latex table generated in R 3.3.2 by xtable 1.8-2 package
% Sun Jan 15 17:52:31 2017
\begin{tabular}{rrrrrrr}
\toprule & spiky & near-normal & flat-top & skew & big-normal & bimodal \\
\midrule ash.n & 0.10 & 0.06 & 0.05 & 0.06 & 0.00 & 0.04 \\
ash.u & 0.13 & 0.07 & 0.06 & 0.07 & 0.00 & 0.04 \\
ash.hu & 0.12 & 0.07 & 0.06 & 0.06 & 0.00 & 0.04 \\
\bottomrule \end{tabular}save_latex_coverage_table(reslong.neg %>% group_by(scenario,method) %>%
summarize(lt = mean(value < coverthresh)),
c("ash.n.s","ash.u.s","ash.hu.s"),
"table/coverage_neg_nopen.tex")% latex table generated in R 3.3.2 by xtable 1.8-2 package
% Sun Jan 15 17:52:31 2017
\begin{tabular}{rrrrrrr}
\toprule & spiky & near-normal & flat-top & skew & big-normal & bimodal \\
\midrule ash.n.s & 0.05 & 0.05 & 0.02 & 0.07 & 0.00 & 0.03 \\
ash.u.s & 0.11 & 0.08 & 0.10 & 0.08 & 0.00 & 0.06 \\
ash.hu.s & 0.11 & 0.08 & 0.09 & 0.06 & 0.00 & 0.06 \\
\bottomrule \end{tabular}save_latex_coverage_table(reslong.pos %>% group_by(scenario,method) %>%
summarize(lt = mean(value < coverthresh)),
c("ash.n.s","ash.u.s","ash.hu.s"),
"table/coverage_pos_nopen.tex")% latex table generated in R 3.3.2 by xtable 1.8-2 package
% Sun Jan 15 17:52:31 2017
\begin{tabular}{rrrrrrr}
\toprule & spiky & near-normal & flat-top & skew & big-normal & bimodal \\
\midrule ash.n.s & 0.06 & 0.06 & 0.08 & 0.12 & 0.00 & 0.06 \\
ash.u.s & 0.07 & 0.07 & 0.08 & 0.12 & 0.00 & 0.05 \\
ash.hu.s & 0.66 & 0.40 & 0.48 & 0.46 & 0.00 & 0.18 \\
\bottomrule \end{tabular}save_latex_coverage_table(reslong %>% group_by(scenario,method) %>%
summarize(lt = mean(value < coverthresh)),
c("ash.n.s","ash.u.s","ash.hu.s"),
"table/coverage_all_nopen.tex")% latex table generated in R 3.3.2 by xtable 1.8-2 package
% Sun Jan 15 17:52:34 2017
\begin{tabular}{rrrrrrr}
\toprule & spiky & near-normal & flat-top & skew & big-normal & bimodal \\
\midrule ash.n.s & 0.05 & 0.05 & 0.05 & 0.05 & 0.00 & 0.04 \\
ash.u.s & 0.06 & 0.05 & 0.05 & 0.06 & 0.00 & 0.04 \\
ash.hu.s & 0.12 & 0.08 & 0.08 & 0.08 & 0.00 & 0.07 \\
\bottomrule \end{tabular}These tables show right tail instead of the left tail.
save_latex_coverage_table(reslong.neg %>% group_by(scenario,method) %>%
summarize(lt = mean(value < coverthresh)),
c("ash.n","ash.u","ash.hu"),
"table/scoverage_neg.tex",switch = TRUE)% latex table generated in R 3.3.2 by xtable 1.8-2 package
% Sun Jan 15 17:52:34 2017
\begin{tabular}{rrrrrrr}
\toprule & spiky & near-normal & flat-top & skew & big-normal & bimodal \\
\midrule ash.n & 0.93 & 0.94 & 1.00 & 0.94 & 1.00 & 0.98 \\
ash.u & 0.88 & 0.90 & 0.93 & 0.91 & 1.00 & 0.94 \\
ash.hu & 0.87 & 0.87 & 0.92 & 0.93 & 1.00 & 0.94 \\
\bottomrule \end{tabular}save_latex_coverage_table(reslong.pos %>% group_by(scenario,method) %>%
summarize(lt = mean(value < coverthresh)),
c("ash.n","ash.u","ash.hu"),
"table/scoverage_pos.tex",switch = TRUE)% latex table generated in R 3.3.2 by xtable 1.8-2 package
% Sun Jan 15 17:52:34 2017
\begin{tabular}{rrrrrrr}
\toprule & spiky & near-normal & flat-top & skew & big-normal & bimodal \\
\midrule ash.n & 0.94 & 0.94 & 0.94 & 0.86 & 1.00 & 0.96 \\
ash.u & 0.93 & 0.93 & 0.93 & 0.84 & 1.00 & 0.95 \\
ash.hu & 0.92 & 0.92 & 0.93 & 0.92 & 1.00 & 0.95 \\
\bottomrule \end{tabular}save_latex_coverage_table(reslong %>% group_by(scenario,method) %>%
summarize(lt = mean(value < coverthresh)),
c("ash.n","ash.u","ash.hu"),
"table/scoverage_all.tex",switch=TRUE)% latex table generated in R 3.3.2 by xtable 1.8-2 package
% Sun Jan 15 17:52:37 2017
\begin{tabular}{rrrrrrr}
\toprule & spiky & near-normal & flat-top & skew & big-normal & bimodal \\
\midrule ash.n & 0.90 & 0.94 & 0.95 & 0.94 & 1.00 & 0.96 \\
ash.u & 0.87 & 0.93 & 0.94 & 0.93 & 1.00 & 0.96 \\
ash.hu & 0.88 & 0.93 & 0.94 & 0.94 & 1.00 & 0.96 \\
\bottomrule \end{tabular}save_latex_coverage_table(reslong.neg %>% group_by(scenario,method) %>%
summarize(lt = mean(value < coverthresh)),
c("ash.n.s","ash.u.s","ash.hu.s"),
"table/scoverage_neg_nopen.tex",switch = TRUE)% latex table generated in R 3.3.2 by xtable 1.8-2 package
% Sun Jan 15 17:52:37 2017
\begin{tabular}{rrrrrrr}
\toprule & spiky & near-normal & flat-top & skew & big-normal & bimodal \\
\midrule ash.n.s & 0.95 & 0.95 & 0.98 & 0.93 & 1.00 & 0.97 \\
ash.u.s & 0.89 & 0.92 & 0.90 & 0.92 & 1.00 & 0.94 \\
ash.hu.s & 0.89 & 0.92 & 0.91 & 0.94 & 1.00 & 0.94 \\
\bottomrule \end{tabular}save_latex_coverage_table(reslong.pos %>% group_by(scenario,method) %>%
summarize(lt = mean(value < coverthresh)),
c("ash.n.s","ash.u.s","ash.hu.s"),
"table/scoverage_pos_nopen.tex",switch = TRUE)% latex table generated in R 3.3.2 by xtable 1.8-2 package
% Sun Jan 15 17:52:37 2017
\begin{tabular}{rrrrrrr}
\toprule & spiky & near-normal & flat-top & skew & big-normal & bimodal \\
\midrule ash.n.s & 0.94 & 0.94 & 0.92 & 0.88 & 1.00 & 0.94 \\
ash.u.s & 0.93 & 0.93 & 0.92 & 0.88 & 1.00 & 0.95 \\
ash.hu.s & 0.34 & 0.60 & 0.52 & 0.54 & 1.00 & 0.82 \\
\bottomrule \end{tabular}save_latex_coverage_table(reslong %>% group_by(scenario,method) %>%
summarize(lt = mean(value < coverthresh)),
c("ash.n.s","ash.u.s","ash.hu.s"),
"table/scoverage_all_nopen.tex",switch = TRUE)% latex table generated in R 3.3.2 by xtable 1.8-2 package
% Sun Jan 15 17:52:40 2017
\begin{tabular}{rrrrrrr}
\toprule & spiky & near-normal & flat-top & skew & big-normal & bimodal \\
\midrule ash.n.s & 0.95 & 0.95 & 0.95 & 0.95 & 1.00 & 0.96 \\
ash.u.s & 0.94 & 0.95 & 0.95 & 0.94 & 1.00 & 0.96 \\
ash.hu.s & 0.88 & 0.92 & 0.92 & 0.92 & 1.00 & 0.93 \\
\bottomrule \end{tabular}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 xtable_1.8-2 dplyr_0.5.0 magrittr_1.5
[5] ggplot2_2.2.0 reshape2_1.4.2 ashr_2.0.5 dscr_0.1.1
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 DBI_0.5-1 shiny_0.14.2
[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 backports_1.0.4 scales_0.4.1
[22] codetools_0.2-15 htmltools_0.3.5 MASS_7.3-45
[25] assertthat_0.1 mime_0.5 colorspace_1.2-6
[28] httpuv_1.3.3 stringi_1.1.2 lazyeval_0.2.0
[31] munsell_0.4.3 doParallel_1.0.10 pscl_1.4.9
[34] truncnorm_1.0-7 SQUAREM_2016.8-2