Skip to content

Instantly share code, notes, and snippets.

@agisga
Created November 30, 2016 19:58
Show Gist options
  • Save agisga/db980b5c7e2d5676d7b1911d8a404222 to your computer and use it in GitHub Desktop.
Save agisga/db980b5c7e2d5676d7b1911d8a404222 to your computer and use it in GitHub Desktop.
An iteractive visualization comparing 6 methods for multiple test corrections
#---
#
# 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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment