1 Summary

In this set of simulations, we consider settings with both null and non-null tests with informative and non-informative covariates as described in simulations-informative-step.Rmd.

This set of simulations differs from simulations-informative-step.Rmd only in the implementation of the AdaPT method for multiple testing correction. In the primary set of simulation results (simulations-informative-step.Rmd), we observed that the AdaPT method can lose power when the informative covariate is only weakly informative. While the method still controls the FDR, we observe that the loss in power is substantial, with AdaPT achieving a lower mean TPR across 100 replications when using a weakly informative covariate rather than a completely uninformative (random) covariate.

In this set of simulations, we assess whether this loss of power with a weakly informative covariate can be mitigated by allowing the AdaPT model selection procedure to decide between the same collection of spline models used before, and an intercept-only (null) model, which makes no use of the covariate. This approach was suggested by the authors of the AdaPT method through personal communication.

Two additional implementations of the AdaPT model are included in this analysis: adapt-withnull, the approach described above, and adapt-onlynull, where the method is forced to ignore the informative covariate. The latter (onlynull) implementation provides a baseline for comparing the two other AdaPT parameter settings.

2 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, "nullAdaPT-informative-step-benchmark-gaussian.rds")
tdist_file <- file.path(resdir, "nullAdaPT-informative-step-benchmark-t5.rds")
tdist11_file <- file.path(resdir, "nullAdaPT-informative-step-benchmark-t11.rds")
chisq_file <- file.path(resdir, "nullAdaPT-informative-step-benchmark-chisq4.rds")

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

## define bechmarking design
bd <- initializeBenchDesign()

As described in simulations-null.Rmd, we include Scott’s FDR Regression in the analysis for simulations with Gaussian or t-distributed test statistics. Again, we include both nulltype = "empirical" and nulltype = "theoretical".

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))

We add a couple alternative calls to AdaPT to test whether the method performance improves in cases of low “informativeness” if given the option to not use the informative covariate. This option is coded by the null model, "pi0 ~ 1", where only an intercept term is included in the model. Specifically, we add two alternative calls to the benchmark - one call with the null model added as one of the candidate models to pi_formulas= and mu_formulas=, and a second call with the null model given as the only candidate model to pi_formulas= and mu_formulas=.

bdplus <- addBMethod(bdplus, "adapt-withnull",
                     adaptMT::adapt_glm,
                     function(x) { x$qvals },
                     pvals = pval,
                     x = data.frame(icov = ind_covariate),
                     pi_formulas = c("~ 1", paste0("splines::ns(icov, df = ", seq(2, 10, 2), ")")),
                     mu_formulas = c("~ 1", paste0("splines::ns(icov, df = ", seq(2, 10, 2), ")")),
                     alphas = 0)

bdplus <- addBMethod(bdplus, "adapt-onlynull",
                     adaptMT::adapt_glm,
                     function(x) { x$qvals },
                     pvals = pval,
                     x = data.frame(icov = ind_covariate),
                     pi_formulas = "~ 1",
                     mu_formulas = "~ 1",
                     alphas = 0)

candycols2 <- rbind(candycols,
                    data.frame(Method = c("adapt-withnull", "adapt-onlynull"),
                               col = c("darkgoldenrod3", "darkgoldenrod3"),
                               lty = c("dashed", "dotted"),
                               stringsAsFactors = FALSE))

All simulation settings will share the following parameters.

m <- 20000                        # integer: number of hypothesis tests
pi0 <- pi0_step(0.90)             # numeric: proportion of null hypotheses
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

3 Gaussian Setting

First, we consider the setting with Gaussian test statistics.

3.1 Data Simulation

es_dist <- rnorm_generator(3)       # functional: dist of alternative test stats
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")

3.2 Covariate Diagnostics

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

onerun <- simIteration(1, 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 = 10, numQ = 3)
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggplot2':
## 
##     ggsave

3.3 Benchmark Metrics

We plot the averaged results across 100 replications.

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,
                palette = candycols2)

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

plotsim_average(resdf, met="TPR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE,
                palette = candycols2)

plotsim_average(resdf, met="TNR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE,
                palette = candycols2)

Since we are primarily interested in comparing just the AdaPT methods, we plot these metrics excluding all other methods.

plotsim_average(resdf, met="rejections",
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE,
                palette = candycols2)

plotsim_average(resdf, met="FDR", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE,
                palette = candycols2)

plotsim_average(resdf, met="TPR", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE,
                palette = candycols2)

plotsim_average(resdf, met="TNR", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE,
                palette = candycols2)

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 = ualpha, covname = "effect_size",
                  palette = candycols2)

covariateLinePlot(res_i, alpha = ualpha, covname = "ind_covariate",
                  palette = candycols2)

We also look at the FDR as a function of the independent covariate.

covariateLinePlot(res_i, alpha = ualpha, covname = "ind_covariate", metric = "FDR",
                  palette = candycols2)
## Warning: Removed 4 rows containing missing values (geom_path).
## Warning: Removed 6 rows containing missing values (geom_errorbar).

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,
                palette = candycols2)

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

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

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

Again, we recreate the plots after subsetting on AdaPT methods.

plotsim_average(resdfiu, met="rejections", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE,
                palette = candycols2)

plotsim_average(resdfiu, met="FDR", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE,
                palette = candycols2)

plotsim_average(resdfiu, met="TPR", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE,
                palette = candycols2)

plotsim_average(resdfiu, met="TNR", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE,
                palette = candycols2)

4 Student’s t Setting (df = 5)

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

4.1 Data Simulation

es_dist <- rnorm_generator(6)    # functional: dist of alternative test stats
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")

4.2 Covariate Diagnostics

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

onerun <- simIteration(1, bdplust, 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 = 10, numQ = 3)

4.3 Benchmark Metrics

We plot the averaged results across 100 replications.

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,
                palette = candycols2)

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

plotsim_average(resdf, met="TPR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE,
                palette = candycols2)

plotsim_average(resdf, met="TNR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE,
                palette = candycols2)

Since we are primarily interested in comparing just the AdaPT methods, we plot these metrics excluding all other methods.

plotsim_average(resdf, met="rejections",
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE,
                palette = candycols2)

plotsim_average(resdf, met="FDR", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE,
                palette = candycols2)

plotsim_average(resdf, met="TPR", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE,
                palette = candycols2)

plotsim_average(resdf, met="TNR", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE,
                palette = candycols2)

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 = ualpha, covname = "effect_size",
                  palette = candycols2)

covariateLinePlot(res_i, alpha = ualpha, covname = "ind_covariate",
                  palette = candycols2)

We also look at the FDR as a function of the independent covariate.

covariateLinePlot(res_i, alpha = ualpha, covname = "ind_covariate", metric = "FDR",
                  palette = candycols2)
## Warning: Removed 31 rows containing missing values (geom_path).
## Warning: Removed 31 rows containing missing values (geom_errorbar).

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,
                palette = candycols2)

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

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

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

Again, we recreate the plots after subsetting on AdaPT methods.

plotsim_average(resdfiu, met="rejections", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE,
                palette = candycols2)

plotsim_average(resdfiu, met="FDR", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE,
                palette = candycols2)

plotsim_average(resdfiu, met="TPR", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE,
                palette = candycols2)

plotsim_average(resdfiu, met="TNR", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE,
                palette = candycols2)

5 Student’s t Setting (df = 11)

Next, we consider a second setting with t-distributed test statistics.

5.1 Data Simulation

es_dist <- rnorm_generator(3)     # functional: dist of alternative test stats
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")

5.2 Covariate Diagnostics

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

onerun <- simIteration(1, bdplust, 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 = 10, numQ = 3)

5.3 Benchmark Metrics

We plot the averaged results across 100 replications.

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,
                palette = candycols2)

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

plotsim_average(resdf, met="TPR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE,
                palette = candycols2)

plotsim_average(resdf, met="TNR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE,
                palette = candycols2)

Since we are primarily interested in comparing just the AdaPT methods, we plot these metrics excluding all other methods.

plotsim_average(resdf, met="rejections",
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE,
                palette = candycols2)

plotsim_average(resdf, met="FDR", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE,
                palette = candycols2)

plotsim_average(resdf, met="TPR", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE,
                palette = candycols2)

plotsim_average(resdf, met="TNR", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE,
                palette = candycols2)

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 = ualpha, covname = "effect_size",
                  palette = candycols2)

covariateLinePlot(res_i, alpha = ualpha, covname = "ind_covariate",
                  palette = candycols2)

We also look at the FDR as a function of the independent covariate.

covariateLinePlot(res_i, alpha = ualpha, covname = "ind_covariate", metric = "FDR",
                  palette = candycols2)
## Warning: Removed 37 rows containing missing values (geom_path).
## Warning: Removed 39 rows containing missing values (geom_errorbar).

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,
                palette = candycols2)

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

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

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

Again, we recreate the plots after subsetting on AdaPT methods.

plotsim_average(resdfiu, met="rejections", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE,
                palette = candycols2)

plotsim_average(resdfiu, met="FDR", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE,
                palette = candycols2)

plotsim_average(resdfiu, met="TPR", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE,
                palette = candycols2)

plotsim_average(resdfiu, met="TNR", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE,
                palette = candycols2)

6 Chi-Squared Setting

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

6.1 Data Simulation

es_dist <- rnorm_generator(15)       # functional: dist of alternative test stats
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(bdplus, "ashq", mode = "empirical")

We also drop the FDRreg methods from the chi-squared setting.

bdchi <- dropBMethod(bdchi, "fdrreg-t")
bdchi <- dropBMethod(bdchi, "fdrreg-e")

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")

6.2 Covariate Diagnostics

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

onerun <- simIteration(1, 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 = 10, numQ = 3)

6.3 Benchmark Metrics

We plot the averaged results across 100 replications.

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,
                palette = candycols2)

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

plotsim_average(resdf, met = "TPR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE,
                palette = candycols2)

plotsim_average(resdf, met = "TNR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE,
                palette = candycols2)

Since we are primarily interested in comparing just the AdaPT methods, we plot these metrics excluding all other methods.

plotsim_average(resdf, met="rejections",
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE,
                palette = candycols2)

plotsim_average(resdf, met="FDR", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE,
                palette = candycols2)

plotsim_average(resdf, met="TPR", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE,
                palette = candycols2)

plotsim_average(resdf, met="TNR", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE,
                palette = candycols2)

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",
                  palette = candycols2)

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

We also look at the FDR as a function of the independent covariate.

covariateLinePlot(res_i, alpha = ualpha, covname = "ind_covariate", metric = "FDR",
                  palette = candycols2)
## Warning: Removed 32 rows containing missing values (geom_path).
## Warning: Removed 37 rows containing missing values (geom_errorbar).

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,
                palette = candycols2)

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

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

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

Again, we recreate the plots after subsetting on AdaPT methods.

plotsim_average(resdfiu, met="rejections", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE,
                palette = candycols2)

plotsim_average(resdfiu, met="FDR", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE,
                palette = candycols2)

plotsim_average(resdfiu, met="TPR", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE,
                palette = candycols2)

plotsim_average(resdfiu, met="TNR", 
                filter_set = grep("adapt", unique(resdf$blabel),
                                  invert = TRUE, value = TRUE), 
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE,
                palette = candycols2)

7 Session Info

sessionInfo()
## R version 3.5.0 (2018-04-23)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: CentOS Linux 7 (Core)
## 
## Matrix products: default
## BLAS: /usr/lib64/libblas.so.3.4.2
## LAPACK: /usr/lib64/liblapack.so.3.4.2
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
##  [1] splines   parallel  stats4    stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] bindrcpp_0.2.2              cowplot_0.9.3              
##  [3] hexbin_1.27.2               adaptMT_1.0.0              
##  [5] FDRreg_0.1                  fda_2.4.8                  
##  [7] Matrix_1.2-14               fdrtool_1.2.15             
##  [9] swfdr_1.6.0                 qvalue_2.12.0              
## [11] ashr_2.2-7                  IHW_1.8.0                  
## [13] magrittr_1.5                SummarizedBenchmark_0.99.2 
## [15] mclust_5.4.1                stringr_1.3.1              
## [17] rlang_0.2.2                 UpSetR_1.3.3               
## [19] SummarizedExperiment_1.10.1 DelayedArray_0.6.6         
## [21] BiocParallel_1.14.2         matrixStats_0.54.0         
## [23] Biobase_2.40.0              GenomicRanges_1.32.7       
## [25] GenomeInfoDb_1.16.0         IRanges_2.14.12            
## [27] S4Vectors_0.18.3            BiocGenerics_0.26.0        
## [29] tidyr_0.8.1                 ggplot2_3.0.0              
## [31] dplyr_0.7.7                
## 
## loaded via a namespace (and not attached):
##  [1] ggdendro_0.1-20        foreach_1.4.4          assertthat_0.2.0      
##  [4] ggstance_0.3.1         GenomeInfoDbData_1.1.0 ggrepel_0.8.0         
##  [7] mosaic_1.4.0           yaml_2.2.0             slam_0.1-43           
## [10] pillar_1.3.0           backports_1.1.2        lattice_0.20-35       
## [13] glue_1.3.0             digest_0.6.18          XVector_0.20.0        
## [16] colorspace_1.3-2       htmltools_0.3.6        plyr_1.8.4            
## [19] pkgconfig_2.0.2        broom_0.5.0            zlibbioc_1.26.0       
## [22] purrr_0.2.5            scales_1.0.0           mosaicData_0.17.0     
## [25] tibble_1.4.2           withr_2.1.2            lazyeval_0.2.1        
## [28] crayon_1.3.4           evaluate_0.12          nlme_3.1-137          
## [31] doParallel_1.0.14      MASS_7.3-51            truncnorm_1.0-8       
## [34] tools_3.5.0            munsell_0.5.0          compiler_3.5.0        
## [37] grid_3.5.0             RCurl_1.95-4.11        iterators_1.0.10      
## [40] labeling_0.3           mosaicCore_0.6.0       bitops_1.0-6          
## [43] rmarkdown_1.10         gtable_0.2.0           codetools_0.2-15      
## [46] reshape2_1.4.3         R6_2.3.0               gridExtra_2.3         
## [49] knitr_1.20             ggformula_0.9.0        bindr_0.1.1           
## [52] rprojroot_1.3-2        lpsymphony_1.8.0       stringi_1.2.4         
## [55] pscl_1.5.2             SQUAREM_2017.10-1      Rcpp_0.12.19          
## [58] tidyselect_0.2.5