From 296ce3cbed5b670ffd2bffdbc3fbcfbba68e4753 Mon Sep 17 00:00:00 2001
From: Sebastian Henz <sebastian.henz@ufz.de>
Date: Sat, 15 Feb 2020 14:29:49 +0100
Subject: [PATCH] Major rework of the source code and documentation. Refactor
everything to improve usability and to make debugging and extending the
package easier. See NEWS.md for details.
---
DESCRIPTION | 4 +-
NAMESPACE | 3 +-
NEWS.md | 17 +
R/ec.R | 104 ++--
R/ecxsys.R | 521 +++++++-----------
R/helpers.R | 59 ++
R/moving_weighted_mean.R | 16 -
R/plot_ecxsys.R | 67 +--
R/plot_effect.R | 102 ++--
R/{plot_system_stress.R => plot_stress.R} | 36 +-
R/predict_ecxsys.R | 122 ++++
R/stress_effect_conversion.R | 15 +-
R/stressaddition-package.R | 9 +-
README.md | 9 +-
man/Stressconversion.Rd | 11 +-
man/ec.Rd | 49 +-
man/ecxsys.Rd | 75 +--
man/plot_ecxsys.Rd | 21 +-
man/predict_ecxsys.Rd | 59 ++
tests/testthat/test-ec.R | 75 ++-
tests/testthat/test-ecxsys.R | 250 ++++-----
.../testthat/test-stress_effect_conversion.R | 8 +-
22 files changed, 856 insertions(+), 776 deletions(-)
create mode 100644 R/helpers.R
delete mode 100644 R/moving_weighted_mean.R
rename R/{plot_system_stress.R => plot_stress.R} (62%)
create mode 100644 R/predict_ecxsys.R
create mode 100644 man/predict_ecxsys.Rd
diff --git a/DESCRIPTION b/DESCRIPTION
index 45b797f..6ca5088 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
Package: stressaddition
Type: Package
Title: Modeling Tri-Phasic Concentration-Response Relationships
-Version: 1.11.1
-Date: 2020-02-13
+Version: 2.0.0
+Date: 2020-02-15
Authors@R: c(person("Sebastian",
"Henz",
role = c("aut", "cre"),
diff --git a/NAMESPACE b/NAMESPACE
index 569bc2a..06741e7 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -4,7 +4,8 @@ export(ec)
export(ecxsys)
export(effect_to_stress)
export(plot_effect)
-export(plot_system_stress)
+export(plot_stress)
+export(predict_ecxsys)
export(stress_to_effect)
import(grDevices)
import(graphics)
diff --git a/NEWS.md b/NEWS.md
index e7d2e27..76e2ec2 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,20 @@
+# stressaddition 2.0.0
+
+* Changed the order of arguments in `ecxsys()`.
+* Removed `hormesis_index` argument from `ecxsys()`. Use `hormesis_concentration` instead.
+* New function `predict_ecxsys()` replaces `fn()` from the `ecxsys()` output.
+* Renamed the arguments in `ec()`.
+* Made `ec()` more flexible. It now also accepts a data.frame with a concentration column and a column of response values.
+* Added LL5 curves to the legend of `plot_effect()`.
+* Replaced every occurrence of "simple" in variable names with "LL5".
+* Replaced every occurrence of "sys_stress" in variable names with "sys" because the extra "stress" was redundant.
+* Renamed `plot_system_stress()` to `plot_stress()` because it is planned to plot more stresses with this function in a future update.
+* Changed the order of the columns in the output of `predict_ecxsys()`.
+* Improved the internal structure of the package.
+* Improved the tests.
+* Improved the documentation.
+
+
# stressaddition 1.11.1
* First public version.
diff --git a/R/ec.R b/R/ec.R
index 650e249..de29753 100644
--- a/R/ec.R
+++ b/R/ec.R
@@ -1,56 +1,94 @@
#' Effect Concentrations
#'
-#' Estimate the concentration to reach a certain effect relative to the control.
+#' Estimate the concentration to reach a certain effect or stress level relative
+#' to the control.
#'
#' If the response level occurs multiple times because of hormesis, which may
-#' happen for low values of \code{target_effect}, then the occurrence with the
+#' happen for low values of \code{response_level}, then the occurrence with the
#' smallest concentration is returned.
#'
-#' @param model The object returned from \code{ecxsys()}.
-#' @param effect_name The name of the effect for which you want to calculate the
-#' EC. Must be the name of a column in \code{model$curves}.
-#' @param target_effect The desired effect percentage between 0 and 100. For
-#' example with the value 10 the function will return the EC10, i.e. the
-#' concentration where the response falls below 90 \% of the maximum possible
+#' @param model This can be one of three types of objects: Either the output of
+#' \code{\link{ecxsys}} or the output of \code{\link{predict_ecxsys}} or a
+#' data frame with a "concentration" column and a \code{response_name} column.
+#' See the examples.
+#' @param response_name The name of the effect or stress for which you want to
+#' calculate the EC. Must be one of \code{colnames(model$curves)}.
+#' @param response_level The desired response level as a percentage between 0
+#' and 100. For example with the value 10 the function will return the EC10,
+#' i.e. the concentration where the response falls below 90 \% of the maximum
#' response.
#'
-#' @return A list containing the effect concentration and the corresponding
-#' effect.
+#' @return A list containing the response concentration and the corresponding
+#' response value.
#'
-#' @examples model <- ecxsys(
+#' @examples # Calculate the EC_10, the concentration where the effect falls
+#' # 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),
-#' effect_tox_env_observed = c(24, 23, 32, 0, 0),
#' hormesis_concentration = 0.3
#' )
-#' # Calculate the EC_10, the concentration where the effect falls
-#' # below 90 % of the effect in the control:
-#' ec_10 <- ec(model, "effect_tox_sys", 10)
+#'
+#' # using the ecxsys() output or the curves therein directly:
+#' ec(model, "effect_tox_sys", 10)
+#' ec(model$curves, "effect_tox_sys", 10)
+#'
+#' # using the output of predict_ecxsys() with custom concentrations:
+#' conc <- 10^seq(-9, 1, length.out = 1000)
+#' curves <- predict_ecxsys(model, conc)
+#' ec(curves, "effect_tox_sys", 10)
+#'
+#' # using a custom data frame:
+#' df_custom <- data.frame(
+#' concentration = curves$concentration,
+#' foo = curves$effect_tox_sys
+#' )
+#' ec(df_custom, "foo", 10)
#'
#' @export
-ec <- function(model, effect_name, target_effect) {
+ec <- function(model, response_name, response_level) {
+ if (inherits(model, "drc")) {
+ stop("Please use drc::ED for drc objects.")
+ }
+
stopifnot(
- is.character(effect_name),
- length(effect_name) == 1,
- effect_name %in% names(model$curves),
- effect_name != "concentration",
- target_effect > 0,
- target_effect < 100
+ is.character(response_name),
+ length(response_name) == 1,
+ response_name != "concentration",
+ response_level > 0,
+ response_level < 100
)
- effect <- model$curves[, effect_name]
- concentration <- model$curves$concentration
- control_effect <- effect[1]
- if (control_effect == 0) {
- stop("Reference effect is zero, calculation of EC not possible.")
+
+ if (!inherits(model, c("ecxsys", "ecxsys_predicted"))) {
+ if (is.data.frame(model)) {
+ concentration <- model$concentration
+ response <- model[, response_name]
+ } else {
+ stop("Invalid first argument.")
+ }
+ } else if (inherits(model, "ecxsys")) {
+ concentration <- model$curves$concentration
+ response <- model$curves[, response_name]
+ } else if (inherits(model, "ecxsys_predicted")) {
+ concentration <- model$concentration
+ response <- model[, response_name]
+ }
+
+ reference <- response[1]
+ if (reference == 0) {
+ stop("Reference value is zero, calculation of EC not possible.")
}
- target_effect <- (1 - target_effect / 100) * control_effect
- output <- list(effect = target_effect)
- # Get the index of where the effect changes from above to below
- # target_effect:
- below <- which(effect < target_effect)[1]
+ response_level <- (1 - response_level / 100) * reference
+ output <- list(response_value = response_level)
+
+ # Get the index of where the response changes from above to below
+ # response_level:
+ below <- which(response < response_level)[1]
above <- below - 1
+
# linear interpolation
- dist <- (target_effect - effect[below]) / (effect[above] - effect[below])
+ dist <- (response_level - response[below]) / (response[above] - response[below])
output$concentration <- dist *
(concentration[above] - concentration[below]) + concentration[below]
output
diff --git a/R/ecxsys.R b/R/ecxsys.R
index 4273586..5919120 100644
--- a/R/ecxsys.R
+++ b/R/ecxsys.R
@@ -17,106 +17,58 @@
#' of 0 at ten times the maximum observed concentration.
#'
#' The vectors \code{concentration}, \code{effect_tox_observed} and
-#' \code{effect_tox_env_observed} must be of equal length and should be sorted
-#' by increasing concentration.
+#' \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 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
#' concentrations and in absence of environmental stress. Values must be
#' between 0 and \code{effect_max}.
#' @param effect_tox_env_observed Effect values observed in the presence of
-#' environmental stress. This argument is optional and can be left out to
-#' model without environmental stress. Values must be between 0 and
-#' \code{effect_max}.
-#' @param hormesis_concentration The concentration where the hormesis occurs.
-#' This is usually the concentration of the highest effect after the control.
-#' @param hormesis_index \strong{Deprecated, will be removed soon.} A single
-#' integer specifying the index of the hormesis concentration in the
-#' concentration vector. This argument exists for compatibility with older
-#' versions of this function.
+#' environmental stress. Must be between 0 and \code{effect_max}.
#' @param effect_max The maximum value the effect could possibly reach. For
#' 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 The result is a list containing many different objects with the most
-#' important being \code{curves} and \code{fn}. You can use \code{fn()} to
-#' calculate the curves at your concentrations of choice, see examples.
-#' \code{curves} is a data frame 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_simple}{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{effect_tox_sys}{Modeled effect resulting from toxicant and system
-#' stress.}
-#' \item{stress_tox}{The toxicant stress.}
-#' \item{sys_stress_tox}{System stress under toxicant stress conditions
-#' without environmental stress.}
-#' \item{stress_tox_sys}{The sum of \code{stress_tox} and
-#' \code{sys_stress_tox}.}
-#' \item{effect_tox_env_simple}{The five-parameter log-logistic model of the
-#' effect 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,
-#' environmental and system stress.}
-#' \item{stress_env}{Environmental stress.}
-#' \item{stress_tox_env}{The sum of toxicant and environmental stress.}
-#' \item{sys_stress_tox_env}{System stress under toxicant and
-#' environmental stress conditions.}
-#' \item{stress_tox_env_sys}{The sum of \code{stress_tox_env} and
-#' \code{sys_stress_tox_env}.}
-#' \item{use_for_plotting}{A boolean vector which is used in the plotting
-#' functions. It controls which parts of the curves are removed for the
-#' broken concentration axis.}
-#' }
+#' @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.
#'
#' @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),
-#' hormesis_concentration = 0.3
+#' effect_tox_env_observed = c(24, 23, 32, 0, 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:
#' model <- ecxsys(
#' concentration = c(0, 0.03, 0.3, 3, 10),
+#' 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),
-#' hormesis_concentration = 0.3,
#' effect_max = 20
#' )
#'
-#' # The returned object contains a function which is useful for calculating
-#' # effect and stress values at specific concentrations:
-#' model$fn(c(0, 0.01, 0.1, 1))
-#'
#' @export
ecxsys <- function(concentration,
- effect_tox_observed,
- effect_tox_env_observed,
hormesis_concentration,
- hormesis_index,
+ effect_tox_observed,
+ effect_tox_env_observed = NULL,
effect_max = 100,
- #stress_tox_at_hormesis = NULL,
p = 3.2,
q = 3.2) {
output <- list(args = as.list(environment()))
+ class(output) <- c("ecxsys", class(output))
original_options <- options()
on.exit(reset_options(original_options))
- warn_error_original <- list(
- warn = getOption("warn"),
- show.error.messages = getOption("show.error.messages")
- )
# input validation ----------------------------------------------------
@@ -129,28 +81,18 @@ ecxsys <- function(concentration,
if (length(concentration) > length(unique(concentration))) {
stop("Concentrations must be unique.")
}
-
- m_hc <- missing(hormesis_concentration)
- m_hi <- missing(hormesis_index)
- if (m_hc && m_hi) {
- stop("Pleace specify either hormesis_concentration or hormesis_index.")
- } else if (!m_hi && !m_hc) {
- stop("Use either hormesis_concentration or hormesis_index but not both.")
- } else if (!m_hc) {
- if (!hormesis_concentration %in% concentration) {
- stop("hormesis_concentration must be one of the values ",
- "in concentration.")
- }
- hormesis_index = which(hormesis_concentration == concentration)
+ if (length(hormesis_concentration) != 1) {
+ stop("hormesis_concentration must be a single number.")
}
-
- if (length(hormesis_index) != 1) {
- stop("hormesis_index must be a single integer.")
- } else if (hormesis_index <= 2 || hormesis_index >= length(concentration)) {
- stop("hormesis_index must be greater than 2 and less than ",
- "(length(concentration)).")
+ if (!hormesis_concentration %in% concentration) {
+ stop("hormesis_concentration must equal one of the concentration values.")
+ }
+ hormesis_index = which(hormesis_concentration == concentration)
+ if (hormesis_index <= 2 || hormesis_index >= length(concentration)) {
+ stop("hormesis_concentration must be greater than the lowest ",
+ "non-control concentration and less than the highest concentration")
}
- if (missing(effect_tox_env_observed)) {
+ if (is.null(effect_tox_env_observed)) {
with_env <- FALSE
effect_tox_env_observed <- rep(NA, length(concentration))
} else {
@@ -176,16 +118,16 @@ ecxsys <- function(concentration,
conc_shift <- 2 # Powers of ten to shift the control downwards from the
# second lowest concentration. This is required to approximate 0 because
# of the logarithmic axis.
+ if (is.unsorted(concentration)) {
+ stop("The values must be sorted by increasing concentration.")
+ }
if (any(concentration < 0)) {
stop("Concentrations must be >= 0")
- } else if (min(concentration) > 0) {
- stop("No control is given. The first concentration must be 0.")
- } else {
- min_conc <- 10 ^ floor(log10(concentration[2]) - conc_shift)
}
- if (is.unsorted(concentration)) {
- stop("The values must be sorted by increasing concentration.")
+ if (min(concentration) > 0) {
+ stop("No control is given. The first concentration must be 0.")
}
+ min_conc <- 10 ^ floor(log10(concentration[2]) - conc_shift)
# scale observed effect -----------------------------------------------
@@ -195,299 +137,128 @@ ecxsys <- function(concentration,
effect_tox_env_observed <- effect_tox_env_observed / effect_max
- # prepare adjusted control concentration ------------------------------
- conc_adjust_factor <- 10^-5
- output$conc_adjust_factor <- conc_adjust_factor
-
-
- # traditional simple model (LL.5) -------------------------------------
- conc_with_control_shifted <- c(min_conc, concentration[-1])
- effect_tox_observed_averaged <- moving_weighted_mean(effect_tox_observed)
- effect_tox_env_observed_averaged <- moving_weighted_mean(effect_tox_env_observed)
- temp <- approx(
- log10(conc_with_control_shifted),
- effect_tox_observed_averaged,
- n = 10
- )
- conc_interpolated <- 10^temp$x
- effect_tox_observed_interpolated_simple_model <- temp$y
- effect_tox_mod_simple <- drc::drm(
- effect_tox_observed_interpolated_simple_model ~ conc_interpolated,
- fct = drc::LL.5(fixed = c(
- NA, 0, effect_tox_observed_averaged[1], NA, NA
- ))
- )
- options(warn_error_original)
- effect_tox_simple <- predict(effect_tox_mod_simple, data.frame(concentration))
- output$effect_tox_mod_simple <- effect_tox_mod_simple
- output$effect_tox_simple <- effect_tox_simple * effect_max
+ # traditional log-logistic model (LL.5) -------------------------------
+ LL5_tox <- fit_LL5_model(min_conc, concentration,
+ effect_tox_observed,
+ original_options)
+ output$effect_tox_LL5_mod <- LL5_tox$effect_LL5_mod
+ output$effect_tox_LL5 <- LL5_tox$effect_LL5 * effect_max
if (with_env) {
- effect_tox_env_observed_interpolated_simple_model <- approx(
- log10(conc_with_control_shifted),
- effect_tox_env_observed_averaged,
- xout = temp$x
- )$y
- effect_tox_env_mod_simple <- drc::drm(
- effect_tox_env_observed_interpolated_simple_model ~ conc_interpolated,
- fct = drc::LL.5(fixed = c(
- NA, 0, effect_tox_env_observed_averaged[1], NA, NA
- ))
- )
- options(warn_error_original)
- effect_tox_env_simple <- predict(
- effect_tox_env_mod_simple,
- data.frame(concentration)
- )
- output$effect_tox_env_mod_simple <- effect_tox_env_mod_simple
- output$effect_tox_env_simple <- effect_tox_env_simple * effect_max
+ LL5_tox_env <- fit_LL5_model(min_conc, concentration,
+ effect_tox_env_observed,
+ original_options)
+ output$effect_tox_env_LL5_mod <- LL5_tox_env$effect_LL5_mod
+ output$effect_tox_env_LL5 <- LL5_tox_env$effect_LL5 * effect_max
}
# interpolation between subhormesis and hormesis ----------------------
n_new <- 3 # number of new points
- len <- n_new + 2 # Add 2 because seq() includes the left and right end.
- subhormesis_index <- hormesis_index - 1
-
- concentration_interpolated <- 10^seq(
- log10(concentration[subhormesis_index]),
- log10(concentration[hormesis_index]),
- length.out = len
- )
- concentration <- append(
- concentration,
- concentration_interpolated[-c(1, len)],
- subhormesis_index
- )
-
- effect_tox_observed_interpolated <- seq(
- effect_tox_observed[subhormesis_index],
- effect_tox_observed[hormesis_index],
- length.out = len
- )
- effect_tox_observed <- append(
- effect_tox_observed,
- effect_tox_observed_interpolated[-c(1, len)],
- subhormesis_index
- )
-
+ concentration <- interpolate(concentration, hormesis_index, n_new, TRUE)
+ effect_tox_observed <- interpolate(effect_tox_observed, hormesis_index, n_new)
if (with_env) {
- effect_tox_env_observed_interpolated <- seq(
- effect_tox_env_observed[subhormesis_index],
- effect_tox_env_observed[hormesis_index],
- length.out = len
- )
- effect_tox_env_observed <- append(
- effect_tox_env_observed,
- effect_tox_env_observed_interpolated[-c(1, len)],
- subhormesis_index
- )
+ effect_tox_env_observed <- interpolate(effect_tox_env_observed,
+ hormesis_index, n_new)
}
-
hormesis_index <- hormesis_index + n_new
-
# In the output return only the values at the original concentrations
# and exclude those at the interpolated concentrations:
- exclude <- seq(subhormesis_index + 1, hormesis_index - 1)
- keep <- !seq_along(concentration) %in% exclude
+ keep <- !seq_along(concentration) %in% seq(hormesis_index - n_new, hormesis_index - 1)
# effect_tox ----------------------------------------------------------
- effect_tox <- effect_tox_observed
- effect_tox[1] <- 1
- effect_to_fit_idx <- 2:(hormesis_index - 1)
- effect_tox[effect_to_fit_idx] <- NA
- effect_tox_mod <- drc::drm(
- effect_tox ~ concentration,
- fct = drc::W1.2()
- )
- options(warn_error_original)
- effect_tox <- predict(
- effect_tox_mod,
- data.frame(concentration = concentration)
- )
+ effect_tox <- c(1, effect_tox_observed[-1])
+ effect_tox[2:(hormesis_index - 1)] <- NA
+ effect_tox_mod <- drc::drm(effect_tox ~ concentration, fct = drc::W1.2())
+ options(original_options["warn"])
+ effect_tox <- predict(effect_tox_mod, data.frame(concentration = concentration))
output$effect_tox_mod <- effect_tox_mod
output$effect_tox <- effect_tox[keep] * effect_max
-
- # system stress without environmental stress --------------------------
- stress_tox_observed <- effect_to_stress(
- effect_tox_observed, p, q
- )
stress_tox <- effect_to_stress(effect_tox, p, q)
output$stress_tox <- stress_tox[keep]
- sys_stress_tox <- stress_tox_observed - stress_tox
- sys_stress_tox <- pmin(pmax(sys_stress_tox, 0), 1)
- sys_stress_tox[hormesis_index:length(sys_stress_tox)] <- 0
- # Add sys_stress_tox to the output before it is fitted:
- output$sys_stress_tox <- sys_stress_tox[keep]
- sys_stress_tox_mod <- tryCatch(
- {
- # There is no other way to suppress that one error message
- # except by changing the options temporarily.
- options(show.error.messages = FALSE)
- drc::drm(sys_stress_tox ~ stress_tox, fct = drc::W1.3())
- },
- error = function(e) {
- options(warn_error_original)
- warning(
- "Using a horizontal linear model for sys_stress_tox_mod ",
- "because the Weibull model did not converge.",
- call. = FALSE
- )
- # Failure to converge often happens when all or almost all sys
- # stress values are zero. Fitting a linear model in this case seems
- # to be the most appropriate remedy.
- stress_tox <- range(stress_tox)
- sys_stress_tox <- c(0, 0)
- return(lm(sys_stress_tox ~ stress_tox))
- },
- finally = options(warn_error_original)
- )
- output$sys_stress_tox_mod <- sys_stress_tox_mod
- sys_stress_tox <- unname(
- predict(sys_stress_tox_mod, data.frame(stress_tox))
+
+
+ # system stress without environmental stress --------------------------
+ fit_sys_output <- fit_sys(
+ effect_to_stress(effect_tox_observed, p, q),
+ stress_tox,
+ stress_tox,
+ hormesis_index,
+ original_options
)
+ output$sys_tox_not_fitted <- fit_sys_output$sys[keep]
+ output$sys_tox_mod <- fit_sys_output$sys_mod
+ if (inherits(fit_sys_output$sys_mod, "lm")) {
+ warning(
+ "Using a horizontal linear model for sys_tox_mod ",
+ "because the Weibull model did not converge."
+ )
+ }
+ sys_tox <- fit_sys_output$sys_modeled
+ output$sys_tox <- sys_tox
# modeled effect without environmental stress -------------------------
- stress_tox_sys <- stress_tox + sys_stress_tox
+ stress_tox_sys <- stress_tox + sys_tox
effect_tox_sys <- stress_to_effect(stress_tox_sys, p, q)
output$effect_tox_sys <- effect_tox_sys[keep] * effect_max
if (with_env) {
# env stress ------------------------------------------------------
- stress_tox_env_observed <- effect_to_stress(
- effect_tox_env_observed, p, q
- )
+ stress_tox_env_observed <- effect_to_stress(effect_tox_env_observed, p, q)
stress_env <- (stress_tox_env_observed - stress_tox)[hormesis_index]
- stress_env <- pmax(stress_env, 0)
+ stress_env <- clamp(stress_env)
output$stress_env <- stress_env
-
- # system stress with environmental stress -------------------------
stress_tox_env <- stress_tox + stress_env
- effect_tox_env <- stress_to_effect(stress_tox_env, p, q)
output$stress_tox_env <- stress_tox_env[keep]
+ effect_tox_env <- stress_to_effect(stress_tox_env, p, q)
output$effect_tox_env <- effect_tox_env[keep] * effect_max
- sys_stress_tox_env <- stress_tox_env_observed - stress_tox_env
- sys_stress_tox_env <- pmin(pmax(sys_stress_tox_env, 0), 1)
- sys_stress_tox_env[hormesis_index:length(sys_stress_tox_env)] <- 0
- # Add sys_stress_tox_env to the output before it is fitted:
- output$sys_stress_tox_env <- sys_stress_tox_env[keep]
- sys_stress_tox_env_mod <- tryCatch(
- {
- # There is no other way to suppress that one error message
- # except by changing the options temporarily.
- options(show.error.messages = FALSE)
- drc::drm(sys_stress_tox_env ~ stress_tox, fct = drc::W1.3())
- },
- error = function(e) {
- options(warn_error_original)
- warning(
- "Using a horizontal linear model for ",
- "sys_stress_tox_env_mod because the Weibull model did ",
- "not converge.",
- call. = FALSE
- )
- # Failure to converge often happens when all or almost all sys
- # stress values are zero. Fitting a linear model in this case
- # seems to be the most appropriate remedy.
- stress_tox <- range(stress_tox)
- sys_stress_tox_env <- c(0, 0)
- return(lm(sys_stress_tox_env ~ stress_tox))
- },
- finally = options(warn_error_original)
- )
- output$sys_stress_tox_env_mod <- sys_stress_tox_env_mod
- sys_stress_tox_env <- unname(
- predict(sys_stress_tox_env_mod, data.frame(stress_tox))
+
+
+ # system stress with environmental stress -------------------------
+ fit_sys_output <- fit_sys(
+ stress_tox_env_observed,
+ stress_tox_env,
+ stress_tox,
+ hormesis_index,
+ original_options
)
+ output$sys_tox_env_not_fitted <- fit_sys_output$sys[keep]
+ output$sys_tox_env_mod <- fit_sys_output$sys_mod
+ if (inherits(fit_sys_output$sys_mod, "lm")) {
+ warning(
+ "Using a horizontal linear model for sys_tox_env_mod ",
+ "because the Weibull model did not converge."
+ )
+ }
+ sys_tox_env <- fit_sys_output$sys_modeled
+ output$sys_tox_env <- sys_tox_env
# modeled effect with environmental stress ------------------------
- stress_tox_env_sys <- stress_tox_env + sys_stress_tox_env
+ stress_tox_env_sys <- stress_tox_env + sys_tox_env
effect_tox_env_sys <- stress_to_effect(stress_tox_env_sys, p, q)
output$effect_tox_env_sys <- effect_tox_env_sys[keep] * effect_max
}
- # building the function -----------------------------------------------
- fn <- function(conc) {
- # This function returns all modeled values at the provided
- # concentrations. conc = a vector of concentrations
- stopifnot(is.numeric(conc))
- effect_tox_simple_fn <- predict(
- effect_tox_mod_simple,
- data.frame(concentration = conc)
- )
- effect_tox_fn <- predict(
- effect_tox_mod,
- data.frame(concentration = conc)
- )
- stress_tox_fn <- effect_to_stress(effect_tox_fn, p, q)
- sys_stress_tox_fn <- predict(
- sys_stress_tox_mod,
- data.frame(stress_tox = stress_tox_fn)
- )
- stress_tox_sys_fn <- stress_tox_fn + sys_stress_tox_fn
- effect_tox_sys_fn <- stress_to_effect(stress_tox_sys_fn, p, q)
-
- out_df <- data.frame(
- concentration = conc,
- effect_tox_simple = effect_tox_simple_fn * effect_max,
- effect_tox = effect_tox_fn * effect_max,
- effect_tox_sys = effect_tox_sys_fn * effect_max,
- stress_tox = stress_tox_fn,
- sys_stress_tox = sys_stress_tox_fn,
- stress_tox_sys = stress_tox_sys_fn
- )
- if (with_env) {
- effect_tox_env_simple_fn = predict(
- effect_tox_env_mod_simple,
- data.frame(concentration = conc)
- )
- stress_tox_env_fn <- stress_tox_fn + stress_env
- effect_tox_env_fn <- stress_to_effect(
- stress_tox_env_fn, p, q
- )
- sys_stress_tox_env_fn <- predict(
- sys_stress_tox_env_mod,
- data.frame(stress_tox = stress_tox_fn)
- )
- stress_tox_env_sys_fn <- stress_tox_env_fn +
- sys_stress_tox_env_fn
- effect_tox_env_sys_fn <- stress_to_effect(
- stress_tox_env_sys_fn, p, q
- )
- out_df <- cbind(out_df, data.frame(
- effect_tox_env_simple = effect_tox_env_simple_fn * effect_max,
- effect_tox_env = effect_tox_env_fn * effect_max,
- effect_tox_env_sys = effect_tox_env_sys_fn * effect_max,
- stress_env = stress_env,
- stress_tox_env = stress_tox_env_fn,
- sys_stress_tox_env = sys_stress_tox_env_fn,
- stress_tox_env_sys = stress_tox_env_sys_fn
- ))
- }
- out_df
- }
-
- output$fn <- fn
-
-
# smooth curves -------------------------------------------------------
# In order to generate a broken x-axis the concentration vector must
# also be broken in two. The left part of the axis is supposed to be at
# 0 but because it's a log axis I have to make the values just really
# small. The concentrations in the gap won't be used for plotting later.
n_smooth <- 1000 # number of points to approximate the curves
+ conc_adjust_factor <- 10^-5
+ output$conc_adjust_factor <- conc_adjust_factor
concentration_smooth <- 10 ^ seq(
log10(min_conc * conc_adjust_factor),
log10(max(concentration)),
length.out = n_smooth
)
- output$curves <- fn(concentration_smooth)
+ output$curves <- predict_ecxsys(output, concentration_smooth)
output$curves$use_for_plotting <-
concentration_smooth < min_conc * conc_adjust_factor * 1.5 |
concentration_smooth > min_conc * 1.5
@@ -513,3 +284,99 @@ reset_options <- function(original_options) {
}
options(changed)
}
+
+
+fit_sys <- function(stress_external_observed,
+ stress_external_modeled,
+ stress_tox,
+ hormesis_index,
+ original_options) {
+ sys <- stress_external_observed - stress_external_modeled
+ sys <- clamp(sys)
+ sys[hormesis_index:length(sys)] <- 0
+ # Add sys to the output before it is fitted:
+ out <- list(sys = sys)
+ sys_mod <- tryCatch(
+ {
+ # There is no other way to suppress that one error message
+ # from optim except by changing the options temporarily.
+ warn_error_original <- original_options[c("warn", "show.error.messages")]
+ options(show.error.messages = FALSE)
+ drc::drm(sys ~ stress_tox, fct = drc::W1.3())
+ },
+ error = function(e) {
+ # Failure to converge often happens when all or almost all sys
+ # stress values are zero. Fitting a linear model in this case seems
+ # to be the most appropriate remedy.
+ stress_tox <- range(stress_tox)
+ sys <- c(0, 0)
+ lm(sys ~ stress_tox)
+ },
+ finally = options(warn_error_original)
+ )
+ out$sys_mod <- sys_mod
+ out$sys_modeled <- unname(
+ predict(sys_mod, data.frame(stress_tox = stress_tox))
+ )
+ out
+}
+
+
+moving_weighted_mean <- function(x) {
+ # This is used to smooth out points which are jumping up and down.
+ count <- rep(1, length(x))
+ x_diff <- diff(x) * -1
+ while (any(x_diff < 0, na.rm = TRUE)) {
+ i <- which(x_diff < 0)[1]
+ j <- i + 1
+ x[i] <- weighted.mean(x[i:j], count[i:j])
+ x <- x[-j]
+ count[i] <- count[i] + count[j]
+ count <- count[-j]
+ x_diff <- diff(x) * -1
+ }
+ rep(x, count)
+}
+
+
+fit_LL5_model <- function(min_conc,
+ concentration,
+ effect_observed,
+ original_options) {
+ # The traditional log-logistic model, here using the five-parameter
+ # log-logistic function drc::LL.5.
+ conc_with_control_shifted <- c(min_conc, concentration[-1])
+ effect_observed_averaged <- moving_weighted_mean(effect_observed)
+ interpolated <- approx(
+ log10(conc_with_control_shifted),
+ effect_observed_averaged,
+ n = 10
+ )
+ conc_interpolated <- 10^interpolated$x
+ effect_interpolated <- interpolated$y
+ effect_LL5_mod <- drc::drm(
+ effect_interpolated ~ conc_interpolated,
+ fct = drc::LL.5(fixed = c(NA, 0, effect_observed_averaged[1], NA, NA))
+ )
+ options(original_options["warn"])
+ effect_LL5 <- predict(
+ effect_LL5_mod,
+ data.frame(conc_interpolated = concentration)
+ )
+ list(
+ effect_LL5_mod = effect_LL5_mod,
+ effect_LL5 = effect_LL5
+ )
+}
+
+
+interpolate <- function(x, to_index, n_new, conc = FALSE) {
+ from_index <- to_index - 1 # subhormesis_index
+ len <- n_new + 2 # Add 2 because seq() includes the left and right end.
+ if (conc) {
+ x_new <- 10^seq(log10(x[from_index]), log10(x[to_index]), length.out = len)
+ } else {
+ x_new <- seq(x[from_index], x[to_index], length.out = len)
+ }
+ append(x, x_new[-c(1, len)], from_index)
+}
diff --git a/R/helpers.R b/R/helpers.R
new file mode 100644
index 0000000..f13f0d0
--- /dev/null
+++ b/R/helpers.R
@@ -0,0 +1,59 @@
+# A collection of internal helper functions which get used in more than one place.
+
+
+clamp <- function(x, lower = 0, upper = 1) {
+ # Returns lower if x < lower, returns upper if x > upper, else returns x
+ pmin(pmax(x, lower), upper)
+}
+
+
+get_log_ticks <- function(x) {
+ # Calculate the positions of major and minor tick marks on a base 10
+ # logarithmic axis.
+ stopifnot(min(x, na.rm = TRUE) > 0)
+ x <- log10(x)
+ major <- 10 ^ seq(
+ floor(min(x, na.rm = TRUE)),
+ ceiling(max(x, na.rm = TRUE))
+ )
+ n_between <- length(major) - 1
+ minor <- integer(n_between * 8)
+ for (i in 1:n_between) {
+ a <- major[i]
+ b <- major[i + 1]
+ minor[seq(i * 8 - 7, i * 8)] <- seq(a + a, b - a, a)
+ }
+ major_tick_labels <- formatC(major, format = "fg")
+ major_tick_labels[1] <- "0"
+ list(
+ major = major,
+ minor = minor,
+ major_labels = major_tick_labels
+ )
+}
+
+
+adjust_smooth_concentrations <- function(curves, conc_adjust_factor) {
+ # 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.
+
+ gap_idx <- min(which(!curves$use_for_plotting))
+
+ # Keep only the values to the left and right of the axis break:
+ curves <- curves[curves$use_for_plotting, ]
+
+ # Add NAs to force breaks in the lines:
+ axis_break_conc <- curves$concentration[gap_idx]
+ curves[gap_idx, ] <- NA
+
+ # 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
+
+ list(
+ curves = curves,
+ axis_break_conc = axis_break_conc
+ )
+}
diff --git a/R/moving_weighted_mean.R b/R/moving_weighted_mean.R
deleted file mode 100644
index a9e1e8e..0000000
--- a/R/moving_weighted_mean.R
+++ /dev/null
@@ -1,16 +0,0 @@
-moving_weighted_mean <- function(x) {
- # Helper function to calculate the moving weighted mean. This is used to
- # smooth out points which are jumping up and down.
- count <- rep(1, length(x))
- x_diff <- diff(x) * -1
- while (any(x_diff < 0, na.rm = TRUE)) {
- i <- which(x_diff < 0)[1]
- j <- i + 1
- x[i] <- weighted.mean(x[i:j], count[i:j])
- x <- x[-j]
- count[i] <- count[i] + count[j]
- count <- count[-j]
- x_diff <- diff(x) * -1
- }
- rep(x, count)
-}
diff --git a/R/plot_ecxsys.R b/R/plot_ecxsys.R
index cf3f28e..0ea5e49 100644
--- a/R/plot_ecxsys.R
+++ b/R/plot_ecxsys.R
@@ -8,74 +8,19 @@
#'
#' @name plot_ecxsys
#'
-#' @param model The list returned from ecxsys().
-#' @param show_simple_model Should the log-logistic models be plotted for
+#' @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.
#'
#' @examples model <- ecxsys(
-#' effect_tox_observed = c(85, 76, 94, 35, 0),
-#' effect_tox_env_observed = c(24, 23, 32, 0, 0),
#' concentration = c(0, 0.03, 0.3, 3, 10),
-#' hormesis_concentration = 0.3
+#' hormesis_concentration = 0.3,
+#' effect_tox_observed = c(85, 76, 94, 35, 0),
+#' effect_tox_env_observed = c(24, 23, 32, 0, 0)
#' )
#' plot_effect(model)
-#' plot_system_stress(model)
+#' plot_stress(model)
NULL
-
-
-
-# internal helper functions for plotting ----------------------------------
-
-get_log_ticks <- function(x) {
- # Calculate the positions of major and minor tick marks on a base 10
- # logarithmic axis.
- stopifnot(min(x, na.rm = TRUE) > 0)
- x <- log10(x)
- major <- 10 ^ seq(
- floor(min(x, na.rm = TRUE)),
- ceiling(max(x, na.rm = TRUE))
- )
- n_between <- length(major) - 1
- minor <- integer(n_between * 8)
- for (i in 1:n_between) {
- a <- major[i]
- b <- major[i+1]
- minor[seq(i * 8 - 7, i * 8)] <- seq(a + a, b - a, a)
- }
- major_tick_labels <- formatC(major, format = "fg")
- major_tick_labels[1] <- "0"
- list(
- major = major,
- minor = minor,
- major_labels = major_tick_labels
- )
-}
-
-
-adjust_smooth_concentrations <- function(curves, conc_adjust_factor) {
- # 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.
-
- gap_idx <- min(which(!curves$use_for_plotting))
-
- # Keep only the values to the left and right of the axis break:
- curves <- curves[curves$use_for_plotting, ]
-
- # Add NAs to force breaks in the lines:
- axis_break_conc <- curves$concentration[gap_idx]
- curves[gap_idx, ] <- NA
-
- # 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
-
- list(
- curves = curves,
- axis_break_conc = axis_break_conc
- )
-}
diff --git a/R/plot_effect.R b/R/plot_effect.R
index cc4995c..6b944e3 100644
--- a/R/plot_effect.R
+++ b/R/plot_effect.R
@@ -1,17 +1,18 @@
#' @rdname plot_ecxsys
#' @export
plot_effect <- function(model,
- show_simple_model = FALSE,
+ show_LL5_model = FALSE,
show_legend = FALSE,
xlab = "concentration",
ylab = "effect") {
+ stopifnot(inherits(model, "ecxsys"))
temp <- adjust_smooth_concentrations(
model$curves,
model$conc_adjust_factor
)
curves <- temp$curves
log_ticks <- get_log_ticks(curves$concentration)
- model$args$concentration[1] <- curves$concentration[1]
+ concentration <- c(curves$concentration[1], model$args$concentration[-1])
plot(
NA,
@@ -25,85 +26,82 @@ plot_effect <- function(model,
bty = "L",
las = 1
)
- if (show_simple_model) {
- lines(
- curves$concentration,
- curves$effect_tox_simple,
- col = rgb(112, 112, 207, maxColorValue = 255),
- lty = 4
- )
- if (model$with_env) {
- lines(
- curves$concentration,
- curves$effect_tox_env_simple,
- col = rgb(207, 112, 112, maxColorValue = 255),
- lty = 4
- )
- }
- }
lines(
curves$concentration,
curves$effect_tox_sys,
col = "blue"
)
- if (model$with_env) {
- lines(
- curves$concentration,
- curves$effect_tox_env_sys,
- col = "red"
- )
- }
lines(
curves$concentration,
curves$effect_tox,
col = "deepskyblue",
lty = 2
)
+ points(
+ concentration,
+ model$args$effect_tox_observed,
+ pch = 16,
+ col = "blue"
+ )
+
if (model$with_env) {
+ lines(
+ curves$concentration,
+ curves$effect_tox_env_sys,
+ col = "red"
+ )
lines(
curves$concentration,
curves$effect_tox_env,
col = "orange",
lty = 2
)
- }
- points(
- model$args$concentration,
- model$args$effect_tox_observed,
- pch = 16,
- col = "blue"
- )
- if (model$with_env) {
- # TODO: merge this with the same condition above.
points(
- model$args$concentration,
+ concentration,
model$args$effect_tox_env_observed,
pch = 16,
col = "red"
)
}
+
+ if (show_LL5_model) {
+ lines(
+ curves$concentration,
+ curves$effect_tox_LL5,
+ col = "darkblue",
+ lty = 4
+ )
+ if (model$with_env) {
+ lines(
+ curves$concentration,
+ curves$effect_tox_env_LL5,
+ col = "darkred",
+ lty = 4
+ )
+ }
+ }
+
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 = temp$axis_break_conc)
+
if (show_legend) {
- # TODO: Maybe do proper subscript in legend
- # TODO: Add simple models to the legend
- legend_text <- c(
- "effect (tox+sys)",
- "effect (tox+env+sys)",
- "toxicant effect",
- "tox_env effect"
- )
- legend_pch <- c(16, 16, NA, NA)
- legend_lty <- c(1, 1, 2, 2)
- legend_col <- c("blue", "red", "deepskyblue", "orange")
- if (!model$with_env) {
- keep <- c(1, 3)
- legend_text <- legend_text[keep]
- legend_pch <- legend_pch[keep]
- legend_lty <- legend_lty[keep]
- legend_col <- legend_col[keep]
+ 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(
"topright",
diff --git a/R/plot_system_stress.R b/R/plot_stress.R
similarity index 62%
rename from R/plot_system_stress.R
rename to R/plot_stress.R
index 72c7326..cccdcbc 100644
--- a/R/plot_system_stress.R
+++ b/R/plot_stress.R
@@ -1,6 +1,7 @@
#' @rdname plot_ecxsys
#' @export
-plot_system_stress <- function(model, show_legend = FALSE) {
+plot_stress <- function(model, show_legend = FALSE) {
+ stopifnot(inherits(model, "ecxsys"))
temp <- adjust_smooth_concentrations(
model$curves,
model$conc_adjust_factor
@@ -8,8 +9,7 @@ plot_system_stress <- function(model, show_legend = FALSE) {
curves <- temp$curves
axis_break_conc <- temp$axis_break_conc
log_ticks <- get_log_ticks(curves$concentration)
- # FIXME: Don't overwrite this arg. Make a new variable instead.
- model$args$concentration[1] <- curves$concentration[1]
+ concentration <- c(curves$concentration[1], model$args$concentration[-1])
plot(
NA,
@@ -23,43 +23,43 @@ plot_system_stress <- function(model, show_legend = FALSE) {
las = 1,
bty = "L"
)
+
lines(
curves$concentration,
- curves$sys_stress_tox,
+ curves$sys_tox,
col = "blue"
)
points(
- model$args$concentration,
- model$sys_stress_tox,
+ concentration,
+ model$sys_tox_not_fitted,
pch = 16,
col = "blue"
)
+
if (model$with_env) {
lines(
curves$concentration,
- curves$sys_stress_tox_env,
+ curves$sys_tox_env,
col = "red"
)
points(
- model$args$concentration,
- model$sys_stress_tox_env,
+ concentration,
+ model$sys_tox_env_not_fitted,
pch = 16,
col = "red"
)
}
+
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)
+
if (show_legend) {
- legend_text <- c(
- # TODO: Maybe do proper subscript in legend
- "system stress (tox)",
- "system stress (tox+env)"
- )
- legend_col <- c("blue", "red")
- if (!model$with_env) {
- legend_text <- legend_text[1]
- legend_col <- legend_col[1]
+ 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(
"topright",
diff --git a/R/predict_ecxsys.R b/R/predict_ecxsys.R
new file mode 100644
index 0000000..604dce4
--- /dev/null
+++ b/R/predict_ecxsys.R
@@ -0,0 +1,122 @@
+#' Predict ECxSyS at various concentrations
+#'
+#' @param model The output of a call to \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{effect_tox_sys}{Modeled effect resulting from toxicant and system
+#' stress.}
+#' \item{stress_tox}{The toxicant stress.}
+#' \item{sys_tox}{System stress under toxicant stress conditions
+#' 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}{Modeled effect resulting from toxicant and
+#' environmental stress.}
+#' \item{effect_tox_env_sys}{Modeled effect resulting from toxicant,
+#' environmental and system stress.}
+#' \item{stress_env}{Environmental stress.}
+#' \item{stress_tox_env}{The sum of toxicant and environmental stress.}
+#' \item{sys_tox_env}{System stress under toxicant and
+#' environmental stress conditions.}
+#' \item{stress_tox_env_sys}{The sum of \code{stress_tox_env} and
+#' \code{sys_tox_env}.}
+#' }
+#'
+#' @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)
+#' )
+#' p <- predict_ecxsys(model, c(0.001, 0.01, 0.1, 1, 10))
+#'
+#' @export
+predict_ecxsys <- function(model, concentration) {
+ # This function returns all modeled values at the provided
+ # concentrations.
+ stopifnot(
+ inherits(model, "ecxsys"),
+ is.numeric(concentration)
+ )
+ p <- model$args$p
+ q <- model$args$q
+ effect_max <- model$args$effect_max
+ out_df <- data.frame(
+ concentration = concentration
+ )
+
+ out_df$effect_tox_LL5 <- predict(
+ model$effect_tox_LL5_mod,
+ data.frame(concentration = concentration)
+ ) * effect_max
+
+ if (model$with_env) {
+ out_df$effect_tox_env_LL5 <- predict(
+ model$effect_tox_env_LL5_mod,
+ data.frame(concentration = concentration)
+ ) * effect_max
+ }
+
+ effect_tox <- predict(
+ model$effect_tox_mod,
+ data.frame(concentration = concentration)
+ )
+ out_df$effect_tox <- effect_tox * effect_max
+
+ stress_tox <- effect_to_stress(effect_tox, p, q)
+ out_df$stress_tox <- stress_tox
+
+ sys_tox <- predict(
+ model$sys_tox_mod,
+ data.frame(stress_tox = stress_tox)
+ )
+ out_df$sys_tox <- sys_tox
+
+ stress_tox_sys <- stress_tox + sys_tox
+ out_df$stress_tox_sys <- stress_tox_sys
+
+ effect_tox_sys <- stress_to_effect(stress_tox_sys, p, q)
+ out_df$effect_tox_sys <- effect_tox_sys * effect_max
+
+ if (model$with_env) {
+ out_df$stress_env <- model$stress_env
+
+ stress_tox_env <- stress_tox + model$stress_env
+ out_df$stress_tox_env <- stress_tox_env
+
+ out_df$effect_tox_env <- stress_to_effect(
+ stress_tox_env, p, q
+ ) * effect_max
+
+ sys_tox_env <- predict(
+ model$sys_tox_env_mod,
+ data.frame(stress_tox = stress_tox)
+ )
+ out_df$sys_tox_env <- sys_tox_env
+
+ stress_tox_env_sys <- stress_tox_env +
+ sys_tox_env
+ out_df$stress_tox_env_sys <- stress_tox_env_sys
+
+ effect_tox_env_sys <- stress_to_effect(
+ stress_tox_env_sys, p, q
+ )
+ out_df$effect_tox_env_sys <- effect_tox_env_sys * effect_max
+ }
+
+ class(out_df) <- c("ecxsys_predicted", class(out_df))
+ out_df
+}
diff --git a/R/stress_effect_conversion.R b/R/stress_effect_conversion.R
index 512f86a..ef6c252 100644
--- a/R/stress_effect_conversion.R
+++ b/R/stress_effect_conversion.R
@@ -6,7 +6,7 @@
#' These are simple wrappers around the beta distribution function
#' \code{\link[stats:Beta]{pbeta}} and the beta quantile function
#' \code{\link[stats:Beta]{qbeta}}. \code{stress_to_effect} returns
-#' \code{1 - pbeta(stress, p, q)}. \code{effect_to_stress} first clips the
+#' \code{1 - pbeta(stress, p, q)}. \code{effect_to_stress} first clamps the
#' effect to the interval [0, 1] and then returns
#' \code{qbeta(1 - effect, p, q)}.
#'
@@ -15,11 +15,12 @@
#'
#' @name Stressconversion
#'
-#' @param effect One or more effect values to convert to general stress.
-#' Should be a value between 0 and 1. Smaller or bigger values are treated as
-#' 0 or 1 respectively.
+#' @param effect One or more effect values to convert to general stress. Should
+#' be a value between 0 and 1. Smaller or bigger values are treated as 0 or 1
+#' respectively.
#' @param stress One or more stress values to convert to effect.
-#' @param p,q The shape parameters of the beta distribution. Default is 3.2.
+#' @param p,q The shape parameters of the \code{\link[stats:Beta]{beta}}
+#' distribution. Default is 3.2.
#'
#' @examples
#' stress <- 0.3
@@ -32,7 +33,8 @@ NULL
#' @rdname Stressconversion
#' @export
effect_to_stress <- function(effect, p = 3.2, q = 3.2) {
- effect <- pmin(pmax(effect, 0), 1)
+ stopifnot(p >= 0, q >= 0)
+ effect <- clamp(effect)
qbeta(1 - effect, p, q)
}
@@ -40,5 +42,6 @@ effect_to_stress <- function(effect, p = 3.2, q = 3.2) {
#' @rdname Stressconversion
#' @export
stress_to_effect <- function(stress, p = 3.2, q = 3.2) {
+ stopifnot(p >= 0, q >= 0)
1 - pbeta(stress, p, q)
}
diff --git a/R/stressaddition-package.R b/R/stressaddition-package.R
index 649790b..16dd5a0 100644
--- a/R/stressaddition-package.R
+++ b/R/stressaddition-package.R
@@ -2,6 +2,10 @@
# ?`stressaddition-package`. Some people recommend using "@keywords internal"
# here to exclude this page from the help index but I'm not going to do that.
+# "_PACKAGE" automatically populates most sections in this package
+# documentation, so I don't need to manually fill in authors, title or
+# description.
+
#' @details See the publication linked below for more information including
#' equations.
@@ -23,8 +27,3 @@
#' @import graphics
#' @import grDevices
"_PACKAGE"
-
-
-# "_PACKAGE" automatically populates most sections in this package
-# documentation, so I don't need to manually fill in authors, title or
-# description.
diff --git a/README.md b/README.md
index 7bc2506..d52dcf3 100644
--- a/README.md
+++ b/README.md
@@ -14,8 +14,7 @@ Alternatively there are binary and source builds downloadable from the [releases
## Updating
RStudio's integrated package updater won't detect updates in packages installed from GitHub or GitLab. I recommend running
```r
-# install.packages("remotes")
-remotes::update_packages()
+devtools::update_packages()
```
in regular intervals to check for updates from those sources.
@@ -28,15 +27,15 @@ In the paper we use the model in the context of survival experiments. However, i
library(stressaddition)
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),
- hormesis_concentration = 0.3
+ effect_tox_env_observed = c(24, 23, 32, 0, 0)
)
# Plot the effect and the system stress:
par(mfrow = c(2, 1))
plot_effect(model)
-plot_system_stress(model)
+plot_stress(model)
# The LC50 of the effect under the influence of toxicant and system tress:
ec(model, "effect_tox_sys", 50)
diff --git a/man/Stressconversion.Rd b/man/Stressconversion.Rd
index b84b0f9..2762fee 100644
--- a/man/Stressconversion.Rd
+++ b/man/Stressconversion.Rd
@@ -11,11 +11,12 @@ effect_to_stress(effect, p = 3.2, q = 3.2)
stress_to_effect(stress, p = 3.2, q = 3.2)
}
\arguments{
-\item{effect}{One or more effect values to convert to general stress.
-Should be a value between 0 and 1. Smaller or bigger values are treated as
-0 or 1 respectively.}
+\item{effect}{One or more effect values to convert to general stress. Should
+be a value between 0 and 1. Smaller or bigger values are treated as 0 or 1
+respectively.}
-\item{p, q}{The shape parameters of the beta distribution. Default is 3.2.}
+\item{p, q}{The shape parameters of the \code{\link[stats:Beta]{beta}}
+distribution. Default is 3.2.}
\item{stress}{One or more stress values to convert to effect.}
}
@@ -27,7 +28,7 @@ distribution.
These are simple wrappers around the beta distribution function
\code{\link[stats:Beta]{pbeta}} and the beta quantile function
\code{\link[stats:Beta]{qbeta}}. \code{stress_to_effect} returns
-\code{1 - pbeta(stress, p, q)}. \code{effect_to_stress} first clips the
+\code{1 - pbeta(stress, p, q)}. \code{effect_to_stress} first clamps the
effect to the interval [0, 1] and then returns
\code{qbeta(1 - effect, p, q)}.
diff --git a/man/ec.Rd b/man/ec.Rd
index f5dd80b..1a9c828 100644
--- a/man/ec.Rd
+++ b/man/ec.Rd
@@ -4,40 +4,59 @@
\alias{ec}
\title{Effect Concentrations}
\usage{
-ec(model, effect_name, target_effect)
+ec(model, response_name, response_level)
}
\arguments{
-\item{model}{The object returned from \code{ecxsys()}.}
+\item{model}{This can be one of three types of objects: Either the output of
+\code{\link{ecxsys}} or the output of \code{\link{predict_ecxsys}} or a
+data frame with a "concentration" column and a \code{response_name} column.
+See the examples.}
-\item{effect_name}{The name of the effect for which you want to calculate the
-EC. Must be the name of a column in \code{model$curves}.}
+\item{response_name}{The name of the effect or stress for which you want to
+calculate the EC. Must be one of \code{colnames(model$curves)}.}
-\item{target_effect}{The desired effect percentage between 0 and 100. For
-example with the value 10 the function will return the EC10, i.e. the
-concentration where the response falls below 90 \% of the maximum possible
+\item{response_level}{The desired response level as a percentage between 0
+and 100. For example with the value 10 the function will return the EC10,
+i.e. the concentration where the response falls below 90 \% of the maximum
response.}
}
\value{
-A list containing the effect concentration and the corresponding
- effect.
+A list containing the response concentration and the corresponding
+ response value.
}
\description{
-Estimate the concentration to reach a certain effect relative to the control.
+Estimate the concentration to reach a certain effect or stress level relative
+to the control.
}
\details{
If the response level occurs multiple times because of hormesis, which may
-happen for low values of \code{target_effect}, then the occurrence with the
+happen for low values of \code{response_level}, then the occurrence with the
smallest concentration is returned.
}
\examples{
+# Calculate the EC_10, the concentration where the effect falls
+# 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),
- effect_tox_env_observed = c(24, 23, 32, 0, 0),
hormesis_concentration = 0.3
)
-# Calculate the EC_10, the concentration where the effect falls
-# below 90 \% of the effect in the control:
-ec_10 <- ec(model, "effect_tox_sys", 10)
+
+# using the ecxsys() output or the curves therein directly:
+ec(model, "effect_tox_sys", 10)
+ec(model$curves, "effect_tox_sys", 10)
+
+# using the output of predict_ecxsys() with custom concentrations:
+conc <- 10^seq(-9, 1, length.out = 1000)
+curves <- predict_ecxsys(model, conc)
+ec(curves, "effect_tox_sys", 10)
+
+# using a custom data frame:
+df_custom <- data.frame(
+ concentration = curves$concentration,
+ foo = curves$effect_tox_sys
+)
+ec(df_custom, "foo", 10)
}
diff --git a/man/ecxsys.Rd b/man/ecxsys.Rd
index e9ce10d..ca2cabc 100644
--- a/man/ecxsys.Rd
+++ b/man/ecxsys.Rd
@@ -6,10 +6,9 @@
\usage{
ecxsys(
concentration,
- effect_tox_observed,
- effect_tox_env_observed,
hormesis_concentration,
- hormesis_index,
+ effect_tox_observed,
+ effect_tox_env_observed = NULL,
effect_max = 100,
p = 3.2,
q = 3.2
@@ -20,22 +19,15 @@ ecxsys(
indicate the control. Should be sorted in ascending order, otherwise it
will be sorted automatically.}
+\item{hormesis_concentration}{The concentration where the hormesis occurs.
+This is usually the concentration of the highest effect after the control.}
+
\item{effect_tox_observed}{A vector of effect values observed at the given
concentrations and in absence of environmental stress. Values must be
between 0 and \code{effect_max}.}
\item{effect_tox_env_observed}{Effect values observed in the presence of
-environmental stress. This argument is optional and can be left out to
-model without environmental stress. Values must be between 0 and
-\code{effect_max}.}
-
-\item{hormesis_concentration}{The concentration where the hormesis occurs.
-This is usually the concentration of the highest effect after the control.}
-
-\item{hormesis_index}{\strong{Deprecated, will be removed soon.} A single
-integer specifying the index of the hormesis concentration in the
-concentration vector. This argument exists for compatibility with older
-versions of this function.}
+environmental stress. Must be between 0 and \code{effect_max}.}
\item{effect_max}{The maximum value the effect could possibly reach. For
survival data in percent this should be 100 (the default).}
@@ -43,43 +35,10 @@ 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{
-The result is a list containing many different objects with the most
- important being \code{curves} and \code{fn}. You can use \code{fn()} to
- calculate the curves at your concentrations of choice, see examples.
- \code{curves} is a data frame 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_simple}{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{effect_tox_sys}{Modeled effect resulting from toxicant and system
- stress.}
- \item{stress_tox}{The toxicant stress.}
- \item{sys_stress_tox}{System stress under toxicant stress conditions
- without environmental stress.}
- \item{stress_tox_sys}{The sum of \code{stress_tox} and
- \code{sys_stress_tox}.}
- \item{effect_tox_env_simple}{The five-parameter log-logistic model of the
- effect 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,
- environmental and system stress.}
- \item{stress_env}{Environmental stress.}
- \item{stress_tox_env}{The sum of toxicant and environmental stress.}
- \item{sys_stress_tox_env}{System stress under toxicant and
- environmental stress conditions.}
- \item{stress_tox_env_sys}{The sum of \code{stress_tox_env} and
- \code{sys_stress_tox_env}.}
- \item{use_for_plotting}{A boolean vector which is used in the plotting
- functions. It controls which parts of the curves are removed for the
- broken concentration axis.}
- }
+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.
}
\description{
The ECx-SyS model for modeling concentration-effect relationships which
@@ -92,29 +51,25 @@ at or close to zero. If the model does not fit properly try adding an effect
of 0 at ten times the maximum observed concentration.
The vectors \code{concentration}, \code{effect_tox_observed} and
-\code{effect_tox_env_observed} must be of equal length and should be sorted
-by increasing concentration.
+\code{effect_tox_env_observed} (if provided) must be of equal length and
+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),
- hormesis_concentration = 0.3
+ effect_tox_env_observed = c(24, 23, 32, 0, 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:
model <- ecxsys(
concentration = c(0, 0.03, 0.3, 3, 10),
+ 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),
- hormesis_concentration = 0.3,
effect_max = 20
)
-# The returned object contains a function which is useful for calculating
-# effect and stress values at specific concentrations:
-model$fn(c(0, 0.01, 0.1, 1))
-
}
diff --git a/man/plot_ecxsys.Rd b/man/plot_ecxsys.Rd
index 6af7b2e..b5d862e 100644
--- a/man/plot_ecxsys.Rd
+++ b/man/plot_ecxsys.Rd
@@ -1,26 +1,25 @@
% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/plot_ecxsys.R, R/plot_effect.R,
-% R/plot_system_stress.R
+% Please edit documentation in R/plot_ecxsys.R, R/plot_effect.R, R/plot_stress.R
\name{plot_ecxsys}
\alias{plot_ecxsys}
\alias{plot_effect}
-\alias{plot_system_stress}
+\alias{plot_stress}
\title{Plot the results of the ECx-SyS model}
\usage{
plot_effect(
model,
- show_simple_model = FALSE,
+ show_LL5_model = FALSE,
show_legend = FALSE,
xlab = "concentration",
ylab = "effect"
)
-plot_system_stress(model, show_legend = FALSE)
+plot_stress(model, show_legend = FALSE)
}
\arguments{
-\item{model}{The list returned from ecxsys().}
+\item{model}{The list returned from \code{\link{ecxsys}}.}
-\item{show_simple_model}{Should the log-logistic models be plotted for
+\item{show_LL5_model}{Should the log-logistic models be plotted for
comparison? Defaults to \code{FALSE}.}
\item{show_legend}{Should the plot include a legend? Defaults to FALSE
@@ -34,11 +33,11 @@ system stress with and without environmental stress.
}
\examples{
model <- ecxsys(
- effect_tox_observed = c(85, 76, 94, 35, 0),
- effect_tox_env_observed = c(24, 23, 32, 0, 0),
concentration = c(0, 0.03, 0.3, 3, 10),
- hormesis_concentration = 0.3
+ hormesis_concentration = 0.3,
+ effect_tox_observed = c(85, 76, 94, 35, 0),
+ effect_tox_env_observed = c(24, 23, 32, 0, 0)
)
plot_effect(model)
-plot_system_stress(model)
+plot_stress(model)
}
diff --git a/man/predict_ecxsys.Rd b/man/predict_ecxsys.Rd
new file mode 100644
index 0000000..194496d
--- /dev/null
+++ b/man/predict_ecxsys.Rd
@@ -0,0 +1,59 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/predict_ecxsys.R
+\name{predict_ecxsys}
+\alias{predict_ecxsys}
+\title{Predict ECxSyS at various concentrations}
+\usage{
+predict_ecxsys(model, concentration)
+}
+\arguments{
+\item{model}{The output of a call to \code{\link{ecxsys}}.}
+
+\item{concentration}{A numeric vector of concentrations.}
+}
+\value{
+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{effect_tox_sys}{Modeled effect resulting from toxicant and system
+ stress.}
+ \item{stress_tox}{The toxicant stress.}
+ \item{sys_tox}{System stress under toxicant stress conditions
+ 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}{Modeled effect resulting from toxicant and
+ environmental stress.}
+ \item{effect_tox_env_sys}{Modeled effect resulting from toxicant,
+ environmental and system stress.}
+ \item{stress_env}{Environmental stress.}
+ \item{stress_tox_env}{The sum of toxicant and environmental stress.}
+ \item{sys_tox_env}{System stress under toxicant and
+ environmental stress conditions.}
+ \item{stress_tox_env_sys}{The sum of \code{stress_tox_env} and
+ \code{sys_tox_env}.}
+ }
+}
+\description{
+Predict ECxSyS at various 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)
+)
+p <- predict_ecxsys(model, c(0.001, 0.01, 0.1, 1, 10))
+
+}
diff --git a/tests/testthat/test-ec.R b/tests/testthat/test-ec.R
index dc1e3ae..f5deb23 100644
--- a/tests/testthat/test-ec.R
+++ b/tests/testthat/test-ec.R
@@ -1,43 +1,80 @@
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)
+ )
+
+ # 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)
+
+ # using the output of predict_ecxsys() with custom concentrations:
+ conc <- 10^seq(-9, 1, length.out = 1000)
+ curves <- predict_ecxsys(model, conc)
+ ec10_c <- ec(curves, "effect_tox_sys", 10)
+
+ # using a custom data frame:
+ df_custom <- data.frame(
+ concentration = curves$concentration,
+ foo = curves$effect_tox_sys
+ )
+ 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)
+})
+
+
test_that("ec values have not changed", {
- result <- ecxsys(
+ 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),
- hormesis_index = 3
+ effect_tox_env_observed = c(24, 23, 32, 0, 0)
)
expect_equal(
- ec(result, "effect_tox_sys",50),
- list(effect = 42.38984413, concentration = 2.547711893)
+ ec(model, "effect_tox_sys", 50),
+ list(response_value = 42.389844, concentration = 2.547712),
+ tolerance = 1e-4
)
expect_equal(
- ec(result, "effect_tox_sys", 10),
- list(effect = 76.30171944, concentration = 1.015730561)
+ ec(model, "effect_tox_sys", 10),
+ list(response_value = 76.301719, concentration = 1.015731),
+ tolerance = 1e-4
)
expect_equal(
- ec(result, "effect_tox_sys", 5),
- list(effect = 80.54070385, concentration = 0.003172831366)
+ ec(model, "effect_tox_sys", 5),
+ list(response_value = 80.540705, concentration = 0.003173),
+ tolerance = 1e-4
)
expect_equal(
- ec(result, "effect_tox_sys", 1),
- list(effect = 83.93189138, concentration = 0.00005727730324)
+ ec(model, "effect_tox_sys", 1),
+ list(response_value = 83.931891, concentration = 0.000057),
+ tolerance = 1e-4
)
expect_equal(
- ec(result, "effect_tox_env_sys", 50),
- list(effect = 12.1555922, concentration = 0.9058480238)
+ ec(model, "effect_tox_env_sys", 50),
+ list(response_value = 12.155592, concentration = 0.905848),
+ tolerance = 1e-4
)
expect_equal(
- ec(result, "effect_tox_env_sys", 10),
- list(effect = 21.88006596, concentration = 0.0007604409081)
+ ec(model, "effect_tox_env_sys", 10),
+ list(response_value = 21.880066, concentration = 0.0007604),
+ tolerance = 1e-4
)
expect_equal(
- ec(result, "effect_tox_env_sys", 5),
- list(effect = 23.09562518, concentration = 0.0001004406979)
+ ec(model, "effect_tox_env_sys", 5),
+ list(response_value = 23.095625, concentration = 0.00010044),
+ tolerance = 1e-4
)
expect_equal(
- ec(result, "effect_tox_env_sys", 1),
- list(effect = 24.06807255, concentration = 0.000001640834495)
+ ec(model, "effect_tox_env_sys", 1),
+ list(response_value = 24.068073, concentration = 0.000002),
+ tolerance = 1e-4
)
})
diff --git a/tests/testthat/test-ecxsys.R b/tests/testthat/test-ecxsys.R
index 3195874..1469b47 100644
--- a/tests/testthat/test-ecxsys.R
+++ b/tests/testthat/test-ecxsys.R
@@ -3,65 +3,67 @@ context("ecxsys")
# Create model here for use in the various tests to save typing:
mod <- ecxsys(
- effect_tox_observed = c(85, 76, 94, 35, 0),
- effect_tox_env_observed = c(24, 23, 32, 0, 0),
concentration = c(0, 0.03, 0.3, 3, 10),
- hormesis_index = 3
- # hormesis_concentration = 0.3
+ hormesis_concentration = 0.3,
+ effect_tox_observed = c(85, 76, 94, 35, 0),
+ effect_tox_env_observed = c(24, 23, 32, 0, 0)
)
-test_that("error when both hormesis arguments passed to ecxsys", {
+test_that("error when hormesis_concentration not in concentration", {
+ errorm <- "hormesis_concentration must equal one of the concentration values."
expect_error(
ecxsys(
- effect_tox_observed = c(85, 76, 94, 35, 0),
concentration = c(0, 0.03, 0.3, 3, 10),
- hormesis_index = 3,
- hormesis_concentration = 0.3
- )
+ hormesis_concentration = 0.4,
+ effect_tox_observed = c(85, 76, 94, 35, 0)
+ ),
+ errorm,
+ fixed = TRUE
)
-})
-
-
-test_that("error when hormesis_concentration not in concentration", {
expect_error(
ecxsys(
- effect_tox_observed = c(85, 76, 94, 35, 0),
concentration = c(0, 0.03, 0.3, 3, 10),
- hormesis_concentration = 0.4
- )
+ hormesis_concentration = 30,
+ effect_tox_observed = c(85, 76, 94, 35, 0)
+ ),
+ errorm,
+ fixed = TRUE
)
})
test_that("error when hormesis_index <= 2 or >= (length(concentration))", {
- expect_error(
- ecxsys(
- effect_tox_observed = c(85, 76, 94, 35, 0),
- concentration = c(0, 0.03, 0.3, 3, 10),
- hormesis_index = 1
- )
+ errorm <- paste(
+ "hormesis_concentration must be greater than the lowest",
+ "non-control concentration and less than the highest concentration"
)
expect_error(
ecxsys(
- effect_tox_observed = c(85, 76, 94, 35, 0),
concentration = c(0, 0.03, 0.3, 3, 10),
- hormesis_index = 2
- )
+ hormesis_concentration = 0,
+ effect_tox_observed = c(85, 76, 94, 35, 0)
+ ),
+ errorm,
+ fixed = TRUE
)
expect_error(
ecxsys(
- effect_tox_observed = c(85, 76, 94, 35, 0),
concentration = c(0, 0.03, 0.3, 3, 10),
- hormesis_index = 5
- )
+ hormesis_concentration = 0.03,
+ effect_tox_observed = c(85, 76, 94, 35, 0)
+ ),
+ errorm,
+ fixed = TRUE
)
expect_error(
ecxsys(
- effect_tox_observed = c(85, 76, 94, 35, 0),
concentration = c(0, 0.03, 0.3, 3, 10),
- hormesis_index = 6
- )
+ hormesis_concentration = 10,
+ effect_tox_observed = c(85, 76, 94, 35, 0)
+ ),
+ errorm,
+ fixed = TRUE
)
})
@@ -73,112 +75,81 @@ test_that("min(concentration) == 0 is shifted the correct amount", {
test_that("the discrete results have not changed", {
expect_equal(
- round(mod$effect_tox_simple, 3),
- c(85.000, 84.542, 78.946, 31.959, 2.610)
+ mod$effect_tox_LL5,
+ c(85, 84.5420, 78.9464, 31.9586, 2.6096),
+ tolerance = 1e-4
)
expect_equal(
- round(mod$effect_tox, 3),
- c(100.000, 99.676, 94.320, 34.860, 0.840)
+ mod$effect_tox,
+ c(100, 99.6760, 94.3198, 34.8601, 0.8400),
+ tolerance = 1e-4
)
expect_equal(
- round(mod$stress_tox, 3),
- c(0.000, 0.078, 0.207, 0.579, 0.893)
+ mod$stress_tox,
+ c(0, 0.0784, 0.2067, 0.5794, 0.8927),
+ tolerance = 1e-4
)
expect_equal(
- round(mod$sys_stress_tox, 3),
- c(0.296, 0.280, 0.000, 0.000, 0.000)
+ mod$sys_tox_not_fitted,
+ c(0.2965, 0.2795, 0, 0, 0),
+ tolerance = 1e-4
)
expect_equal(
- round(mod$effect_tox_sys, 3),
- c(84.791, 77.481, 93.138, 34.86, 0.84)
+ mod$effect_tox_sys,
+ c(84.7912, 77.4811, 93.1380, 34.8602, 0.8400),
+ tolerance = 1e-4
)
- expect_equal(round(mod$stress_env, 3), 0.388)
+ expect_equal(mod$stress_env, 0.3885, tolerance = 1e-3)
expect_equal(
- round(mod$effect_tox_env_simple, 3),
- c(26.333, 26.333, 26.062, 0.710, 0.037)
+ mod$effect_tox_env_LL5,
+ c(26.3333, 26.3333, 26.0620, 0.7104, 0.0374),
+ tolerance = 1e-4
)
expect_equal(
- round(mod$stress_tox_env, 3),
- c(0.388, 0.467, 0.595, 0.968, 1.281)
+ mod$stress_tox_env,
+ c(0.3884, 0.4669, 0.5952, 0.9679, 1.2812),
+ tolerance = 1e-4
)
expect_equal(
- round(mod$effect_tox_env, 3),
- c(70.879, 56.414, 32.000, 0.020, 0.000)
+ mod$effect_tox_env,
+ c(70.8786, 56.4138, 32, 0.0202, 0),
+ tolerance = 1e-4
)
expect_equal(
- round(mod$sys_stress_tox_env, 3),
- c(0.254, 0.181, 0.000, 0.000, 0.000)
+ mod$sys_tox_env_not_fitted,
+ c(0.2537, 0.1815, 0, 0, 0),
+ tolerance = 1e-4
)
expect_equal(
- round(mod$effect_tox_env_sys, 3),
- c(24.325, 22.323, 29.386, 0.020, 0.000)
+ mod$effect_tox_env_sys,
+ c(24.3254, 22.3232, 29.3860, 0.0202, 0),
+ tolerance = 1e-4
)
})
-test_that("the smooth curves have not changed", {
- curves <- mod$curves
- i <- c(1, 714, 810, 905, 1000)
- expect_equal(
- round(curves$concentration[i], 3),
- c(0.000, 0.014, 0.125, 1.120, 10.000)
- )
- expect_equal(
- round(curves$effect_tox_simple[i], 3),
- c(85.000, 84.812, 82.700, 61.288, 2.610)
- )
- expect_equal(
- round(curves$effect_tox[i], 3),
- c(100.000, 99.879, 98.064, 73.665, 0.840)
- )
- expect_equal(
- round(curves$effect_tox_sys[i], 3),
- c(84.780, 77.932, 87.580, 73.665, 0.840)
- )
- expect_equal(
- round(curves$stress_tox[i], 3),
- c(0.000, 0.057, 0.142, 0.372, 0.893)
- )
- expect_equal(
- round(curves$sys_stress_tox[i], 3),
- c(0.298, 0.289, 0.134, 0.000, 0.000)
- )
- expect_equal(
- round(curves$stress_tox_sys[i], 3),
- c(0.298, 0.346, 0.276, 0.372, 0.893)
- )
- expect_equal(
- round(curves$effect_tox_env[i], 3),
- c(70.863, 60.496, 44.083, 8.514, 0)
- )
- expect_equal(
- round(curves$effect_tox_env_sys[i], 3),
- c(24.311, 20.727, 29.948, 8.513, 0.000)
- )
- expect_equal(
- round(curves$stress_env[i], 3),
- c(0.388, 0.388, 0.388, 0.388, 0.388)
- )
- expect_equal(
- round(curves$effect_tox_env_simple[i], 3),
- c(26.333, 26.333, 26.329, 7.538, 0.037)
- )
- expect_equal(
- round(curves$stress_tox_env[i], 3),
- c(0.389, 0.445, 0.531, 0.761, 1.281)
- )
- expect_equal(
- round(curves$sys_stress_tox_env[i], 3),
- c(0.252, 0.218, 0.076, 0.000, 0.000)
- )
- expect_equal(
- round(curves$stress_tox_env_sys[i], 3),
- c(0.640, 0.663, 0.607, 0.761, 1.281)
- )
- expect_equal(
- curves$use_for_plotting[i],
- c(TRUE, TRUE, TRUE, TRUE, TRUE)
- )
+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),
+ use_for_plotting = c(TRUE, TRUE, TRUE, TRUE, TRUE)
+ )
+ class(new_curves) <- class(reference_curves)
+ expect_equal(new_curves, reference_curves, tolerance = 1e-3)
})
@@ -186,36 +157,38 @@ test_that("the returned fn works the same way as internally", {
# I don't know why it would fail but it doesn't hurt to test it.
curves <- mod$curves
curves$use_for_plotting <- NULL # remove column "use_for_plotting"
- expect_identical(curves, mod$fn(curves$concentration))
+ expect_identical(curves, predict_ecxsys(mod, curves$concentration))
})
test_that("function arguments are returned unchanged", {
- expect_equal(mod$args$effect_tox_observed, c(85, 76, 94, 35, 0))
- expect_equal(mod$args$effect_tox_env_observed, c(24, 23, 32, 0, 0))
- expect_equal(mod$args$concentration, c(0, 0.03, 0.3, 3, 10))
- expect_equal(mod$args$hormesis_index, 3)
- # expect_equal(mod$args$hormesis_concentration, 0.3)
- expect_equal(mod$args$effect_max, 100)
- expect_equal(mod$args$p, 3.2)
- expect_equal(mod$args$q, 3.2)
+ 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),
+ effect_max = 100,
+ p = 3.2,
+ q = 3.2
+ )
+ expect_identical(args_reference, mod$args)
})
test_that("results are independent of concentration shift", {
mod_2 <- ecxsys(
- effect_tox_observed = c(85, 76, 94, 35, 0),
- effect_tox_env_observed = c(24, 23, 32, 0, 0),
concentration = c(0, 0.03, 0.3, 3, 10) * 2,
- hormesis_index = 3
+ hormesis_concentration = 0.3 * 2,
+ effect_tox_observed = c(85, 76, 94, 35, 0),
+ effect_tox_env_observed = c(24, 23, 32, 0, 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(
- effect_tox_observed = c(85, 76, 94, 35, 0),
- effect_tox_env_observed = c(24, 23, 32, 0, 0),
concentration = c(0, 0.03, 0.3, 3, 10) * 10,
- hormesis_index = 3
+ 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 shifts by factors other than powers of 10 will affect
# the result because of the way the zero concentration is "corrected".
@@ -230,9 +203,9 @@ test_that("results are independent of concentration shift", {
test_that("effect_tox_env_observed can be left out", {
mod_without_env <- ecxsys(
- effect_tox_observed = c(85, 76, 94, 35, 0),
concentration = c(0, 0.03, 0.3, 3, 10),
- hormesis_index = 3
+ hormesis_concentration = 0.3,
+ effect_tox_observed = c(85, 76, 94, 35, 0)
)
expect_equal(mod$effect_tox_sys, mod_without_env$effect_tox_sys)
expect_equal(mod$curves$effect_tox_sys,
@@ -240,12 +213,17 @@ test_that("effect_tox_env_observed can be left out", {
})
-test_that("model not converging produces warnings", {
+test_that("sys model not converging produces a warning, not an error", {
expect_warning(
ecxsys(
concentration = c(0, 0.1, 0.5, 1, 10, 33),
- effect_tox_observed = c(98, 98, 96, 76, 26, 0),
- hormesis_index = 3
- )
+ hormesis_concentration = 0.5,
+ effect_tox_observed = c(98, 98, 96, 76, 26, 0)
+ ),
+ paste(
+ "Using a horizontal linear model for sys_tox_mod because the",
+ "Weibull model did not converge."
+ ),
+ fixed = TRUE
)
})
diff --git a/tests/testthat/test-stress_effect_conversion.R b/tests/testthat/test-stress_effect_conversion.R
index abd1317..506491d 100644
--- a/tests/testthat/test-stress_effect_conversion.R
+++ b/tests/testthat/test-stress_effect_conversion.R
@@ -5,16 +5,16 @@ test_that("stress_to_effect handles bad input", {
expect_equal(stress_to_effect(-3), 1)
expect_equal(stress_to_effect(10), 0)
expect_error(stress_to_effect("a"))
- expect_warning(stress_to_effect(0.5, -3, 10))
- expect_warning(stress_to_effect(0.5, 3, -10))
+ expect_error(stress_to_effect(0.5, -3, 10))
+ expect_error(stress_to_effect(0.5, 3, -10))
})
test_that("effect_to_stress handles bad input", {
expect_equal(effect_to_stress(-15), 1)
expect_equal(effect_to_stress(20), 0)
expect_error(effect_to_stress("a"))
- expect_warning(effect_to_stress(0.5, -3, 10))
- expect_warning(effect_to_stress(0.5, 3, -10))
+ expect_error(effect_to_stress(0.5, -3, 10))
+ expect_error(effect_to_stress(0.5, 3, -10))
})
test_that("stress_to_effect is correct", {
--
GitLab