diff --git a/NAMESPACE b/NAMESPACE
index 06741e7eab9a60266f9993496c37c8b444cd9238..58db063863ead9fe8cb2bc7ceebf64bc51392656 100644
--- a/NAMESPACE
+++ b/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)
diff --git a/NEWS.md b/NEWS.md
index 22c5af0c2f0954f1d37c7daa13eaea418127bf86..21adea503dc8fa4e56ec37900e266b814ea88ee5 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,6 +1,7 @@
# 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
diff --git a/R/predict_mixture.R b/R/predict_mixture.R
new file mode 100644
index 0000000000000000000000000000000000000000..7dc546fa72ad078e44f60140d711075b42e4db82
--- /dev/null
+++ b/R/predict_mixture.R
@@ -0,0 +1,104 @@
+# 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)
+}
diff --git a/man/predict_mixture.Rd b/man/predict_mixture.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..901eba64e27c25c5a0d715960038d66fce0efe95
--- /dev/null
+++ b/man/predict_mixture.Rd
@@ -0,0 +1,53 @@
+% 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
+)
+}
diff --git a/tests/testthat/test-predict_mixture.R b/tests/testthat/test-predict_mixture.R
new file mode 100644
index 0000000000000000000000000000000000000000..f482f5c5e3dd0e97052390678fac0ad4d9cf126e
--- /dev/null
+++ b/tests/testthat/test-predict_mixture.R
@@ -0,0 +1,34 @@
+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)
+})