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


#' Linear object
#'
#' A MicrobialGrowth object specialized for the linear 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 linear 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.linear.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.linear.core = function(...) {
  env = .new.MicrobialGrowth.core(...)
  class(env) = c("linear", 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 = NULL) {
    if (!is.null(base)) { stop("confint.lower is only supported 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).")
    }

    cond = do.call(getFormula("linear"), coefficients<-list(N0=NULL, Nmax=NULL, 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 only supported for base=NULL") }

    conf = env$f$confint(level)
    cond = do.call(getFormula("linear"), coefficients<-list(N0=NULL, Nmax=NULL, 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) {
    # le doubling time n a pas de sens sur une moisissure dont on regarde le rayon
    stop("Doubling time cannot be calculated for linear model")
  }

  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( y/env$coefficients$mu + env$coefficients$lambda )
      }
      else {
        if (y <= 0) {
          warning("NaNs produced")
          return(NaN)
        }
        return(NaN)
      }
    }
  }

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

  env
}


#' Linear regression function
#'
#' Regression function for Linear'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{.linear.formula}
.MicrobialGrowth.linear = 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)
  clip[1] = 0

  { # 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.linear(x, y, start, lower, upper)
  start = default.values$start
  lower = default.values$lower
  upper = default.values$upper

  # Creation of the return value
  linear.out = .new.linear.core()
  linear.out$data = list(x=x, y=y)
  linear.out$call = match.call()
  linear.out$call$clip=clip
  linear.out$call$start=start
  linear.out$call$lower=lower
  linear.out$call$upper=upper
  linear.out$call$nls.args=nls.args
  linear.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 linear-object with `NULL` values rather than stopping execution.
    # The `callbackError=stop` parameter can be used to make execution blocking.
    {
      # Linear
      nls.args = modifyList(list(), nls.args)
      reg = do.call(stats::nls, append(
        list(formula = y ~ .linear.formula(N0=NULL, Nmax=NULL, 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"))])
      )

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

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

  linear.out
}


.getDefaultNlsValues.linear = function(x, y, start = list(), lower = list(), upper = list()) {
  .start = start
  startN0 = min(y)
  start = list()
  start$lambda = if("lambda" %in% names(.start)) .start[["lambda"]] else x[tail(which(y <= startN0), 1)]
  start$mu = if("mu" %in% names(.start)) .start[["mu"]] else tryCatch(
      {
        cond = x >= start$lambda
        r = lm(y[cond] ~ x[cond])$coefficients[[2]]
        if (is.na(r) | is.null(r)) stop("lm fail")
        r
      },
      error = function(e) { return(max( (y[2:length(y)] - y[1:(length(y)-1)]) / (x[2:length(x)] - x[1:(length(x)-1)]) )) }
    )

  lower = list(
    mu = if("mu" %in% names(lower)) lower[["mu"]] else ((max(y)-min(y))/(max(x)-min(x))),
    lambda = if("lambda" %in% names(lower)) lower[["lambda"]] else 0
  )

  upper = list(
    mu = if("mu" %in% names(upper)) upper[["mu"]] else max( (y[2:length(y)] - y[1:(length(y)-1)]) / (x[2:length(x)] - x[1:(length(x)-1)])),
    lambda = if("lambda" %in% names(upper)) upper[["lambda"]] else x[which(y > startN0)[1]]    #max(x)
  )

  if (start$mu < lower$mu) start$mu = lower$mu
  if (start$mu > upper$mu) start$mu = upper$mu

  return(list(start=start[c( "mu", "lambda")], lower=lower[c("mu", "lambda")], upper=upper[c("mu", "lambda")]))
}


#' Linear equation.
#'
#' @param N0 initial radius.
#' @param Nmax final/maximum radius.
#' @param mu growth rate.
#' @param lambda latency time.
#' @param base decimal base used for plot y-scaling.
#'
#' @return a function taking as input x (the time) and outputting the value of the linear equation.
#'
#' @export
#'
#' @examples
#' f <- .linear.formula(0.1, 2, 0.2, 5)
#' f(4)
#' ## [1] 0
#' f(20)
#' ## [1] 3
.linear.formula = function(N0, Nmax, mu, lambda, base = NULL) {
  if (is.null(base)) {
    return(
      function(x) {
        y = rep(NA, length(x))

        y[x <= lambda] = 0
        y[x > lambda] = mu*(x[x > lambda]-lambda)

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

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

      return (y)
    }
  )
}


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



#' Linear create function
#'
#' Linear-object creator from the 4 biological meaning parameters.
#'
#' @param N0 initial radius.
#' @param Nmax final/maximum radius.
#' @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 Linear-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}
linear.create = function(N0, Nmax, mu, lambda, xlim, n = 101) {
  linear.out = .new.linear.core()
  linear.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)

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

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

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

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

  linear.out$isValid = TRUE

  linear.out
}


#' Linear plot function
#'
#' Default plot.MicrobialGrowth function can be overwritten with the following function
#'
#' @param x linear-object.
#' @param base base used for plot y-scaling.
#' @param ... further arguments passed to or from other methods.
#'
#' @return No return value, called to plot a MicrobialGrowth-object based on the linear model.
#'
#' @exportS3Method
#'
#' @seealso \link{plot.MicrobialGrowth}
plot.linear = function(x, base=NULL, ...) {
  if (!is.null(base)) {stop('linear model only supports base=NULL')}
  plot.MicrobialGrowth(x, base=base, ...)
}