An iteractive visualization comparing 6 methods for multiple test corrections

ggvis_multiple_hypotheses_testing.R

#---
#
#    Copyright (C) 2016 Alexej Gossmann
#
#    This program is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program.  If not, see <http://www.gnu.org/licenses/>.
#
#---

#---
# This code uses ggvis to produce an interactive visualization of several
# methods for multiple test corrections. The available methods are:
# Holm, Hochberg, Hommel, Bonferroni, Benjamini-Hochberg, 
# Benjamini-Yekutieli, and None.
#---

library(ggvis)
library(dplyr)
library(tidyr)

set.seed(20161129)
# generate an i.i.d. Gaussian sample
n.samples <- 100
x <- rnorm(n.samples)
# make the first n.signif variables significant
n.signif <- 20
x[1:n.signif] <- x[1:n.signif] + 2.5
# get p-values for the test that mean=0
p <- 2*pnorm(abs(x), lower.tail = FALSE)
# put everything together in a data frame
pvalues <- data.frame(significant = c(rep(TRUE, n.signif), 
                                      rep(FALSE, n.samples - n.signif)),
                      pvalues = p) %>% 
  tbl_df %>% arrange(pvalues) %>% mutate(index = 1:n.samples)

# apply mulitple correction methods
correction.methods <- p.adjust.methods[p.adjust.methods != "fdr"]
pvalues.adj <- sapply(correction.methods, 
                      function(meth) p.adjust(pvalues$pvalues, meth)) %>%
  data.frame() %>% tbl_df() %>% mutate(index = 1:n.samples)
pvalues <- left_join(pvalues, pvalues.adj)

# assign TP, FP, TN, FN labels
get_status <- function(p, significant, threshold = 0.05) {
  result <- rep(NA, length(p))
  for (i in 1:length(p)) {
    if (p[i] <= threshold & significant[i]) {
      result[i] <- "TP"
    } else if (p[i] <= threshold & !significant[i]) {
      result[i] <- "FP"
    } else if (p[i] > threshold & significant[i]) {
      result[i] <- "FN"
    } else {
      result[i] <- "TN"
    }
  }
  return(result)
}
pvalues.result <- transmute(pvalues, 
                            holm = get_status(holm, significant), 
                            hochberg = get_status(hochberg, significant),
                            hommel = get_status(hommel, significant),
                            bonferroni = get_status(bonferroni, 
                                                    significant),
                            BH = get_status(BH, significant),
                            BY = get_status(BY, significant),
                            none = get_status(none, significant),
                            pvalues, index, significant)

# ggvis
map_to_method <- function(x) {
  methods.list <- list("Holm" = pvalues.result$holm, 
                       "Hochberg" = pvalues.result$hochberg, 
                       "Hommel" = pvalues.result$hommel,
                       "Bonferroni" = pvalues.result$bonferroni,
                       "Benjamini-Hochberg" = pvalues.result$BH,
                       "Benjamini-Yekutieli" = pvalues.result$BY,
                       "None" = pvalues.result$none)
  return(methods.list[[x]])
}

pvalues.result %>%
  ggvis(x = ~index, y = ~pvalues, 
        fill = input_radiobuttons(choices = c("Holm", "Hochberg", 
                                        "Hommel", "Bonferroni", 
                                        "Benjamini-Hochberg", 
                                        "Benjamini-Yekutieli", 
                                        "None"), 
                            selected = "None",
                            map = map_to_method,
                            label = "Multiple Correction Method")) %>%
  layer_points() %>%
  add_legend("fill", title = "Outcome",
             properties = legend_props(labels = list(fontSize = 16))) %>%
  add_axis("x", title = "Index") %>%
  add_axis("y", title = "p-values") %>%
  set_options(height = 400, width = 900)