diff --git a/.gitignore b/.gitignore
index 065f56108d9a9934746a7a4d890cdb3ff7ab55ef..eddaf3d24df4c6520b1ef99aa062c58441c54acb 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,3 +4,4 @@
.Ruserdata
Notizen.md
+*.png
diff --git a/.gitlab/merge_request_templates/merge_into_master.md b/.gitlab/merge_request_templates/merge_into_master.md
index 021f1edaab0642d0ad93098425c6951f99f151d8..670d525aab37918df1c66d7cd3c46122f9aa5e09 100644
--- a/.gitlab/merge_request_templates/merge_into_master.md
+++ b/.gitlab/merge_request_templates/merge_into_master.md
@@ -20,7 +20,7 @@ Comment on any notes that come up during devtools::check().
<!-- Please check `[x]` the applicable boxes. -->
* [ ] This merge includes bug fixes or changes in the functionality of the package (whether they are user-visible or not).
* [ ] update documentation including examples (if necessary)
- * [ ] update NEWS.md
+ * [ ] update and polish NEWS.md
* [ ] update version in DESCRIPTION
* [ ] update date in DESCRIPTION
* [ ] This merge introduces new features.
@@ -35,5 +35,5 @@ If this merge introduces a new version (major, minor or patch) to the
master branch do these steps after the branch has been successfully merged:
- git tag
- build binary and source package
-- gitlab release including the builds
+- gitlab release and upload the builds
-->
diff --git a/DESCRIPTION b/DESCRIPTION
index ff3def4338af703ffda55cd527818ea844937aa5..86af1897ed88671bb33a486fb22337d35b85a230 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
Package: stressaddition
Type: Package
Title: Modeling Tri-Phasic Concentration-Response Relationships
-Version: 2.0.0
-Date: 2020-02-17
+Version: 2.1.0
+Date: 2020-02-25
Authors@R: person("Sebastian", "Henz", role = c("aut", "cre"),
email = "sebastian.henz@ufz.de",
comment = c(ORCID = "0000-0001-8299-8852"))
@@ -19,10 +19,7 @@ Depends:
R (>= 3.5)
Imports:
drc (>= 3.0),
- plotrix,
- stats,
- graphics,
- grDevices
+ plotrix
Suggests:
testthat (>= 2.1.0)
RoxygenNote: 7.0.2
diff --git a/NAMESPACE b/NAMESPACE
index 06741e7eab9a60266f9993496c37c8b444cd9238..58db063863ead9fe8cb2bc7ceebf64bc51392656 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -6,6 +6,7 @@ export(effect_to_stress)
export(plot_effect)
export(plot_stress)
export(predict_ecxsys)
+export(predict_mixture)
export(stress_to_effect)
import(grDevices)
import(graphics)
diff --git a/NEWS.md b/NEWS.md
index 76e2ec2ef128af6718962136f839c8ac760bcf9c..012fd68e8f66b0d08fb57dc73bc41c32a2487e0a 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,12 @@
+# stressaddition 2.1.0
+
+* The functions `plot_effect()` and `plot_stress()` gain a `which` argument that controls which curves are plotted. Consequently, the `show_LL5_model` argument of `plot_effect()` was removed.
+* Added arguments `xlab` and `ylab` to `plot_stress`.
+* Added argument `main` to both plot functions.
+* Changed some colors of the stress curves so they better match with the colors of related effect curves.
+* Added `predict_mixture()` for the prediction of the effects of mixtures of two toxicants.
+* Fixed documentation of `ecxsys()` and `predict_ecxsys()`.
+
# stressaddition 2.0.0
* Changed the order of arguments in `ecxsys()`.
diff --git a/R/ec.R b/R/ec.R
index de29753db77fea5338f70e4ccad1b632553a3008..66d25927c69b34c14f5e498baed86cd72e971461 100644
--- a/R/ec.R
+++ b/R/ec.R
@@ -25,9 +25,9 @@
#' # below 90 % of the effect in the control.
#'
#' model <- ecxsys(
-#' concentration = c(0, 0.03, 0.3, 3, 10),
-#' effect_tox_observed = c(85, 76, 94, 35, 0),
-#' hormesis_concentration = 0.3
+#' concentration = c(0, 0.05, 0.5, 5, 30),
+#' hormesis_concentration = 0.5,
+#' effect_tox_observed = c(90, 81, 92, 28, 0)
#' )
#'
#' # using the ecxsys() output or the curves therein directly:
diff --git a/R/ecxsys.R b/R/ecxsys.R
index 591912098baf039efdcc9aab4b868e7ff707c784..95e090d4362f864d1a4b15a6a50c08bf55e87f1f 100644
--- a/R/ecxsys.R
+++ b/R/ecxsys.R
@@ -8,8 +8,8 @@
#' ECx-SyS
#'
-#' The ECx-SyS model for modeling concentration-effect relationships which
-#' indicate signs of hormesis.
+#' The ECx-SyS model for modeling concentration-effect relationships whith
+#' hormesis.
#'
#' It is advised to complete the curve down to zero for optimal prediction.
#' Therefore \code{effect_tox_observed} in the highest concentration should be
@@ -20,9 +20,8 @@
#' \code{effect_tox_env_observed} (if provided) must be of equal length and
#' sorted by increasing concentration.
#'
-#' @param concentration A vector of concentrations, one of which must be 0 to
-#' indicate the control. Should be sorted in ascending order, otherwise it
-#' will be sorted automatically.
+#' @param concentration A vector of concentrations. Must be sorted in ascending
+#' order and the first element must be 0 to indicate the control.
#' @param hormesis_concentration The concentration where the hormesis occurs.
#' This is usually the concentration of the highest effect after the control.
#' @param effect_tox_observed A vector of effect values observed at the given
@@ -34,28 +33,53 @@
#' survival data in percent this should be 100 (the default).
#' @param p,q The shape parameters of the beta distribution. Default is 3.2.
#'
-#' @return A list (of class ecxsys) containing many different objects with the
-#' most important being \code{curves}, a data frame containing effect and
-#' stress values at different concentrations. See \code{\link{predict_ecxsys}}
-#' for details.
+#' @return A list (of class ecxsys) containing many different objects of which
+#' the most important are listed below. The effect and stress vectors
+#' correspond to the provided concentrations.
+#' \describe{
+#' \item{effect_tox}{Modeled effect resulting from toxicant stress.}
+#' \item{effect_tox_sys}{Modeled effect resulting from toxicant and system
+#' stress.}
+#' \item{effect_tox_env}{Modeled effect resulting from toxicant and
+#' environmental stress.}
+#' \item{effect_tox_env_sys}{Modeled effect resulting from toxicant,
+#' environmental and system stress.}
+#' \item{effect_tox_LL5}{The effect predicted by the five-parameter
+#' log-logistic model derived from the observations under toxicant stress
+#' but without environmental stress.}
+#' \item{effect_tox_env_LL5}{The effect predicted by the five-parameter
+#' log-logistic model derived from the observations under toxicant stress
+#' with environmental stress.}
+#' \item{curves}{A data frame containing effect and stress values as
+#' returned by \code{\link{predict_ecxsys}}. The concentrations are
+#' regularly spaced on a logarithmic scale in the given concentration range.
+#' The control is approximated by the lowest non-control concentration times
+#' 1e-7. The additional column \code{use_for_plotting} is used by the
+#' plotting functions of this package to approximate the control and
+#' generate a break in the concentration axis.}
+#' }
#'
#' @examples model <- ecxsys(
-#' concentration = c(0, 0.03, 0.3, 3, 10),
-#' hormesis_concentration = 0.3,
-#' effect_tox_observed = c(85, 76, 94, 35, 0),
-#' effect_tox_env_observed = c(24, 23, 32, 0, 0)
+#' concentration = c(0, 0.05, 0.5, 5, 30),
+#' hormesis_concentration = 0.5,
+#' effect_tox_observed = c(90, 81, 92, 28, 0),
+#' effect_tox_env_observed = c(29, 27, 33, 5, 0)
#' )
#'
-#' # Use effect_max if for example the effect is given as the number of
-#' # surviving animals and the initial number of animals is 20:
+#' # Use effect_max if for example the effect is given as the average number of
+#' # surviving animals and the initial number of animals is 21:
#' model <- ecxsys(
-#' concentration = c(0, 0.03, 0.3, 3, 10),
+#' concentration = c(0, 0.03, 0.3, 3, 30),
#' hormesis_concentration = 0.3,
#' effect_tox_observed = c(17, 15.2, 18.8, 7, 0),
#' effect_tox_env_observed = c(4.8, 4.6, 6.4, 0, 0),
-#' effect_max = 20
+#' effect_max = 21
#' )
#'
+#' @references \href{https://doi.org/10.1038/s41598-019-51645-4}{Liess, M.,
+#' Henz, S. & Knillmann, S. Predicting low-concentration effects of
+#' pesticides. Sci Rep 9, 15248 (2019).}
+#'
#' @export
ecxsys <- function(concentration,
hormesis_concentration,
@@ -111,8 +135,7 @@ ecxsys <- function(concentration,
if (any(is.na(c(all_observations, concentration)))) {
stop("Values containing NA are not supported.")
}
- if (any(all_observations > effect_max) ||
- any(all_observations < 0)) {
+ if (any(all_observations > effect_max) || any(all_observations < 0)) {
stop("Observed effect must be between 0 and effect_max.")
}
conc_shift <- 2 # Powers of ten to shift the control downwards from the
diff --git a/R/helpers.R b/R/helpers.R
index f13f0d09f1c7da52f091e498dda5960e0bb3ea95..1b9fea35ab1ac236d2b0d9d296a7d483263247dd 100644
--- a/R/helpers.R
+++ b/R/helpers.R
@@ -33,12 +33,12 @@ get_log_ticks <- function(x) {
}
-adjust_smooth_concentrations <- function(curves, conc_adjust_factor) {
+adjust_smooth_concentrations <- function(model) {
# Helper for the plot functions, not exported.
# Deals with the concentrations which are unnecessary for plotting.
# This means it removes the concentrations in the gap and scales the
# concentrations left of the gap up.
-
+ curves <- model$curves
gap_idx <- min(which(!curves$use_for_plotting))
# Keep only the values to the left and right of the axis break:
@@ -50,7 +50,7 @@ adjust_smooth_concentrations <- function(curves, conc_adjust_factor) {
# Shift the small concentrations upwards so the plot has a nice x-axis:
curves$concentration[1:gap_idx] <-
- curves$concentration[1:gap_idx] / conc_adjust_factor
+ curves$concentration[1:gap_idx] / model$conc_adjust_factor
list(
curves = curves,
diff --git a/R/plot_ecxsys.R b/R/plot_ecxsys.R
index 0ea5e493068b8d9d98eb359c181b357f546ef6ea..18a31e67edf2e11215e40f180023411f7498de84 100644
--- a/R/plot_ecxsys.R
+++ b/R/plot_ecxsys.R
@@ -3,24 +3,34 @@
#' Plot the results of the ECx-SyS model
#'
-#' Convenience functions to plot the observed and modeled effect and the
-#' system stress with and without environmental stress.
+#' Plot the observed and modeled effects and stresses.
#'
#' @name plot_ecxsys
#'
#' @param model The list returned from \code{\link{ecxsys}}.
-#' @param show_LL5_model Should the log-logistic models be plotted for
-#' comparison? Defaults to \code{FALSE}.
-#' @param show_legend Should the plot include a legend? Defaults to FALSE
-#' because it may cover some points or lines depending on the plot size.
-#' @param xlab,ylab Axis labels.
+#' @param which A vector of curve names to plot. Allowed values are the column
+#' names of the \code{model$curves} data frame. The default \code{NULL} only
+#' plots the most important curves. Use \code{which = "all"} to display all
+#' curves.
+#' @param show_legend Should the plot include a legend? Defaults to \code{FALSE}
+#' because it may cover some parts of the plot depending on the plot size and
+#' the number of elements shown.
+#' @param xlab,ylab,main Axis labels and title.
#'
#' @examples model <- ecxsys(
-#' concentration = c(0, 0.03, 0.3, 3, 10),
-#' hormesis_concentration = 0.3,
-#' effect_tox_observed = c(85, 76, 94, 35, 0),
-#' effect_tox_env_observed = c(24, 23, 32, 0, 0)
+#' concentration = c(0, 0.05, 0.5, 5, 30),
+#' hormesis_concentration = 0.5,
+#' effect_tox_observed = c(90, 81, 92, 28, 0),
+#' effect_tox_env_observed = c(29, 27, 33, 5, 0)
#' )
#' plot_effect(model)
#' plot_stress(model)
+#'
+#' # Plot all curves:
+#' plot_effect(model, which = "all")
+#' plot_stress(model, which = "all")
+#'
+#' # Plot only some selected curves:
+#' plot_effect(model, which = c("effect_tox", "effect_tox_env"))
+#' plot_stress(model, which = c("stress_tox", "stress_tox_env"))
NULL
diff --git a/R/plot_effect.R b/R/plot_effect.R
index 6b944e3a46ab9ea9c3eccb9ba87317c5c3fac931..7fa995dd952e3042e9e0a95f07dde5bdaf64f19c 100644
--- a/R/plot_effect.R
+++ b/R/plot_effect.R
@@ -1,19 +1,49 @@
#' @rdname plot_ecxsys
#' @export
plot_effect <- function(model,
- show_LL5_model = FALSE,
+ which = NULL,
show_legend = FALSE,
xlab = "concentration",
- ylab = "effect") {
+ ylab = "effect",
+ main = NULL) {
stopifnot(inherits(model, "ecxsys"))
- temp <- adjust_smooth_concentrations(
- model$curves,
- model$conc_adjust_factor
- )
+
+ curve_names <- names(model$curves)
+ valid_names <- curve_names[startsWith(curve_names, "effect")]
+ if (is.null(which)) {
+ which <- c("effect_tox_sys")
+ if (model$with_env) {
+ which <- c(which, "effect_tox_env_sys")
+ }
+ } else if ("all" %in% which) {
+ if (length(which) == 1) {
+ which <- valid_names
+ } else {
+ stop("'all' must not be combined with other curve names.")
+ }
+ } else if (!all(which %in% valid_names)) {
+ warning("Argument 'which' contains invalid names.")
+ if (!model$with_env && any(grepl("env", which, fixed = TRUE))) {
+ warning("'which' contains names with 'env' but the model was",
+ " built without environmental effects.")
+ }
+ which <- which[which %in% valid_names]
+ }
+
+ temp <- adjust_smooth_concentrations(model)
curves <- temp$curves
log_ticks <- get_log_ticks(curves$concentration)
concentration <- c(curves$concentration[1], model$args$concentration[-1])
+ legend_df <- data.frame(
+ text = character(),
+ pch = numeric(),
+ lty = numeric(),
+ col = character(),
+ order = numeric(), # controls sorting of legend elements
+ stringsAsFactors = FALSE
+ )
+
plot(
NA,
NA,
@@ -22,93 +52,103 @@ plot_effect <- function(model,
log = "x",
xlab = xlab,
ylab = ylab,
+ main = main,
xaxt = "n",
- bty = "L",
- las = 1
+ yaxt = "n",
+ bty = "L"
)
- lines(
- curves$concentration,
- curves$effect_tox_sys,
- col = "blue"
- )
- lines(
- curves$concentration,
- curves$effect_tox,
- col = "deepskyblue",
- lty = 2
- )
points(
concentration,
model$args$effect_tox_observed,
pch = 16,
col = "blue"
)
-
+ legend_df[nrow(legend_df) + 1, ] <- list("tox (observed)", 16, 0, "blue", 1)
if (model$with_env) {
+ points(
+ concentration,
+ model$args$effect_tox_env_observed,
+ pch = 16,
+ col = "red"
+ )
+ legend_df[nrow(legend_df) + 1, ] <- list("tox + env (observed)", 16, 0, "red", 2)
+ }
+ # The lines are drawn in this order to ensure that dotted and dashed lines
+ # are on top of solid lines for better visibility.
+ if ("effect_tox_sys" %in% which) {
lines(
curves$concentration,
- curves$effect_tox_env_sys,
- col = "red"
+ curves$effect_tox_sys,
+ col = "blue"
)
+ legend_df[nrow(legend_df) + 1, ] <- list("tox + sys", NA, 1, "blue", 5)
+ }
+ if ("effect_tox" %in% which) {
lines(
curves$concentration,
- curves$effect_tox_env,
- col = "orange",
+ curves$effect_tox,
+ col = "deepskyblue",
lty = 2
)
- points(
- concentration,
- model$args$effect_tox_env_observed,
- pch = 16,
- col = "red"
- )
+ legend_df[nrow(legend_df) + 1, ] <- list("tox", NA, 2, "deepskyblue", 4)
}
-
- if (show_LL5_model) {
+ if ("effect_tox_LL5" %in% which) {
lines(
curves$concentration,
curves$effect_tox_LL5,
col = "darkblue",
- lty = 4
+ lty = 3
)
- if (model$with_env) {
+ legend_df[nrow(legend_df) + 1, ] <- list("tox (LL5)", NA, 3, "darkblue", 3)
+ }
+ if (model$with_env) {
+ if ("effect_tox_env_sys" %in% which) {
+ lines(
+ curves$concentration,
+ curves$effect_tox_env_sys,
+ col = "red"
+ )
+ legend_df[nrow(legend_df) + 1, ] <- list("tox + env + sys", NA, 1, "red", 8)
+ }
+ if ("effect_tox_env" %in% which) {
+ lines(
+ curves$concentration,
+ curves$effect_tox_env,
+ col = "orange",
+ lty = 2
+ )
+ legend_df[nrow(legend_df) + 1, ] <- list("tox + env", NA, 2, "orange", 7)
+ }
+ if ("effect_tox_env_LL5" %in% which) {
lines(
curves$concentration,
curves$effect_tox_env_LL5,
col = "darkred",
- lty = 4
+ lty = 3
)
+ legend_df[nrow(legend_df) + 1, ] <- list("tox + env (LL5)", NA, 3, "darkred", 6)
}
}
- axis(1, at = log_ticks$major, labels = log_ticks$major_labels)
- axis(1, at = log_ticks$minor, labels = FALSE, tcl = -0.25)
+ # The setting of col = NA and col.ticks = par("fg") is to prevent ugly line
+ # thickness issues when plotting as a png with type = "cairo" and at a low
+ # resolution.
+ axis(1, at = log_ticks$major, labels = log_ticks$major_labels,
+ col = NA, col.ticks = par("fg"))
+ axis(1, at = log_ticks$minor, labels = FALSE, tcl = -0.25,
+ col = NA, col.ticks = par("fg"))
plotrix::axis.break(1, breakpos = temp$axis_break_conc)
+ axis(2, col = NA, col.ticks = par("fg"), las = 1)
if (show_legend) {
- legend_text <- c("tox", "tox + sys")
- legend_pch <- c(NA, 16)
- legend_lty <- c(2, 1)
- legend_col <- c("deepskyblue", "blue")
- if (model$with_env) {
- legend_text <- c(legend_text, "tox + env", "tox + env + sys")
- legend_pch <- c(legend_pch, NA, 16)
- legend_lty <- c(legend_lty, 2, 1)
- legend_col <- c(legend_col, "orange", "red")
- }
- if (show_LL5_model) {
- legend_text <- c(legend_text, "tox (LL5)", "tox + env (LL5)")
- legend_pch <- c(legend_pch, NA, NA)
- legend_lty <- c(legend_lty, 4, 4)
- legend_col <- c(legend_col, "darkblue", "darkred")
- }
+ legend_df <- legend_df[order(legend_df$order), ]
legend(
"topright",
- legend = legend_text,
- pch = legend_pch,
- lty = legend_lty,
- col = legend_col
+ legend = legend_df$text,
+ pch = legend_df$pch,
+ lty = legend_df$lty,
+ col = legend_df$col
)
}
}
diff --git a/R/plot_stress.R b/R/plot_stress.R
index cccdcbcb0e7176262087fa331d74806514e90c63..7fc70f8f5e4970c0620c7ea15b87ab9a4be61fd3 100644
--- a/R/plot_stress.R
+++ b/R/plot_stress.R
@@ -1,72 +1,162 @@
#' @rdname plot_ecxsys
#' @export
-plot_stress <- function(model, show_legend = FALSE) {
+plot_stress <- function(model,
+ which = NULL,
+ show_legend = FALSE,
+ xlab = "concentration",
+ ylab = "stress",
+ main = NULL) {
stopifnot(inherits(model, "ecxsys"))
- temp <- adjust_smooth_concentrations(
- model$curves,
- model$conc_adjust_factor
- )
+
+ curve_names <- names(model$curves)
+ valid_names <- curve_names[
+ startsWith(curve_names, "stress") | startsWith(curve_names, "sys")
+ ]
+ if (is.null(which)) {
+ which <- c("sys_tox")
+ if (model$with_env) {
+ which <- c(which, "sys_tox_env")
+ }
+ } else if ("all" %in% which) {
+ if (length(which) == 1) {
+ which <- valid_names
+ } else {
+ stop("'all' must not be combined with other curve names.")
+ }
+ } else if (!all(which %in% valid_names)) {
+ warning("Argument 'which' contains invalid names.")
+ if (!model$with_env && any(grepl("env", which, fixed = TRUE))) {
+ warning("'which' contains names with 'env' but the model was",
+ " built without environmental effects.")
+ }
+ which <- which[which %in% valid_names]
+ if (length(which) == 0) {
+ stop("No curves to display.")
+ }
+ }
+
+ temp <- adjust_smooth_concentrations(model)
curves <- temp$curves
- axis_break_conc <- temp$axis_break_conc
log_ticks <- get_log_ticks(curves$concentration)
concentration <- c(curves$concentration[1], model$args$concentration[-1])
+ legend_df <- data.frame(
+ text = character(),
+ pch = numeric(),
+ lty = numeric(),
+ col = character(),
+ order = numeric(), # controls sorting of legend elements
+ stringsAsFactors = FALSE
+ )
+
plot(
NA,
NA,
xlim = range(curves$concentration, na.rm = TRUE),
- ylim = c(0, 1),
+ ylim = c(0, max(curves[, which], 1, na.rm = TRUE)),
log = "x",
- xlab = "concentration",
- ylab = "system stress",
+ xlab = xlab,
+ ylab = ylab,
+ main = main,
xaxt = "n",
- las = 1,
+ yaxt = "n",
bty = "L"
)
- lines(
- curves$concentration,
- curves$sys_tox,
- col = "blue"
- )
- points(
- concentration,
- model$sys_tox_not_fitted,
- pch = 16,
- col = "blue"
- )
-
- if (model$with_env) {
+ # The lines are drawn in this order to ensure that dotted and dashed lines
+ # are on top of solid lines for better visibility.
+ if ("stress_tox_sys" %in% which) {
+ lines(
+ curves$concentration,
+ curves$stress_tox_sys,
+ col = "blue"
+ )
+ legend_df[nrow(legend_df) + 1, ] <- list("tox + sys", NA, 1, "blue", 3)
+ }
+ if ("stress_tox" %in% which) {
lines(
curves$concentration,
- curves$sys_tox_env,
- col = "red"
+ curves$stress_tox,
+ col = "deepskyblue",
+ lty = 2
)
+ legend_df[nrow(legend_df) + 1, ] <- list("tox", NA, 2, "deepskyblue", 1)
+ }
+ if ("sys_tox" %in% which) {
points(
concentration,
- model$sys_tox_env_not_fitted,
+ model$sys_tox_not_fitted,
pch = 16,
- col = "red"
+ col = "steelblue3"
)
+ lines(
+ curves$concentration,
+ curves$sys_tox,
+ col = "steelblue3"
+ )
+ legend_df[nrow(legend_df) + 1, ] <- list("sys (tox)", 16, 1, "steelblue3", 2)
+ }
+ if (model$with_env) {
+ if ("stress_tox_env_sys" %in% which) {
+ lines(
+ curves$concentration,
+ curves$stress_tox_env_sys,
+ col = "red"
+ )
+ legend_df[nrow(legend_df) + 1, ] <- list("tox + env + sys", NA, 1, "red", 7)
+ }
+ if ("stress_env" %in% which) {
+ lines(
+ curves$concentration,
+ curves$stress_env,
+ col = "forestgreen",
+ lty = 3
+ )
+ legend_df[nrow(legend_df) + 1, ] <- list("env", NA, 3, "forestgreen", 4)
+ }
+ if ("stress_tox_env" %in% which) {
+ lines(
+ curves$concentration,
+ curves$stress_tox_env,
+ col = "orange",
+ lty = 2
+ )
+ legend_df[nrow(legend_df) + 1, ] <- list("tox + env", NA, 2, "orange", 5)
+ }
+ if ("sys_tox_env" %in% which) {
+ points(
+ concentration,
+ model$sys_tox_env_not_fitted,
+ pch = 16,
+ col = "violetred"
+ )
+ lines(
+ curves$concentration,
+ curves$sys_tox_env,
+ col = "violetred"
+ )
+ legend_df[nrow(legend_df) + 1, ] <- list("sys (tox + env)", 16, 1, "violetred", 6)
+ }
}
- axis(1, at = log_ticks$major, labels = log_ticks$major_labels)
- axis(1, at = log_ticks$minor, labels = FALSE, tcl = -0.25)
- plotrix::axis.break(1, breakpos = axis_break_conc)
+ # The setting of col = NA and col.ticks = par("fg") is to prevent ugly line
+ # thickness issues when plotting as a png with type = "cairo" and at a low
+ # resolution.
+ axis(1, at = log_ticks$major, labels = log_ticks$major_labels,
+ col = NA, col.ticks = par("fg"))
+ axis(1, at = log_ticks$minor, labels = FALSE, tcl = -0.25,
+ col = NA, col.ticks = par("fg"))
+ plotrix::axis.break(1, breakpos = temp$axis_break_conc)
+ axis(2, col = NA, col.ticks = par("fg"), las = 1)
if (show_legend) {
- legend_text <- c("tox")
- legend_col <- c("blue")
- if (model$with_env) {
- legend_text <- c(legend_text, "tox + env")
- legend_col <- c(legend_col, "red")
- }
+ legend_df <- legend_df[order(legend_df$order), ]
legend(
"topright",
- legend = legend_text,
- col = legend_col,
- pch = 16,
- lty = 1
+ legend = legend_df$text,
+ pch = legend_df$pch,
+ lty = legend_df$lty,
+ col = legend_df$col
)
}
}
diff --git a/R/predict_ecxsys.R b/R/predict_ecxsys.R
index 604dce43729b25ca640cd239027df06e92c6167b..7d3d82a5a48eac83aff288cc68771a56ba007272 100644
--- a/R/predict_ecxsys.R
+++ b/R/predict_ecxsys.R
@@ -1,19 +1,19 @@
-#' Predict ECxSyS at various concentrations
+#' Predict effects and stresses
#'
-#' @param model The output of a call to \code{\link{ecxsys}}.
+#' Calculate the effects and stresses of an ECx-SyS model at arbitrary
+#' concentrations.
+#'
+#' @param model An ECx-SyS model as returned by \code{\link{ecxsys}}.
#' @param concentration A numeric vector of concentrations.
#'
#' @return A data frame (of class "ecxsys_predicted") with the following
#' columns:
#' \describe{
-#' \item{concentration}{Concentrations regularly spaced on a logarithmic
-#' scale in the given concentration range. The control is approximated by
-#' the lowest non-control concentration times 1e-7.}
-#' \item{effect_tox_LL5}{The five-parameter log-logistic model of the
-#' effect derived from the observations under toxicant stress but without
-#' environmental stress.}
-#' \item{effect_tox}{Modeled effect resulting from toxicant and system
-#' stress.}
+#' \item{concentration}{The supplied concentrations.}
+#' \item{effect_tox_LL5}{The effect predicted by the five-parameter
+#' log-logistic model derived from the observations under toxicant stress
+#' but without environmental stress.}
+#' \item{effect_tox}{Modeled effect resulting from toxicant stress.}
#' \item{effect_tox_sys}{Modeled effect resulting from toxicant and system
#' stress.}
#' \item{stress_tox}{The toxicant stress.}
@@ -21,9 +21,9 @@
#' without environmental stress.}
#' \item{stress_tox_sys}{The sum of \code{stress_tox} and
#' \code{sys_tox}.}
-#' \item{effect_tox_env_LL5}{The five-parameter log-logistic model of the
-#' effect derived from the observations under toxicant stress with
-#' environmental stress.}
+#' \item{effect_tox_env_LL5}{The effect predicted by the five-parameter
+#' log-logistic model derived from the observations under toxicant stress
+#' with environmental stress.}
#' \item{effect_tox_env}{Modeled effect resulting from toxicant and
#' environmental stress.}
#' \item{effect_tox_env_sys}{Modeled effect resulting from toxicant,
@@ -37,9 +37,10 @@
#' }
#'
#' @examples model <- ecxsys(
-#' concentration = c(0, 0.03, 0.3, 3, 10),
-#' hormesis_concentration = 0.3,
-#' effect_tox_observed = c(85, 76, 94, 35, 0)
+#' concentration = c(0, 0.05, 0.5, 5, 30),
+#' hormesis_concentration = 0.5,
+#' effect_tox_observed = c(90, 81, 92, 28, 0),
+#' effect_tox_env_observed = c(29, 27, 33, 5, 0)
#' )
#' p <- predict_ecxsys(model, c(0.001, 0.01, 0.1, 1, 10))
#'
diff --git a/R/predict_mixture.R b/R/predict_mixture.R
new file mode 100644
index 0000000000000000000000000000000000000000..253503e9f7cbfa3f036508471a34493fd557e5c5
--- /dev/null
+++ b/R/predict_mixture.R
@@ -0,0 +1,98 @@
+#' Predict the effect of a mixture of two toxicants
+#'
+#' This method is only suitable for experiments without environmental stress.
+#' Any environmental stress in \code{model_1} or \code{model_2} is ignored.
+#'
+#' @param model_1 The ecxsys model of the toxicant that varies in concentration.
+#' @param model_2 The ecxsys model of the toxicant with a fixed concentration.
+#' @param concentration_1 The concentration of toxicant 1 in the mixture.
+#' @param concentration_2 The concentration of toxicant 2 in the mixture.
+#' @param proportion_ca How much of the combined toxicant stress is determined
+#' by concentration addition. Must be between 0 and 1.
+#'
+#' @return A vector of the effects of the mixture, scaled to the effect_max
+#' value of model_1.
+#'
+#' @examples toxicant_model_1 <- ecxsys(
+#' concentration = c(0, 0.05, 0.5, 5, 30),
+#' hormesis_concentration = 0.5,
+#' effect_tox_observed = c(90, 81, 92, 28, 0),
+#' )
+#' toxicant_model_2 <- ecxsys(
+#' concentration = c(0, 0.1, 1, 10, 100, 1000),
+#' hormesis_concentration = 10,
+#' effect_tox_observed = c(26, 25, 24, 27, 5, 0),
+#' effect_max = 30
+#' )
+#' predict_mixture(
+#' toxicant_model_1,
+#' toxicant_model_2,
+#' c(0, 0.02, 0.2, 2, 20),
+#' 3
+#' )
+#' @export
+predict_mixture <- function(model_1,
+ model_2,
+ concentration_1,
+ concentration_2,
+ proportion_ca = 0.5) {
+ stopifnot(
+ inherits(model_1, "ecxsys"),
+ inherits(model_2, "ecxsys"),
+ is.numeric(concentration_1),
+ is.numeric(concentration_2),
+ length(concentration_2) == 1,
+ !is.na(concentration_2),
+ proportion_ca >= 0,
+ proportion_ca <= 1
+ )
+
+ predicted_model_1 <- predict_ecxsys(model_1, concentration_1)
+ predicted_model_2 <- predict_ecxsys(model_2, concentration_2)
+
+ # tox stress ----------------------------------------------------------
+ stress_tox_sam <- predicted_model_1$stress_tox + predicted_model_2$stress_tox
+
+ # Convert the model_2 concentration into an equivalent model_1 concentration
+ # anc vice versa. Clamp the response levels because drc::ED can't deal with
+ # 0 and 100.
+ response_level_2 <- 100 - predicted_model_2$effect_tox / model_2$args$effect_max * 100
+ response_level_2 <- clamp(response_level_2, 1e-10, 100 - 1e-10)
+ concentration_2_equivalent <- drc::ED(
+ model_1$effect_tox_mod,
+ response_level_2,
+ display = FALSE
+ )[, "Estimate"]
+ effect_tox_ca_1 <- predict(
+ model_1$effect_tox_mod,
+ data.frame(concentration = concentration_1 + concentration_2_equivalent)
+ )
+
+ response_level_1 <- 100 - predicted_model_1$effect_tox / model_1$args$effect_max * 100
+ response_level_1 <- clamp(response_level_1, 1e-10, 100 - 1e-10)
+ concentration_1_equivalent <- drc::ED(
+ model_2$effect_tox_mod,
+ response_level_1,
+ display = FALSE
+ )[, "Estimate"]
+ effect_tox_ca_2 <- predict(
+ model_2$effect_tox_mod,
+ data.frame(concentration = concentration_2 + concentration_1_equivalent)
+ )
+
+ stress_tox_ca <- effect_to_stress(effect_tox_ca_1 * 0.5 + effect_tox_ca_2 * 0.5)
+
+ # sys -----------------------------------------------------------------
+ sys_1 <- predict(model_1$sys_tox_mod, data.frame(stress_tox = stress_tox_ca))
+ sys_2 <- predict(model_2$sys_tox_mod, data.frame(stress_tox = stress_tox_ca))
+ sys_total <- sys_1 * 0.5 + sys_2 * 0.5
+
+ # combined stress and result ------------------------------------------
+ proportion_sam <- 1 - proportion_ca
+ stress_tox_total <- stress_tox_ca * proportion_ca + stress_tox_sam * proportion_sam
+ stress_total <- stress_tox_total + sys_total
+ effect_total <- stress_to_effect(stress_total) * model_1$args$effect_max
+
+ # unname() to remove the name when concentration_1 is a single number.
+ unname(effect_total)
+}
diff --git a/man/ec.Rd b/man/ec.Rd
index 1a9c8289114d58707f4386f1edbbaddeb3da18a0..dea65818c30d40ae7981cfa4cc8ba61084de4267 100644
--- a/man/ec.Rd
+++ b/man/ec.Rd
@@ -38,9 +38,9 @@ smallest concentration is returned.
# below 90 \% of the effect in the control.
model <- ecxsys(
- concentration = c(0, 0.03, 0.3, 3, 10),
- effect_tox_observed = c(85, 76, 94, 35, 0),
- hormesis_concentration = 0.3
+ concentration = c(0, 0.05, 0.5, 5, 30),
+ hormesis_concentration = 0.5,
+ effect_tox_observed = c(90, 81, 92, 28, 0)
)
# using the ecxsys() output or the curves therein directly:
diff --git a/man/ecxsys.Rd b/man/ecxsys.Rd
index ca2cabc4ec802efb1ca64fa7e8d25a6530c0d351..620adb3536863dfb892d7411497aa117e38ae0fa 100644
--- a/man/ecxsys.Rd
+++ b/man/ecxsys.Rd
@@ -15,9 +15,8 @@ ecxsys(
)
}
\arguments{
-\item{concentration}{A vector of concentrations, one of which must be 0 to
-indicate the control. Should be sorted in ascending order, otherwise it
-will be sorted automatically.}
+\item{concentration}{A vector of concentrations. Must be sorted in ascending
+order and the first element must be 0 to indicate the control.}
\item{hormesis_concentration}{The concentration where the hormesis occurs.
This is usually the concentration of the highest effect after the control.}
@@ -35,14 +34,35 @@ survival data in percent this should be 100 (the default).}
\item{p, q}{The shape parameters of the beta distribution. Default is 3.2.}
}
\value{
-A list (of class ecxsys) containing many different objects with the
- most important being \code{curves}, a data frame containing effect and
- stress values at different concentrations. See \code{\link{predict_ecxsys}}
- for details.
+A list (of class ecxsys) containing many different objects of which
+ the most important are listed below. The effect and stress vectors
+ correspond to the provided concentrations.
+ \describe{
+ \item{effect_tox}{Modeled effect resulting from toxicant stress.}
+ \item{effect_tox_sys}{Modeled effect resulting from toxicant and system
+ stress.}
+ \item{effect_tox_env}{Modeled effect resulting from toxicant and
+ environmental stress.}
+ \item{effect_tox_env_sys}{Modeled effect resulting from toxicant,
+ environmental and system stress.}
+ \item{effect_tox_LL5}{The effect predicted by the five-parameter
+ log-logistic model derived from the observations under toxicant stress
+ but without environmental stress.}
+ \item{effect_tox_env_LL5}{The effect predicted by the five-parameter
+ log-logistic model derived from the observations under toxicant stress
+ with environmental stress.}
+ \item{curves}{A data frame containing effect and stress values as
+ returned by \code{\link{predict_ecxsys}}. The concentrations are
+ regularly spaced on a logarithmic scale in the given concentration range.
+ The control is approximated by the lowest non-control concentration times
+ 1e-7. The additional column \code{use_for_plotting} is used by the
+ plotting functions of this package to approximate the control and
+ generate a break in the concentration axis.}
+ }
}
\description{
-The ECx-SyS model for modeling concentration-effect relationships which
-indicate signs of hormesis.
+The ECx-SyS model for modeling concentration-effect relationships whith
+hormesis.
}
\details{
It is advised to complete the curve down to zero for optimal prediction.
@@ -56,20 +76,25 @@ sorted by increasing concentration.
}
\examples{
model <- ecxsys(
- concentration = c(0, 0.03, 0.3, 3, 10),
- hormesis_concentration = 0.3,
- effect_tox_observed = c(85, 76, 94, 35, 0),
- effect_tox_env_observed = c(24, 23, 32, 0, 0)
+ concentration = c(0, 0.05, 0.5, 5, 30),
+ hormesis_concentration = 0.5,
+ effect_tox_observed = c(90, 81, 92, 28, 0),
+ effect_tox_env_observed = c(29, 27, 33, 5, 0)
)
-# Use effect_max if for example the effect is given as the number of
-# surviving animals and the initial number of animals is 20:
+# Use effect_max if for example the effect is given as the average number of
+# surviving animals and the initial number of animals is 21:
model <- ecxsys(
- concentration = c(0, 0.03, 0.3, 3, 10),
+ concentration = c(0, 0.03, 0.3, 3, 30),
hormesis_concentration = 0.3,
effect_tox_observed = c(17, 15.2, 18.8, 7, 0),
effect_tox_env_observed = c(4.8, 4.6, 6.4, 0, 0),
- effect_max = 20
+ effect_max = 21
)
}
+\references{
+\href{https://doi.org/10.1038/s41598-019-51645-4}{Liess, M.,
+ Henz, S. & Knillmann, S. Predicting low-concentration effects of
+ pesticides. Sci Rep 9, 15248 (2019).}
+}
diff --git a/man/plot_ecxsys.Rd b/man/plot_ecxsys.Rd
index b5d862e04a97a426af65ef99066a053c7d4c27ac..19d5c61ed5b184d8a9d0dd157f877806c3c8b0b5 100644
--- a/man/plot_ecxsys.Rd
+++ b/man/plot_ecxsys.Rd
@@ -8,36 +8,54 @@
\usage{
plot_effect(
model,
- show_LL5_model = FALSE,
+ which = NULL,
show_legend = FALSE,
xlab = "concentration",
- ylab = "effect"
+ ylab = "effect",
+ main = NULL
)
-plot_stress(model, show_legend = FALSE)
+plot_stress(
+ model,
+ which = NULL,
+ show_legend = FALSE,
+ xlab = "concentration",
+ ylab = "stress",
+ main = NULL
+)
}
\arguments{
\item{model}{The list returned from \code{\link{ecxsys}}.}
-\item{show_LL5_model}{Should the log-logistic models be plotted for
-comparison? Defaults to \code{FALSE}.}
+\item{which}{A vector of curve names to plot. Allowed values are the column
+names of the \code{model$curves} data frame. The default \code{NULL} only
+plots the most important curves. Use \code{which = "all"} to display all
+curves.}
-\item{show_legend}{Should the plot include a legend? Defaults to FALSE
-because it may cover some points or lines depending on the plot size.}
+\item{show_legend}{Should the plot include a legend? Defaults to \code{FALSE}
+because it may cover some parts of the plot depending on the plot size and
+the number of elements shown.}
-\item{xlab, ylab}{Axis labels.}
+\item{xlab, ylab, main}{Axis labels and title.}
}
\description{
-Convenience functions to plot the observed and modeled effect and the
-system stress with and without environmental stress.
+Plot the observed and modeled effects and stresses.
}
\examples{
model <- ecxsys(
- concentration = c(0, 0.03, 0.3, 3, 10),
- hormesis_concentration = 0.3,
- effect_tox_observed = c(85, 76, 94, 35, 0),
- effect_tox_env_observed = c(24, 23, 32, 0, 0)
+ concentration = c(0, 0.05, 0.5, 5, 30),
+ hormesis_concentration = 0.5,
+ effect_tox_observed = c(90, 81, 92, 28, 0),
+ effect_tox_env_observed = c(29, 27, 33, 5, 0)
)
plot_effect(model)
plot_stress(model)
+
+# Plot all curves:
+plot_effect(model, which = "all")
+plot_stress(model, which = "all")
+
+# Plot only some selected curves:
+plot_effect(model, which = c("effect_tox", "effect_tox_env"))
+plot_stress(model, which = c("stress_tox", "stress_tox_env"))
}
diff --git a/man/predict_ecxsys.Rd b/man/predict_ecxsys.Rd
index 194496d1710b8c30498a83f92e7690adec1dcc6b..5e86de00b1dfa2f0a3fbc7f5bf0ad940391833aa 100644
--- a/man/predict_ecxsys.Rd
+++ b/man/predict_ecxsys.Rd
@@ -2,12 +2,12 @@
% Please edit documentation in R/predict_ecxsys.R
\name{predict_ecxsys}
\alias{predict_ecxsys}
-\title{Predict ECxSyS at various concentrations}
+\title{Predict effects and stresses}
\usage{
predict_ecxsys(model, concentration)
}
\arguments{
-\item{model}{The output of a call to \code{\link{ecxsys}}.}
+\item{model}{An ECx-SyS model as returned by \code{\link{ecxsys}}.}
\item{concentration}{A numeric vector of concentrations.}
}
@@ -15,14 +15,11 @@ predict_ecxsys(model, concentration)
A data frame (of class "ecxsys_predicted") with the following
columns:
\describe{
- \item{concentration}{Concentrations regularly spaced on a logarithmic
- scale in the given concentration range. The control is approximated by
- the lowest non-control concentration times 1e-7.}
- \item{effect_tox_LL5}{The five-parameter log-logistic model of the
- effect derived from the observations under toxicant stress but without
- environmental stress.}
- \item{effect_tox}{Modeled effect resulting from toxicant and system
- stress.}
+ \item{concentration}{The supplied concentrations.}
+ \item{effect_tox_LL5}{The effect predicted by the five-parameter
+ log-logistic model derived from the observations under toxicant stress
+ but without environmental stress.}
+ \item{effect_tox}{Modeled effect resulting from toxicant stress.}
\item{effect_tox_sys}{Modeled effect resulting from toxicant and system
stress.}
\item{stress_tox}{The toxicant stress.}
@@ -30,9 +27,9 @@ A data frame (of class "ecxsys_predicted") with the following
without environmental stress.}
\item{stress_tox_sys}{The sum of \code{stress_tox} and
\code{sys_tox}.}
- \item{effect_tox_env_LL5}{The five-parameter log-logistic model of the
- effect derived from the observations under toxicant stress with
- environmental stress.}
+ \item{effect_tox_env_LL5}{The effect predicted by the five-parameter
+ log-logistic model derived from the observations under toxicant stress
+ with environmental stress.}
\item{effect_tox_env}{Modeled effect resulting from toxicant and
environmental stress.}
\item{effect_tox_env_sys}{Modeled effect resulting from toxicant,
@@ -46,13 +43,15 @@ A data frame (of class "ecxsys_predicted") with the following
}
}
\description{
-Predict ECxSyS at various concentrations
+Calculate the effects and stresses of an ECx-SyS model at arbitrary
+concentrations.
}
\examples{
model <- ecxsys(
- concentration = c(0, 0.03, 0.3, 3, 10),
- hormesis_concentration = 0.3,
- effect_tox_observed = c(85, 76, 94, 35, 0)
+ concentration = c(0, 0.05, 0.5, 5, 30),
+ hormesis_concentration = 0.5,
+ effect_tox_observed = c(90, 81, 92, 28, 0),
+ effect_tox_env_observed = c(29, 27, 33, 5, 0)
)
p <- predict_ecxsys(model, c(0.001, 0.01, 0.1, 1, 10))
diff --git a/man/predict_mixture.Rd b/man/predict_mixture.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..bc43850edd096f54a6841cf3d35623ca6920fc86
--- /dev/null
+++ b/man/predict_mixture.Rd
@@ -0,0 +1,53 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/predict_mixture.R
+\name{predict_mixture}
+\alias{predict_mixture}
+\title{Predict the effect of a mixture of two toxicants}
+\usage{
+predict_mixture(
+ model_1,
+ model_2,
+ concentration_1,
+ concentration_2,
+ proportion_ca = 0.5
+)
+}
+\arguments{
+\item{model_1}{The ecxsys model of the toxicant that varies in concentration.}
+
+\item{model_2}{The ecxsys model of the toxicant with a fixed concentration.}
+
+\item{concentration_1}{The concentration of toxicant 1 in the mixture.}
+
+\item{concentration_2}{The concentration of toxicant 2 in the mixture.}
+
+\item{proportion_ca}{How much of the combined toxicant stress is determined
+by concentration addition. Must be between 0 and 1.}
+}
+\value{
+A vector of the effects of the mixture, scaled to the effect_max
+ value of model_1.
+}
+\description{
+This method is only suitable for experiments without environmental stress.
+Any environmental stress in \code{model_1} or \code{model_2} is ignored.
+}
+\examples{
+toxicant_model_1 <- ecxsys(
+ concentration = c(0, 0.05, 0.5, 5, 30),
+ hormesis_concentration = 0.5,
+ effect_tox_observed = c(90, 81, 92, 28, 0),
+)
+toxicant_model_2 <- ecxsys(
+ concentration = c(0, 0.1, 1, 10, 100, 1000),
+ hormesis_concentration = 10,
+ effect_tox_observed = c(26, 25, 24, 27, 5, 0),
+ effect_max = 30
+)
+predict_mixture(
+ toxicant_model_1,
+ toxicant_model_2,
+ c(0, 0.02, 0.2, 2, 20),
+ 3
+)
+}
diff --git a/tests/testthat/test-ec.R b/tests/testthat/test-ec.R
index f5deb23002d3f9a41c442050332ee4c942af94cf..a4cb7402f2422f5d3d4001b4857def212c0f7a27 100644
--- a/tests/testthat/test-ec.R
+++ b/tests/testthat/test-ec.R
@@ -1,13 +1,9 @@
-context("effect concentration")
-
-
test_that("all input formats produce identical models", {
model <- ecxsys(
- concentration = c(0, 0.03, 0.3, 3, 10),
- hormesis_concentration = 0.3,
- effect_tox_observed = c(85, 76, 94, 35, 0)
+ concentration = c(0, 0.05, 0.5, 5, 30),
+ hormesis_concentration = 0.5,
+ effect_tox_observed = c(90, 81, 92, 28, 0)
)
-
# using the ecxsys() output or the curves therein directly:
ec10_a <- ec(model, "effect_tox_sys", 10)
ec10_b <- ec(model$curves, "effect_tox_sys", 10)
@@ -24,57 +20,57 @@ test_that("all input formats produce identical models", {
)
ec10_d <- ec(df_custom, "foo", 10)
- expect_identical(ec10_a, ec10_b)
- expect_identical(ec10_b, ec10_c)
- expect_identical(ec10_c, ec10_d)
+ expect_equal(ec10_a, ec10_b, tolerance = 1e-5)
+ expect_equal(ec10_b, ec10_c, tolerance = 1e-5)
+ expect_equal(ec10_c, ec10_d, tolerance = 1e-5)
})
test_that("ec values have not changed", {
model <- ecxsys(
- concentration = c(0, 0.03, 0.3, 3, 10),
- hormesis_concentration = 0.3,
- effect_tox_observed = c(85, 76, 94, 35, 0),
- effect_tox_env_observed = c(24, 23, 32, 0, 0)
+ concentration = c(0, 0.05, 0.5, 5, 30),
+ hormesis_concentration = 0.5,
+ effect_tox_observed = c(90, 81, 92, 28, 0),
+ effect_tox_env_observed = c(29, 27, 33, 5, 0)
)
expect_equal(
ec(model, "effect_tox_sys", 50),
- list(response_value = 42.389844, concentration = 2.547712),
+ list(response_value = 44.95368, concentration = 3.375735),
tolerance = 1e-4
)
expect_equal(
ec(model, "effect_tox_sys", 10),
- list(response_value = 76.301719, concentration = 1.015731),
+ list(response_value = 80.91662, concentration = 1.098648),
tolerance = 1e-4
)
expect_equal(
ec(model, "effect_tox_sys", 5),
- list(response_value = 80.540705, concentration = 0.003173),
+ list(response_value = 85.41198, concentration = 0.005502182),
tolerance = 1e-4
)
expect_equal(
ec(model, "effect_tox_sys", 1),
- list(response_value = 83.931891, concentration = 0.000057),
+ list(response_value = 89.00828, concentration = 8.53288e-05),
tolerance = 1e-4
)
expect_equal(
ec(model, "effect_tox_env_sys", 50),
- list(response_value = 12.155592, concentration = 0.905848),
+ list(response_value = 14.67725, concentration = 1.299516),
tolerance = 1e-4
)
expect_equal(
ec(model, "effect_tox_env_sys", 10),
- list(response_value = 21.880066, concentration = 0.0007604),
+ list(response_value = 26.41904, concentration = 0.0008571244),
tolerance = 1e-4
)
expect_equal(
ec(model, "effect_tox_env_sys", 5),
- list(response_value = 23.095625, concentration = 0.00010044),
+ list(response_value = 27.88677, concentration = 0.0001044245),
tolerance = 1e-4
)
expect_equal(
ec(model, "effect_tox_env_sys", 1),
- list(response_value = 24.068073, concentration = 0.000002),
+ list(response_value = 29.06095, concentration = 1.388112e-06),
tolerance = 1e-4
)
})
diff --git a/tests/testthat/test-ecxsys.R b/tests/testthat/test-ecxsys.R
index 1469b4755e4f95cd966ee871feba1ac7510b22c4..0120163817762a933746afd5ef24f755108c76ab 100644
--- a/tests/testthat/test-ecxsys.R
+++ b/tests/testthat/test-ecxsys.R
@@ -1,12 +1,8 @@
-context("ecxsys")
-
-
-# Create model here for use in the various tests to save typing:
mod <- ecxsys(
- concentration = c(0, 0.03, 0.3, 3, 10),
- hormesis_concentration = 0.3,
- effect_tox_observed = c(85, 76, 94, 35, 0),
- effect_tox_env_observed = c(24, 23, 32, 0, 0)
+ concentration = c(0, 0.05, 0.5, 5, 30),
+ hormesis_concentration = 0.5,
+ effect_tox_observed = c(90, 81, 92, 28, 0),
+ effect_tox_env_observed = c(29, 27, 33, 5, 0)
)
@@ -14,18 +10,18 @@ test_that("error when hormesis_concentration not in concentration", {
errorm <- "hormesis_concentration must equal one of the concentration values."
expect_error(
ecxsys(
- concentration = c(0, 0.03, 0.3, 3, 10),
+ concentration = c(0, 0.05, 0.5, 5, 30),
hormesis_concentration = 0.4,
- effect_tox_observed = c(85, 76, 94, 35, 0)
+ effect_tox_observed = c(90, 81, 92, 28, 0)
),
errorm,
fixed = TRUE
)
expect_error(
ecxsys(
- concentration = c(0, 0.03, 0.3, 3, 10),
- hormesis_concentration = 30,
- effect_tox_observed = c(85, 76, 94, 35, 0)
+ concentration = c(0, 0.05, 0.5, 5, 30),
+ hormesis_concentration = 20,
+ effect_tox_observed = c(90, 81, 92, 28, 0)
),
errorm,
fixed = TRUE
@@ -40,7 +36,7 @@ test_that("error when hormesis_index <= 2 or >= (length(concentration))", {
)
expect_error(
ecxsys(
- concentration = c(0, 0.03, 0.3, 3, 10),
+ concentration = c(0, 0.05, 0.5, 5, 30),
hormesis_concentration = 0,
effect_tox_observed = c(85, 76, 94, 35, 0)
),
@@ -49,8 +45,8 @@ test_that("error when hormesis_index <= 2 or >= (length(concentration))", {
)
expect_error(
ecxsys(
- concentration = c(0, 0.03, 0.3, 3, 10),
- hormesis_concentration = 0.03,
+ concentration = c(0, 0.05, 0.5, 5, 30),
+ hormesis_concentration = 0.05,
effect_tox_observed = c(85, 76, 94, 35, 0)
),
errorm,
@@ -58,8 +54,8 @@ test_that("error when hormesis_index <= 2 or >= (length(concentration))", {
)
expect_error(
ecxsys(
- concentration = c(0, 0.03, 0.3, 3, 10),
- hormesis_concentration = 10,
+ concentration = c(0, 0.05, 0.5, 5, 30),
+ hormesis_concentration = 30,
effect_tox_observed = c(85, 76, 94, 35, 0)
),
errorm,
@@ -76,53 +72,53 @@ test_that("min(concentration) == 0 is shifted the correct amount", {
test_that("the discrete results have not changed", {
expect_equal(
mod$effect_tox_LL5,
- c(85, 84.5420, 78.9464, 31.9586, 2.6096),
+ c(90, 89.745092, 82.340325, 26.787104, 4.306868),
tolerance = 1e-4
)
expect_equal(
mod$effect_tox,
- c(100, 99.6760, 94.3198, 34.8601, 0.8400),
+ c(100, 99.455327884, 92.000028403, 27.999995346, 0.002452598),
tolerance = 1e-4
)
expect_equal(
mod$stress_tox,
- c(0, 0.0784, 0.2067, 0.5794, 0.8927),
+ c(0, 0.09296422, 0.23401450, 0.61804415, 0.98352101),
tolerance = 1e-4
)
expect_equal(
mod$sys_tox_not_fitted,
- c(0.2965, 0.2795, 0, 0, 0),
+ c(0.2541156, 0.2324304, 0, 0, 0),
tolerance = 1e-4
)
expect_equal(
mod$effect_tox_sys,
- c(84.7912, 77.4811, 93.1380, 34.8602, 0.8400),
+ c(89.920669034, 81.890300846, 90.863799559, 27.999995346, 0.002452598),
tolerance = 1e-4
)
- expect_equal(mod$stress_env, 0.3885, tolerance = 1e-3)
+ expect_equal(mod$stress_env, 0.3556369, tolerance = 1e-4)
expect_equal(
mod$effect_tox_env_LL5,
- c(26.3333, 26.3333, 26.0620, 0.7104, 0.0374),
+ c(29.6666667, 29.6666657, 29.5214959, 5.4076411, 0.7871201),
tolerance = 1e-4
)
expect_equal(
mod$stress_tox_env,
- c(0.3884, 0.4669, 0.5952, 0.9679, 1.2812),
+ c(0.3556369, 0.4486011, 0.5896514, 0.9736810, 1.3391579),
tolerance = 1e-4
)
expect_equal(
mod$effect_tox_env,
- c(70.8786, 56.4138, 32, 0.0202, 0),
+ c(76.36558967, 59.90195211, 33, 0.01079038, 0),
tolerance = 1e-4
)
expect_equal(
mod$sys_tox_env_not_fitted,
- c(0.2537, 0.1815, 0, 0, 0),
+ c(0.2566005, 0.1753250, 0, 0, 0),
tolerance = 1e-4
)
expect_equal(
mod$effect_tox_env_sys,
- c(24.3254, 22.3232, 29.3860, 0.0202, 0),
+ c(29.37633973, 25.99861674, 30.16676004, 0.01079038, 0),
tolerance = 1e-4
)
})
@@ -132,20 +128,20 @@ test_that("the curves have not changed", {
new_curves <- mod$curves[c(1, 714, 810, 905, 1000), ] # random indices
rownames(new_curves) <- NULL
reference_curves <- data.frame(
- concentration = c(0.0000, 0.0137, 0.1253, 1.1196, 10.0000),
- effect_tox_LL5 = c(85.0000, 84.8116, 82.6996, 61.2882, 2.6096),
- effect_tox_env_LL5 = c(26.3333, 26.3333, 26.3289, 7.5384, 0.0374),
- effect_tox = c(100.0000, 99.8787, 98.0645, 73.6652, 0.8400),
- stress_tox = c(0.0001, 0.0570, 0.1421, 0.3721, 0.8927),
- sys_tox = c(0.2980, 0.2887, 0.1336, 0.0000, 0.0000),
- stress_tox_sys = c(0.2981, 0.3457, 0.2757, 0.3721, 0.8927),
- effect_tox_sys = c(84.7797, 77.9323, 87.5799, 73.6652, 0.8400),
- stress_env = c(0.3885, 0.3885, 0.3885, 0.3885, 0.3885),
- stress_tox_env = c(0.3886, 0.4455, 0.5306, 0.7606, 1.2812),
- effect_tox_env = c(70.8634, 60.4957, 44.0834, 8.5137, 0.0000),
- sys_tox_env = c(0.2516, 0.2175, 0.0762, 0.0000, 0.0000),
- stress_tox_env_sys = c(0.6402, 0.6630, 0.6068, 0.7606, 1.2812),
- effect_tox_env_sys = c(24.3112, 20.7272, 29.9481, 8.5133, 0.0000),
+ concentration = c(0.00000, 0.03004, 0.30512, 3.02551, 30.00000),
+ effect_tox_LL5 = c(90.00000, 89.88245, 86.18591, 40.42560, 4.30687),
+ effect_tox_env_LL5 = c(29.66667, 29.66667, 29.65530, 9.27616, 0.78712),
+ effect_tox = c(100.00000, 99.70167, 95.45944, 49.54173, 0.00245),
+ stress_tox = c(0.00013, 0.07631, 0.19093, 0.50236, 0.98352),
+ sys_tox = c(0.25487, 0.24017, 0.05667, 0.00000, 0.00000),
+ stress_tox_sys = c(0.25499, 0.31648, 0.24760, 0.50236, 0.98352),
+ effect_tox_sys = c(89.90735, 82.27932, 90.67515, 49.54173, 0.00245),
+ stress_env = c(0.35564, 0.35564, 0.35564, 0.35564, 0.35564),
+ stress_tox_env = c(0.35576, 0.43195, 0.54657, 0.85800, 1.33916),
+ effect_tox_env = c(76.34546, 63.03379, 41.01651, 1.93183, 0.00000),
+ sys_tox_env = c(0.25443, 0.20496, 0.04445, 0.00000, 0.00000),
+ stress_tox_env_sys = c(0.61020, 0.63691, 0.59102, 0.85800, 1.33916),
+ effect_tox_env_sys = c(29.35449, 24.84147, 32.75385, 1.93182, 0),
use_for_plotting = c(TRUE, TRUE, TRUE, TRUE, TRUE)
)
class(new_curves) <- class(reference_curves)
@@ -163,10 +159,10 @@ test_that("the returned fn works the same way as internally", {
test_that("function arguments are returned unchanged", {
args_reference <- list(
- concentration = c(0, 0.03, 0.3, 3, 10),
- hormesis_concentration = 0.3,
- effect_tox_observed = c(85, 76, 94, 35, 0),
- effect_tox_env_observed = c(24, 23, 32, 0, 0),
+ concentration = c(0, 0.05, 0.5, 5, 30),
+ hormesis_concentration = 0.5,
+ effect_tox_observed = c(90, 81, 92, 28, 0),
+ effect_tox_env_observed = c(29, 27, 33, 5, 0),
effect_max = 100,
p = 3.2,
q = 3.2
@@ -175,22 +171,29 @@ test_that("function arguments are returned unchanged", {
})
+ecxsys(
+ concentration = c(0, 0.05, 0.5, 5, 30),
+ hormesis_concentration = 0.5,
+ effect_tox_observed = c(90, 81, 92, 28, 0),
+ effect_tox_env_observed = c(29, 27, 33, 5, 0)
+)
+
test_that("results are independent of concentration shift", {
mod_2 <- ecxsys(
- concentration = c(0, 0.03, 0.3, 3, 10) * 2,
- hormesis_concentration = 0.3 * 2,
- effect_tox_observed = c(85, 76, 94, 35, 0),
- effect_tox_env_observed = c(24, 23, 32, 0, 0)
+ concentration = c(0, 0.05, 0.5, 5, 30) * 2,
+ hormesis_concentration = 0.5 * 2,
+ effect_tox_observed = c(90, 81, 92, 28, 0),
+ effect_tox_env_observed = c(29, 27, 33, 5, 0)
)
expect_equal(mod$effect_tox_sys, mod_2$effect_tox_sys)
expect_equal(mod$effect_tox_env_sys, mod_2$effect_tox_env_sys)
mod_10 <- ecxsys(
- concentration = c(0, 0.03, 0.3, 3, 10) * 10,
- hormesis_concentration = 0.3 * 10,
- effect_tox_observed = c(85, 76, 94, 35, 0),
- effect_tox_env_observed = c(24, 23, 32, 0, 0)
+ concentration = c(0, 0.05, 0.5, 5, 30) * 10,
+ hormesis_concentration = 0.5 * 10,
+ effect_tox_observed = c(90, 81, 92, 28, 0),
+ effect_tox_env_observed = c(29, 27, 33, 5, 0)
)
- # Concentration shifts by factors other than powers of 10 will affect
+ # Concentration shifts by factors other than powers of 10 may affect
# the result because of the way the zero concentration is "corrected".
expect_equal(mod$curves$effect_tox,
mod_10$curves$effect_tox)
@@ -203,9 +206,9 @@ test_that("results are independent of concentration shift", {
test_that("effect_tox_env_observed can be left out", {
mod_without_env <- ecxsys(
- concentration = c(0, 0.03, 0.3, 3, 10),
- hormesis_concentration = 0.3,
- effect_tox_observed = c(85, 76, 94, 35, 0)
+ concentration = c(0, 0.05, 0.5, 5, 30),
+ hormesis_concentration = 0.5,
+ effect_tox_observed = c(90, 81, 92, 28, 0)
)
expect_equal(mod$effect_tox_sys, mod_without_env$effect_tox_sys)
expect_equal(mod$curves$effect_tox_sys,
@@ -226,4 +229,17 @@ test_that("sys model not converging produces a warning, not an error", {
),
fixed = TRUE
)
+ expect_warning(
+ ecxsys(
+ concentration = c(0, 0.1, 0.5, 1, 10, 33),
+ effect_tox_observed = c(90, 81, 92, 50, 18, 0),
+ effect_tox_env_observed = c(75, 89, 54, 7, 0, 0),
+ hormesis_concentration = 0.5
+ ),
+ paste(
+ "Using a horizontal linear model for sys_tox_env_mod because the",
+ "Weibull model did not converge."
+ ),
+ fixed = TRUE
+ )
})
diff --git a/tests/testthat/test-options.R b/tests/testthat/test-options.R
index dc5945baeefcc9cf31aa32e387837fcba6f6bea3..6d7c85d1db31a21e9bd235a7c173fcb26fd6b8e1 100644
--- a/tests/testthat/test-options.R
+++ b/tests/testthat/test-options.R
@@ -1,6 +1,3 @@
-context("options")
-
-
test_that("user options are not permanently changed by ecxsys()", {
# drc::drm() messes with the options and fails to return them to their
# previous values. And options are temporarily modified in ecxsys(). This
@@ -13,10 +10,10 @@ test_that("user options are not permanently changed by ecxsys()", {
original_options <- options()
model <- ecxsys(
- concentration = c(0, 0.03, 0.3, 3, 10),
- effect_tox_observed = c(85, 76, 94, 35, 0),
- effect_tox_env_observed = c(24, 23, 32, 0, 0),
- hormesis_concentration = 0.3
+ concentration = c(0, 0.05, 0.5, 5, 30),
+ hormesis_concentration = 0.5,
+ effect_tox_observed = c(90, 81, 92, 28, 0),
+ effect_tox_env_observed = c(29, 27, 33, 5, 0)
)
new_options <- options()
diff --git a/tests/testthat/test-plot_ecxsys.R b/tests/testthat/test-plot_ecxsys.R
index 541028b5a21b5770e6e417476fac6ffc6e3391b6..4a61b99d27c535b04eef909476a1bf6e128d6232 100644
--- a/tests/testthat/test-plot_ecxsys.R
+++ b/tests/testthat/test-plot_ecxsys.R
@@ -1,6 +1,3 @@
-context("plots")
-
-
test_that("log ticks are correct", {
x <- c(0.03, 0.3, 3, 30)
ticks <- get_log_ticks(x)
diff --git a/tests/testthat/test-predict_mixture.R b/tests/testthat/test-predict_mixture.R
new file mode 100644
index 0000000000000000000000000000000000000000..3471aeb429a0e7cff5f1496d08445e36253eb99f
--- /dev/null
+++ b/tests/testthat/test-predict_mixture.R
@@ -0,0 +1,34 @@
+model_1 <- ecxsys(
+ concentration = c(0, 0.05, 0.5, 5, 30),
+ hormesis_concentration = 0.5,
+ effect_tox_observed = c(90, 81, 92, 28, 0),
+ effect_tox_env_observed = c(29, 27, 33, 5, 0),
+ effect_max = 101
+)
+model_2 <- ecxsys(
+ concentration = c(0, 0.1, 1, 10, 100, 1000),
+ hormesis_concentration = 10,
+ effect_tox_observed = c(26, 25, 24, 27, 5, 0),
+ effect_max = 30
+)
+
+
+test_that("results have not changed", {
+ new <- predict_mixture(model_1, model_2, c(0, 0.01, 0.1, 1, 7, 15), 5, 0.3)
+ reference <- c(89.460324, 85.205168, 81.440100, 57.297855, 2.911545, 0)
+ expect_equal(new, reference, tolerance = 1e-5)
+})
+
+
+test_that("predictions are symmetric", {
+ conc_1 <- c(0, 10^seq(log10(0.001), log10(40), length.out = 50))
+ conc_2 <- 3.5
+ prop_ca <- 0.8
+ effect_12 <- predict_mixture(model_1, model_2, conc_1, conc_2, prop_ca)
+ effect_21 <- sapply(
+ conc_1,
+ function(x) predict_mixture(model_2, model_1, conc_2, x, prop_ca)
+ )
+ effect_21 <- effect_21 / model_2$args$effect_max * model_1$args$effect_max
+ expect_equal(effect_12, effect_21)
+})
diff --git a/tests/testthat/test-stress_effect_conversion.R b/tests/testthat/test-stress_effect_conversion.R
index 506491d38c69d6d246d45c549ae8ab2d86b625b4..b583ece438293fb19df092d467a5b016054c6c98 100644
--- a/tests/testthat/test-stress_effect_conversion.R
+++ b/tests/testthat/test-stress_effect_conversion.R
@@ -1,6 +1,3 @@
-context("stress effect conversion")
-
-
test_that("stress_to_effect handles bad input", {
expect_equal(stress_to_effect(-3), 1)
expect_equal(stress_to_effect(10), 0)