In mshasan/OPWpaper1: Reproduce figures and simulations of OPWeight

knitr::opts_chunk$set(fig.width = 8.5, fig.height = 4)
knitr::opts_chunk$set(tidy = FALSE, cache = FALSE, autodep = TRUE)

Below we generate the necessary plots to demonstrate that the Covariate Rank Weighting (CRW) method controls the Family Wise Error Rate (FWER) and the False Discovery Rate (FDR).

Load the necessary libraries

library(OPWpaper)       
library(ggplot2)
library(reshape2)       # library for the melt function
library(cowplot)        # plot_grid function

Load the data stored in r Biocpkg("OPWpaper")

fwer_dat <- system.file("simulations/results", package = "OPWpaper")
setwd(fwer_dat)
load("simu_fwer.RDATA")
#load(system.file("simulations/results", "simu_fwer.RDATA",
#                 package = "OPWpaper"), envir = environment())

FWER when all tests are from the true null models

This simulation is conducted to verify that the proposed method controls the Family Wise Error Rate (FWER). To conduct this simulation, we assumed that all test statistics are from the true null models.

fwer_by_alpha <- matrix(rowMeans(fwer_mat, na.rm = TRUE), nrow = 4, byrow = FALSE)

alphaVal = seq(.01, .1, .02)
datError <- data.frame(alphaVal, t(fwer_by_alpha))
colnames(datError) <- c("alpha","BON","CRW_bin","CRW_cont", "IHW")
datError2 <- melt(datError, id.var="alpha")

ggplot(datError2, aes(x = alpha, y = value, col=variable)) +
    geom_line(size=1.5) +
    geom_abline(linetype="dashed") +
    xlab(expression(bold(paste("Nominal ",alpha)))) +
    ylab("FWER")+
    scale_x_continuous(limits = c(0.01,0.1), breaks=seq(0.01,0.09,length=5)) +
    #ylim(0,0.9) +
    theme(legend.title = element_blank())+
    theme(axis.title = element_text(face="bold"))+
    theme(panel.background = element_rect(fill = 'white', colour = 'black'))

Figure 1: The FWER for the different significance levels of $\alpha$. In the legend, the representation is: BON = bonferroni, CRW_bin = CRW binary, CRW_cont = CRW continuous, and IHW = Independent Hypothesis Weighting [@ignatiadis2016natmeth].

Continuous effects

This part show the FWER and FDR control when the effect sizes are continuous. We only show the results for the continuous effect sizes. For the binary effects, one just needs to load $.Rdata$ for the binary effects and then apply the same R-code.

Load data stored in r Biocpkg("OPWpaper")

load(system.file("simulations/results", "simu_fwerPowerFdrPower_cont.RDATA",
                 package = "OPWpaper"), envir = environment())

Extract legend to use for the following plots

# this part is for legend------------------------------------------------------
ey_vec <- c(seq(0, 1, .2), 2, 3, 5, 8)

dat_99 <- data.frame(ey_vec, t(FwerPowerFdrPower5e1[13:16,]))
colnames(dat_99) <- c("effectSize", "CRW", "BH", "RDW", "IHW")
dat_99_par <- melt(dat_99[1:6,], id.var = "effectSize")

p_99_par <- ggplot(dat_99_par, aes(x = effectSize, y = value,group = variable, col=variable)) +
    geom_line(aes(linetype = variable), size = 1.5) +
    labs(x = "ey", y = "power", title = "null = 99%") +
    theme(legend.title = element_blank())

legend <- get_legend(p_99_par + theme(legend.direction = "horizontal", legend.position = "bottom"))

FWER for the different effect sizes

This simulation is conducted to see the FWER control for the different effect sizes. The simulation procedures are divided into three groups based on the proportion of the true null hypothesis. Three groups are composed of $50\%,90\%,$ and $99\%$ true null tests. For each group of simulations, we considered mean effect sizes ${et=ey$ or $et \sim Normal(ey,CV.ey)}$, where $et,ey,$ and $CV$ refers to the mean test-effect, mean covariate-effect, and coefficient of variations. For the mean effects we considered a vector of ${.2,.4,.6,.8,1,2,3,5,8}$ and $CV = 1/2$. We compare the Covariate Rank Weighting (CRW) method with three other methods which are Bonferroni (BON), Roeder and Wasserman (RDW) [@roeder2009genome], and Independent Hypotheis Weighting (IHW).

Plot to observe the Family Wise Error Rate (FWER) for the different effect sizes

# plots FWER et = ey (i.e cv =0)
#-------------------------------------------------
p_.5_eq_fwer <- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower2e1, fdr = FALSE, power = FALSE, null = 50, figure = "effectVsFPFP")
p_.9_eq_fwer <- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower4e1, fdr = FALSE, power = FALSE, null = 90, figure = "effectVsFPFP")
p_.99_eq_fwer<- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower5e1, fdr = FALSE, power = FALSE, null = 99, figure = "effectVsFPFP")

p_eq_fwer = plot_grid(p_.5_eq_fwer, p_.9_eq_fwer, p_.99_eq_fwer, ncol = 3, labels = letters[1:3], align = 'hv')
title <- ggdraw() + draw_label(expression(paste("FWER: et = ey, ", alpha, " = .05")))
plot_grid(title, p_eq_fwer, legend, ncol = 1, rel_heights=c(.1, .5, .1))

Figure 2: The FWER of four methods for the different mean effects when the mean test-effect ($et$) is equal to mean covariate-effect ($ey$). Three columns are based on three groups composed of $50\%,90\%$, and $99\%$ true null hypothesis.

# plots FWER et ~ Normal(ey, ey/2) (i.e cv = 1/2)
#-------------------------------------------------
p_.5_uneq_fwer <- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower2e2, fdr = FALSE, power = FALSE, null = 50, figure = "effectVsFPFP")
p_.9_uneq_fwer <- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower4e2, fdr = FALSE, power = FALSE, null = 90, figure = "effectVsFPFP")
p_.99_uneq_fwer<- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower5e2, fdr = FALSE, power = FALSE, null = 99, figure = "effectVsFPFP")

p_uneq_fwer = plot_grid(p_.5_uneq_fwer, p_.9_uneq_fwer, p_.99_uneq_fwer, ncol = 3, labels = letters[1:3], align = 'hv')
title <- ggdraw() + draw_label(expression(paste("FWER: et ~ Normal(ey, ey/2), ", alpha, " = .05")))
plot_grid(title, p_uneq_fwer, legend, ncol = 1, rel_heights=c(.1, .5, .1))

Figure 3: The FWER when the mean test effect ($et$) is not equal to mean covariate effect ($ey$); rather $et \sim Normal(ey, CV.ey )$, where $CV =1/2$.

FDR for the different effect sizes

Plot to observe the False Discovery Rate (FDR) for the different effect sizes. For this simulation, we kept the same set-up described in the FWER section. We compare the CRW method with the three other methods which are Benjamini and Hochberg (BH) [@benjamini1997false], Roeder and Wasserman (RDW), and Independent Hypotheis Weighting (IHW).

# plots FDR et = ey (i.e cv =0)
#-------------------------------------------------
p_.5_eq_fdr <- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower2e1, fdr = TRUE, power = FALSE, null = 50, figure = "effectVsFPFP")
p_.9_eq_fdr <- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower4e1, fdr = TRUE, power = FALSE, null = 90, figure = "effectVsFPFP")
p_.99_eq_fdr<- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower5e1, fdr = TRUE, power = FALSE, null = 99, figure = "effectVsFPFP")

p_eq_fdr = plot_grid(p_.5_eq_fdr, p_.9_eq_fdr, p_.99_eq_fdr, ncol = 3, labels = letters[1:3], align = 'hv')
title <- ggdraw() + draw_label(expression(paste("FDR: et = ey, ", alpha, " = .05")))
plot_grid(title, p_eq_fdr, legend, ncol = 1, rel_heights=c(.1, .5, .1))

**Figure 4:** The FDR when the mean test effect ($et$) is equal to mean covariate effect ($ey$)$.

# plots FDR et ~ Normal(ey, ey/2) (i.e cv = 1/2)
#-------------------------------------------------
p_.5_uneq_fdr <- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower2e2, fdr = TRUE, power = FALSE, null = 50, figure = "effectVsFPFP")
p_.9_uneq_fdr <- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower4e2, fdr = TRUE, power = FALSE, null = 90, figure = "effectVsFPFP")
p_.99_uneq_fdr<- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower5e2, fdr = TRUE, power = FALSE, null = 99, figure = "effectVsFPFP")

p_uneq_fdr = plot_grid(p_.5_uneq_fdr, p_.9_uneq_fdr, p_.99_uneq_fdr, ncol = 3, labels = letters[1:3], align = 'hv')
title <- ggdraw() + draw_label(expression(paste("FDR: et ~ Normal(ey, ey/2), ", alpha, " = .05")))
plot_grid(title, p_uneq_fdr, legend, ncol = 1, rel_heights=c(.1, .5, .1))

Figure 5: The FDR when the mean test effect ($et$) is not equal to mean covariate effect ($ey$); rather $et \sim Normal(ey, CV.ey )$, where $CV =1/2$.