- Journal List
- HHS Author Manuscripts
- PMC9997717

As a library, NLM provides access to scientific literature. Inclusion in an NLM database does not imply endorsem*nt of, or agreement with, the contents by NLM or the National Institutes of Health.

Learn more: PMC Disclaimer | PMC Copyright Notice

Comput Toxicol. Author manuscript; available in PMC 2024 Feb 1.

*Published in final edited form as:*

Comput Toxicol. 2023 Feb; 25: 100259.

Published online 2022 Dec 27. doi:10.1016/j.comtox.2022.100259

PMCID: PMC9997717

NIHMSID: NIHMS1865839

PMID: 36909352

Matthew W. Wheeler,^{†,}^{1} Sooyeong Lim,^{2} John House,^{1} Keith Shockley,^{1} A. John Bailer,^{2} Jennifer Fostel,^{3} Longlong Yang,^{3} Dawan Talley,^{3} Ashwin Raghuraman,^{3} Jeffery S. Gift,^{4} J. Allen Davis,^{5} Scott S. Auerbach,^{6} and Alison A. Motsinger-Reif^{1}

Author information Copyright and License information PMC Disclaimer

## Associated Data

- Supplementary Materials

## Abstract

The need to analyze the complex relationships observed in high-throughput toxicogenomic and other omic platforms has resulted in an explosion of methodological advances in computational toxicology. However, advancements in the literature often outpace the development of software researchers can implement in their pipelines, and existing software is frequently based on pre-specified workflows built from well-vetted assumptions that may not be optimal for novel research questions. Accordingly, there is a need for a stable platform and open-source codebase attached to a programming language that allows users to program new algorithms. To fill this gap, the Biostatistics and Computational Biology Branch of the National Institute of Environmental Health Sciences, in cooperation with the National Toxicology Program (NTP) and US Environmental Protection Agency (EPA), developed ToxicR, an open-source R programming package. The ToxicR platform implements many of the standard analyses used by the NTP and EPA, including dose-response analyses for continuous and dichotomous data that employ Bayesian, maximum likelihood, and model averaging methods, as well as many standard tests the NTP uses in rodent toxicology and carcinogenicity studies, such as the poly-K and Jonckheere trend tests. ToxicR is built on the same codebase as current versions of the EPA’s Benchmark Dose software and NTP’s BMDExpress software but has increased flexibility because it directly accesses this software. To demonstrate ToxicR, we developed a custom workflow to illustrate its capabilities for analyzing toxicogenomic data. The unique features of ToxicR will allow researchers in other fields to add modules, increasing its functionality in the future.

**Keywords: **Benchmark Dose, High-Throughput Pathway Analysis, Model Averaged Dose-Response

## 1. Introduction

The development of high-throughput bioassays to understand cellular-level genetic expression and molecular-level responses to chemical exposures has exploded over the past 15 years. Bioassays enable novel toxicology experiments that illuminate a chemical’s effects on biological processes and enhance mechanistic understanding of toxicological processes. To analyze the large amounts of data produced by these assays, scientists employ various computational and data analytic tools to extract biological signals while removing spurious noise. For example, the National Toxicology Program (NTP) developed BMDExpress 2 and 3 (Phillips et al., 2019) based on BMDExpress () to analyze toxicogenomic data, and the US Environmental Protection Agency (EPA) uses its Benchmark Dose Software (United States Environmental Protection Agency, 2012) to examine dose-response relationships.

Some R programs for dose-response modeling using high-throughput data, such as the recently released tclpfit2 R package (), have limitations. Such programs consider only continuous responses and use a single definition of adverse events (e.g., tclpfit2 defaults to the 1-SD definition of the benchmark dose (Crump, 1995) and requires effort beyond a simple function argument to modify this default). Additionally, many of these tools execute a specific workflow and are thus built based on a one-size-fits-all model. Accordingly, when implementing the tools, it is often challenging for labs to reproduce methodologies proposed by other researchers or modify them to a different domain. Further, in many cases, the methods are written as research-level code that is not sufficiently accessible to a broader audience or has undergone limited testing beyond a given study’s purposes. This multi-faceted domain transfer issue limits the adoption of new methodologies by the broader scientific community, ultimately impeding scientific discovery.

Considering these limitations, we developed ToxicR, an R programming language package, as a community-contributed but centrally managed repository for computational toxicology methodologies with the aim of allowing access to state-of-the-art methods for dose-response and high-throughput analysis so that researchers can create novel analysis pipelines that are fit for purpose. While ToxicR consists of methodologies pioneered by researchers at the National Institute of Environmental Health Sciences, NTP, and EPA, we anticipate that additional modules will be added from the broader community of researchers.

Existing tools for high-throughput dose-response analyses have limitations. For individual dose-response analyses, risk assessors use the EPA’s Benchmark Dose Software (BMDS) (US EPA,United States Environmental Protection Agency, 2022) and the PROAST (Slob, 2018) software platforms. While well-vetted, these packages are designed for single or, at most, a handful of dose-response analyses simultaneously. Consequently, they do not natively support the throughput needs of modern toxicology studies. BMDExpress (Phillips et al., 2019; Yang et al., 2007) and FastBMD () are well-developed state-of-the-art platforms for toxicogenomic dose-response analyses in high-throughput pipelines. Although the tools are feature-rich and provide various analysis options, they have limited flexibility because they are not programming interfaces.

On the other side of the spectrum are R package libraries for dose-response modeling, which also have limitations. A notable example is the drc R package (), which provides dose-response modeling tools in R based on maximum likelihood estimation. This package provides a comprehensive hierarchical modeling environment but has a disadvantage in that it does not implement benchmark dose analyses, which are standard in regulatory settings when estimating the point of departure. The bmd R package was developed on top of the drc R package (). This package performs similarly to BMDS and PROAST but does not compute profile likelihood-based confidence intervals, which are recommended for use in regulatory toxicology (see the US EPA Technical Guidance Document (US United States Environmental Protection Agency, 2012)). Further, the bmd R package allows frequentist inference and does not permit Bayesian approaches, which are preferred by the World Health Organization ().

Additionally, programming access to methodologies such as model averaging (; Simmons et al., 2015; ; Wheeler et al., 2020; ) and Bayesian inference is even more limited. Tools such as Bayesian BMD (BBMD) () provide access to model averaging and Bayesian methods, but it is currently with a pay-for-access model for full functionality. Further, the service is web-based and does not allow access to source code or use through a customizable programming interface.

To address the deficits of existing tools, we designed ToxicR to provide a programming interface similar to the drc and bmd R packages. ToxicR uses the same model source libraries as the EPA’s BMDS and the NTP’s BMDExpress. Because ToxicR is a programming interface, it offers increased modeling flexibility and customizability compared to other programs. The source code is distributed under the MIT license and thus all the algorithms are open source, allowing transparency. In addition to dose-response analyses, ToxicR provides access to various statistical analyses routinely used by the NTP to analyze bioassay data in addition to multiple trend tests (e.g., the Williams trend (Williams, 1971) and the poly-K test for trend in carcinogenicity studies ()). While some modifications were made to ensure the methods are accessible to a general audience, the core of the code is unchanged. Consequently, ToxicR provides community access to the standard analyses used by the NTP, ensuring transparency, reproducibility, and confidence when developing analysis pipelines.

## 2. Materials and Methods

There are numerous tools available to ToxicR users. Table 1 compares the functionality of ToxicR to that of similar R packages (i.e., drc/bmd) designed as dose-response modeling tools. ToxicR is currently the most complete of these packages, and the gap is expected to widen as functionality is added to ToxicR.

### Table 1:

Comparison of the dose-response functionality between ToxicR, drc/bmd, and other dose-response modeling platforms. The table shows that ToxicR has comparable or more features than most other platforms.

ToxicR Compared with Other Software | |||||
---|---|---|---|---|---|

R Package Dose-Response Capability | Other Software | ||||

ToxicR | drc/bmd | PROAST | BBMDS | EPA BMDS | |

Data Type | |||||

Continuous data | Yes | Yes | Yes | Yes | Yes |

Dichotomous data | Yes | Yes | Yes | Yes | Yes |

Count data | No | Yes | No | No^{2} | No^{2} |

DR Estimation Method | |||||

Maximum likelihood | Yes | Yes | Yes | No | Yes |

Bayesian | Yes | No | No | Yes | Yes^{3} |

Identical to EPA algorithms | Yes | No | No | No | Yes |

BMDL Calculation | |||||

Profile likelihood | Yes | No | Yes | No | Yes |

Wald-based | Yes | Yes | No | No | No |

Bootstrap | No^{1} | Yes | Yes | No | No |

MCMC | Yes | No | No | Yes | No |

MA Estimation | |||||

Frequentist | No | Yes | Yes | No | Yes |

Bayesian | Yes | No | No | Yes | Yes^{3} |

Open in a separate window

^{1}Although not available directly, this method can be implemented with ToxicR using the bootstrap package.

^{2}Though count data cannot be modeled, it can be categorized and one can model it as categorical.

^{3}Currently available only for dichotomous data.

To install the package, users type the command

$$\u2018\text{install}.\text{packages}(\u201c\text{ToxicR}\u201d)\u2019$$

inside the R prompt. The developmental source code can be downloaded from https://github.com/NIEHS/ToxicR, and the manual and compiled package are available at https://cran.r-project.org/web/packages/ToxicR/index.html. Occasionally, R’s evolving CRAN checks flag new issues with previously approved code (e.g., new compilers may deprecate C++ functionality currently used in ToxicR). As CRAN seeks to be compatible with a wide variety of computer systems, this may result in ToxicR being removed from CRAN until the code is changed to adhere to new requirements. If this happens, we recommend going to the GitHub repository where the latest release and pre-compiled binaries for most systems (e.g., Windows and macOS) are available with installation instructions.

### 2.1. Dose-Response Modeling

Dose-response modeling assumes a statistical model that is dependent upon observed data. Given a data type, ToxicR makes assumptions on the distribution of the response and the mean and, when applicable, variance in that distribution. Dichotomous data are assumed to have a binomial distribution, and continuous data are assumed to have either a normal or log-normal distribution, with the normal distribution having an additional specification of the variance proportional to a power of the mean. In the following, dose-response represents the expected response and is estimated using either maximum likelihood estimation or Bayesian approaches. Tables 2 and and33 outline the dose-response models currently available in ToxicR for dichotomous and continuous data, respectively.

### Table 2:

Dichotomous dose-response models available in ToxicR. Here, *φ*[·] is the cumulative standard normal distribution, and the quantal linear model is the multistage model of degree 1.

Model | Form |
---|---|

Weibull | $f[dose\text{}\mid \theta =(a,b,g)]=a+(1-a)\left[1-\mathrm{exp}\left(-b\times dose{\text{}}^{g}\right)\right]$ |

Probit | $f[dose\text{}\mid \theta =(b,g)]=\text{\Phi}(b+g\times dose\text{})$ |

Log-Probit | $f[dose\text{}\mid \theta =(a,b,g)]=a+(1-a)\text{\Phi}(b+g\times dose\text{})$ |

Logistic | $f[dose\text{}\mid \theta =(a,b,g)]=\frac{1}{1+\mathrm{exp}(-a-g\times dose\text{})}$ |

Log-Logistic | $f[dose\text{}\mid \theta =(a,b,g)]=a+\frac{1-a}{1+\mathrm{exp}(-b-g\times \mathrm{log}[dose])}$ |

Gamma | $f[dose\mid \theta =(a,b,g)]=a+(1-a){\int}_{0}^{b\times dose\text{}}\frac{1}{\text{\Gamma}(g)}{t}^{g-1}\mathrm{exp}(-t)dt$ |

Hill | $f[dose\text{}\mid \theta =(a,b,c,g)]=a+\frac{c\times (1-a)}{1+\mathrm{exp}(-b-g\times \mathrm{log}[dose\text{}])}$ |

Multistage | $f\left[dose\mid \theta =\left(a,{\left\{{b}_{i}\right\}}_{i=1}^{k}\right)\right]=a+(1-a)[1-\mathrm{exp}\left(-\sum _{i=1}^{k}{b}_{i}\times dos{e}^{i}\right)$ |

Open in a separate window

### Table 3:

Continuous dose-response models available in ToxicR. These are the same models available in the EPA’s BMDS 3.3 software.

Model | Form |
---|---|

Hill | $f[dose\text{}\mid \theta =(a,b,c,g)]=a+\frac{b\times dos{e}^{g}}{{c}^{g}+dos{e}^{g}}$ |

Exponential – 3 | $f[dose\text{}\mid \theta =(a,b,g)]=a\times \mathrm{exp}\left(\pm \text{b}\times {\text{dose}}^{\text{g}}\right)$ |

Exponential – 5 | $f[dose\text{}\theta =(a,b,c,g)]=a\times (c-(c-1)\mathrm{exp}\left[-{(b\times \text{dose})}^{\text{g}}\right]$ |

Power | $f[dose\text{}\mid \theta =(a,b,g)]=a+b\times dos{e}^{g}$ |

Polynomial | $f\left[dose\mid \theta =\left(a,{\left\{{b}_{i}\right\}}_{i=1}^{k}\right)\right]=a+\sum _{i=1}^{k}{b}_{i}\times dos{e}^{i}$ |

Open in a separate window

Bayesian estimation requires prior information. In ToxicR, by default, the model parameters’ prior distributions are specified as in Wheeler et al. (2020) for dichotomous data or Wheeler et al. (2022) for continuous data by default. Additionally, users can modify the model’s parameter prior distributions based on their needs. Bayesian estimation is performed using either the Laplace maximum a posteriori (MAP) (Gelman et al., 2013) approach, which is the default, or Markov chain Monte Carlo (MCMC) simulation (; Gelman et al., 2013).

The ToxicR dose-response modeling suite accesses the same code repository utilized by the EPA for the BMDS 3.x software and BMDExpress 3(https://github.com/auerbachs/BMDExpress-3/wiki), but ToxicR provides additional functionality. In contrast to these programs, which specify default behavior and are difficult to customize, ToxicR allows for customizable dose-response analyses. As an example, Bayesian analyses enable user-modifiable priors and MCMC analyses (Brooks et al., 2011). ToxicR also provides the option of rapid BMD analyses using Wald-based intervals (Buse, 1982) vs. the traditional and computationally more demanding likelihood ratio approach. The use of Wald-based confidence intervals follows the work of Ewald et al. (2021). The supplement provides further information on the implementation of fitting algorithms.

#### 2.1.2. Model Averaging

ToxicR conducts model averaging using methodologies defined in Wheeler et al. (2020) and Wheeler et al. (2022). The underlying algorithms are improved versions of the code developed in these works. The model averaging algorithms take advantage of multicore computers using the openMP library (Chandra et al., 2001) to increase computational speed. When possible, multiple models can be run simultaneously using different CPU cores. For computers with eight or more cores, this enables a model averaging fit in roughly the same time required for an individual model. For macOS users who want to control multi-threading, ToxicR must be recompiled. By default, openMP is not supported on macOS, and we supply only a single-threaded precompiled binary on this platform.

### 2.2. NTP Testing Routines

The NTP, which sets the standard for toxicological testing, has developed source code for analyzing numerous toxicology experiments, including long-term bioassays, that represents a specific implementation of well-known statistical tests for ordered alternatives. For example, the NTP uses the Jonckheere trend test (Jonckheere, 1954), a non-parametric global test for ordered trends, to determine whether response changes monotonically with dose. While there are many implementations of this trend test, ToxicR employs essentially the same algorithm used by the NTP with the same default decision rules.

NTP’s trend testing strategy first globally employs the Jonckheere trend test (the function ‘ntp_ ntp_jonckeere()’ in ToxicR), with the assumption that either Williams’ trend test (Williams, 1971) or Shirley’s trend test (Shirley, 1977) will be used to test individual monotonic deviations from the background (functions ‘ntp_williams()’ and ‘ntp_shirley’). If, however, the significance level given by Jonkheere’s test is less than 0.01, NTP defaults to Dunnett’s multiple comparison test (Dunnett, 1955) (function ‘ntp_dunnett()’) when Williams’ trend test is specified or Dunn’s multiple comparison test (Dunn, 1961) to maximize power (function ‘ntp_dunn()’). The NTP’s decision logic is maintained in ToxicR.

Finally, the NTP uses the poly-K test (; ) to adjust for exposure-induced mortality that occurs prior to terminal sacrifice, which may interfere with understanding the tumorgenicity of a given chemical. The original code, which we moved into C++, is available in ToxicR as the function ‘ntp_polyk()’. In the future, we plan to add additional routines to the package.

### 2.3. Dose-Response Plotting

ToxicR provides a unified platform for plotting graphics in R. All model fits are returned as objects that can be graphically plotted using the ‘plot’ command. A ggplot2 (Wickham, 2011) object is produced, allowing user customizability beyond the default settings. For example, by default, plots are made on the original dose/data scale; however, users may rescale the doses to the log scale by applying the ‘scale_x_continuous()’ function. For example, for a fit object named ‘ma.fit’, the command

$$\u2018\text{plot}(\text{ma}.\text{fit})+\text{scale}\_\text{x}\_\text{continuous}(\text{trans}=\u201d\text{pseudo}\_\text{log}\u201d)\u2019$$

produces a plot of the fit at log scale. Additional plotting functions that are unique to ToxicR help users understand the distribution of the model average fit. The functions ‘ma_continuous_fit()’ and ‘cleveland_plot()’ sort the BMD results by posterior probability in the model average for MCMC (ma_continuous_fit()) and Laplace (cleveland_plot()) methodologies.

## 3. Example Workflow

To demonstrate the use of ToxicR, we developed a custom workflow using a pipeline to analyze toxicogenomic data observed by Ramaiahgari et al. (2019). In that study, differentiated HepaRG cells (a human liver hepatocyte-like cell) plated in a 96-well format were exposed to one of 24 chemicals for 96 hours. Chemicals were administered at 10 dose levels (half-log spacing), and three biological replicates per dose level were generated except in the untreated control, for which there were 12 biological replicates. We investigate three of the chemicals and analyze perfluorooctanoic acid (PFOA) data on the CYP3A1 gene receptor using ToxicR 22.8.1.0.2.

The supplement provides code and data files for the analysis and illustrates how we used R to select data and model dose-response with a Bayesian Exponential-5 model using the ‘single_continuous_fit()’ function. As ToxicR was developed to enable analyses beyond the standards recommended in other programs, we conducted a nonstandard analysis. First, we estimated the dose that would cause a 1.5 standard-deviation shift in the mean expression data as the benchmark dose (Crump, 1984). We also changed the default model prior distributions in ToxicR, which are set for gross pathology data not expected to change more than several standard deviations from the background. As gene expression can change by multiple orders of magnitude, we relaxed the prior on the slope term ‘b’ to be more diffuse. Instead of having a default standard deviation of 1, we set this parameter to have a standard deviation of 100. The transparency achieved by providing the code and data used for our analysis in the supplemental material demonstrates the full transparency achievable when using ToxicR with R Markdown (). From a perspective of reproducibility of research and transparency, this provides a significant benefit for toxicology research.

An individual fit takes 0.1 seconds on a MacBook Pro with an eight-core 2.3 ghz processor. For each fit, a fit object is returned and plotted using the default ‘plot’ function. If additional functionality is desired, the output can be modified using ‘ggplot’ commands or by writing a custom function. We chose the latter route and wrote the ‘rev_toxicRplot()’ function provided in the supplement. Figure 1 shows these two plots and the differences between the default and user-defined settings, highlighting the utility of the ToxicR interface. By default, users have access to high-quality graphic functionality and can augment this functionality to suit their needs. The supplement shows how the interface can be used to perform model averaging and analyze the cytochrome p450 response to PFOA.

We extended this to compare the cytochrome p450 response for two similar drugs—rosiglitazone and troglitazone. We examined p450 genes because their expression is commonly affected following chemical exposure through activation of the pregane X receptor (PXR). In ToxicR, this comparison is relatively easy. In R, we selected genes containing CYP in their name from our sample datasets, which essentially allowed us to select all phase I CyP450 probes to determine if they are transcriptionally altered in response to rosiglitazone and troglitazone. Next, we modeled only the genes with a significant Williams’ trend test, a function given in ToxicR, to remove genes exhibiting no response. We then selected the model average of each chemical/gene probe using the side-by-side chemical/gene response plot given for each gene probe.

Figure 2 plots an active probe and compares response to the two drugs. Table 4 gives the BMD estimates for the selected genes. The BMD estimates are qualitatively similar. Consistent with what would be expected for chemicals with similar potency on PXR, P450 exhibits roughly equivalent BMD values for both drugs. The similarity of responses between the two drugs is unsurprising. At high doses (above the therapeutic range), both drugs likely activate the nuclear receptors PXR/CAR, which would explain the activation of CYP2C8, CYP3A4, and CYP3A5. While this workflow did not utilize every feature of ToxicR, it is illustrative of the features and functionality of ToxicR.

Figure 2:

Comparison of troglitazone and rosiglitazone dose-response curves for the CYP2C8_15146 TempOseqO probe.

### Table 4:

Benchmark dose calculations for troglitazone and rosiglitazone for all active CYP 450 probes. Credible intervals are computed using the Laplacian approximation described in Wheeler et al. (2020).

Benchmark Dose | ||
---|---|---|

Probe | Troglitazone | Rosiglitazone |

CYP51A1_1742 | 61560.8 (39631.8, 94151.6) | 89327.8 (56464.0,347674.8) |

CYP2C8_15146 | 4993.0 (3061.0, 8798.1) | 4112.375(2340.6,7946.3) |

CYP3A4_20661 | 2181.1 (1263.4, 3632.3) | 2141.3 (1314.1,3435.6) |

CYP3A5_27090 | 1986.2 (1082.8, 3537.4) | 1742.9 (1030.8,2935.9) |

CYP3A4_28338 | 2186.9 (1195.8, 3884.8) | 3145.8 (2132.7, 4466.0) |

CYP3A5_28493 | 20564.3 (11668.0,38450.6) | 8306.4 (4018.1, 22295.2) |

Open in a separate window

The above functionality uses the R programming interface, and its use assumes a more advanced skillset. We are currently developing an R Shiny interface (Sievert, 2020) to provide basic functionality for users with less expertise in the R programming language. This interface would provide standard workflows allowing users to use ToxicR with typical options in an easier to use point-and-click interface.

## 4. Discussion

The open-source ToxicR R package can be used as an analysis tool for data arising in toxicology experiments. Users can deploy the package with other analysis tools in the R programming environment and develop personalized workflows. We released ToxicR under the MIT license with publicly available source code, which allows a high level of scrutiny. If an analysis is combined with R Markdown, researchers and risk assessors will be able to provide code and documentation of their analysis, which will enable an extra level of transparency. While this package was designed to develop novel omic dose-response analyses, it can easily be used to model data from more traditional toxicology experiments. The same underlying code is used in ToxicR as in the US EPA’s Benchmark Dose Software, so both researchers and regulators can use the platform.

As new research methods become available, we will add functionality to ToxicR. For example, we plan to add models described in Aerts, Wheeler, and Abrahantes (2020), which will harmonize ToxicR with the dose-response recommendations of the European Food Safety Authority (John-More et al., 2022). As ToxicR is under open development at NIEHS and NTP in partnership with other researchers, we will offer rolling releases that implement improvements and bug fixes for existing features. Finally, although the examples shown here using ToxicR 22.8.1.0.2, we expect to add additional continuous modeling features in future releases (e.g., using the suggested EFSA suite of continuous models). As these developments become available, we will place new releases on CRAN.

Figure 1:

Default (top) and user-defined (bottom) plotting for continuous model fits using ggplot2 in ToxicR.

### Highlights:

Open-source Toxicology dose-response routines available in R.

Same underlying code base as used in BMDExpress and BMDS.

Allows for Bayesian models and model averaging.

Updates planned to include state-of-the-science developments.

## Supplementary Material

#### 4

Click here to view.^{(220K, zip)}

#### 3

Click here to view.^{(229K, zip)}

#### 2

Click here to view.^{(309K, zip)}

#### 5

Click here to view.^{(850K, zip)}

#### 1

Click here to view.^{(136K, docx)}

## Footnotes

CRediT author statement

**Matthew W. Wheeler:** conceptualization, methodology, software, writing – original draft, and reviewing editing **Sooyeong Lim:** software, and visualization, **John House:** formal analysis, visualization, writing – review & editing, **Keith Shockley:** supervision and software **A. John Bailer:** supervision, conceptualization, **Jennifer Fostel:** project administration, software **Dawan Talley:** software, **Ashwin Raghuraman:** software **Cari Martini:** software, **Jeffery S. Gift:** conceptualization, writing -review & editing, **J. Allen Davis:** conceptualization, writing -review editing **Scott S. Auerbach:** conceptualization, writing -review & editing and **Alison. A. Motsinger-Reif:** conceptualization, supervision, writing -review & editing.

Declaration of interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

**Publisher's Disclaimer: **This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

## References

- Aerts M, Wheeler MW, & Abrahantes JC (2020). An extended and unified modeling framework for benchmark dose estimation for both continuous and binary data. Environmetrics, 31(7), e2630. doi: 10.1002/env.2630 [PMC free article] [PubMed] [CrossRef] [Google Scholar]
- Bailer AJ., & Portier CJ. (1988). Effects of treatment-induced mortality and tumor-induced mortality on tests for carcinogenicity in small samples. Biometrics, 44(2), 417–431. doi: 10.2307/2531856 [PubMed] [CrossRef] [Google Scholar]
- Bieler GS, & Williams RL (1993). Ratio estimates, the delta method, and quantal response tests for increased carcinogenicity. Biometrics, 49(3), 793–801. Retrieved from https://www.ncbi.nlm.nih.gov/pubmed/8241374 [PubMed] [Google Scholar]
- Brooks S, Gelman A, Jones G, & Meng X-L (2011). Handbook of Markov Chain Monte Carlo: CRC press. [Google Scholar]
- Buse A. (1982). The likelihood ratio, Wald, and Lagrange multiplier tests: An expository note. The American Statistician, 36(3a), 153–157. [Google Scholar]
- Chandra R, Dagum L, Kohr D, Menon R, Maydan D, & McDonald J. (2001). Parallel programming in OpenMP: Morgan kaufmann. [Google Scholar]
- Crump KS (1984). A new method for determining allowable daily intakes. Toxicological Sciences, 4(5), 854–871. doi: 10.1016/0272-0590(84)90107-6 [PubMed] [CrossRef] [Google Scholar]
- Crump KS (1995). Calculation of benchmark doses from continuous data. Risk Analysis, 15(1), 79–89. [Google Scholar]
- Dunn OJ (1961). Multiple comparisons among means. Journal of the American Statistical Association, 56(293), 52–64. [Google Scholar]
- Dunnett CW (1955). A multiple comparison procedure for comparing several treatments with a control. Journal of the American Statistical Association, 50(272), 1096–1121. [Google Scholar]
- Ewald J, Soufan O, Xia J, & Basu N. (2021). FastBMD: an online tool for rapid benchmark dose–response analysis of transcriptomics data. Bioinformatics, 37(7), 1035–1036. doi: 10.1093/bioinformatics/btaa700 [PMC free article] [PubMed] [CrossRef] [Google Scholar]
- Gelman A, Carlin JB, Stern HS, Dunson DB, Vehtari A, & Rubin DB (2013). Bayesian Data Analysis: CRC Press. [Google Scholar]
- Jensen SM, Kluxen FM, Streibig JC, Cedergreen N, & Ritz C. (2020). bmd: an R package for benchmark dose estimation. PeerJ, 8, e10557. doi: 10.7717/peerj.10557 [PMC free article] [PubMed] [CrossRef] [Google Scholar]
- John-More S, Bampidis V, Benford D, Bragard CL, Halldorsson TI, F H-JA, . . . Schlatter J. (2022). Draft Guidance on the use of the Benchmark Dose approach in risk assessment. The EFSA Journal(Draft). [PMC free article] [PubMed] [Google Scholar]
- Jonckheere AR (1954). A distribution-free k-sample test against ordered alternatives. Biometrika, 41(1/2), 133–145. [Google Scholar]
- Phillips JR, Svoboda DL, Tandon A, Patel S, Sedykh A, Mav D, . . . Thomas RS (2019). BMDExpress 2: enhanced transcriptomic dose-response analysis workflow. Bioinformatics, 35(10), 1780–1782. doi: 10.1093/bioinformatics/bty878 [PMC free article] [PubMed] [CrossRef] [Google Scholar]
- Ramaiahgari SC, Auerbach SS, Saddler TO, Rice JR, Dunlap PE, Sipes NS, . . . Merrick BA (2019). The power of resolution: contextualized understanding of biological responses to liver injury chemicals using high-throughput transcriptomics and benchmark concentration modeling. Toxicological Sciences, 169(2), 553–566. [PMC free article] [PubMed] [Google Scholar]
- Ritz C, Baty F, Streibig JC, & Gerhard D. (2015). Dose-response analysis using R. PloS one, 10(12), e0146021. [PMC free article] [PubMed] [Google Scholar]
- Shao K, & Gift JS (2014). Model uncertainty and Bayesian model averaged benchmark dose estimation for continuous data. Risk Analysis, 34(1), 101–120. doi: 10.1111/risa.12078 [PubMed] [CrossRef] [Google Scholar]
- Shao K, & Shapiro AJ (2018). A web-based system for Bayesian benchmark dose estimation. Environmental health perspectives, 126(1), 017002. doi: 10.1289/EHP1289 [PMC free article] [PubMed] [CrossRef] [Google Scholar]
- Sheffield T, Brown J, Davidson S, Friedman KP, & Judson R. (2022). tcplfit2: an R-language general purpose concentration–response modeling package. Bioinformatics, 38(4), 1157–1158. [PMC free article] [PubMed] [Google Scholar]
- Shirley E. (1977). A non-parametric equivalent of Williams’ test for contrasting increasing dose levels of a treatment. Biometrics, 33(2), 386–389. Retrieved from https://www.ncbi.nlm.nih.gov/pubmed/884197 [PubMed] [Google Scholar]
- Sievert C. (2020). Interactive web-based data visualization with R, plotly, and shiny: CRC Press. [Google Scholar]
- Simmons SJ, Chen C, Li X, Wang Y, Piegorsch WW, Fang Q, . . . Dunn GE (2015).Bayesian model averaging for benchmark dose estimation. Environmental and Ecological Statistics, 22(1), 5–16. [Google Scholar]
- Slob W. (2018). Joint project on benchmark dose modelling with RIVM. EFSA Supporting Publications, 15(12), 1497E. [Google Scholar]
- United States Environmental Protection Agency. (2012). Benchmark Dose Technical Guidance. (EPA/100/R-12/001). US EPA [Google Scholar]
- United States Environmental Protection Agency. (2022). Benchmark Dose Software (Build 3.3; Model Library Version 2021.09). (Available from https://www.epa.gov/bmds/benchmark-dose-software-bmds-version-33-download). [Google Scholar]
- Wheeler MW, & Bailer AJ (2007). Properties of model-averaged BMDLs: A study of model averaging in dichotomous response risk estimation. Risk Analysis: An International Journal, 27(3), 659–670. doi: 10.1111/j.1539-6924.2007.00920.x [PubMed] [CrossRef] [Google Scholar]
- Wheeler MW, Blessinger T, Shao K, Allen BC, Olszyk L, Davis JA, & Gift JS (2020). Quantitative Risk Assessment: Developing a Bayesian Approach to Dichotomous Dose–Response Uncertainty. Risk Analysis, 40(9), 1706–1722. doi: 10.1111/risa.13537 [PMC free article] [PubMed] [CrossRef] [Google Scholar]
- Wheeler MW, Cortiñas Abrahantes J, Aerts M, Gift JS, & Allen Davis J. (2022). Continuous model averaging for benchmark dose analysis: Averaging over distributional forms. Environmetrics, e2728. doi: 10.1002/env.2728 [PMC free article] [PubMed] [CrossRef] [Google Scholar]
- WHO, & FAO. (2020). Principles and methods for the risk assessment of chemicals in food. Environmental Health Criteria, 240. [Google Scholar]
- Williams D. (1971). A test for differences between treatment means when several dose levels are compared with a zero dose control. Biometrics, 103–117. [PubMed] [Google Scholar]
- Xie Y, Dervieux C, & Riederer E. (2020). R markdown cookbook: Chapman and Hall/CRC. [Google Scholar]
- Yang L, Allen BC, & Thomas RS (2007). BMDExpress: a software tool for the benchmark dose analyses of genomic data. BMC genomics, 8(1), 1–8. doi: 10.1186/1471-2164-8-387 [PMC free article] [PubMed] [CrossRef] [Google Scholar]