The inflated vs deflated noise figure in cashr

Introduction

Generate the inflated vs in-between vs deflated noise figure in the paper.

library(ggplot2)

FDP.list <- readRDS("../output/paper/simulation/FDP.list.rds")
TDP.list <- readRDS("../output/paper/simulation/TDP.list.rds")
z.list <- readRDS("../output/paper/simulation/z.list.rds")
pi0.list <- readRDS("../output/paper/simulation/pi0.list.rds")

pi0.vec <- c(0.5, 0.9, 0.99)
q.vec <- seq(0.001, 0.20, by = 0.001)
q <- 0.1
method.name.FDR <- c("cashr", "BH", "qvalue", "ashr", "locfdr")
method.col.FDR <- scales::hue_pal()(length(method.name.FDR))[c(5, 1, 2, 4, 3)]

#####################################

boxplot.quantile <- function(x) {
  r <- quantile(x, probs = c(0.10, 0.25, 0.5, 0.75, 0.90))
  names(r) <- c("ymin", "lower", "middle", "upper", "ymax")
  return(r)
}

boxplot.quantile.sq <- function (x) {
  r <- sqrt(quantile(x^2, probs = c(0.10, 0.25, 0.5, 0.75, 0.90)))
  names(r) <- c("ymin", "lower", "middle", "upper", "ymax")
  return(r)
}

mean.sq <- function (x) {
  r <- sqrt(mean(x^2))
  return(r)
}

mysqrt_trans <- function() {
  scales::trans_new("mysqrt", 
                    transform = base::sqrt,
                    inverse = function(x) ifelse(x<0, 0, x^2),
                    domain = c(0, Inf))
}

###############################

sd.z <- sapply(z.list, sd)
Noise <- cut(sd.z, breaks = c(0, quantile(sd.z, probs = 1 : 2 / 3), Inf), labels = c("Pseudo Deflation", "In-between", "Pseudo Inflation"))

###############################

pi0.0.9 <- which(pi0.list == 0.9)
sd.z.0.9 <- sd.z[pi0.0.9]
typical.noise <- pi0.0.9[order(sd.z.0.9)[floor(quantile(seq(sd.z.0.9), c(0.1, 0.5, 0.9)))]]

z.list.sel <- z.list[typical.noise]
names(z.list.sel) <- c("Pseudo Deflation", "In-between", "Pseudo Inflation")
z.sep.ggdata <- reshape2::melt(z.list.sel, value.name = "z")
z.sep.ggdata$L1 <- factor(z.sep.ggdata$L1, levels = c("Pseudo Deflation", "In-between", "Pseudo Inflation"))
z.sep.plot <- ggplot(data = z.sep.ggdata, aes(x = z)) +
  geom_histogram(aes(y = ..density..), binwidth = 0.2) +
  facet_wrap(~L1, nrow = 1) +
  stat_function(fun = dnorm, aes(color = "N(0,1)"), lwd = 1.5, show.legend = TRUE) +
#  scale_color_manual(values = "blue") +
  theme(# axis.title.x = element_text(size = 15),
        # axis.text.x = element_text(size = 10),
        axis.title.y = element_text(size = 15),
        axis.text.y = element_text(size = 10),
        strip.text = element_text(size = 15),
      # strip.text = element_blank(),
        legend.position = "left",
      # legend.title = element_blank(),
        legend.text = element_text(size = 12),
        legend.key = element_blank()
        )

#########################################

FDP.q <- FDP.list[[which(q.vec == q)]]
FDP.q.noise.mat <- rbind.data.frame(
  cbind.data.frame(Noise = rep("All", length(Noise)),
                   pi0 = factor(do.call(rbind, pi0.list)),
                   FDP.q),
  cbind.data.frame(Noise,
                   pi0 = factor(do.call(rbind, pi0.list)),
                   FDP.q)
)
FDP.q.ggdata <- reshape2::melt(FDP.q.noise.mat, id.vars = c("Noise", "pi0"), variable.name = "Method", value.name = "FDP")

FDP.q.ggdata$Method <- factor(FDP.q.ggdata$Method, levels = c("CASH", "BHq", "qvalue", "ASH", "locfdr"))
FDP.q.ggdata$Method <- plyr::mapvalues(FDP.q.ggdata$Method, from = c("CASH", "BHq", "qvalue", "ASH", "locfdr"), to = c("cashr", "BH", "qvalue", "ashr", "locfdr"))

FDP.q.ggdata.0.9 <- FDP.q.ggdata[FDP.q.ggdata$pi0 == 0.9, ]

################################################

TDP.q <- TDP.list[[which(q.vec == q)]]
TDP.q.noise.mat <- rbind.data.frame(
  cbind.data.frame(Noise = rep("All", length(Noise)),
                   pi0 = factor(do.call(rbind, pi0.list)),
                   TDP.q),
  cbind.data.frame(Noise,
                   pi0 = factor(do.call(rbind, pi0.list)),
                   TDP.q)
)
TDP.q.ggdata <- reshape2::melt(TDP.q.noise.mat, id.vars = c("Noise", "pi0"), variable.name = "Method", value.name = "TDP")

TDP.q.ggdata$Method <- factor(TDP.q.ggdata$Method, levels = c("CASH", "BHq", "qvalue", "ASH", "locfdr"))
TDP.q.ggdata$Method <- plyr::mapvalues(TDP.q.ggdata$Method, from = c("CASH", "BHq", "qvalue", "ASH", "locfdr"), to = c("cashr", "BH", "qvalue", "ashr", "locfdr"))

TDP.q.ggdata.0.9 <- TDP.q.ggdata[TDP.q.ggdata$pi0 == 0.9, ]

##############################################################

# FDP.q.all.sep.plot <- ggplot(data = FDP.q.ggdata, aes(x = pi0, y = FDP, fill = Method, color = Method)) +
#   stat_summary(fun.data = boxplot.quantile, geom = "boxplot", position = "dodge") +
#   stat_summary(fun.y = mean, geom = "point", position = position_dodge(width = 0.9), show.legend = TRUE) +
#   scale_color_manual(labels = method.name.FDR, values = method.col.FDR) +
#   scale_fill_manual(labels = method.name.FDR, values = alpha(method.col.FDR, 0.35)) +
#   facet_wrap(~Noise, nrow = 1) +
#   geom_hline(yintercept = q, col = "black", linetype = "dashed", size = 1) +
#   labs(x = expression(pi[0]), y = "FDP", title = bquote(paste("Nominal FDR = ", .(q)))) +
#   theme(plot.title = element_text(size = 12, hjust = 0),
#         axis.title.x = element_text(size = 15),
#         axis.text.x = element_text(size = 10),
#         axis.title.y = element_text(size = 15),
#         axis.text.y = element_text(size = 10),
#         strip.text = element_text(size = 15),
#         legend.position = "bottom",
#         legend.background = element_rect(color = "grey"),
#         legend.text = element_text(size = 12)
#         )

FDP.q.all.sep.plot <- ggplot(data = FDP.q.ggdata.0.9, aes(x = Method, y = FDP, fill = Method, color = Method)) +
  stat_summary(fun.data = boxplot.quantile, geom = "boxplot", aes(width = 0.5), position = position_dodge(), show.legend = FALSE) +
  stat_summary(fun.y = mean, geom = "point", position = position_dodge(width = 0.9), show.legend = FALSE) +
  scale_x_discrete(limits = rev(levels(FDP.q.ggdata.0.9$Method))) +
  coord_flip() +
  scale_color_manual(labels = method.name.FDR, values = method.col.FDR) +
  scale_fill_manual(labels = method.name.FDR, values = alpha(method.col.FDR, 0.35)) +
  facet_wrap(~Noise, nrow = 1) +
  geom_hline(yintercept = q, col = "black", linetype = "dashed", size = 1) +
  labs(y = "FDP", title = bquote(paste("Nominal FDR = ", .(q), " (", g[1], " is Gaussian; ", pi[0] == 0.9, ")"))) +
  expand_limits(y = 0) +
  theme(plot.title = element_text(size = 12, hjust = 0),
        axis.title.y = element_blank(),
        axis.text.y = element_text(size = 15),
        axis.title.x = element_text(size = 15),
        axis.text.x = element_text(size = 10),
        strip.text = element_text(size = 15),
        legend.position = "bottom",
        legend.background = element_rect(color = "grey"),
        legend.text = element_text(size = 12)
        )

Warning: Ignoring unknown aesthetics: width

FDP.sqrt.q.all.sep.plot <- ggplot(data = FDP.q.ggdata.0.9, aes(x = Method, y = FDP, fill = Method, color = Method)) +
  stat_summary(fun.data = boxplot.quantile.sq, geom = "boxplot", aes(width = 0.75), position = position_dodge(), show.legend = FALSE) +
  stat_summary(fun.y = mean.sq, geom = "point", position = position_dodge(width = 0.9), show.legend = FALSE, shape = 13, size = 3) +
  scale_x_discrete(limits = rev(levels(FDP.q.ggdata.0.9$Method))) +
  scale_y_continuous(trans = "mysqrt", breaks = c(0, 0.1, 0.2, 0.4, 0.6, 0.8)) +
  coord_flip() +
  scale_color_manual(labels = method.name.FDR, values = method.col.FDR) +
  scale_fill_manual(labels = method.name.FDR, values = alpha(method.col.FDR, 0.35)) +
  facet_wrap(~Noise, nrow = 1) +
  geom_hline(yintercept = q, col = "black", linetype = "dashed", size = 1) +
  labs(x = expression(pi[0]), y = "FDP", title = bquote(paste("Nominal FDR = ", .(q), " (", g[1], " is Gaussian; ", pi[0] == 0.9, ")"))) +
  expand_limits(y = 0) +
  theme(plot.title = element_text(size = 12, hjust = 0),
        axis.title.y = element_blank(),
        axis.text.y = element_text(size = 15),
        axis.title.x = element_text(size = 15),
        axis.text.x = element_text(size = 10),
        strip.text = element_text(size = 15),
        legend.position = "bottom",
        legend.background = element_rect(color = "grey"),
        legend.text = element_text(size = 12)
        )

Warning: Ignoring unknown aesthetics: width

###############################################

blank.ggdata <- data.frame()
blank.plot <- ggplot(data = blank.ggdata) +
  theme(panel.background = element_blank(),
        panel.border = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.background = element_blank()) +
  geom_blank()

z.sep.plot.save <- gridExtra::arrangeGrob(blank.plot, z.sep.plot, nrow = 1, widths = c(0.55, 3.5))

FDP.q.sep.plot.save <- gridExtra::arrangeGrob(z.sep.plot.save, FDP.q.all.sep.plot, heights = c(1, 1.2))

FDP.sqrt.q.sep.plot.save <- gridExtra::arrangeGrob(z.sep.plot.save, FDP.sqrt.q.all.sep.plot, heights = c(1, 1.2))

ggsave("../output/paper/FDP.q.sep.pdf", FDP.q.sep.plot.save, height = 4.5, width = 10)

ggsave("../output/paper/FDP.sqrt.q.sep.pdf", FDP.sqrt.q.sep.plot.save, height = 4.5, width = 10)

##############################################

TDP.q.all.sep.plot <- ggplot(data = TDP.q.ggdata.0.9, aes(x = Method, y = TDP, fill = Method, color = Method)) +
  stat_summary(fun.data = boxplot.quantile, geom = "boxplot", aes(width = 0.75), position = position_dodge(), show.legend = FALSE) +
  coord_flip() +
  stat_summary(fun.y = mean, geom = "point", position = position_dodge(width = 0.9), show.legend = FALSE, shape = 13, size = 3) +
  scale_x_discrete(limits = rev(levels(FDP.q.ggdata.0.9$Method))) +
  scale_color_manual(labels = method.name.FDR, values = method.col.FDR) +
  scale_fill_manual(labels = method.name.FDR, values = alpha(method.col.FDR, 0.35)) +
  facet_wrap(~Noise, nrow = 1) +
  labs(x = expression(pi[0]), y = "TDP") +
  expand_limits(y = 0) +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_text(size = 15),
        axis.title.x = element_text(size = 15),
        axis.text.x = element_text(size = 10),
        strip.text = element_text(size = 15),
        legend.position = "bottom",
        legend.background = element_rect(color = "grey"),
        legend.text = element_text(size = 12)
        )

Warning: Ignoring unknown aesthetics: width

FDP.TDP.q.sep.plot.save <- gridExtra::arrangeGrob(z.sep.plot.save, FDP.q.all.sep.plot + theme(plot.title = element_text(size = 15, hjust = 0)), TDP.q.all.sep.plot, heights = c(1, 1.2, 1.1))

ggsave("../output/paper/FDP.TDP.q.sep.pdf", FDP.TDP.q.sep.plot.save, height = 6.5, width = 10)

FDP.sqrt.TDP.q.sep.plot.save <- gridExtra::arrangeGrob(z.sep.plot.save, FDP.sqrt.q.all.sep.plot + theme(plot.title = element_text(size = 15, hjust = 0)), TDP.q.all.sep.plot, heights = c(1, 1.2, 1.1))

FDP.sqrt.TDP.q.sep.plot.save <- gridExtra::arrangeGrob(
  z.sep.plot +
    labs(title = "Typical Examples of Realized Correlated N(0,1) Noise") +
    scale_color_manual(values = "blue") +
    theme(legend.title = element_blank(),
         # legend.title = element_text(size = 15),
         plot.margin = grid::unit(c(5.5, 5.5, 5.5, 67.5), "points"),
         plot.title = element_text(size = 12, hjust = 0.5),
         strip.text = element_text(size = 12),
         axis.title.x = element_blank(),
         axis.text.x = element_text(size = 10),
         axis.title.y = element_text(size = 12)
    ),
  gridExtra::arrangeGrob(
  FDP.sqrt.q.all.sep.plot +
    labs(title = "FDP") +
    theme(plot.title = element_text(size = 12, hjust = 0.5),
          axis.title.x = element_blank(),
          strip.text = element_text(size = 12),
          plot.margin = grid::unit(c(5.5, 5.5, 5.5, 5.5), "points"),
          axis.text.y = element_text(size = 12)
    ),
  TDP.q.all.sep.plot +
    labs(title = "TDP") +
    theme(plot.margin = grid::unit(c(5.5, 5.5, 5.5, 5.5), "points"),
          axis.title.x = element_blank(),
          strip.text = element_text(size = 12),
          plot.title = element_text(size = 12, hjust = 0.5),
          axis.text.y = element_text(size = 12)
    ),
  heights = c(1, 1),
  top = grid::textGrob(bquote(paste("Nominal FDR = ", .(q), " (", g[1], " is Gaussian; ", pi[0] == 0.9, ")")), gp = grid::gpar(fontsize = 12), hjust = 1.15)
  ),
  heights = c(1, 2)
)

ggsave("../output/paper/FDP.sqrt.TDP.q.sep.pdf", FDP.sqrt.TDP.q.sep.plot.save, height = 6, width = 8)

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.6

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] ggplot2_2.2.1

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.16      knitr_1.20        whisker_0.3-2    
 [4] magrittr_1.5      workflowr_1.1.1   munsell_0.4.3    
 [7] colorspace_1.3-2  rlang_0.2.0       stringr_1.3.1    
[10] plyr_1.8.4        tools_3.4.3       grid_3.4.3       
[13] gtable_0.2.0      R.oo_1.22.0       git2r_0.21.0     
[16] htmltools_0.3.6   yaml_2.1.19       lazyeval_0.2.1   
[19] rprojroot_1.3-2   digest_0.6.15     tibble_1.4.2     
[22] gridExtra_2.3     reshape2_1.4.3    R.utils_2.6.0    
[25] evaluate_0.10.1   rmarkdown_1.9     labeling_0.3     
[28] stringi_1.2.2     pillar_1.2.2      compiler_3.4.3   
[31] scales_0.5.0      backports_1.1.2   R.methodsS3_1.7.1