Merge branch 'develop' into 'master'

Improve predict_mixture

See merge request !25

DESCRIPTION

 Package: stressaddition
 Type: Package
 Title: Modeling Tri-Phasic Concentration-Response Relationships
-Version: 2.2.1
-Date: 2020-03-03
+Version: 2.3.0
+Date: 2020-03-11
 Authors@R: c(person("Sebastian", 
                     "Henz", 
                     role = c("aut", "cre"), 

NEWS.md

+# stressaddition 2.3.0
+
+* `predict_mixture()` accepts multiple values for the concentration of the second toxicant. Both concentration vectors must be the same length.
+* `predict_mixture()` returns a data frame with the concentrations and effects. Previously it was only a vector of effects.
+* `predict_mixture()` received a new argument "effect_max" which scales the returned effect values.
+* Renamed the arguments of `predict_mixture()` to use underscore letters a and b instad of 1 and 2. For example model_1 is now model_a.
+
 # stressaddition 2.2.1
 
 * Improve documentation of `predict_mixture()` and include example of symmetry.

R/predict_mixture.R

@@ -22,127 +22,126 @@
 #' Given the ecxsys models of two toxicants this method predicts the effects of
 #' different mixtures of both.
 #'
-#' The prediction happens at one or multiple concentrations of toxicant 1 and
-#' one concentration of toxicant 2. This allows for easy plotting over the
-#' concentrations of toxicant 1.
-#'
 #' The predictions are symmetric, i.e. it does not matter which of the toxicant
-#' models is 1 or 2 as long as the concentration arguments are supplied
-#' in the right order. See the example below.
+#' models is 1 or 2 as long as the concentration arguments are supplied in the
+#' right order. See the example below.
 #'
 #' This method is only suitable for experiments without environmental stress.
-#' Any environmental stress in \code{model_1} or \code{model_2} is ignored.
+#' Any environmental stress in \code{model_a} or \code{model_b} is ignored.
 #'
-#' @param model_1 The ecxsys model of the toxicant that varies in concentration.
-#' @param model_2 The ecxsys model of the toxicant with a fixed concentration.
-#' @param concentration_1 The concentration of toxicant 1 in the mixture.
-#' @param concentration_2 The concentration of toxicant 2 in the mixture.
-#' @param proportion_ca How much of the combined toxicant stress is determined
-#'   by concentration addition. Must be between 0 and 1.
+#' @param model_a,model_b The ecxsys models of the toxicants.
+#' @param concentration_a,concentration_b The concentrations of the toxicants in
+#'   the mixture. Both vectors must be the same length.
+#' @param proportion_ca The proportion of concentration addition in the
+#'   calculation of the toxicant stress of the mixture. Must be between 0 and 1.
+#' @param effect_max Controls the scaling of the result. This represents the
+#'   maximum value the effect could possibly reach. For survival data in percent
+#'   this should be 100 (the default).
 #'
-#' @return A vector of the effects of the mixture, scaled to the effect_max
-#'   value of model_1.
+#' @return A data frame with columns of the supplied concentrations and the
+#'   corresponding mixture effects.
 #'
-#' @examples toxicant_1  <- ecxsys(
+#' @examples toxicant_a  <- ecxsys(
 #'     concentration = c(0, 0.05, 0.5, 5, 30),
 #'     hormesis_concentration = 0.5,
 #'     effect_tox_observed = c(90, 81, 92, 28, 0),
 #' )
-#' toxicant_2  <- ecxsys(
+#' toxicant_b  <- ecxsys(
 #'     concentration = c(0, 0.1, 1, 10, 100, 1000),
 #'     hormesis_concentration = 10,
 #'     effect_tox_observed = c(26, 25, 24, 27, 5, 0),
 #'     effect_max = 30
 #' )
 #' predict_mixture(
-#'     toxicant_1 ,
-#'     toxicant_2 ,
+#'     toxicant_a ,
+#'     toxicant_b ,
 #'     c(0, 0.02, 0.2, 2, 20),
-#'     3
+#'     c(0, 0.03, 0.4, 5, 10)
 #' )
 #'
 #' # Example of symmetric prediction:
-#' conc_1 <- c(0, 0.03, 0.3, 3)
-#' conc_2 <- 5.5
+#' conc_a <- c(0, 0.03, 0.3, 3)
+#' conc_b <- rep(5.5, 4)
 #' prop_ca <- 0.75
-#' effect_a <- predict_mixture(toxicant_1 , toxicant_2 , conc_1, conc_2, prop_ca)
-#' effect_b <- sapply(
-#'     conc_1,
-#'     function(x) predict_mixture(toxicant_2 , toxicant_1 , conc_2, x, prop_ca)
-#' )
-#' # The effect_max values of the models are different. effect_b is scaled to
-#' # the one from toxicant 2 but to compare the results effect_b must be scaled
-#' # to the effect_max of toxicant 1:
-#' effect_b <- effect_b / toxicant_2 $args$effect_max * toxicant_1$args$effect_max
-#' identical(effect_a, effect_b)
+#' effect_a <- predict_mixture(toxicant_a , toxicant_b , conc_a, conc_b, prop_ca)
+#' effect_b <- predict_mixture(toxicant_b , toxicant_a , conc_b, conc_a, prop_ca)
+#' identical(effect_a$mixture_effect, effect_b$mixture_effect)
 #'
 #' @export
-predict_mixture <- function(model_1,
-                            model_2,
-                            concentration_1,
-                            concentration_2,
-                            proportion_ca = 0.5) {
+predict_mixture <- function(model_a,
+                            model_b,
+                            concentration_a,
+                            concentration_b,
+                            proportion_ca = 0.5,
+                            effect_max = 100) {
     stopifnot(
-        inherits(model_1, "ecxsys"),
-        inherits(model_2, "ecxsys"),
-        is.numeric(concentration_1),
-        is.numeric(concentration_2),
-        length(concentration_2) == 1,
-        !is.na(concentration_2),
+        inherits(model_a, "ecxsys"),
+        inherits(model_b, "ecxsys"),
+        is.numeric(concentration_a),
+        is.numeric(concentration_b),
+        length(concentration_a) > 0,
+        length(concentration_b) > 0,
+        length(concentration_a) == length(concentration_b),
+        all(!is.na(concentration_a)),
+        all(!is.na(concentration_b)),
         proportion_ca >= 0,
         proportion_ca <= 1,
-        model_1$args$p == model_2$args$p,
-        model_1$args$q == model_2$args$q
+        model_a$args$p == model_b$args$p,
+        model_a$args$q == model_b$args$q
     )
 
-    predicted_model_1 <- predict_ecxsys(model_1, concentration_1)
-    predicted_model_2 <- predict_ecxsys(model_2, concentration_2)
+    predicted_model_a <- predict_ecxsys(model_a, concentration_a)
+    predicted_model_b <- predict_ecxsys(model_b, concentration_b)
 
     # tox stress ----------------------------------------------------------
-    stress_tox_sam <- predicted_model_1$stress_tox + predicted_model_2$stress_tox
+    stress_tox_sam <- predicted_model_a$stress_tox + predicted_model_b$stress_tox
 
-    # Convert the model_2 concentration into an equivalent model_1 concentration
-    # anc vice versa. Clamp the response levels because drc::ED can't deal with
-    # 0 and 100.
-    response_level_2 <- 100 - predicted_model_2$effect_tox / model_2$args$effect_max * 100
-    response_level_2 <- clamp(response_level_2, 1e-10, 100 - 1e-10)
-    concentration_2_equivalent <- drc::ED(
-        model_1$effect_tox_mod,
-        response_level_2,
-        display = FALSE
-    )[, "Estimate"]
-    effect_tox_ca_1 <- predict(
-        model_1$effect_tox_mod,
-        data.frame(concentration = concentration_1 + concentration_2_equivalent)
+    # Convert the model_b concentration into an equivalent model_a concentration
+    # and vice versa.
+    concentration_b_equivalent <- W1.2_inverse(
+        model_a$effect_tox_mod,
+        predicted_model_b$effect_tox / model_b$args$effect_max
     )
-    stress_tox_ca_1 <- effect_to_stress(effect_tox_ca_1)
+    effect_tox_ca_a <- predict(
+        model_a$effect_tox_mod,
+        data.frame(concentration = concentration_a + concentration_b_equivalent)
+    )
+    stress_tox_ca_a <- effect_to_stress(effect_tox_ca_a)
 
-    response_level_1 <- 100 - predicted_model_1$effect_tox / model_1$args$effect_max * 100
-    response_level_1 <- clamp(response_level_1, 1e-10, 100 - 1e-10)
-    concentration_1_equivalent <- drc::ED(
-        model_2$effect_tox_mod,
-        response_level_1,
-        display = FALSE
-    )[, "Estimate"]
-    effect_tox_ca_2 <- predict(
-        model_2$effect_tox_mod,
-        data.frame(concentration = concentration_2 + concentration_1_equivalent)
+    concentration_a_equivalent <- W1.2_inverse(
+        model_b$effect_tox_mod,
+        predicted_model_a$effect_tox / model_a$args$effect_max
+    )
+    effect_tox_ca_b <- predict(
+        model_b$effect_tox_mod,
+        data.frame(concentration = concentration_b + concentration_a_equivalent)
     )
-    stress_tox_ca_2 <- effect_to_stress(effect_tox_ca_2)
+    stress_tox_ca_b <- effect_to_stress(effect_tox_ca_b)
 
-    stress_tox_ca <- (stress_tox_ca_1 + stress_tox_ca_2) / 2
+    stress_tox_ca <- (stress_tox_ca_a + stress_tox_ca_b) / 2
 
     # sys -----------------------------------------------------------------
-    sys_1 <- predict(model_1$sys_tox_mod, data.frame(stress_tox = stress_tox_ca))
-    sys_2 <- predict(model_2$sys_tox_mod, data.frame(stress_tox = stress_tox_ca))
-    sys_total <- sys_1 * 0.5 + sys_2 * 0.5
+    sys_a <- predict(model_a$sys_tox_mod, data.frame(stress_tox = stress_tox_ca))
+    sys_b <- predict(model_b$sys_tox_mod, data.frame(stress_tox = stress_tox_ca))
+    sys_total <- (sys_a + sys_b) / 2
 
     # combined stress and result ------------------------------------------
     proportion_sam <- 1 - proportion_ca
     stress_tox_total <- stress_tox_ca * proportion_ca + stress_tox_sam * proportion_sam
     stress_total <- stress_tox_total + sys_total
-    effect_total <- stress_to_effect(stress_total) * model_1$args$effect_max
+    mixture_effect <- stress_to_effect(stress_total) * effect_max
+
+    # unname() to remove the name when concentration_a is a single number.
+    mixture_effect <- unname(mixture_effect)
+    data.frame(concentration_a, concentration_b, mixture_effect)
+}
+
 
-    # unname() to remove the name when concentration_1 is a single number.
-    unname(effect_total)
+W1.2_inverse <- function(mod, x) {
+    # Using drc::ED() with respLev = 0 does not work, which is why I use
+    # this function instead.
+    # x = 1 - respLev / 100
+    b <- mod$fit$par[1]
+    e <- mod$fit$par[2]
+    exp(log(-log(x)) / b + log(e))
 }

man/predict_mixture.Rd

@@ -5,77 +5,68 @@
 \title{Predict the effect of a mixture of two toxicants}
 \usage{
 predict_mixture(
-  model_1,
-  model_2,
-  concentration_1,
-  concentration_2,
-  proportion_ca = 0.5
+  model_a,
+  model_b,
+  concentration_a,
+  concentration_b,
+  proportion_ca = 0.5,
+  effect_max = 100
 )
 }
 \arguments{
-\item{model_1}{The ecxsys model of the toxicant that varies in concentration.}
+\item{model_a, model_b}{The ecxsys models of the toxicants.}
 
-\item{model_2}{The ecxsys model of the toxicant with a fixed concentration.}
+\item{concentration_a, concentration_b}{The concentrations of the toxicants in
+the mixture. Both vectors must be the same length.}
 
-\item{concentration_1}{The concentration of toxicant 1 in the mixture.}
+\item{proportion_ca}{The proportion of concentration addition in the
+calculation of the toxicant stress of the mixture. Must be between 0 and 1.}
 
-\item{concentration_2}{The concentration of toxicant 2 in the mixture.}
-
-\item{proportion_ca}{How much of the combined toxicant stress is determined
-by concentration addition. Must be between 0 and 1.}
+\item{effect_max}{Controls the scaling of the result. This represents the
+maximum value the effect could possibly reach. For survival data in percent
+this should be 100 (the default).}
 }
 \value{
-A vector of the effects of the mixture, scaled to the effect_max
-  value of model_1.
+A data frame with columns of the supplied concentrations and the
+  corresponding mixture effects.
 }
 \description{
 Given the ecxsys models of two toxicants this method predicts the effects of
 different mixtures of both.
 }
 \details{
-The prediction happens at one or multiple concentrations of toxicant 1 and
-one concentration of toxicant 2. This allows for easy plotting over the
-concentrations of toxicant 1.
-
 The predictions are symmetric, i.e. it does not matter which of the toxicant
-models is 1 or 2 as long as the concentration arguments are supplied
-in the right order. See the example below.
+models is 1 or 2 as long as the concentration arguments are supplied in the
+right order. See the example below.
 
 This method is only suitable for experiments without environmental stress.
-Any environmental stress in \code{model_1} or \code{model_2} is ignored.
+Any environmental stress in \code{model_a} or \code{model_b} is ignored.
 }
 \examples{
-toxicant_1  <- ecxsys(
+toxicant_a  <- ecxsys(
     concentration = c(0, 0.05, 0.5, 5, 30),
     hormesis_concentration = 0.5,
     effect_tox_observed = c(90, 81, 92, 28, 0),
 )
-toxicant_2  <- ecxsys(
+toxicant_b  <- ecxsys(
     concentration = c(0, 0.1, 1, 10, 100, 1000),
     hormesis_concentration = 10,
     effect_tox_observed = c(26, 25, 24, 27, 5, 0),
     effect_max = 30
 )
 predict_mixture(
-    toxicant_1 ,
-    toxicant_2 ,
+    toxicant_a ,
+    toxicant_b ,
     c(0, 0.02, 0.2, 2, 20),
-    3
+    c(0, 0.03, 0.4, 5, 10)
 )
 
 # Example of symmetric prediction:
-conc_1 <- c(0, 0.03, 0.3, 3)
-conc_2 <- 5.5
+conc_a <- c(0, 0.03, 0.3, 3)
+conc_b <- rep(5.5, 4)
 prop_ca <- 0.75
-effect_a <- predict_mixture(toxicant_1 , toxicant_2 , conc_1, conc_2, prop_ca)
-effect_b <- sapply(
-    conc_1,
-    function(x) predict_mixture(toxicant_2 , toxicant_1 , conc_2, x, prop_ca)
-)
-# The effect_max values of the models are different. effect_b is scaled to
-# the one from toxicant 2 but to compare the results effect_b must be scaled
-# to the effect_max of toxicant 1:
-effect_b <- effect_b / toxicant_2 $args$effect_max * toxicant_1$args$effect_max
-identical(effect_a, effect_b)
+effect_a <- predict_mixture(toxicant_a , toxicant_b , conc_a, conc_b, prop_ca)
+effect_b <- predict_mixture(toxicant_b , toxicant_a , conc_b, conc_a, prop_ca)
+identical(effect_a$mixture_effect, effect_b$mixture_effect)
 
 }

tests/testthat/test-predict_mixture.R

@@ -17,14 +17,14 @@
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 
-model_1 <- ecxsys(
+model_a <- 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),
     effect_max = 101
 )
-model_2 <- ecxsys(
+model_b <- ecxsys(
     concentration = c(0, 0.1, 1, 10, 100, 1000),
     hormesis_concentration = 10,
     effect_tox_observed = c(26, 25, 24, 27, 5, 0),
@@ -33,21 +33,58 @@ model_2 <- ecxsys(
 
 
 test_that("results have not changed", {
-    new <- predict_mixture(model_1, model_2, c(0, 0.01, 0.1, 1, 7, 15), 5, 0.3)
-    reference <- c(89.460324, 85.205168, 81.440100, 57.297855, 2.911545, 0)
+    # one concentration_b
+    new <- predict_mixture(
+        model_a,
+        model_b,
+        c(0, 0.01, 0.1, 1, 7, 15),
+        rep(5, 6),
+        0.3
+    )$mixture_effect
+    reference <- c(88.574578, 84.361552, 80.633762, 56.730550, 2.882718, 0)
+    expect_equal(new, reference, tolerance = 1e-5)
+
+    # diverse concentration_b
+    new <- predict_mixture(
+        model_a,
+        model_b,
+        c(0, 0.01, 0.1, 1, 7, 15),
+        c(0, 0.02, 0.2, 2, 14, 30),
+        0.3
+    )$mixture_effect
+    reference <- c(88.2698383, 79.9617127, 78.1574808, 65.7999834, 0.3861678, 0)
+    expect_equal(new, reference, tolerance = 1e-5)
+
+    # diverse concentration_b and custom effect_max
+    new <- predict_mixture(
+        model_a,
+        model_b,
+        c(0, 0.01, 0.1, 1, 7, 15),
+        c(0, 0.02, 0.2, 2, 14, 30),
+        0.3,
+        42
+    )$mixture_effect
+    reference <- c(88.2698383, 79.9617127, 78.1574808, 65.7999834, 0.3861678, 0) * 0.42
     expect_equal(new, reference, tolerance = 1e-5)
 })
 
 
 test_that("predictions are symmetric", {
-    conc_1 <- c(0, 10^seq(log10(0.001), log10(40), length.out = 50))
-    conc_2 <- 3.5
+    conc_a <- c(0, 10^seq(log10(0.001), log10(40), length.out = 50))
+    conc_b <- rep(3.5, length(conc_a))
     prop_ca <- 0.8
-    effect_12 <- predict_mixture(model_1, model_2, conc_1, conc_2, prop_ca)
-    effect_21 <- sapply(
-        conc_1,
-        function(x) predict_mixture(model_2, model_1, conc_2, x, prop_ca)
-    )
-    effect_21 <- effect_21 / model_2$args$effect_max * model_1$args$effect_max
+    effect_12 <- predict_mixture(model_a, model_b, conc_a, conc_b, prop_ca)$mixture_effect
+    effect_21 <- predict_mixture(model_b, model_a, conc_b, conc_a, prop_ca)$mixture_effect
     expect_equal(effect_12, effect_21)
 })
+
+
+test_that("effect_tox_mod is a W1.2 model", {
+    # This is a safeguard for if I decide to use something else than drc::W1.2
+    # for effect_tox_mod. It is extremely unlikely that this ever happens but
+    # better safe than sorry. Currently, the inverse only works for W1.2 models.
+    # If this test throws errors then I probably forget to adjust the inverse
+    # function accordingly.
+    expect_s3_class(model_a$effect_tox_mod$fct, "Weibull-1")
+    expect_equal(model_a$effect_tox_mod$fct$noParm, 2)
+})