R/extcoef_logistic.R

In growthrates: Estimate Growth Rates from Experimental Data

#' Extended Set of Coefficients of a Logistic Growth Model
#'
#' Estimate model-specific derived parameters of the logistic growth
#' model
#'
#' @param object model object fited by \code{fit_growthmodel}
#' @param quantile fraction of the capacity parameter (\code{K}) for the quantile method
#' @param time 2-valued vector of the search interval for the independent
#'   variable (\code{time}).
#'   Note: this needs to be set this manually if saturation is not
#'   reached within the observation time period taken from the data.
#' @param ... reserved for future extensions
#'
#' @return vector that contains the fitted parameters and some
#'   derived characteristics (extended parameters) of the logistic
#'   function.
#'
#' @details This function returns the estimated parameters of a logistic growth model
#'  (\code{y0}, \code{mumax}, \code{K}) and a series of estimates for the time
#'  of approximate saturation.
#'  The estimates are defined as follows:
#'  \itemize{
#'    \item \code{turnpoint}: time of turnpoint (50\% saturation)
#'    \item \code{sat1}: time of the minimum of the 2nd derivative
#'    \item \code{sat2}: time of the intercept between the steepest increase
#'      (the tangent at \code{mumax}) and the carrying capacity \code{K}
#'    \item \code{sat3}: time when a quantile of \code{K} (default 0.95)
#'      is reached
#'  }
#'
#'  This function is normally not directly called by the user.
#'  It is usually called indirectly from \code{coef} or \code{results} if
#'    \code{extended=TRUE}.
#'
#' @note
#'  The estimates for the turnpoint and the time of approximate saturation
#'    (\code{sat1}, \code{sat2}, \code{sat3}) may be unreliable, if saturation
#'    is not reached within the observation time period. See example below.
#'  A set of extended parameters exists currently only for the standard logistic
#'    growth model (\code{grow_logistic}).
#'  The code and naming of the parameters is preliminary and may change in
#'    future versions.
#'
#'
#' @examples
#'
#' ## =========================================================================
#' ## The 'extended parameters' are usually derived
#' ## =========================================================================
#'
#' data(antibiotic)
#'
#' ## fit a logistic model to a single data set
#' dat <- subset(antibiotic, conc==0.078 & repl=="R4")
#'
#' parms <- c(y0=0.01, mumax=0.2, K=0.5)
#' fit <- fit_growthmodel(grow_logistic, parms, dat$time, dat$value)
#' coef(fit, extended=TRUE)
#'
#' ## fit the logistic to all data sets
#' myData <- subset(antibiotic, repl=="R3")
#' parms <- c(y0=0.01, mumax=0.2, K=0.5)
#' all <- all_growthmodels(value ~ time | conc,
#'                          data = myData, FUN=grow_logistic,
#'                          p = parms, ncores = 2)
#'
#'
#' par(mfrow=c(3,4))
#' plot(all)
#' results(all, extended=TRUE)
#' ## we see that the the last 3 series (10...12) do not go into saturation
#' ## within the observation time period.
#'
#' ## We can try to extend the search range:
#' results(all[10:12], extended=TRUE, time=c(0, 5000))
#'
#'
#' ## =========================================================================
#' ## visualisation how the 'extended parameters' are derived
#' ## =========================================================================
#'
#' # Derivatives of the logistic:
#' #   The 1st and 2nd derivatives are internal functions of the package.
#' #   They are used here for the visualisation of the algorithm.
#'
#' deriv1 <- function(time, y0, mumax, K) {
#'   ret <- (K*mumax*y0*(K - y0)*exp(mumax * time))/
#'     ((K + y0 * (exp(mumax * time) - 1))^2)
#'   unname(ret)
#' }
#'
#' deriv2 <- function(time, y0, mumax, K) {
#'   ret <- -(K * mumax^2 * y0 * (K - y0) * exp(mumax * time) *
#'              (-K + y0 * exp(mumax * time) + y0))/
#'     (K + y0 * (exp(mumax * time) - 1))^3
#'   unname(ret)
#' }
#' ## =========================================================================
#'
#' data(bactgrowth)
#' ## extract one growth experiment by name
#' dat <- multisplit(bactgrowth, c("strain", "conc", "replicate"))[["D:0:1"]]
#'
#'
#' ## unconstraied fitting
#' p <- c(y0 = 0.01, mumax = 0.2, K = 0.1) # start parameters
#' fit1 <- fit_growthmodel(FUN = grow_logistic, p = p, dat$time, dat$value)
#' summary(fit1)
#' p <- coef(fit1, extended=TRUE)
#'
#' ## copy parameters to separate variables to improve readability ------------
#' y0 <-    p["y0"]
#' mumax <- p["mumax"]
#' K  <-    p["K"]
#' turnpoint <- p["turnpoint"]
#' sat1 <-  p["sat1"]  # 2nd derivative
#' sat2 <-  p["sat2"]  # intercept between steepest increase and K
#' sat3 <-  p["sat3"]  # a given quantile of K, default 95\%
#'
#' ## show saturation values in growth curve and 1st and 2nd derivatives ------
#' opar <- par(no.readonly=TRUE)
#' par(mfrow=c(3, 1), mar=c(4,4,0.2,0))
#' plot(fit1)
#'
#' ## 95% saturation
#' abline(h=0.95*K, col="magenta", lty="dashed")
#'
#' ## Intercept between steepest increase and 100% saturation
#' b <- deriv1(turnpoint, y0, mumax, K)
#' a <- K/2 - b*turnpoint
#' abline(a=a, b=b, col="orange", lty="dashed")
#' abline(h=K, col="orange", lty="dashed")
#' points(sat2, K, pch=16, col="orange")
#' points(turnpoint, K/2, pch=16, col="blue")
#'
#' ## sat2 is the minimum of the 2nd derivative
#' abline(v=c(turnpoint, sat1, sat2, sat3),
#'        col=c("blue", "grey", "orange", "magenta"), lty="dashed")
#'
#' ## plot the derivatives
#' with(dat, plot(time, deriv1(time, y0, mumax, K), type="l", ylab="y'"))
#' abline(v=c(turnpoint, sat1), col=c("blue", "grey"), lty="dashed")
#'
#' with(dat, plot(time, deriv2(time, y0, mumax, K), type="l",  ylab="y''"))
#' abline(v=sat1, col="grey", lty="dashed")
#' par(opar)
#'
#' @export
#'
extcoef_logistic <- function(object, quantile=0.95, time=NULL, ...) {

  ## analytical derivatives of the logistic as local functions
  deriv1 <- function(time, y0, mumax, K) {
    (K*mumax*y0*(K - y0)*exp(mumax * time))/((K + y0 * (exp(mumax * time) - 1))^2)
  }

  deriv2 <- function(time, y0, mumax, K) {
    -(K * mumax^2 * y0 * (K - y0) * exp(mumax * time) *
        (-K + y0 * exp(mumax * time) + y0))/(K + y0 * (exp(mumax * time) - 1))^3
  }


  ## argument checks
  if (is.null(time)) time <- object@obs$time
  if (is.null(quantile)) quantile <- 0.95
  ## todo: check if fitted function was grow_logistic

  p  <- coef(object)
  r2 <- rsquared(object)

  y0     <- p["y0"]
  K      <- p["K"]
  mumax  <- p["mumax"]
  trange <- range(time)

  ## indentify max/min by numerical search
  time_turn1 <- optimize(deriv1, trange, y0=y0, mumax=mumax, K=K, maximum=TRUE)$maximum
  time_turn2 <- optimize(deriv2, trange, y0=y0, mumax=mumax, K=K)$minimum

  time_quantile <- (log((quantile * (y0 - K))/((quantile - 1) * y0)))/mumax

  ## intercept between steepest increase and saturation
  y_turn1 <- deriv1(time_turn1, y0, mumax, K)
  b <- y_turn1
  a <- K/2 - b * time_turn1

  time_sat <- (K-a)/b

  c(y0 = unname(y0),
    mumax = unname(mumax),
    K = unname(K),
    turnpoint = time_turn1,
    sat1 = unname(time_turn2),
    sat2 = unname(time_sat),
    sat3 = unname(time_quantile))
}