From 6e2fa522fd71cf06edea512519d5acaeb743cea5 Mon Sep 17 00:00:00 2001
From: Sebastian Henz <sebastian.henz@ufz.de>
Date: Thu, 27 Feb 2020 13:41:19 +0100
Subject: [PATCH] * `ecxsys()` gains a new argument `curves_concentration_max`
which sets the maximum concentration for the curves. * `ec()` raises an error
if the requested curve is too short to calculate the result. * Improve merge
request template. * Update readme example.
---
.../merge_into_master.md | 9 +++---
DESCRIPTION | 4 +--
NEWS.md | 5 +++
R/ec.R | 26 +++++++++++----
R/ecxsys.R | 32 ++++++++++++-------
README.md | 8 ++---
man/ec.Rd | 7 +++-
man/ecxsys.Rd | 5 +++
tests/testthat/test-ec.R | 14 ++++++++
tests/testthat/test-ecxsys.R | 19 +++++++++--
10 files changed, 97 insertions(+), 32 deletions(-)
diff --git a/.gitlab/merge_request_templates/merge_into_master.md b/.gitlab/merge_request_templates/merge_into_master.md
index 670d525..a8b3a26 100644
--- a/.gitlab/merge_request_templates/merge_into_master.md
+++ b/.gitlab/merge_request_templates/merge_into_master.md
@@ -18,14 +18,13 @@ Comment on any notes that come up during devtools::check().
### checklist
<!-- Please check `[x]` the applicable boxes. -->
-* [ ] This merge includes bug fixes or changes in the functionality of the package (whether they are user-visible or not).
- * [ ] update documentation including examples (if necessary)
- * [ ] update and polish NEWS.md
- * [ ] update version in DESCRIPTION
- * [ ] update date in DESCRIPTION
* [ ] This merge introduces new features.
* [ ] new features are documented
* [ ] new features have tests
+* [ ] update documentation including examples (if necessary)
+* [ ] update and polish NEWS.md
+* [ ] update version in DESCRIPTION
+* [ ] update date in DESCRIPTION
* [ ] no remaining TODO, FIXME, or debug prints anywhere in the source files
* [ ] `devtools::document()`
* [ ] `devtools::check()` without errors or warnings
diff --git a/DESCRIPTION b/DESCRIPTION
index e96c337..ae78b57 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
Package: stressaddition
Type: Package
Title: Modeling Tri-Phasic Concentration-Response Relationships
-Version: 2.1.1
-Date: 2020-02-26
+Version: 2.2.0
+Date: 2020-02-27
Authors@R: c(person("Sebastian",
"Henz",
role = c("aut", "cre"),
diff --git a/NEWS.md b/NEWS.md
index 53a54a4..1379a56 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,8 @@
+# stressaddition 2.2.0
+
+* `ec()` raises an error if the curve does not cross the desired response level.
+* `ecxsys()` gains a new argument `curves_concentration_max` which allows setting the maximum concentration of the predicted curves.
+
# stressaddition 2.1.1
* Restore the default behaviour of `plot_effect()` to also show `effect_tox` and `effect_tox_env`.
diff --git a/R/ec.R b/R/ec.R
index 66d2592..5fb444b 100644
--- a/R/ec.R
+++ b/R/ec.R
@@ -7,6 +7,11 @@
#' 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
+#' \code{sys_tox} and \code{sys_tox_env}. Others are untested and may give
+#' unexpected results, if any.
+#'
#' @param model This can be one of three types of objects: Either the output of
#' \code{\link{ecxsys}} or the output of \code{\link{predict_ecxsys}} or a
#' data frame with a "concentration" column and a \code{response_name} column.
@@ -15,7 +20,7 @@
#' calculate the EC. Must be one of \code{colnames(model$curves)}.
#' @param response_level The desired response level as a percentage between 0
#' and 100. For example with the value 10 the function will return the EC10,
-#' i.e. the concentration where the response falls below 90 \% of the maximum
+#' i.e. the concentration where the response falls below 90 \% of the control
#' response.
#'
#' @return A list containing the response concentration and the corresponding
@@ -77,18 +82,25 @@ ec <- function(model, response_name, response_level) {
reference <- response[1]
if (reference == 0) {
- stop("Reference value is zero, calculation of EC not possible.")
+ stop("Reference value is zero, calculation of EC is impossible.")
}
- response_level <- (1 - response_level / 100) * reference
- output <- list(response_value = response_level)
+ response_value <- (1 - response_level / 100) * reference
+ output <- list(response_value = response_value)
# Get the index of where the response changes from above to below
- # response_level:
- below <- which(response < response_level)[1]
+ # response_value:
+ below <- which(response < response_value)[1]
+ if (is.na(below)) {
+ stop("The curve '", response_name, "' does not fall below ",
+ 100 - response_level, "% of the control, which makes ",
+ "determining the EC", response_level, " impossible.\n You could ",
+ "try using predict_ecxsys() to predict more values in a wider ",
+ "concentration range.")
+ }
above <- below - 1
# linear interpolation
- dist <- (response_level - response[below]) / (response[above] - response[below])
+ dist <- (response_value - response[below]) / (response[above] - response[below])
output$concentration <- dist *
(concentration[above] - concentration[below]) + concentration[below]
output
diff --git a/R/ecxsys.R b/R/ecxsys.R
index 95e090d..ccab533 100644
--- a/R/ecxsys.R
+++ b/R/ecxsys.R
@@ -1,9 +1,10 @@
# Note about the setting and resetting of options:
# 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 ecxsys function can generate some warnings. So every time drc::drm is
-# called this option is cached beforehand and reset afterwards. Additionally,
-# all options are cached at the beginning of ecxsys and reset on exit.
+# 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
+# beforehand and reset immediately afterwards. Additionally, all options are
+# cached at the beginning of ecxsys and reset on exit.
#' ECx-SyS
@@ -31,6 +32,9 @@
#' environmental stress. Must be between 0 and \code{effect_max}.
#' @param effect_max The maximum value the effect could possibly reach. For
#' survival data in percent this should be 100 (the default).
+#' @param 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.
#'
#' @return A list (of class ecxsys) containing many different objects of which
@@ -86,6 +90,7 @@ ecxsys <- function(concentration,
effect_tox_observed,
effect_tox_env_observed = NULL,
effect_max = 100,
+ curves_concentration_max = NULL,
p = 3.2,
q = 3.2) {
output <- list(args = as.list(environment()))
@@ -268,23 +273,28 @@ ecxsys <- function(concentration,
}
- # smooth curves -------------------------------------------------------
+ # curves -------------------------------------------------------
# In order to generate a broken x-axis the concentration vector must
# also be broken in two. The left part of the axis is supposed to be at
# 0 but because it's a log axis I have to make the values just really
# small. The concentrations in the gap won't be used for plotting later.
- n_smooth <- 1000 # number of points to approximate the curves
+ len_curves <- 1000 # number of points to approximate the curves
conc_adjust_factor <- 10^-5
output$conc_adjust_factor <- conc_adjust_factor
- concentration_smooth <- 10 ^ seq(
+ if (is.null(curves_concentration_max)) {
+ curves_concentration_max <- max(concentration)
+ } else if (curves_concentration_max < min(concentration[concentration > 0])) {
+ stop("'curves_concentration_max' is too low.")
+ }
+ curves_concentration <- 10 ^ seq(
log10(min_conc * conc_adjust_factor),
- log10(max(concentration)),
- length.out = n_smooth
+ log10(curves_concentration_max),
+ length.out = len_curves
)
- output$curves <- predict_ecxsys(output, concentration_smooth)
+ output$curves <- predict_ecxsys(output, curves_concentration)
output$curves$use_for_plotting <-
- concentration_smooth < min_conc * conc_adjust_factor * 1.5 |
- concentration_smooth > min_conc * 1.5
+ curves_concentration < min_conc * conc_adjust_factor * 1.5 |
+ curves_concentration > min_conc * 1.5
output
}
diff --git a/README.md b/README.md
index 6522297..7ffe064 100644
--- a/README.md
+++ b/README.md
@@ -26,10 +26,10 @@ In the paper we use the model in the context of survival experiments. However, i
```r
library(stressaddition)
model <- ecxsys(
- concentration = c(0, 0.03, 0.3, 3, 10),
- hormesis_concentration = 0.3,
- effect_tox_observed = c(85, 76, 94, 35, 0),
- effect_tox_env_observed = c(24, 23, 32, 0, 0)
+ concentration = c(0, 0.05, 0.5, 5, 30),
+ hormesis_concentration = 0.5,
+ effect_tox_observed = c(90, 81, 92, 28, 0),
+ effect_tox_env_observed = c(29, 27, 33, 5, 0)
)
# Plot the effect and the system stress:
diff --git a/man/ec.Rd b/man/ec.Rd
index dea6581..0cacd6a 100644
--- a/man/ec.Rd
+++ b/man/ec.Rd
@@ -17,7 +17,7 @@ calculate the EC. Must be one of \code{colnames(model$curves)}.}
\item{response_level}{The desired response level as a percentage between 0
and 100. For example with the value 10 the function will return the EC10,
-i.e. the concentration where the response falls below 90 \% of the maximum
+i.e. the concentration where the response falls below 90 \% of the control
response.}
}
\value{
@@ -32,6 +32,11 @@ to the control.
If the response level occurs multiple times because of hormesis, which may
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
+\code{sys_tox} and \code{sys_tox_env}. Others are untested and may give
+unexpected results, if any.
}
\examples{
# Calculate the EC_10, the concentration where the effect falls
diff --git a/man/ecxsys.Rd b/man/ecxsys.Rd
index 620adb3..4d5bed1 100644
--- a/man/ecxsys.Rd
+++ b/man/ecxsys.Rd
@@ -10,6 +10,7 @@ ecxsys(
effect_tox_observed,
effect_tox_env_observed = NULL,
effect_max = 100,
+ curves_concentration_max = NULL,
p = 3.2,
q = 3.2
)
@@ -31,6 +32,10 @@ environmental stress. Must be between 0 and \code{effect_max}.}
\item{effect_max}{The maximum value the effect could possibly reach. For
survival data in percent this should be 100 (the default).}
+\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.}
+
\item{p, q}{The shape parameters of the beta distribution. Default is 3.2.}
}
\value{
diff --git a/tests/testthat/test-ec.R b/tests/testthat/test-ec.R
index a4cb740..138fbce 100644
--- a/tests/testthat/test-ec.R
+++ b/tests/testthat/test-ec.R
@@ -74,3 +74,17 @@ test_that("ec values have not changed", {
tolerance = 1e-4
)
})
+
+
+test_that("ec fails 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)
+ )
+ expect_error(
+ ec(model, "effect_tox_sys", 90),
+ "You could try using predict_ecxsys() to predict more values",
+ fixed = TRUE
+ )
+})
diff --git a/tests/testthat/test-ecxsys.R b/tests/testthat/test-ecxsys.R
index 0120163..9911f1d 100644
--- a/tests/testthat/test-ecxsys.R
+++ b/tests/testthat/test-ecxsys.R
@@ -164,6 +164,7 @@ test_that("function arguments are returned unchanged", {
effect_tox_observed = c(90, 81, 92, 28, 0),
effect_tox_env_observed = c(29, 27, 33, 5, 0),
effect_max = 100,
+ curves_concentration_max = NULL,
p = 3.2,
q = 3.2
)
@@ -232,9 +233,9 @@ test_that("sys model not converging produces a warning, not an error", {
expect_warning(
ecxsys(
concentration = c(0, 0.1, 0.5, 1, 10, 33),
+ 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),
- hormesis_concentration = 0.5
+ effect_tox_env_observed = c(75, 89, 54, 7, 0, 0)
),
paste(
"Using a horizontal linear model for sys_tox_env_mod because the",
@@ -243,3 +244,17 @@ test_that("sys model not converging produces a warning, not an error", {
fixed = TRUE
)
})
+
+
+test_that("error when curves_concentration_max is too low", {
+ expect_error(
+ ecxsys(
+ concentration = c(0, 0.05, 0.5, 5, 30),
+ hormesis_concentration = 0.5,
+ effect_tox_observed = c(90, 81, 92, 28, 0),
+ curves_concentration_max = 0.01
+ ),
+ "'curves_concentration_max' is too low.",
+ fixed = TRUE
+ )
+})
--
GitLab