bmdcalc: Computation of benchmark doses for responsive items
In aursiber/DRomics: Dose Response for Omics

bmdcalc

R Documentation

Computation of benchmark doses for responsive items

Description

Computes x-fold and z-SD benchmark doses for each responsive item using the best fit dose-reponse model.

Usage

bmdcalc(f, z = 1, x = 10, minBMD, ratio2switchinlog = 100)

## S3 method for class 'bmdcalc'
print(x, ...)

## S3 method for class 'bmdcalc'
plot(x, BMDtype = c("zSD", "xfold"), 
            plottype = c("ecdf", "hist", "density"), 
            by = c("none", "trend", "model", "typology"),
            hist.bins = 30, BMD_log_transfo = TRUE, ...)

Arguments

`f`	An object of class `"drcfit"` returned by the function `drcfit`.
`z`	Value of z defining the BMD-zSD as the dose at which the response is reaching y0 +/- z * SD, with y0 the level at the control given by the dose-response fitted model and SD the residual standard deviation of the dose-response fitted model.
`x`	Value of x given as a percentage and defining the BMD-xfold as the dose at which the response is reaching y0 +/- (x/100) * y0, with y0 the level at the control given by the dose-response fitted model. For `print` and `plot` functions, an object of class `"bmdcalc"`.
`minBMD`	minimal value for calculated BMDs, so a value considered negligible compared to the tested doses. If not given by the user this argument is fixed at the minimal non null tested dose divided by 100.
`ratio2switchinlog`	ratio between maximal and minimal tested doses above which the numerical computation (when the use of `uniroot` is necessary) of the BMD is performed on a log scale of dose.
`BMDtype`	The type of BMD to plot, `"zSD"` (default choice) or `"xfold"`.
`plottype`	The type plot, `"ecdf"` for an empirical cumulative distribution plot (default choice), `"hist"` for a histogram or `"density"` for a density plot.
`by`	If different from `"none"` the plot is split by trend (if `"trend"`), by model (if `"model"`) or by typology (if `"typology"`).
`hist.bins`	The number of bins, only used for histogram(s).
`BMD_log_transfo`	If TRUE, default option, a log transformation of the BMD is used in the plot.
`...`	further arguments passed to graphical or print functions.

Details

The two types of benchmark doses (BMD) proposed by the EFSA (2017) were computed for each responsive item using the best fit dose-reponse model previously obtained using the drcfit function (see Larras et al. 2018 for details):

the BMD-zSD defined as the dose at which the response is reaching y0 +/- z * SD, with y0 the level at the control given by the dose-response model, SD the residual standard deviation of the dose response model fit and z given as an input (z fixed to 1 by default),
the BMD-xfold defined as the dose at which the response is reaching y0 +/- (x/100) * y0, with y0 the level at the control given by the dose-response fitted model and x the percentage given as an input (x fixed at 10 by default.)

When there is no analytical solution for the BMD, it is numerically searched along the fitted curve using the uniroot function.

In cases where the BMD cannot be reached due to the asymptote at high doses, NaN is returned. In cases where the BMD is not reached at the highest tested dose, NA is returned. Very low BMD values obtained by extrapolation between 0 and the smallest non null tested dose, that correspond to very sensitive items (that we do not want to exclude), are thresholded at minBMD, an argument by default fixed at the smallest non null tested dose divided by 100, but that can be fixed by the user as what he considers to be a negligible dose.

Value

bmdcalc returns an object of class "bmdcalc", a list with 4 components:

`res`	a data frame reporting the results of the fit and BMD computation on each selected item sorted in the ascending order of the adjusted p-values returned by function `itemselect`. The different columns correspond to the identifier of each item (`id`), the row number of this item in the initial data set (`irow`), the adjusted p-value of the selection step (`adjpvalue`), the name of the best fit model (`model`), the number of fitted parameters (`nbpar`), the values of the parameters `b`, `c`, `d`, `e` and `f`, (`NA` for non used parameters), the residual standard deviation (`SDres`), the typology of the curve (`typology`, (16 class typology described in the help of the `drcfit` function)), the rough trend of the curve (`trend`) defined with four classes (U, bell, increasing or decreasing shape), the theoretical y value at the control (`y0`), the theoretical y value at the maximal dose `yatdosemax`), the theoretical y range for x within the range of tested doses (`yrange`), the maximal absolute y change (up or down) from the control(`maxychange`) and for biphasic curves the x value at which their extremum is reached (`xextrem`) and the corresponding y value (`yextrem`), the BMD-zSD value (`BMD.zSD`) with the corresponding BMR-zSD value (reached or not, `BMR.zSD`) and the BMD-xfold value (`BMD.xfold`) with the corresponding BMR-xfold value (reached or not, `BMR.xfold`).
`z`	Value of z given in input to define the BMD-zSD.
`x`	Value of x given in input as a percentage to define the BMD-xfold.
`minBMD`	minimal value for calculated BMDs given in input or fixed at the minimal non null tested dose divided by 100.
`ratio2switchinlog`	ratio between maximal and minimal tested doses above which the numerical computations are performed in a log scale (as given in input).
`omicdata`	The corresponding object given in input (component of itemselect).

Author(s)

Marie-Laure Delignette-Muller and Elise Billoir

References

EFSA Scientific Committee, Hardy A, Benford D, Halldorsson T, Jeger MJ, Knutsen KH, ... & Schlatter JR (2017). Update: use of the benchmark dose approach in risk assessment. EFSA Journal, 15(1), e04658.

Larras F, Billoir E, Baillard V, Siberchicot A, Scholz S, Wubet T, Tarkka M, Schmitt-Jansen M and Delignette-Muller ML (2018). DRomics: a turnkey tool to support the use of the dose-response framework for omics data in ecological risk assessment. Environmental science & technology.\Sexpr[results=rd]{tools:::Rd_expr_doi("10.1021/acs.est.8b04752")}

Examples


# (1) a toy example (a very small subsample of a microarray data set) 
#
datafilename <- system.file("extdata", "transcripto_very_small_sample.txt", package="DRomics")

# to test the package on a small (for a quick calculation) but not very small data set
# use the following commented line
# datafilename <- system.file("extdata", "transcripto_sample.txt", package="DRomics")

(o <- microarraydata(datafilename, check = TRUE, norm.method = "cyclicloess"))
(s_quad <- itemselect(o, select.method = "quadratic", FDR = 0.01))
(f <- drcfit(s_quad, progressbar = TRUE))
(r <- bmdcalc(f))
plot(r) 


# same plot in raw scale of BMD (without log transformation of BMD values)
plot(r, BMD_log_transfo = FALSE) 

# changing the values of z and x for BMD calculation
(rb <- bmdcalc(f, z = 2, x = 50))
plot(rb)



# Plot of fits with BMD values 

# example with the BMD-1SD
plot(f, BMDoutput = r, BMDtype = "zSD")

# example with the BMD-2SD
plot(f, BMDoutput = rb, BMDtype = "zSD")

# example with the BMD-xfold with x = 10 percent
plot(f, BMDoutput = r, BMDtype = "xfold")


# (2) an example on a microarray data set (a subsample of a greater data set) 
#

datafilename <- system.file("extdata", "transcripto_sample.txt", package="DRomics")

# to test the package on a small (for a quick calculation) but not very small data set
# use the following commented line
# datafilename <- system.file("extdata", "transcripto_sample.txt", package="DRomics")

(o <- microarraydata(datafilename, check = TRUE, norm.method = "cyclicloess"))
(s_quad <- itemselect(o, select.method = "quadratic", FDR = 0.01))
(f <- drcfit(s_quad, progressbar = TRUE))
(r <- bmdcalc(f))
plot(r) 

# different plots of BMD-zSD

plot(r, plottype = "hist") 
plot(r, plottype = "density") 
plot(r, plottype = "density", by = "trend") 
plot(r, plottype = "ecdf", by = "trend") 
plot(r, plottype = "ecdf", by = "model") 
plot(r, plottype = "ecdf", by = "typology") 

# a plot of BMD-xfold (by default BMD-zSD is plotted)
plot(r, BMDtype = "xfold", plottype = "hist", by = "typology", hist.bins = 10)

aursiber/DRomics documentation built on Feb. 6, 2024, 4:28 p.m.