The inflated vs deflated noise figure in cashr

Introduction

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

library(ggplot2)

Warning: package 'ggplot2' was built under R version 3.4.4

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.05, 0.25, 0.5, 0.75, 0.95))
  names(r) <- c("ymin", "lower", "middle", "upper", "ymax")
  return(r)
}

boxplot.quantile.sq <- function (x) {
  r <- sqrt(quantile(x^2, probs = c(0.05, 0.25, 0.5, 0.75, 0.95)))
  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) <- levels(Noise)
z.sep.ggdata <- reshape2::melt(z.list.sel, value.name = "z")
z.sep.ggdata$L1 <- factor(z.sep.ggdata$L1, levels = levels(Noise))
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 = "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)

FDP_noise.save <- gridExtra::arrangeGrob(
  z.sep.plot +
    labs(y = 'Density', title = "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, 110), "points"),
         plot.title = element_text(size = 15),
         strip.text = element_text(size = 15),
         axis.title.x = element_blank(),
         axis.text.x = element_text(size = 12),
         axis.title.y = element_text(size = 15)
    ),
  FDP.sqrt.q.all.sep.plot +
    labs(x = 'FDP', title = bquote(paste('Distribution of empirical FDP at nominal FDR = ', .(q)))) +
    theme(plot.title = element_text(size = 15),
          axis.title.x = element_text(size = 15),
          strip.text = element_text(size = 15),
          plot.margin = grid::unit(c(15, 5.5, 5.5, 5.5), "points"),
          axis.text.y = element_text(size = 15),
          axis.text.x = element_text(size = 12)
    ),
  TDP.q.all.sep.plot +
    labs(x = 'TDP', title = bquote(paste('Distribution of empirical TDP at nominal FDR = ', .(q)))) +
    theme(plot.margin = grid::unit(c(15, 5.5, 5.5, 5.5), "points"),
          axis.title.x = element_text(size = 15),
          strip.text = element_text(size = 15),
          plot.title = element_text(size = 15),
          axis.text.y = element_text(size = 15),
          axis.text.x = element_text(size = 12)
    ),
  heights = c(1.75, 2, 2)
)

ggsave("../output/paper/FDP_noise.pdf", FDP_noise.save, height = 7, width = 10)

FDP.q.ggdata.0.9$Noise <- factor(FDP.q.ggdata.0.9$Noise, levels = levels(FDP.q.ggdata.0.9$Noise)[c(2, 3, 4, 1)])

TDP.q.ggdata.0.9$Noise <- factor(TDP.q.ggdata.0.9$Noise, levels = levels(TDP.q.ggdata.0.9$Noise)[c(2, 3, 4, 1)])

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

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_noise_new.save <- gridExtra::arrangeGrob(
  z.sep.plot +
    labs(y = 'Density', title = '(a): Examples of realized correlated N(0,1) noise') +
    scale_color_manual(values = "blue") +
    expand_limits(x = -4) +
    theme(legend.title = element_blank(),
         legend.position = 'right',
         plot.margin = grid::unit(c(10, 95.5, 10, 17.5), "points"),
         plot.title = element_text(size = 15),
         strip.text = element_text(size = 15),
         axis.title.x = element_blank(),
         axis.text.x = element_text(size = 12),
         axis.title.y = element_text(size = 15)
    ),
  FDP.sqrt.q.all.sep.plot +
    labs(x = 'FDP', title = bquote(paste('(b): Distribution of FDP'))) +
    theme(plot.title = element_text(size = 15),
          axis.title.x = element_text(size = 15),
          strip.text = element_text(size = 15),
          plot.margin = grid::unit(c(5.5, 5.5, 5.5, 5.5), "points"),
          axis.text.y = element_text(size = 15),
          axis.text.x = element_text(size = 12)
    ),
  TDP.q.all.sep.plot +
    labs(x = 'TDP', title = bquote(paste('(c): Distribution of TDP'))) +
    theme(plot.margin = grid::unit(c(5.5, 5.5, 5.5, 5.5), "points"),
          axis.title.x = element_text(size = 15),
          strip.text = element_text(size = 15),
          plot.title = element_text(size = 15),
          axis.text.y = element_text(size = 15),
          axis.text.x = element_text(size = 12)
    ),
  heights = c(1.1, 1, 1),
  top = grid::textGrob(label = bquote(paste(
    'Nominal FDR = ', .(q), ' (', g[1], ' is Guassian; ', pi[0] == 0.9, ')')), gp = grid::gpar(fontsize = 15, fontface = 'bold'), hjust = 1.125)
)

ggsave("../output/paper/FDP_noise_new.pdf", FDP_noise_new.save, height = 7, width = 10)

Session information

sessionInfo()

R version 3.4.3 (2017-11-30)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS  10.14.1

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_3.1.0

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