path: root/R/twa.R

# Copyright (C) 2016,2017  Johannes Ranke

# Contact: jranke@uni-bremen.de
# This file is part of the R package pfm

# This program 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 <http://www.gnu.org/licenses/>

#' Create a time series of decline data
#'
#' @param x When numeric, this is the half-life to be used for an exponential
#'   decline. When a character string specifying a parent decline model is given
#'   e.g. \code{FOMC}, \code{parms} must contain the corresponding parameters.
#'   If x is an \code{\link{mkinfit}} object, the decline is calculated from this
#'   object.  
#' @param ini The initial amount. If x is an \code{\link{mkinfit}} object, and
#'   ini is 'model', the fitted initial concentrations are used. Otherwise, ini
#'   must be numeric. If it has length one, it is used for the parent and
#'   initial values of metabolites are zero, otherwise, it must give values for
#'   all observed variables.
#' @param t_end End of the time series
#' @param res Resolution of the time series
#' @param ... Further arguments passed to methods
#' @return An object of class \code{one_box}, inheriting from \code{\link{ts}}.
#' @importFrom stats filter frequency time ts
#' @export
#' @examples
#' # Only use a half-life
#' pred_0 <- one_box(10)
#' plot(pred_0)
#'
#' # Use a fitted mkinfit model
#' require(mkin)
#' fit <- mkinfit("FOMC", FOCUS_2006_C, quiet = TRUE)
#' pred_1 <- one_box(fit)
#' plot(pred_1)
#'
#' # Use a model with more than one observed variable
#' m_2 <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO"))
#' fit_2 <- mkinfit(m_2, FOCUS_2006_D, quiet = TRUE)
#' pred_2 <- one_box(fit_2, ini = "model")
#' plot(pred_2)
one_box <- function(x, ini, ...,
  t_end = 100, res = 0.01)
{
  UseMethod("one_box")
}

#' @rdname one_box
#' @export
one_box.numeric <- function(x, ini = 1, ...,
  t_end = 100, res = 0.01)
{
  half_life = x
  k = log(2)/half_life
  t_out <- seq(0, t_end, by = res)
  raw <- matrix(ini * exp( - k * t_out), ncol = 1)
  dimnames(raw) <- list(NULL, "parent")
  result <- ts(raw, 0, t_end, frequency = 1/res)
  class(result) <- c("one_box", "ts")
  return(result)
}

#' @rdname one_box
#' @param parms A named numeric vector containing the model parameters
#' @export
one_box.character <- function(x, ini = 1, parms, ...,
  t_end = 100, res = 0.01)
{
  parent_models_available = c("SFO", "FOMC", "DFOP", "HS", "SFORB", "IORE")
  if (length(x) == 1 & x %in% parent_models_available) {
    m <- mkinmod(parent = mkinsub(x))
  } else {
    stop("If you specify the decline model using a character string, ",
         "x has to be one of\n ",
          paste(parent_models_available, collapse = ", "))
  }

  if (!setequal(names(parms), m$par)) {
    stop("Please supply the parameters\n",
         paste(m$par, collapse = ", "))
  }

  t_out <- seq(0, t_end, by = res)
  pred <- mkinpredict(m, odeparms = parms, odeini = c(parent = ini),
                      outtimes = t_out, solution_type = "analytical")[-1]
  result <- ts(pred, 0, t_end, frequency = 1/res)
  class(result) <- c("one_box", "ts")
  return(result)
}

#' @rdname one_box
#' @importFrom mkin mkinpredict
#' @export
one_box.mkinfit <- function(x, ini = "model", ..., t_end = 100, res = 0.01) {
  fit <- x
  if (ini[1] == "model") {
    odeini = x$bparms.state
  } else {
    if (!is.numeric(ini[1])) stop ("Argument ini can only be 'model' or numeric")
    if (length(ini) == 1) odeini <- c(ini[1], rep(0, length(fit$mkinmod$spec) - 1))
    else odeini = ini
    names(odeini) <- names(fit$mkinmod$spec)
  }

  t_out = seq(0, t_end, by = res)
  if (length(fit$mkinmod$spec) == 1) solution_type = "analytical"
  else solution_type = "deSolve"

  tmp <- mkinpredict(fit$mkinmod, odeparms = fit$bparms.ode, odeini = odeini,
    outtimes = t_out, solution_type = solution_type)[-1]
  result <- ts(tmp, 0, t_end, frequency = 1/res)
  class(result) <- c("one_box", "ts")
  return(result)
}

#' Plot time series of decline data
#'
#' @param x The object of type \code{\link{one_box}} to be plotted
#' @param xlim Limits for the x axis
#' @param ylim Limits for the y axis
#' @param xlab Label for the x axis
#' @param ylab Label for the y axis
#' @param max_twa If a numeric value is given, the maximum time weighted
#'   average concentration(s) is/are shown in the graph.
#' @param max_twa_var Variable for which the maximum time weighted average should
#'   be shown if max_twa is not NULL.
#' @param ... Further arguments passed to methods
#' @importFrom stats plot.ts
#' @seealso \code{\link{sawtooth}}
#' @export
#' @examples
#' dfop_pred <- one_box("DFOP", parms = c(k1 = 0.2, k2 = 0.02, g = 0.7))
#' plot(dfop_pred)
#' plot(sawtooth(dfop_pred, 3, 7), max_twa = 21)
#'
#' # Use a fitted mkinfit model
#' m_2 <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO"))
#' fit_2 <- mkinfit(m_2, FOCUS_2006_D, quiet = TRUE)
#' pred_2 <- one_box(fit_2, ini = 1)
#' pred_2_saw <- sawtooth(pred_2, 2, 7)
#' plot(pred_2_saw, max_twa = 21, max_twa_var = "m1")
plot.one_box <- function(x,
                         xlim = range(time(x)), ylim = c(0, max(x)),
                         xlab = "Time", ylab = "Residue",
                         max_twa = NULL, max_twa_var = dimnames(x)[[2]][1], ...)
{
  obs_vars <- dimnames(x)[[2]]
  plot.ts(x, plot.type = "single", xlab = xlab, ylab = ylab,
          lty = 1:length(obs_vars), col = 1:length(obs_vars),
          las = 1, xlim = xlim, ylim = ylim)
  if (!is.null(max_twa)) {
    x_twa <- max_twa(x, window = max_twa)
    value <- x_twa$max[max_twa_var]
    rect(x_twa$window_start[max_twa_var], 0,
         x_twa$window_end[max_twa_var], value, col = "grey")
    text(x_twa$window_end[max_twa_var], value, paste("Maximum:", signif(value, 3)), pos = 4)
    # Plot a second time to cover the grey rectangle
    matlines(time(x), as.matrix(x), lty = 1:length(obs_vars), col = 1:length(obs_vars))
  }
}

#' Create decline time series for multiple applications
#'
#' If the application pattern is specified in \code{applications}, 
#' \code{n} and \code{i} are disregarded.
#' @param x A \code{\link{one_box}} object
#' @param n The number of applications. If \code{applications} is specified, \code{n} is ignored
#' @param i The interval between applications. If \code{applications} is specified, \code{i}
#'   is ignored
#' @param applications A data frame holding the application times in the first column and
#'   the corresponding amounts applied in the second column.
#' @export
#' @examples
#' applications = data.frame(time = seq(0, 14, by = 7), amount = c(1, 2, 3))
#' pred <- one_box(10)
#' plot(sawtooth(pred, applications = applications))
#'
#' m_2 <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO"))
#' fit_2 <- mkinfit(m_2, FOCUS_2006_D, quiet = TRUE)
#' pred_2 <- one_box(fit_2, ini = 1)
#' pred_2_saw <- sawtooth(pred_2, 2, 7)
#' plot(pred_2_saw, max_twa = 21, max_twa_var = "m1")
#'
#' max_twa(pred_2_saw)
sawtooth <- function(x, n = 1, i = 365,
                     applications = data.frame(time = seq(0, (n - 1) * i, length.out = n),
                                               amount = 1))
{
  n_obs = ncol(as.matrix(x))
  t_end = max(time(x))
  freq = frequency(x)
  empty <- ts(matrix(0, nrow = t_end * freq, ncol = n_obs), 0, t_end, freq)
  result <- empty
  for (i_app in 1:nrow(applications)) {
    t_app <- applications[i_app, "time"]
    amount_app <- applications[i_app, "amount"]
    if (t_app == 0) {
      result <- result + x * amount_app
    } else {
      lag_phase <- as.matrix(empty)[1:(t_app * freq), , drop = FALSE]
      app_phase <- amount_app * as.matrix(x)[1:((t_end - t_app) * freq + 1), , drop = FALSE]
      app_ts <- ts(rbind(lag_phase, app_phase), 0, t_end, frequency = freq)
      result <- result + app_ts
    }
  }
  class(result) = c("one_box", "ts")
  dimnames(result) <- dimnames(x)
  return(result)
}

#' Calculate a time weighted average concentration
#'
#' The moving average is built only using the values in the past, so
#' the earliest possible time for the maximum in the time series returned
#' is after one window has passed.
#'
#' @param x An object of type \code{\link{one_box}}
#' @param window The size of the moving window
#' @seealso \code{\link{max_twa}}
#' @importFrom stats start end
#' @export
#' @examples
#' pred <- sawtooth(one_box(10),
#'   applications = data.frame(time = c(0, 7), amount = c(1, 1)))
#' max_twa(pred)
twa <- function(x, window = 21) UseMethod("twa")

#' @rdname twa
#' @export
twa.one_box <- function(x, window = 21)
{
  length_ts <- end(x) - start(x)
  if (window >= length_ts[1]) {
    stop("The window must be smaller than the length of the time series")
  }

  resolution = 1/frequency(x)
  n_filter = window/resolution
  result = filter(x, rep(1/n_filter, n_filter), method = "convolution", sides = 1)
  class(result) = c("one_box", "ts")
  dimnames(result) <- dimnames(x)
  return(result)
}

#' The maximum time weighted average concentration for a moving window
#'
#' If you generate your time series using \code{\link{sawtooth}},
#' you need to make sure that the length of the time series allows
#' for finding the maximum. It is therefore recommended to check this using
#' \code{\link{plot.one_box}} using the window size for the argument
#' \code{max_twa}.
#' 
#' The method working directly on fitted \code{\link{mkinfit}} objects uses the
#' equations given in the PEC soil section of the FOCUS guidance and is restricted
#' SFO, FOMC and DFOP models and to the parent compound
#' @references FOCUS (2006) \dQuote{Guidance Document on Estimating Persistence and
#'   Degradation Kinetics from Environmental Fate Studies on Pesticides in EU
#'   Registration} Report of the FOCUS Work Group on Degradation Kinetics,
#'   EC Document Reference Sanco/10058/2005 version 2.0, 434 pp,
#'   \url{http://esdac.jrc.ec.europa.eu/projects/degradation-kinetics}
#' @seealso \code{\link{twa}}
#' @inheritParams twa
#' @export
#' @examples
#' pred <- sawtooth(one_box(10),
#'   applications = data.frame(time = c(0, 7), amount = c(1, 1)))
#' max_twa(pred)
#' pred_FOMC <- mkinfit("FOMC", FOCUS_2006_C, quiet = TRUE)
#' max_twa(pred_FOMC)
max_twa <- function(x, window = 21) UseMethod("max_twa")

#' @export
max_twa.mkinfit <- function(x, window = 21) {
  fit <- x
  parms.all <- c(fit$bparms.optim, fit$bparms.fixed)
  obs_vars <- fit$obs_vars
  if (length(obs_vars) > 1) {
    warning("Calculation of maximum time weighted average concentrations is",
            "currently only implemented for the parent compound using",
            "analytical solutions")
  }
  obs_var <- obs_vars[1]
  spec = fit$mkinmod$spec
  type = spec[[1]]$type

  M0 <- parms.all[paste0(obs_var, "_0")]

  if (type == "SFO") {
    k_name <- paste0("k_", obs_var)
    if (fit$mkinmod$use_of_ff == "min") {
      k_name <- paste0(k_name, "_sink")
    }
    k <- parms.all[k_name]
    twafunc <- function(t) {
      M0 * (1 - exp(- k * t)) / (k * t)
    }
  }
  if (type == "FOMC") {
    alpha <- parms.all["alpha"]
    beta <- parms.all["beta"]
    twafunc <- function(t) {
      M0 * (beta)/(t * (1 - alpha)) * ((t/beta + 1)^(1 - alpha) - 1)
    }
  }
  if (type == "DFOP") {
    k1 <- parms.all["k1"]
    k2 <- parms.all["k2"]
    g <- parms.all["g"]
    twafunc <- function(t) {
      M0/t * ((g/k1) * (1 - exp(- k1 * t)) + ((1 - g)/k2) * (1 - exp(- k2 * t)))
    }
  }
  if (type %in% c("HS", "IORE", "SFORB")) {
    stop("Calculation of maximum time weighted average concentrations is currently ",
         "not implemented for the ", type, " model.")
  }
  res <- twafunc(t = window)
  names(res) <- window
  return(res)
}

#' @export
max_twa.one_box <- function(x, window = 21)
{
  freq = frequency(x)

  twa_ts <- twa(x, window = window)
  window_end <- apply(twa_ts, 2, which.max) / freq
  result <- list()
  result$max <- apply(twa_ts, 2, max, na.rm = TRUE)
  result$window_start <- window_end - window
  result$window_end <- window_end
  return(result)
}