R/ecxsys.R

In stressaddition: Modelling Tri-Phasic Concentration-Response Relationships

# 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/>.


# Note about the setting and resetting of options:
# The drc package changes some of the options which some users may have
# modified. Of particular interest is options("warn") because it is important
# that the user is able to control the visibility of warnings generated by
# ecxsys(). Therefore every time drc::drm() is called this option is cached
# beforehand and reset immediately afterwards. Additionally, all options are
# cached at the beginning of ecxsys() and reset on exit.


#' ECx-SyS
#'
#' ECx-SyS is a model for tri-phasic concentration-response relationships where
#' 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.
#'
#' It is advised to complete the curve down to zero for optimal prediction.
#' 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 a survival
#' of 0 at ten times the maximum observed concentration.
#'
#' The vectors \code{concentration}, \code{survival_tox_observed} and
#' \code{survival_tox_env_observed} (if provided) must be of equal length and
#' sorted by increasing concentration.
#'
#' @param concentration A vector of concentrations. Must be sorted in ascending
#'   order and the first element must be 0 to indicate the control.
#' @param hormesis_concentration The concentration where the hormesis occurs.
#'   This is usually the concentration of the highest survival after the
#'   control.
#' @param survival_tox_observed A vector of survival values observed at the given
#'   concentrations and in absence of environmental stress. Values must be
#'   between 0 and \code{survival_max}.
#' @param survival_tox_env_observed Survival values observed in the presence of
#'   environmental stress. Must be between 0 and \code{survival_max}.
#' @param survival_max The maximum value the survival could possibly reach. For
#'   survival data in percent this should be 100 (the default).
#' @param curves_concentration_max The maximum concentration of the predicted
#'   curves. This might be useful if for example your highest observed
#'   concentration is 30 but you would like to know the predicted values on a
#'   scale between 0 and 100.
#' @param p,q The shape parameters of the beta distribution. Default is
#'   \code{p = q = 3.2}.
#'
#' @return A list (of class ecxsys) containing many different objects of which
#'   the most important are listed below. The survival and stress vectors
#'   correspond to the provided concentrations.
#'   \describe{
#'     \item{survival_tox}{Modeled survival resulting from toxicant stress.}
#'     \item{survival_tox_sys}{Modeled survival resulting from toxicant and
#'     system stress.}
#'     \item{survival_tox_env}{Modeled survival resulting from toxicant and
#'     environmental stress.}
#'     \item{survival_tox_env_sys}{Modeled survival resulting from toxicant,
#'     environmental and system stress.}
#'     \item{survival_tox_LL5}{The survival predicted by the five-parameter
#'     log-logistic model derived from the observations under toxicant stress
#'     but without environmental stress.}
#'     \item{survival_tox_env_LL5}{The survival predicted by the five-parameter
#'     log-logistic model derived from the observations under toxicant stress
#'     with environmental stress.}
#'     \item{curves}{A data frame containing survival and stress values as
#'     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{concentration_for_plots} is used by the
#'     plotting functions of this package to approximate the control and
#'     generate a nice concentration axis.}
#'   }
#'
#' @examples model <- ecxsys(
#'     concentration = c(0, 0.05, 0.5, 5, 30),
#'     hormesis_concentration = 0.5,
#'     survival_tox_observed = c(90, 81, 92, 28, 0),
#'     survival_tox_env_observed = c(29, 27, 33, 5, 0)
#' )
#'
#' # Use survival_max if for example the survival is given as the average number
#' # of surviving animals and the initial number of animals is 21:
#' model <- ecxsys(
#'     concentration = c(0, 0.03, 0.3, 3, 30),
#'     hormesis_concentration = 0.3,
#'     survival_tox_observed = c(17, 15.2, 18.8, 7, 0),
#'     survival_tox_env_observed = c(4.8, 4.6, 6.4, 0, 0),
#'     survival_max = 21
#' )
#'
#' @references \href{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).}
#'
#' @export
ecxsys <- function(concentration,
                   hormesis_concentration,
                   survival_tox_observed,
                   survival_tox_env_observed = NULL,
                   survival_max = 100,
                   curves_concentration_max = NULL,
                   p = 3.2,
                   q = 3.2) {
    output <- list(args = as.list(environment()))
    class(output) <- c("ecxsys", class(output))

    original_options <- options()
    on.exit(reset_options(original_options))


    # input validation ----------------------------------------------------
    if (survival_max <= 0) {
        stop("survival_max must be >= 0")
    }
    if (length(concentration) != length(survival_tox_observed)) {
        stop("concentration and survival_tox_observed must have the same length.")
    }
    if (length(concentration) > length(unique(concentration))) {
        stop("Concentrations must be unique.")
    }
    if (length(hormesis_concentration) != 1) {
        stop("hormesis_concentration must be a single number.")
    }
    if (!hormesis_concentration %in% concentration) {
        stop("hormesis_concentration must equal one of the concentration values.")
    }
    hormesis_index = which(hormesis_concentration == concentration)
    if (hormesis_index <= 2 || hormesis_index >= length(concentration)) {
        stop("hormesis_concentration must be greater than the lowest ",
             "non-control concentration and less than the highest concentration")
    }
    if (is.null(survival_tox_env_observed)) {
        with_env <- FALSE
        survival_tox_env_observed <- rep(NA, length(concentration))
    } else {
        with_env <- TRUE
    }
    output$with_env <- with_env
    # Creating all_observations makes it easier to test some assumptions.
    all_observations <- survival_tox_observed
    if (with_env) {
        if (length(survival_tox_observed) != length(survival_tox_env_observed)) {
            stop("survival_tox_observed and survival_tox_env_observed ",
                 "must have the same length.")
        }
        all_observations <- c(all_observations, survival_tox_env_observed)
    }
    if (any(is.na(c(all_observations, concentration)))) {
        stop("Values containing NA are not supported.")
    }
    if (any(all_observations > survival_max | all_observations < 0)) {
        stop("Observed survival must be between 0 and survival_max.")
    }
    conc_shift <- 2  # Powers of ten to shift the control downwards from the
    # second lowest concentration. This is required to approximate 0 because
    # of the logarithmic axis.
    if (is.unsorted(concentration)) {
        stop("The values must be sorted by increasing concentration.")
    }
    if (any(concentration < 0)) {
        stop("Concentrations must be >= 0")
    }
    if (min(concentration) > 0) {
        stop("No control is given. The first concentration must be 0.")
    }
    min_conc <- 10 ^ floor(log10(concentration[2]) - conc_shift)


    # scale observed survival ---------------------------------------------
    # scale the observed survival to [0,1] to make calculations independent of
    # value of the theoretical maximum survival
    survival_tox_observed <- survival_tox_observed / survival_max
    survival_tox_env_observed <- survival_tox_env_observed / survival_max


    # traditional log-logistic model (LL.5) -------------------------------
    LL5_tox <- fit_LL5_model(min_conc, concentration,
                             survival_tox_observed,
                             original_options)
    output$survival_tox_LL5_mod <- LL5_tox$survival_LL5_mod
    output$survival_tox_LL5 <- LL5_tox$survival_LL5 * survival_max
    if (with_env) {
        LL5_tox_env <- fit_LL5_model(min_conc, concentration,
                                     survival_tox_env_observed,
                                     original_options)
        output$survival_tox_env_LL5_mod <- LL5_tox_env$survival_LL5_mod
        output$survival_tox_env_LL5 <- LL5_tox_env$survival_LL5 * survival_max
    }


    # interpolation between subhormesis and hormesis ----------------------
    n_new <- 3  # number of new points
    concentration <- interpolate_sub_horm(
        concentration,
        hormesis_index - 1,
        hormesis_index,
        n_new,
        TRUE
    )
    survival_tox_observed <- interpolate_sub_horm(
        survival_tox_observed,
        hormesis_index - 1,
        hormesis_index,
        n_new
    )
    if (with_env) {
        survival_tox_env_observed <- interpolate_sub_horm(
            survival_tox_env_observed,
            hormesis_index - 1,
            hormesis_index,
            n_new
        )
    }
    hormesis_index <- hormesis_index + n_new
    # In the output return only the values at the original concentrations
    # and exclude those at the interpolated concentrations:
    keep <- !seq_along(concentration) %in% seq(hormesis_index - n_new, hormesis_index - 1)


    # survival_tox --------------------------------------------------------
    survival_tox <- c(1, survival_tox_observed[-1])
    survival_tox[2:(hormesis_index - 1)] <- NA
    survival_tox_mod <- drc::drm(survival_tox ~ concentration, fct = drc::W1.2())
    options(original_options["warn"])
    survival_tox <- predict(survival_tox_mod, data.frame(concentration = concentration))
    output$survival_tox_mod <- survival_tox_mod
    output$survival_tox <- survival_tox[keep] * survival_max

    stress_tox <- survival_to_stress(survival_tox, p, q)
    output$stress_tox <- stress_tox[keep]


    # system stress without environmental stress --------------------------
    fit_sys_output <- fit_sys(
        survival_to_stress(survival_tox_observed, p, q),
        stress_tox,
        stress_tox,
        hormesis_index,
        original_options
    )
    output$sys_tox_observed <- fit_sys_output$sys[keep]
    output$sys_tox_mod <- fit_sys_output$sys_mod
    if (inherits(fit_sys_output$sys_mod, "lm")) {
        warning(
            "Using a horizontal linear model for sys_tox_mod ",
            "because the Weibull model did not converge."
        )
    }
    sys_tox <- fit_sys_output$sys_modeled
    output$sys_tox <- sys_tox


    # modeled survival without environmental stress -----------------------
    stress_tox_sys <- stress_tox + sys_tox
    survival_tox_sys <- stress_to_survival(stress_tox_sys, p, q)
    output$survival_tox_sys <- survival_tox_sys[keep] * survival_max


    if (with_env) {
        # env stress ------------------------------------------------------
        stress_tox_env_observed <- survival_to_stress(survival_tox_env_observed, p, q)
        stress_env <- (stress_tox_env_observed - stress_tox)[hormesis_index]
        stress_env <- clamp(stress_env)
        output$stress_env <- stress_env

        stress_tox_env <- stress_tox + stress_env
        output$stress_tox_env <- stress_tox_env[keep]
        survival_tox_env <- stress_to_survival(stress_tox_env, p, q)
        output$survival_tox_env <- survival_tox_env[keep] * survival_max


        # system stress with environmental stress -------------------------
        fit_sys_output <- fit_sys(
            stress_tox_env_observed,
            stress_tox_env,
            stress_tox,
            hormesis_index,
            original_options
        )
        output$sys_tox_env_observed <- fit_sys_output$sys[keep]
        output$sys_tox_env_mod <- fit_sys_output$sys_mod
        if (inherits(fit_sys_output$sys_mod, "lm")) {
            warning(
                "Using a horizontal linear model for sys_tox_env_mod ",
                "because the Weibull model did not converge."
            )
        }
        sys_tox_env <- fit_sys_output$sys_modeled
        output$sys_tox_env <- sys_tox_env


        # modeled survival with environmental stress ----------------------
        stress_tox_env_sys <- stress_tox_env + sys_tox_env
        survival_tox_env_sys <- stress_to_survival(stress_tox_env_sys, p, q)
        output$survival_tox_env_sys <- survival_tox_env_sys[keep] * survival_max
    }


    # curves -------------------------------------------------------
    # In order to generate a break in the x-axis the concentration vector must
    # also be broken in two. The left part of the axis is supposed to be at
    # 0 but because it's a log axis I have to make the values just really
    # small. The concentrations in the gap won't be used for plotting later.
    len_curves <- 1000  #  number of points to approximate the curves
    conc_adjust_factor <- 10^-5
    output$conc_adjust_factor <- conc_adjust_factor
    if (is.null(curves_concentration_max)) {
        curves_concentration_max <- max(concentration)
    } else if (curves_concentration_max < min(concentration[concentration > 0])) {
        stop("'curves_concentration_max' is too low.")
    }
    curves_concentration <- 10 ^ seq(
        log10(min_conc * conc_adjust_factor),
        log10(curves_concentration_max),
        length.out = len_curves
    )
    output$curves <- predict_ecxsys(output, curves_concentration)

    # Add a column of concentrations which helps in plotting. It does not make
    # sense to show concentrations many orders of magnitude lower than the
    # lowest measured concentration. So this cuts out a large chunk and raises
    # the concentrations left of the cut so they make a nicer axis.
    conc_axis <- make_axis_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
}


reset_options <- function(original_options) {
    # Reset all the options which have changed.

    # You may ask why I don't just reset the options using
    # options(original_options). The reason is that when I do this and ecxsys
    # generates warnings then those warnings don't show up in the console. I
    # don't know why, but resetting only the options which have changed
    # alleviates that problem.

    changed <- list()
    for (n in names(original_options)) {
        orig_opt <- original_options[[n]]
        if (!identical(orig_opt, getOption(n))) {
            changed[n] <- orig_opt
        }
    }
    options(changed)
}


fit_sys <- function(stress_external_observed,
                    stress_external_modeled,
                    stress_tox,
                    hormesis_index,
                    original_options) {
    sys <- stress_external_observed - stress_external_modeled
    sys <- clamp(sys)
    sys[hormesis_index:length(sys)] <- 0
    # Add sys to the output before it is fitted:
    out <- list(sys = sys)
    sys_mod <- tryCatch(
        {
            # There is no other way to suppress that one error message
            # from optim except by changing the options temporarily.
            warn_error_original <- original_options[c("warn", "show.error.messages")]
            options(show.error.messages = FALSE)
            drc::drm(sys ~ stress_tox, fct = drc::W1.3())
        },
        error = function(e) {
            # Failure to converge often happens when all or almost all sys
            # stress values are zero. Fitting a linear model in this case seems
            # to be the most appropriate remedy.
            stress_tox <- range(stress_tox)
            sys <- c(0, 0)
            lm(sys ~ stress_tox)
        },
        finally = options(warn_error_original)
    )
    out$sys_mod <- sys_mod
    out$sys_modeled <- unname(
        predict(sys_mod, data.frame(stress_tox = stress_tox))
    )
    out
}


moving_weighted_mean <- function(x) {
    # This is used to smooth out points which are jumping up and down.
    count <- rep(1, length(x))
    x_diff <- diff(x) * -1
    while (any(x_diff < 0, na.rm = TRUE)) {
        i <- which(x_diff < 0)[1]
        j <- i + 1
        x[i] <- weighted.mean(x[i:j], count[i:j])
        x <- x[-j]
        count[i] <- count[i] + count[j]
        count <- count[-j]
        x_diff <- diff(x) * -1
    }
    rep(x, count)
}


fit_LL5_model <- function(min_conc,
                          concentration,
                          survival_observed,
                          original_options) {
    # The traditional log-logistic model, here using the five-parameter
    # log-logistic function drc::LL.5.
    conc_with_control_shifted <- c(min_conc, concentration[-1])
    survival_observed_averaged <- moving_weighted_mean(survival_observed)
    interpolated <- approx(
        log10(conc_with_control_shifted),
        survival_observed_averaged,
        n = 10
    )
    conc_interpolated <- 10^interpolated$x
    survival_interpolated <- interpolated$y
    survival_LL5_mod <- drc::drm(
        survival_interpolated ~ conc_interpolated,
        fct = drc::LL.5(fixed = c(NA, 0, survival_observed_averaged[1], NA, NA))
    )
    options(original_options["warn"])
    survival_LL5 <- predict(
        survival_LL5_mod,
        data.frame(conc_interpolated = concentration)
    )
    list(
        survival_LL5_mod = survival_LL5_mod,
        survival_LL5 = survival_LL5
    )
}


interpolate_sub_horm <- function(x,
                                 from_index,
                                 to_index,
                                 n_new,
                                 logarithmic = FALSE) {
    # Interpolate between subhormesis and hormesis. This makes the curves
    # smoother, less extreme. Without it the slope between subhormesis and
    # hormesis would be much steeper in many cases.
    len <- n_new + 2  # Add 2 because seq() includes the left and right end.
    if (logarithmic) {
        x_new <- 10^seq(log10(x[from_index]), log10(x[to_index]), length.out = len)
    } else {
        x_new <- seq(x[from_index], x[to_index], length.out = len)
    }
    append(x, x_new[-c(1, len)], from_index)
}


make_axis_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
    )
}