From fdca761396edbf01e32921e81fb243795ec0664e Mon Sep 17 00:00:00 2001 From: Sebastian Henz <sebastian.henz@ufz.de> Date: Thu, 7 May 2020 19:03:40 +0200 Subject: [PATCH] Enable check on the previous minor R release closes #38 --- .gitlab-ci.yml | 26 +- .../merge_into_master.md | 2 +- DESCRIPTION | 6 +- NAMESPACE | 11 +- NEWS.md | 15 ++ R/ecxsys.R | 233 ++++++++++-------- R/helpers.R | 30 +-- R/{ec.R => lc.R} | 71 +++--- R/log10_ticks.R | 72 ++++++ R/{predict_mixture.R => multi_tox.R} | 84 ++++--- R/plot_ecxsys.R | 26 +- R/plot_stress.R | 23 +- R/{plot_effect.R => plot_survival.R} | 70 +++--- R/predict_ecxsys.R | 56 ++--- ...version.R => stress_survival_conversion.R} | 34 ++- R/stressaddition-package.R | 5 +- README.md | 30 +-- man/Stressconversion.Rd | 36 ++- man/ecxsys.Rd | 74 +++--- man/{ec.Rd => lc.Rd} | 51 ++-- man/log10_ticks.Rd | 34 +++ man/{predict_mixture.Rd => multi_tox.Rd} | 48 ++-- man/plot_ecxsys.Rd | 39 +-- man/predict_ecxsys.Rd | 22 +- man/stressaddition-package.Rd | 5 +- tests/testthat/test-ecxsys.R | 90 +++---- tests/testthat/{test-ec.R => test-lc.R} | 62 ++--- ...{test-plot_ecxsys.R => test-log10_ticks.R} | 30 ++- tests/testthat/test-multi_tox.R | 117 +++++++++ tests/testthat/test-options.R | 5 +- tests/testthat/test-predict_mixture.R | 90 ------- .../testthat/test-stress_effect_conversion.R | 44 ---- .../test-stress_survival_conversion.R | 44 ++++ 33 files changed, 889 insertions(+), 696 deletions(-) rename R/{ec.R => lc.R} (75%) create mode 100644 R/log10_ticks.R rename R/{predict_mixture.R => multi_tox.R} (65%) rename R/{plot_effect.R => plot_survival.R} (70%) rename R/{stress_effect_conversion.R => stress_survival_conversion.R} (57%) rename man/{ec.Rd => lc.Rd} (63%) create mode 100644 man/log10_ticks.Rd rename man/{predict_mixture.Rd => multi_tox.Rd} (53%) rename tests/testthat/{test-ec.R => test-lc.R} (58%) rename tests/testthat/{test-plot_ecxsys.R => test-log10_ticks.R} (50%) create mode 100644 tests/testthat/test-multi_tox.R delete mode 100644 tests/testthat/test-predict_mixture.R delete mode 100644 tests/testthat/test-stress_effect_conversion.R create mode 100644 tests/testthat/test-stress_survival_conversion.R diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 6215d20..67e72b9 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -1,4 +1,9 @@ -# Using the rocker/tidyverse image because it has devtools installed. +# Using the rocker/tidyverse image because it has devtools and testthat. + +# Use "document = FALSE" to catch potential mismatches between the roxygen +# comments and the .Rd files. This could happen when checking locally +# with "document = TRUE" (the default), which updates the .Rd files, but then +# forgetting to commit and push those updated .Rd files. .check_template: &check_job_template only: @@ -6,15 +11,20 @@ - merge_requests - tags script: - - R -e 'devtools::install_deps()' - - R -e 'devtools::check(error_on = "warning")' + - R -e 'devtools::install_deps(quiet = TRUE)' + - R -e 'devtools::check(error_on = "note", document = FALSE)' + +# Earliest supported version +check-3.5.1: + image: rocker/tidyverse:3.5.1 + <<: *check_job_template -# # latest version of the previous major R release -# check-previous: -# image: rocker/tidyverse:3.5.3 -# <<: *check_job_template +# Last version before 4.0.0 +check-3.6.3: + image: rocker/tidyverse:3.6.3 + <<: *check_job_template -# latest R release +# Latest R release check-latest: image: rocker/tidyverse:latest <<: *check_job_template diff --git a/.gitlab/merge_request_templates/merge_into_master.md b/.gitlab/merge_request_templates/merge_into_master.md index a8b3a26..7790e4c 100644 --- a/.gitlab/merge_request_templates/merge_into_master.md +++ b/.gitlab/merge_request_templates/merge_into_master.md @@ -26,7 +26,7 @@ Comment on any notes that come up during devtools::check(). * [ ] update version in DESCRIPTION * [ ] update date in DESCRIPTION * [ ] no remaining TODO, FIXME, or debug prints anywhere in the source files -* [ ] `devtools::document()` +* [ ] `devtools::document()` and commit the updated .Rd files * [ ] `devtools::check()` without errors or warnings <!-- diff --git a/DESCRIPTION b/DESCRIPTION index e6e6b2c..a607fb5 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: stressaddition Type: Package Title: Modeling Tri-Phasic Concentration-Response Relationships -Version: 2.7.0 -Date: 2020-04-17 +Version: 3.0.0 +Date: 2020-05-07 Authors@R: c(person("Sebastian", "Henz", role = c("aut", "cre"), @@ -25,7 +25,7 @@ URL: https://git.ufz.de/oekotox/stressaddition Encoding: UTF-8 LazyData: true Depends: - R (>= 3.5) + R (> 3.5) Imports: drc (>= 3.0), plotrix diff --git a/NAMESPACE b/NAMESPACE index 58db063..a0b1a5e 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,13 +1,14 @@ # Generated by roxygen2: do not edit by hand -export(ec) export(ecxsys) -export(effect_to_stress) -export(plot_effect) +export(lc) +export(log10_ticks) +export(multi_tox) export(plot_stress) +export(plot_survival) export(predict_ecxsys) -export(predict_mixture) -export(stress_to_effect) +export(stress_to_survival) +export(survival_to_stress) import(grDevices) import(graphics) import(stats) diff --git a/NEWS.md b/NEWS.md index 3dcbf77..413ddf9 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,18 @@ +# stressaddition 3.0.0 + +## Breaking changes +* Renamed all instances of "effect" to "survival". +* Renamed all instances of "ec" to "lc". +* Renamed `predict_mixture()`, which was a temporary development name, to `multi_tox()`. +* The argument `proportion_ca` in the mixture model `multi_tox()` was renamed and its value reversed. It is now called `sa_contribution` and specifies the proportion of stress addition in the calculation of toxicant stress. To convert your code from the old version use this equation: sa_contribution = 1 - proportion_ca. +* Renamed `stress_tox_sam` to `stress_tox_sa` in the output of `multi_tox()`. + +## Bugfixes +* Fixed a bug where `plot_stress()` with argument `which = NULL` would result in an error. Now it correctly draws the axes without data. + +## New +* Exporte function `log10_ticks()` for calculating tick mark labels and positions on a base 10 logarithmic axis. + # stressaddition 2.7.0 * Fixed some spelling mistakes. diff --git a/R/ecxsys.R b/R/ecxsys.R index 5eaed94..43e4448 100644 --- a/R/ecxsys.R +++ b/R/ecxsys.R @@ -21,82 +21,88 @@ # The drc package changes some of the options which some users may have # modified. Of particular interest is options("warn") because it is important # that the user is able to control the visibility of warnings generated by -# ecxsys. Therefore every time drc::drm is called this option is cached +# ecxsys(). Therefore every time drc::drm() is called this option is cached # beforehand and reset immediately afterwards. Additionally, all options are -# cached at the beginning of ecxsys and reset on exit. +# cached at the beginning of ecxsys() and reset on exit. #' ECx-SyS #' -#' The ECx-SyS model for modeling concentration-effect relationships with -#' hormesis. +#' ECx-SyS is a model for tri-phasic concentration-response relationships where +#' hormetic and subhormetic effects are observed at low concentrations. It +#' expands the Stress Addition Model (SAM) by introducing system stress (SyS) +#' which is negatively correlated with toxicant stress. A constant environmental +#' stress can be included. See the publication for details. #' #' It is advised to complete the curve down to zero for optimal prediction. -#' Therefore \code{effect_tox_observed} in the highest concentration should be -#' at or close to zero. If the model does not fit properly try adding an effect +#' Therefore \code{survival_tox_observed} in the highest concentration should be +#' at or close to zero. If the model does not fit properly try adding a survival #' of 0 at ten times the maximum observed concentration. #' -#' The vectors \code{concentration}, \code{effect_tox_observed} and -#' \code{effect_tox_env_observed} (if provided) must be of equal length and +#' The vectors \code{concentration}, \code{survival_tox_observed} and +#' \code{survival_tox_env_observed} (if provided) must be of equal length and #' sorted by increasing concentration. #' #' @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 +#' This is usually the concentration of the highest survival after the +#' control. +#' @param survival_tox_observed A vector of survival 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. Must be between 0 and \code{effect_max}. -#' @param effect_max The maximum value the effect could possibly reach. For +#' between 0 and \code{survival_max}. +#' @param survival_tox_env_observed Survival values observed in the presence of +#' environmental stress. Must be between 0 and \code{survival_max}. +#' @param survival_max The maximum value the survival could possibly reach. For #' survival data in percent this should be 100 (the default). #' @param curves_concentration_max The maximum concentration of the predicted #' curves. This might be useful if for example your highest observed -#' concentration is 30 but you would like to know the predicted values at 100. -#' @param p,q The shape parameters of the beta distribution. Default is 3.2. +#' concentration is 30 but you would like to know the predicted values on a +#' scale between 0 and 100. +#' @param p,q The shape parameters of the beta distribution. Default is +#' \code{p = q = 3.2}. #' #' @return A list (of class ecxsys) containing many different objects of which -#' the most important are listed below. The effect and stress vectors +#' the most important are listed below. The survival 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 +#' \item{survival_tox}{Modeled survival resulting from toxicant stress.} +#' \item{survival_tox_sys}{Modeled survival resulting from toxicant and +#' system stress.} +#' \item{survival_tox_env}{Modeled survival resulting from toxicant and #' environmental stress.} -#' \item{effect_tox_env_sys}{Modeled effect resulting from toxicant, +#' \item{survival_tox_env_sys}{Modeled survival resulting from toxicant, #' environmental and system stress.} -#' \item{effect_tox_LL5}{The effect predicted by the five-parameter +#' \item{survival_tox_LL5}{The survival 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 +#' \item{survival_tox_env_LL5}{The survival 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 +#' \item{curves}{A data frame containing survival 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{concentration_for_plots} is used by the #' plotting functions of this package to approximate the control and -#' generate a break in the concentration axis.} +#' generate a nice concentration axis.} #' } #' #' @examples model <- 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) +#' survival_tox_observed = c(90, 81, 92, 28, 0), +#' survival_tox_env_observed = c(29, 27, 33, 5, 0) #' ) #' -#' # 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: +#' # Use survival_max if for example the survival 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, 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 = 21 +#' survival_tox_observed = c(17, 15.2, 18.8, 7, 0), +#' survival_tox_env_observed = c(4.8, 4.6, 6.4, 0, 0), +#' survival_max = 21 #' ) #' #' @references \href{https://doi.org/10.1038/s41598-019-51645-4}{Liess, M., @@ -106,9 +112,9 @@ #' @export ecxsys <- function(concentration, hormesis_concentration, - effect_tox_observed, - effect_tox_env_observed = NULL, - effect_max = 100, + survival_tox_observed, + survival_tox_env_observed = NULL, + survival_max = 100, curves_concentration_max = NULL, p = 3.2, q = 3.2) { @@ -120,11 +126,11 @@ ecxsys <- function(concentration, # input validation ---------------------------------------------------- - if (effect_max <= 0) { - stop("effect_max must be >= 0") + if (survival_max <= 0) { + stop("survival_max must be >= 0") } - if (length(concentration) != length(effect_tox_observed)) { - stop("concentration and effect_tox_observed must have the same length.") + if (length(concentration) != length(survival_tox_observed)) { + stop("concentration and survival_tox_observed must have the same length.") } if (length(concentration) > length(unique(concentration))) { stop("Concentrations must be unique.") @@ -140,27 +146,27 @@ ecxsys <- function(concentration, stop("hormesis_concentration must be greater than the lowest ", "non-control concentration and less than the highest concentration") } - if (is.null(effect_tox_env_observed)) { + if (is.null(survival_tox_env_observed)) { with_env <- FALSE - effect_tox_env_observed <- rep(NA, length(concentration)) + survival_tox_env_observed <- rep(NA, length(concentration)) } else { with_env <- TRUE } output$with_env <- with_env # Creating all_observations makes it easier to test some assumptions. - all_observations <- effect_tox_observed + all_observations <- survival_tox_observed if (with_env) { - if (length(effect_tox_observed) != length(effect_tox_env_observed)) { - stop("effect_tox_observed and effect_tox_env_observed must have ", - "the same length.") + if (length(survival_tox_observed) != length(survival_tox_env_observed)) { + stop("survival_tox_observed and survival_tox_env_observed ", + "must have the same length.") } - all_observations <- c(all_observations, effect_tox_env_observed) + all_observations <- c(all_observations, survival_tox_env_observed) } 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)) { - stop("Observed effect must be between 0 and effect_max.") + if (any(all_observations > survival_max | all_observations < 0)) { + stop("Observed survival must be between 0 and survival_max.") } conc_shift <- 2 # Powers of ten to shift the control downwards from the # second lowest concentration. This is required to approximate 0 because @@ -177,35 +183,50 @@ ecxsys <- function(concentration, min_conc <- 10 ^ floor(log10(concentration[2]) - conc_shift) - # scale observed effect ----------------------------------------------- - # scale the observed effect to [0,1] to make calculations independent of - # value of the theoretical maximum effect - effect_tox_observed <- effect_tox_observed / effect_max - effect_tox_env_observed <- effect_tox_env_observed / effect_max + # scale observed survival --------------------------------------------- + # scale the observed survival to [0,1] to make calculations independent of + # value of the theoretical maximum survival + survival_tox_observed <- survival_tox_observed / survival_max + survival_tox_env_observed <- survival_tox_env_observed / survival_max # traditional log-logistic model (LL.5) ------------------------------- LL5_tox <- fit_LL5_model(min_conc, concentration, - effect_tox_observed, + survival_tox_observed, original_options) - output$effect_tox_LL5_mod <- LL5_tox$effect_LL5_mod - output$effect_tox_LL5 <- LL5_tox$effect_LL5 * effect_max + output$survival_tox_LL5_mod <- LL5_tox$survival_LL5_mod + output$survival_tox_LL5 <- LL5_tox$survival_LL5 * survival_max if (with_env) { LL5_tox_env <- fit_LL5_model(min_conc, concentration, - effect_tox_env_observed, + survival_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 + output$survival_tox_env_LL5_mod <- LL5_tox_env$survival_LL5_mod + output$survival_tox_env_LL5 <- LL5_tox_env$survival_LL5 * survival_max } # interpolation between subhormesis and hormesis ---------------------- n_new <- 3 # number of new points - concentration <- interpolate(concentration, hormesis_index, n_new, TRUE) - effect_tox_observed <- interpolate(effect_tox_observed, hormesis_index, n_new) + concentration <- interpolate_sub_horm( + concentration, + hormesis_index - 1, + hormesis_index, + n_new, + TRUE + ) + survival_tox_observed <- interpolate_sub_horm( + survival_tox_observed, + hormesis_index - 1, + hormesis_index, + n_new + ) if (with_env) { - effect_tox_env_observed <- interpolate(effect_tox_env_observed, - hormesis_index, n_new) + survival_tox_env_observed <- interpolate_sub_horm( + survival_tox_env_observed, + hormesis_index - 1, + hormesis_index, + n_new + ) } hormesis_index <- hormesis_index + n_new # In the output return only the values at the original concentrations @@ -213,22 +234,22 @@ ecxsys <- function(concentration, keep <- !seq_along(concentration) %in% seq(hormesis_index - n_new, hormesis_index - 1) - # effect_tox ---------------------------------------------------------- - 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()) + # survival_tox -------------------------------------------------------- + survival_tox <- c(1, survival_tox_observed[-1]) + survival_tox[2:(hormesis_index - 1)] <- NA + survival_tox_mod <- drc::drm(survival_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 + survival_tox <- predict(survival_tox_mod, data.frame(concentration = concentration)) + output$survival_tox_mod <- survival_tox_mod + output$survival_tox <- survival_tox[keep] * survival_max - stress_tox <- effect_to_stress(effect_tox, p, q) + stress_tox <- survival_to_stress(survival_tox, p, q) output$stress_tox <- stress_tox[keep] # system stress without environmental stress -------------------------- fit_sys_output <- fit_sys( - effect_to_stress(effect_tox_observed, p, q), + survival_to_stress(survival_tox_observed, p, q), stress_tox, stress_tox, hormesis_index, @@ -246,23 +267,23 @@ ecxsys <- function(concentration, output$sys_tox <- sys_tox - # modeled effect without environmental stress ------------------------- + # modeled survival without environmental stress ----------------------- 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 + survival_tox_sys <- stress_to_survival(stress_tox_sys, p, q) + output$survival_tox_sys <- survival_tox_sys[keep] * survival_max if (with_env) { # env stress ------------------------------------------------------ - stress_tox_env_observed <- effect_to_stress(effect_tox_env_observed, p, q) + stress_tox_env_observed <- survival_to_stress(survival_tox_env_observed, p, q) stress_env <- (stress_tox_env_observed - stress_tox)[hormesis_index] stress_env <- clamp(stress_env) output$stress_env <- stress_env stress_tox_env <- stress_tox + stress_env 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 + survival_tox_env <- stress_to_survival(stress_tox_env, p, q) + output$survival_tox_env <- survival_tox_env[keep] * survival_max # system stress with environmental stress ------------------------- @@ -285,15 +306,15 @@ ecxsys <- function(concentration, output$sys_tox_env <- sys_tox_env - # modeled effect with environmental stress ------------------------ + # modeled survival with environmental stress ---------------------- 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 + survival_tox_env_sys <- stress_to_survival(stress_tox_env_sys, p, q) + output$survival_tox_env_sys <- survival_tox_env_sys[keep] * survival_max } # curves ------------------------------------------------------- - # In order to generate a broken x-axis the concentration vector must + # In order to generate a break in the 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. @@ -312,7 +333,11 @@ ecxsys <- function(concentration, ) output$curves <- predict_ecxsys(output, curves_concentration) - conc_axis <- adjust_plot_concentrations( + # Add a column of concentrations which helps in plotting. It does not make + # sense to show concentrations many orders of magnitude lower than the + # lowest measured concentration. So this cuts out a large chunk and raises + # the concentrations left of the cut so they make a nicer axis. + conc_axis <- make_axis_concentrations( curves_concentration, min_conc, conc_adjust_factor @@ -399,39 +424,45 @@ moving_weighted_mean <- function(x) { fit_LL5_model <- function(min_conc, concentration, - effect_observed, + survival_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) + survival_observed_averaged <- moving_weighted_mean(survival_observed) interpolated <- approx( log10(conc_with_control_shifted), - effect_observed_averaged, + survival_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)) + survival_interpolated <- interpolated$y + survival_LL5_mod <- drc::drm( + survival_interpolated ~ conc_interpolated, + fct = drc::LL.5(fixed = c(NA, 0, survival_observed_averaged[1], NA, NA)) ) options(original_options["warn"]) - effect_LL5 <- predict( - effect_LL5_mod, + survival_LL5 <- predict( + survival_LL5_mod, data.frame(conc_interpolated = concentration) ) list( - effect_LL5_mod = effect_LL5_mod, - effect_LL5 = effect_LL5 + survival_LL5_mod = survival_LL5_mod, + survival_LL5 = survival_LL5 ) } -interpolate <- function(x, to_index, n_new, conc = FALSE) { - from_index <- to_index - 1 # subhormesis_index +interpolate_sub_horm <- function(x, + from_index, + to_index, + n_new, + logarithmic = FALSE) { + # Interpolate between subhormesis and hormesis. This makes the curves + # smoother, less extreme. Without it the slope between subhormesis and + # hormesis would be much steeper in many cases. len <- n_new + 2 # Add 2 because seq() includes the left and right end. - if (conc) { + if (logarithmic) { 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) @@ -440,9 +471,9 @@ interpolate <- function(x, to_index, n_new, conc = FALSE) { } -adjust_plot_concentrations <- function(concentration, - min_conc, - conc_adjust_factor) { +make_axis_concentrations <- function(concentration, + min_conc, + conc_adjust_factor) { # Deals with the concentrations which are unnecessary for plotting. This # means it removes the concentrations in the gap and increases the # concentrations below the gap. diff --git a/R/helpers.R b/R/helpers.R index fcc3052..c56ee33 100644 --- a/R/helpers.R +++ b/R/helpers.R @@ -17,36 +17,12 @@ # along with this program. If not, see <https://www.gnu.org/licenses/>. -# A collection of some internal helper functions. +# A collection of internal helper functions which are small enough that they +# don't desere their own scripts. And which are also not exclusively used by one +# function because then they would be in the script of that function. 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 - ) -} diff --git a/R/ec.R b/R/lc.R similarity index 75% rename from R/ec.R rename to R/lc.R index 9f89821..784ceae 100644 --- a/R/ec.R +++ b/R/lc.R @@ -17,9 +17,9 @@ # along with this program. If not, see <https://www.gnu.org/licenses/>. -#' Effect Concentrations +#' Lethal Concentrations #' -#' Estimate the concentration to reach a certain effect or stress level relative +#' Estimate the concentration to reach a certain mortality relative #' to the control. #' #' If the response level occurs multiple times because of hormesis, which may @@ -27,7 +27,7 @@ #' smallest concentration is returned. #' #' This function only makes sense for curves which generally go down with -#' increasing concentration, i.e. all \code{effect_*} curves and also +#' increasing concentration, i.e. all \code{survival_*} curves and also #' \code{sys_tox} and \code{sys_tox_env}. Others are untested and may give #' unexpected results, if any. #' @@ -36,12 +36,13 @@ #' data frame with a "concentration" column and a \code{response_name} column. #' In the case of the data frame the first row is assumed to be the control. #' 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_name The name of the survival or stress for which you want to +#' calculate the LC. #' @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, +#' and 100. For example with the value 10 the function will return the LC10, #' i.e. the concentration where the response falls below 90 \% of the control -#' response. +#' response. In other words: where a mortality of 10 \% relative to the control +#' is reached. #' @param reference The reference value of the response, usually the control at #' concentration 0. This argument is optional. If it is missing, the first #' value of the response is used as control. This value determines what number @@ -49,49 +50,49 @@ #' response_level is 10 and reference is 45, then the function returns the #' concentration where the curve is reduced by 10\% relative to 45 = 40.5. #' @param warn Logical. Should the function emit a warning if the calculation of -#' the effect concentration is not possible? +#' the lethal concentration is not possible? #' -#' @return A list containing the effect concentration and the corresponding -#' effect. The effect will be \code{NA} if its calculation is impossible from -#' the supplied data. +#' @return A list containing the lethal concentration and the corresponding +#' survival. The survival will be \code{NA} if its calculation is impossible +#' using the supplied data. #' -#' @examples # Calculate the EC10, the concentration where the effect falls -#' # below 90% of the effect in the control. +#' @examples # Calculate the LC10, the concentration where the survival falls +#' # below 90% of the survival in the control. #' #' model <- ecxsys( #' concentration = c(0, 0.05, 0.5, 5, 30), #' hormesis_concentration = 0.5, -#' effect_tox_observed = c(90, 81, 92, 28, 0) +#' survival_tox_observed = c(90, 81, 92, 28, 0) #' ) #' #' # using the ecxsys() output or the curves therein directly: -#' ec(model, "effect_tox_sys", 10) -#' ec(model$curves, "effect_tox_sys", 10) +#' lc(model, "survival_tox_sys", 10) +#' lc(model$curves, "survival_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) +#' lc(curves, "survival_tox_sys", 10) #' #' # using a custom data frame: #' df_custom <- data.frame( #' concentration = curves$concentration, -#' foo = curves$effect_tox_sys +#' foo = curves$survival_tox_sys #' ) -#' ec(df_custom, "foo", 10) +#' lc(df_custom, "foo", 10) #' -#' # Calculate the EC50 relative to an effect of 100 +#' # Calculate the LC50 relative to an survival of 100 #' # instead of relative to the control: -#' ec(model, "effect_tox_sys", 50, reference = 100) +#' lc(model, "survival_tox_sys", 50, reference = 100) #' #' @export -ec <- function(model, +lc <- function(model, response_name, response_level, reference, warn = TRUE) { if (inherits(model, "drc")) { - stop("Please use drc::ED for drc objects.") + stop("Please use drc::ED() for drc objects.") } stopifnot( @@ -122,7 +123,7 @@ ec <- function(model, } else { stopifnot(reference > 0) if (inherits(model, "ecxsys")) { - stopifnot(reference <= model$args$effect_max) + stopifnot(reference <= model$args$survival_max) } reference_value <- reference } @@ -130,15 +131,15 @@ ec <- function(model, if (reference_value == 0) { # May happen with horizontal Sys curves. if (warn) { - warning("Reference value is zero, calculation of EC is impossible.") + warning("Reference value is zero, calculation of LC is impossible.") } - return(list(effect = 0, concentration = NA)) + return(list(response = 0, concentration = NA)) } if (response_level == 0) { - return(list(effect = reference_value, concentration = 0)) + return(list(response = reference_value, concentration = 0)) } else if (response_level == 100) { - return(list(effect = 0, concentration = Inf)) + return(list(response = 0, concentration = Inf)) } response_value <- (1 - response_level / 100) * reference_value @@ -151,12 +152,12 @@ ec <- function(model, response_name, "' are smaller than ", 100 - response_level, - "% of the reference value, which makes determining the EC", + "% of the reference value, which makes determining the LC", response_level, " impossible." ) } - return(list(effect = response_value, concentration = NA)) + return(list(response = response_value, concentration = NA)) } if (!any(response_smaller)) { if (warn) { @@ -165,16 +166,16 @@ ec <- function(model, response_name, "' does not fall below ", 100 - response_level, - "% of the reference, which makes determining the EC", + "% of the reference, which makes determining the LC", response_level, " impossible. You could try using predict_ecxsys() to predict ", "more values in a wider concentration range." ) } - return(list(effect = response_value, concentration = NA)) + return(list(response = response_value, concentration = NA)) } if (response[1] == response_value) { - return(list(effect = response_value, concentration = 0)) + return(list(response = response_value, concentration = 0)) } below <- which(response < response_value)[1] @@ -182,8 +183,8 @@ ec <- function(model, # linear interpolation dist <- (response_value - response[below]) / (response[above] - response[below]) - effect_concentration <- dist * + survival_concentration <- dist * (concentration[above] - concentration[below]) + concentration[below] - list(effect = response_value, concentration = effect_concentration) + list(response = response_value, concentration = survival_concentration) } diff --git a/R/log10_ticks.R b/R/log10_ticks.R new file mode 100644 index 0000000..3f5fd58 --- /dev/null +++ b/R/log10_ticks.R @@ -0,0 +1,72 @@ +# Copyright (C) 2020 Helmholtz-Zentrum fuer Umweltforschung GmbH - UFZ +# See file inst/COPYRIGHTS for details. +# +# This file is part of the R package stressaddition. +# +# stressaddition is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see <https://www.gnu.org/licenses/>. + + +#' Logarithmic axis tick marks +#' +#' Calculate the positions and labels of major and minor tick marks for a base +#' 10 logarithmic axis. +#' +#' @param x A vector of axis values. Can be arbitrarily long but only the +#' minimum and maximum are necessary. +#' @param label_zero Whether or not to replace the smallest major label with +#' "0". This defaults to \code{TRUE} and is useful for some types of plots +#' used to display concentration-response data where the leftmost data point +#' represents the control. +#' +#' @return A list with the positions and labels of the major and minor tick +#' marks. The labels are formatted without trailing zeros using +#' \code{formatC(labels, format = "fg")}. +#' +#' @examples +#' x <- c(0.01, 0.2, 3, 10, 50) +#' plot(x, c(5, 4, 2.5, 1, 0), xaxt = "n", log = "x") +#' ticks <- log10_ticks(x) +#' axis(1, at = ticks$major, labels = ticks$major_labels) +#' axis(1, at = ticks$minor, labels = FALSE, tcl = -0.25) +#' +#' @export +log10_ticks <- function(x, label_zero = TRUE) { + stopifnot(min(x, na.rm = TRUE) > 0) + x <- log10(x) + major <- seq( + floor(min(x, na.rm = TRUE)), + ceiling(max(x, na.rm = TRUE)) + ) + major <- 10 ^ major + 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_labels <- formatC(major, format = "fg") + if (label_zero) { + major_labels[1] <- "0" + } + minor_labels <- formatC(minor, format = "fg") + + list( + major = major, + minor = minor, + major_labels = major_labels, + minor_labels = minor_labels + ) +} diff --git a/R/predict_mixture.R b/R/multi_tox.R similarity index 65% rename from R/predict_mixture.R rename to R/multi_tox.R index ccd9e35..882c069 100644 --- a/R/predict_mixture.R +++ b/R/multi_tox.R @@ -17,10 +17,11 @@ # along with this program. If not, see <https://www.gnu.org/licenses/>. -#' Predict the effect of a mixture of two toxicants +#' Predict the survival of binary toxicant mixtures #' -#' Given the ecxsys models of two toxicants this method predicts the effects of -#' different mixtures of both. +#' The Multi-TOX model predicts the effects of binary toxicant mixtures based on +#' three-phasic concentration-response relationships. See the publication for +#' details. #' #' The predictions are symmetric, i.e. it does not matter which of the toxicant #' models is 1 or 2 as long as the concentration arguments are supplied in the @@ -34,27 +35,28 @@ #' the mixture. Both vectors must either be the same length or the longer #' length must be a multiple of the shorter length. That's because the shorter #' concentration vector gets recycled to the length of the longer one. -#' @param proportion_ca The proportion of concentration addition in the -#' calculation of the toxicant stress of the mixture. Must be between 0 and 1. -#' @param effect_max Controls the scaling of the result. This represents the -#' maximum value the effect could possibly reach. For survival data in percent -#' this should be 100 (the default). +#' @param sa_contribution The proportion of stress addition contributing to the +#' calculation of the toxicant stress in the mixture. Must be between 0 and 1 +#' where 1 stands for 100 \% stress addition. +#' @param survival_max Controls the scaling of the result. This represents the +#' maximum value the survival could possibly reach. For survival data in +#' percent this should be 100 (the default). #' #' @return A data frame with columns of the supplied concentrations and the -#' corresponding mixture effects and stresses. +#' corresponding mixture survival and stresses. #' #' @examples toxicant_a <- ecxsys( #' concentration = c(0, 0.05, 0.5, 5, 30), #' hormesis_concentration = 0.5, -#' effect_tox_observed = c(90, 81, 92, 28, 0), +#' survival_tox_observed = c(90, 81, 92, 28, 0), #' ) #' toxicant_b <- 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 +#' survival_tox_observed = c(26, 25, 24, 27, 5, 0), +#' survival_max = 30 #' ) -#' predict_mixture( +#' multi_tox( #' toxicant_a , #' toxicant_b , #' c(0, 0.02, 0.2, 2, 20), @@ -64,18 +66,18 @@ #' # Example of symmetric prediction: #' conc_a <- c(0, 0.03, 0.3, 3) #' conc_b <- 5.5 -#' prop_ca <- 0.75 -#' mix_a <- predict_mixture(toxicant_a , toxicant_b , conc_a, conc_b, prop_ca) -#' mix_b <- predict_mixture(toxicant_b , toxicant_a , conc_b, conc_a, prop_ca) -#' identical(mix_a$effect, mix_b$effect) +#' sa_contrib <- 0.75 +#' mix_a <- multi_tox(toxicant_a , toxicant_b , conc_a, conc_b, sa_contrib) +#' mix_b <- multi_tox(toxicant_b , toxicant_a , conc_b, conc_a, sa_contrib) +#' identical(mix_a$survival, mix_b$survival) #' #' @export -predict_mixture <- function(model_a, - model_b, - concentration_a, - concentration_b, - proportion_ca = 0.5, - effect_max = 100) { +multi_tox <- function(model_a, + model_b, + concentration_a, + concentration_b, + sa_contribution = 0.5, + survival_max = 100) { stopifnot( inherits(model_a, "ecxsys"), inherits(model_b, "ecxsys"), @@ -85,8 +87,8 @@ predict_mixture <- function(model_a, length(concentration_b) > 0, all(!is.na(concentration_a)), all(!is.na(concentration_b)), - proportion_ca >= 0, - proportion_ca <= 1, + sa_contribution >= 0, + sa_contribution <= 1, model_a$args$p == model_b$args$p, model_a$args$q == model_b$args$q ) @@ -95,29 +97,29 @@ predict_mixture <- function(model_a, predicted_model_b <- predict_ecxsys(model_b, concentration_b) # tox stress ---------------------------------------------------------- - stress_tox_sam <- predicted_model_a$stress_tox + predicted_model_b$stress_tox + stress_tox_sa <- predicted_model_a$stress_tox + predicted_model_b$stress_tox # Convert the model_b concentration into an equivalent model_a concentration # and vice versa. concentration_b_equivalent <- W1.2_inverse( - model_a$effect_tox_mod, - predicted_model_b$effect_tox / model_b$args$effect_max + model_a$survival_tox_mod, + predicted_model_b$survival_tox / model_b$args$survival_max ) - effect_tox_ca_a <- predict( - model_a$effect_tox_mod, + survival_tox_ca_a <- predict( + model_a$survival_tox_mod, data.frame(concentration = concentration_a + concentration_b_equivalent) ) - stress_tox_ca_a <- effect_to_stress(effect_tox_ca_a) + stress_tox_ca_a <- survival_to_stress(survival_tox_ca_a) concentration_a_equivalent <- W1.2_inverse( - model_b$effect_tox_mod, - predicted_model_a$effect_tox / model_a$args$effect_max + model_b$survival_tox_mod, + predicted_model_a$survival_tox / model_a$args$survival_max ) - effect_tox_ca_b <- predict( - model_b$effect_tox_mod, + survival_tox_ca_b <- predict( + model_b$survival_tox_mod, data.frame(concentration = concentration_b + concentration_a_equivalent) ) - stress_tox_ca_b <- effect_to_stress(effect_tox_ca_b) + stress_tox_ca_b <- survival_to_stress(survival_tox_ca_b) stress_tox_ca <- (stress_tox_ca_a + stress_tox_ca_b) / 2 @@ -127,16 +129,16 @@ predict_mixture <- function(model_a, sys_total <- (sys_a + sys_b) / 2 # combined stress and result ------------------------------------------ - proportion_sam <- 1 - proportion_ca - stress_tox_total <- stress_tox_ca * proportion_ca + stress_tox_sam * proportion_sam + ca_contribution <- 1 - sa_contribution + stress_tox_total <- stress_tox_sa * sa_contribution + stress_tox_ca * ca_contribution stress_total <- stress_tox_total + sys_total - effect <- stress_to_effect(stress_total) * effect_max + survival <- stress_to_survival(stress_total) * survival_max data.frame( concentration_a = concentration_a, concentration_b = concentration_b, - effect = effect, - stress_tox_sam = stress_tox_sam, + survival = survival, + stress_tox_sa = stress_tox_sa, stress_tox_ca = stress_tox_ca, stress_tox = stress_tox_total, sys = sys_total, diff --git a/R/plot_ecxsys.R b/R/plot_ecxsys.R index e3d1f3e..f0d9c0f 100644 --- a/R/plot_ecxsys.R +++ b/R/plot_ecxsys.R @@ -22,18 +22,18 @@ #' Plot the results of the ECx-SyS model #' -#' Plot the observed and modeled effects and stresses. +#' Plot the observed and modeled survivals and stresses. #' #' @name plot_ecxsys #' #' @param model The list returned from \code{\link{ecxsys}}. #' @param which A vector of names to plot. Allowed are the column names of the -#' \code{model$curves} data frame. There is also \code{"effect_tox_observed"} -#' and \code{"effect_tox_env_observed"} for the observed effects and -#' \code{"sys_tox_observed"} and \code{"sys_tox_env_observed"} for the -#' observed Sys. The default \code{NA} only plots the most important curves. -#' Use \code{which = "all"} to display all curves. An empty vector or -#' \code{NULL} creates just the axes. +#' \code{model$curves} data frame. There is also +#' \code{"survival_tox_observed"} and \code{"survival_tox_env_observed"} for +#' the observed survival and \code{"sys_tox_observed"} and +#' \code{"sys_tox_env_observed"} for the observed Sys. The default \code{NA} +#' only plots the most important curves. Use \code{which = "all"} to display +#' all curves. An empty vector or \code{NULL} creates just the axes. #' #' @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 @@ -43,21 +43,21 @@ #' @examples model <- 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) +#' survival_tox_observed = c(90, 81, 92, 28, 0), +#' survival_tox_env_observed = c(29, 27, 33, 5, 0) #' ) -#' plot_effect(model, show_legend = TRUE) +#' plot_survival(model, show_legend = TRUE) #' plot_stress(model, show_legend = TRUE) #' #' # Plot all curves: -#' plot_effect(model, which = "all") +#' plot_survival(model, which = "all") #' plot_stress(model, which = "all") #' #' # Plot only some selected curves: -#' plot_effect(model, which = c("effect_tox_sys", "effect_tox_env_sys")) +#' plot_survival(model, which = c("survival_tox_sys", "survival_tox_env_sys")) #' plot_stress(model, which = c("sys_tox", "sys_tox_env")) #' #' # Plot only the observed values: -#' plot_effect(model, which = c("effect_tox_observed", "effect_tox_env_observed")) +#' plot_survival(model, which = c("survival_tox_observed", "survival_tox_env_observed")) #' plot_stress(model, which = c("sys_tox_observed", "sys_tox_env_observed")) NULL diff --git a/R/plot_stress.R b/R/plot_stress.R index 7c2f5da..1d702a9 100644 --- a/R/plot_stress.R +++ b/R/plot_stress.R @@ -44,7 +44,7 @@ plot_stress <- function(model, which <- valid_names } else if (!model$with_env && any(grepl("env", which, fixed = TRUE))) { warning("'which' contains names with 'env' but the model was built ", - "without environmental effects.") + "without environmental stress.") which <- which[which %in% valid_names] } else if (any(!which %in% valid_names)) { warning("Argument 'which' contains invalid names.") @@ -52,14 +52,25 @@ plot_stress <- function(model, } curves <- model$curves - log_ticks <- get_log_ticks(curves$concentration_for_plots) + ticks <- log10_ticks(curves$concentration_for_plots) point_concentration <- c( curves$concentration_for_plots[1], model$args$concentration[-1] ) - curves_w <- curves[, which[!endsWith(which, "observed")]] - ymax <- if (NCOL(curves_w) == 0) 1 else max(curves_w, 1, na.rm = TRUE) + if (is.null(which)) { + ymax = 1 + } else { + which_lines <- which[!endsWith(which, "observed")] + if (length(which_lines) == 0) { + ymax <- 1 + } else { + ymax <- max(curves[, which_lines], 1, na.rm = TRUE) + # No need to include the observed stress in the call to max() no + # matter if those are in "which" or not because these vectors are + # clamped to [0, 1] anyway. + } + } plot( NA, @@ -151,9 +162,9 @@ plot_stress <- function(model, # 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, + axis(1, at = ticks$major, labels = ticks$major_labels, col = NA, col.ticks = par("fg")) - axis(1, at = log_ticks$minor, labels = FALSE, tcl = -0.25, + axis(1, at = ticks$minor, labels = FALSE, tcl = -0.25, col = NA, col.ticks = par("fg")) plotrix::axis.break(1, breakpos = model$axis_break_conc) axis(2, col = NA, col.ticks = par("fg"), las = 1) diff --git a/R/plot_effect.R b/R/plot_survival.R similarity index 70% rename from R/plot_effect.R rename to R/plot_survival.R index 0fb0293..7963504 100644 --- a/R/plot_effect.R +++ b/R/plot_survival.R @@ -19,24 +19,24 @@ #' @rdname plot_ecxsys #' @export -plot_effect <- function(model, - which = NA, - show_legend = FALSE, - xlab = "concentration", - ylab = "effect", - main = NULL) { +plot_survival <- function(model, + which = NA, + show_legend = FALSE, + xlab = "concentration", + ylab = "survival", + main = NULL) { stopifnot(inherits(model, "ecxsys")) curve_names <- names(model$curves) valid_names <- c( - curve_names[startsWith(curve_names, "effect")], - "effect_tox_observed", "effect_tox_env_observed" # the observed points + curve_names[startsWith(curve_names, "survival")], + "survival_tox_observed", "survival_tox_env_observed" # observed points ) if (length(which) == 1 && is.na(which)) { - which <- c("effect_tox", "effect_tox_sys", "effect_tox_observed") + which <- c("survival_tox", "survival_tox_sys", "survival_tox_observed") if (model$with_env) { - which <- c(which, "effect_tox_env", "effect_tox_env_sys", - "effect_tox_env_observed") + which <- c(which, "survival_tox_env", "survival_tox_env_sys", + "survival_tox_env_observed") } } else if ("all" %in% which) { if (length(which) > 1) { @@ -45,7 +45,7 @@ plot_effect <- function(model, which <- valid_names } else if (!model$with_env && any(grepl("env", which, fixed = TRUE))) { warning("'which' contains names with 'env' but the model was built ", - "without environmental effects.") + "without environmental stress.") which <- which[which %in% valid_names] } else if (any(!which %in% valid_names)) { warning("Argument 'which' contains invalid names.") @@ -53,7 +53,7 @@ plot_effect <- function(model, } curves <- model$curves - log_ticks <- get_log_ticks(curves$concentration_for_plots) + ticks <- log10_ticks(curves$concentration_for_plots) point_concentration <- c( curves$concentration_for_plots[1], model$args$concentration[-1] @@ -63,7 +63,7 @@ plot_effect <- function(model, NA, NA, xlim = range(curves$concentration_for_plots, na.rm = TRUE), - ylim = c(0, model$args$effect_max), + ylim = c(0, model$args$survival_max), log = "x", xlab = xlab, ylab = ylab, @@ -75,65 +75,65 @@ plot_effect <- function(model, # 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_observed" %in% which) { + if ("survival_tox_observed" %in% which) { points( point_concentration, - model$args$effect_tox_observed, + model$args$survival_tox_observed, pch = 16, col = "blue" ) } - if ("effect_tox_sys" %in% which) { + if ("survival_tox_sys" %in% which) { lines( curves$concentration_for_plots, - curves$effect_tox_sys, + curves$survival_tox_sys, col = "blue" ) } - if ("effect_tox" %in% which) { + if ("survival_tox" %in% which) { lines( curves$concentration_for_plots, - curves$effect_tox, + curves$survival_tox, col = "deepskyblue", lty = 2 ) } - if ("effect_tox_LL5" %in% which) { + if ("survival_tox_LL5" %in% which) { lines( curves$concentration_for_plots, - curves$effect_tox_LL5, + curves$survival_tox_LL5, col = "darkblue", lty = 3 ) } if (model$with_env) { - if ("effect_tox_env_observed" %in% which) { + if ("survival_tox_env_observed" %in% which) { points( point_concentration, - model$args$effect_tox_env_observed, + model$args$survival_tox_env_observed, pch = 16, col = "red" ) } - if ("effect_tox_env_sys" %in% which) { + if ("survival_tox_env_sys" %in% which) { lines( curves$concentration_for_plots, - curves$effect_tox_env_sys, + curves$survival_tox_env_sys, col = "red" ) } - if ("effect_tox_env" %in% which) { + if ("survival_tox_env" %in% which) { lines( curves$concentration_for_plots, - curves$effect_tox_env, + curves$survival_tox_env, col = "orange", lty = 2 ) } - if ("effect_tox_env_LL5" %in% which) { + if ("survival_tox_env_LL5" %in% which) { lines( curves$concentration_for_plots, - curves$effect_tox_env_LL5, + curves$survival_tox_env_LL5, col = "darkred", lty = 3 ) @@ -143,18 +143,18 @@ plot_effect <- function(model, # 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, + axis(1, at = ticks$major, labels = ticks$major_labels, col = NA, col.ticks = par("fg")) - axis(1, at = log_ticks$minor, labels = FALSE, tcl = -0.25, + axis(1, at = ticks$minor, labels = FALSE, tcl = -0.25, col = NA, col.ticks = par("fg")) plotrix::axis.break(1, breakpos = model$axis_break_conc) axis(2, col = NA, col.ticks = par("fg"), las = 1) if (show_legend) { legend_df <- data.frame( - name = c("effect_tox_observed", "effect_tox", "effect_tox_sys", - "effect_tox_LL5", "effect_tox_env_observed", - "effect_tox_env", "effect_tox_env_sys", "effect_tox_env_LL5"), + name = c("survival_tox_observed", "survival_tox", "survival_tox_sys", + "survival_tox_LL5", "survival_tox_env_observed", + "survival_tox_env", "survival_tox_env_sys", "survival_tox_env_LL5"), text = c("tox (observed)", "tox", "tox + sys", "tox (LL5)", "tox + env (observed)", "tox + env", "tox + env + sys", "tox + env (LL5)"), diff --git a/R/predict_ecxsys.R b/R/predict_ecxsys.R index 567ac67..28e06da 100644 --- a/R/predict_ecxsys.R +++ b/R/predict_ecxsys.R @@ -17,9 +17,9 @@ # along with this program. If not, see <https://www.gnu.org/licenses/>. -#' Predict effects and stresses +#' Predict survival and stress #' -#' Calculate the effects and stresses of an ECx-SyS model at arbitrary +#' Calculate the survivals and stresses of an ECx-SyS model at arbitrary #' concentrations. #' #' @param model An ECx-SyS model as returned by \code{\link{ecxsys}}. @@ -29,23 +29,23 @@ #' columns: #' \describe{ #' \item{concentration}{The supplied concentrations.} -#' \item{effect_tox_LL5}{The effect predicted by the five-parameter +#' \item{survival_tox_LL5}{The survival 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{survival_tox}{Modeled survival resulting from toxicant stress.} +#' \item{survival_tox_sys}{Modeled survival 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 effect predicted by the five-parameter +#' \item{survival_tox_env_LL5}{The survival 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 +#' \item{survival_tox_env}{Modeled survival resulting from toxicant and #' environmental stress.} -#' \item{effect_tox_env_sys}{Modeled effect resulting from toxicant, +#' \item{survival_tox_env_sys}{Modeled survival resulting from toxicant, #' environmental and system stress.} #' \item{stress_env}{Environmental stress.} #' \item{stress_tox_env}{The sum of toxicant and environmental stress.} @@ -58,8 +58,8 @@ #' @examples model <- 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) +#' survival_tox_observed = c(90, 81, 92, 28, 0), +#' survival_tox_env_observed = c(29, 27, 33, 5, 0) #' ) #' p <- predict_ecxsys(model, c(0.001, 0.01, 0.1, 1, 10)) #' @@ -73,30 +73,30 @@ predict_ecxsys <- function(model, concentration) { ) p <- model$args$p q <- model$args$q - effect_max <- model$args$effect_max + survival_max <- model$args$survival_max out_df <- data.frame( concentration = concentration ) - out_df$effect_tox_LL5 <- predict( - model$effect_tox_LL5_mod, + out_df$survival_tox_LL5 <- predict( + model$survival_tox_LL5_mod, data.frame(concentration = concentration) - ) * effect_max + ) * survival_max if (model$with_env) { - out_df$effect_tox_env_LL5 <- predict( - model$effect_tox_env_LL5_mod, + out_df$survival_tox_env_LL5 <- predict( + model$survival_tox_env_LL5_mod, data.frame(concentration = concentration) - ) * effect_max + ) * survival_max } - effect_tox <- predict( - model$effect_tox_mod, + survival_tox <- predict( + model$survival_tox_mod, data.frame(concentration = concentration) ) - out_df$effect_tox <- effect_tox * effect_max + out_df$survival_tox <- survival_tox * survival_max - stress_tox <- effect_to_stress(effect_tox, p, q) + stress_tox <- survival_to_stress(survival_tox, p, q) out_df$stress_tox <- stress_tox sys_tox <- predict( @@ -108,8 +108,8 @@ predict_ecxsys <- function(model, concentration) { 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 + survival_tox_sys <- stress_to_survival(stress_tox_sys, p, q) + out_df$survival_tox_sys <- survival_tox_sys * survival_max if (model$with_env) { out_df$stress_env <- model$stress_env @@ -117,9 +117,9 @@ predict_ecxsys <- function(model, concentration) { stress_tox_env <- stress_tox + model$stress_env out_df$stress_tox_env <- stress_tox_env - out_df$effect_tox_env <- stress_to_effect( + out_df$survival_tox_env <- stress_to_survival( stress_tox_env, p, q - ) * effect_max + ) * survival_max sys_tox_env <- predict( model$sys_tox_env_mod, @@ -131,10 +131,10 @@ predict_ecxsys <- function(model, concentration) { sys_tox_env out_df$stress_tox_env_sys <- stress_tox_env_sys - effect_tox_env_sys <- stress_to_effect( + survival_tox_env_sys <- stress_to_survival( stress_tox_env_sys, p, q ) - out_df$effect_tox_env_sys <- effect_tox_env_sys * effect_max + out_df$survival_tox_env_sys <- survival_tox_env_sys * survival_max } class(out_df) <- c("ecxsys_predicted", class(out_df)) diff --git a/R/stress_effect_conversion.R b/R/stress_survival_conversion.R similarity index 57% rename from R/stress_effect_conversion.R rename to R/stress_survival_conversion.R index 49dcc5b..50e614e 100644 --- a/R/stress_effect_conversion.R +++ b/R/stress_survival_conversion.R @@ -17,50 +17,46 @@ # along with this program. If not, see <https://www.gnu.org/licenses/>. -#' Convert Between Stress and Effect +#' Convert Between Stress and Survival #' -#' Functions to convert effect to general stress or vice versa using the beta +#' Functions to convert survival to general stress or vice versa using the beta #' 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 clamps the -#' effect to the interval [0, 1] and then returns -#' \code{qbeta(1 - effect, p, q)}. -#' -#' "Effect" is a response which is negatively correlated with stress. For -#' example increasing toxicant concentration (stress) reduces survival (effect). +#' \code{\link[stats:Beta]{qbeta}}. \code{stress_to_survival} returns \code{1 - +#' pbeta(stress, p, q)}. \code{survival_to_stress} first clamps the survival to +#' the interval [0, 1] and then returns \code{qbeta(1 - survival, p, q)}. #' #' @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 stress One or more stress values to convert to effect. +#' @param survival One or more survival 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 survival. #' @param p,q The shape parameters of the \code{\link[stats:Beta]{beta}} #' distribution. Default is 3.2. #' #' @examples #' stress <- 0.3 -#' effect <- stress_to_effect(stress) -#' effect_to_stress(effect) +#' survival <- stress_to_survival(stress) +#' survival_to_stress(survival) #' NULL #' @rdname Stressconversion #' @export -effect_to_stress <- function(effect, p = 3.2, q = 3.2) { +survival_to_stress <- function(survival, p = 3.2, q = 3.2) { stopifnot(p >= 0, q >= 0) - effect <- clamp(effect) - qbeta(1 - effect, p, q) + survival <- clamp(survival) + qbeta(1 - survival, p, q) } #' @rdname Stressconversion #' @export -stress_to_effect <- function(stress, p = 3.2, q = 3.2) { +stress_to_survival <- 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 abd7efa..e36e1f8 100644 --- a/R/stressaddition-package.R +++ b/R/stressaddition-package.R @@ -31,9 +31,8 @@ #' equations. #' #' In the paper the model is introduced in the context of survival -#' experiments. However, it can also be used to model other concentration -#' dependent responses. For this reason the more general term "effect" instead -#' of "survival" is used throughout the package. +#' experiments. However, we conjecture that it can also be used to model +#' other concentration dependent responses. #' #' This package is not on CRAN. It is hosted on the UFZ GitLab server. Visit #' the repository linked below for the newest version. See the readme file diff --git a/README.md b/README.md index 426439e..0a1c02a 100644 --- a/README.md +++ b/README.md @@ -1,20 +1,21 @@ # stressaddition This is the R implementation of the tri-phasic concentration-response model introduced in -[Liess, M., Henz, S. & Knillmann, S. Predicting low-concentration effects of pesticides. Sci Rep 9, 15248 (2019).](https://doi.org/10.1038/s41598-019-51645-4) -It allows modeling of ecotoxicological experiments where the response shows signs of a hormesis effect. +[Liess, M., Henz, S. & Knillmann, S. Predicting low-concentration effects of pesticides. Sci Rep 9, 15248 (2019)](https://doi.org/10.1038/s41598-019-51645-4). It allows modeling of ecotoxicological experiments where the response shows signs of a hormesis effect. + +The EC<sub>x-SyS</sub> and Multi-TOX models from this package are also available as part of the [Indicate app](http://www.systemecology.eu/indicate) which offers a graphical user interface. ## Installation -Stressaddition is not on CRAN. You can install the most recent stable version from GitLab using the devtools package: +Stressaddition is not yet on CRAN. You can install the most recent stable version from GitLab using the remotes package: ``` r -# install.packages("devtools") -devtools::install_gitlab("oekotox/stressaddition", host = "git.ufz.de") +install.packages("remotes") +remotes::install_gitlab("oekotox/stressaddition", host = "git.ufz.de") ``` -Alternatively there are binary and source builds downloadable from the [releases page](https://git.ufz.de/oekotox/stressaddition/-/releases). +Alternatively there are binary and source builds of various versions downloadable from the [releases page](https://git.ufz.de/oekotox/stressaddition/-/releases). ## Updating RStudio's integrated package updater won't detect updates in packages installed from GitHub or GitLab. I recommend running ```r -devtools::update_packages() +remotes::update_packages() ``` in regular intervals to check for updates from those sources. @@ -22,26 +23,25 @@ in regular intervals to check for updates from those sources. Please cite this package if you use it in your analysis. See `citation("stressaddition")` for details. ## Example -In the paper we use the model in the context of survival experiments. However, it can also be applicable in modeling other concentration or dose dependent responses. For this reason the more general term "effect" instead of "survival" is used throughout the package. ```r library(stressaddition) model <- 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) + survival_tox_observed = c(90, 81, 92, 28, 0), + survival_tox_env_observed = c(29, 27, 33, 5, 0) ) # Plot the effect and the system stress: par(mfrow = c(2, 1)) -plot_effect(model) +plot_survival(model) plot_stress(model) -# The LC50 of the effect under the influence of toxicant and system tress: -ec(model, "effect_tox_sys", 50) +# The LC50 under the influence of toxicant and system tress: +lc(model, "survival_tox_sys", 50) -# The LC10 of the effect under the influence of toxicant, system and environmental tress: -ec(model, "effect_tox_env_sys", 10) +# The LC10 under the influence of toxicant, system and environmental tress: +lc(model, "survival_tox_env_sys", 10) ``` ## Copyright and License diff --git a/man/Stressconversion.Rd b/man/Stressconversion.Rd index 2762fee..6ffa537 100644 --- a/man/Stressconversion.Rd +++ b/man/Stressconversion.Rd @@ -1,43 +1,39 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/stress_effect_conversion.R +% Please edit documentation in R/stress_survival_conversion.R \name{Stressconversion} \alias{Stressconversion} -\alias{effect_to_stress} -\alias{stress_to_effect} -\title{Convert Between Stress and Effect} +\alias{survival_to_stress} +\alias{stress_to_survival} +\title{Convert Between Stress and Survival} \usage{ -effect_to_stress(effect, p = 3.2, q = 3.2) +survival_to_stress(survival, p = 3.2, q = 3.2) -stress_to_effect(stress, p = 3.2, q = 3.2) +stress_to_survival(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{survival}{One or more survival 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 \code{\link[stats:Beta]{beta}} distribution. Default is 3.2.} -\item{stress}{One or more stress values to convert to effect.} +\item{stress}{One or more stress values to convert to survival.} } \description{ -Functions to convert effect to general stress or vice versa using the beta +Functions to convert survival to general stress or vice versa using the beta distribution. } \details{ 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 clamps the -effect to the interval [0, 1] and then returns -\code{qbeta(1 - effect, p, q)}. - -"Effect" is a response which is negatively correlated with stress. For -example increasing toxicant concentration (stress) reduces survival (effect). +\code{\link[stats:Beta]{qbeta}}. \code{stress_to_survival} returns \code{1 - +pbeta(stress, p, q)}. \code{survival_to_stress} first clamps the survival to +the interval [0, 1] and then returns \code{qbeta(1 - survival, p, q)}. } \examples{ stress <- 0.3 -effect <- stress_to_effect(stress) -effect_to_stress(effect) +survival <- stress_to_survival(stress) +survival_to_stress(survival) } diff --git a/man/ecxsys.Rd b/man/ecxsys.Rd index 42d14bd..71f2568 100644 --- a/man/ecxsys.Rd +++ b/man/ecxsys.Rd @@ -7,9 +7,9 @@ ecxsys( concentration, hormesis_concentration, - effect_tox_observed, - effect_tox_env_observed = NULL, - effect_max = 100, + survival_tox_observed, + survival_tox_env_observed = NULL, + survival_max = 100, curves_concentration_max = NULL, p = 3.2, q = 3.2 @@ -20,81 +20,87 @@ ecxsys( 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.} +This is usually the concentration of the highest survival after the +control.} -\item{effect_tox_observed}{A vector of effect values observed at the given +\item{survival_tox_observed}{A vector of survival values observed at the given concentrations and in absence of environmental stress. Values must be -between 0 and \code{effect_max}.} +between 0 and \code{survival_max}.} -\item{effect_tox_env_observed}{Effect values observed in the presence of -environmental stress. Must be between 0 and \code{effect_max}.} +\item{survival_tox_env_observed}{Survival values observed in the presence of +environmental stress. Must be between 0 and \code{survival_max}.} -\item{effect_max}{The maximum value the effect could possibly reach. For +\item{survival_max}{The maximum value the survival could possibly reach. For survival data in percent this should be 100 (the default).} \item{curves_concentration_max}{The maximum concentration of the predicted curves. This might be useful if for example your highest observed -concentration is 30 but you would like to know the predicted values at 100.} +concentration is 30 but you would like to know the predicted values on a +scale between 0 and 100.} -\item{p, q}{The shape parameters of the beta distribution. Default is 3.2.} +\item{p, q}{The shape parameters of the beta distribution. Default is +\code{p = q = 3.2}.} } \value{ A list (of class ecxsys) containing many different objects of which - the most important are listed below. The effect and stress vectors + the most important are listed below. The survival 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 + \item{survival_tox}{Modeled survival resulting from toxicant stress.} + \item{survival_tox_sys}{Modeled survival resulting from toxicant and + system stress.} + \item{survival_tox_env}{Modeled survival resulting from toxicant and environmental stress.} - \item{effect_tox_env_sys}{Modeled effect resulting from toxicant, + \item{survival_tox_env_sys}{Modeled survival resulting from toxicant, environmental and system stress.} - \item{effect_tox_LL5}{The effect predicted by the five-parameter + \item{survival_tox_LL5}{The survival 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 + \item{survival_tox_env_LL5}{The survival 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 + \item{curves}{A data frame containing survival 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{concentration_for_plots} is used by the plotting functions of this package to approximate the control and - generate a break in the concentration axis.} + generate a nice concentration axis.} } } \description{ -The ECx-SyS model for modeling concentration-effect relationships with -hormesis. +ECx-SyS is a model for tri-phasic concentration-response relationships where +hormetic and subhormetic effects are observed at low concentrations. It +expands the Stress Addition Model (SAM) by introducing system stress (SyS) +which is negatively correlated with toxicant stress. A constant environmental +stress can be included. See the publication for details. } \details{ It is advised to complete the curve down to zero for optimal prediction. -Therefore \code{effect_tox_observed} in the highest concentration should be -at or close to zero. If the model does not fit properly try adding an effect +Therefore \code{survival_tox_observed} in the highest concentration should be +at or close to zero. If the model does not fit properly try adding a survival of 0 at ten times the maximum observed concentration. -The vectors \code{concentration}, \code{effect_tox_observed} and -\code{effect_tox_env_observed} (if provided) must be of equal length and +The vectors \code{concentration}, \code{survival_tox_observed} and +\code{survival_tox_env_observed} (if provided) must be of equal length and sorted by increasing concentration. } \examples{ model <- 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) + survival_tox_observed = c(90, 81, 92, 28, 0), + survival_tox_env_observed = c(29, 27, 33, 5, 0) ) -# 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: +# Use survival_max if for example the survival 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, 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 = 21 + survival_tox_observed = c(17, 15.2, 18.8, 7, 0), + survival_tox_env_observed = c(4.8, 4.6, 6.4, 0, 0), + survival_max = 21 ) } diff --git a/man/ec.Rd b/man/lc.Rd similarity index 63% rename from man/ec.Rd rename to man/lc.Rd index 9ddf03f..6d1cc1f 100644 --- a/man/ec.Rd +++ b/man/lc.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/ec.R -\name{ec} -\alias{ec} -\title{Effect Concentrations} +% Please edit documentation in R/lc.R +\name{lc} +\alias{lc} +\title{Lethal Concentrations} \usage{ -ec(model, response_name, response_level, reference, warn = TRUE) +lc(model, response_name, response_level, reference, warn = TRUE) } \arguments{ \item{model}{This can be one of three types of objects: Either the output of @@ -13,13 +13,14 @@ data frame with a "concentration" column and a \code{response_name} column. In the case of the data frame the first row is assumed to be the control. See the examples.} -\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{response_name}{The name of the survival or stress for which you want to +calculate the LC.} \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, +and 100. For example with the value 10 the function will return the LC10, i.e. the concentration where the response falls below 90 \% of the control -response.} +response. In other words: where a mortality of 10 \% relative to the control +is reached.} \item{reference}{The reference value of the response, usually the control at concentration 0. This argument is optional. If it is missing, the first @@ -29,15 +30,15 @@ response_level is 10 and reference is 45, then the function returns the concentration where the curve is reduced by 10\% relative to 45 = 40.5.} \item{warn}{Logical. Should the function emit a warning if the calculation of -the effect concentration is not possible?} +the lethal concentration is not possible?} } \value{ -A list containing the effect concentration and the corresponding - effect. The effect will be \code{NA} if its calculation is impossible from - the supplied data. +A list containing the lethal concentration and the corresponding + survival. The survival will be \code{NA} if its calculation is impossible + using the supplied data. } \description{ -Estimate the concentration to reach a certain effect or stress level relative +Estimate the concentration to reach a certain mortality relative to the control. } \details{ @@ -46,38 +47,38 @@ happen for low values of \code{response_level}, then the occurrence with the smallest concentration is returned. This function only makes sense for curves which generally go down with -increasing concentration, i.e. all \code{effect_*} curves and also +increasing concentration, i.e. all \code{survival_*} curves and also \code{sys_tox} and \code{sys_tox_env}. Others are untested and may give unexpected results, if any. } \examples{ -# Calculate the EC10, the concentration where the effect falls -# below 90\% of the effect in the control. +# Calculate the LC10, the concentration where the survival falls +# below 90\% of the survival in the control. model <- ecxsys( concentration = c(0, 0.05, 0.5, 5, 30), hormesis_concentration = 0.5, - effect_tox_observed = c(90, 81, 92, 28, 0) + survival_tox_observed = c(90, 81, 92, 28, 0) ) # using the ecxsys() output or the curves therein directly: -ec(model, "effect_tox_sys", 10) -ec(model$curves, "effect_tox_sys", 10) +lc(model, "survival_tox_sys", 10) +lc(model$curves, "survival_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) +lc(curves, "survival_tox_sys", 10) # using a custom data frame: df_custom <- data.frame( concentration = curves$concentration, - foo = curves$effect_tox_sys + foo = curves$survival_tox_sys ) -ec(df_custom, "foo", 10) +lc(df_custom, "foo", 10) -# Calculate the EC50 relative to an effect of 100 +# Calculate the LC50 relative to an survival of 100 # instead of relative to the control: -ec(model, "effect_tox_sys", 50, reference = 100) +lc(model, "survival_tox_sys", 50, reference = 100) } diff --git a/man/log10_ticks.Rd b/man/log10_ticks.Rd new file mode 100644 index 0000000..e83aa0b --- /dev/null +++ b/man/log10_ticks.Rd @@ -0,0 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/log10_ticks.R +\name{log10_ticks} +\alias{log10_ticks} +\title{Logarithmic axis tick marks} +\usage{ +log10_ticks(x, label_zero = TRUE) +} +\arguments{ +\item{x}{A vector of axis values. Can be arbitrarily long but only the +minimum and maximum are necessary.} + +\item{label_zero}{Whether or not to replace the smallest major label with +"0". This defaults to \code{TRUE} and is useful for some types of plots +used to display concentration-response data where the leftmost data point +represents the control.} +} +\value{ +A list with the positions and labels of the major and minor tick + marks. The labels are formatted without trailing zeros using + \code{formatC(labels, format = "fg")}. +} +\description{ +Calculate the positions and labels of major and minor tick marks for a base +10 logarithmic axis. +} +\examples{ +x <- c(0.01, 0.2, 3, 10, 50) +plot(x, c(5, 4, 2.5, 1, 0), xaxt = "n", log = "x") +ticks <- log10_ticks(x) +axis(1, at = ticks$major, labels = ticks$major_labels) +axis(1, at = ticks$minor, labels = FALSE, tcl = -0.25) + +} diff --git a/man/predict_mixture.Rd b/man/multi_tox.Rd similarity index 53% rename from man/predict_mixture.Rd rename to man/multi_tox.Rd index 7c1302d..ebd17c4 100644 --- a/man/predict_mixture.Rd +++ b/man/multi_tox.Rd @@ -1,16 +1,16 @@ % 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} +% Please edit documentation in R/multi_tox.R +\name{multi_tox} +\alias{multi_tox} +\title{Predict the survival of binary toxicant mixtures} \usage{ -predict_mixture( +multi_tox( model_a, model_b, concentration_a, concentration_b, - proportion_ca = 0.5, - effect_max = 100 + sa_contribution = 0.5, + survival_max = 100 ) } \arguments{ @@ -21,20 +21,22 @@ the mixture. Both vectors must either be the same length or the longer length must be a multiple of the shorter length. That's because the shorter concentration vector gets recycled to the length of the longer one.} -\item{proportion_ca}{The proportion of concentration addition in the -calculation of the toxicant stress of the mixture. Must be between 0 and 1.} +\item{sa_contribution}{The proportion of stress addition contributing to the +calculation of the toxicant stress in the mixture. Must be between 0 and 1 +where 1 stands for 100 \% stress addition.} -\item{effect_max}{Controls the scaling of the result. This represents the -maximum value the effect could possibly reach. For survival data in percent -this should be 100 (the default).} +\item{survival_max}{Controls the scaling of the result. This represents the +maximum value the survival could possibly reach. For survival data in +percent this should be 100 (the default).} } \value{ A data frame with columns of the supplied concentrations and the - corresponding mixture effects and stresses. + corresponding mixture survival and stresses. } \description{ -Given the ecxsys models of two toxicants this method predicts the effects of -different mixtures of both. +The Multi-TOX model predicts the effects of binary toxicant mixtures based on +three-phasic concentration-response relationships. See the publication for +details. } \details{ The predictions are symmetric, i.e. it does not matter which of the toxicant @@ -48,15 +50,15 @@ Any environmental stress in \code{model_a} or \code{model_b} is ignored. toxicant_a <- ecxsys( concentration = c(0, 0.05, 0.5, 5, 30), hormesis_concentration = 0.5, - effect_tox_observed = c(90, 81, 92, 28, 0), + survival_tox_observed = c(90, 81, 92, 28, 0), ) toxicant_b <- 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 + survival_tox_observed = c(26, 25, 24, 27, 5, 0), + survival_max = 30 ) -predict_mixture( +multi_tox( toxicant_a , toxicant_b , c(0, 0.02, 0.2, 2, 20), @@ -66,9 +68,9 @@ predict_mixture( # Example of symmetric prediction: conc_a <- c(0, 0.03, 0.3, 3) conc_b <- 5.5 -prop_ca <- 0.75 -mix_a <- predict_mixture(toxicant_a , toxicant_b , conc_a, conc_b, prop_ca) -mix_b <- predict_mixture(toxicant_b , toxicant_a , conc_b, conc_a, prop_ca) -identical(mix_a$effect, mix_b$effect) +sa_contrib <- 0.75 +mix_a <- multi_tox(toxicant_a , toxicant_b , conc_a, conc_b, sa_contrib) +mix_b <- multi_tox(toxicant_b , toxicant_a , conc_b, conc_a, sa_contrib) +identical(mix_a$survival, mix_b$survival) } diff --git a/man/plot_ecxsys.Rd b/man/plot_ecxsys.Rd index 1b5e970..0cb3187 100644 --- a/man/plot_ecxsys.Rd +++ b/man/plot_ecxsys.Rd @@ -1,26 +1,27 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/plot_ecxsys.R, R/plot_effect.R, R/plot_stress.R +% Please edit documentation in R/plot_ecxsys.R, R/plot_stress.R, +% R/plot_survival.R \name{plot_ecxsys} \alias{plot_ecxsys} -\alias{plot_effect} \alias{plot_stress} +\alias{plot_survival} \title{Plot the results of the ECx-SyS model} \usage{ -plot_effect( +plot_stress( model, which = NA, show_legend = FALSE, xlab = "concentration", - ylab = "effect", + ylab = "stress", main = NULL ) -plot_stress( +plot_survival( model, which = NA, show_legend = FALSE, xlab = "concentration", - ylab = "stress", + ylab = "survival", main = NULL ) } @@ -28,12 +29,12 @@ plot_stress( \item{model}{The list returned from \code{\link{ecxsys}}.} \item{which}{A vector of names to plot. Allowed are the column names of the -\code{model$curves} data frame. There is also \code{"effect_tox_observed"} -and \code{"effect_tox_env_observed"} for the observed effects and -\code{"sys_tox_observed"} and \code{"sys_tox_env_observed"} for the -observed Sys. The default \code{NA} only plots the most important curves. -Use \code{which = "all"} to display all curves. An empty vector or -\code{NULL} creates just the axes.} +\code{model$curves} data frame. There is also +\code{"survival_tox_observed"} and \code{"survival_tox_env_observed"} for +the observed survival and \code{"sys_tox_observed"} and +\code{"sys_tox_env_observed"} for the observed Sys. The default \code{NA} +only plots the most important curves. Use \code{which = "all"} to display +all curves. An empty vector or \code{NULL} creates just the axes.} \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 @@ -42,27 +43,27 @@ the number of elements shown.} \item{xlab, ylab, main}{Axis labels and title.} } \description{ -Plot the observed and modeled effects and stresses. +Plot the observed and modeled survivals and stresses. } \examples{ model <- 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) + survival_tox_observed = c(90, 81, 92, 28, 0), + survival_tox_env_observed = c(29, 27, 33, 5, 0) ) -plot_effect(model, show_legend = TRUE) +plot_survival(model, show_legend = TRUE) plot_stress(model, show_legend = TRUE) # Plot all curves: -plot_effect(model, which = "all") +plot_survival(model, which = "all") plot_stress(model, which = "all") # Plot only some selected curves: -plot_effect(model, which = c("effect_tox_sys", "effect_tox_env_sys")) +plot_survival(model, which = c("survival_tox_sys", "survival_tox_env_sys")) plot_stress(model, which = c("sys_tox", "sys_tox_env")) # Plot only the observed values: -plot_effect(model, which = c("effect_tox_observed", "effect_tox_env_observed")) +plot_survival(model, which = c("survival_tox_observed", "survival_tox_env_observed")) plot_stress(model, which = c("sys_tox_observed", "sys_tox_env_observed")) } diff --git a/man/predict_ecxsys.Rd b/man/predict_ecxsys.Rd index 5e86de0..2039c7c 100644 --- a/man/predict_ecxsys.Rd +++ b/man/predict_ecxsys.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/predict_ecxsys.R \name{predict_ecxsys} \alias{predict_ecxsys} -\title{Predict effects and stresses} +\title{Predict survival and stress} \usage{ predict_ecxsys(model, concentration) } @@ -16,23 +16,23 @@ A data frame (of class "ecxsys_predicted") with the following columns: \describe{ \item{concentration}{The supplied concentrations.} - \item{effect_tox_LL5}{The effect predicted by the five-parameter + \item{survival_tox_LL5}{The survival 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{survival_tox}{Modeled survival resulting from toxicant stress.} + \item{survival_tox_sys}{Modeled survival 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 effect predicted by the five-parameter + \item{survival_tox_env_LL5}{The survival 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 + \item{survival_tox_env}{Modeled survival resulting from toxicant and environmental stress.} - \item{effect_tox_env_sys}{Modeled effect resulting from toxicant, + \item{survival_tox_env_sys}{Modeled survival resulting from toxicant, environmental and system stress.} \item{stress_env}{Environmental stress.} \item{stress_tox_env}{The sum of toxicant and environmental stress.} @@ -43,15 +43,15 @@ A data frame (of class "ecxsys_predicted") with the following } } \description{ -Calculate the effects and stresses of an ECx-SyS model at arbitrary +Calculate the survivals and stresses of an ECx-SyS model at arbitrary concentrations. } \examples{ model <- 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) + survival_tox_observed = c(90, 81, 92, 28, 0), + survival_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/stressaddition-package.Rd b/man/stressaddition-package.Rd index 5a4ca3c..9908cea 100644 --- a/man/stressaddition-package.Rd +++ b/man/stressaddition-package.Rd @@ -17,9 +17,8 @@ See the publication linked below for more information including equations. In the paper the model is introduced in the context of survival - experiments. However, it can also be used to model other concentration - dependent responses. For this reason the more general term "effect" instead - of "survival" is used throughout the package. + experiments. However, we conjecture that it can also be used to model + other concentration dependent responses. This package is not on CRAN. It is hosted on the UFZ GitLab server. Visit the repository linked below for the newest version. See the readme file diff --git a/tests/testthat/test-ecxsys.R b/tests/testthat/test-ecxsys.R index ffc2608..cd880a7 100644 --- a/tests/testthat/test-ecxsys.R +++ b/tests/testthat/test-ecxsys.R @@ -20,8 +20,8 @@ mod <- 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) + survival_tox_observed = c(90, 81, 92, 28, 0), + survival_tox_env_observed = c(29, 27, 33, 5, 0) ) @@ -31,7 +31,7 @@ test_that("error when hormesis_concentration not in concentration", { ecxsys( concentration = c(0, 0.05, 0.5, 5, 30), hormesis_concentration = 0.4, - effect_tox_observed = c(90, 81, 92, 28, 0) + survival_tox_observed = c(90, 81, 92, 28, 0) ), errorm, fixed = TRUE @@ -40,7 +40,7 @@ test_that("error when hormesis_concentration not in concentration", { ecxsys( concentration = c(0, 0.05, 0.5, 5, 30), hormesis_concentration = 20, - effect_tox_observed = c(90, 81, 92, 28, 0) + survival_tox_observed = c(90, 81, 92, 28, 0) ), errorm, fixed = TRUE @@ -57,7 +57,7 @@ test_that("error when hormesis_index <= 2 or >= (length(concentration))", { ecxsys( concentration = c(0, 0.05, 0.5, 5, 30), hormesis_concentration = 0, - effect_tox_observed = c(85, 76, 94, 35, 0) + survival_tox_observed = c(85, 76, 94, 35, 0) ), errorm, fixed = TRUE @@ -66,7 +66,7 @@ test_that("error when hormesis_index <= 2 or >= (length(concentration))", { ecxsys( concentration = c(0, 0.05, 0.5, 5, 30), hormesis_concentration = 0.05, - effect_tox_observed = c(85, 76, 94, 35, 0) + survival_tox_observed = c(85, 76, 94, 35, 0) ), errorm, fixed = TRUE @@ -75,7 +75,7 @@ test_that("error when hormesis_index <= 2 or >= (length(concentration))", { ecxsys( concentration = c(0, 0.05, 0.5, 5, 30), hormesis_concentration = 30, - effect_tox_observed = c(85, 76, 94, 35, 0) + survival_tox_observed = c(85, 76, 94, 35, 0) ), errorm, fixed = TRUE @@ -90,12 +90,12 @@ test_that("min(concentration) == 0 is shifted the correct amount", { test_that("the discrete results have not changed", { expect_equal( - mod$effect_tox_LL5, + mod$survival_tox_LL5, c(90, 89.745092, 82.340325, 26.787104, 4.306868), tolerance = 1e-4 ) expect_equal( - mod$effect_tox, + mod$survival_tox, c(100, 99.455327884, 92.000028403, 27.999995346, 0.002452598), tolerance = 1e-4 ) @@ -110,13 +110,13 @@ test_that("the discrete results have not changed", { tolerance = 1e-4 ) expect_equal( - mod$effect_tox_sys, + mod$survival_tox_sys, c(89.920669034, 81.890300846, 90.863799559, 27.999995346, 0.002452598), tolerance = 1e-4 ) expect_equal(mod$stress_env, 0.3556369, tolerance = 1e-4) expect_equal( - mod$effect_tox_env_LL5, + mod$survival_tox_env_LL5, c(29.6666667, 29.6666657, 29.5214959, 5.4076411, 0.7871201), tolerance = 1e-4 ) @@ -126,7 +126,7 @@ test_that("the discrete results have not changed", { tolerance = 1e-4 ) expect_equal( - mod$effect_tox_env, + mod$survival_tox_env, c(76.36558967, 59.90195211, 33, 0.01079038, 0), tolerance = 1e-4 ) @@ -136,7 +136,7 @@ test_that("the discrete results have not changed", { tolerance = 1e-4 ) expect_equal( - mod$effect_tox_env_sys, + mod$survival_tox_env_sys, c(29.37633973, 25.99861674, 30.16676004, 0.01079038, 0), tolerance = 1e-4 ) @@ -148,19 +148,19 @@ test_that("the curves have not changed", { rownames(new_curves) <- NULL reference_curves <- data.frame( 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), + survival_tox_LL5 = c(90.00000, 89.88245, 86.18591, 40.42560, 4.30687), + survival_tox_env_LL5 = c(29.66667, 29.66667, 29.65530, 9.27616, 0.78712), + survival_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), + survival_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), + survival_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), + survival_tox_env_sys = c(29.35449, 24.84147, 32.75385, 1.93182, 0), concentration_for_plots = c(0.0001, 0.03004, 0.30512, 3.02551, 30) ) class(new_curves) <- class(reference_curves) @@ -180,9 +180,9 @@ test_that("function arguments are returned unchanged", { args_reference <- list( 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, + survival_tox_observed = c(90, 81, 92, 28, 0), + survival_tox_env_observed = c(29, 27, 33, 5, 0), + survival_max = 100, curves_concentration_max = NULL, p = 3.2, q = 3.2 @@ -194,45 +194,45 @@ 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) + survival_tox_observed = c(90, 81, 92, 28, 0), + survival_tox_env_observed = c(29, 27, 33, 5, 0) ) test_that("results are independent of concentration shift", { mod_2 <- ecxsys( 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) + survival_tox_observed = c(90, 81, 92, 28, 0), + survival_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) + expect_equal(mod$survival_tox_sys, mod_2$survival_tox_sys) + expect_equal(mod$survival_tox_env_sys, mod_2$survival_tox_env_sys) mod_10 <- ecxsys( 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) + survival_tox_observed = c(90, 81, 92, 28, 0), + survival_tox_env_observed = c(29, 27, 33, 5, 0) ) # 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) - expect_equal(mod$curves$effect_tox_sys, - mod_10$curves$effect_tox_sys) - expect_equal(mod$curves$effect_tox_env_sys, - mod_10$curves$effect_tox_env_sys) + expect_equal(mod$curves$survival_tox, + mod_10$curves$survival_tox) + expect_equal(mod$curves$survival_tox_sys, + mod_10$curves$survival_tox_sys) + expect_equal(mod$curves$survival_tox_env_sys, + mod_10$curves$survival_tox_env_sys) }) -test_that("effect_tox_env_observed can be left out", { +test_that("survival_tox_env_observed can be left out", { mod_without_env <- ecxsys( concentration = c(0, 0.05, 0.5, 5, 30), hormesis_concentration = 0.5, - effect_tox_observed = c(90, 81, 92, 28, 0) + survival_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, - mod_without_env$curves$effect_tox_sys) + expect_equal(mod$survival_tox_sys, mod_without_env$survival_tox_sys) + expect_equal(mod$curves$survival_tox_sys, + mod_without_env$curves$survival_tox_sys) }) @@ -241,7 +241,7 @@ test_that("sys model not converging produces a warning, not an error", { ecxsys( concentration = c(0, 0.1, 0.5, 1, 10, 33), hormesis_concentration = 0.5, - effect_tox_observed = c(98, 98, 96, 76, 26, 0) + survival_tox_observed = c(98, 98, 96, 76, 26, 0) ), paste( "Using a horizontal linear model for sys_tox_mod because the", @@ -253,8 +253,8 @@ test_that("sys model not converging produces a warning, not an error", { ecxsys( concentration = c(0, 0.1, 0.5, 1, 10, 33), hormesis_concentration = 0.5, - effect_tox_observed = c(90, 81, 92, 50, 18, 0), - effect_tox_env_observed = c(75, 89, 54, 7, 0, 0) + survival_tox_observed = c(90, 81, 92, 50, 18, 0), + survival_tox_env_observed = c(75, 89, 54, 7, 0, 0) ), paste( "Using a horizontal linear model for sys_tox_env_mod because the", @@ -270,7 +270,7 @@ test_that("error when curves_concentration_max is too low", { ecxsys( concentration = c(0, 0.05, 0.5, 5, 30), hormesis_concentration = 0.5, - effect_tox_observed = c(90, 81, 92, 28, 0), + survival_tox_observed = c(90, 81, 92, 28, 0), curves_concentration_max = 0.01 ), "'curves_concentration_max' is too low.", diff --git a/tests/testthat/test-ec.R b/tests/testthat/test-lc.R similarity index 58% rename from tests/testthat/test-ec.R rename to tests/testthat/test-lc.R index fe5bbfc..b77738e 100644 --- a/tests/testthat/test-ec.R +++ b/tests/testthat/test-lc.R @@ -21,27 +21,27 @@ model <- 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) + survival_tox_observed = c(90, 81, 92, 28, 0), + survival_tox_env_observed = c(29, 27, 33, 5, 0) ) test_that("all input formats produce identical models", { # 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) + ec10_a <- lc(model, "survival_tox_sys", 10) + ec10_b <- lc(model$curves, "survival_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) + ec10_c <- lc(curves, "survival_tox_sys", 10) # using a custom data frame: df_custom <- data.frame( concentration = curves$concentration, - foo = curves$effect_tox_sys + foo = curves$survival_tox_sys ) - ec10_d <- ec(df_custom, "foo", 10) + ec10_d <- lc(df_custom, "foo", 10) expect_equal(ec10_a, ec10_b, tolerance = 1e-5) expect_equal(ec10_b, ec10_c, tolerance = 1e-5) @@ -49,58 +49,58 @@ test_that("all input formats produce identical models", { }) -test_that("ec values have not changed", { +test_that("lc values have not changed", { expect_equal( - ec(model, "effect_tox_sys", 50), - list(effect = 44.95368, concentration = 3.375735), + lc(model, "survival_tox_sys", 50), + list(response = 44.95368, concentration = 3.375735), tolerance = 1e-4 ) expect_equal( - ec(model, "effect_tox_sys", 10), - list(effect = 80.91662, concentration = 1.098648), + lc(model, "survival_tox_sys", 10), + list(response = 80.91662, concentration = 1.098648), tolerance = 1e-4 ) expect_equal( - ec(model, "effect_tox", 100/3), - list(effect = 66.66667, concentration = 1.902125), + lc(model, "survival_tox", 100/3), + list(response = 66.66667, concentration = 1.902125), tolerance = 1e-4 ) expect_equal( - ec(model, "effect_tox_LL5", 42), - list(effect = 52.2, concentration = 2.054426), + lc(model, "survival_tox_LL5", 42), + list(response = 52.2, concentration = 2.054426), tolerance = 1e-4 ) expect_equal( - ec(model, "effect_tox_env_sys", 50), - list(effect = 14.67725, concentration = 1.299516), + lc(model, "survival_tox_env_sys", 50), + list(response = 14.67725, concentration = 1.299516), tolerance = 1e-4 ) expect_equal( - ec(model, "effect_tox_env_sys", 10), - list(effect = 26.41904, concentration = 0.0008571244), + lc(model, "survival_tox_env_sys", 10), + list(response = 26.41904, concentration = 0.0008571244), tolerance = 1e-4 ) expect_equal( - ec(model, "effect_tox_env", 67.89), - list(effect = 24.51453, concentration = 0.7890892), + lc(model, "survival_tox_env", 67.89), + list(response = 24.51453, concentration = 0.7890892), tolerance = 1e-4 ) expect_equal( - ec(model, "effect_tox_env_LL5", 3.14159), - list(effect = 28.73466, concentration = 0.7267524), + lc(model, "survival_tox_env_LL5", 3.14159), + list(response = 28.73466, concentration = 0.7267524), tolerance = 1e-4 ) }) -test_that("ec warns when response_level is outside the range of the curve", { +test_that("lc warns when response_level is outside the range of the curve", { model <- ecxsys( concentration = c(0, 0.05, 0.5, 5), hormesis_concentration = 0.5, - effect_tox_observed = c(90, 81, 92, 28) + survival_tox_observed = c(90, 81, 92, 28) ) expect_warning( - ec(model, "effect_tox_sys", 90), + lc(model, "survival_tox_sys", 90), "You could try using predict_ecxsys() to predict more values", fixed = TRUE ) @@ -109,13 +109,13 @@ test_that("ec warns when response_level is outside the range of the curve", { test_that("reference argument works", { expect_equal( - ec(model, "effect_tox_LL5", 50, reference = 100), - list(effect = 50, concentration = 2.208119), + lc(model, "survival_tox_LL5", 50, reference = 100), + list(response = 50, concentration = 2.208119), tolerance = 1e-4 ) expect_equal( - ec(model, "effect_tox_LL5", 50, reference = 75), - list(effect = 37.5, concentration = 3.342715), + lc(model, "survival_tox_LL5", 50, reference = 75), + list(response = 37.5, concentration = 3.342715), tolerance = 1e-4 ) }) diff --git a/tests/testthat/test-plot_ecxsys.R b/tests/testthat/test-log10_ticks.R similarity index 50% rename from tests/testthat/test-plot_ecxsys.R rename to tests/testthat/test-log10_ticks.R index a0dea0a..3f6f93d 100644 --- a/tests/testthat/test-plot_ecxsys.R +++ b/tests/testthat/test-log10_ticks.R @@ -18,13 +18,25 @@ test_that("log ticks are correct", { - x <- c(0.03, 0.3, 3, 30) - ticks <- get_log_ticks(x) - expect_equal(ticks$major, c(0.01, 0.10, 1.00, 10.00, 100.00)) - expect_equal(ticks$major_labels, c("0", "0.1", "1", "10", "100")) - expect_equal( - ticks$minor, - c(0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.20, 0.30, 0.40, 0.50, - 0.60, 0.70, 0.80, 0.90, 2.00, 3.00, 4.00, 5.00, 6.00, 7.00, 8.00, 9.00, - 20.00, 30.00, 40.00, 50.00, 60.00, 70.00, 80.00, 90.00)) + x <- c(0.03, 0.3, 3, 30) + reference <- list( + major = c(0.01, 0.1, 1, 10, 100), + minor = c( + 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, + 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, + 2, 3, 4, 5, 6, 7, 8, 9, + 20, 30, 40, 50, 60, 70, 80, 90 + ), + major_labels = c("0", "0.1", "1", "10", "100"), + minor_labels = c( + "0.02", "0.03", "0.04", "0.05", "0.06", "0.07", "0.08", "0.09", + "0.2", "0.3", "0.4", "0.5", "0.6", "0.7", "0.8", "0.9", + "2", "3", "4", "5", "6", "7", "8", "9", + "20", "30", "40", "50", "60", "70", "80", "90" + ) + ) + expect_equal(reference, log10_ticks(x)) + + reference$major_labels[1] <- "0.01" + expect_equal(reference, log10_ticks(x, label_zero = FALSE)) }) diff --git a/tests/testthat/test-multi_tox.R b/tests/testthat/test-multi_tox.R new file mode 100644 index 0000000..0ac7256 --- /dev/null +++ b/tests/testthat/test-multi_tox.R @@ -0,0 +1,117 @@ +# Copyright (C) 2020 Helmholtz-Zentrum fuer Umweltforschung GmbH - UFZ +# See file inst/COPYRIGHTS for details. +# +# This file is part of the R package stressaddition. +# +# stressaddition is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see <https://www.gnu.org/licenses/>. + + +model_a <- ecxsys( + concentration = c(0, 0.05, 0.5, 5, 30), + hormesis_concentration = 0.5, + survival_tox_observed = c(90, 81, 92, 28, 0), + survival_tox_env_observed = c(29, 27, 33, 5, 0), + survival_max = 101 +) +model_b <- ecxsys( + concentration = c(0, 0.1, 1, 10, 100, 1000), + hormesis_concentration = 10, + survival_tox_observed = c(26, 25, 24, 27, 5, 0), + survival_max = 30 +) + + +test_that("results have not changed", { + # one concentration_b + new <- multi_tox( + model_a, + model_b, + c(0, 0.01, 0.1, 1, 7, 15), + 5, + 0.7 + ) + reference <- data.frame( + concentration_a = c(0, 0.01, 0.1, 1, 7, 15), + concentration_b = 5, + survival = c(88.574578, 84.361552, 80.633762, 56.730550, 2.882718, 0), + stress_tox_sa = c(0.1886787, 0.2436883, 0.3185115, 0.5111081, 0.8904295, 1.0700716), + stress_tox_ca = c(0.18867873, 0.19199327, 0.21634282, 3.579444e-01, 7.126441e-01, 8.856914e-01), + stress_tox = c(0.18867873, 0.22817978, 0.28786093, 4.651590e-01, 8.370939e-01, 1.014758), + sys = c(0.07845536, 0.07312129, 0.04004696, 4.663355e-05, 0, 0), + stress_total = c(0.26713409, 0.30130107, 0.32790789, 4.652056e-01, 8.370939e-01, 1.014758) + ) + expect_equal(new, reference, tolerance = 1e-5) + + # diverse concentration_b + new <- multi_tox( + model_a, + model_b, + c(0, 0.01, 0.1, 1, 7, 15), + c(0, 0.02, 0.2, 2, 14, 30), + 0.7 + ) + reference <- data.frame( + concentration_a = c(0, 0.01, 0.1, 1, 7, 15), + concentration_b = c(0, 0.02, 0.2, 2, 14, 30), + survival = c(88.2698383, 79.96171270, 78.1574808, 6.579998e+01, 3.861678e-01, 0), + stress_tox_sa = c(0, 0.07567745, 0.1807501, 4.510696e-01, 9.961584e-01, 1.294227), + stress_tox_ca = c(0, 0.05654834, 0.1343429, 3.367461e-01, 7.321990e-01, 9.063311e-01), + stress_tox = c(0, 0.06993872, 0.1668279, 4.167725e-01, 9.169706e-01, 1.177858), + sys = c(0.2698033, 0.26248944, 0.1774490, 1.829786e-04, 0, 0), + stress_total = c(0.2698033, 0.33242815, 0.3442769, 4.169555e-01, 9.169706e-01, 1.177858) + ) + expect_equal(new, reference, tolerance = 1e-5) + + # diverse concentration_b and custom survival_max + new <- multi_tox( + model_a, + model_b, + c(0, 0.01, 0.1, 1, 7, 15), + c(0, 0.02, 0.2, 2, 14, 30), + 0.7, + 42 + ) + reference <- data.frame( + concentration_a = c(0, 0.01, 0.1, 1, 7, 15), + concentration_b = c(0, 0.02, 0.2, 2, 14, 30), + survival = c(37.0733321, 33.58391933, 32.8261419, 2.763599e+01, 1.621905e-01, 0), + stress_tox_sa = c(0, 0.07567745, 0.1807501, 4.510696e-01, 9.961584e-01, 1.294227), + stress_tox_ca = c(0, 0.05654834, 0.1343429, 3.367461e-01, 7.321990e-01, 9.063311e-01), + stress_tox = c(0, 0.06993872, 0.1668279, 4.167725e-01, 9.169706e-01, 1.177858), + sys = c(0.2698033, 0.26248944, 0.1774490, 1.829786e-04, 0, 0), + stress_total = c(0.2698033, 0.33242815, 0.3442769, 4.169555e-01, 9.169706e-01, 1.177858) + ) + expect_equal(new, reference, tolerance = 1e-5) +}) + + +test_that("predictions are symmetric", { + conc_a <- c(0, 10^seq(log10(0.001), log10(40), length.out = 50)) + conc_b <- 3.5 + sa_contrib <- 0.8 + survival_12 <- multi_tox(model_a, model_b, conc_a, conc_b, sa_contrib)$survival + survival_21 <- multi_tox(model_b, model_a, conc_b, conc_a, sa_contrib)$survival + expect_equal(survival_12, survival_21) +}) + + +test_that("survival_tox_mod is a W1.2 model", { + # This is a safeguard for if I decide to use something else than drc::W1.2 + # for survival_tox_mod. It is extremely unlikely that this ever happens but + # better safe than sorry. Currently, the inverse only works for W1.2 models. + # If this test throws errors then I probably forget to adjust the inverse + # function accordingly. + expect_s3_class(model_a$survival_tox_mod$fct, "Weibull-1") + expect_equal(model_a$survival_tox_mod$fct$noParm, 2) +}) diff --git a/tests/testthat/test-options.R b/tests/testthat/test-options.R index 2ff4a81..2a4f03c 100644 --- a/tests/testthat/test-options.R +++ b/tests/testthat/test-options.R @@ -26,13 +26,14 @@ test_that("user options are not permanently changed by ecxsys()", { # This problem becomes visible only if the user changes the default options: options(show.error.messages = FALSE) # default is TRUE options(warn = 1) # default is 0 + on.exit(options(warn = 0)) original_options <- options() model <- 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) + survival_tox_observed = c(90, 81, 92, 28, 0), + survival_tox_env_observed = c(29, 27, 33, 5, 0) ) new_options <- options() diff --git a/tests/testthat/test-predict_mixture.R b/tests/testthat/test-predict_mixture.R deleted file mode 100644 index a1974bb..0000000 --- a/tests/testthat/test-predict_mixture.R +++ /dev/null @@ -1,90 +0,0 @@ -# Copyright (C) 2020 Helmholtz-Zentrum fuer Umweltforschung GmbH - UFZ -# See file inst/COPYRIGHTS for details. -# -# This file is part of the R package stressaddition. -# -# stressaddition is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see <https://www.gnu.org/licenses/>. - - -model_a <- 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_b <- 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", { - # one concentration_b - new <- predict_mixture( - model_a, - model_b, - c(0, 0.01, 0.1, 1, 7, 15), - 5, - 0.3 - )$effect - reference <- c(88.574578, 84.361552, 80.633762, 56.730550, 2.882718, 0) - expect_equal(new, reference, tolerance = 1e-5) - - # diverse concentration_b - new <- predict_mixture( - model_a, - model_b, - c(0, 0.01, 0.1, 1, 7, 15), - c(0, 0.02, 0.2, 2, 14, 30), - 0.3 - )$effect - reference <- c(88.2698383, 79.9617127, 78.1574808, 65.7999834, 0.3861678, 0) - expect_equal(new, reference, tolerance = 1e-5) - - # diverse concentration_b and custom effect_max - new <- predict_mixture( - model_a, - model_b, - c(0, 0.01, 0.1, 1, 7, 15), - c(0, 0.02, 0.2, 2, 14, 30), - 0.3, - 42 - )$effect - reference <- c(88.2698383, 79.9617127, 78.1574808, 65.7999834, 0.3861678, 0) * 0.42 - expect_equal(new, reference, tolerance = 1e-5) -}) - - -test_that("predictions are symmetric", { - conc_a <- c(0, 10^seq(log10(0.001), log10(40), length.out = 50)) - conc_b <- 3.5 - prop_ca <- 0.8 - effect_12 <- predict_mixture(model_a, model_b, conc_a, conc_b, prop_ca)$effect - effect_21 <- predict_mixture(model_b, model_a, conc_b, conc_a, prop_ca)$effect - expect_equal(effect_12, effect_21) -}) - - -test_that("effect_tox_mod is a W1.2 model", { - # This is a safeguard for if I decide to use something else than drc::W1.2 - # for effect_tox_mod. It is extremely unlikely that this ever happens but - # better safe than sorry. Currently, the inverse only works for W1.2 models. - # If this test throws errors then I probably forget to adjust the inverse - # function accordingly. - expect_s3_class(model_a$effect_tox_mod$fct, "Weibull-1") - expect_equal(model_a$effect_tox_mod$fct$noParm, 2) -}) diff --git a/tests/testthat/test-stress_effect_conversion.R b/tests/testthat/test-stress_effect_conversion.R deleted file mode 100644 index e17a2ae..0000000 --- a/tests/testthat/test-stress_effect_conversion.R +++ /dev/null @@ -1,44 +0,0 @@ -# Copyright (C) 2020 Helmholtz-Zentrum fuer Umweltforschung GmbH - UFZ -# See file inst/COPYRIGHTS for details. -# -# This file is part of the R package stressaddition. -# -# stressaddition is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see <https://www.gnu.org/licenses/>. - - -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_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_error(effect_to_stress(0.5, -3, 10)) - expect_error(effect_to_stress(0.5, 3, -10)) -}) - -test_that("stress_to_effect is correct", { - expect_equal(stress_to_effect(0.5), 0.5) - expect_equal(round(stress_to_effect(0.1), 5), 0.99321) -}) - -test_that("effect_to_stress is correct", { - expect_equal(effect_to_stress(0.5), 0.5) - expect_equal(round(effect_to_stress(0.1), 5), 0.74588) -}) diff --git a/tests/testthat/test-stress_survival_conversion.R b/tests/testthat/test-stress_survival_conversion.R new file mode 100644 index 0000000..76aef7e --- /dev/null +++ b/tests/testthat/test-stress_survival_conversion.R @@ -0,0 +1,44 @@ +# Copyright (C) 2020 Helmholtz-Zentrum fuer Umweltforschung GmbH - UFZ +# See file inst/COPYRIGHTS for details. +# +# This file is part of the R package stressaddition. +# +# stressaddition is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see <https://www.gnu.org/licenses/>. + + +test_that("stress_to_survival handles bad input", { + expect_equal(stress_to_survival(-3), 1) + expect_equal(stress_to_survival(10), 0) + expect_error(stress_to_survival("a")) + expect_error(stress_to_survival(0.5, -3, 10)) + expect_error(stress_to_survival(0.5, 3, -10)) +}) + +test_that("survival_to_stress handles bad input", { + expect_equal(survival_to_stress(-15), 1) + expect_equal(survival_to_stress(20), 0) + expect_error(survival_to_stress("a")) + expect_error(survival_to_stress(0.5, -3, 10)) + expect_error(survival_to_stress(0.5, 3, -10)) +}) + +test_that("stress_to_survival is correct", { + expect_equal(stress_to_survival(0.5), 0.5) + expect_equal(round(stress_to_survival(0.1), 5), 0.99321) +}) + +test_that("survival_to_stress is correct", { + expect_equal(survival_to_stress(0.5), 0.5) + expect_equal(round(survival_to_stress(0.1), 5), 0.74588) +}) -- GitLab