diff --git a/DESCRIPTION b/DESCRIPTION
index 8610a2b08012c934edfcb5b64b782a98315e9f1d..422ecaa2ca550bd37402085b183fab19ea0aa1be 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
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"),
diff --git a/NEWS.md b/NEWS.md
index b954e4d6be1cda2075ff18af1db1fa5baebe4fde..2b8b99f441931ee20633ddb96ff63ca1c671240e 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,10 @@
+# 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.
diff --git a/R/predict_mixture.R b/R/predict_mixture.R
index f70f6a60e381d8cdd0de15759dd3ac24a83fe7e3..c79582c6ef3dea42682f09cc94e9dc105fe8fd0f 100644
--- a/R/predict_mixture.R
+++ b/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))
}
diff --git a/man/predict_mixture.Rd b/man/predict_mixture.Rd
index 8595ac80da97ffe27d7fcc38bd4e9361212800cc..297c291ba1662e2ff5c3320bb3396fed24e650e6 100644
--- a/man/predict_mixture.Rd
+++ b/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)
}
diff --git a/tests/testthat/test-predict_mixture.R b/tests/testthat/test-predict_mixture.R
index 715697f791ee352f83cbf77a79c15fbd8ff880d1..2ad99b6f8d49429b6b8943e14ba3e1fbdb410c5e 100644
--- a/tests/testthat/test-predict_mixture.R
+++ b/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)
+})