Improve ec() and axis concentration handling

description

• The curves data frame in the output of ecxsys() now contains a column with the concentrations which are used for the plot functions in this package. This is useful for generating a nicer concentration axis.
• Changes to ec():
  • Renamed response_value to effect in the output list.
  • response_level of 0 or 100 is now allowed. 0 returns the concentration 0 and 100 returns the concentration Inf. Previously this resulted in an error.
  • It is now possible to set the reference to a custom value, for example 100.

system details

• RStudio Version: 1.3.937
• R Version: 3.6.3
• OS Version: Windows 10 x64 (build 18363)

checklist

• 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

R/ec.R

@@ -34,6 +34,7 @@
 #' @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.
+#'   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)}.
@@ -41,12 +42,21 @@
 #'   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 control
 #'   response.
+#' @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
+#'   \code{response_level} actually corresponds to. For example if
+#'   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?
 #'
-#' @return A list containing the response concentration and the corresponding
-#'   response value.
+#' @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.
 #'
-#' @examples # Calculate the EC_10, the concentration where the effect falls
-#' # below 90 % of the effect in the control.
+#' @examples # Calculate the EC10, the concentration where the effect falls
+#' # below 90% of the effect in the control.
 #'
 #' model <- ecxsys(
 #'     concentration = c(0, 0.05, 0.5, 5, 30),
@@ -70,8 +80,16 @@
 #' )
 #' ec(df_custom, "foo", 10)
 #'
+#' # Calculate the EC50 relative to an effect of 100
+#' # instead of relative to the control:
+#' ec(model, "effect_tox_sys", 50, reference = 100)
+#'
 #' @export
-ec <- function(model, response_name, response_level) {
+ec <- function(model,
+               response_name,
+               response_level,
+               reference,
+               warn = TRUE) {
     if (inherits(model, "drc")) {
         stop("Please use drc::ED for drc objects.")
     }
@@ -80,8 +98,8 @@ ec <- function(model, response_name, response_level) {
         is.character(response_name),
         length(response_name) == 1,
         response_name != "concentration",
-        response_level > 0,
-        response_level < 100
+        response_level >= 0,
+        response_level <= 100
     )
 
     if (!inherits(model, c("ecxsys", "ecxsys_predicted"))) {
@@ -99,28 +117,73 @@ ec <- function(model, response_name, response_level) {
         response <- model[, response_name]
     }
 
-    reference <- response[1]
-    if (reference == 0) {
-        stop("Reference value is zero, calculation of EC is impossible.")
+    if (missing(reference)) {
+        reference_value <- response[1]
+    } else {
+        stopifnot(reference > 0)
+        if (inherits(model, "ecxsys")) {
+            stopifnot(reference <= model$args$effect_max)
+        }
+        reference_value <- reference
     }
-    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_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.")
+    if (reference_value == 0) {
+        # May happen with horizontal Sys curves.
+        if (warn) {
+            warning("Reference value is zero, calculation of EC is impossible.")
+        }
+        return(list(effect = 0, concentration = NA))
+    }
+
+    if (response_level == 0) {
+        return(list(effect = reference_value, concentration = 0))
+    } else if (response_level == 100) {
+        return(list(effect = 0, concentration = Inf))
+    }
+
+    response_value <- (1 - response_level / 100) * reference_value
+
+    response_smaller <- response < response_value
+    if (response_smaller[1]) {
+        if (warn) {
+            warning(
+                "The first values of '",
+                response_name,
+                "' are smaller than ",
+                100 - response_level,
+                "% of the reference value, which makes determining the EC",
+                response_level,
+                " impossible."
+            )
+        }
+        return(list(effect = response_value, concentration = NA))
+    }
+    if (!any(response_smaller)) {
+        if (warn) {
+            warning(
+                "The curve '",
+                response_name,
+                "' does not fall below ",
+                100 - response_level,
+                "% of the reference, which makes determining the EC",
+                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))
     }
+    if (response[1] == response_value) {
+        return(list(effect = response_value, concentration = 0))
+    }
+
+    below <- which(response < response_value)[1]
     above <- below - 1
 
     # linear interpolation
     dist <- (response_value - response[below]) / (response[above] - response[below])
-    output$concentration <- dist *
+    effect_concentration <- dist *
         (concentration[above] - concentration[below]) + concentration[below]
-    output
+
+    list(effect = response_value, concentration = effect_concentration)
 }