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


#' Rosso object
#'
#' A MicrobialGrowth object specialized for the rosso 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 Rosso 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.rosso.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.rosso.core = function(...) {
  env = .new.MicrobialGrowth.core(...)
  class(env) = c("rosso", 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["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).")
    }
    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.")
    }

    cond = do.call(getFormula("rosso"), coefficients<-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 = cond(x)

      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).")
    }

    cond = do.call(getFormula("rosso"), coefficients<-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 = cond(x)

      return(y)
    }

    return(f)
  }

  env$f$doublingTime = function(level = 0.95) {
    # TODO : To rework or validate for Rosso -----------------------------------
    #
    # # TODO : To check N0 sur le confint ? => risque de N0 négatif et donc de faire échouer log(Nmax/N0)
    # if (!rosso.out$isValid) stop("Coefficient mu not defined, impossible to calculate the doubling time")
    #
    # confint = rosso.out$f$confint(level)
    # N0 = append(confint["N0",], c(coefficient=rosso.out$coefficient$N0), 1)
    # Nmax = append(confint["Nmax",], c(coefficient=rosso.out$coefficient$Nmax), 1)
    # mu = append(confint["mu",], c(coefficient=rosso.out$coefficient$mu), 1)
    # lambda = append(confint["lambda",], c(coefficient=rosso.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)) {
        if (y <= log(env$coefficients$N0)) {
          warning("NaNs produced")
          return(NaN)
        }
        return( env$coefficients$lambda - ((log( (exp(-log(y/env$coefficients$Nmax))-1)/((env$coefficients$Nmax/env$coefficients$N0)-1) ))/env$coefficients$mu) )
      }
      else {
        if (y <= 0) {
          warning("NaNs produced")
          return(NaN)
        }
        return( env$coefficients$lambda - ((log( (exp(log(env$coefficients$Nmax/env$coefficients$N0)-y*log(base))-1)/((env$coefficients$Nmax/env$coefficients$N0)-1) ))/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("rosso"), append(env$coefficients, list(base = base))) )
  }

  env
}


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

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

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

  rosso.out
}



#' Rosso 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 Rosso 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 <- .rosso.formula(0.1, 2, 0.2, 5)
#' f(4)
#' ## [1] 0
#' f(20)
#' ## [1] 2.32998
.rosso.formula = function(N0, Nmax, mu, lambda, base = exp(1)) {
  if (is.null(base)) {
    return(
      function(x) {
        y = rep(NA, length(x))

        y[x <= lambda] = N0
        y[x > lambda] = exp( log(Nmax) - log( 1+((Nmax/N0)-1)*exp(-mu*(x[x > lambda]-lambda)) ) )

        return (y)
      }
    )
  }
  return(
    function(x) {
      y = rep(NA, length(x))

      y[x <= lambda] = 0
      y[x > lambda] = (log(Nmax/N0) - log( 1+((Nmax/N0)-1)*exp(-mu*(x[x > lambda]-lambda)) )) * log(exp(1), base = base)

      return (y)
    }
  )
}



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



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

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

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

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

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

  rosso.out$isValid = TRUE

  rosso.out
}