In stephaniehicks/benchmarkfdrData2019: Data and Benchmarking Results from Korthauer and Kimes et al. (2019)

knitr::opts_chunk$set(echo = TRUE)

Summary

In this set of simulations, we only consider settings where all tests are null.

Workspace Setup

library(dplyr)
library(ggplot2)
library(SummarizedBenchmark)
library(parallel)

## load helper functions
for (f in list.files("../R", "\\.(r|R)$", full.names = TRUE)) {
    source(f)
}

## project data/results folders
resdir <- "results"
dir.create(resdir, showWarnings = FALSE, recursive = TRUE)

## intermediary files we create below
gauss_file <- file.path(resdir, "null-benchmark-gaussian.rds")
tdist_file <- file.path(resdir, "null-benchmark-t5.rds")
tdist11_file <- file.path(resdir, "null-benchmark-t11.rds")
chisq_file <- file.path(resdir, "null-benchmark-chisq4.rds")

## number of cores for parallelization
cores <- 20
B <- 100

## define bechmarking design
bd <- initializeBenchDesign()

We include Scott's FDR Regression in the analysis. We include both nulltype = "empirical" and nulltype = "theoretical". In both cases, test statistics are modeled as being sampled from mean-shifted normal distributions with equal variance under the null and alternative. While nulltype = "theoretical" makes the additional assumption that the null is standard normal, nulltype = "empirical" attempts to estimate the mean and variance under the null empirically. While the assumptions of nulltype = "theoretical" are met for a subset of the simulations, we include nulltype = "empirical" as this is not true for all simulations, and not necessarily reflective of all real data sets.

bdplus <- bd
bdplus <- addBMethod(bdplus, "fdrreg-t",
                     FDRreg::FDRreg,
                     function(x) { x$FDR },
                     z = test_statistic,
                     features = model.matrix( ~  splines::bs(ind_covariate, df = 3) - 1),
                     nulltype = 'theoretical',
                     control = list(lambda = 0.01))
bdplus <- addBMethod(bdplus, "fdrreg-e",
                     FDRreg::FDRreg,
                     function(x) { x$FDR },
                     z = test_statistic,
                     features = model.matrix( ~  splines::bs(ind_covariate, df = 3) - 1),
                     nulltype = 'empirical',
                     control = list(lambda = 0.01))

All simulation settings will share the following parameters.

m <- 20000                        # integer: number of hypothesis tests
pi0 <- 1                          # numeric: proportion of null hypotheses
es_dist <- function(x) { 0 }        # functional: dist of alternative test stats
icovariate <- runif               # functional: independent covariate

Simulation results will be presented excluding a subset of methods, and for certain plots (upset plots), a single alpha cutoff will be used.

excludeSet <- c("unadjusted", "bl-df02", "bl-df04", "bl-df05")
ualpha <- 0.05

Gaussian Setting

First, we consider the null setting with gaussian test statistics.

Data Simulation

ts_dist <- rnorm_perturber(1)  # functional: sampling dist/noise for test stats
null_dist <- rnorm_2pvaluer(1)    # functional: dist to calc p-values
seed <- 608

We next run the simulations (including Scott's FDR Regression).

if (file.exists(gauss_file)) {
    res <- readRDS(gauss_file)
} else {
    res <- mclapply(X = 1:B, FUN = simIteration, bench = bdplus, m = m,
                    pi0 = pi0, es_dist = es_dist, icovariate = icovariate,
                    ts_dist = ts_dist, null_dist = null_dist,
                    seed = seed, mc.cores = cores)
    saveRDS(res, file = gauss_file)
}
res_i <- lapply(res, `[[`, "informative")
res_u <- lapply(res, `[[`, "uninformative")

Covariate Diagnostics

Here, we show the relationship between the independent covariate and p-values for a single replication of the experiment.

onerun <- simIteration(bdplus, m = m, pi0 = pi0, es_dist = es_dist, ts_dist = ts_dist,
                       icovariate = icovariate, null_dist = null_dist, execute = FALSE)

rank_scatter(onerun, pvalue = "pval", covariate = "ind_covariate")

strat_hist(onerun, pvalue = "pval", covariate = "ind_covariate", maxy = 5, numQ = 3)

Benchmark Metrics

We plot the average results across r B simulations.

resdf <- plotsim_standardize(res_i, alpha = seq(0.01, 0.10, 0.01))

plotsim_average(resdf, met="rejections", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE) 

plotsim_average(resdf, met="FDR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE) 

We also take a look at the distribution of rejects for each method as a function of the effect size and independent covariate.

covariateLinePlot(res_i, alpha = 0.05, covname = "effect_size", trans = "log1p")

covariateLinePlot(res_i, alpha = 0.05, covname = "ind_covariate", trans = "log1p")

Finally, (if enough methods produce rejections at r ualpha) we take a look at the overlap of rejections between methods.

if (numberMethodsReject(resdf, alphacutoff = ualpha, filterSet = excludeSet) >= 3) {
    aggupset(res_i, alpha = ualpha, supplementary = FALSE, return_list = FALSE)
} else {
    message("Not enough methods found rejections at alpha ", ualpha, 
            "; skipping upset plot")
}

We also compare the simulation results with and without an informative covariate.

resdfu <- plotsim_standardize(res_u, alpha = seq(0.01, 0.10, 0.01))

resdfiu <- dplyr::full_join(select(resdf, rep, blabel, param.alpha, key,
                                   performanceMetric, alpha, value),
                            select(resdfu, rep, blabel, param.alpha, key,
                                   performanceMetric, alpha, value),
                            by = c("rep", "blabel", "param.alpha", "key",
                                   "performanceMetric", "alpha"),
                            suffix = c(".info", ".uninfo"))
resdfiu <- dplyr::mutate(resdfiu, value = value.info - value.uninfo)

plotsim_average(resdfiu, met="rejections", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE)

plotsim_average(resdfiu, met="FDR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE)

plotsim_average(resdfiu, met="TPR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE)

plotsim_average(resdfiu, met="TNR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE)

Student's t Setting (df = 5)

Next, we consider the null setting with t-distributed test statistics.

Data Simulation

ts_dist <- rt_perturber(5)  # functional: sampling dist/noise for test stats
null_dist <- rt_2pvaluer(5)    # functional: dist to calc p-values
seed <- 815

For the t-distributed setting, we must specify the number of degrees of freedom for ASH. We add an additional parameter to the ashq method with the corresponding degrees of freedom of the test statistic distribution.

bdplust <- modifyBMethod(bdplus, "ashq", df = 5)

We next run the simulations (including Scott's FDR Regression and ASH with degrees of freedom specified).

if (file.exists(tdist_file)) {
    res <- readRDS(tdist_file)
} else {
    res <- mclapply(X = 1:B, FUN = simIteration, bench = bdplust, m = m,
                    pi0 = pi0, es_dist = es_dist, icovariate = icovariate,
                    ts_dist = ts_dist, null_dist = null_dist,
                    seed = seed, mc.cores = cores)
    saveRDS(res, file = tdist_file)
}
res_i <- lapply(res, `[[`, "informative")
res_u <- lapply(res, `[[`, "uninformative")

Covariate Diagnostics

Here, we show the relationship between the independent covariate and p-values for a single replication of the experiment.

onerun <- simIteration(bdplus, m = m, pi0 = pi0, es_dist = es_dist, ts_dist = ts_dist,
                       icovariate = icovariate, null_dist = null_dist, execute = FALSE)

rank_scatter(onerun, pvalue = "pval", covariate = "ind_covariate")

strat_hist(onerun, pvalue = "pval", covariate = "ind_covariate", maxy = 5, numQ = 3)

Benchmark Metrics

We plot the average results across r B simulations.

resdf <- plotsim_standardize(res_i, alpha = seq(0.01, 0.10, 0.01))

plotsim_average(resdf, met="rejections", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE) 

plotsim_average(resdf, met="FDR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE) 

We also take a look at the distribution of rejects for each method as a function of the effect size and independent covariate.

covariateLinePlot(res_i, alpha = 0.05, covname = "effect_size", trans = "log1p")

covariateLinePlot(res_i, alpha = 0.05, covname = "ind_covariate", trans = "log1p")

Finally, (if enough methods produce rejections at r ualpha) we take a look at the overlap of rejections between methods.

if (numberMethodsReject(resdf, alphacutoff = ualpha, filterSet = excludeSet) >= 3) {
    aggupset(res_i, alpha = ualpha, supplementary = FALSE, return_list = FALSE) 
} else {
    message("Not enough methods found rejections at alpha ", ualpha, 
            "; skipping upset plot")
}

We also compare the simulation results with and without an informative covariate.

resdfu <- plotsim_standardize(res_u, alpha = seq(0.01, 0.10, 0.01))

resdfiu <- dplyr::full_join(select(resdf, rep, blabel, param.alpha, key,
                                   performanceMetric, alpha, value),
                            select(resdfu, rep, blabel, param.alpha, key,
                                   performanceMetric, alpha, value),
                            by = c("rep", "blabel", "param.alpha", "key",
                                   "performanceMetric", "alpha"),
                            suffix = c(".info", ".uninfo"))
resdfiu <- dplyr::mutate(resdfiu, value = value.info - value.uninfo)

plotsim_average(resdfiu, met="rejections", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE)

plotsim_average(resdfiu, met="FDR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE)

plotsim_average(resdfiu, met="TPR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE)

plotsim_average(resdfiu, met="TNR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE)

Student's t Setting (df = 11)

Next, we consider the null setting with t-distributed test statistics.

Data Simulation

ts_dist <- rt_perturber(11)  # functional: sampling dist/noise for test stats
null_dist <- rt_2pvaluer(11)    # functional: dist to calc p-values
seed <- 9158

For the t-distributed setting, we must specify the number of degrees of freedom for ASH. We add an additional parameter to the ashq method with the corresponding degrees of freedom of the test statistic distribution.

bdplust <- modifyBMethod(bdplus, "ashq", df = 11)

We next run the simulations (including Scott's FDR Regression and ASH with degrees of freedom specified).

if (file.exists(tdist11_file)) {
    res <- readRDS(tdist11_file)
} else {
    res <- mclapply(X = 1:B, FUN = simIteration, bench = bdplust, m = m,
                    pi0 = pi0, es_dist = es_dist, icovariate = icovariate,
                    ts_dist = ts_dist, null_dist = null_dist,
                    seed = seed, mc.cores = cores)
    saveRDS(res, file = tdist11_file)
}
res_i <- lapply(res, `[[`, "informative")
res_u <- lapply(res, `[[`, "uninformative")

Covariate Diagnostics

Here, we show the relationship between the independent covariate and p-values for a single replication of the experiment.

onerun <- simIteration(bdplus, m = m, pi0 = pi0, es_dist = es_dist, ts_dist = ts_dist,
                       icovariate = icovariate, null_dist = null_dist, execute = FALSE)

rank_scatter(onerun, pvalue = "pval", covariate = "ind_covariate")

strat_hist(onerun, pvalue = "pval", covariate = "ind_covariate", maxy = 5, numQ = 3)

Benchmark Metrics

We plot the average results across r B simulations.

resdf <- plotsim_standardize(res_i, alpha = seq(0.01, 0.10, 0.01))

plotsim_average(resdf, met="rejections", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE) 

plotsim_average(resdf, met="FDR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE) 

We also take a look at the distribution of rejects for each method as a function of the effect size and independent covariate.

covariateLinePlot(res_i, alpha = 0.05, covname = "effect_size", trans = "log1p")

covariateLinePlot(res_i, alpha = 0.05, covname = "ind_covariate", trans = "log1p")

Finally, (if enough methods produce rejections at r ualpha) we take a look at the overlap of rejections between methods.

if (numberMethodsReject(resdf, alphacutoff = ualpha, filterSet = excludeSet) >= 3) {
    aggupset(res_i, alpha = ualpha, supplementary = FALSE, return_list = FALSE) 
} else {
    message("Not enough methods found rejections at alpha ", ualpha, 
            "; skipping upset plot")
}

We also compare the simulation results with and without an informative covariate.

resdfu <- plotsim_standardize(res_u, alpha = seq(0.01, 0.10, 0.01))

resdfiu <- dplyr::full_join(select(resdf, rep, blabel, param.alpha, key,
                                   performanceMetric, alpha, value),
                            select(resdfu, rep, blabel, param.alpha, key,
                                   performanceMetric, alpha, value),
                            by = c("rep", "blabel", "param.alpha", "key",
                                   "performanceMetric", "alpha"),
                            suffix = c(".info", ".uninfo"))
resdfiu <- dplyr::mutate(resdfiu, value = value.info - value.uninfo)

plotsim_average(resdfiu, met="rejections", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE)

plotsim_average(resdfiu, met="FDR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE)

plotsim_average(resdfiu, met="TPR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE)

plotsim_average(resdfiu, met="TNR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE)

Chi-Squared Setting

Finally, we consider the null setting with chi-squared distributed test statistics.

Data Simulation

ts_dist <- rchisq_perturber(4)  # functional: sampling dist/noise for test stats
null_dist <- rchisq_pvaluer(4)     # functional: dist to calc p-values
seed <- 1023

For the chi-squared distributed setting, we must change the "mode" setting for ASH from the default of 0 to "estimate" since the mode of null and alternative test statistics are no longer centered at 0. While both approximate Normality and unimodality of effects are violated in this simulation setting, by allowing the mode to be estimated, rather than forced to 0, should return more comparable results.

bdchi <- modifyBMethod(bd, "ashq", mode = "empirical")

We next run the simulations. We do not include FDR Regression because the test statistics are not approximately normally distributed.

if (file.exists(chisq_file)) {
    res <- readRDS(chisq_file)
} else {
    res <- mclapply(X = 1:B, FUN = simIteration, bench = bdchi, m = m,
                    pi0 = pi0, es_dist = es_dist, icovariate = icovariate,
                    ts_dist = ts_dist, null_dist = null_dist,
                    seed = seed, mc.cores = cores)
    saveRDS(res, file = chisq_file)
}
res_i <- lapply(res, `[[`, "informative")
res_u <- lapply(res, `[[`, "uninformative")

Covariate Diagnostics

Here, we show the relationship between the independent covariate and p-values for a single replication of the experiment.

onerun <- simIteration(bdchi, m = m, pi0 = pi0, es_dist = es_dist, ts_dist = ts_dist,
                       icovariate = icovariate, null_dist = null_dist, execute = FALSE)

rank_scatter(onerun, pvalue = "pval", covariate = "ind_covariate")

strat_hist(onerun, pvalue = "pval", covariate = "ind_covariate", maxy = 5, numQ = 3)

Benchmark Metrics

We plot the average results across r B simulations.

resdf <- plotsim_standardize(res_i, alpha = seq(0.01, 0.10, 0.01))

plotsim_average(resdf, met="rejections", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE) 

plotsim_average(resdf, met="FDR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE) 

Since ashq is not appropriate when the test statistics are not approximately normal or t, we also plot the results without this method.

plotsim_average(resdf, met="rejections", filter_set = c(excludeSet, "ashq"),
                merge_ihw = TRUE, errorBars = TRUE) 

plotsim_average(resdf, met="FDR", filter_set = c(excludeSet, "ashq"),
                merge_ihw = TRUE, errorBars = TRUE) 

We also take a look at the distribution of rejects for each method as a function of the effect size and independent covariate. Again, we plot these without ashq.

covariateLinePlot(lapply(res_i, function(x) { x[, -which(colnames(x) == "ashq")] }),
                  alpha = ualpha, covname = "effect_size")

covariateLinePlot(lapply(res_i, function(x) { x[, -which(colnames(x) == "ashq")] }),
                  alpha = ualpha, covname = "ind_covariate")

Finally, (if enough methods produce rejections at r ualpha) we take a look at the overlap of rejections between methods. Again, without ashq.

if (numberMethodsReject(resdf, alphacutoff = ualpha, filterSet = c(excludeSet, "ashq")) >= 3) {
    aggupset(lapply(res_i, function(x) { x[, -which(colnames(x) == "ashq")] }),
             alpha = ualpha, supplementary = FALSE, return_list = FALSE)
} else {
    message("Not enough methods found rejections at alpha ", ualpha, 
            "; skipping upset plot")
}

We also compare the simulation results with and without an informative covariate.

resdfu <- plotsim_standardize(res_u, alpha = seq(0.01, 0.10, 0.01))

resdfiu <- dplyr::full_join(select(resdf, rep, blabel, param.alpha, key,
                                   performanceMetric, alpha, value),
                            select(resdfu, rep, blabel, param.alpha, key,
                                   performanceMetric, alpha, value),
                            by = c("rep", "blabel", "param.alpha", "key",
                                   "performanceMetric", "alpha"),
                            suffix = c(".info", ".uninfo"))
resdfiu <- dplyr::mutate(resdfiu, value = value.info - value.uninfo)

plotsim_average(resdfiu, met="rejections", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE)

plotsim_average(resdfiu, met="FDR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE)

plotsim_average(resdfiu, met="TPR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE)

plotsim_average(resdfiu, met="TNR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE)

Session Info

sessionInfo()

stephaniehicks/benchmarkfdrData2019 documentation built on June 20, 2021, 10 a.m.