R/baranyi.R

In MicrobialGrowth: Estimates Growth Parameters from Models and Plots the Curve

# Copyright 2023, Florentin L'Homme, Clarisse Breard
#
# 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 <https://www.gnu.org/licenses/>


#' Baranyi object
#'
#' A MicrobialGrowth object specialized for the baranyi model.
#' Most of the methods are pre-implemented (some of these can be overwritten for a specific regression/create function).
#' Must be completed for `data`, `isValid` (regression successful), etc.
#'
#' @param ... further arguments passed to or from other methods.
#'
#' @return a Baranyi object skeleton.
#' @importFrom stats integrate
#' @export
#'
#' @details the three dots `...` are passed to the \link{.new.MicrobialGrowth.core} function.
#'
#' @examples
#' # First, create the skeleton.
#' model.object = .new.baranyi.core()
#'
#' # Then complete with data, functions, etc.
#' model.object$data$x = c(1,2,3)
#' model.object$data$y = c(1,2,3)
#' model.object$coefficients = list(N0 = 0, Nmax=0, mu=0, lambda=0)
#'
#' # You can print, plot, etc., with the generic functions of MicrobialGrowth.
#' print(model.object)
#' ##MicrobialGrowth, model specialized.model:
#' ##    N0   Nmax     mu lambda
#' ##     0      0      0      0
#' plot(model.object)
#'
#' # Don't forget to change `isValid` to TRUE to confirm the success of the regression.
#' model.object$isValid = TRUE # Not a good idea here, since we have no `reg` value.
.new.baranyi.core = function(...) {
  env = .new.MicrobialGrowth.core(...)
  class(env) = c("baranyi", class(env))

  env$f$confint = function(level = 0.95) {
    if (!is.null(env$reg)) {
      return( suppressMessages(nlstools::confint2(env$reg, level = level)) )
    }
    else {
      stop("Since the regression did not pass, no confidence interval can be inferred.")
    }
  }

  env$f$confint.lower = function(level = 0.95, base = exp(1)) {
    if (is.null(base)) { stop("confint.lower is not available for base=NULL") }

    conf = env$f$confint(level)

    if (conf["mu", 1] <= 0) {
      conf["mu",1]  = 1/.Machine$double.xmax
      warning("Lower bound of the confidence interval for mu was below or equal to 0. Its value was replaced by a value close to 0 (=1/.Machine$double.xmax).")
    }

    cond1 = do.call(getFormula("baranyi"), list(N0=conf["N0",2], Nmax=conf["Nmax",1], mu=conf["mu",1], lambda=conf["lambda",2], base = base))
    cond2 = function(x) {
      # x is a vector of values
      y = rep(NA, length(x))

      alpha = log(conf["lambda",2]/(conf["lambda",2]-x))/x

      y[conf["mu",2] < alpha] = do.call(getFormula("baranyi"), list(N0=conf["N0",2], Nmax=conf["Nmax",1], mu=conf["mu",1], lambda=conf["lambda",2], base = base))(x[conf["mu",2] < alpha])
      y[conf["mu",1] > alpha] = do.call(getFormula("baranyi"), list(N0=conf["N0",2], Nmax=conf["Nmax",1], mu=conf["mu",2], lambda=conf["lambda",2], base = base))(x[conf["mu",1] > alpha])
      y[!(conf["mu",2] < alpha) & !(conf["mu",1] > alpha) & (abs(exp(x*alpha)-exp(x*conf["mu",1]))>abs(exp(x*alpha)-exp(x*conf["mu",2])))] = do.call(getFormula("baranyi"), list(N0=conf["N0",2], Nmax=conf["Nmax",1], mu=conf["mu",1], lambda=conf["lambda",2], base = base))(x[!(conf["mu",2] < alpha) & !(conf["mu",1] > alpha) & (abs(exp(x*alpha)-exp(x*conf["mu",1]))>abs(exp(x*alpha)-exp(x*conf["mu",2])))])
      y[!(conf["mu",2] < alpha) & !(conf["mu",1] > alpha) & !(abs(exp(x*alpha)-exp(x*conf["mu",1]))>abs(exp(x*alpha)-exp(x*conf["mu",2])))] = do.call(getFormula("baranyi"), list(N0=conf["N0",2], Nmax=conf["Nmax",1], mu=conf["mu",2], lambda=conf["lambda",2], base = base))(x[!(conf["mu",2] < alpha) & !(conf["mu",1] > alpha) & !(abs(exp(x*alpha)-exp(x*conf["mu",1]))>abs(exp(x*alpha)-exp(x*conf["mu",2])))])

      return(y)
    }

    f = function(x) {
      y = rep(NA, length(x))
      y[x >= (conf["lambda",2])] = cond1(x[x >= (conf["lambda",2])])
      y[!(x >= (conf["lambda",2]))] = cond2(x[!(x >= (conf["lambda",2]))])

      y[x < 0] = NA

      return(y)
    }

    return(f)
  }

  env$f$confint.upper = function(level = 0.95, base = exp(1)) {
    if (is.null(base)) { stop("confint.upper is not available for base=NULL") }

    conf = env$f$confint(level)

    if (conf["mu", 1] <= 0) {
      conf["mu",1]  = 1/.Machine$double.xmax
      warning("Lower bound of the confidence interval for mu was below or equal to 0. Its value was replaced by a value close to 0 (=1/.Machine$double.xmax).")
    }

    cond1 = do.call(getFormula("baranyi"), list(N0=conf["N0",1], Nmax=conf["Nmax",2], mu=conf["mu",2], lambda=conf["lambda",1], base = base))
    cond2 = function(x) {
      # x is a vector of values
      y = rep(NA, length(x))

      alpha = log(conf["lambda",1]/(conf["lambda",1]-x))/x

      y[conf["mu",2] < alpha] = do.call(getFormula("baranyi"), list(N0=conf["N0",1], Nmax=conf["Nmax",2], mu=conf["mu",2], lambda=conf["lambda",1], base = base))(x[conf["mu",2] < alpha])
      y[conf["mu",1] > alpha] = do.call(getFormula("baranyi"), list(N0=conf["N0",1], Nmax=conf["Nmax",2], mu=conf["mu",1], lambda=conf["lambda",1], base = base))(x[conf["mu",1] > alpha])
      y[!(conf["mu",2] < alpha) & !(conf["mu",1] > alpha)] = do.call(getFormula("baranyi"), list(N0=conf["N0",1], Nmax=conf["Nmax",2], mu=alpha[!(conf["mu",2] < alpha) & !(conf["mu",1] > alpha)] , lambda=conf["lambda",1], base = base))(x[!(conf["mu",2] < alpha) & !(conf["mu",1] > alpha)])

      return(y)
    }

    f = function(x) {
      y = rep(NA, length(x))
      y[x >= (conf["lambda",1])] = cond1(x[x >= (conf["lambda",1])])
      y[!(x >= (conf["lambda",1]))] = cond2(x[!(x >= (conf["lambda",1]))])

      y[x < 0] = NA

      return(y)
    }

    return(f)
  }

  env$f$doublingTime = function(level = 0.95) {
    # TODO : To rework or validate for Baranyi -----------------------------------
    #
    # # TODO : To check N0 sur le confint ? => risque de N0 négatif et donc de faire échouer log(Nmax/N0)
    # if (!baranyi.out$isValid) stop("Coefficient mu not defined, impossible to calculate the doubling time")
    #
    # confint = baranyi.out$f$confint(level)
    # N0 = append(confint["N0",], c(coefficient=baranyi.out$coefficient$N0), 1)
    # Nmax = append(confint["Nmax",], c(coefficient=baranyi.out$coefficient$Nmax), 1)
    # mu = append(confint["mu",], c(coefficient=baranyi.out$coefficient$mu), 1)
    # lambda = append(confint["lambda",], c(coefficient=baranyi.out$coefficient$lambda), 1)
    #
    # doublingTime = log(2)/mu
    # A = log(Nmax/N0)
    # start = lambda - (A*log(-log(log(2)/A)) - A)/(mu * exp(1))
    #
    # return( as.matrix(t(data.frame(doublingTime = doublingTime, start = start))) )
  }

  env$f$intersection = function(x = NULL, y = NULL, base = exp(1)) {
    if (!is.null(x) & !is.null(y)) { stop("Only one of x or y must be specified") }
    if (!is.null(x)) {
      return( env$f$formula(base = base)(x) )
    }
    else { # !is.null(y)
      if (is.null(base)) {
        return( log(1+(env$coefficients$Nmax/y - env$coefficients$Nmax/env$coefficients$N0)/(exp(-env$coefficients$mu*env$coefficients$lambda)*(1-env$coefficients$Nmax/y)))/env$coefficients$mu )
      }
      else {
        return( log(1+((env$coefficients$Nmax/env$coefficients$N0)*(1-base^y))/(exp(-env$coefficients$mu*env$coefficients$lambda)*(base^y-(env$coefficients$Nmax/env$coefficients$N0))))/env$coefficients$mu )
      }
    }
  }

  env$f$auc = function(from = NULL, to = NULL, curve = c("main", "lower", "upper"), level = 0.95, base = exp(1)) {
    if (is.null(from)) from = min(env$data$x)
    if (is.null(to)) to = max(env$data$x)
    curve = match.arg(curve)

    expr = switch (curve,
                   "main" = env$f$formula(base=base),
                   "lower" = env$f$confint.lower(level=level, base=base),
                   "upper" = env$f$confint.upper(level=level, base=base)
    )

    return(integrate(expr, from, to))
  }

  env$f$formula = function(base = exp(1)) {
    if (!env$isValid) { stop("Regression failed, no formula available") }
    return( do.call(getFormula("baranyi"), append(env$coefficients, list(base = base))) )
  }

  env
}



#' Baranyi regression function
#'
#' Regression function for Baranyi's model
#'
#' @param x index series or time series.
#' @param y values or list of values to regress (should not be logged).
#' @param clip a pair of values indicating in which interval to clip the data `y`. When `clip` is missing, default values are used.
#' @param start a named list of starting estimates. When `start` is missing, default values are used.
#' @param lower a named list of lower bounds. When `lower` is missing, default values are used.
#' @param upper a named list of upper bounds. When `upper` is missing, default values are used.
#' @param nls.args additional parameters to use when calling \link[stats]{nls}.
#' @param callbackError function to call on error during regression.
#'
#' @returns a MicrobialGrowth-object composed of
#' \item{call}{the matched call with several components.}
#' \item{coefficients}{coefficients obtained by regression.}
#' \item{data}{data used for regression, once the y values are clipped}
#' \item{f}{a list of functions such as `formula` to retrieve the function of the model with the coefficients obtained by regression, `confint` to retrieve the confidence intervals, etc.}
#' \item{isValid}{a boolean indicating whether the regression was successful or not.}
#' \item{message}{contains the error message if the regression fails, `NULL` otherwise.}
#' \item{reg}{the `nls` object returned by the `nls` function.}
#'
#' @include utils.R
#' @importFrom nlstools confint2
#' @importFrom utils modifyList tail
#' @importFrom stats lm nls
#' @export
#'
#' @details
#' The default values for `clip`, `start`, `lower` and `upper` are calculated based on the given data. These default values can be known through the `call` member of the returned value.
#'
#' The `nls.args` argument is a list that can contain any \link[stats]{nls} function argument except `formula`, `algorithm`, `start`, `lower` and `upper` which are already fixed (via a homonymous or hard-coded argument).
#'
#' For the `callbackError` argument, prefer the `stop` function to block or `warning` to not be blocking.
#'
#' @seealso \link{MicrobialGrowth}, \link{.baranyi.formula}
.MicrobialGrowth.baranyi = function(x, y, clip = c(-Inf, Inf), start = list(), lower = list(), upper = list(), nls.args = list(), callbackError = NULL) {
  clip = .checkMicrobialGrowthArgs(x=x, y=y, clip=clip, start=start, lower=lower, upper=upper, nls.args=nls.args, callbackError=callbackError)

  { # Removal of NA values
    idx = !(is.na(x) | is.na(y))
    x = x[idx]
    y = y[idx]
  }

  # Clip y-values to be log function friendly
  y[y < clip[1]] = clip[1]
  y[y > clip[2]] = clip[2]

  # Default values
  default.values = .getDefaultNlsValues(x, y, start, lower, upper)
  start = default.values$start
  lower = default.values$lower
  upper = default.values$upper


  # Creation of the return value
  baranyi.out = .new.baranyi.core()
  baranyi.out$data = list(x=x, y=y)
  baranyi.out$call = match.call()
  baranyi.out$call$clip=clip
  baranyi.out$call$start=start
  baranyi.out$call$lower=lower
  baranyi.out$call$upper=upper
  baranyi.out$call$nls.args=nls.args
  baranyi.out$call$callbackError=callbackError

  tryCatch(
    # The `nls` function may return an error if the regression fails.
    # In order to take a non-blocking approach (especially when dealing with multiple data series),
    # we apply a `tryCatch` to return a baranyi-object with `NULL` values rather than stopping execution.
    # The `callbackError=stop` parameter can be used to make execution blocking.
    {
      # Baranyi
      nls.args = modifyList(list(), nls.args)
      reg = do.call(stats::nls, append(
        list(formula = log(y) ~ log(.baranyi.formula(N0, Nmax, mu, lambda, base = NULL)(x)),
             start = start,
             lower = lower,
             upper = upper,
             algorithm = "port",
             data = data.frame(x = x, y = y)),
        nls.args[!(names(nls.args) %in% c("formula", "algorithm", "start", "lower", "upper"))])
      )

      baranyi.out$reg <- reg
      baranyi.out$coefficients <- as.list(summary(reg)$coefficients[,1])
      baranyi.out$isValid <- TRUE
    },
    error = function(e) {
      if (!is.null(callbackError)) {
        callbackError(e)
      }

      baranyi.out$message <<- e$message
    }
  )

  baranyi.out
}



#' Baranyi equation.
#'
#' @param N0 initial population.
#' @param Nmax final/maximum population.
#' @param mu growth rate.
#' @param lambda latency time.
#' @param base the logarithm base used for plot y-scaling. By default, the natural logarithm is used. Set `NULL` to not scale.
#'
#' @returns a function taking as input x (the time) and outputting the value of the Baranyi equation.
#'
#' @export
#'
#' @details
#' The output result is by default in the form `ln(N_t/N0)` (with `N_t` the population at time `t`).
#' The base used can be modified by specifying the desired base in the `base` argument.
#' For example, specifying `base=10` corresponds to output in the form `log_{10}(N_t/N0)`.
#' It is possible to specify `base = NULL` to retrieve the normal `N_t` output.
#'
#' @examples
#' f <- .baranyi.formula(0.1, 2, 0.2, 5)
#' f(4)
#' ## [1] 0.3498583
#' f(20)
#' ## [1] 2.344923
.baranyi.formula = function(N0, Nmax, mu, lambda, base = exp(1)) {
  if (is.null(base)) {
    return(
      function(x) {
        A = x + 1/mu * log(exp(-mu * x) + exp(-mu * lambda) - exp(-mu*(x+lambda)))
        return( exp(log(N0) + mu*A - log(1+( (exp(mu*A)-1)/(Nmax/N0) ))) )
      }
    )
  }
  return(
    function(x) {
      A = x + 1/mu * log(exp(-mu * x) + exp(-mu * lambda) - exp(-mu*(x+lambda)))
      return( (mu*A - log(1+( (exp(mu*A)-1)/(Nmax/N0) ))) * log(exp(1), base = base) )
    }
  )
}



#' Baranyi class
#'
#' Test if a variable or variable list is/are baranyi-object(s).
#'
#' @param x variable or list.
#'
#' @return `TRUE` if the object or all objects are of class `baranyi`.
#'
#' @export
#'
#' @examples
#' # TRUE return
#' r1 <- MicrobialGrowth(example_data$time, example_data$y1, model="baranyi")
#' r2 <- MicrobialGrowth.create(N0 = 0.14, Nmax = 1.43, mu = 0.07, lambda = 45,
#'                                    xlim = c(0, 100), model="baranyi")
#' is.baranyi(r1)
#' is.baranyi(r2)
#' is.baranyi(c(r1, r2))
#' is.baranyi(list(r1,r2))
#'
#' # FALSE return
#' is.baranyi(1)
#' is.baranyi(list())
#' is.baranyi(c(r1, r2, 1))
#' is.baranyi(list(r1, r2, 1))
is.baranyi <- function(x) .is.baranyi(x)
.is.baranyi <- function(x) UseMethod(".is.baranyi")
.is.baranyi.default <- function(x) "baranyi" %in% class(x)
.is.baranyi.baranyi <- function(x) TRUE
.is.baranyi.list <- function(x) length(x) && all(sapply(x, is.baranyi))



#' Baranyi create function
#'
#' Baranyi-object creator from the 4 biological meaning parameters.
#'
#' @param N0 initial population.
#' @param Nmax final/maximum population.
#' @param mu growth rate.
#' @param lambda latency time.
#' @param xlim range of values to simulate `x` and `y` (and hence plotting, *etc.*)
#' @param n number of points to simulate in the interval `xlim`.
#'
#' @returns a Baranyi-object composed of
#' \item{call}{the matched call with several components.}
#' \item{coefficients}{coefficients obtained by regression.}
#' \item{data}{data used for regression, once the y values are clipped}
#' \item{f}{a list of functions such as `formula` to retrieve the function of the model with the coefficients obtained by regression, `confint` to retrieve the confidence intervals, etc.}
#' \item{isValid}{a boolean indicating whether the regression was successful or not.}
#' \item{message}{always with this method.}
#' \item{reg}{always with this method.}
#'
#' @include utils.R
#'
#' @seealso \link{MicrobialGrowth.create}
baranyi.create = function(N0, Nmax, mu, lambda, xlim, n = 101) {
  baranyi.out = .new.baranyi.core()
  baranyi.out$call = match.call()

  N0 = .parseMicrobialGrowthCreateArgs(N0)
  Nmax = .parseMicrobialGrowthCreateArgs(Nmax)
  mu = .parseMicrobialGrowthCreateArgs(mu)
  lambda = .parseMicrobialGrowthCreateArgs(lambda)

  .checkMicrobialGrowthCreateArgs(N0=N0, Nmax=Nmax, mu=mu, lambda=lambda, xlim=xlim, n=n)

  baranyi.out$data$x = seq(xlim[1], xlim[2], length.out = n)
  baranyi.out$coefficients = list(N0 = N0[2], Nmax = Nmax[2], mu = mu[2], lambda = lambda[2])

  baranyi.out$data$y = .baranyi.formula(N0[2], Nmax[2], mu[2], lambda[2], base = NULL)(baranyi.out$data$x)

  .confint = matrix(c(N0[1], N0[3], Nmax[1], Nmax[3], mu[1], mu[3], lambda[1], lambda[3]),
                    byrow = TRUE, ncol = 2,
                    dimnames = list(c("N0", "Nmax", "mu", "lambda"), c("min", "max")))

  baranyi.out$f$confint = function(level = NULL) {
    if (!is.null(level)) { message("In the context of a user-created baranyi-object, the `level` argument is ignored") }
    return(.confint)
  }

  default.doublingTime = baranyi.out$f$doublingTime
  baranyi.out$f$doublingTime = function(level = NULL) {
    default.doublingTime(level)
  }

  baranyi.out$isValid = TRUE

  baranyi.out
}