From 1aae72c3b2aa64bf5b70978ef1474bff72f61f93 Mon Sep 17 00:00:00 2001
From: Sebastian Henz <sebastian.henz@ufz.de>
Date: Tue, 7 Apr 2020 15:08:26 +0200
Subject: [PATCH] Version 2.6.0: Improve ec() and axis concentration handling
 Replace the "use_for_plotting" column in mod$curves with the plot axis
 concentrations, which is far more useful. Add "refence" argument to ec() and
 some minor improvements.

---
 DESCRIPTION                  |   4 +-
 NEWS.md                      |   8 +++
 R/ec.R                       | 109 +++++++++++++++++++++++++++--------
 R/ecxsys.R                   |  40 +++++++++++--
 R/helpers.R                  |  26 ---------
 R/plot_effect.R              |  30 +++++-----
 R/plot_stress.R              |  32 +++++-----
 man/ec.Rd                    |  26 +++++++--
 man/ecxsys.Rd                |   2 +-
 tests/testthat/test-ec.R     |  62 ++++++++++++--------
 tests/testthat/test-ecxsys.R |   4 +-
 11 files changed, 226 insertions(+), 117 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index f327c51..74f812b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: stressaddition
 Type: Package
 Title: Modeling Tri-Phasic Concentration-Response Relationships
-Version: 2.5.0
-Date: 2020-03-16
+Version: 2.6.0
+Date: 2020-04-07
 Authors@R: c(person("Sebastian", 
                     "Henz", 
                     role = c("aut", "cre"), 
diff --git a/NEWS.md b/NEWS.md
index 1270c2a..3c1cfec 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,11 @@
+# stressaddition 2.6.0
+
+* 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.
+
 # stressaddition 2.5.0
 
 * Fixed unintended behaviour in `plot_effect()` and `plot_stress()` where supplying an empty vector caused the four standard curves to show. Now setting `which` to an empty vector or `NULL` shows just the axes. The default value is NA.
diff --git a/R/ec.R b/R/ec.R
index 69c9da6..9f89821 100644
--- a/R/ec.R
+++ b/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)
 }
diff --git a/R/ecxsys.R b/R/ecxsys.R
index 80e05b4..a257531 100644
--- a/R/ecxsys.R
+++ b/R/ecxsys.R
@@ -77,7 +77,7 @@
 #'     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{use_for_plotting} is used by the
+#'     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.}
 #'   }
@@ -311,9 +311,15 @@ ecxsys <- function(concentration,
         length.out = len_curves
     )
     output$curves <- predict_ecxsys(output, curves_concentration)
-    output$curves$use_for_plotting <-
-        curves_concentration < min_conc * conc_adjust_factor * 1.5 |
-        curves_concentration > min_conc * 1.5
+
+    conc_axis <- adjust_plot_concentrations(
+        curves_concentration,
+        min_conc,
+        conc_adjust_factor
+    )
+    output$curves$concentration_for_plots <- conc_axis$concentration
+    output$axis_break_conc <- conc_axis$axis_break_conc
+
     output
 }
 
@@ -432,3 +438,29 @@ interpolate <- function(x, to_index, n_new, conc = FALSE) {
     }
     append(x, x_new[-c(1, len)], from_index)
 }
+
+
+adjust_plot_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.
+    use_for_plotting <- (
+        concentration < min_conc * conc_adjust_factor * 1.5 |
+        concentration > min_conc * 1.5
+    )
+    gap_idx <- min(which(!use_for_plotting))
+
+    # Add NAs to force breaks in the lines:
+    axis_break_conc <- concentration[use_for_plotting][gap_idx]
+    concentration[!use_for_plotting] <- NA
+
+    # Shift the small concentrations upwards so the plot has a nice x-axis:
+    concentration[1:gap_idx] <- concentration[1:gap_idx] / conc_adjust_factor
+
+    list(
+        concentration = concentration,
+        axis_break_conc = axis_break_conc
+    )
+}
diff --git a/R/helpers.R b/R/helpers.R
index acdfaff..fcc3052 100644
--- a/R/helpers.R
+++ b/R/helpers.R
@@ -50,29 +50,3 @@ get_log_ticks <- function(x) {
         major_labels = major_tick_labels
     )
 }
-
-
-adjust_plot_concentrations <- function(model) {
-    # Helper for the plot functions, not exported.
-    # Deals with the concentrations which are unnecessary for plotting. This
-    # means it removes the concentrations in the gap and increases the
-    # concentrations left of the gap.
-    curves <- model$curves
-    gap_idx <- min(which(!curves$use_for_plotting))
-
-    # Keep only the values to the left and right of the axis break:
-    curves <- curves[curves$use_for_plotting, ]
-
-    # Add NAs to force breaks in the lines:
-    axis_break_conc <- curves$concentration[gap_idx]
-    curves[gap_idx, ] <- NA
-
-    # Shift the small concentrations upwards so the plot has a nice x-axis:
-    curves$concentration[1:gap_idx] <-
-        curves$concentration[1:gap_idx] / model$conc_adjust_factor
-
-    list(
-        curves = curves,
-        axis_break_conc = axis_break_conc
-    )
-}
diff --git a/R/plot_effect.R b/R/plot_effect.R
index 942b3cc..0fb0293 100644
--- a/R/plot_effect.R
+++ b/R/plot_effect.R
@@ -52,15 +52,17 @@ plot_effect <- function(model,
         which <- which[which %in% valid_names]
     }
 
-    temp <- adjust_plot_concentrations(model)
-    curves <- temp$curves
-    log_ticks <- get_log_ticks(curves$concentration)
-    concentration <- c(curves$concentration[1], model$args$concentration[-1])
+    curves <- model$curves
+    log_ticks <- get_log_ticks(curves$concentration_for_plots)
+    point_concentration <- c(
+        curves$concentration_for_plots[1],
+        model$args$concentration[-1]
+    )
 
     plot(
         NA,
         NA,
-        xlim = range(curves$concentration, na.rm = TRUE),
+        xlim = range(curves$concentration_for_plots, na.rm = TRUE),
         ylim = c(0, model$args$effect_max),
         log = "x",
         xlab = xlab,
@@ -75,7 +77,7 @@ plot_effect <- function(model,
     # are on top of solid lines for better visibility.
     if ("effect_tox_observed" %in% which) {
         points(
-            concentration,
+            point_concentration,
             model$args$effect_tox_observed,
             pch = 16,
             col = "blue"
@@ -83,14 +85,14 @@ plot_effect <- function(model,
     }
     if ("effect_tox_sys" %in% which) {
         lines(
-            curves$concentration,
+            curves$concentration_for_plots,
             curves$effect_tox_sys,
             col = "blue"
         )
     }
     if ("effect_tox" %in% which) {
         lines(
-            curves$concentration,
+            curves$concentration_for_plots,
             curves$effect_tox,
             col = "deepskyblue",
             lty = 2
@@ -98,7 +100,7 @@ plot_effect <- function(model,
     }
     if ("effect_tox_LL5" %in% which) {
         lines(
-            curves$concentration,
+            curves$concentration_for_plots,
             curves$effect_tox_LL5,
             col = "darkblue",
             lty = 3
@@ -107,7 +109,7 @@ plot_effect <- function(model,
     if (model$with_env) {
         if ("effect_tox_env_observed" %in% which) {
             points(
-                concentration,
+                point_concentration,
                 model$args$effect_tox_env_observed,
                 pch = 16,
                 col = "red"
@@ -115,14 +117,14 @@ plot_effect <- function(model,
         }
         if ("effect_tox_env_sys" %in% which) {
             lines(
-                curves$concentration,
+                curves$concentration_for_plots,
                 curves$effect_tox_env_sys,
                 col = "red"
             )
         }
         if ("effect_tox_env" %in% which) {
             lines(
-                curves$concentration,
+                curves$concentration_for_plots,
                 curves$effect_tox_env,
                 col = "orange",
                 lty = 2
@@ -130,7 +132,7 @@ plot_effect <- function(model,
         }
         if ("effect_tox_env_LL5" %in% which) {
             lines(
-                curves$concentration,
+                curves$concentration_for_plots,
                 curves$effect_tox_env_LL5,
                 col = "darkred",
                 lty = 3
@@ -145,7 +147,7 @@ plot_effect <- function(model,
          col = NA, col.ticks = par("fg"))
     axis(1, at = log_ticks$minor, labels = FALSE, tcl = -0.25,
          col = NA, col.ticks = par("fg"))
-    plotrix::axis.break(1, breakpos = temp$axis_break_conc)
+    plotrix::axis.break(1, breakpos = model$axis_break_conc)
     axis(2, col = NA, col.ticks = par("fg"), las = 1)
 
     if (show_legend) {
diff --git a/R/plot_stress.R b/R/plot_stress.R
index 152b0a3..7c2f5da 100644
--- a/R/plot_stress.R
+++ b/R/plot_stress.R
@@ -51,10 +51,12 @@ plot_stress <- function(model,
         which <- which[which %in% valid_names]
     }
 
-    temp <- adjust_plot_concentrations(model)
-    curves <- temp$curves
-    log_ticks <- get_log_ticks(curves$concentration)
-    concentration <- c(curves$concentration[1], model$args$concentration[-1])
+    curves <- model$curves
+    log_ticks <- get_log_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)
@@ -62,7 +64,7 @@ plot_stress <- function(model,
     plot(
         NA,
         NA,
-        xlim = range(curves$concentration, na.rm = TRUE),
+        xlim = range(curves$concentration_for_plots, na.rm = TRUE),
         ylim = c(0, ymax),
         log = "x",
         xlab = xlab,
@@ -77,7 +79,7 @@ plot_stress <- function(model,
     # are on top of solid lines for better visibility.
     if ("sys_tox_observed" %in% which) {
         points(
-            concentration,
+            point_concentration,
             model$sys_tox_observed,
             pch = 16,
             col = "steelblue3"
@@ -85,14 +87,14 @@ plot_stress <- function(model,
     }
     if ("stress_tox_sys" %in% which) {
         lines(
-            curves$concentration,
+            curves$concentration_for_plots,
             curves$stress_tox_sys,
             col = "blue"
         )
     }
     if ("stress_tox" %in% which) {
         lines(
-            curves$concentration,
+            curves$concentration_for_plots,
             curves$stress_tox,
             col = "deepskyblue",
             lty = 2
@@ -100,7 +102,7 @@ plot_stress <- function(model,
     }
     if ("sys_tox" %in% which) {
         lines(
-            curves$concentration,
+            curves$concentration_for_plots,
             curves$sys_tox,
             col = "steelblue3"
         )
@@ -108,7 +110,7 @@ plot_stress <- function(model,
     if (model$with_env) {
         if ("sys_tox_env_observed" %in% which) {
             points(
-                concentration,
+                point_concentration,
                 model$sys_tox_env_observed,
                 pch = 16,
                 col = "violetred"
@@ -116,14 +118,14 @@ plot_stress <- function(model,
         }
         if ("stress_tox_env_sys" %in% which) {
             lines(
-                curves$concentration,
+                curves$concentration_for_plots,
                 curves$stress_tox_env_sys,
                 col = "red"
             )
         }
         if ("stress_env" %in% which) {
             lines(
-                curves$concentration,
+                curves$concentration_for_plots,
                 curves$stress_env,
                 col = "forestgreen",
                 lty = 3
@@ -131,7 +133,7 @@ plot_stress <- function(model,
         }
         if ("stress_tox_env" %in% which) {
             lines(
-                curves$concentration,
+                curves$concentration_for_plots,
                 curves$stress_tox_env,
                 col = "orange",
                 lty = 2
@@ -139,7 +141,7 @@ plot_stress <- function(model,
         }
         if ("sys_tox_env" %in% which) {
             lines(
-                curves$concentration,
+                curves$concentration_for_plots,
                 curves$sys_tox_env,
                 col = "violetred"
             )
@@ -153,7 +155,7 @@ plot_stress <- function(model,
          col = NA, col.ticks = par("fg"))
     axis(1, at = log_ticks$minor, labels = FALSE, tcl = -0.25,
          col = NA, col.ticks = par("fg"))
-    plotrix::axis.break(1, breakpos = temp$axis_break_conc)
+    plotrix::axis.break(1, breakpos = model$axis_break_conc)
     axis(2, col = NA, col.ticks = par("fg"), las = 1)
 
     if (show_legend) {
diff --git a/man/ec.Rd b/man/ec.Rd
index 0cacd6a..9ddf03f 100644
--- a/man/ec.Rd
+++ b/man/ec.Rd
@@ -4,12 +4,13 @@
 \alias{ec}
 \title{Effect Concentrations}
 \usage{
-ec(model, response_name, response_level)
+ec(model, response_name, response_level, reference, warn = TRUE)
 }
 \arguments{
 \item{model}{This can be one of three types of objects: Either the output of
 \code{\link{ecxsys}} or the output of \code{\link{predict_ecxsys}} or a
 data frame with a "concentration" column and a \code{response_name} column.
+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
@@ -19,10 +20,21 @@ calculate the EC. Must be one of \code{colnames(model$curves)}.}
 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.}
+
+\item{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.}
+
+\item{warn}{Logical. Should the function emit a warning if the calculation of
+the effect concentration is not possible?}
 }
 \value{
-A list containing the response concentration and the corresponding
-  response 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.
 }
 \description{
 Estimate the concentration to reach a certain effect or stress level relative
@@ -39,8 +51,8 @@ increasing concentration, i.e. all \code{effect_*} curves and also
 unexpected results, if any.
 }
 \examples{
-# Calculate the EC_10, the concentration where the effect falls
-# below 90 \% of the effect in the control.
+# 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),
@@ -64,4 +76,8 @@ df_custom <- data.frame(
 )
 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)
+
 }
diff --git a/man/ecxsys.Rd b/man/ecxsys.Rd
index 4d5bed1..f8a83b8 100644
--- a/man/ecxsys.Rd
+++ b/man/ecxsys.Rd
@@ -60,7 +60,7 @@ A list (of class ecxsys) containing many different objects of which
     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{use_for_plotting} is used by the
+    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.}
   }
diff --git a/tests/testthat/test-ec.R b/tests/testthat/test-ec.R
index 0371c7d..fe5bbfc 100644
--- a/tests/testthat/test-ec.R
+++ b/tests/testthat/test-ec.R
@@ -17,12 +17,16 @@
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 
+# shared model for some of the tests
+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)
+)
+
+
 test_that("all input formats produce identical models", {
-    model <- ecxsys(
-        concentration = c(0, 0.05, 0.5, 5, 30),
-        hormesis_concentration = 0.5,
-        effect_tox_observed = c(90, 81, 92, 28, 0)
-    )
     # using the ecxsys() output or the curves therein directly:
     ec10_a <- ec(model, "effect_tox_sys", 10)
     ec10_b <- ec(model$curves, "effect_tox_sys", 10)
@@ -46,64 +50,72 @@ test_that("all input formats produce identical models", {
 
 
 test_that("ec values have not changed", {
-    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)
-    )
     expect_equal(
         ec(model, "effect_tox_sys", 50),
-        list(response_value = 44.95368, concentration = 3.375735),
+        list(effect = 44.95368, concentration = 3.375735),
         tolerance = 1e-4
     )
     expect_equal(
         ec(model, "effect_tox_sys", 10),
-        list(response_value = 80.91662, concentration = 1.098648),
+        list(effect = 80.91662, concentration = 1.098648),
         tolerance = 1e-4
     )
     expect_equal(
-        ec(model, "effect_tox_sys", 5),
-        list(response_value = 85.41198, concentration = 0.005502182),
+        ec(model, "effect_tox", 100/3),
+        list(effect = 66.66667, concentration = 1.902125),
         tolerance = 1e-4
     )
     expect_equal(
-        ec(model, "effect_tox_sys", 1),
-        list(response_value = 89.00828, concentration = 8.53288e-05),
+        ec(model, "effect_tox_LL5", 42),
+        list(effect = 52.2, concentration = 2.054426),
         tolerance = 1e-4
     )
     expect_equal(
         ec(model, "effect_tox_env_sys", 50),
-        list(response_value = 14.67725, concentration = 1.299516),
+        list(effect = 14.67725, concentration = 1.299516),
         tolerance = 1e-4
     )
     expect_equal(
         ec(model, "effect_tox_env_sys", 10),
-        list(response_value = 26.41904, concentration = 0.0008571244),
+        list(effect = 26.41904, concentration = 0.0008571244),
         tolerance = 1e-4
     )
     expect_equal(
-        ec(model, "effect_tox_env_sys", 5),
-        list(response_value = 27.88677, concentration = 0.0001044245),
+        ec(model, "effect_tox_env", 67.89),
+        list(effect = 24.51453, concentration = 0.7890892),
         tolerance = 1e-4
     )
     expect_equal(
-        ec(model, "effect_tox_env_sys", 1),
-        list(response_value = 29.06095, concentration = 1.388112e-06),
+        ec(model, "effect_tox_env_LL5", 3.14159),
+        list(effect = 28.73466, concentration = 0.7267524),
         tolerance = 1e-4
     )
 })
 
 
-test_that("ec fails when response_level is outside the range of the curve", {
+test_that("ec 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)
     )
-    expect_error(
+    expect_warning(
         ec(model, "effect_tox_sys", 90),
         "You could try using predict_ecxsys() to predict more values",
         fixed = TRUE
     )
 })
+
+
+test_that("reference argument works", {
+    expect_equal(
+        ec(model, "effect_tox_LL5", 50, reference = 100),
+        list(effect = 50, concentration = 2.208119),
+        tolerance = 1e-4
+    )
+    expect_equal(
+        ec(model, "effect_tox_LL5", 50, reference = 75),
+        list(effect = 37.5, concentration = 3.342715),
+        tolerance = 1e-4
+    )
+})
diff --git a/tests/testthat/test-ecxsys.R b/tests/testthat/test-ecxsys.R
index d04e16f..ffc2608 100644
--- a/tests/testthat/test-ecxsys.R
+++ b/tests/testthat/test-ecxsys.R
@@ -161,7 +161,7 @@ test_that("the curves have not changed", {
         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),
-        use_for_plotting = c(TRUE, TRUE, TRUE, TRUE, TRUE)
+        concentration_for_plots = c(0.0001, 0.03004, 0.30512, 3.02551, 30)
     )
     class(new_curves) <- class(reference_curves)
     expect_equal(new_curves, reference_curves, tolerance = 1e-3)
@@ -171,7 +171,7 @@ test_that("the curves have not changed", {
 test_that("the returned fn works the same way as internally", {
     # I don't know why it would fail but it doesn't hurt to test it.
     curves <- mod$curves
-    curves$use_for_plotting <- NULL  # remove column "use_for_plotting"
+    curves$concentration_for_plots <- NULL
     expect_identical(curves, predict_ecxsys(mod, curves$concentration))
 })
 
-- 
GitLab