Merge branch 'develop' into 'master'

curves_concentration_max and ec error

Closes #25 and #29

See merge request !21

.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

DESCRIPTION

 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: person("Sebastian", "Henz", role = c("aut", "cre"), 
                   email = "sebastian.henz@ufz.de", 
                   comment = c(ORCID = "0000-0001-8299-8852"))

NEWS.md

+# 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`.

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

R/ecxsys.R

 # 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
 }
 

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:

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

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{

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
+    )
+})

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
+    )
+})