1 Summary

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

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

3 Gaussian Setting

First, we consider the null setting with gaussian noise.

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

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(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)
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggplot2':
## 
##     ggsave

3.3 Benchmark Metrics

We plot the average results across 100 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) 
## Warning: Removed 45 rows containing missing values (geom_path).
## Warning: Removed 48 rows containing missing values (geom_errorbar).

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 0.05) 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")
}
## Not enough methods found rejections at alpha 0.05; 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)
## Warning: Removed 45 rows containing missing values (geom_path).
## Warning: Removed 50 rows containing missing values (geom_errorbar).

plotsim_average(resdfiu, met="TPR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE)
## Warning: Removed 100 rows containing missing values (geom_path).
## Warning: Removed 100 rows containing missing values (geom_errorbar).

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

4 Student’s t Setting (df = 5)

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

4.1 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 noise 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 noise 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(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)

4.3 Benchmark Metrics

We plot the average results across 100 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) 
## Warning: Removed 12 rows containing missing values (geom_path).
## Warning: Removed 12 rows containing missing values (geom_errorbar).

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")
## Warning: Removed 25 rows containing missing values (geom_path).
## Warning: Removed 25 rows containing missing values (geom_errorbar).

covariateLinePlot(res_i, alpha = 0.05, covname = "ind_covariate", trans = "log1p")
## Warning: Removed 25 rows containing missing values (geom_path).

## Warning: Removed 25 rows containing missing values (geom_errorbar).

Finally, (if enough methods produce rejections at 0.05) 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)
## Warning: Removed 17 rows containing missing values (geom_path).
## Warning: Removed 18 rows containing missing values (geom_errorbar).

plotsim_average(resdfiu, met="TPR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE)
## Warning: Removed 100 rows containing missing values (geom_path).
## Warning: Removed 100 rows containing missing values (geom_errorbar).

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

5 Student’s t Setting (df = 11)

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

5.1 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 noise 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 noise 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(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)

5.3 Benchmark Metrics

We plot the average results across 100 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) 
## Warning: Removed 10 rows containing missing values (geom_path).
## Warning: Removed 29 rows containing missing values (geom_errorbar).

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 0.05) 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)
## Warning: Removed 15 rows containing missing values (geom_path).
## Warning: Removed 30 rows containing missing values (geom_errorbar).

plotsim_average(resdfiu, met="TPR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE)
## Warning: Removed 100 rows containing missing values (geom_path).
## Warning: Removed 100 rows containing missing values (geom_errorbar).

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

6 Chi-Squared Setting

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

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

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

6.3 Benchmark Metrics

We plot the average results across 100 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) 
## Warning: Removed 25 rows containing missing values (geom_path).
## Warning: Removed 27 rows containing missing values (geom_errorbar).

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) 
## Warning: Removed 15 rows containing missing values (geom_path).
## Warning: Removed 17 rows containing missing values (geom_errorbar).

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 0.05) 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")
}
## Not enough methods found rejections at alpha 0.05; 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)
## Warning: Removed 25 rows containing missing values (geom_path).
## Warning: Removed 27 rows containing missing values (geom_errorbar).

plotsim_average(resdfiu, met="TPR", filter_set = excludeSet,
                merge_ihw = TRUE, errorBars = TRUE, diffplot = TRUE)
## Warning: Removed 80 rows containing missing values (geom_path).
## Warning: Removed 80 rows containing missing values (geom_errorbar).

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

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.2             
##  [3] hexbin_1.27.2              adaptMT_1.0.0             
##  [5] FDRreg_0.2-1               mvtnorm_1.0-8             
##  [7] BayesLogit_0.6             fda_2.4.7                 
##  [9] Matrix_1.2-14              fdrtool_1.2.15            
## [11] swfdr_1.4.0                qvalue_2.10.0             
## [13] ashr_2.2-7                 IHW_1.6.0                 
## [15] magrittr_1.5               SummarizedBenchmark_0.99.2
## [17] mclust_5.4.1               BiocParallel_1.12.0       
## [19] stringr_1.3.1              rlang_0.2.1               
## [21] UpSetR_1.3.3               SummarizedExperiment_1.8.1
## [23] DelayedArray_0.4.1         matrixStats_0.53.1        
## [25] Biobase_2.38.0             GenomicRanges_1.30.3      
## [27] GenomeInfoDb_1.14.0        IRanges_2.12.0            
## [29] S4Vectors_0.16.0           BiocGenerics_0.24.0       
## [31] tidyr_0.8.1                ggplot2_3.0.0             
## [33] dplyr_0.7.6               
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-137           bitops_1.0-6           doParallel_1.0.11     
##  [4] rprojroot_1.3-2        tools_3.5.0            backports_1.1.2       
##  [7] R6_2.2.2               lazyeval_0.2.1         colorspace_1.3-2      
## [10] withr_2.1.2            tidyselect_0.2.4       gridExtra_2.3         
## [13] mnormt_1.5-5           compiler_3.5.0         ggdendro_0.1-20       
## [16] labeling_0.3           slam_0.1-43            mosaicCore_0.5.0      
## [19] scales_0.5.0           SQUAREM_2017.10-1      psych_1.8.4           
## [22] digest_0.6.15          ggformula_0.7.0        foreign_0.8-70        
## [25] rmarkdown_1.9          XVector_0.18.0         pscl_1.5.2            
## [28] pkgconfig_2.0.1        htmltools_0.3.6        lpsymphony_1.7.1      
## [31] bindr_0.1.1            RCurl_1.95-4.10        GenomeInfoDbData_1.0.0
## [34] mosaicData_0.16.0      Rcpp_0.12.18           munsell_0.4.3         
## [37] stringi_1.2.2          yaml_2.1.19            MASS_7.3-50           
## [40] zlibbioc_1.24.0        plyr_1.8.4             grid_3.5.0            
## [43] lattice_0.20-35        knitr_1.20             pillar_1.2.3          
## [46] reshape2_1.4.3         codetools_0.2-15       glue_1.3.0            
## [49] evaluate_0.10.1        foreach_1.4.4          gtable_0.2.0          
## [52] purrr_0.2.5            assertthat_0.2.0       broom_0.4.4           
## [55] truncnorm_1.0-8        tibble_1.4.2           iterators_1.0.9       
## [58] mosaic_1.2.0