R/gompertz.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/>


#' Gompertz object
#'
#' A MicrobialGrowth object specialized for the gompertz 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 Gompertz 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.gompertz.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.gompertz.core = function(...) {
  env = .new.MicrobialGrowth.core(...)
  class(env) = c("gompertz", 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["Nmax",1] <= conf["N0",2]) {
      mean_conf = mean(c(env$coefficients$N0, env$coefficients$Nmax))
      noise = 0.01 * (min(conf["Nmax",2]-conf["Nmax",1], conf["N0",2]-conf["N0",1]))
      conf["N0",2]  = mean_conf - noise
      conf["Nmax",1] = mean_conf + noise
      warning("Lower bound of the confidence interval for Nmax overlaps the upper bound of the confidence interval for N0. The values have been replaced by the average between Nmax and N0, and the addition or subtraction of a noise, respectively for Nmax and N0, equal to 1% of the minimum deviation of the limits of N0 or Nmax.")
    }
    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("gompertz"), coefficients1<-list(N0=conf["N0",2], Nmax=conf["Nmax",1], mu=conf["mu",2], lambda=conf["lambda",2], base = base))
    cond2 = do.call(getFormula("gompertz"), coefficients2<-list(N0=conf["N0",2], Nmax=conf["Nmax",1], mu=conf["mu",1], lambda=conf["lambda",2], base = base))

    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]))])

      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["N0",1] <= 0) {
      conf["N0",1]  = 1/.Machine$double.xmax
      warning("Lower bound of the confidence interval for N0 was below or equal to 0. Its value was replaced by a value close to 0 (=1/.Machine$double.xmax).")
    }
    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("gompertz"), coefficients1<-list(N0=conf["N0",1], Nmax=conf["Nmax",2], mu=conf["mu",1], lambda=conf["lambda",1], base = base))
    cond2 = do.call(getFormula("gompertz"), coefficients2<-list(N0=conf["N0",1], Nmax=conf["Nmax",2], mu=conf["mu",2], lambda=conf["lambda",1], base = base))

    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]))])

      return(y)
    }

    return(f)
  }

  env$f$doublingTime = function(level = 0.95) {
    # TODO : To check N0 sur le confint ? => risque de N0 négatif et donc de faire échouer log(Nmax/N0)
    if (!env$isValid) stop("Coefficient mu not defined, impossible to calculate the doubling time")

    confint = env$f$confint(level)
    N0 = append(confint["N0",], c(coefficient=env$coefficient$N0), 1)
    Nmax = append(confint["Nmax",], c(coefficient=env$coefficient$Nmax), 1)
    mu = append(confint["mu",], c(coefficient=env$coefficient$mu), 1)
    lambda = append(confint["lambda",], c(coefficient=env$coefficient$lambda), 1)

    doublingTime = log(2)/mu
    A = log(Nmax/N0)
    start = lambda - (A*log(-log(log(10)/A)) - A)/(mu * exp(1)) # En fait, c'est pas le start du doublingTime mais plutôt la première fois où on multiplie N0 par n=2 (:log(2))

    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)) {
        A = log(env$coefficients$Nmax/env$coefficients$N0)
        return( env$coefficients$lambda - ((log(-log((log(y/env$coefficients$N0))/A))-1) / (env$coefficients$mu*exp(1)/A)) )
      }
      else {
        A = log(env$coefficients$Nmax/env$coefficients$N0)
        return( env$coefficients$lambda - ((log(-log((y*log(base))/A))-1) / (env$coefficients$mu*exp(1)/A)) )
      }
    }
  }

  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("gompertz"), append(env$coefficients, list(base = base))) )
  }

  env
}



#' Gompertz regression function
#'
#' Regression function for Gompertz'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{.gompertz.formula}
.MicrobialGrowth.gompertz = 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
  gompertz.out = .new.gompertz.core()
  gompertz.out$data = list(x=x, y=y)
  gompertz.out$call = match.call()
  gompertz.out$call$clip=clip
  gompertz.out$call$start=start
  gompertz.out$call$lower=lower
  gompertz.out$call$upper=upper
  gompertz.out$call$nls.args=nls.args
  gompertz.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 gompertz-object with `NULL` values rather than stopping execution.
    # The `callbackError=stop` parameter can be used to make execution blocking.
    {
      # Gompertz
      nls.args = modifyList(list(), nls.args)
      reg = do.call(stats::nls, append(
        list(formula = log(y) ~ log(N0) + (log(Nmax)-log(N0)) * exp(-exp(((mu*exp(1))/(log(Nmax)-log(N0))) * (lambda-x) + 1 )),
             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"))])
      )

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

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

  gompertz.out
}



#' Gompertz 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.
#'
#' @return a function taking as input x (the time) and outputting the value of the Gompertz 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 <- .gompertz.formula(0.1, 2, 0.2, 5)
#' f(4)
#' ## [1] 0.1150952
#' f(20)
#' ## [1] 2.505549
.gompertz.formula = function(N0, Nmax, mu, lambda, base = exp(1)) {
  if (is.null(base)) {
    return(
      function(x) {
        A = log(Nmax/N0)
        return ({ exp(A * exp( -exp( (mu*exp(1)/A) * (lambda-x) + 1 ) )) * N0 })
      }
    )
  }
  return(
    function(x) {
      A = log(Nmax/N0)
      return ({ (A * exp( -exp( (mu*exp(1)/A) * (lambda-x) + 1 ) )) * log(exp(1), base) })
    }
  )
}



#' Gompertz print function
#'
#' Default print.MicrobialGrowth function can be overwritten with the following function
#'
#' @param x gompertz-object.
#' @param ... further arguments passed to or from other methods.
#'
#' @return No return value, called to print information about a MicrobialGrowth-object based on the Gompertz model.
#'
#' @exportS3Method
#'
#' @seealso \link{print.MicrobialGrowth}
print.gompertz = function(x, ...) {
  print.MicrobialGrowth(x, ...)
}



#' Gompertz plot function
#'
#' Default plot.MicrobialGrowth function can be overwritten with the following function
#'
#' @param x gompertz-object.
#' @param ... further arguments passed to or from other methods.
#'
#' @return No return value, called to plot a MicrobialGrowth-object based on the Gompertz model.
#'
#' @exportS3Method
#'
#' @seealso \link{plot.MicrobialGrowth}
plot.gompertz = function(x, ...) {
  plot.MicrobialGrowth(x, ...)
}



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



#' Gompertz create function
#'
#' Gompertz-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 Gompertz-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}
gompertz.create = function(N0, Nmax, mu, lambda, xlim, n = 101) {
  gompertz.out = .new.gompertz.core()
  gompertz.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)

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

  gompertz.out$data$y = .gompertz.formula(N0[2], Nmax[2], mu[2], lambda[2], base = NULL)(gompertz.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")))

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

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

  gompertz.out$isValid = TRUE

  gompertz.out
}


#' Graphical Example of Modified Gompertz Equation.
#'
#' @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`.
#'
#' @return No return value, called to plot a MicrobialGrowth object with the Gompertz model to illustrate the different coefficients.
#'
#' @importFrom graphics abline arrows text
#' @export
#'
#' @examples
#' gompertz.explain()
#'
#' gompertz.explain(0.15, 2, 0.1, 40, c(0,100))
gompertz.explain = function(N0 = 0.14, Nmax = 1.43, mu = 0.07, lambda = 40, xlim = c(0, 100)) {
  g = MicrobialGrowth.create(N0, Nmax, mu, lambda, xlim)
  A = log(Nmax/N0)

  plot(g, col = NULL, reg.args = list(col = "black"), ylim = c(log(1), A), xlab ="Time t", main = "Curve example")

  abline(h = A, lty = 2, col = "red")
  text(x = 0, y = A, adj = c(0,1.2), label = expression(ln(N[max] / N[0])), col = "red")
  abline(h = log(N0/N0), lty = 2, col = "red")
  text(x = 0, y = log(N0/N0), adj = c(0,-0.2), label = expression(ln(N[0] / N[0])), col = "red")

  arrows(x0 = 0, y0 = 0.15*A, x1 = lambda, y1 = 0.15*A, lty = 2, col = "red", code = 3, angle = 90, length = 0.1)
  text(x = lambda/2, y = 0.15*A, adj = c(0.5, -0.5), label = expression(lambda), col = "red")

  ti = ((mu*exp(1)/A)*lambda+1) / (mu*exp(1)/A)
  abline(g$f$intersection(x=ti)-mu*ti, mu, lty = 2, col = "red")
  text(x = ti, y = g$f$intersection(x=ti), adj = c(1.2,-0.2), label = expression(mu[m]), col = "red")
}