Improve predict_mixture

description

Multiple changes to predict_mixture.

• The concentration vectors are checked for NAs, which raise an error when found.
• The use of drc::ED() was replaced with a custom inverse weibull function, which does the same but better.
• 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.

system details

• RStudio Version: Version 1.2.5033
• 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/predict_mixture.R

@@ -27,10 +27,10 @@
 #' 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,model_2 The ecxsys models of the toxicants.
-#' @param concentration_1,concentration_2 The concentrations of the toxicants in
+#' @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.
@@ -41,96 +41,89 @@
 #' @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,
+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_1) > 0,
-        length(concentration_2) > 0,
-        length(concentration_1) == length(concentration_2),
-        all(!is.na(concentration_1)),
-        all(!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
+    # Convert the model_b concentration into an equivalent model_a concentration
     # and vice versa.
-    concentration_2_equivalent <- W1.2_inverse(
-        model_1$effect_tox_mod,
-        predicted_model_2$effect_tox / model_2$args$effect_max
+    concentration_b_equivalent <- W1.2_inverse(
+        model_a$effect_tox_mod,
+        predicted_model_b$effect_tox / model_b$args$effect_max
     )
-    effect_tox_ca_1 <- predict(
-        model_1$effect_tox_mod,
-        data.frame(concentration = concentration_1 + concentration_2_equivalent)
+    effect_tox_ca_a <- predict(
+        model_a$effect_tox_mod,
+        data.frame(concentration = concentration_a + concentration_b_equivalent)
     )
-    stress_tox_ca_1 <- effect_to_stress(effect_tox_ca_1)
+    stress_tox_ca_a <- effect_to_stress(effect_tox_ca_a)
 
-    concentration_1_equivalent <- W1.2_inverse(
-        model_2$effect_tox_mod,
-        predicted_model_1$effect_tox / model_1$args$effect_max
+    concentration_a_equivalent <- W1.2_inverse(
+        model_b$effect_tox_mod,
+        predicted_model_a$effect_tox / model_a$args$effect_max
     )
-    effect_tox_ca_2 <- predict(
-        model_2$effect_tox_mod,
-        data.frame(concentration = concentration_2 + concentration_1_equivalent)
+    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 + sys_2) / 2
+    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
@@ -138,9 +131,9 @@ predict_mixture <- function(model_1,
     stress_total <- stress_tox_total + sys_total
     mixture_effect <- stress_to_effect(stress_total) * effect_max
 
-    # unname() to remove the name when concentration_1 is a single number.
+    # unname() to remove the name when concentration_a is a single number.
     mixture_effect <- unname(mixture_effect)
-    data.frame(concentration_1, concentration_2, mixture_effect)
+    data.frame(concentration_a, concentration_b, mixture_effect)
 }