bmdcalc | R Documentation |
Computes x-fold and z-SD benchmark doses for each responsive item using the best fit dose-reponse model.
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, ...)
f |
An object of class |
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 |
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 |
BMDtype |
The type of BMD to plot, |
plottype |
The type plot, |
by |
If different from |
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. |
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.
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 |
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). |
Marie-Laure Delignette-Muller and Elise Billoir
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")}
See uniroot
for details about the function used for the numerical
search of the benchmark dose for cases where there is no analytical solution.
# (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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.