Sort out the options and warnings mess,

closes #13

R/ecxsys.R

-# Note about the setting and resetting of options():
-# In drc version 3.0-1 a call to drm changes the options "warn" and some others.
-# That is why I cache the options and reset them every time after calling drm.
-# I filed a bug report and they fixed it. In the next version of drc I may
-# remove some of that manual option handling. Not, however, when I change them
-# in tryCatch(). I must make sure to require at least a drc version that has
-# that fix.
+# 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.
 
 
 #' ECx-SyS
@@ -19,7 +18,7 @@
 #'
 #' The vectors \code{concentration}, \code{effect_tox_observed} and
 #' \code{effect_tox_env_observed} (if provided) must be of equal length and
-#' should be sorted by increasing concentration.
+#' sorted by increasing concentration.
 #'
 #' @param concentration A vector of concentrations, one of which must be 0 to
 #'   indicate the control. Should be sorted in ascending order, otherwise it
@@ -68,6 +67,10 @@ ecxsys <- function(concentration,
     output <- list(args = as.list(environment()))
     class(output) <- c("ecxsys", class(output))
 
+    original_options <- options()
+    on.exit(reset_options(original_options))
+
+
     # input validation ----------------------------------------------------
     if (effect_max <= 0) {
         stop("effect_max must be >= 0")
@@ -120,26 +123,14 @@ ecxsys <- function(concentration,
     # second lowest concentration. This is required to approximate 0 because
     # of the logarithmic axis.
     if (any(concentration < 0)) {
-        stop("Concentration must be >= 0")
+        stop("Concentrations must be >= 0")
     } else if (min(concentration) > 0) {
-        warning("No control is given and therefore the smallest concentration ",
-                "is assumed to be the control.")
-        min_conc <- min(concentration)
-        concentration[which.min(concentration)] <- 0
+        stop("No control is given. The first concentration must be 0.")
     } else {
         min_conc <- 10 ^ floor(log10(concentration[2]) - conc_shift)
     }
     if (is.unsorted(concentration)) {
-        warning("The concentrations are not sorted in increasing order. The ",
-                "provided effect vectors will be sorted by concentration.")
-        od <- order(concentration)
-        concentration <- concentration[od]
-        effect_tox_env_observed <- effect_tox_env_observed[od]
-        effect_tox_observed <- effect_tox_observed[od]
-    }
-    if (effect_tox_observed[length(effect_tox_observed)] > 0) {
-        warning("It is advised to complete the curve down to zero for ",
-                "optimal prediction.")
+        stop("The values must be sorted by increasing concentration.")
     }
 
 
@@ -166,14 +157,14 @@ ecxsys <- function(concentration,
     )
     conc_interpolated <- 10^temp$x
     effect_tox_observed_interpolated_simple_model <- temp$y
-    original_options <- options()
+    option_warn_original <- getOption("warn")
     effect_tox_mod_simple <- drc::drm(
         effect_tox_observed_interpolated_simple_model ~ conc_interpolated,
         fct = drc::LL.5(fixed = c(
             NA, 0, effect_tox_observed_averaged[1], NA, NA
         ))
     )
-    options(original_options)  #  because drm modifies options
+    options(warn = option_warn_original)
     effect_tox_simple <- predict(effect_tox_mod_simple, data.frame(concentration))
     output$effect_tox_mod_simple <- effect_tox_mod_simple
     output$effect_tox_simple <- effect_tox_simple * effect_max
@@ -183,14 +174,14 @@ ecxsys <- function(concentration,
             effect_tox_env_observed_averaged,
             xout = temp$x
         )$y
-        original_options <- options()
+        option_warn_original <- getOption("warn")
         effect_tox_env_mod_simple <- drc::drm(
             effect_tox_env_observed_interpolated_simple_model ~ conc_interpolated,
             fct = drc::LL.5(fixed = c(
                 NA, 0, effect_tox_env_observed_averaged[1], NA, NA
             ))
         )
-        options(original_options)  #  because drm modifies options
+        options(warn = option_warn_original)
         effect_tox_env_simple <- predict(
             effect_tox_env_mod_simple,
             data.frame(concentration)
@@ -253,12 +244,12 @@ ecxsys <- function(concentration,
     effect_tox[1] <- 1
     effect_to_fit_idx <- 2:(hormesis_index - 1)
     effect_tox[effect_to_fit_idx] <- NA
-    original_options <- options()
+    option_warn_original <- getOption("warn")
     effect_tox_mod <- drc::drm(
         effect_tox ~ concentration,
         fct = drc::W1.2()
     )
-    options(original_options)  #  because drm modifies options
+    options(warn = option_warn_original)
     effect_tox <- predict(
         effect_tox_mod,
         data.frame(concentration = concentration)
@@ -345,28 +336,54 @@ ecxsys <- function(concentration,
 }
 
 
+reset_options <- function(original_options) {
+    # Reset all the options which have changed.
+
+    # You may ask why I don't just cache and reset the options like this:
+    # original_options <- options()
+    # options(original_options)
+    # The reason is that when I do this and ecxsys generates warnings then those
+    # warnings don't show up in the console. I don't know why, but resetting
+    # only the options which have changed alleviates that problem.
+
+    changed <- list()
+    for (n in names(original_options)) {
+        orig_opt <- original_options[[n]]
+        if (!identical(orig_opt, getOption(n))) {
+            changed[n] <- orig_opt
+        }
+    }
+    options(changed)
+}
+
+
 fit_sys <- function(stress_tox, sys) {
+    result_name <- paste0(deparse(substitute(sys)), "_mod")
+    # There is no other way to suppress that one error message
+    # from optim except by changing the options temporarily.
+    warn_error_original <- list(
+        warn = getOption("warn"),
+        show.error.messages = getOption("show.error.messages")
+    )
+    options(show.error.messages = FALSE)
     tryCatch(
-        {
-            # There is no other way to suppress that one error message
-            # from optim except by changing the options temporarily.
-            original_options <- options()
-            options(show.error.messages = FALSE)
-            drc::drm(sys ~ stress_tox, fct = drc::W1.3())
-        },
+        drc::drm(sys ~ stress_tox, fct = drc::W1.3()),
         error = function(e) {
             # Failure to converge often happens when all or almost all sys
             # stress values are zero. Fitting a linear model in this case seems
             # to be the most appropriate remedy.
-            # I would also like to generate a warning in this case but for some
-            # unholy reason it gets printed sometimes and other times it's not.
-            # So for consistency I think it is better to leave it out
-            # completely until I figure out what's going on.
+            options(warn_error_original)
+            warning(
+                "Using a horizontal linear model for ",
+                result_name,
+                " because the Weibull model did not converge.",
+                call. = FALSE
+            )
             stress_tox <- range(stress_tox)
             sys <- c(0, 0)
             lm(sys ~ stress_tox)
         },
-        finally = options(original_options)
+        finally = options(warn_error_original)
     )
 }
 

man/ecxsys.Rd

@@ -52,7 +52,7 @@ 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
-should be sorted by increasing concentration.
+sorted by increasing concentration.
 }
 \examples{
 model <- ecxsys(

tests/testthat/test_ecxsys.R

@@ -70,26 +70,6 @@ test_that("min(concentration) == 0 is shifted the correct amount", {
 })
 
 
-test_that("min(concentration) > 0 is conserved", {
-    expect_warning(
-        ecxsys(
-            concentration = c(0.0005, 0.03, 0.3, 3, 10),
-            hormesis_concentration = 0.3,
-            effect_tox_observed = c(85, 76, 94, 35, 0)
-        )
-    )
-    suppressWarnings({
-        mod <- ecxsys(
-            concentration = c(0.0005, 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)
-        )
-    })
-    expect_equal(mod$curves$concentration[1] * 10^5, 0.0005)
-})
-
-
 test_that("the discrete results have not changed", {
     expect_equal(
         round(mod$effect_tox_simple, 3),
@@ -258,16 +238,12 @@ test_that("effect_tox_env_observed can be left out", {
 })
 
 
-# This is commented out for now until I'm finished refactoring ecxsys and have
-# time to tackle the issue around the warnings which don't show up in the
-# console. See issue #13.
-# Also remember to test if the error produced by optim is suppressed or not.
-# test_that("model not converging produces no errors", {
-#     expect_warning(
-#         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)
-#         )
-#     )
-# })
+test_that("sys model not converging produces a warning", {
+    expect_warning(
+        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)
+        )
+    )
+})