Merge branch 'synergism' into 'develop'

Add predict_mixture() to predict the effects of a mixture of two toxicants.

Closes #23.

See merge request !16

NAMESPACE

@@ -6,6 +6,7 @@ export(effect_to_stress)
 export(plot_effect)
 export(plot_stress)
 export(predict_ecxsys)
+export(predict_mixture)
 export(stress_to_effect)
 import(grDevices)
 import(graphics)

NEWS.md

 # stressaddition (development version)
 
-* Fix documentation of `ecxsys()` and `predict_ecxsys`.
+* Added predict_mixture() for the prediction of the effects of mixtures of two toxicants.
+* Fixed documentation of `ecxsys()` and `predict_ecxsys`.
 
 # stressaddition 2.0.0
 

R/predict_mixture.R

+# TODO: reference to the paper TODO: explanatio how this works (mixture between
+# sam, and ca).
+# TODO: Add example with picture. Three plots: model_1, model_2, model_1 with
+# predicted line.
+
+#' Predict the effect of a mixture of two toxicants
+#'
+#' This method is only suitable for experiments without environmental stress.
+#' Any environmental stress in \code{model_1} or \code{model_2} 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.
+#'
+#' @return A vector of the effects of the mixture, scaled to the effect_max
+#'   value of model_1.
+#'
+#' @examples toxicant_model_1 <- ecxsys(
+#'     concentration = c(0, 0.03, 0.3, 3, 10),
+#'     hormesis_concentration = 0.3,
+#'     effect_tox_observed = c(85, 76, 94, 35, 0)
+#' )
+#' toxicant_model_2 <- 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_model_1,
+#'     toxicant_model_2,
+#'     c(0, 0.02, 0.2, 2, 20),
+#'     5
+#' )
+#' @export
+predict_mixture <- function(model_1,
+                            model_2,
+                            concentration_1,
+                            concentration_2,
+                            proportion_ca = 0.5) {
+    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),
+        proportion_ca >= 0,
+        proportion_ca <= 1
+    )
+
+    predicted_model_1 <- predict_ecxsys(model_1, concentration_1)
+    predicted_model_2 <- predict_ecxsys(model_2, concentration_2)
+
+    # tox stress ----------------------------------------------------------
+    stress_tox_sam <- predicted_model_1$stress_tox + predicted_model_2$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)
+    )
+
+    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)
+    )
+
+    stress_tox_ca <- effect_to_stress(effect_tox_ca_1 * 0.5 + effect_tox_ca_2 * 0.5)
+
+    # sys -----------------------------------------------------------------
+    stress_tox_max <- pmax(predicted_model_1$stress_tox, predicted_model_2$stress_tox)
+    sys_1 <- predict(model_1$sys_tox_mod, data.frame(stress_tox = stress_tox_max))
+    sys_2 <- predict(model_2$sys_tox_mod, data.frame(stress_tox = stress_tox_max))
+    sys_total <- sys_1 * 0.5 + sys_2 * 0.5
+
+    # 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
+
+    # unname() to remove the name when concentration_1 is a single number.
+    unname(effect_total)
+}

man/predict_mixture.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/predict_mixture.R
+\name{predict_mixture}
+\alias{predict_mixture}
+\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
+)
+}
+\arguments{
+\item{model_1}{The ecxsys model of the toxicant that varies in concentration.}
+
+\item{model_2}{The ecxsys model of the toxicant with a fixed concentration.}
+
+\item{concentration_1}{The concentration of toxicant 1 in the mixture.}
+
+\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.}
+}
+\value{
+A vector of the effects of the mixture, scaled to the effect_max
+  value of model_1.
+}
+\description{
+This method is only suitable for experiments without environmental stress.
+Any environmental stress in \code{model_1} or \code{model_2} is ignored.
+}
+\examples{
+toxicant_model_1 <- ecxsys(
+    concentration = c(0, 0.03, 0.3, 3, 10),
+    hormesis_concentration = 0.3,
+    effect_tox_observed = c(85, 76, 94, 35, 0)
+)
+toxicant_model_2 <- 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_model_1,
+    toxicant_model_2,
+    c(0, 0.02, 0.2, 2, 20),
+    5
+)
+}

tests/testthat/test-predict_mixture.R

+model_1 <- ecxsys(
+    concentration = c(0, 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),
+    effect_max = 101
+)
+model_2 <- 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
+)
+
+
+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(90.8690654, 86.1041220, 76.2504530, 49.5585562, 0.1958095, 0)
+    expect_equal(new, reference, tolerance = 1e-5)
+})
+
+
+test_that("predictions are symmetric", {
+    conc_1 <- c(0, 10^seq(log10(0.001), log10(20), length.out = 50))
+    conc_2 <- 3.5
+    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
+    expect_equal(effect_12, effect_21)
+})