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

Merge branch 'develop' into 'master' 

Improve ec() and axis concentration handling

Closes #35

See merge request !28
parents 3013a08c 1aae72c3
No related branches found
No related tags found
1 merge request!28Improve ec() and axis concentration handling
Pipeline #3310 passed with stage
in 8 minutes and 38 seconds
Hide whitespace changes
Inline Side-by-side
DESCRIPTION
+ 2
− 2
View file @ a1f2a052
Package: stressaddition
Type: Package
Title: Modeling Tri-Phasic Concentration-Response Relationships
Version: 2.5.0
Date: 2020-03-16
Version: 2.6.0
Date: 2020-04-07
Authors@R: c(person("Sebastian",
"Henz",
role = c("aut", "cre"),
......
NEWS.md
+ 8
− 0
View file @ a1f2a052
# stressaddition 2.6.0
* The `curves` data frame in the output of `ecxsys()` now contains a column with the concentrations which are used for the plot functions in this package. This is useful for generating a nicer concentration axis.
* Changes to `ec()`:
* Renamed `response_value` to `effect` in the output list.
* `response_level` of 0 or 100 is now allowed. 0 returns the concentration 0 and 100 returns the concentration `Inf`. Previously this resulted in an error.
* It is now possible to set the reference to a custom value, for example 100.
# stressaddition 2.5.0
* Fixed unintended behaviour in `plot_effect()` and `plot_stress()` where supplying an empty vector caused the four standard curves to show. Now setting `which` to an empty vector or `NULL` shows just the axes. The default value is NA.
......
R/ec.R
+ 86
− 23
View file @ a1f2a052
......@@ -34,6 +34,7 @@
#' @param model This can be one of three types of objects: Either the output of
#' \code{\link{ecxsys}} or the output of \code{\link{predict_ecxsys}} or a
#' 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.
#' See the examples.
#' @param response_name The name of the effect or stress for which you want to
#' calculate the EC. Must be one of \code{colnames(model$curves)}.
......@@ -41,12 +42,21 @@
#' and 100. For example with the value 10 the function will return the EC10,
#' i.e. the concentration where the response falls below 90 \% of the control
#' response.
#' @param reference The reference value of the response, usually the control at
#' 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
#' \code{response_level} actually corresponds to. For example if
#' 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.
#' @param warn Logical. Should the function emit a warning if the calculation of
#' the effect concentration is not possible?
#'
#' @return A list containing the response concentration and the corresponding
#' response value.
#' @return A list containing the effect concentration and the corresponding
#' effect. The effect will be \code{NA} if its calculation is impossible from
#' the supplied data.
#'
#' @examples # Calculate the EC_10, the concentration where the effect falls
#' # below 90 % of the effect in the control.
#' @examples # Calculate the EC10, the concentration where the effect falls
#' # below 90% of the effect in the control.
#'
#' model <- ecxsys(
#' concentration = c(0, 0.05, 0.5, 5, 30),
......@@ -70,8 +80,16 @@
#' )
#' ec(df_custom, "foo", 10)
#'
#' # Calculate the EC50 relative to an effect of 100
#' # instead of relative to the control:
#' ec(model, "effect_tox_sys", 50, reference = 100)
#'
#' @export
ec <- function(model, response_name, response_level) {
ec <- function(model,
response_name,
response_level,
reference,
warn = TRUE) {
if (inherits(model, "drc")) {
stop("Please use drc::ED for drc objects.")
}
......@@ -80,8 +98,8 @@ ec <- function(model, response_name, response_level) {
is.character(response_name),
length(response_name) == 1,
response_name != "concentration",
response_level > 0,
response_level < 100
response_level >= 0,
response_level <= 100
)
if (!inherits(model, c("ecxsys", "ecxsys_predicted"))) {
......@@ -99,28 +117,73 @@ ec <- function(model, response_name, response_level) {
response <- model[, response_name]
}
reference <- response[1]
if (reference == 0) {
stop("Reference value is zero, calculation of EC is impossible.")
if (missing(reference)) {
reference_value <- response[1]
} else {
stopifnot(reference > 0)
if (inherits(model, "ecxsys")) {
stopifnot(reference <= model$args$effect_max)
}
reference_value <- reference
}
response_value <- (1 - response_level / 100) * reference
output <- list(response_value = response_value)
# Get the index of where the response changes from above to below
# response_value:
below <- which(response < response_value)[1]
if (is.na(below)) {
stop("The curve '", response_name, "' does not fall below ",
100 - response_level, "% of the control, which makes ",
"determining the EC", response_level, " impossible.\n You could ",
"try using predict_ecxsys() to predict more values in a wider ",
"concentration range.")
if (reference_value == 0) {
# May happen with horizontal Sys curves.
if (warn) {
warning("Reference value is zero, calculation of EC is impossible.")
}
return(list(effect = 0, concentration = NA))
}
if (response_level == 0) {
return(list(effect = reference_value, concentration = 0))
} else if (response_level == 100) {
return(list(effect = 0, concentration = Inf))
}
response_value <- (1 - response_level / 100) * reference_value
response_smaller <- response < response_value
if (response_smaller[1]) {
if (warn) {
warning(
"The first values of '",
response_name,
"' are smaller than ",
100 - response_level,
"% of the reference value, which makes determining the EC",
response_level,
" impossible."
)
}
return(list(effect = response_value, concentration = NA))
}
if (!any(response_smaller)) {
if (warn) {
warning(
"The curve '",
response_name,
"' does not fall below ",
100 - response_level,
"% of the reference, which makes determining the EC",
response_level,
" impossible. You could try using predict_ecxsys() to predict ",
"more values in a wider concentration range."
)
}
return(list(effect = response_value, concentration = NA))
}
if (response[1] == response_value) {
return(list(effect = response_value, concentration = 0))
}
below <- which(response < response_value)[1]
above <- below - 1
# linear interpolation
dist <- (response_value - response[below]) / (response[above] - response[below])
output$concentration <- dist *
effect_concentration <- dist *
(concentration[above] - concentration[below]) + concentration[below]
output
list(effect = response_value, concentration = effect_concentration)
}
R/ecxsys.R
+ 36
− 4
View file @ a1f2a052
......@@ -77,7 +77,7 @@
#' returned by \code{\link{predict_ecxsys}}. The concentrations are
#' regularly spaced on a logarithmic scale in the given concentration range.
#' The control is approximated by the lowest non-control concentration times
#' 1e-7. The additional column \code{use_for_plotting} 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
#' generate a break in the concentration axis.}
#' }
......@@ -311,9 +311,15 @@ ecxsys <- function(concentration,
length.out = len_curves
)
output$curves <- predict_ecxsys(output, curves_concentration)
output$curves$use_for_plotting <-
curves_concentration < min_conc * conc_adjust_factor * 1.5 |
curves_concentration > min_conc * 1.5
conc_axis <- adjust_plot_concentrations(
curves_concentration,
min_conc,
conc_adjust_factor
)
output$curves$concentration_for_plots <- conc_axis$concentration
output$axis_break_conc <- conc_axis$axis_break_conc
output
}
......@@ -432,3 +438,29 @@ interpolate <- function(x, to_index, n_new, conc = FALSE) {
}
append(x, x_new[-c(1, len)], from_index)
}
adjust_plot_concentrations <- function(concentration,
min_conc,
conc_adjust_factor) {
# Deals with the concentrations which are unnecessary for plotting. This
# means it removes the concentrations in the gap and increases the
# concentrations below the gap.
use_for_plotting <- (
concentration < min_conc * conc_adjust_factor * 1.5 |
concentration > min_conc * 1.5
)
gap_idx <- min(which(!use_for_plotting))
# Add NAs to force breaks in the lines:
axis_break_conc <- concentration[use_for_plotting][gap_idx]
concentration[!use_for_plotting] <- NA
# Shift the small concentrations upwards so the plot has a nice x-axis:
concentration[1:gap_idx] <- concentration[1:gap_idx] / conc_adjust_factor
list(
concentration = concentration,
axis_break_conc = axis_break_conc
)
}
R/helpers.R
+ 0
− 26
View file @ a1f2a052
......@@ -50,29 +50,3 @@ get_log_ticks <- function(x) {
major_labels = major_tick_labels
)
}
adjust_plot_concentrations <- function(model) {
# Helper for the plot functions, not exported.
# Deals with the concentrations which are unnecessary for plotting. This
# means it removes the concentrations in the gap and increases the
# concentrations left of the gap.
curves <- model$curves
gap_idx <- min(which(!curves$use_for_plotting))
# Keep only the values to the left and right of the axis break:
curves <- curves[curves$use_for_plotting, ]
# Add NAs to force breaks in the lines:
axis_break_conc <- curves$concentration[gap_idx]
curves[gap_idx, ] <- NA
# Shift the small concentrations upwards so the plot has a nice x-axis:
curves$concentration[1:gap_idx] <-
curves$concentration[1:gap_idx] / model$conc_adjust_factor
list(
curves = curves,
axis_break_conc = axis_break_conc
)
}
R/plot_effect.R
+ 16
− 14
View file @ a1f2a052
......@@ -52,15 +52,17 @@ plot_effect <- function(model,
which <- which[which %in% valid_names]
}
temp <- adjust_plot_concentrations(model)
curves <- temp$curves
log_ticks <- get_log_ticks(curves$concentration)
concentration <- c(curves$concentration[1], model$args$concentration[-1])
curves <- model$curves
log_ticks <- get_log_ticks(curves$concentration_for_plots)
point_concentration <- c(
curves$concentration_for_plots[1],
model$args$concentration[-1]
)
plot(
NA,
NA,
xlim = range(curves$concentration, na.rm = TRUE),
xlim = range(curves$concentration_for_plots, na.rm = TRUE),
ylim = c(0, model$args$effect_max),
log = "x",
xlab = xlab,
......@@ -75,7 +77,7 @@ plot_effect <- function(model,
# are on top of solid lines for better visibility.
if ("effect_tox_observed" %in% which) {
points(
concentration,
point_concentration,
model$args$effect_tox_observed,
pch = 16,
col = "blue"
......@@ -83,14 +85,14 @@ plot_effect <- function(model,
}
if ("effect_tox_sys" %in% which) {
lines(
curves$concentration,
curves$concentration_for_plots,
curves$effect_tox_sys,
col = "blue"
)
}
if ("effect_tox" %in% which) {
lines(
curves$concentration,
curves$concentration_for_plots,
curves$effect_tox,
col = "deepskyblue",
lty = 2
......@@ -98,7 +100,7 @@ plot_effect <- function(model,
}
if ("effect_tox_LL5" %in% which) {
lines(
curves$concentration,
curves$concentration_for_plots,
curves$effect_tox_LL5,
col = "darkblue",
lty = 3
......@@ -107,7 +109,7 @@ plot_effect <- function(model,
if (model$with_env) {
if ("effect_tox_env_observed" %in% which) {
points(
concentration,
point_concentration,
model$args$effect_tox_env_observed,
pch = 16,
col = "red"
......@@ -115,14 +117,14 @@ plot_effect <- function(model,
}
if ("effect_tox_env_sys" %in% which) {
lines(
curves$concentration,
curves$concentration_for_plots,
curves$effect_tox_env_sys,
col = "red"
)
}
if ("effect_tox_env" %in% which) {
lines(
curves$concentration,
curves$concentration_for_plots,
curves$effect_tox_env,
col = "orange",
lty = 2
......@@ -130,7 +132,7 @@ plot_effect <- function(model,
}
if ("effect_tox_env_LL5" %in% which) {
lines(
curves$concentration,
curves$concentration_for_plots,
curves$effect_tox_env_LL5,
col = "darkred",
lty = 3
......@@ -145,7 +147,7 @@ plot_effect <- function(model,
col = NA, col.ticks = par("fg"))
axis(1, at = log_ticks$minor, labels = FALSE, tcl = -0.25,
col = NA, col.ticks = par("fg"))
plotrix::axis.break(1, breakpos = temp$axis_break_conc)
plotrix::axis.break(1, breakpos = model$axis_break_conc)
axis(2, col = NA, col.ticks = par("fg"), las = 1)
if (show_legend) {
......
R/plot_stress.R
+ 17
− 15
View file @ a1f2a052
......@@ -51,10 +51,12 @@ plot_stress <- function(model,
which <- which[which %in% valid_names]
}
temp <- adjust_plot_concentrations(model)
curves <- temp$curves
log_ticks <- get_log_ticks(curves$concentration)
concentration <- c(curves$concentration[1], model$args$concentration[-1])
curves <- model$curves
log_ticks <- get_log_ticks(curves$concentration_for_plots)
point_concentration <- c(
curves$concentration_for_plots[1],
model$args$concentration[-1]
)
curves_w <- curves[, which[!endsWith(which, "observed")]]
ymax <- if (NCOL(curves_w) == 0) 1 else max(curves_w, 1, na.rm = TRUE)
......@@ -62,7 +64,7 @@ plot_stress <- function(model,
plot(
NA,
NA,
xlim = range(curves$concentration, na.rm = TRUE),
xlim = range(curves$concentration_for_plots, na.rm = TRUE),
ylim = c(0, ymax),
log = "x",
xlab = xlab,
......@@ -77,7 +79,7 @@ plot_stress <- function(model,
# are on top of solid lines for better visibility.
if ("sys_tox_observed" %in% which) {
points(
concentration,
point_concentration,
model$sys_tox_observed,
pch = 16,
col = "steelblue3"
......@@ -85,14 +87,14 @@ plot_stress <- function(model,
}
if ("stress_tox_sys" %in% which) {
lines(
curves$concentration,
curves$concentration_for_plots,
curves$stress_tox_sys,
col = "blue"
)
}
if ("stress_tox" %in% which) {
lines(
curves$concentration,
curves$concentration_for_plots,
curves$stress_tox,
col = "deepskyblue",
lty = 2
......@@ -100,7 +102,7 @@ plot_stress <- function(model,
}
if ("sys_tox" %in% which) {
lines(
curves$concentration,
curves$concentration_for_plots,
curves$sys_tox,
col = "steelblue3"
)
......@@ -108,7 +110,7 @@ plot_stress <- function(model,
if (model$with_env) {
if ("sys_tox_env_observed" %in% which) {
points(
concentration,
point_concentration,
model$sys_tox_env_observed,
pch = 16,
col = "violetred"
......@@ -116,14 +118,14 @@ plot_stress <- function(model,
}
if ("stress_tox_env_sys" %in% which) {
lines(
curves$concentration,
curves$concentration_for_plots,
curves$stress_tox_env_sys,
col = "red"
)
}
if ("stress_env" %in% which) {
lines(
curves$concentration,
curves$concentration_for_plots,
curves$stress_env,
col = "forestgreen",
lty = 3
......@@ -131,7 +133,7 @@ plot_stress <- function(model,
}
if ("stress_tox_env" %in% which) {
lines(
curves$concentration,
curves$concentration_for_plots,
curves$stress_tox_env,
col = "orange",
lty = 2
......@@ -139,7 +141,7 @@ plot_stress <- function(model,
}
if ("sys_tox_env" %in% which) {
lines(
curves$concentration,
curves$concentration_for_plots,
curves$sys_tox_env,
col = "violetred"
)
......@@ -153,7 +155,7 @@ plot_stress <- function(model,
col = NA, col.ticks = par("fg"))
axis(1, at = log_ticks$minor, labels = FALSE, tcl = -0.25,
col = NA, col.ticks = par("fg"))
plotrix::axis.break(1, breakpos = temp$axis_break_conc)
plotrix::axis.break(1, breakpos = model$axis_break_conc)
axis(2, col = NA, col.ticks = par("fg"), las = 1)
if (show_legend) {
......
man/ec.Rd
+ 21
− 5
View file @ a1f2a052
......@@ -4,12 +4,13 @@
\alias{ec}
\title{Effect Concentrations}
\usage{
ec(model, response_name, response_level)
ec(model, response_name, response_level, reference, warn = TRUE)
}
\arguments{
\item{model}{This can be one of three types of objects: Either the output of
\code{\link{ecxsys}} or the output of \code{\link{predict_ecxsys}} or a
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.
See the examples.}
\item{response_name}{The name of the effect or stress for which you want to
......@@ -19,10 +20,21 @@ calculate the EC. Must be one of \code{colnames(model$curves)}.}
and 100. For example with the value 10 the function will return the EC10,
i.e. the concentration where the response falls below 90 \% of the control
response.}
\item{reference}{The reference value of the response, usually the control at
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
\code{response_level} actually corresponds to. For example if
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.}
\item{warn}{Logical. Should the function emit a warning if the calculation of
the effect concentration is not possible?}
}
\value{
A list containing the response concentration and the corresponding
response value.
A list containing the effect concentration and the corresponding
effect. The effect will be \code{NA} if its calculation is impossible from
the supplied data.
}
\description{
Estimate the concentration to reach a certain effect or stress level relative
......@@ -39,8 +51,8 @@ increasing concentration, i.e. all \code{effect_*} curves and also
unexpected results, if any.
}
\examples{
# Calculate the EC_10, the concentration where the effect falls
# below 90 \% of the effect in the control.
# Calculate the EC10, the concentration where the effect falls
# below 90\% of the effect in the control.
model <- ecxsys(
concentration = c(0, 0.05, 0.5, 5, 30),
......@@ -64,4 +76,8 @@ df_custom <- data.frame(
)
ec(df_custom, "foo", 10)
# Calculate the EC50 relative to an effect of 100
# instead of relative to the control:
ec(model, "effect_tox_sys", 50, reference = 100)
}
man/ecxsys.Rd
+ 1
− 1
View file @ a1f2a052
......@@ -60,7 +60,7 @@ A list (of class ecxsys) containing many different objects of which
returned by \code{\link{predict_ecxsys}}. The concentrations are
regularly spaced on a logarithmic scale in the given concentration range.
The control is approximated by the lowest non-control concentration times
1e-7. The additional column \code{use_for_plotting} 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
generate a break in the concentration axis.}
}
......
tests/testthat/test-ec.R
+ 37
− 25
View file @ a1f2a052
......@@ -17,12 +17,16 @@
# along with this program. If not, see <https://www.gnu.org/licenses/>.
# shared model for some of the tests
model <- 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)
)
test_that("all input formats produce identical models", {
model <- ecxsys(
concentration = c(0, 0.05, 0.5, 5, 30),
hormesis_concentration = 0.5,
effect_tox_observed = c(90, 81, 92, 28, 0)
)
# using the ecxsys() output or the curves therein directly:
ec10_a <- ec(model, "effect_tox_sys", 10)
ec10_b <- ec(model$curves, "effect_tox_sys", 10)
......@@ -46,64 +50,72 @@ test_that("all input formats produce identical models", {
test_that("ec values have not changed", {
model <- 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)
)
expect_equal(
ec(model, "effect_tox_sys", 50),
list(response_value = 44.95368, concentration = 3.375735),
list(effect = 44.95368, concentration = 3.375735),
tolerance = 1e-4
)
expect_equal(
ec(model, "effect_tox_sys", 10),
list(response_value = 80.91662, concentration = 1.098648),
list(effect = 80.91662, concentration = 1.098648),
tolerance = 1e-4
)
expect_equal(
ec(model, "effect_tox_sys", 5),
list(response_value = 85.41198, concentration = 0.005502182),
ec(model, "effect_tox", 100/3),
list(effect = 66.66667, concentration = 1.902125),
tolerance = 1e-4
)
expect_equal(
ec(model, "effect_tox_sys", 1),
list(response_value = 89.00828, concentration = 8.53288e-05),
ec(model, "effect_tox_LL5", 42),
list(effect = 52.2, concentration = 2.054426),
tolerance = 1e-4
)
expect_equal(
ec(model, "effect_tox_env_sys", 50),
list(response_value = 14.67725, concentration = 1.299516),
list(effect = 14.67725, concentration = 1.299516),
tolerance = 1e-4
)
expect_equal(
ec(model, "effect_tox_env_sys", 10),
list(response_value = 26.41904, concentration = 0.0008571244),
list(effect = 26.41904, concentration = 0.0008571244),
tolerance = 1e-4
)
expect_equal(
ec(model, "effect_tox_env_sys", 5),
list(response_value = 27.88677, concentration = 0.0001044245),
ec(model, "effect_tox_env", 67.89),
list(effect = 24.51453, concentration = 0.7890892),
tolerance = 1e-4
)
expect_equal(
ec(model, "effect_tox_env_sys", 1),
list(response_value = 29.06095, concentration = 1.388112e-06),
ec(model, "effect_tox_env_LL5", 3.14159),
list(effect = 28.73466, concentration = 0.7267524),
tolerance = 1e-4
)
})
test_that("ec fails when response_level is outside the range of the curve", {
test_that("ec warns when response_level is outside the range of the curve", {
model <- ecxsys(
concentration = c(0, 0.05, 0.5, 5),
hormesis_concentration = 0.5,
effect_tox_observed = c(90, 81, 92, 28)
)
expect_error(
expect_warning(
ec(model, "effect_tox_sys", 90),
"You could try using predict_ecxsys() to predict more values",
fixed = TRUE
)
})
test_that("reference argument works", {
expect_equal(
ec(model, "effect_tox_LL5", 50, reference = 100),
list(effect = 50, concentration = 2.208119),
tolerance = 1e-4
)
expect_equal(
ec(model, "effect_tox_LL5", 50, reference = 75),
list(effect = 37.5, concentration = 3.342715),
tolerance = 1e-4
)
})
tests/testthat/test-ecxsys.R
+ 2
− 2
View file @ a1f2a052
......@@ -161,7 +161,7 @@ test_that("the curves have not changed", {
sys_tox_env = c(0.25443, 0.20496, 0.04445, 0.00000, 0.00000),
stress_tox_env_sys = c(0.61020, 0.63691, 0.59102, 0.85800, 1.33916),
effect_tox_env_sys = c(29.35449, 24.84147, 32.75385, 1.93182, 0),
use_for_plotting = c(TRUE, TRUE, TRUE, TRUE, TRUE)
concentration_for_plots = c(0.0001, 0.03004, 0.30512, 3.02551, 30)
)
class(new_curves) <- class(reference_curves)
expect_equal(new_curves, reference_curves, tolerance = 1e-3)
......@@ -171,7 +171,7 @@ test_that("the curves have not changed", {
test_that("the returned fn works the same way as internally", {
# I don't know why it would fail but it doesn't hurt to test it.
curves <- mod$curves
curves$use_for_plotting <- NULL # remove column "use_for_plotting"
curves$concentration_for_plots <- NULL
expect_identical(curves, predict_ecxsys(mod, curves$concentration))
})
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment