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


#' MicrobialGrowth regression function
#'
#' Regression function to different microbial growth models.
#'
#' @param x index series or time series.
#' @param y values or list of values to regress (should not be logged).
#' @param model wanted growth model : "baranyi", "gompertz" or "rosso".
#' @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.
#' @param ... further arguments passed to or from other methods.
#'
#' @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
#' Use `listAvailableModels()` function to see all values accepted by `model` parameter.
#'
#' 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.
#'
#' @examples
#' # Using the embedded data example_data
#' # Simple example
#' g <- MicrobialGrowth(example_data$time, example_data$y1, model="gompertz")
#'
#' # Multiple regression example
#' G <- MicrobialGrowth(example_data$time, example_data[2:ncol(example_data)], model="gompertz")
#'
#' # Example of multiple parameter changes
#' g <- MicrobialGrowth(example_data$time, example_data$y1, model="gompertz",
#'                            clip = c(0.15, Inf), start = list(N0=0.1, Nmax=2,
#'                            mu=0.05, lambda=40), lower = list(lambda = 40))
#'
#' # Example of using `nls.args` to apply weight to some data
#' g <- MicrobialGrowth(example_data$time, example_data$y1, model="gompertz",
#' nls.args = list(weights = (function(x){(x >= 50 & x <= 70)*9 + 1})(example_data$time)))
#'
#' # Example of callbackError (remaining non-blocking)
#' g <- MicrobialGrowth(example_data$time, example_data$y15, model="gompertz",
#'                            callbackError = warning)
#'
#' # Example of callbackError (becoming blocking)
#' try(
#'   g <- MicrobialGrowth(c(1,2,3,4,5),c(1,1,1,1,1), model="gompertz", callbackError = stop)
#' )
MicrobialGrowth = function(x, y, model = "gompertz", clip = c(-Inf, Inf), start = list(), lower = list(), upper = list(), nls.args = list(), callbackError = NULL, ...) {
  functionName = getFunctionName(model)
  if (!exists(functionName)) {
    stop(sprintf("Model %s not found. Available models are : %s", model, paste(listAvailableModels(), collapse = ", ")))
  }

  doCall = function(y) {
    # TODO : Tester de normaliser systématique les valeurs entre 0 et 1 (?) *
    #        Puis réélever les valeurs en fonction de la division/soustraction effectée
    #        Par contre, on aura une certaine disonance entre le `reg` et les données en accès direct,
    #        Le premier étant sur les nombres relatifs et le second retravaillé -> nécessite de prendre en
    #        compte ça dans toutes les fonctions (confint, formula, etc.) !?
    #
    #       * on a pu voir sur des données dont les ordres de grandeurs sont élevés que la régression était
    #         difficile (échec) mais passait régulièrement en appliquant une division de 10E6 par exemple.
    r = do.call(functionName,
                append(list(x=x, y=y, clip=clip, start=start, lower=lower, upper=upper, nls.args=nls.args, callbackError=callbackError), list(...)))
    r
  }

  if (is.list(y)) { lapply(y, function(y2) { doCall(y2) }) }
  else { doCall(y) }
}





#' MicrobialGrowth print function
#'
#' Print function of MicrobialGrowth-objects.
#'
#' @param x MicrobialGrowth-object.
#' @param ... further arguments passed to or from other methods.
#'
#' @return No return value, called to print information about a MicrobialGrowth-object.
#'
#' @exportS3Method
#'
#' @seealso \link{MicrobialGrowth}, \link{MicrobialGrowth.create}
#'
#' @examples
#' # Print from regressed MicrobialGrowth-object
#' g <- MicrobialGrowth(example_data$time, example_data$y1, model="gompertz")
#' print(g) # or just `g` in the console
#'
#' # Print from a user-created MicrobialGrowth-object (via MicrobialGrowth.create)
#' g <- MicrobialGrowth.create(N0 = 0.14, Nmax = 1.43, mu = 0.07, lambda = 45,
#'                                   xlim = c(0, 100), model="gompertz")
#' print(g) # or just `g` in the console
print.MicrobialGrowth = function(x, ...) {
  cat(sprintf("MicrobialGrowth, model %s:", getModelName(x)), "\n")
  print(unlist(x$coefficients), ...)
}


#' MicrobialGrowth summary function
#'
#' Summarizes the regression of an MicrobialGrowth-object.
#'
#' @param object MicrobialGrowth-object.
#' @param ... additional arguments affecting the summary produced.
#'
#' @return The summary of the successful regression, `NULL` otherwise.
#'
#' @include utils.R
#' @exportS3Method
#'
#' @details
#' Equivalent to `summary(MicrobialGrowthObject$reg, ...)` to which we add the corresponding `model` member and the `summary.MicrobialGrowth` class.
#'
#' @examples
#' # Simple example
#' g <- MicrobialGrowth(example_data$time, example_data$y1)
#' summary(g)
#'
#' # Example without summary available
#' g <- MicrobialGrowth(example_data$time, example_data$y15)
#' summary(g)
#'
#' g <- MicrobialGrowth.create(0.14, 1.5, 0.07, 45, c(0,100), model="gompertz")
#' summary(g)
summary.MicrobialGrowth = function(object, ...) {
  if (!is.null(object$reg)) {
    s = summary(object$reg, ...)
    s$model = getModelName(object)
    class(s) = c("summary.MicrobialGrowth", class(s))
    return(s)
  }
  return(NULL)
}




#' MicrobialGrowth plot function
#'
#' Plot function of MicrobialGrowth-objects.
#'
#' @param x MicrobialGrowth-object.
#' @param main main title for the plot.
#' @param xlab title for the x axis.
#' @param ylab title for the y axis.
#' @param n the number of x values at which to evaluate. See details section.
#' @param base the logarithm base used for plot y-scaling. By default, the natural logarithm is used. Set `NULL` to not scale.
#' @param display.coefficients boolean indicating the display or not of the values of coefficients (under the main title).
#' @param display.model boolean indicating the model used for regression (on right side).
#' @param display.confint boolean indicating the display or not of confidence intervals (in the form of curves and area).
#' @param reg.args customization parameters of the curve obtained by regression (see \link[graphics]{curve} for possible parameters).
#' @param title.args title customization parameters `main`, `xlab` and `ylab` (see \link[graphics]{title} for possible parameters).
#' @param model.args model display customization parameters (see \link[graphics]{mtext} for possible parameters).
#' @param coefficients.args coefficient display customization parameters (see \link[graphics]{mtext} for possible parameters).
#' @param confint.args parameters for customizing the plotting of curves and area, corresponding to the confidence interval (see details section).
#' @param ... other graphical parameters (see \link[base]{plot}).
#'
#' @return No return value, called to plot a MicrobialGrowth-object.
#'
#' @include utils.R
#' @import graphics
#' @importFrom utils modifyList
#' @importFrom grDevices col2rgb
#' @exportS3Method
#'
#' @details
#' Similar to the `curve` function, the `n` argument corresponds to the number of points evaluated to draw the curves of regression, confidence bounds and the associated area. Increase its value for a more accurate representation.
#'
#' When `base` is not `NULL`, the plot produced is \eqn{log_n(N/N0)}, where n is the value specified in the `base` argument.
#'
#' @examples
#' # Example plot of a MicrobialGrowth-object obtained by regression
#' g <- MicrobialGrowth(example_data$time, example_data$y1, model="gompertz")
#' plot(g)
#'
#' # Example plot of a user-created MicrobialGrowth-object (via MicrobialGrowth.create)
#' g <- MicrobialGrowth.create(N0 = 0.14, Nmax = 1.43, mu = 0.07, lambda = 45,
#'                                   xlim = c(0, 100), model="gompertz")
#' plot(g)
#'
#' # Example plot with usual graphical parameters
#' plot(g, pch = 4, cex = 2, col = "blue", xlab = "Time (hours)", main = "Gompertz regression")
#'
#' # Example of plot hiding the coefficients and customizing the curve obtained by regression
#' plot(g, display.coefficients = FALSE, reg.args = list(col = "green", lty = 2, lwd = 5))
#'
#' # Example of a plot displaying the curves and area of the confidence interval
#' g <- MicrobialGrowth(example_data$time, example_data$y1, model="gompertz")
#' plot(g, display.confint = TRUE)
#'
#' # Example of a plot customizing the confidence interval
#' plot(g, xlim = c(80, 100), ylim = c(1.8, 2.4), # Zoom in to see the example better
#' display.confint = TRUE,
#' confint.args = list(
#'   lines = list(col = "purple", lty = 2, lwd = 2),
#'   area = list(col = "green", opacity = 0.1)
#' ))
#'
#' # Example of a plot customizing the display of coefficients and titles
#' plot(g, main = "Gompertz",
#' coefficients.args = list(cex = 1.5, side = 4, line = 1),
#' title.args = list(col.main = "blue", col.lab = "red"))
plot.MicrobialGrowth = function(x, main = NULL, xlab = NULL, ylab = NULL, n = 101, base = exp(1), display.coefficients = TRUE, display.model = TRUE, display.confint = FALSE, reg.args = list(col = "red"), title.args = list(line = 2), model.args = list(side = 4, line = 0), coefficients.args = list(cex = 0.9, line = 0.2, side = 3), confint.args = list(), ...) {
  model = getModelName(x)

  if (is.null(main)) { main = "" }
  if (is.null(xlab)) { xlab = "Time" }
  if (is.null(ylab)) { ylab = if(x$isValid & !is.null(base)) (if(base == exp(1)) expression(ln(N/N[0])) else bquote(~log[.(base)]*"("*N*"/"*N[0]*")")) else "y (clipped)" }

  # If regression failed
  if (!x$isValid) {
    plot(x$data$x, x$data$y, main = "", xlab = "", ylab = "", ...)

    if (display.coefficients) {
      coefficients.args = modifyList(list(cex = 0.9, line = 0.2, side = 3, col = "red"), coefficients.args)
      do.call(graphics::mtext, append(
        list(text = "Regression failed"),
        coefficients.args[!(names(coefficients.args) %in% c("text"))])
      )
    }
  }
  # Else, regression successful
  else {
    plot(x$data$x, if (is.null(base)) x$data$y else log(x$data$y/x$coefficients$N0, base=base), main = "", xlab = "", ylab = "", ...)

    reg.args = modifyList(list(col = "red", xlim = c(max(c(0, par("usr")[1])), par("usr")[2])), reg.args)
    expr = function(t){getFormula(model)(x$coefficients$N0, x$coefficients$Nmax, x$coefficients$mu, x$coefficients$lambda, base = base)(t)}
    do.call(graphics::curve, append(
      list(expr = substitute(expr, parent.frame()), add = TRUE, n = n),
      reg.args[!(names(reg.args) %in% c("expr", "add"))])
    )

    if (display.confint) suppressMessages({
      confint.level = 0.95
      if ("level" %in% names(confint.args)) {
        if (is.numeric(confint.args$level) && (confint.args$level >= 0) && (confint.args$level < 1)) { confint.level = confint.args$level }
        else { stop("Sub-parameter \"level\" in \"confint.args\" parameter must be a numeric in [0,1[") }
      }


      # Simulates x values
      x.values = seq(par("usr")[1], par("usr")[2], length.out = n)
      x.values = x.values[x.values >= 0]

      f_inf = x$f$confint.lower(confint.level, base = base)(x.values)
      f_sup = x$f$confint.upper(confint.level, base = base)(x.values)

      # Managing parameters for the confidence interval
      confint.lines.args = modifyList(reg.args, if("lines" %in% names(confint.args)) confint.args$lines else list())
      area.opacity = if("area" %in% names(confint.args) & "opacity" %in% names(confint.args$area)) confint.args$area$opacity else 1/3
      area.color = if("area" %in% names(confint.args) & "col" %in% names(confint.args$area)) confint.args$area$col else reg.args$col
      confint.area.args = modifyList(modifyList(modifyList(reg.args, list(border = NA)), if("area" %in% names(confint.args)) confint.args$area else list()), list(col = do.call(grDevices::rgb, as.list(((col2rgb(area.color, TRUE) * c(1, 1, 1, area.opacity))/255)[,1]))))

      # Confidence interval plot
      do.call(graphics::polygon, append(
        list(x = c(x.values, rev(x.values)), y = c(f_inf, rev(f_sup))),
        confint.area.args[!(names(confint.area.args) %in% c("x", "y", "opacity"))])
      )
      # do.call(graphics::lines, append(
      #   list(x = x.values, y = f_inf(x.values)),
      #   confint.lines.args[!(names(confint.lines.args) %in% c("x", "y"))])
      # )
      do.call(graphics::lines, append(
        list(x = x.values, y = f_inf),
        confint.lines.args[!(names(confint.lines.args) %in% c("x", "y"))])#TODO TODELETE coloration de la borne inf pour debug
      )
      # do.call(graphics::lines, append(
      #   list(x = x.values, y = f_sup(x.values)),
      #   confint.lines.args[!(names(confint.lines.args) %in% c("x", "y"))])
      # )
      do.call(graphics::lines, append(
        list(x = x.values, y = f_sup),
        confint.lines.args[!(names(confint.lines.args) %in% c("x", "y"))])#TODO TODELETE coloration de la borne inf pour debug
      )
    })

    if (display.coefficients) {
      coefficients.args = modifyList(list(cex = 0.9, line = 0.2, side = 3), coefficients.args)
      r = if("round" %in% names(coefficients.args) && is.numeric(coefficients.args$round)) coefficients.args$round else 4
      do.call(graphics::mtext, append(
        list(text = bquote(N[0] ~ "=" ~ .(round(x$coefficients$N0,r)) * "," ~ N[max] ~ "=" ~ .(round(x$coefficients$Nmax,r)) * "," ~ mu[m] ~ "=" ~ .(round(x$coefficients$mu,r)) ~ "and" ~ lambda ~ "=" ~ .(round(x$coefficients$lambda,r)))),
        coefficients.args[!(names(coefficients.args) %in% c("text", "round"))])
      )
    }
  }

  title.args = modifyList(list(line = 2), title.args)
  do.call(graphics::title, append(
    list(main = main, xlab = xlab, ylab = ylab),
    title.args[!(names(title.args) %in% c("main", "xlab", "ylab"))])
  )

  if (display.model) {
    model.args = modifyList(list(side = 4, line = 0), model.args)
    do.call(graphics::mtext, append(
      list(text = sprintf("Model: %s", model)),
      model.args[!(names(model.args) %in% c("text"))])
    )
  }

}


#' MicrobialGrowth class
#'
#' Test if a variable or variable list is/are MicrobialGrowth-object(s).
#'
#' @param x variable or list.
#'
#' @return `TRUE` if the object or all objects are of class `MicrobialGrowth`.
#'
#' @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.MicrobialGrowth(g1)
#' is.MicrobialGrowth(g2)
#' is.MicrobialGrowth(list(g1,g2))
#'
#' # FALSE return
#' is.MicrobialGrowth(1)
#' is.MicrobialGrowth(list())
#' is.MicrobialGrowth(c(g1, g2))
#' is.MicrobialGrowth(list(g1, g2, 1))
is.MicrobialGrowth <- function(x) {
  if ("list" %in% class(x)) { return(length(x) && all(sapply(x, is.MicrobialGrowth))) }
  return("MicrobialGrowth" %in% class(x))
}


#' MicrobialGrowth create function
#'
#' MicrobialGrowth-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 model wanted growth model: "baranyi", "gompertz" or "rosso"
#' @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`.
#' @param ... further arguments passed to or from other methods.
#'
#' @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}{always with this method.}
#' \item{reg}{always with this method.}
#'
#' @include utils.R
#' @export
#'
#' @details
#' The `N0`, `Nmax`, `mu` and `lambda` parameter-coefficients can be given as one to three values.
#'
#' A single value means that the coefficient and the confidence interval values will be identical.
#'
#' Two values means that they will correspond to the confidence interval, and the coefficient will be calculated as the average of these two values.
#'
#' Three values means that each of these values will be associated with the coefficient or the confidence interval.
#'
#' Values are always sorted automatically, which means that `c(2,1,3)` is equivalent to `c(1,2,3)`.
#'
#' @examples
#' # Simple example
#' g <- MicrobialGrowth.create(N0 = 0.14, Nmax = 1.43, mu = 0.07, lambda = 45,
#'      xlim = c(0, 100), model="gompertz")
#'
#' # Example from a regression (whose values can be stored and then reused later)
#' g <- MicrobialGrowth(example_data$time, example_data$y1, model="gompertz")
#' c <- g$coefficients
#' g2 <- MicrobialGrowth.create(c$N0, c$Nmax, c$mu, c$lambda,
#'       xlim = c(min(g$data$x),max(g$data$x)), n = length(g$data$x), model="gompertz")
#'
#' # Example with confidence intervals
#' g <- MicrobialGrowth.create(N0 = c(0.13, 0.15), Nmax = 1.43, mu = c(0.05, 0.07, 0.09),
#' lambda = c(45, 49, 43), xlim = c(0, 100), model="gompertz")
#' # Coefficient N0 is 0.14 and the confidence interval is (0.13, 0.15)
#' # Coefficient Nmax is 1.43 and the confidence interval is (1.43, 1.43)
#' # Coefficient mu is 0.07 and the confidence interval is (0.05, 0.09)
#' # Coefficient lambda is 45 and the confidence interval is (43, 49)
MicrobialGrowth.create = function(N0, Nmax, mu, lambda, xlim, model = "gompertz", n = 101, ...) {
  functionName = getCreateFunctionName(model)
  if (!exists(functionName)) {
    stop(sprintf("Model %s not found. Available models are : %s", model, paste(listAvailableModels(), collapse = ", ")))
  }

  r = do.call(functionName, list(N0=N0, Nmax=Nmax, mu=mu, lambda=lambda, xlim=xlim, n=n, ...))
  r
}