R/parametric_fits.R

In QurvE: Robust and User-Friendly Analysis of Growth and Fluorescence Curves

#' Fit nonlinear growth models to growth data
#'
#' \code{growth.gcFitModel} determines a parametric growth model that best describes the data.
#'
#' @param time Vector of the independent variable (usually time).
#' @param data Vector of dependent variable (usually growth values).
#' @param gcID (Character) The name of the analyzed sample.
#' @param control A \code{grofit.control} object created with \code{\link{growth.control}}, defining relevant fitting options.
#'
#' @return A \code{gcFitModel} object that contains physiological parameters and information about the best fit. Use \code{\link{plot.gcFitModel}} to visualize the parametric fit and growth equation.
#' \item{raw.time}{Raw time values provided to the function as \code{time}.}
#' \item{raw.data}{Raw growth data provided to the function as \code{data}.}
#' \item{gcID}{(Character) Identifies the tested sample.}
#' \item{fit.time}{Fitted time values.}
#' \item{fit.data}{Fitted growth values.}
#' \item{parameters}{List of determined growth parameters.}
#' \itemize{
#' \item \code{A}: Maximum growth.
#' \item \code{dY}: Difference in maximum growth and minimum growth of the fitted model.
#' \item \code{mu}: Maximum growth rate (i.e., maximum in first derivative of the spline).
#' \item \code{lambda}: Lag time.
#' \item \code{b.tangent}: Intersection of the tangent at the maximum growth rate with the abscissa.
#' \item \code{fitpar}: For some models: list of additional parameters used in the equations describing the growth curve.
#' \item \code{integral}: Area under the curve of the parametric fit.
#' }
#' \item{model}{(Character) The model that obtained the fit with the lowest AIC, determined by \code{\link{AIC}}.}
#' \item{nls}{\code{nls} object for the chosen model generated by the \code{\link{nls}} function.}
#' \item{reliable}{(Logical) Indicates whether the performed fit is reliable (to be set manually).}
#' \item{fitFlag}{(Logical) Indicates whether a parametric model was successfully fitted on the data.}
#' \item{control}{Object of class \code{grofit.control} containing list of options passed to the function as \code{control}.}
#'
#' @references Matthias Kahm, Guido Hasenbrink, Hella Lichtenberg-Frate, Jost Ludwig, Maik Kschischo (2010). _grofit: Fitting Biological Growth Curves with R_. Journal of Statistical Software, 33(7), 1-21. DOI: 10.18637/jss.v033.i07
#'
#' @family growth fitting functions
#'
#' @export
#'
#' @examples
#' # Create random growth dataset
#' rnd.dataset <- rdm.data(d = 35, mu = 0.8, A = 5, label = 'Test1')
#'
#' # Extract time and growth data for single sample
#' time <- rnd.dataset$time[1,]
#' data <- rnd.dataset$data[1,-(1:3)] # Remove identifier columns
#'
#' # Perform parametric fit
#' TestFit <- growth.gcFitModel(time, data, gcID = 'TestFit',
#'                  control = growth.control(fit.opt = 'm'))
#'
#' plot(TestFit, basesize = 18, eq.size = 1.5)
#'
growth.gcFitModel <- function(
    time, data, gcID = "undefined", control = growth.control()
)
    {
    # /// check input parameters
    if (methods::is(control) !=
        "grofit.control")
        stop("control must be of class grofit.control!")
    if (!any(
        c("m", "a") %in%
            control$fit.opt
    ))
        stop(
            "Fit option is not set for a model fit. See growth.control()"
        )

    # /// conversion to handle even data.frame
    # inputs
    time <- as.vector(as.numeric(as.matrix(time)))[!is.na(time)][!is.na(data)]
    data <- as.vector(as.numeric(as.matrix(data)))[!is.na(time)][!is.na(data)]

    if (length(data[data < 0]) >
        0)
        {
        data <- data + abs(min(data[data < 0])) +
            0.01  # add the absolute value of the minimum negative growth (+ 0.01) to the data
    }

    # /// check length of input data
    if (length(time) !=
        length(data))
        stop(
            "gcFitModel: length of time and data input vectors differ!"
        )
    if (max(data) <
        control$growth.thresh * data[1])
        {
        if (control$suppress.messages == FALSE)
            message(
                paste0(
                  "Parametric fit: No significant growth detected (with all values below ",
                  control$growth.thresh, " * start_value)."
              )
            )
        gcFitModel <- list(
            time.in = time, data.in = data, raw.time = time,
            raw.data = data, gcID = gcID, fit.time = NA,
            fit.data = NA, parameters = list(
                A = NA, mu = 0, tD = NA, lambda = NA,
                integral = NA, RMSE = NA
            ),
            model = NA, nls = NA, reliable = NULL,
            fitFlag = FALSE, control = control
        )
        class(gcFitModel) <- "gcFitModel"
        return(gcFitModel)
    }
    # /// check if there are enough data points
    if (length(data) <
        5)
        {
        if (control$suppress.messages == FALSE)
            message(
                "gcFitModel: There is not enough valid data. Must have at least 5 unique values!"
            )
        gcFitModel <- list(
            time.in = time, data.in = data, raw.time = time,
            raw.data = data, gcID = gcID, fit.time = NA,
            fit.data = NA, parameters = list(
                A = NA, mu = NA, tD = NA, lambda = NA,
                integral = NA, RMSE = NA
            ),
            model = NA, nls = NA, reliable = NULL,
            fitFlag = FALSE, control = control
        )
        class(gcFitModel) <- "gcFitModel"
        return(gcFitModel)
    }
    # /// check if there are enough growth
    # values above the start value
    data.growth <- data[data > data[1]]
    if (length(data.growth) <
        5)
        {
        if (control$suppress.messages == FALSE)
            message(
                "gcFitModel: No significant amount of growth values above the start value!"
            )
        gcFitModel <- list(
            time.in = time, data.in = data, raw.time = time,
            raw.data = data, gcID = gcID, fit.time = NA,
            fit.data = NA, parameters = list(
                A = NA, mu = NA, tD = NA, lambda = NA,
                integral = NA, RMSE = NA
            ),
            model = NA, nls = NA, reliable = NULL,
            fitFlag = FALSE, control = control
        )
        class(gcFitModel) <- "gcFitModel"
        return(gcFitModel)
    } else
    {
        gcFitModel <- growth.param(time, data, gcID, control)
    }
    invisible(gcFitModel)
}

#' Internal function called within \code{\link{growth.gcFitModel}} to perform growth model fitting.
#'
#' @param time Vector of the independent variable (usually time).
#' @param data Vector of dependent variable (usually growth values).
#' @param gcID (Character) The name of the analyzed sample.
#' @param control A \code{grofit.control} object created with \code{\link{growth.control}}, defining relevant fitting options.
#'
#' @return A \code{gcFitModel} object that contains physiological parameters and information about the best fit. Use \code{\link{plot.gcFitModel}} to visualize the parametric fit and growth equation.
#' \item{raw.time}{Raw time values provided to the function as \code{time}.}
#' \item{raw.data}{Raw growth data provided to the function as \code{data}.}
#' \item{gcID}{(Character) Identifies the tested sample.}
#' \item{fit.time}{Fitted time values.}
#' \item{fit.data}{Fitted growth values.}
#' \item{parameters}{List of determined growth parameters.}
#' \itemize{
#' \item \code{A}: Maximum growth.
#' \item \code{dY}: Difference in maximum growth and minimum growth of the fitted model.
#' \item \code{mu}: Maximum growth rate (i.e., maximum in first derivative of the spline).
#' \item \code{lambda}: Lag time.
#' \item \code{b.tangent}: Intersection of the tangent at the maximum growth rate with the abscissa.
#' \item \code{fitpar}: For some models: list of additional parameters used in the equations describing the growth curve.
#' \item \code{integral}: Area under the curve of the parametric fit.
#' }
#' \item{model}{(Character) The model that obtained the fit with the lowest AIC, determined by \code{\link{AIC}}.}
#' \item{nls}{\code{nls} object for the chosen model generated by the \code{\link{nls}} function.}
#' \item{reliable}{(Logical) Indicates whether the performed fit is reliable (to be set manually).}
#' \item{fitFlag}{(Logical) Indicates whether a parametric model was successfully fitted on the data.}
#' \item{control}{Object of class \code{grofit.control} containing list of options passed to the function as \code{control}.}
#'
#' @keywords internal
#' @noRd
#'
growth.param <- function(time, data, gcID = "undefined", control)
    {
    time.in <- time
    data.in <- data

    if (!is.null(control$t0) &&
        !is.na(control$t0) &&
        control$t0 != "")
        {
        t0 <- as.numeric(control$t0)
    } else
    {
        t0 <- 0
    }

    # Implement min.growth into dataset
    if (!is.null(control$min.growth))
        {
        if (!is.na(control$min.growth) &&
            control$min.growth != 0)
            {
            min.growth <- control$min.growth
            time <- time[max(
                which.min(abs(time.in - t0)),
                which.min(abs(data - min.growth))
            ):length(time)]
            data <- data[max(
                which.min(abs(time.in - t0)),
                which.min(abs(data - min.growth))
            ):length(data)]
        }
    }

    # Perform log transformation of data
    if (control$log.y.model == TRUE)
    {
        data <- log(data/data[1])
    }

    # apply t0 to dataset
    if (is.numeric(t0) &&
        t0 > 0)
        {
        data <- data[which.min(abs(time - t0)):length(data)]
        time <- time[which.min(abs(time - t0)):length(time)]
    }

    # /// determine which values are not valid
    bad.values <- (is.na(time)) |
        (time < 0) | (is.na(data)) |
        (data < 0)

    # /// remove bad values or stop program
    if (TRUE %in% bad.values)
    {
        if (control$neg.nan.act == FALSE)
        {
            time <- time[!bad.values]
            data <- data[!bad.values]
        } else
        {
            stop("Bad values in gcFitModel")
        }
    }

    # fitting parametric growth curves
    y.model <- NULL
    bestAIC <- NULL
    bestRMSE <- NULL
    best <- NULL
    used <- NULL

    # starting values for parametric fitting from
    # spline fit
    control.tmp <- control
    control.tmp$fit.opt <- "s"
    control.tmp$log.y.spline <- control$log.y.model
    nonpara <- growth.gcFitSpline(time.in, data.in, gcID, control.tmp)
    if (nonpara$fitFlag == FALSE)
    {
        gcFitModel <- list(
            time.in = time, data.in = data, raw.time = time,
            raw.data = data, gcID = gcID, fit.time = NA,
            fit.data = NA, parameters = list(
                A = NA, mu = 0, tD = NA, lambda = NA,
                integral = NA
            ),
            model = NA, nls = NA, reliable = NULL,
            fitFlag = FALSE, control = control
        )
        class(gcFitModel) <- "gcFitModel"
        return(gcFitModel)
    }
    mu.start <- nonpara$parameters$mu
    lambda.start <- nonpara$parameters$lambda
    A.start <- nonpara$parameters$A

    # /// determine length of model names
    l <- 10
    for (i in 1:length(control$model.type))
        {
        l[i] <- nchar((control$model.type)[i])
    }
    lmax <- max(l)

    # /// loop over all parametric models
    # requested
    for (i in 1:length(control$model.type))
        {
        if (control$suppress.messages == FALSE)
        {
            cat(
                paste("--> Try to fit model", (control$model.type)[i])
            )
        }
        initmodel <- paste("init", (control$model.type)[i], sep = "")

        formulamodel <- as.formula(
            paste(
                "data ~ ", (control$model.type)[i],
                "(time, A, mu, lambda, addpar)", sep = ""
            )
        )
        if ((exists((control$model.type)[i])) &&
            (exists(initmodel)))
            {
                init.model <- do.call(
                  initmodel, list(
                    y = data, time = time, A = A.start,
                    mu = mu.start, lambda = lambda.start
                )
              )

                y.model <- try(
                  nls(formulamodel, start = init.model),
                  silent = TRUE
              )

                if (methods::is(y.model) ==
                  "nls")
                  {
                  AIC <- stats::AIC(y.model)
                  model_residuals <- stats::residuals(y.model)
                  RMSE <- sqrt(mean(model_residuals^2))
                }

                if (control$suppress.messages == FALSE)
                {
                  if (methods::is(y.model) ==
                    "nls")
                    {
                    if (y.model$convInfo$isConv ==
                      TRUE)
                      {
                      message(
                        paste(rep(".", lmax + 3 - l[i])),
                        " OK"
                    )
                    } else
                    {
                      message(
                        paste(
                          rep(".", lmax + 3 - l[i]),
                          " nls() failed to converge with stopCode ",
                          as.character(y.model$convInfo$stopCode)
                      )
                    )
                    }
                  } else
                  {
                    message(
                      paste(rep(".", lmax + 3 - l[i])),
                      " ERROR in nls(). For further information see help(growth.gcFitModel)"
                  )
                  }
                }
                if (exists("AIC", inherits = FALSE) &&
                  FALSE %in% is.null(AIC))
                    {
                  if (is.null(best) ||
                    AIC < bestAIC)
                    {
                    bestAIC <- AIC
                    bestRMSE <- RMSE
                    best <- y.model
                    used <- (control$model.type)[i]
                  }
                }
            }  # of if ( (exists((control$model.type)[i])) ...
 else
        {
          stop(
            paste0(
              "The model definition '",
              (control$model.type)[i],
              " does not exist! Spelling error?"
            )
          )
        }
        y.model <- NULL
    }

    if (control$suppress.messages == FALSE)
    {
        cat("\n")
        cat(
            paste0(
                "Best fitting model: ", sub(
                  "\\(.+", "", sub(
                    "data", "", paste(
                      formula(best),
                      collapse = ""
                  )
                )
              )
            )
        )
    }
    # /// extract parameters from data fit
    if (is.null(best) ==
        FALSE)
        {
        Abest <- summary(best)$parameters["A",
            1:2]
        mubest <- summary(best)$parameters["mu",
            1:2]

        if (any(
            grepl("addpar", rownames(summary(best)$parameters))
        ))
            {
            fitparbest <- summary(best)$parameters[grep("addpar", rownames(summary(best)$parameters)),
                1:2]
            if (summary(best)[["formula"]][[3]][[1]] ==
                "richards")
                {
                fitparbest <- list(nu = as.data.frame(t(fitparbest)))
            } else if (summary(best)[["formula"]][[3]][[1]] ==
                "gompertz.exp")
                {
                fitparbest <- list(
                  alpha = fitparbest[1, ], t_shift = fitparbest[2,
                    ]
              )
            } else if (summary(best)[["formula"]][[3]][[1]] ==
                "huang")
                {
                fitparbest <- list(y0 = as.data.frame(t(fitparbest)))
            } else if (summary(best)[["formula"]][[3]][[1]] ==
                "baranyi")
                {
                fitparbest <- list(y0 = as.data.frame(t(fitparbest)))
            }
        }
        fitFlag <- TRUE
        lambdabest <- summary(best)$parameters["lambda",
            1:2]

        best.spline <- stats::smooth.spline(
            time, fitted.values(best),
            spar = 0, keep.data = FALSE
        )
        best.deriv1 <- stats::predict(best.spline, deriv = 1)
        mumax.index <- which.max(best.deriv1$y)

        y.max <- fitted.values(best)[mumax.index]
        t.max <- time[mumax.index]
        b.tangent <- y.max - max(best.deriv1$y) *
            t.max

        if (length(time) ==
            length(as.numeric(fitted.values(best))))
                {
            Integralbest <- low.integrate(time, as.numeric(fitted.values(best)))
        } else
        {
            Integralbest <- NA
        }
    } else
    {
        if (control$suppress.messages == FALSE)
        {
            message(
                "gcFitModel: Unable to fit this curve parametrically!"
            )
        }
        Abest <- c(NA, NA)
        mubest <- c(NA, NA)
        lambdabest <- c(NA, NA)
        Integralbest <- NA
        fitFlag <- FALSE
        b.tangent <- NA
        bestRMSE <- NA
    }
    dY <- ifelse(
        is.null(best),
        NA, max(fitted.values(best)) -
            min(fitted.values(best))
    )

    gcFitModel <- list(
        time.in = time.in, data.in = data.in, raw.time = time,
        raw.data = data, gcID = gcID, fit.time = time,
        fit.data = if (is.null(best))
            {
            NA
        } else
        {
            as.numeric(fitted.values(best))
        }, parameters = list(
            A = Abest, A.orig = if (control$log.y.model ==
                TRUE)
                {
                data.in[1] * exp(Abest)
            } else
            {
                Abest
            }, dY = dY, dY.orig = if (control$log.y.model ==
                TRUE)
                {
                data.in[1] * exp(dY)
            } else
            {
                dY
            }, mu = mubest, tD = log(2)/as.numeric(mubest),
            lambda = lambdabest, b.tangent = b.tangent,
            RMSE = bestRMSE,
            fitpar = if (exists("fitparbest"))
                {
                fitparbest
            } else
            {
                NULL
            }, integral = Integralbest
        ),
        model = used, nls = best, reliable = NULL,
        fitFlag = fitFlag, control = control
    )

    class(gcFitModel) <- "gcFitModel"

    invisible(gcFitModel)
}

#' Calculate the values of the exponential growth model for given time points and growth parameters.
#'
#' @param time A vector of time values
#' @param parms A vector of length 2 with parameters: 1. start growth value (y0), 2. exponential growth rate (mumax)
#'
#' @return A dataframe with time and calculated growth columns
#' @keywords internal
#' @noRd
#'
grow_exponential <- function(time, parms)
    {
    if (is.null(names(parms)))
        {
        y0 <- parms[1]
        mumax <- parms[2]
    } else
    {
        y0 <- parms["y0_lm"]
        mumax <- ifelse(
            !is.na(parms["max_slope"]),
            parms["max_slope"], parms["mumax"]
        )
    }
    y <- y0 * exp(mumax * time)
    return(as.matrix(data.frame(time = time, y = y)))
}

#' Calculate the values of the linear growth model for given time points and growth parameters.
#'
#' @param time A vector of time values
#' @param parms A vector of length 2 with parameters: 1. start growth value (y0), 2. growth rate as slope of the linear model (mumax)
#'
#' @return A dataframe with time and calculated growth columns
#' @keywords internal
#' @noRd
#'
grow_linear <- function(time, parms)
    {
    if (length(parms) >
        2)
        {
        y0 <- parms[4]
        mumax <- parms[5]
    } else
    {
        y0 <- parms[1]
        mumax <- parms[2]
    }
    y <- y0 + mumax * time
    return(as.matrix(data.frame(time = time, y = y)))
}

#' Function to estimate the area under a curve given as x and y(x) values
#'
#' @param x Numeric vector of x values.
#' @param y Numeric vector of y values with the same length as \code{x}.
#'
#' @return Numeric value: Area under the smoothed spline.
#' @export
#'
#' @details The function uses the the R internal function \code{\link{smooth.spline}}.
#' @seealso \code{\link{smooth.spline}}
#'
#' @examples
#' # Create random growth dataset
#' rnd.dataset <- rdm.data(d = 35, mu = 0.8, A = 5, label = 'Test1')
#'
#' # Extract time and growth data for single sample
#' time <- rnd.dataset$time[1,]
#' data <- as.numeric(rnd.dataset$data[1,-(1:3)]) # Remove identifier columns
#'
#' plot(time, data)
#'
#' print(low.integrate(time, data))
#'
low.integrate <- function(x, y)
    {
    if (is.vector(x) ==
        FALSE || is.vector(y) ==
        FALSE)
        stop(
            "low.integrate: two vectors x and y are needed !"
        )
    if (length(x) !=
        length(y))
        stop(
            "low.integrate: x and y have to be of same length !"
        )
    spline <- NULL
    try(spline <- smooth.spline(x, y, keep.data = FALSE))
    if (is.null(spline))
        {
        try(
            spline <- smooth.spline(x, y, keep.data = FALSE, spar = 0.1)
        )
    }
    if (is.null(spline))
        {
        try(
            spline <- smooth.spline(x, y, keep.data = FALSE, spar = 0.1)
        )
        warning("Spline could not be fitted to data!")
        stop()
    }
    f <- function(t)
        {
        p <- stats::predict(spline, t)
        f <- p$y
    }
    low.integrate <- integrate(
        f, min(x),
        max(x)
    )$value
}

#' Calculate the values of the logistic growth model for given time points and growth parameters.
#'
#' @param time A vector of time values
#' @param A Maximum growth value
#' @param mu Growth rate
#' @param lambda Lag time
#' @param addpar Additional parameters have no effect in this type of model. They belong to the standard model description in \code{QurvE} and are initialized as \code{addpar=NULL} in the function header.
#'
#' @return A vector of growth values
#' @keywords internal
#' @noRd
#'
logistic <- function(time, A, mu, lambda, addpar = NULL)
    {
    A <- A[1]
    mu <- mu[1]
    lambda <- lambda[1]
    if (is.numeric(time) ==
        FALSE)
        stop("Need numeric vector for: time")
    if (is.numeric(mu) ==
        FALSE)
        stop("Need numeric vector for: mu")
    if (is.numeric(lambda) ==
        FALSE)
        stop("Need numeric vector for: lambda")
    if (is.numeric(A) ==
        FALSE)
        stop("Need numeric vector for: A")
    y <- A/(1 + exp(4 * mu * (lambda - time)/A + 2))
    logistic <- y
}

#' Generate initial values for parameter estimation with the logistic growth model
#'
#'\code{initlogistic} receives time points, growth data, values for \eqn{A}, \eqn{\mu} and \eqn{\lambda} and returns a list object which entries are used as initial values in the nonlinear fit procedure \code{\link{nls}} with the logistic growth model.
#'
#' @param time A vector of time values.
#' @param y A vector with growth values (e.g., growth).
#' @param A Maximum of the curve. If a vector is provided only the first entry is used.
#' @param mu Maximum slope. If a vector is provided only the first entry is used.
#' @param lambda Lag time. If a vector is provided only the first entry is used.
#'
#' @return The parameters input as arguments in a list:
#' \item{A}{Maximum of the curve.}
#' \item{mu}{Maximum slope.}
#' \item{lambda}{Lag-phase.}
#' \item{addpar}{additional parameters are not used in this type of model.}
#'
#' @keywords internal
#' @noRd
#'
initlogistic <- function(time, y, A, mu, lambda)
    {
    if (is.numeric(time) ==
        FALSE)
        stop("Need numeric vector for: time")
    if (is.numeric(y) ==
        FALSE)
        stop("Need numeric vector for: y")
    if (is.numeric(mu) ==
        FALSE)
        stop("Need numeric vector for: mu")
    if (is.numeric(lambda) ==
        FALSE)
        stop("Need numeric vector for: lambda")
    if (is.numeric(A) ==
        FALSE)
        stop("Need numeric vector for: A")
    A <- max(y)
    mu <- mu[1]
    lambda <- lambda[1]
    initlogistic <- list(A = A, mu = mu, lambda = lambda, addpar = NULL)
}

#' Generate initial values for parameter estimation with the Richard's growth model
#'
#'\code{initlogistic} receives time points, growth data, values for \eqn{A}, \eqn{\mu} and \eqn{\lambda} and returns a list object which entries are used as initial values in the nonlinear fit procedure \code{\link{nls}} with the Richard's growth model.
#'
#' @param time A vector of time values.
#' @param y A vector with growth values (e.g., growth).
#' @param A Maximum of the curve. If a vector is provided only the first entry is used.
#' @param mu Maximum slope. If a vector is provided only the first entry is used.
#' @param lambda Lag time. If a vector is provided only the first entry is used.
#'
#' @return The parameters input as arguments in a list:
#' \item{A}{Maximum of the curve.}
#' \item{mu}{Maximum slope.}
#' \item{lambda}{Lag time.}
#' \item{addpar}{Shape exponent \eqn{\nu}.}
#'
#' @keywords internal
#' @noRd
#'
initrichards <- function(time, y, A, mu, lambda)
    {
    if (is.numeric(time) ==
        FALSE)
        stop("Need numeric vector for: time")
    if (is.numeric(y) ==
        FALSE)
        stop("Need numeric vector for: y")
    if (is.numeric(mu) ==
        FALSE)
        stop("Need numeric vector for: mu")
    if (is.numeric(lambda) ==
        FALSE)
        stop("Need numeric vector for: lambda")
    if (is.numeric(A) ==
        FALSE)
        stop("Need numeric vector for: A")
    nu <- 0.1
    A <- max(y)
    mu <- mu[1]
    lambda <- lambda[1]
    initrichards <- list(A = A, mu = mu, lambda = lambda, addpar = nu)
}

#' Generate initial values for parameter estimation with the Huang's growth model
#'
#'\code{initlogistic} receives time points, growth data, values for \eqn{A}, \eqn{\mu} and \eqn{\lambda} and returns a list object which entries are used as initial values in the nonlinear fit procedure \code{\link{nls}} with the Huang's growth model.
#'
#' @param time A vector of time values.
#' @param y A vector with growth values (e.g., growth).
#' @param A Maximum of the curve. If a vector is provided only the first entry is used.
#' @param mu Maximum slope. If a vector is provided only the first entry is used.
#' @param lambda Lag time. If a vector is provided only the first entry is used.
#'
#' @return The parameters input as arguments in a list:
#' \item{A}{Maximum of the curve.}
#' \item{mu}{Maximum slope.}
#' \item{lambda}{Lag time.}
#' \item{addpar}{Minimum of the curve \eqn{y0}.}
#'
#' @keywords internal
#' @noRd
#'
inithuang <- function(time, y, A, mu, lambda)
    {
    if (is.numeric(time) ==
        FALSE)
        stop("Need numeric vector for: time")
    if (is.numeric(y) ==
        FALSE)
        stop("Need numeric vector for: y")
    if (is.numeric(mu) ==
        FALSE)
        stop("Need numeric vector for: mu")
    if (is.numeric(lambda) ==
        FALSE)
        stop("Need numeric vector for: lambda")
    if (is.numeric(A) ==
        FALSE)
        stop("Need numeric vector for: A")
    y0 <- y[1]
    A <- max(y)
    mu <- mu[1]
    lambda <- lambda[1]
    inithuang <- list(A = A, mu = mu, lambda = lambda, addpar = y0)
}

#' Calculate the values of the Huang growth model for given time points and growth parameters.
#'
#' @param time A vector of time values
#' @param A Maximum growth value
#' @param mu Growth rate
#' @param lambda Lag time
#' @param addpar additional parameters (here: minimum growth)
#'
#' @return A vector of growth values
#' @keywords internal
#' @noRd
#'
huang <- function(time, A, mu, lambda, addpar)
    {
    A <- A[1]
    mu <- mu[1]
    lambda <- lambda[1]
    y0 <- addpar[1]
    if (is.numeric(time) ==
        FALSE)
        stop("Need numeric vector for: time")
    if (is.numeric(mu) ==
        FALSE)
        stop("Need numeric vector for: mu")
    if (is.numeric(lambda) ==
        FALSE)
        stop("Need numeric vector for: lambda")
    if (is.numeric(A) ==
        FALSE)
        stop("Need numeric vector for: A")
    suppressWarnings(
        y <- y0 + A - log(
            exp(y0) +
                (exp(A) -
                  exp(y0)) *
                  exp(
                    -mu * (time + 0.25 * log(
                      (1 + exp(-4 * (time - lambda)))/(1 +
                        exp(4 * lambda))
                  ))
                )
        )
    )
    huang <- y
}

#' Calculate the values of the Baranyi growth model for given time points and growth parameters.
#'
#' @param time A vector of time values
#' @param A Maximum growth value
#' @param mu Growth rate
#' @param lambda Lag time
#' @param addpar additional parameters (here: minimum growth)
#'
#' @return A vector of growth values
#' @keywords internal
#' @noRd
#'
baranyi <- function(time, A, mu, lambda, addpar)
    {
    A <- A[1]
    mu <- mu[1]
    lambda <- lambda[1]
    y0 <- addpar[1]
    if (is.numeric(time) ==
        FALSE)
        stop("Need numeric vector for: time")
    if (is.numeric(mu) ==
        FALSE)
        stop("Need numeric vector for: mu")
    if (is.numeric(lambda) ==
        FALSE)
        stop("Need numeric vector for: lambda")
    if (is.numeric(A) ==
        FALSE)
        stop("Need numeric vector for: A")
    suppressWarnings(
        B <- time + 1/mu * log(
            exp(-mu * time) +
                exp(-mu * lambda) -
                exp(-mu * (time + lambda))
        )
    )
    suppressWarnings(
        y <- y0 + mu * B - log(
            1 + (exp(mu * B) -
                1)/exp(A - y0)
        )
    )
    baranyi <- y
}

#' Generate initial values for parameter estimation with the Baranyi's growth model
#'
#'\code{initlogistic} receives time points, growth data, values for \eqn{A}, \eqn{\mu} and \eqn{\lambda} and returns a list object which entries are used as initial values in the nonlinear fit procedure \code{\link{nls}} with the Baranyi's growth model.
#'
#' @param time A vector of time values.
#' @param y A vector with growth values (e.g., growth).
#' @param A Maximum of the curve. If a vector is provided only the first entry is used.
#' @param mu Maximum slope. If a vector is provided only the first entry is used.
#' @param lambda Lag time. If a vector is provided only the first entry is used.
#'
#' @return The parameters input as arguments in a list:
#' \item{A}{Maximum of the curve.}
#' \item{mu}{Maximum slope.}
#' \item{lambda}{Lag time.}
#' \item{addpar}{Minimum of the curve \eqn{y0}.}
#'
#' @keywords internal
#' @noRd
#'
initbaranyi <- function(time, y, A, mu, lambda)
    {
    if (is.numeric(time) ==
        FALSE)
        stop("Need numeric vector for: time")
    if (is.numeric(y) ==
        FALSE)
        stop("Need numeric vector for: y")
    if (is.numeric(mu) ==
        FALSE)
        stop("Need numeric vector for: mu")
    if (is.numeric(lambda) ==
        FALSE)
        stop("Need numeric vector for: lambda")
    if (is.numeric(A) ==
        FALSE)
        stop("Need numeric vector for: A")
    y0 <- y[1]
    A <- max(y)
    mu <- mu[1]
    lambda <- lambda[1]
    initbaranyi <- list(A = A, mu = mu, lambda = lambda, addpar = y0)
}


# liquori <- function (time, A, mu, addpar) { A
# <- A[1] mu <- mu[1] y0 <- addpar[1] t_a1 <-
# addpar[2] t_a2 <- addpar[3] t_b1 <- addpar[4]
# t_b2 <- addpar[5] x <- addpar[6] if
# (is.numeric(time) == FALSE) stop('Need numeric
# vector for: time') if (is.numeric(mu) == FALSE)
# stop('Need numeric vector for: mu') if
# (is.numeric(y0) == FALSE) stop('Need numeric
# vector for: y0') if (is.numeric(A) == FALSE)
# stop('Need numeric vector for: A') if
# (is.numeric(t_a1) == FALSE) stop('Need numeric
# vector for: addpar[2]') if (is.numeric(t_a2) ==
# FALSE) stop('Need numeric vector for:
# addpar[3]') if (is.numeric(t_b1) == FALSE)
# stop('Need numeric vector for: addpar[4]') if
# (is.numeric(t_b2) == FALSE) stop('Need numeric
# vector for: addpar[5]') t <- time y1 <-
# (1-exp(-(t/t_a1)))/(1-exp(-(t/t_a1))+exp(-(t/t_a2)))
# y2 <-
# (1-exp(-(t/t_b1)))/(1-exp(-(t/t_b1))+exp(-(t/t_b2)))
# y <- y0 + A*x*y1 + A*(1-x)*y2 liquori <- y }

#' Calculate the values of the Richards growth model for given time points and growth parameters.
#'
#' @param time A vector of time values
#' @param A Maximum growth value
#' @param mu Growth rate
#' @param lambda Lag time
#' @param addpar additional parameters (here: Shape exponent \eqn{\nu}.)
#'
#' @return A vector of growth values
#' @keywords internal
#' @noRd
#'
richards <- function(time, A, mu, lambda, addpar)
    {
    A <- A[1]
    mu <- mu[1]
    lambda <- lambda[1]
    nu <- addpar[1]
    if (is.numeric(time) ==
        FALSE)
        stop("Need numeric vector for: time")
    if (is.numeric(mu) ==
        FALSE)
        stop("Need numeric vector for: mu")
    if (is.numeric(lambda) ==
        FALSE)
        stop("Need numeric vector for: lambda")
    if (is.numeric(A) ==
        FALSE)
        stop("Need numeric vector for: A")
    if (is.numeric(nu) ==
        FALSE)
        stop("Need numeric vector for: addpar[1]")
    y <- A * (1 + nu * exp(1 + nu) *
        exp(mu * (1 + nu)^(1 + 1/nu) * (lambda - time)/A))^(-1/nu)
    richards <- y
}

#' Calculate the values of the Gompertz growth model for given time points and growth parameters.
#'
#' @param time A vector of time values
#' @param A Maximum growth value
#' @param mu Growth rate
#' @param lambda Lag time
#' @param addpar Additional parameters have no effect in this type of model. They belong to the standard model description in \code{grofit} and are initialized as \code{addpar=NULL} in the function header.
#'
#' @return A vector of growth values
#' @keywords internal
#' @noRd
#'
gompertz <- function(time, A, mu, lambda, addpar = NULL)
    {
    A <- A[1]
    mu <- mu[1]
    lambda <- lambda[1]
    if (is.numeric(time) ==
        FALSE)
        stop("Need numeric vector for: time")
    if (is.numeric(mu) ==
        FALSE)
        stop("Need numeric vector for: mu")
    if (is.numeric(lambda) ==
        FALSE)
        stop("Need numeric vector for: lambda")
    if (is.numeric(A) ==
        FALSE)
        stop("Need numeric vector for: A")
    e <- exp(1)
    y <- A * exp(-exp(mu * e * (lambda - time)/A + 1))
    gompertz <- y
}

#' Generate initial values for parameter estimation with the Gompertz's growth model
#'
#'\code{initlogistic} receives time points, growth data, values for \eqn{A}, \eqn{\mu} and \eqn{\lambda} and returns a list object which entries are used as initial values in the nonlinear fit procedure \code{\link{nls}} with the Gompertz's growth model.
#'
#' @param time A vector of time values.
#' @param y A vector with growth values (e.g., growth).
#' @param A Maximum of the curve. If a vector is provided only the first entry is used.
#' @param mu Maximum slope. If a vector is provided only the first entry is used.
#' @param lambda Lag time. If a vector is provided only the first entry is used.
#'
#' @return The parameters input as arguments in a list:
#' \item{A}{Maximum of the curve.}
#' \item{mu}{Maximum slope.}
#' \item{lambda}{Lag time.}
#' \item{addpar}{Additional parameters are not used in this type of model.}
#'
#' @keywords internal
#' @noRd
#'
initgompertz <- function(time, y, A, mu, lambda)
    {
    if (is.numeric(time) ==
        FALSE)
        stop("Need numeric vector for: time")
    if (is.numeric(y) ==
        FALSE)
        stop("Need numeric vector for: y")
    if (is.numeric(mu) ==
        FALSE)
        stop("Need numeric vector for: mu")
    if (is.numeric(lambda) ==
        FALSE)
        stop("Need numeric vector for: lambda")
    if (is.numeric(A) ==
        FALSE)
        stop("Need numeric vector for: A")
    A <- max(y)
    mu <- mu[1]
    lambda <- lambda[1]
    initgompertz <- list(A = A, mu = mu, lambda = lambda, addpar = NULL)
}

#' Generate initial values for parameter estimation with the modified Gompertz growth model
#'
#'\code{initlogistic} receives time points, growth data, values for \eqn{A}, \eqn{\mu} and \eqn{\lambda} and returns a list object which entries are used as initial values in the nonlinear fit procedure \code{\link{nls}} with the Gompertz's growth model.
#'
#' @param time A vector of time values.
#' @param y A vector with growth values (e.g., growth).
#' @param A Maximum of the curve. If a vector is provided only the first entry is used.
#' @param mu Maximum slope. If a vector is provided only the first entry is used.
#' @param lambda Lag time. If a vector is provided only the first entry is used.
#'
#' @return The parameters input as arguments in a list:
#' \item{A}{Maximum of the curve.}
#' \item{mu}{Maximum slope.}
#' \item{lambda}{Lag time.}
#' \item{addpar}{Two element vector defining scaling parameter \eqn{\alpha} and shifting parameter \eqn{t_{shift}}.}
#'
#' @keywords internal
#' @noRd
#'
initgompertz.exp <- function(time, y, A, mu, lambda)
    {
    if (is.numeric(time) ==
        FALSE)
        stop("Need numeric vector for: time")
    if (is.numeric(y) ==
        FALSE)
        stop("Need numeric vector for: y")
    if (is.numeric(mu) ==
        FALSE)
        stop("Need numeric vector for: mu")
    if (is.numeric(lambda) ==
        FALSE)
        stop("Need numeric vector for: lambda")
    if (is.numeric(A) ==
        FALSE)
        stop("Need numeric vector for: A")
    alpha <- 0.1
    t_shift <- max(time)/10
    A <- max(y)
    mu <- mu[1]
    lambda <- lambda[1]
    initgompertz.exp <- list(
        A = A, mu = mu, lambda = lambda, addpar = c(alpha, t_shift)
    )
}

#' Calculate the values of the modified Gompertz growth model for given time points and growth parameters.
#'
#' @param time A vector of time values
#' @param A Maximum growth value
#' @param mu Growth rate
#' @param lambda Lag time
#' @param addpar Numeric vector of size two, addpar`[`1`]` corresponds to scaling parameter \eqn{\alpha} and addpar`[`2`]` corresponds to shifting parameter \eqn{t_{shift}}.
#'
#' @return A vector of growth values
#' @keywords internal
#' @noRd
#'
gompertz.exp <- function(time, A, mu, lambda, addpar)
    {
    A <- A[1]
    mu <- mu[1]
    lambda <- lambda[1]
    alpha <- addpar[1]
    t_shift <- addpar[2]
    if (is.numeric(time) ==
        FALSE)
        stop("Need numeric vector for: time")
    if (is.numeric(mu) ==
        FALSE)
        stop("Need numeric vector for: mu")
    if (is.numeric(lambda) ==
        FALSE)
        stop("Need numeric vector for: lambda")
    if (is.numeric(A) ==
        FALSE)
        stop("Need numeric vector for: A")
    if (is.numeric(alpha) ==
        FALSE)
        stop("Need numeric vector for: addpar[1]")
    if (is.numeric(t_shift) ==
        FALSE)
        stop("Need numeric vector for: addpar[2]")
    e <- exp(1)
    y <- A * exp(-exp(mu * e * (lambda - time)/A + 1)) +
        A * exp(alpha * (time - t_shift))
    gompertz.exp <- y
}