In mshasan/OPWpaper: Reproduce figures and simulations of OPWeight

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

Below we generate the necessary plots to demonstrate the behavior of the Power of the Covariate Rank Weighting (CRW) method as well as the other methods under the different circumstances.

Load the necessary libraries

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

Continuous effects

This part show the power 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 the 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 the 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"))

Power vs. effect size

The simulation procedures for the power is divided into three groups based on the proportion of the true null hypothesis. Three groups composed of $50\%, 90\%,$ and $99\%$ true null tests. For each group of simulations, we considered the effect sizes ${et=ey,$ or $et \sim N(ey,CV.ey)}$, where $et,ey,$ and $CV$ refers to the mean test-effect, mean covariate-effect, and coefficient of variations. For the mean effect sizes we considered a vector of ${.2, .4, .6, .8, 1, 2, 3, 5, 8}$ and $CV = 1/2$. We compared the results to that of using the BH [@benjamini1997false], RDW [@roeder2009genome], and IHW [@ignatiadis2016natmeth] methods.

# plots of power for the mean covariate-effect(ey) = mean test-effect(et)
#-------------------------------------------------------------------------------
p_.5_eq_power <- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower2e1, null = 50, figure = "effectVsFPFP")
p_.9_eq_power <- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower4e1, null = 90, figure = "effectVsFPFP")
p_.99_eq_power<- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower5e1, null = 99, figure = "effectVsFPFP")

p_.5_low_ef_eq_power <- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower2e1, null = 50, low_eff_plot = TRUE, figure = "effectVsFPFP")
p_.9_low_ef_eq_power <- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower4e1, null = 90, low_eff_plot = TRUE, figure = "effectVsFPFP")
p_.99_low_ef_eq_power<- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower5e1, null = 99, low_eff_plot = TRUE, figure = "effectVsFPFP")

p_eq_power = plot_grid(p_.5_eq_power, p_.9_eq_power, p_.99_eq_power,
                       p_.5_low_ef_eq_power, p_.9_low_ef_eq_power, p_.99_low_ef_eq_power,
                       ncol = 3, labels = letters[1:3], align = 'hv')
title <- ggdraw() + draw_label("Power: et = ey")
plot_grid(title, p_eq_power, legend, ncol = 1, rel_heights=c(.1, 1, .1))

Figure 1: The Power of the four methods when the mean test-effect ($et$) is equal to mean covariate-effect ($ey$). Each plot consists of the four curves of CRW, BH, RDW, and IHW methods. The first row shows the power for the low to high effect sizes, and the second row shows $log(power)$ for the low effect sizes. Three columns represent three groups of $50\%, 90\%$, and $99\%$ true null hypothesis.

# plots of power for
# mean test-effect(et) ~ Normal (mean covariate-effect, mean covariate- effect/2) (i.e cv = 1/2)
#-----------------------------------------------------------------------------------------
p_.5_uneq_power <- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower2e2, null = 50, figure = "effectVsFPFP")
p_.9_uneq_power <- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower4e2, null = 90, figure = "effectVsFPFP")
p_.99_uneq_power<- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower5e2, null = 99, figure = "effectVsFPFP")

p_.5_low_ef_uneq_power <- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower2e2, null = 50, low_eff_plot = TRUE, figure = "effectVsFPFP")
p_.9_low_ef_uneq_power <- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower4e2, null = 90, low_eff_plot = TRUE, figure = "effectVsFPFP")
p_.99_low_ef_uneq_power<- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower5e2, null = 99, low_eff_plot = TRUE, figure = "effectVsFPFP")

p_uneq_power = plot_grid(p_.5_uneq_power, p_.9_uneq_power, p_.99_uneq_power,
                         p_.5_low_ef_uneq_power, p_.9_low_ef_uneq_power, p_.99_low_ef_uneq_power,
                         ncol = 3, labels = letters[1:3], align = 'hv')
title <- ggdraw() + draw_label("Power: et ~ Normal(ey, ey/2)")
plot_grid(title, p_uneq_power, legend, ncol = 1, rel_heights=c(.1, 1, .1))

Figure 2: The power for the same parameters described in Figure 1 except the mean test-effect ($et$) is not equal to the mean covariate-effect ($ey$); rather $et \sim Normal(ey, ey/2)$, where $CV=1/2$.

Power vs. proportion of the true null hypothesis

This section shows the impact of the proportion of the true null hypothesis tests for the various mean covariate-effects and the correlations. For the simulation procedures, we kept the similar set-up described in the Power vs. effect size section except here we considered five groups of true null hypothesis composed of $20\%, 50\%, 75\%, 90\%,$ and $99\%$ true null tests.

# see the influence of the null proportion
#------------------------------------------------
# this code is to load the saved workspace from parallel computing
nullProp <- c(20, 50, 75, 90, 99)

# corr = .3-------------
mat_ef.6<- rbind(FwerPowerFdrPower1f1[13:16, 4], FwerPowerFdrPower2f1[13:16, 4],
                 FwerPowerFdrPower3f1[13:16, 4], FwerPowerFdrPower4f1[13:16, 4],
                 FwerPowerFdrPower5f1[13:16, 4])
p_ef.6 <- nice_plots(x_vec = nullProp, y_matrix = mat_ef.6, ey = 0.6, figure = "nullPropVsPower")


mat_ef1 <- rbind(FwerPowerFdrPower1f1[13:16, 6], FwerPowerFdrPower2f1[13:16, 6],
                 FwerPowerFdrPower3f1[13:16, 6], FwerPowerFdrPower4f1[13:16, 6],
                 FwerPowerFdrPower5f1[13:16, 6])
p_ef1 <- nice_plots(x_vec = nullProp, y_matrix = mat_ef1, ey = 1.0, figure = "nullPropVsPower")


mat_ef3 <- rbind(FwerPowerFdrPower1f1[13:16, 8], FwerPowerFdrPower2f1[13:16, 8],
                 FwerPowerFdrPower3f1[13:16, 8], FwerPowerFdrPower4f1[13:16, 8],
                 FwerPowerFdrPower5f1[13:16, 8])
p_ef3 <- nice_plots(x_vec = nullProp, y_matrix = mat_ef3, ey = 3.0, figure = "nullPropVsPower")


plot_cor.3 <- plot_grid(p_ef.6, p_ef1, p_ef3, align = 'hv', ncol = 3, labels=letters[1:3])
title <- ggdraw() + draw_label("corr = 0.3", fontface='bold')
plot_cor.3_title <- plot_grid(title, plot_cor.3, ncol=1, rel_heights=c(0.1, 1), align = 'hv')



# corr = .7 -------------
mat_ef.6.7 <- rbind(FwerPowerFdrPower1h1[13:16, 4], FwerPowerFdrPower2h1[13:16, 4],
                 FwerPowerFdrPower3h1[13:16, 4], FwerPowerFdrPower4h1[13:16, 4],
                 FwerPowerFdrPower5h1[13:16, 4])
p_ef.6.7 <- nice_plots(x_vec = nullProp, y_matrix = mat_ef.6.7, ey = 0.6, figure = "nullPropVsPower")


mat_ef1.7 <- rbind(FwerPowerFdrPower1h1[13:16, 6], FwerPowerFdrPower2h1[13:16, 6],
                 FwerPowerFdrPower3h1[13:16, 6], FwerPowerFdrPower4h1[13:16, 6],
                 FwerPowerFdrPower5h1[13:16, 6])
p_ef1.7 <- nice_plots(x_vec = nullProp, y_matrix = mat_ef1.7, ey = 1.0, figure = "nullPropVsPower")


mat_ef3.7 <- rbind(FwerPowerFdrPower1h1[13:16, 8], FwerPowerFdrPower2h1[13:16, 8],
                 FwerPowerFdrPower3h1[13:16, 8], FwerPowerFdrPower4h1[13:16, 8],
                 FwerPowerFdrPower5h1[13:16, 8])
p_ef3.7 <- nice_plots(x_vec = nullProp, y_matrix = mat_ef3.7, ey = 3.0, figure = "nullPropVsPower")


plot_cor.7 <- plot_grid(p_ef.6.7, p_ef1.7, p_ef3.7, align = 'hv', ncol = 3, labels=letters[4:6])
title <- ggdraw() + draw_label("corr = 0.7", fontface='bold')
plot_cor.7_title <- plot_grid(title, plot_cor.7, ncol=1, rel_heights=c(0.1, 1), align = 'hv')


# for the main title------------
p_prop = plot_grid(plot_cor.3_title, plot_cor.7_title, ncol = 1, align = 'hv')
title_main <- ggdraw() + draw_label("Continuous: Power vs. prop. of null")
# make sure get legend----------
plot_grid(title_main, p_prop, legend, ncol = 1, rel_heights=c(.1, 1, .1))

**Figure 3:** The power across the proportion of the true null hypothesis for the different test correlations.

Power vs. test correlation

This section shows the impact of the correlations between the test statistics for the various mean effect sizes. We considered the combination of correlations $\times$ effect sizes = ${0, .3, .5, .7, .9} \times {.2, .4, .6, .8, 1, 2, 3, 5, 8}$.

covariateEffectVec <- c(seq(0,1,.2),2,3,5,8)

# use 2 for 50% and 4 for 90% nulls-----------
E = FwerPowerFdrPower2e1
FF = FwerPowerFdrPower2f1
G = FwerPowerFdrPower2g1
H = FwerPowerFdrPower2h1
I = FwerPowerFdrPower2i1

corr = c(0,.3,.5,.7,.9)       # correlations
r = 13                        # row starts for fdr based power

gplots <- list()
for(e in 3:8)               # effect size index
{
    CRW = c(E[r,e],    FF[r,e],    G[r,e],    H[r,e],    I[r,e])
    BH = c(E[(r+1),e],FF[(r+1),e],G[(r+1),e],H[(r+1),e],I[(r+1),e])
    RDW = c(E[(r+2),e],FF[(r+2),e],G[(r+2),e],H[(r+2),e],I[(r+2),e])
    IHW = c(E[(r+3),e],FF[(r+3),e],G[(r+3),e],H[(r+3),e],I[(r+3),e])
    dat = data.frame(corr, CRW, BH, RDW, IHW)
    dat2 = melt(dat, id.var = "corr")
    gplots[[e]] <- ggplot(dat2, aes(x = corr, y = value, group = variable,
                                    col = variable)) +
        geom_line(aes(linetype = variable), size = 1.5) +
        labs(x = "corr", y = "power", title = paste("et = ", covariateEffectVec[e])) +
        #theme(legend.title = element_blank())
        theme(legend.position="none")
}

gplots[[8]] <- gplots[[8]] + theme(legend.position="bottom", legend.title = element_blank())

legend_corr <- get_legend(gplots[[8]])

gplots[[8]] <- gplots[[8]] + theme(legend.position="none")


p = plot_grid(gplots[[3]],gplots[[4]],gplots[[5]],gplots[[6]],gplots[[7]],gplots[[8]],
              ncol = 3, labels = letters[1:6], align = 'hv')
title <- ggdraw() + draw_label("Continuous: null = 50%, et = ey")
plot_grid(title, p, legend_corr, ncol = 1, rel_heights=c(.1, 1, .1))

**Figure 4:** Power across different correlations between the test statistics for the different effect sizes when $50\%$ tests are true null

# 90% nulls-----------
E = FwerPowerFdrPower4e1
FF = FwerPowerFdrPower4f1
G = FwerPowerFdrPower4g1
H = FwerPowerFdrPower4h1
I = FwerPowerFdrPower4i1

corr = c(0,.3,.5,.7,.9)       # correlations
r = 13                        # row starts for fdr based power

gplots <- list()
for(e in 3:8)               # effect size index
{
    CRW = c(E[r,e],    FF[r,e],    G[r,e],    H[r,e],    I[r,e])
    BH = c(E[(r+1),e],FF[(r+1),e],G[(r+1),e],H[(r+1),e],I[(r+1),e])
    RDW = c(E[(r+2),e],FF[(r+2),e],G[(r+2),e],H[(r+2),e],I[(r+2),e])
    IHW = c(E[(r+3),e],FF[(r+3),e],G[(r+3),e],H[(r+3),e],I[(r+3),e])
    dat = data.frame(corr, CRW, BH, RDW, IHW)
    dat2 = melt(dat, id.var = "corr")
    gplots[[e]] <- ggplot(dat2, aes(x = corr, y = value, group = variable,
                                    col = variable)) +
        geom_line(aes(linetype = variable), size = 1.5) +
        labs(x = "corr", y = "power", title = paste("et = ", covariateEffectVec[e])) +
        #theme(legend.title = element_blank())
        theme(legend.position="none")
}

gplots[[8]] <- gplots[[8]] + theme(legend.position="bottom", legend.title = element_blank())

legend_corr <- get_legend(gplots[[8]])

gplots[[8]] <- gplots[[8]] + theme(legend.position="none")


p = plot_grid(gplots[[3]],gplots[[4]],gplots[[5]],gplots[[6]],gplots[[7]],gplots[[8]],
              ncol = 3, labels = letters[1:6], align = 'hv')
title <- ggdraw() + draw_label("Continuous: null = 90%, et = ey")
plot_grid(title, p, legend_corr, ncol = 1, rel_heights=c(.1, 1, .1))

**Figure 5:** Power across different correlations between the test statistics for the different effect sizes when $90\%$ tests are true null.

Power vs. variance of the test effect

This section shows the impact of the variation of the mean test-effect ($et$) for the different mean covariate-effect ($ey$) via simulation. For the simulation procedures, we kept the similar set-up described in the Power vs. effect size section except here we considered one group of the true null proportion at a time for the different coefficient of variations ($CV$). That is, we considered that the mean test-effect $et \sim N(ey, CV.ey)$ and $CV={1,3,10}$.

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

# plots of power for miss variance of the test effect size; et ~ normal(ey, CV*ey)
# CV = coefficient of variance (i.e cv = 1, 3, 10)
# 50% null case
#----------------------------------------------------------------------------
p_cv1 <- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower2a2, cv = 1, figure = "CV")
p_cv3 <- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower2a3, cv = 3, figure = "CV")
p_cv10<- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower2a10,cv = 10,figure = "CV")

p_cv1_low_ef <- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower2a2, cv = 1, low_eff_plot = TRUE, figure = "CV")
p_cv3_low_ef <- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower2a3, cv = 3, low_eff_plot = TRUE, figure = "CV")
p_cv10_low_ef<- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower2a10,cv = 10,low_eff_plot = TRUE, figure = "CV")

p_cv_power_50 = plot_grid(p_cv1, p_cv3, p_cv10, p_cv1_low_ef, p_cv3_low_ef, p_cv10_low_ef,
                          ncol = 3, labels = letters[1:3], align = 'hv')
title <- ggdraw() + draw_label("Power: et ~ Normal(ey, CV*ey), null = 50%")
plot_grid(title, p_cv_power_50, legend, ncol = 1, rel_heights=c(.1, 1, .1))

Figure 6: The power when the test effect ($et$) of the true alternative hypothesis is not equal to mean covariate effect ($ey$); rather $et \sim Normal(ey, CV.ey)$ if $50\%$ tests are from the true null models. The first row shows the power for low to high effect sizes, and the second row shows power for the low effect sizes. Three columns are based on three groups composed of $CV = 1,3,$ or $10$.

# CV = coefficient of variance (i.e cv = 1, 3, 10)
# 90% null case------------
#--------------------------------------------------------------------
p_cv1 <- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower4a2, cv = 1, figure = "CV")
p_cv3 <- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower4a3, cv = 3, figure = "CV")
p_cv10<- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower4a10,cv = 10,figure = "CV")

p_cv1_low_ef <- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower4a2, cv = 1, low_eff_plot = TRUE, figure = "CV")
p_cv3_low_ef <- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower4a3, cv = 3, low_eff_plot = TRUE, figure = "CV")
p_cv10_low_ef<- nice_plots(x_vec = ey_vec, y_matrix = FwerPowerFdrPower4a10,cv = 10,low_eff_plot = TRUE, figure = "CV")

p_cv_power_90 = plot_grid(p_cv1, p_cv3, p_cv10, p_cv1_low_ef, p_cv3_low_ef, p_cv10_low_ef,
                          ncol = 3, labels = letters[1:3], align = 'hv')
title <- ggdraw() + draw_label("Power: et ~ Normal(ey, CV*ey), null = 90%")
plot_grid(title, p_cv_power_90, legend, ncol = 1, rel_heights=c(.1, 1, .1))