## ----warning=F, message=F-------------------------------------------------- library("ggplot2") library("grid") library("dplyr") library("tidyr") library("cowplot") library("IHWpaper") ## -------------------------------------------------------------------------- # color from https://github.com/dill/beyonce # beyonce::beyonce_palette(30,5) colors <- c("#040320", "#352140", "#871951", "#EB4A60", "#CFAB7A") ## -------------------------------------------------------------------------- sim_folder <- system.file("simulation_benchmarks/theory_paper", package = "IHWpaper") sim_df <- readRDS(file.path(sim_folder, "ihwc_wasserman_normal_sim.Rds")) ## -------------------------------------------------------------------------- taus_df <- expand.grid(fdr_pars = unique(sim_df$fdr_pars), fdr_method = c("BH","IHW","IHW-Storey"), sim_pars = unique(sim_df$sim_pars)) %>% na.exclude() df1 <- left_join(taus_df, select(sim_df, -fdr_pars), by=c("fdr_method","sim_pars")) %>% select(-fdr_pars.y) %>% rename(fdr_pars = fdr_pars.x) df <- full_join(df1, filter(sim_df, fdr_method %in% c("IHWc_Storey","IHWc"))) %>% mutate( fdr_method = ifelse(fdr_method=="IHWc_Storey","IHWc-Storey", fdr_method)) ## -------------------------------------------------------------------------- pi0s <- sapply(strsplit(df$sim_pars, split="[[:punct:]]"), function(x) paste0(x[[2]],".",x[[3]])) df$pi0s <- pi0s ## -------------------------------------------------------------------------- fdr_plot <- ggplot(df, aes(x=fdr_pars, y=FDR, col=fdr_method)) + geom_line() + facet_grid(.~pi0s, labeller=label_bquote(cols=pi[0] == .(pi0s))) + scale_x_log10() + scale_color_manual(values=colors) + theme(legend.title=element_blank()) + xlab(expression(tau~"(censoring threshold)")) + ylab("FDR") fdr_plot ## -------------------------------------------------------------------------- power_plot <- ggplot(df, aes(x=fdr_pars, y=power, col=fdr_method)) + geom_line() + facet_grid(.~pi0s, labeller=label_bquote(cols=pi[0] == .(pi0s))) + scale_x_log10() + scale_color_manual(values=colors) + theme(legend.title=element_blank()) + xlab(expression(tau~"(censoring threshold)")) + ylab("Power") power_plot