In mshasan/OPWpaper: Reproduce figures and simulations of OPWeight

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

Below we generate the necessary plots to demonstrate the ranks probability for both the continuous effects and the binary effects.

Load necessary libraries

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

Ranks proabbility comparison

Load data stored in r Biocpkg("OPWpaper")

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

To verify the proposed ranks probability, we conducted simulations. Figure 1 shows the ranks probability of a test from three different approaches: 1) Simulation approach and 2) Exact numerical solution and 3) Normal approximation of the proposed method. All plots of the simulation suggest nearly perfect alignment with the proposed (exact and approximate) methods.

dat00 <- prob_50_0_cont[[1]]
colnames(dat00) <- c("ranks", "SH0","SH1","EH0","EH1","AH0","AH1")
dat00 <- melt(dat00, id.var="ranks")

dat01 <- prob_50_1_cont[[2]]
colnames(dat01) <- c("ranks", "SH0","SH1","EH0","EH1","AH0","AH1")
dat01 <- melt(dat01, id.var="ranks")

dat12 <- prob_50_1_cont[[3]]
colnames(dat12) <- c("ranks", "SH0","SH1","EH0","EH1","AH0","AH1")
dat12 <- melt(dat12, id.var="ranks")

dat02 <- prob_50_2_cont[[2]]
colnames(dat02) <- c("ranks", "SH0","SH1","EH0","EH1","AH0","AH1")
dat02 <- melt(dat02,id.var="ranks")


p00 <- ggplot(dat00, aes(x = ranks, y = value, group = variable,
                       colour = variable)) +
          geom_line(aes(linetype = variable), size = 1.5) +
          labs(x = "Ranks", y = "p(rank | effect)", title = "ey = 0, e.one = 0") +
          theme(legend.position="none") +
          annotate("text", x=50, y=.011, label="P(rank | effect = 0)")


p01 <- ggplot(dat01, aes(x = ranks, y = value, group = variable,
                       colour = variable)) +
          geom_line(aes(linetype = variable), size = 1.5) +
          labs(x = "Ranks", y = "p(rank | effect)", title = "ey ~ U(0, 1), e.one = 1") +
          theme(legend.position="none") +
          annotate("text", x = c(50, 50), y = c(.005, .018),
                   label = c(paste(sprintf('\u2190'),"P(rank | effect = 0)"),
                    paste("P(rank | effect = e.one)", sprintf('\u2192'))))

p12 <- ggplot(dat12, aes(x = ranks, y = value, group = variable,
                       colour = variable)) +
          geom_line(aes(linetype = variable), size = 1.5) +
          labs(x = "Ranks", y = "p(rank | effect)", title = "ey ~ U(1, 2), e.one = 1") +
          theme(legend.position="none") +
          annotate("text", x = c(60, 40), y=c(.001, .02),
                   label = c(paste(sprintf('\u2190'),"P(rank | effect = 0)"),
                    paste("P(rank | effect = e.one)", sprintf('\u2192'))))

p02 <- ggplot(dat02, aes(x = ranks, y = value, group = variable,
                       colour = variable)) +
          geom_line(aes(linetype = variable), size = 1.5) +
          labs(x = "Ranks", y = "p(rank | effect)", title = "ey ~ U(0, 1), e.one = 2") +
          theme(legend.title = element_blank(), legend.position="bottom") +
          annotate("text", x = c(50, 40), y=c(.04, .15),
                  label = c("P(rank | effect = 0)",paste(sprintf('\u2190'),"P(rank | effect = e.one)")))


# extract the legend from one of the plots
legend_art <- get_legend(p02 + theme(legend.direction="horizontal",
                                 legend.position="bottom"))

p02 = p02 +    theme(legend.position="none")


p = plot_grid(p00, p01, p12, p02, nrow = 2, labels = letters[1:4], align = 'hv')
# now add the title
title <- ggdraw() + draw_label("Continuous: m0 = 50, m1 = 50")
plot_grid(title, p, legend_art, ncol = 1, rel_heights=c(0.1, 1, .1))

Figure 1: Ranks probability, $p(r_i \mid \varepsilon_i)$, of the continuous effects for three different scenarios. For these plots, we assumed that there are $m = 100$ tests of which $m_0 = 50$ are true nulls and $m_1 = 50$ are true alternatives. Each plot consists of six curves, SH0, SH1, EH0, EH1, AH0, and AH1, where the first letter represents the method (S = simulated, E = exact, and A = approximate), and H0 and H1 represent the hypothesis type.

Ranks probability of the continuous effects

To perform the simulations, we generated $10,000,000$ samples; each sample is composed of $m = 100$ test statistics from the normal distribution for a particular combination of the effect sizes ${0}\times Uniform(a,b)$, where ${0}$ represents the effect size of the null model, and $Uniform(a,b)$ represents the effect sizes of the alternative models. We consider the combination of $a = {0,1 }$ and $b = {1,2}$, respectively. For each sample, we assumed that there are $m=100$ hypotheses tests of which $m_0={20,50,75,90,99}$ tests are from the null, and $m_1={80,50,25,10,1}$ tests are from the alternative models, respectively. For example, if we are interested in testing the hypothesis $H_0: \varepsilon_i = 0$ vs. $H_1: \varepsilon_i > 0; \varepsilon_i \sim uniform(0,1); i = 1,...,100$, and we assumed that $50$ tests are from the null, and $50$ tests are from the alternative models, then the simulation parameters will be $a=0,b=1,m=100,m_0=50$ and $m_1=50$. That is, $50$ test statistics are from the normal distribution (null model) with mean $0$ (effect size 0) and standard deviation $1$, and the remaining $50$ test statistics are from the normal distribution (alternative model) with mean $uniorm(0,1)$ and standard deviation $1$.

The following plots shows the ranks probabilities for the different combination of the true null and the alternative hypotheses.

nullSize <- c(20, 50, 75, 90, 99)    # proportion of the true null tests
lapply(nullSize, ranksProb_compare_plots, effectType = "continuous")

Ranks probability of the binary effects

Consider the similar set-up described for the continuous effects case except the effect sizes combinations. For the binary case, we considered the combination of the effect sizes ${0} \times {0,1,2}$, where ${0}$ represent the effect size of the null model and ${0,1,2}$ represents the effect sizes of the alternative models. Note that, for the binary case the alternative effects of the tests are the same.

nullSize <- c(20, 50, 75, 90, 99)
lapply(nullSize, ranksProb_compare_plots, effectType = "binary")

mshasan/OPWpaper documentation built on March 3, 2021, 7:02 a.m.