Skip to content
Snippets Groups Projects
Commit fdca7613 authored by Sebastian Henz's avatar Sebastian Henz
Browse files

Enable check on the previous minor R release

closes #38
parent e3b4ac94
No related branches found
No related tags found
No related merge requests found
# Using the rocker/tidyverse image because it has devtools installed. # Using the rocker/tidyverse image because it has devtools and testthat.
# Use "document = FALSE" to catch potential mismatches between the roxygen
# comments and the .Rd files. This could happen when checking locally
# with "document = TRUE" (the default), which updates the .Rd files, but then
# forgetting to commit and push those updated .Rd files.
.check_template: &check_job_template .check_template: &check_job_template
only: only:
...@@ -6,15 +11,20 @@ ...@@ -6,15 +11,20 @@
- merge_requests - merge_requests
- tags - tags
script: script:
- R -e 'devtools::install_deps()' - R -e 'devtools::install_deps(quiet = TRUE)'
- R -e 'devtools::check(error_on = "warning")' - R -e 'devtools::check(error_on = "note", document = FALSE)'
# Earliest supported version
check-3.5.1:
image: rocker/tidyverse:3.5.1
<<: *check_job_template
# # latest version of the previous major R release # Last version before 4.0.0
# check-previous: check-3.6.3:
# image: rocker/tidyverse:3.5.3 image: rocker/tidyverse:3.6.3
# <<: *check_job_template <<: *check_job_template
# latest R release # Latest R release
check-latest: check-latest:
image: rocker/tidyverse:latest image: rocker/tidyverse:latest
<<: *check_job_template <<: *check_job_template
......
...@@ -26,7 +26,7 @@ Comment on any notes that come up during devtools::check(). ...@@ -26,7 +26,7 @@ Comment on any notes that come up during devtools::check().
* [ ] update version in DESCRIPTION * [ ] update version in DESCRIPTION
* [ ] update date in DESCRIPTION * [ ] update date in DESCRIPTION
* [ ] no remaining TODO, FIXME, or debug prints anywhere in the source files * [ ] no remaining TODO, FIXME, or debug prints anywhere in the source files
* [ ] `devtools::document()` * [ ] `devtools::document()` and commit the updated .Rd files
* [ ] `devtools::check()` without errors or warnings * [ ] `devtools::check()` without errors or warnings
<!-- <!--
......
Package: stressaddition Package: stressaddition
Type: Package Type: Package
Title: Modeling Tri-Phasic Concentration-Response Relationships Title: Modeling Tri-Phasic Concentration-Response Relationships
Version: 2.7.0 Version: 3.0.0
Date: 2020-04-17 Date: 2020-05-07
Authors@R: c(person("Sebastian", Authors@R: c(person("Sebastian",
"Henz", "Henz",
role = c("aut", "cre"), role = c("aut", "cre"),
...@@ -25,7 +25,7 @@ URL: https://git.ufz.de/oekotox/stressaddition ...@@ -25,7 +25,7 @@ URL: https://git.ufz.de/oekotox/stressaddition
Encoding: UTF-8 Encoding: UTF-8
LazyData: true LazyData: true
Depends: Depends:
R (>= 3.5) R (> 3.5)
Imports: Imports:
drc (>= 3.0), drc (>= 3.0),
plotrix plotrix
......
# Generated by roxygen2: do not edit by hand # Generated by roxygen2: do not edit by hand
export(ec)
export(ecxsys) export(ecxsys)
export(effect_to_stress) export(lc)
export(plot_effect) export(log10_ticks)
export(multi_tox)
export(plot_stress) export(plot_stress)
export(plot_survival)
export(predict_ecxsys) export(predict_ecxsys)
export(predict_mixture) export(stress_to_survival)
export(stress_to_effect) export(survival_to_stress)
import(grDevices) import(grDevices)
import(graphics) import(graphics)
import(stats) import(stats)
# stressaddition 3.0.0
## Breaking changes
* Renamed all instances of "effect" to "survival".
* Renamed all instances of "ec" to "lc".
* Renamed `predict_mixture()`, which was a temporary development name, to `multi_tox()`.
* The argument `proportion_ca` in the mixture model `multi_tox()` was renamed and its value reversed. It is now called `sa_contribution` and specifies the proportion of stress addition in the calculation of toxicant stress. To convert your code from the old version use this equation: sa_contribution = 1 - proportion_ca.
* Renamed `stress_tox_sam` to `stress_tox_sa` in the output of `multi_tox()`.
## Bugfixes
* Fixed a bug where `plot_stress()` with argument `which = NULL` would result in an error. Now it correctly draws the axes without data.
## New
* Exporte function `log10_ticks()` for calculating tick mark labels and positions on a base 10 logarithmic axis.
# stressaddition 2.7.0 # stressaddition 2.7.0
* Fixed some spelling mistakes. * Fixed some spelling mistakes.
......
This diff is collapsed.
...@@ -17,36 +17,12 @@ ...@@ -17,36 +17,12 @@
# along with this program. If not, see <https://www.gnu.org/licenses/>. # along with this program. If not, see <https://www.gnu.org/licenses/>.
# A collection of some internal helper functions. # A collection of internal helper functions which are small enough that they
# don't desere their own scripts. And which are also not exclusively used by one
# function because then they would be in the script of that function.
clamp <- function(x, lower = 0, upper = 1) { clamp <- function(x, lower = 0, upper = 1) {
# Returns lower if x < lower, returns upper if x > upper, else returns x # Returns lower if x < lower, returns upper if x > upper, else returns x
pmin(pmax(x, lower), upper) pmin(pmax(x, lower), upper)
} }
get_log_ticks <- function(x) {
# Calculate the positions of major and minor tick marks on a base 10
# logarithmic axis.
stopifnot(min(x, na.rm = TRUE) > 0)
x <- log10(x)
major <- 10 ^ seq(
floor(min(x, na.rm = TRUE)),
ceiling(max(x, na.rm = TRUE))
)
n_between <- length(major) - 1
minor <- integer(n_between * 8)
for (i in 1:n_between) {
a <- major[i]
b <- major[i + 1]
minor[seq(i * 8 - 7, i * 8)] <- seq(a + a, b - a, a)
}
major_tick_labels <- formatC(major, format = "fg")
major_tick_labels[1] <- "0"
list(
major = major,
minor = minor,
major_labels = major_tick_labels
)
}
...@@ -17,9 +17,9 @@ ...@@ -17,9 +17,9 @@
# along with this program. If not, see <https://www.gnu.org/licenses/>. # along with this program. If not, see <https://www.gnu.org/licenses/>.
#' Effect Concentrations #' Lethal Concentrations
#' #'
#' Estimate the concentration to reach a certain effect or stress level relative #' Estimate the concentration to reach a certain mortality relative
#' to the control. #' to the control.
#' #'
#' If the response level occurs multiple times because of hormesis, which may #' If the response level occurs multiple times because of hormesis, which may
...@@ -27,7 +27,7 @@ ...@@ -27,7 +27,7 @@
#' smallest concentration is returned. #' smallest concentration is returned.
#' #'
#' This function only makes sense for curves which generally go down with #' This function only makes sense for curves which generally go down with
#' increasing concentration, i.e. all \code{effect_*} curves and also #' increasing concentration, i.e. all \code{survival_*} curves and also
#' \code{sys_tox} and \code{sys_tox_env}. Others are untested and may give #' \code{sys_tox} and \code{sys_tox_env}. Others are untested and may give
#' unexpected results, if any. #' unexpected results, if any.
#' #'
...@@ -36,12 +36,13 @@ ...@@ -36,12 +36,13 @@
#' data frame with a "concentration" column and a \code{response_name} column. #' data frame with a "concentration" column and a \code{response_name} column.
#' In the case of the data frame the first row is assumed to be the control. #' In the case of the data frame the first row is assumed to be the control.
#' See the examples. #' See the examples.
#' @param response_name The name of the effect or stress for which you want to #' @param response_name The name of the survival or stress for which you want to
#' calculate the EC. Must be one of \code{colnames(model$curves)}. #' calculate the LC.
#' @param response_level The desired response level as a percentage between 0 #' @param response_level The desired response level as a percentage between 0
#' and 100. For example with the value 10 the function will return the EC10, #' and 100. For example with the value 10 the function will return the LC10,
#' i.e. the concentration where the response falls below 90 \% of the control #' i.e. the concentration where the response falls below 90 \% of the control
#' response. #' response. In other words: where a mortality of 10 \% relative to the control
#' is reached.
#' @param reference The reference value of the response, usually the control at #' @param reference The reference value of the response, usually the control at
#' concentration 0. This argument is optional. If it is missing, the first #' concentration 0. This argument is optional. If it is missing, the first
#' value of the response is used as control. This value determines what number #' value of the response is used as control. This value determines what number
...@@ -49,49 +50,49 @@ ...@@ -49,49 +50,49 @@
#' response_level is 10 and reference is 45, then the function returns the #' response_level is 10 and reference is 45, then the function returns the
#' concentration where the curve is reduced by 10\% relative to 45 = 40.5. #' concentration where the curve is reduced by 10\% relative to 45 = 40.5.
#' @param warn Logical. Should the function emit a warning if the calculation of #' @param warn Logical. Should the function emit a warning if the calculation of
#' the effect concentration is not possible? #' the lethal concentration is not possible?
#' #'
#' @return A list containing the effect concentration and the corresponding #' @return A list containing the lethal concentration and the corresponding
#' effect. The effect will be \code{NA} if its calculation is impossible from #' survival. The survival will be \code{NA} if its calculation is impossible
#' the supplied data. #' using the supplied data.
#' #'
#' @examples # Calculate the EC10, the concentration where the effect falls #' @examples # Calculate the LC10, the concentration where the survival falls
#' # below 90% of the effect in the control. #' # below 90% of the survival in the control.
#' #'
#' model <- ecxsys( #' model <- ecxsys(
#' concentration = c(0, 0.05, 0.5, 5, 30), #' concentration = c(0, 0.05, 0.5, 5, 30),
#' hormesis_concentration = 0.5, #' hormesis_concentration = 0.5,
#' effect_tox_observed = c(90, 81, 92, 28, 0) #' survival_tox_observed = c(90, 81, 92, 28, 0)
#' ) #' )
#' #'
#' # using the ecxsys() output or the curves therein directly: #' # using the ecxsys() output or the curves therein directly:
#' ec(model, "effect_tox_sys", 10) #' lc(model, "survival_tox_sys", 10)
#' ec(model$curves, "effect_tox_sys", 10) #' lc(model$curves, "survival_tox_sys", 10)
#' #'
#' # using the output of predict_ecxsys() with custom concentrations: #' # using the output of predict_ecxsys() with custom concentrations:
#' conc <- 10^seq(-9, 1, length.out = 1000) #' conc <- 10^seq(-9, 1, length.out = 1000)
#' curves <- predict_ecxsys(model, conc) #' curves <- predict_ecxsys(model, conc)
#' ec(curves, "effect_tox_sys", 10) #' lc(curves, "survival_tox_sys", 10)
#' #'
#' # using a custom data frame: #' # using a custom data frame:
#' df_custom <- data.frame( #' df_custom <- data.frame(
#' concentration = curves$concentration, #' concentration = curves$concentration,
#' foo = curves$effect_tox_sys #' foo = curves$survival_tox_sys
#' ) #' )
#' ec(df_custom, "foo", 10) #' lc(df_custom, "foo", 10)
#' #'
#' # Calculate the EC50 relative to an effect of 100 #' # Calculate the LC50 relative to an survival of 100
#' # instead of relative to the control: #' # instead of relative to the control:
#' ec(model, "effect_tox_sys", 50, reference = 100) #' lc(model, "survival_tox_sys", 50, reference = 100)
#' #'
#' @export #' @export
ec <- function(model, lc <- function(model,
response_name, response_name,
response_level, response_level,
reference, reference,
warn = TRUE) { warn = TRUE) {
if (inherits(model, "drc")) { if (inherits(model, "drc")) {
stop("Please use drc::ED for drc objects.") stop("Please use drc::ED() for drc objects.")
} }
stopifnot( stopifnot(
...@@ -122,7 +123,7 @@ ec <- function(model, ...@@ -122,7 +123,7 @@ ec <- function(model,
} else { } else {
stopifnot(reference > 0) stopifnot(reference > 0)
if (inherits(model, "ecxsys")) { if (inherits(model, "ecxsys")) {
stopifnot(reference <= model$args$effect_max) stopifnot(reference <= model$args$survival_max)
} }
reference_value <- reference reference_value <- reference
} }
...@@ -130,15 +131,15 @@ ec <- function(model, ...@@ -130,15 +131,15 @@ ec <- function(model,
if (reference_value == 0) { if (reference_value == 0) {
# May happen with horizontal Sys curves. # May happen with horizontal Sys curves.
if (warn) { if (warn) {
warning("Reference value is zero, calculation of EC is impossible.") warning("Reference value is zero, calculation of LC is impossible.")
} }
return(list(effect = 0, concentration = NA)) return(list(response = 0, concentration = NA))
} }
if (response_level == 0) { if (response_level == 0) {
return(list(effect = reference_value, concentration = 0)) return(list(response = reference_value, concentration = 0))
} else if (response_level == 100) { } else if (response_level == 100) {
return(list(effect = 0, concentration = Inf)) return(list(response = 0, concentration = Inf))
} }
response_value <- (1 - response_level / 100) * reference_value response_value <- (1 - response_level / 100) * reference_value
...@@ -151,12 +152,12 @@ ec <- function(model, ...@@ -151,12 +152,12 @@ ec <- function(model,
response_name, response_name,
"' are smaller than ", "' are smaller than ",
100 - response_level, 100 - response_level,
"% of the reference value, which makes determining the EC", "% of the reference value, which makes determining the LC",
response_level, response_level,
" impossible." " impossible."
) )
} }
return(list(effect = response_value, concentration = NA)) return(list(response = response_value, concentration = NA))
} }
if (!any(response_smaller)) { if (!any(response_smaller)) {
if (warn) { if (warn) {
...@@ -165,16 +166,16 @@ ec <- function(model, ...@@ -165,16 +166,16 @@ ec <- function(model,
response_name, response_name,
"' does not fall below ", "' does not fall below ",
100 - response_level, 100 - response_level,
"% of the reference, which makes determining the EC", "% of the reference, which makes determining the LC",
response_level, response_level,
" impossible. You could try using predict_ecxsys() to predict ", " impossible. You could try using predict_ecxsys() to predict ",
"more values in a wider concentration range." "more values in a wider concentration range."
) )
} }
return(list(effect = response_value, concentration = NA)) return(list(response = response_value, concentration = NA))
} }
if (response[1] == response_value) { if (response[1] == response_value) {
return(list(effect = response_value, concentration = 0)) return(list(response = response_value, concentration = 0))
} }
below <- which(response < response_value)[1] below <- which(response < response_value)[1]
...@@ -182,8 +183,8 @@ ec <- function(model, ...@@ -182,8 +183,8 @@ ec <- function(model,
# linear interpolation # linear interpolation
dist <- (response_value - response[below]) / (response[above] - response[below]) dist <- (response_value - response[below]) / (response[above] - response[below])
effect_concentration <- dist * survival_concentration <- dist *
(concentration[above] - concentration[below]) + concentration[below] (concentration[above] - concentration[below]) + concentration[below]
list(effect = response_value, concentration = effect_concentration) list(response = response_value, concentration = survival_concentration)
} }
# Copyright (C) 2020 Helmholtz-Zentrum fuer Umweltforschung GmbH - UFZ
# See file inst/COPYRIGHTS for details.
#
# This file is part of the R package stressaddition.
#
# stressaddition is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
#' Logarithmic axis tick marks
#'
#' Calculate the positions and labels of major and minor tick marks for a base
#' 10 logarithmic axis.
#'
#' @param x A vector of axis values. Can be arbitrarily long but only the
#' minimum and maximum are necessary.
#' @param label_zero Whether or not to replace the smallest major label with
#' "0". This defaults to \code{TRUE} and is useful for some types of plots
#' used to display concentration-response data where the leftmost data point
#' represents the control.
#'
#' @return A list with the positions and labels of the major and minor tick
#' marks. The labels are formatted without trailing zeros using
#' \code{formatC(labels, format = "fg")}.
#'
#' @examples
#' x <- c(0.01, 0.2, 3, 10, 50)
#' plot(x, c(5, 4, 2.5, 1, 0), xaxt = "n", log = "x")
#' ticks <- log10_ticks(x)
#' axis(1, at = ticks$major, labels = ticks$major_labels)
#' axis(1, at = ticks$minor, labels = FALSE, tcl = -0.25)
#'
#' @export
log10_ticks <- function(x, label_zero = TRUE) {
stopifnot(min(x, na.rm = TRUE) > 0)
x <- log10(x)
major <- seq(
floor(min(x, na.rm = TRUE)),
ceiling(max(x, na.rm = TRUE))
)
major <- 10 ^ major
n_between <- length(major) - 1
minor <- integer(n_between * 8)
for (i in 1:n_between) {
a <- major[i]
b <- major[i + 1]
minor[seq(i * 8 - 7, i * 8)] <- seq(a + a, b - a, a)
}
major_labels <- formatC(major, format = "fg")
if (label_zero) {
major_labels[1] <- "0"
}
minor_labels <- formatC(minor, format = "fg")
list(
major = major,
minor = minor,
major_labels = major_labels,
minor_labels = minor_labels
)
}
...@@ -17,10 +17,11 @@ ...@@ -17,10 +17,11 @@
# along with this program. If not, see <https://www.gnu.org/licenses/>. # along with this program. If not, see <https://www.gnu.org/licenses/>.
#' Predict the effect of a mixture of two toxicants #' Predict the survival of binary toxicant mixtures
#' #'
#' Given the ecxsys models of two toxicants this method predicts the effects of #' The Multi-TOX model predicts the effects of binary toxicant mixtures based on
#' different mixtures of both. #' three-phasic concentration-response relationships. See the publication for
#' details.
#' #'
#' The predictions are symmetric, i.e. it does not matter which of the toxicant #' 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 #' models is 1 or 2 as long as the concentration arguments are supplied in the
...@@ -34,27 +35,28 @@ ...@@ -34,27 +35,28 @@
#' the mixture. Both vectors must either be the same length or the longer #' the mixture. Both vectors must either be the same length or the longer
#' length must be a multiple of the shorter length. That's because the shorter #' length must be a multiple of the shorter length. That's because the shorter
#' concentration vector gets recycled to the length of the longer one. #' concentration vector gets recycled to the length of the longer one.
#' @param proportion_ca The proportion of concentration addition in the #' @param sa_contribution The proportion of stress addition contributing to the
#' calculation of the toxicant stress of the mixture. Must be between 0 and 1. #' calculation of the toxicant stress in the mixture. Must be between 0 and 1
#' @param effect_max Controls the scaling of the result. This represents the #' where 1 stands for 100 \% stress addition.
#' maximum value the effect could possibly reach. For survival data in percent #' @param survival_max Controls the scaling of the result. This represents the
#' this should be 100 (the default). #' maximum value the survival could possibly reach. For survival data in
#' percent this should be 100 (the default).
#' #'
#' @return A data frame with columns of the supplied concentrations and the #' @return A data frame with columns of the supplied concentrations and the
#' corresponding mixture effects and stresses. #' corresponding mixture survival and stresses.
#' #'
#' @examples toxicant_a <- ecxsys( #' @examples toxicant_a <- ecxsys(
#' concentration = c(0, 0.05, 0.5, 5, 30), #' concentration = c(0, 0.05, 0.5, 5, 30),
#' hormesis_concentration = 0.5, #' hormesis_concentration = 0.5,
#' effect_tox_observed = c(90, 81, 92, 28, 0), #' survival_tox_observed = c(90, 81, 92, 28, 0),
#' ) #' )
#' toxicant_b <- ecxsys( #' toxicant_b <- ecxsys(
#' concentration = c(0, 0.1, 1, 10, 100, 1000), #' concentration = c(0, 0.1, 1, 10, 100, 1000),
#' hormesis_concentration = 10, #' hormesis_concentration = 10,
#' effect_tox_observed = c(26, 25, 24, 27, 5, 0), #' survival_tox_observed = c(26, 25, 24, 27, 5, 0),
#' effect_max = 30 #' survival_max = 30
#' ) #' )
#' predict_mixture( #' multi_tox(
#' toxicant_a , #' toxicant_a ,
#' toxicant_b , #' toxicant_b ,
#' c(0, 0.02, 0.2, 2, 20), #' c(0, 0.02, 0.2, 2, 20),
...@@ -64,18 +66,18 @@ ...@@ -64,18 +66,18 @@
#' # Example of symmetric prediction: #' # Example of symmetric prediction:
#' conc_a <- c(0, 0.03, 0.3, 3) #' conc_a <- c(0, 0.03, 0.3, 3)
#' conc_b <- 5.5 #' conc_b <- 5.5
#' prop_ca <- 0.75 #' sa_contrib <- 0.75
#' mix_a <- predict_mixture(toxicant_a , toxicant_b , conc_a, conc_b, prop_ca) #' mix_a <- multi_tox(toxicant_a , toxicant_b , conc_a, conc_b, sa_contrib)
#' mix_b <- predict_mixture(toxicant_b , toxicant_a , conc_b, conc_a, prop_ca) #' mix_b <- multi_tox(toxicant_b , toxicant_a , conc_b, conc_a, sa_contrib)
#' identical(mix_a$effect, mix_b$effect) #' identical(mix_a$survival, mix_b$survival)
#' #'
#' @export #' @export
predict_mixture <- function(model_a, multi_tox <- function(model_a,
model_b, model_b,
concentration_a, concentration_a,
concentration_b, concentration_b,
proportion_ca = 0.5, sa_contribution = 0.5,
effect_max = 100) { survival_max = 100) {
stopifnot( stopifnot(
inherits(model_a, "ecxsys"), inherits(model_a, "ecxsys"),
inherits(model_b, "ecxsys"), inherits(model_b, "ecxsys"),
...@@ -85,8 +87,8 @@ predict_mixture <- function(model_a, ...@@ -85,8 +87,8 @@ predict_mixture <- function(model_a,
length(concentration_b) > 0, length(concentration_b) > 0,
all(!is.na(concentration_a)), all(!is.na(concentration_a)),
all(!is.na(concentration_b)), all(!is.na(concentration_b)),
proportion_ca >= 0, sa_contribution >= 0,
proportion_ca <= 1, sa_contribution <= 1,
model_a$args$p == model_b$args$p, model_a$args$p == model_b$args$p,
model_a$args$q == model_b$args$q model_a$args$q == model_b$args$q
) )
...@@ -95,29 +97,29 @@ predict_mixture <- function(model_a, ...@@ -95,29 +97,29 @@ predict_mixture <- function(model_a,
predicted_model_b <- predict_ecxsys(model_b, concentration_b) predicted_model_b <- predict_ecxsys(model_b, concentration_b)
# tox stress ---------------------------------------------------------- # tox stress ----------------------------------------------------------
stress_tox_sam <- predicted_model_a$stress_tox + predicted_model_b$stress_tox stress_tox_sa <- predicted_model_a$stress_tox + predicted_model_b$stress_tox
# Convert the model_b concentration into an equivalent model_a concentration # Convert the model_b concentration into an equivalent model_a concentration
# and vice versa. # and vice versa.
concentration_b_equivalent <- W1.2_inverse( concentration_b_equivalent <- W1.2_inverse(
model_a$effect_tox_mod, model_a$survival_tox_mod,
predicted_model_b$effect_tox / model_b$args$effect_max predicted_model_b$survival_tox / model_b$args$survival_max
) )
effect_tox_ca_a <- predict( survival_tox_ca_a <- predict(
model_a$effect_tox_mod, model_a$survival_tox_mod,
data.frame(concentration = concentration_a + concentration_b_equivalent) data.frame(concentration = concentration_a + concentration_b_equivalent)
) )
stress_tox_ca_a <- effect_to_stress(effect_tox_ca_a) stress_tox_ca_a <- survival_to_stress(survival_tox_ca_a)
concentration_a_equivalent <- W1.2_inverse( concentration_a_equivalent <- W1.2_inverse(
model_b$effect_tox_mod, model_b$survival_tox_mod,
predicted_model_a$effect_tox / model_a$args$effect_max predicted_model_a$survival_tox / model_a$args$survival_max
) )
effect_tox_ca_b <- predict( survival_tox_ca_b <- predict(
model_b$effect_tox_mod, model_b$survival_tox_mod,
data.frame(concentration = concentration_b + concentration_a_equivalent) data.frame(concentration = concentration_b + concentration_a_equivalent)
) )
stress_tox_ca_b <- effect_to_stress(effect_tox_ca_b) stress_tox_ca_b <- survival_to_stress(survival_tox_ca_b)
stress_tox_ca <- (stress_tox_ca_a + stress_tox_ca_b) / 2 stress_tox_ca <- (stress_tox_ca_a + stress_tox_ca_b) / 2
...@@ -127,16 +129,16 @@ predict_mixture <- function(model_a, ...@@ -127,16 +129,16 @@ predict_mixture <- function(model_a,
sys_total <- (sys_a + sys_b) / 2 sys_total <- (sys_a + sys_b) / 2
# combined stress and result ------------------------------------------ # combined stress and result ------------------------------------------
proportion_sam <- 1 - proportion_ca ca_contribution <- 1 - sa_contribution
stress_tox_total <- stress_tox_ca * proportion_ca + stress_tox_sam * proportion_sam stress_tox_total <- stress_tox_sa * sa_contribution + stress_tox_ca * ca_contribution
stress_total <- stress_tox_total + sys_total stress_total <- stress_tox_total + sys_total
effect <- stress_to_effect(stress_total) * effect_max survival <- stress_to_survival(stress_total) * survival_max
data.frame( data.frame(
concentration_a = concentration_a, concentration_a = concentration_a,
concentration_b = concentration_b, concentration_b = concentration_b,
effect = effect, survival = survival,
stress_tox_sam = stress_tox_sam, stress_tox_sa = stress_tox_sa,
stress_tox_ca = stress_tox_ca, stress_tox_ca = stress_tox_ca,
stress_tox = stress_tox_total, stress_tox = stress_tox_total,
sys = sys_total, sys = sys_total,
......
...@@ -22,18 +22,18 @@ ...@@ -22,18 +22,18 @@
#' Plot the results of the ECx-SyS model #' Plot the results of the ECx-SyS model
#' #'
#' Plot the observed and modeled effects and stresses. #' Plot the observed and modeled survivals and stresses.
#' #'
#' @name plot_ecxsys #' @name plot_ecxsys
#' #'
#' @param model The list returned from \code{\link{ecxsys}}. #' @param model The list returned from \code{\link{ecxsys}}.
#' @param which A vector of names to plot. Allowed are the column names of the #' @param which A vector of names to plot. Allowed are the column names of the
#' \code{model$curves} data frame. There is also \code{"effect_tox_observed"} #' \code{model$curves} data frame. There is also
#' and \code{"effect_tox_env_observed"} for the observed effects and #' \code{"survival_tox_observed"} and \code{"survival_tox_env_observed"} for
#' \code{"sys_tox_observed"} and \code{"sys_tox_env_observed"} for the #' the observed survival and \code{"sys_tox_observed"} and
#' observed Sys. The default \code{NA} only plots the most important curves. #' \code{"sys_tox_env_observed"} for the observed Sys. The default \code{NA}
#' Use \code{which = "all"} to display all curves. An empty vector or #' only plots the most important curves. Use \code{which = "all"} to display
#' \code{NULL} creates just the axes. #' all curves. An empty vector or \code{NULL} creates just the axes.
#' #'
#' @param show_legend Should the plot include a legend? Defaults to \code{FALSE} #' @param show_legend Should the plot include a legend? Defaults to \code{FALSE}
#' because it may cover some parts of the plot depending on the plot size and #' because it may cover some parts of the plot depending on the plot size and
...@@ -43,21 +43,21 @@ ...@@ -43,21 +43,21 @@
#' @examples model <- ecxsys( #' @examples model <- ecxsys(
#' concentration = c(0, 0.05, 0.5, 5, 30), #' concentration = c(0, 0.05, 0.5, 5, 30),
#' hormesis_concentration = 0.5, #' hormesis_concentration = 0.5,
#' effect_tox_observed = c(90, 81, 92, 28, 0), #' survival_tox_observed = c(90, 81, 92, 28, 0),
#' effect_tox_env_observed = c(29, 27, 33, 5, 0) #' survival_tox_env_observed = c(29, 27, 33, 5, 0)
#' ) #' )
#' plot_effect(model, show_legend = TRUE) #' plot_survival(model, show_legend = TRUE)
#' plot_stress(model, show_legend = TRUE) #' plot_stress(model, show_legend = TRUE)
#' #'
#' # Plot all curves: #' # Plot all curves:
#' plot_effect(model, which = "all") #' plot_survival(model, which = "all")
#' plot_stress(model, which = "all") #' plot_stress(model, which = "all")
#' #'
#' # Plot only some selected curves: #' # Plot only some selected curves:
#' plot_effect(model, which = c("effect_tox_sys", "effect_tox_env_sys")) #' plot_survival(model, which = c("survival_tox_sys", "survival_tox_env_sys"))
#' plot_stress(model, which = c("sys_tox", "sys_tox_env")) #' plot_stress(model, which = c("sys_tox", "sys_tox_env"))
#' #'
#' # Plot only the observed values: #' # Plot only the observed values:
#' plot_effect(model, which = c("effect_tox_observed", "effect_tox_env_observed")) #' plot_survival(model, which = c("survival_tox_observed", "survival_tox_env_observed"))
#' plot_stress(model, which = c("sys_tox_observed", "sys_tox_env_observed")) #' plot_stress(model, which = c("sys_tox_observed", "sys_tox_env_observed"))
NULL NULL
...@@ -44,7 +44,7 @@ plot_stress <- function(model, ...@@ -44,7 +44,7 @@ plot_stress <- function(model,
which <- valid_names which <- valid_names
} else if (!model$with_env && any(grepl("env", which, fixed = TRUE))) { } else if (!model$with_env && any(grepl("env", which, fixed = TRUE))) {
warning("'which' contains names with 'env' but the model was built ", warning("'which' contains names with 'env' but the model was built ",
"without environmental effects.") "without environmental stress.")
which <- which[which %in% valid_names] which <- which[which %in% valid_names]
} else if (any(!which %in% valid_names)) { } else if (any(!which %in% valid_names)) {
warning("Argument 'which' contains invalid names.") warning("Argument 'which' contains invalid names.")
...@@ -52,14 +52,25 @@ plot_stress <- function(model, ...@@ -52,14 +52,25 @@ plot_stress <- function(model,
} }
curves <- model$curves curves <- model$curves
log_ticks <- get_log_ticks(curves$concentration_for_plots) ticks <- log10_ticks(curves$concentration_for_plots)
point_concentration <- c( point_concentration <- c(
curves$concentration_for_plots[1], curves$concentration_for_plots[1],
model$args$concentration[-1] model$args$concentration[-1]
) )
curves_w <- curves[, which[!endsWith(which, "observed")]] if (is.null(which)) {
ymax <- if (NCOL(curves_w) == 0) 1 else max(curves_w, 1, na.rm = TRUE) ymax = 1
} else {
which_lines <- which[!endsWith(which, "observed")]
if (length(which_lines) == 0) {
ymax <- 1
} else {
ymax <- max(curves[, which_lines], 1, na.rm = TRUE)
# No need to include the observed stress in the call to max() no
# matter if those are in "which" or not because these vectors are
# clamped to [0, 1] anyway.
}
}
plot( plot(
NA, NA,
...@@ -151,9 +162,9 @@ plot_stress <- function(model, ...@@ -151,9 +162,9 @@ plot_stress <- function(model,
# The setting of col = NA and col.ticks = par("fg") is to prevent ugly line # The setting of col = NA and col.ticks = par("fg") is to prevent ugly line
# thickness issues when plotting as a png with type = "cairo" and at a low # thickness issues when plotting as a png with type = "cairo" and at a low
# resolution. # resolution.
axis(1, at = log_ticks$major, labels = log_ticks$major_labels, axis(1, at = ticks$major, labels = ticks$major_labels,
col = NA, col.ticks = par("fg")) col = NA, col.ticks = par("fg"))
axis(1, at = log_ticks$minor, labels = FALSE, tcl = -0.25, axis(1, at = ticks$minor, labels = FALSE, tcl = -0.25,
col = NA, col.ticks = par("fg")) col = NA, col.ticks = par("fg"))
plotrix::axis.break(1, breakpos = model$axis_break_conc) plotrix::axis.break(1, breakpos = model$axis_break_conc)
axis(2, col = NA, col.ticks = par("fg"), las = 1) axis(2, col = NA, col.ticks = par("fg"), las = 1)
......
...@@ -19,24 +19,24 @@ ...@@ -19,24 +19,24 @@
#' @rdname plot_ecxsys #' @rdname plot_ecxsys
#' @export #' @export
plot_effect <- function(model, plot_survival <- function(model,
which = NA, which = NA,
show_legend = FALSE, show_legend = FALSE,
xlab = "concentration", xlab = "concentration",
ylab = "effect", ylab = "survival",
main = NULL) { main = NULL) {
stopifnot(inherits(model, "ecxsys")) stopifnot(inherits(model, "ecxsys"))
curve_names <- names(model$curves) curve_names <- names(model$curves)
valid_names <- c( valid_names <- c(
curve_names[startsWith(curve_names, "effect")], curve_names[startsWith(curve_names, "survival")],
"effect_tox_observed", "effect_tox_env_observed" # the observed points "survival_tox_observed", "survival_tox_env_observed" # observed points
) )
if (length(which) == 1 && is.na(which)) { if (length(which) == 1 && is.na(which)) {
which <- c("effect_tox", "effect_tox_sys", "effect_tox_observed") which <- c("survival_tox", "survival_tox_sys", "survival_tox_observed")
if (model$with_env) { if (model$with_env) {
which <- c(which, "effect_tox_env", "effect_tox_env_sys", which <- c(which, "survival_tox_env", "survival_tox_env_sys",
"effect_tox_env_observed") "survival_tox_env_observed")
} }
} else if ("all" %in% which) { } else if ("all" %in% which) {
if (length(which) > 1) { if (length(which) > 1) {
...@@ -45,7 +45,7 @@ plot_effect <- function(model, ...@@ -45,7 +45,7 @@ plot_effect <- function(model,
which <- valid_names which <- valid_names
} else if (!model$with_env && any(grepl("env", which, fixed = TRUE))) { } else if (!model$with_env && any(grepl("env", which, fixed = TRUE))) {
warning("'which' contains names with 'env' but the model was built ", warning("'which' contains names with 'env' but the model was built ",
"without environmental effects.") "without environmental stress.")
which <- which[which %in% valid_names] which <- which[which %in% valid_names]
} else if (any(!which %in% valid_names)) { } else if (any(!which %in% valid_names)) {
warning("Argument 'which' contains invalid names.") warning("Argument 'which' contains invalid names.")
...@@ -53,7 +53,7 @@ plot_effect <- function(model, ...@@ -53,7 +53,7 @@ plot_effect <- function(model,
} }
curves <- model$curves curves <- model$curves
log_ticks <- get_log_ticks(curves$concentration_for_plots) ticks <- log10_ticks(curves$concentration_for_plots)
point_concentration <- c( point_concentration <- c(
curves$concentration_for_plots[1], curves$concentration_for_plots[1],
model$args$concentration[-1] model$args$concentration[-1]
...@@ -63,7 +63,7 @@ plot_effect <- function(model, ...@@ -63,7 +63,7 @@ plot_effect <- function(model,
NA, NA,
NA, NA,
xlim = range(curves$concentration_for_plots, na.rm = TRUE), xlim = range(curves$concentration_for_plots, na.rm = TRUE),
ylim = c(0, model$args$effect_max), ylim = c(0, model$args$survival_max),
log = "x", log = "x",
xlab = xlab, xlab = xlab,
ylab = ylab, ylab = ylab,
...@@ -75,65 +75,65 @@ plot_effect <- function(model, ...@@ -75,65 +75,65 @@ plot_effect <- function(model,
# The lines are drawn in this order to ensure that dotted and dashed lines # The lines are drawn in this order to ensure that dotted and dashed lines
# are on top of solid lines for better visibility. # are on top of solid lines for better visibility.
if ("effect_tox_observed" %in% which) { if ("survival_tox_observed" %in% which) {
points( points(
point_concentration, point_concentration,
model$args$effect_tox_observed, model$args$survival_tox_observed,
pch = 16, pch = 16,
col = "blue" col = "blue"
) )
} }
if ("effect_tox_sys" %in% which) { if ("survival_tox_sys" %in% which) {
lines( lines(
curves$concentration_for_plots, curves$concentration_for_plots,
curves$effect_tox_sys, curves$survival_tox_sys,
col = "blue" col = "blue"
) )
} }
if ("effect_tox" %in% which) { if ("survival_tox" %in% which) {
lines( lines(
curves$concentration_for_plots, curves$concentration_for_plots,
curves$effect_tox, curves$survival_tox,
col = "deepskyblue", col = "deepskyblue",
lty = 2 lty = 2
) )
} }
if ("effect_tox_LL5" %in% which) { if ("survival_tox_LL5" %in% which) {
lines( lines(
curves$concentration_for_plots, curves$concentration_for_plots,
curves$effect_tox_LL5, curves$survival_tox_LL5,
col = "darkblue", col = "darkblue",
lty = 3 lty = 3
) )
} }
if (model$with_env) { if (model$with_env) {
if ("effect_tox_env_observed" %in% which) { if ("survival_tox_env_observed" %in% which) {
points( points(
point_concentration, point_concentration,
model$args$effect_tox_env_observed, model$args$survival_tox_env_observed,
pch = 16, pch = 16,
col = "red" col = "red"
) )
} }
if ("effect_tox_env_sys" %in% which) { if ("survival_tox_env_sys" %in% which) {
lines( lines(
curves$concentration_for_plots, curves$concentration_for_plots,
curves$effect_tox_env_sys, curves$survival_tox_env_sys,
col = "red" col = "red"
) )
} }
if ("effect_tox_env" %in% which) { if ("survival_tox_env" %in% which) {
lines( lines(
curves$concentration_for_plots, curves$concentration_for_plots,
curves$effect_tox_env, curves$survival_tox_env,
col = "orange", col = "orange",
lty = 2 lty = 2
) )
} }
if ("effect_tox_env_LL5" %in% which) { if ("survival_tox_env_LL5" %in% which) {
lines( lines(
curves$concentration_for_plots, curves$concentration_for_plots,
curves$effect_tox_env_LL5, curves$survival_tox_env_LL5,
col = "darkred", col = "darkred",
lty = 3 lty = 3
) )
...@@ -143,18 +143,18 @@ plot_effect <- function(model, ...@@ -143,18 +143,18 @@ plot_effect <- function(model,
# The setting of col = NA and col.ticks = par("fg") is to prevent ugly line # The setting of col = NA and col.ticks = par("fg") is to prevent ugly line
# thickness issues when plotting as a png with type = "cairo" and at a low # thickness issues when plotting as a png with type = "cairo" and at a low
# resolution. # resolution.
axis(1, at = log_ticks$major, labels = log_ticks$major_labels, axis(1, at = ticks$major, labels = ticks$major_labels,
col = NA, col.ticks = par("fg")) col = NA, col.ticks = par("fg"))
axis(1, at = log_ticks$minor, labels = FALSE, tcl = -0.25, axis(1, at = ticks$minor, labels = FALSE, tcl = -0.25,
col = NA, col.ticks = par("fg")) col = NA, col.ticks = par("fg"))
plotrix::axis.break(1, breakpos = model$axis_break_conc) plotrix::axis.break(1, breakpos = model$axis_break_conc)
axis(2, col = NA, col.ticks = par("fg"), las = 1) axis(2, col = NA, col.ticks = par("fg"), las = 1)
if (show_legend) { if (show_legend) {
legend_df <- data.frame( legend_df <- data.frame(
name = c("effect_tox_observed", "effect_tox", "effect_tox_sys", name = c("survival_tox_observed", "survival_tox", "survival_tox_sys",
"effect_tox_LL5", "effect_tox_env_observed", "survival_tox_LL5", "survival_tox_env_observed",
"effect_tox_env", "effect_tox_env_sys", "effect_tox_env_LL5"), "survival_tox_env", "survival_tox_env_sys", "survival_tox_env_LL5"),
text = c("tox (observed)", "tox", "tox + sys", "tox (LL5)", text = c("tox (observed)", "tox", "tox + sys", "tox (LL5)",
"tox + env (observed)", "tox + env", "tox + env + sys", "tox + env (observed)", "tox + env", "tox + env + sys",
"tox + env (LL5)"), "tox + env (LL5)"),
......
...@@ -17,9 +17,9 @@ ...@@ -17,9 +17,9 @@
# along with this program. If not, see <https://www.gnu.org/licenses/>. # along with this program. If not, see <https://www.gnu.org/licenses/>.
#' Predict effects and stresses #' Predict survival and stress
#' #'
#' Calculate the effects and stresses of an ECx-SyS model at arbitrary #' Calculate the survivals and stresses of an ECx-SyS model at arbitrary
#' concentrations. #' concentrations.
#' #'
#' @param model An ECx-SyS model as returned by \code{\link{ecxsys}}. #' @param model An ECx-SyS model as returned by \code{\link{ecxsys}}.
...@@ -29,23 +29,23 @@ ...@@ -29,23 +29,23 @@
#' columns: #' columns:
#' \describe{ #' \describe{
#' \item{concentration}{The supplied concentrations.} #' \item{concentration}{The supplied concentrations.}
#' \item{effect_tox_LL5}{The effect predicted by the five-parameter #' \item{survival_tox_LL5}{The survival predicted by the five-parameter
#' log-logistic model derived from the observations under toxicant stress #' log-logistic model derived from the observations under toxicant stress
#' but without environmental stress.} #' but without environmental stress.}
#' \item{effect_tox}{Modeled effect resulting from toxicant stress.} #' \item{survival_tox}{Modeled survival resulting from toxicant stress.}
#' \item{effect_tox_sys}{Modeled effect resulting from toxicant and system #' \item{survival_tox_sys}{Modeled survival resulting from toxicant
#' stress.} #' and system stress.}
#' \item{stress_tox}{The toxicant stress.} #' \item{stress_tox}{The toxicant stress.}
#' \item{sys_tox}{System stress under toxicant stress conditions #' \item{sys_tox}{System stress under toxicant stress conditions
#' without environmental stress.} #' without environmental stress.}
#' \item{stress_tox_sys}{The sum of \code{stress_tox} and #' \item{stress_tox_sys}{The sum of \code{stress_tox} and
#' \code{sys_tox}.} #' \code{sys_tox}.}
#' \item{effect_tox_env_LL5}{The effect predicted by the five-parameter #' \item{survival_tox_env_LL5}{The survival predicted by the five-parameter
#' log-logistic model derived from the observations under toxicant stress #' log-logistic model derived from the observations under toxicant stress
#' with environmental stress.} #' with environmental stress.}
#' \item{effect_tox_env}{Modeled effect resulting from toxicant and #' \item{survival_tox_env}{Modeled survival resulting from toxicant and
#' environmental stress.} #' environmental stress.}
#' \item{effect_tox_env_sys}{Modeled effect resulting from toxicant, #' \item{survival_tox_env_sys}{Modeled survival resulting from toxicant,
#' environmental and system stress.} #' environmental and system stress.}
#' \item{stress_env}{Environmental stress.} #' \item{stress_env}{Environmental stress.}
#' \item{stress_tox_env}{The sum of toxicant and environmental stress.} #' \item{stress_tox_env}{The sum of toxicant and environmental stress.}
...@@ -58,8 +58,8 @@ ...@@ -58,8 +58,8 @@
#' @examples model <- ecxsys( #' @examples model <- ecxsys(
#' concentration = c(0, 0.05, 0.5, 5, 30), #' concentration = c(0, 0.05, 0.5, 5, 30),
#' hormesis_concentration = 0.5, #' hormesis_concentration = 0.5,
#' effect_tox_observed = c(90, 81, 92, 28, 0), #' survival_tox_observed = c(90, 81, 92, 28, 0),
#' effect_tox_env_observed = c(29, 27, 33, 5, 0) #' survival_tox_env_observed = c(29, 27, 33, 5, 0)
#' ) #' )
#' p <- predict_ecxsys(model, c(0.001, 0.01, 0.1, 1, 10)) #' p <- predict_ecxsys(model, c(0.001, 0.01, 0.1, 1, 10))
#' #'
...@@ -73,30 +73,30 @@ predict_ecxsys <- function(model, concentration) { ...@@ -73,30 +73,30 @@ predict_ecxsys <- function(model, concentration) {
) )
p <- model$args$p p <- model$args$p
q <- model$args$q q <- model$args$q
effect_max <- model$args$effect_max survival_max <- model$args$survival_max
out_df <- data.frame( out_df <- data.frame(
concentration = concentration concentration = concentration
) )
out_df$effect_tox_LL5 <- predict( out_df$survival_tox_LL5 <- predict(
model$effect_tox_LL5_mod, model$survival_tox_LL5_mod,
data.frame(concentration = concentration) data.frame(concentration = concentration)
) * effect_max ) * survival_max
if (model$with_env) { if (model$with_env) {
out_df$effect_tox_env_LL5 <- predict( out_df$survival_tox_env_LL5 <- predict(
model$effect_tox_env_LL5_mod, model$survival_tox_env_LL5_mod,
data.frame(concentration = concentration) data.frame(concentration = concentration)
) * effect_max ) * survival_max
} }
effect_tox <- predict( survival_tox <- predict(
model$effect_tox_mod, model$survival_tox_mod,
data.frame(concentration = concentration) data.frame(concentration = concentration)
) )
out_df$effect_tox <- effect_tox * effect_max out_df$survival_tox <- survival_tox * survival_max
stress_tox <- effect_to_stress(effect_tox, p, q) stress_tox <- survival_to_stress(survival_tox, p, q)
out_df$stress_tox <- stress_tox out_df$stress_tox <- stress_tox
sys_tox <- predict( sys_tox <- predict(
...@@ -108,8 +108,8 @@ predict_ecxsys <- function(model, concentration) { ...@@ -108,8 +108,8 @@ predict_ecxsys <- function(model, concentration) {
stress_tox_sys <- stress_tox + sys_tox stress_tox_sys <- stress_tox + sys_tox
out_df$stress_tox_sys <- stress_tox_sys out_df$stress_tox_sys <- stress_tox_sys
effect_tox_sys <- stress_to_effect(stress_tox_sys, p, q) survival_tox_sys <- stress_to_survival(stress_tox_sys, p, q)
out_df$effect_tox_sys <- effect_tox_sys * effect_max out_df$survival_tox_sys <- survival_tox_sys * survival_max
if (model$with_env) { if (model$with_env) {
out_df$stress_env <- model$stress_env out_df$stress_env <- model$stress_env
...@@ -117,9 +117,9 @@ predict_ecxsys <- function(model, concentration) { ...@@ -117,9 +117,9 @@ predict_ecxsys <- function(model, concentration) {
stress_tox_env <- stress_tox + model$stress_env stress_tox_env <- stress_tox + model$stress_env
out_df$stress_tox_env <- stress_tox_env out_df$stress_tox_env <- stress_tox_env
out_df$effect_tox_env <- stress_to_effect( out_df$survival_tox_env <- stress_to_survival(
stress_tox_env, p, q stress_tox_env, p, q
) * effect_max ) * survival_max
sys_tox_env <- predict( sys_tox_env <- predict(
model$sys_tox_env_mod, model$sys_tox_env_mod,
...@@ -131,10 +131,10 @@ predict_ecxsys <- function(model, concentration) { ...@@ -131,10 +131,10 @@ predict_ecxsys <- function(model, concentration) {
sys_tox_env sys_tox_env
out_df$stress_tox_env_sys <- stress_tox_env_sys out_df$stress_tox_env_sys <- stress_tox_env_sys
effect_tox_env_sys <- stress_to_effect( survival_tox_env_sys <- stress_to_survival(
stress_tox_env_sys, p, q stress_tox_env_sys, p, q
) )
out_df$effect_tox_env_sys <- effect_tox_env_sys * effect_max out_df$survival_tox_env_sys <- survival_tox_env_sys * survival_max
} }
class(out_df) <- c("ecxsys_predicted", class(out_df)) class(out_df) <- c("ecxsys_predicted", class(out_df))
......
...@@ -17,50 +17,46 @@ ...@@ -17,50 +17,46 @@
# along with this program. If not, see <https://www.gnu.org/licenses/>. # along with this program. If not, see <https://www.gnu.org/licenses/>.
#' Convert Between Stress and Effect #' Convert Between Stress and Survival
#' #'
#' Functions to convert effect to general stress or vice versa using the beta #' Functions to convert survival to general stress or vice versa using the beta
#' distribution. #' distribution.
#' #'
#' These are simple wrappers around the beta distribution function #' These are simple wrappers around the beta distribution function
#' \code{\link[stats:Beta]{pbeta}} and the beta quantile function #' \code{\link[stats:Beta]{pbeta}} and the beta quantile function
#' \code{\link[stats:Beta]{qbeta}}. \code{stress_to_effect} returns #' \code{\link[stats:Beta]{qbeta}}. \code{stress_to_survival} returns \code{1 -
#' \code{1 - pbeta(stress, p, q)}. \code{effect_to_stress} first clamps the #' pbeta(stress, p, q)}. \code{survival_to_stress} first clamps the survival to
#' effect to the interval [0, 1] and then returns #' the interval [0, 1] and then returns \code{qbeta(1 - survival, p, q)}.
#' \code{qbeta(1 - effect, p, q)}.
#'
#' "Effect" is a response which is negatively correlated with stress. For
#' example increasing toxicant concentration (stress) reduces survival (effect).
#' #'
#' @name Stressconversion #' @name Stressconversion
#' #'
#' @param effect One or more effect values to convert to general stress. Should #' @param survival One or more survival values to convert to general stress.
#' be a value between 0 and 1. Smaller or bigger values are treated as 0 or 1 #' Should be a value between 0 and 1. Smaller or bigger values are treated as
#' respectively. #' 0 or 1 respectively.
#' @param stress One or more stress values to convert to effect. #' @param stress One or more stress values to convert to survival.
#' @param p,q The shape parameters of the \code{\link[stats:Beta]{beta}} #' @param p,q The shape parameters of the \code{\link[stats:Beta]{beta}}
#' distribution. Default is 3.2. #' distribution. Default is 3.2.
#' #'
#' @examples #' @examples
#' stress <- 0.3 #' stress <- 0.3
#' effect <- stress_to_effect(stress) #' survival <- stress_to_survival(stress)
#' effect_to_stress(effect) #' survival_to_stress(survival)
#' #'
NULL NULL
#' @rdname Stressconversion #' @rdname Stressconversion
#' @export #' @export
effect_to_stress <- function(effect, p = 3.2, q = 3.2) { survival_to_stress <- function(survival, p = 3.2, q = 3.2) {
stopifnot(p >= 0, q >= 0) stopifnot(p >= 0, q >= 0)
effect <- clamp(effect) survival <- clamp(survival)
qbeta(1 - effect, p, q) qbeta(1 - survival, p, q)
} }
#' @rdname Stressconversion #' @rdname Stressconversion
#' @export #' @export
stress_to_effect <- function(stress, p = 3.2, q = 3.2) { stress_to_survival <- function(stress, p = 3.2, q = 3.2) {
stopifnot(p >= 0, q >= 0) stopifnot(p >= 0, q >= 0)
1 - pbeta(stress, p, q) 1 - pbeta(stress, p, q)
} }
...@@ -31,9 +31,8 @@ ...@@ -31,9 +31,8 @@
#' equations. #' equations.
#' #'
#' In the paper the model is introduced in the context of survival #' In the paper the model is introduced in the context of survival
#' experiments. However, it can also be used to model other concentration #' experiments. However, we conjecture that it can also be used to model
#' dependent responses. For this reason the more general term "effect" instead #' other concentration dependent responses.
#' of "survival" is used throughout the package.
#' #'
#' This package is not on CRAN. It is hosted on the UFZ GitLab server. Visit #' This package is not on CRAN. It is hosted on the UFZ GitLab server. Visit
#' the repository linked below for the newest version. See the readme file #' the repository linked below for the newest version. See the readme file
......
# stressaddition # stressaddition
This is the R implementation of the tri-phasic concentration-response model introduced in This is the R implementation of the tri-phasic concentration-response model introduced in
[Liess, M., Henz, S. & Knillmann, S. Predicting low-concentration effects of pesticides. Sci Rep 9, 15248 (2019).](https://doi.org/10.1038/s41598-019-51645-4) [Liess, M., Henz, S. & Knillmann, S. Predicting low-concentration effects of pesticides. Sci Rep 9, 15248 (2019)](https://doi.org/10.1038/s41598-019-51645-4). It allows modeling of ecotoxicological experiments where the response shows signs of a hormesis effect.
It allows modeling of ecotoxicological experiments where the response shows signs of a hormesis effect.
The EC<sub>x-SyS</sub> and Multi-TOX models from this package are also available as part of the [Indicate app](http://www.systemecology.eu/indicate) which offers a graphical user interface.
## Installation ## Installation
Stressaddition is not on CRAN. You can install the most recent stable version from GitLab using the devtools package: Stressaddition is not yet on CRAN. You can install the most recent stable version from GitLab using the remotes package:
``` r ``` r
# install.packages("devtools") install.packages("remotes")
devtools::install_gitlab("oekotox/stressaddition", host = "git.ufz.de") remotes::install_gitlab("oekotox/stressaddition", host = "git.ufz.de")
``` ```
Alternatively there are binary and source builds downloadable from the [releases page](https://git.ufz.de/oekotox/stressaddition/-/releases). Alternatively there are binary and source builds of various versions downloadable from the [releases page](https://git.ufz.de/oekotox/stressaddition/-/releases).
## Updating ## Updating
RStudio's integrated package updater won't detect updates in packages installed from GitHub or GitLab. I recommend running RStudio's integrated package updater won't detect updates in packages installed from GitHub or GitLab. I recommend running
```r ```r
devtools::update_packages() remotes::update_packages()
``` ```
in regular intervals to check for updates from those sources. in regular intervals to check for updates from those sources.
...@@ -22,26 +23,25 @@ in regular intervals to check for updates from those sources. ...@@ -22,26 +23,25 @@ in regular intervals to check for updates from those sources.
Please cite this package if you use it in your analysis. See `citation("stressaddition")` for details. Please cite this package if you use it in your analysis. See `citation("stressaddition")` for details.
## Example ## Example
In the paper we use the model in the context of survival experiments. However, it can also be applicable in modeling other concentration or dose dependent responses. For this reason the more general term "effect" instead of "survival" is used throughout the package.
```r ```r
library(stressaddition) library(stressaddition)
model <- ecxsys( model <- ecxsys(
concentration = c(0, 0.05, 0.5, 5, 30), concentration = c(0, 0.05, 0.5, 5, 30),
hormesis_concentration = 0.5, hormesis_concentration = 0.5,
effect_tox_observed = c(90, 81, 92, 28, 0), survival_tox_observed = c(90, 81, 92, 28, 0),
effect_tox_env_observed = c(29, 27, 33, 5, 0) survival_tox_env_observed = c(29, 27, 33, 5, 0)
) )
# Plot the effect and the system stress: # Plot the effect and the system stress:
par(mfrow = c(2, 1)) par(mfrow = c(2, 1))
plot_effect(model) plot_survival(model)
plot_stress(model) plot_stress(model)
# The LC50 of the effect under the influence of toxicant and system tress: # The LC50 under the influence of toxicant and system tress:
ec(model, "effect_tox_sys", 50) lc(model, "survival_tox_sys", 50)
# The LC10 of the effect under the influence of toxicant, system and environmental tress: # The LC10 under the influence of toxicant, system and environmental tress:
ec(model, "effect_tox_env_sys", 10) lc(model, "survival_tox_env_sys", 10)
``` ```
## Copyright and License ## Copyright and License
......
% Generated by roxygen2: do not edit by hand % Generated by roxygen2: do not edit by hand
% Please edit documentation in R/stress_effect_conversion.R % Please edit documentation in R/stress_survival_conversion.R
\name{Stressconversion} \name{Stressconversion}
\alias{Stressconversion} \alias{Stressconversion}
\alias{effect_to_stress} \alias{survival_to_stress}
\alias{stress_to_effect} \alias{stress_to_survival}
\title{Convert Between Stress and Effect} \title{Convert Between Stress and Survival}
\usage{ \usage{
effect_to_stress(effect, p = 3.2, q = 3.2) survival_to_stress(survival, p = 3.2, q = 3.2)
stress_to_effect(stress, p = 3.2, q = 3.2) stress_to_survival(stress, p = 3.2, q = 3.2)
} }
\arguments{ \arguments{
\item{effect}{One or more effect values to convert to general stress. Should \item{survival}{One or more survival values to convert to general stress.
be a value between 0 and 1. Smaller or bigger values are treated as 0 or 1 Should be a value between 0 and 1. Smaller or bigger values are treated as
respectively.} 0 or 1 respectively.}
\item{p, q}{The shape parameters of the \code{\link[stats:Beta]{beta}} \item{p, q}{The shape parameters of the \code{\link[stats:Beta]{beta}}
distribution. Default is 3.2.} distribution. Default is 3.2.}
\item{stress}{One or more stress values to convert to effect.} \item{stress}{One or more stress values to convert to survival.}
} }
\description{ \description{
Functions to convert effect to general stress or vice versa using the beta Functions to convert survival to general stress or vice versa using the beta
distribution. distribution.
} }
\details{ \details{
These are simple wrappers around the beta distribution function These are simple wrappers around the beta distribution function
\code{\link[stats:Beta]{pbeta}} and the beta quantile function \code{\link[stats:Beta]{pbeta}} and the beta quantile function
\code{\link[stats:Beta]{qbeta}}. \code{stress_to_effect} returns \code{\link[stats:Beta]{qbeta}}. \code{stress_to_survival} returns \code{1 -
\code{1 - pbeta(stress, p, q)}. \code{effect_to_stress} first clamps the pbeta(stress, p, q)}. \code{survival_to_stress} first clamps the survival to
effect to the interval [0, 1] and then returns the interval [0, 1] and then returns \code{qbeta(1 - survival, p, q)}.
\code{qbeta(1 - effect, p, q)}.
"Effect" is a response which is negatively correlated with stress. For
example increasing toxicant concentration (stress) reduces survival (effect).
} }
\examples{ \examples{
stress <- 0.3 stress <- 0.3
effect <- stress_to_effect(stress) survival <- stress_to_survival(stress)
effect_to_stress(effect) survival_to_stress(survival)
} }
...@@ -7,9 +7,9 @@ ...@@ -7,9 +7,9 @@
ecxsys( ecxsys(
concentration, concentration,
hormesis_concentration, hormesis_concentration,
effect_tox_observed, survival_tox_observed,
effect_tox_env_observed = NULL, survival_tox_env_observed = NULL,
effect_max = 100, survival_max = 100,
curves_concentration_max = NULL, curves_concentration_max = NULL,
p = 3.2, p = 3.2,
q = 3.2 q = 3.2
...@@ -20,81 +20,87 @@ ecxsys( ...@@ -20,81 +20,87 @@ ecxsys(
order and the first element must be 0 to indicate the control.} order and the first element must be 0 to indicate the control.}
\item{hormesis_concentration}{The concentration where the hormesis occurs. \item{hormesis_concentration}{The concentration where the hormesis occurs.
This is usually the concentration of the highest effect after the control.} This is usually the concentration of the highest survival after the
control.}
\item{effect_tox_observed}{A vector of effect values observed at the given \item{survival_tox_observed}{A vector of survival values observed at the given
concentrations and in absence of environmental stress. Values must be concentrations and in absence of environmental stress. Values must be
between 0 and \code{effect_max}.} between 0 and \code{survival_max}.}
\item{effect_tox_env_observed}{Effect values observed in the presence of \item{survival_tox_env_observed}{Survival values observed in the presence of
environmental stress. Must be between 0 and \code{effect_max}.} environmental stress. Must be between 0 and \code{survival_max}.}
\item{effect_max}{The maximum value the effect could possibly reach. For \item{survival_max}{The maximum value the survival could possibly reach. For
survival data in percent this should be 100 (the default).} survival data in percent this should be 100 (the default).}
\item{curves_concentration_max}{The maximum concentration of the predicted \item{curves_concentration_max}{The maximum concentration of the predicted
curves. This might be useful if for example your highest observed curves. This might be useful if for example your highest observed
concentration is 30 but you would like to know the predicted values at 100.} concentration is 30 but you would like to know the predicted values on a
scale between 0 and 100.}
\item{p, q}{The shape parameters of the beta distribution. Default is 3.2.} \item{p, q}{The shape parameters of the beta distribution. Default is
\code{p = q = 3.2}.}
} }
\value{ \value{
A list (of class ecxsys) containing many different objects of which A list (of class ecxsys) containing many different objects of which
the most important are listed below. The effect and stress vectors the most important are listed below. The survival and stress vectors
correspond to the provided concentrations. correspond to the provided concentrations.
\describe{ \describe{
\item{effect_tox}{Modeled effect resulting from toxicant stress.} \item{survival_tox}{Modeled survival resulting from toxicant stress.}
\item{effect_tox_sys}{Modeled effect resulting from toxicant and system \item{survival_tox_sys}{Modeled survival resulting from toxicant and
stress.} system stress.}
\item{effect_tox_env}{Modeled effect resulting from toxicant and \item{survival_tox_env}{Modeled survival resulting from toxicant and
environmental stress.} environmental stress.}
\item{effect_tox_env_sys}{Modeled effect resulting from toxicant, \item{survival_tox_env_sys}{Modeled survival resulting from toxicant,
environmental and system stress.} environmental and system stress.}
\item{effect_tox_LL5}{The effect predicted by the five-parameter \item{survival_tox_LL5}{The survival predicted by the five-parameter
log-logistic model derived from the observations under toxicant stress log-logistic model derived from the observations under toxicant stress
but without environmental stress.} but without environmental stress.}
\item{effect_tox_env_LL5}{The effect predicted by the five-parameter \item{survival_tox_env_LL5}{The survival predicted by the five-parameter
log-logistic model derived from the observations under toxicant stress log-logistic model derived from the observations under toxicant stress
with environmental stress.} with environmental stress.}
\item{curves}{A data frame containing effect and stress values as \item{curves}{A data frame containing survival and stress values as
returned by \code{\link{predict_ecxsys}}. The concentrations are returned by \code{\link{predict_ecxsys}}. The concentrations are
regularly spaced on a logarithmic scale in the given concentration range. regularly spaced on a logarithmic scale in the given concentration range.
The control is approximated by the lowest non-control concentration times The control is approximated by the lowest non-control concentration times
1e-7. The additional column \code{concentration_for_plots} is used by the 1e-7. The additional column \code{concentration_for_plots} is used by the
plotting functions of this package to approximate the control and plotting functions of this package to approximate the control and
generate a break in the concentration axis.} generate a nice concentration axis.}
} }
} }
\description{ \description{
The ECx-SyS model for modeling concentration-effect relationships with ECx-SyS is a model for tri-phasic concentration-response relationships where
hormesis. hormetic and subhormetic effects are observed at low concentrations. It
expands the Stress Addition Model (SAM) by introducing system stress (SyS)
which is negatively correlated with toxicant stress. A constant environmental
stress can be included. See the publication for details.
} }
\details{ \details{
It is advised to complete the curve down to zero for optimal prediction. It is advised to complete the curve down to zero for optimal prediction.
Therefore \code{effect_tox_observed} in the highest concentration should be Therefore \code{survival_tox_observed} in the highest concentration should be
at or close to zero. If the model does not fit properly try adding an effect at or close to zero. If the model does not fit properly try adding a survival
of 0 at ten times the maximum observed concentration. of 0 at ten times the maximum observed concentration.
The vectors \code{concentration}, \code{effect_tox_observed} and The vectors \code{concentration}, \code{survival_tox_observed} and
\code{effect_tox_env_observed} (if provided) must be of equal length and \code{survival_tox_env_observed} (if provided) must be of equal length and
sorted by increasing concentration. sorted by increasing concentration.
} }
\examples{ \examples{
model <- ecxsys( model <- ecxsys(
concentration = c(0, 0.05, 0.5, 5, 30), concentration = c(0, 0.05, 0.5, 5, 30),
hormesis_concentration = 0.5, hormesis_concentration = 0.5,
effect_tox_observed = c(90, 81, 92, 28, 0), survival_tox_observed = c(90, 81, 92, 28, 0),
effect_tox_env_observed = c(29, 27, 33, 5, 0) survival_tox_env_observed = c(29, 27, 33, 5, 0)
) )
# Use effect_max if for example the effect is given as the average number of # Use survival_max if for example the survival is given as the average number
# surviving animals and the initial number of animals is 21: # of surviving animals and the initial number of animals is 21:
model <- ecxsys( model <- ecxsys(
concentration = c(0, 0.03, 0.3, 3, 30), concentration = c(0, 0.03, 0.3, 3, 30),
hormesis_concentration = 0.3, hormesis_concentration = 0.3,
effect_tox_observed = c(17, 15.2, 18.8, 7, 0), survival_tox_observed = c(17, 15.2, 18.8, 7, 0),
effect_tox_env_observed = c(4.8, 4.6, 6.4, 0, 0), survival_tox_env_observed = c(4.8, 4.6, 6.4, 0, 0),
effect_max = 21 survival_max = 21
) )
} }
......
% Generated by roxygen2: do not edit by hand % Generated by roxygen2: do not edit by hand
% Please edit documentation in R/ec.R % Please edit documentation in R/lc.R
\name{ec} \name{lc}
\alias{ec} \alias{lc}
\title{Effect Concentrations} \title{Lethal Concentrations}
\usage{ \usage{
ec(model, response_name, response_level, reference, warn = TRUE) lc(model, response_name, response_level, reference, warn = TRUE)
} }
\arguments{ \arguments{
\item{model}{This can be one of three types of objects: Either the output of \item{model}{This can be one of three types of objects: Either the output of
...@@ -13,13 +13,14 @@ data frame with a "concentration" column and a \code{response_name} column. ...@@ -13,13 +13,14 @@ data frame with a "concentration" column and a \code{response_name} column.
In the case of the data frame the first row is assumed to be the control. In the case of the data frame the first row is assumed to be the control.
See the examples.} See the examples.}
\item{response_name}{The name of the effect or stress for which you want to \item{response_name}{The name of the survival or stress for which you want to
calculate the EC. Must be one of \code{colnames(model$curves)}.} calculate the LC.}
\item{response_level}{The desired response level as a percentage between 0 \item{response_level}{The desired response level as a percentage between 0
and 100. For example with the value 10 the function will return the EC10, and 100. For example with the value 10 the function will return the LC10,
i.e. the concentration where the response falls below 90 \% of the control i.e. the concentration where the response falls below 90 \% of the control
response.} response. In other words: where a mortality of 10 \% relative to the control
is reached.}
\item{reference}{The reference value of the response, usually the control at \item{reference}{The reference value of the response, usually the control at
concentration 0. This argument is optional. If it is missing, the first concentration 0. This argument is optional. If it is missing, the first
...@@ -29,15 +30,15 @@ response_level is 10 and reference is 45, then the function returns the ...@@ -29,15 +30,15 @@ response_level is 10 and reference is 45, then the function returns the
concentration where the curve is reduced by 10\% relative to 45 = 40.5.} concentration where the curve is reduced by 10\% relative to 45 = 40.5.}
\item{warn}{Logical. Should the function emit a warning if the calculation of \item{warn}{Logical. Should the function emit a warning if the calculation of
the effect concentration is not possible?} the lethal concentration is not possible?}
} }
\value{ \value{
A list containing the effect concentration and the corresponding A list containing the lethal concentration and the corresponding
effect. The effect will be \code{NA} if its calculation is impossible from survival. The survival will be \code{NA} if its calculation is impossible
the supplied data. using the supplied data.
} }
\description{ \description{
Estimate the concentration to reach a certain effect or stress level relative Estimate the concentration to reach a certain mortality relative
to the control. to the control.
} }
\details{ \details{
...@@ -46,38 +47,38 @@ happen for low values of \code{response_level}, then the occurrence with the ...@@ -46,38 +47,38 @@ happen for low values of \code{response_level}, then the occurrence with the
smallest concentration is returned. smallest concentration is returned.
This function only makes sense for curves which generally go down with This function only makes sense for curves which generally go down with
increasing concentration, i.e. all \code{effect_*} curves and also increasing concentration, i.e. all \code{survival_*} curves and also
\code{sys_tox} and \code{sys_tox_env}. Others are untested and may give \code{sys_tox} and \code{sys_tox_env}. Others are untested and may give
unexpected results, if any. unexpected results, if any.
} }
\examples{ \examples{
# Calculate the EC10, the concentration where the effect falls # Calculate the LC10, the concentration where the survival falls
# below 90\% of the effect in the control. # below 90\% of the survival in the control.
model <- ecxsys( model <- ecxsys(
concentration = c(0, 0.05, 0.5, 5, 30), concentration = c(0, 0.05, 0.5, 5, 30),
hormesis_concentration = 0.5, hormesis_concentration = 0.5,
effect_tox_observed = c(90, 81, 92, 28, 0) survival_tox_observed = c(90, 81, 92, 28, 0)
) )
# using the ecxsys() output or the curves therein directly: # using the ecxsys() output or the curves therein directly:
ec(model, "effect_tox_sys", 10) lc(model, "survival_tox_sys", 10)
ec(model$curves, "effect_tox_sys", 10) lc(model$curves, "survival_tox_sys", 10)
# using the output of predict_ecxsys() with custom concentrations: # using the output of predict_ecxsys() with custom concentrations:
conc <- 10^seq(-9, 1, length.out = 1000) conc <- 10^seq(-9, 1, length.out = 1000)
curves <- predict_ecxsys(model, conc) curves <- predict_ecxsys(model, conc)
ec(curves, "effect_tox_sys", 10) lc(curves, "survival_tox_sys", 10)
# using a custom data frame: # using a custom data frame:
df_custom <- data.frame( df_custom <- data.frame(
concentration = curves$concentration, concentration = curves$concentration,
foo = curves$effect_tox_sys foo = curves$survival_tox_sys
) )
ec(df_custom, "foo", 10) lc(df_custom, "foo", 10)
# Calculate the EC50 relative to an effect of 100 # Calculate the LC50 relative to an survival of 100
# instead of relative to the control: # instead of relative to the control:
ec(model, "effect_tox_sys", 50, reference = 100) lc(model, "survival_tox_sys", 50, reference = 100)
} }
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment