R/fitGrowth.R

In pcvr: Plant Phenotyping and Bayesian Statistics

#' Ease of use wrapper function for fitting various growth models specified by \link{growthSS}
#'
#' @param ss A list generated by \code{growthSS}.
#' @param ... Additional arguments passed to model fitting functions determined by \code{ss$type}.
#' For type = "nlme" these are passed to \code{nlme::nlmeControl}, not \code{nlme::nlme}.
#' Additional arguments are documented in \link{fitGrowthbrms}, \link{fitGrowthnlme},
#' \link{fitGrowthnls}, \link{fitGrowthnlrq}, \link{fitGrowthmgcvgam}, \link{fitGrowthsurvreg},
#' \link{fitGrowthflexsurv}.
#' @keywords Bayesian brms nlme nls nlrq
#' @return A fit model from the selected type.
#' @seealso \link{growthPlot} for model visualization, \link{testGrowth} for hypothesis testing,
#' \link{barg} for Bayesian model reporting metrics.
#' @examples
#'
#' simdf <- growthSim("logistic",
#'   n = 20, t = 25,
#'   params = list("A" = c(200, 160), "B" = c(13, 11), "C" = c(3, 3.5))
#' )
#' ss <- growthSS(
#'   model = "logistic", form = y ~ time | group,
#'   df = simdf, type = "nls"
#' )
#' fitGrowth(ss)
#' ss <- growthSS(
#'   model = "gam", form = y ~ time | group,
#'   df = simdf, type = "nls"
#' )
#' fitGrowth(ss)
#'
#' @export

fitGrowth <- function(ss, ...) {
  if (ss$model == "gam") {
    suffix <- "gam"
  } else {
    suffix <- NULL
  }
  fit_function <- match.fun(paste0("fitGrowth", ss$type, suffix))
  fit <- fit_function(ss, ...)
  return(fit)
}

#' Ease of use brms wrapper function for fitting growth models specified by \code{growthSS}
#'
#' Helper function generally called from \link{fitGrowth}.
#'
#' @param ss A list generated by \code{growthSS}.
#' @param iter A number of iterations to sample for each chain.
#' By default half this length is taken as warm-up for the MCMC algorithm.
#' This defaults to 2000.
#' @param cores A number of cores to run in parallel.
#' This defaults to 1 if the "mc.cores" option is not set.
#' Generally this is specified as one core per chain so that the model is fit in parallel.
#' @param chains A number of markov chains to use, this defaults to 4.
#' @param prior A \code{brmsprior} object if \code{growthSS} did not have priors specified.
#' If left NULL (the default) and ss does not contain priors then a warning is
#' issued but the model will still attempt to fit.
#' @param backend A backend for brms to use Stan through.
#'   This defaults to use "cmdstanr".
#' @param silent Passed to \code{brms::brm} to control verbosity.
#' This defaults to 0, the most verbose option so that messages and progress are printed.
#' With changes to \code{cmdstanr} and \code{brms} this may be removed, but the option
#' will be available through \code{...}. Note that this is likely to print lots of
#' messages during warmup iterations as the MCMC gets started.
#' @param ... Additional arguments passed to \code{brms::brm}.
#' @keywords Bayesian brms
#' @return A \code{brmsfit} object, see \code{?`brmsfit-class`} for details.
#' @export

fitGrowthbrms <- function(ss, iter = 2000, cores = getOption("mc.cores", 1), chains = 4, prior = NULL,
                          backend = "cmdstanr", silent = 0, ...) {
  if (!"prior" %in% names(ss) && is.null(prior)) {
    warning(
      paste0(
        "No prior was specified. Flat priors will be used, this is likely to cause problems in model ",
        "fitting and yield less accurate results. If you are fitting a gam then ignore this warning."
      )
    )
  } else if ("prior" %in% names(ss)) {
    prior <- ss$prior
  }
  fit <- brms::brm(
    formula = ss$formula, prior = prior, data = ss$df, family = ss$family,
    iter = iter, init = ss$initfun, cores = cores, chains = chains,
    backend = backend, silent = silent, ...
  )
  return(fit)
}

#' @rdname fitGrowthbrms
#' @export
fitGrowthbrmsgam <- fitGrowthbrms


#' Ease of use nlme wrapper function for fitting growth models specified by \code{growthSS}
#'
#' Helper function generally called from \link{fitGrowth}.
#'
#' @param ss A list generated by \code{growthSS}.
#' @param ... Additional arguments passed to \code{nlme::nlmeControl}.
#' @keywords nlme
#' @return An \code{nlme} object, see \code{?nlme} for details.
#' @importFrom nlme nlme
#' @export

fitGrowthnlme <- function(ss, ...) {
  fit <- do.call(nlme,
    args = list(
      model = ss[["formula"]][["model"]],
      data = quote(ss[["df"]]),
      fixed = ss[["formula"]][["fixed"]],
      random = ss[["formula"]][["random"]],
      groups = ss[["formula"]][["groups"]],
      weights = ss[["formula"]][["weights"]],
      correlation = ss[["formula"]][["cor_form"]],
      start = ss[["start"]],
      control = nlme::nlmeControl(returnObject = TRUE, ...)
    )
  )
  return(fit)
}

#' Ease of use lme wrapper function for fitting gams specified by \code{growthSS}
#'
#' Helper function generally called from \link{fitGrowth}.
#'
#' @param ss A list generated by \code{growthSS}.
#' @param ... Additional arguments passed to \code{nlme::lmeControl}.
#' @keywords nlme gam
#' @return An \code{lme} object, see \code{?lme} for details.
#' @importFrom nlme lme
#' @export

fitGrowthnlmegam <- function(ss, ...) {
  fit <- do.call(lme,
    args = list(
      fixed = ss[["formula"]][["model"]],
      data = quote(ss[["df"]]),
      random = ss[["formula"]][["random"]],
      weights = ss[["formula"]][["weights"]],
      correlation = ss[["formula"]][["cor_form"]],
      control = nlme::lmeControl(returnObject = TRUE, ...)
    )
  )
  return(fit)
}


#' Ease of use nls wrapper function for fitting growth models specified by \code{growthSS}
#'
#' Helper function generally called from \link{fitGrowth}.
#'
#' @param ss A list generated by \code{growthSS}.
#' @param ... Additional arguments passed to \code{stats::nls}.
#' @keywords nls
#' @return An \code{nls} object, see \code{?nls} for details.
#' @importFrom stats nls
#' @export

fitGrowthnls <- function(ss, ...) {
  fit <- stats::nls(
    formula = ss[["formula"]],
    data = ss[["df"]],
    start = ss[["start"]], ...
  )
  return(fit)
}

#' Ease of use lm wrapper function for fitting gams specified by \code{growthSS}
#'
#' Helper function generally called from \link{fitGrowth}.
#'
#' @param ss A list generated by \code{growthSS}.
#' @param ... Additional arguments passed to \code{stats::lm}.
#' @keywords nls
#' @return An \code{lm} object, see \code{?lm} for details.
#' @importFrom stats lm
#' @export

fitGrowthnlsgam <- function(ss, ...) {
  fit <- stats::lm(
    formula = ss[["formula"]],
    data = ss[["df"]],
    ...
  )
  return(fit)
}

#' Ease of use lm wrapper function for fitting growth models specified by \code{mvSS}
#'
#' Helper function generally called from \link{fitGrowth}.
#'
#' @param ss A list generated by \code{mvSS}.
#' @param ... Additional arguments passed to \code{stats::lm}.
#' @keywords nls
#' @return An \code{lm} object, see \code{?lm} for details.
#' @importFrom stats lm
#' @export

fitGrowthlm <- function(ss, ...) {
  fit <- do.call("lm", args = list(
    formula = ss[["formula"]],
    data = quote(ss[["df"]]),
    weights = ss[["weights"]], ...
  ))
  ocall <- as.character(fit$call)
  ocall[4] <- "weights"
  fit$call <- as.call(str2expression(ocall))
  return(fit)
}

#' Ease of use nlrq wrapper function for fitting growth models specified by \code{growthSS}
#'
#' Helper function generally called from \link{fitGrowth}.
#'
#' @param ss A list generated by \code{growthSS}.
#' @param cores Optionally specify how many cores to run in parallel if ss$taus is >1L.
#' Defaults to 1 if mc.cores option is not set globally.
#' @param ... Additional arguments passed to \code{quantreg::nlrq}.
#' @keywords nls
#' @return An \code{nlrqModel} object if tau is length of 1 or a list of such models for multiple taus,
#'   see \code{?quantreg::nlrq} or \code{?nls::nlsModel} for details.
#' @importFrom parallel mclapply
#' @importFrom quantreg nlrq
#' @export

fitGrowthnlrq <- function(ss, cores = getOption("mc.cores", 1), ...) {
  if (length(ss[["taus"]]) > 1) {
    fits <- parallel::mclapply(ss[["taus"]], function(tau) {
      fit <- do.call("nlrq", args = list(
        formula = ss[["formula"]],
        data = quote(ss[["df"]]),
        tau = tau,
        start = ss[["start"]], ...
      ))
      return(fit)
    }, mc.cores = cores)
    names(fits) <- ss[["taus"]]
  } else {
    fits <- do.call("nlrq", args = list(
      formula = ss[["formula"]],
      data = quote(ss[["df"]]),
      tau = ss[["taus"]],
      start = ss[["start"]], ...
    ))
  }
  return(fits)
}

#' Ease of use rq wrapper function for fitting gams specified by \code{growthSS}
#'
#' Helper function generally called from \link{fitGrowth}.
#'
#' @param ss A list generated by \code{growthSS}.
#' @param cores number of cores to run in parallel
#' @param ... Additional arguments passed to \code{quantreg::rq}.
#' @keywords nls
#' @return An \code{rq} object, see \code{?rq} for details.
#' @importFrom quantreg rq
#' @export

fitGrowthnlrqgam <- function(ss, cores = getOption("mc.cores", 1), ...) {
  if (length(ss[["taus"]]) > 1) {
    fits <- parallel::mclapply(ss[["taus"]], function(tau) {
      fit <- do.call("rq", args = list(
        formula = ss[["formula"]],
        data = quote(ss[["df"]]),
        tau = tau,
        ...
      ))
      return(fit)
    }, mc.cores = cores)
    names(fits) <- ss[["taus"]]
  } else {
    fits <- do.call("rq", args = list(
      formula = ss[["formula"]],
      data = quote(ss[["df"]]),
      tau = ss[["taus"]],
      ...
    ))
  }
  return(fits)
}

#' Ease of use rq wrapper function for fitting models specified by \code{mvSS}
#'
#' Helper function generally called from \link{fitGrowth}.
#'
#' @param ss A list generated by \code{mvSS}.
#' @param cores number of cores to run in parallel
#' @param ... Additional arguments passed to \code{quantreg::rq}.
#' @keywords nls
#' @return An \code{rq} object, see \code{?rq} for details.
#' @importFrom quantreg rq
#' @export

fitGrowthrq <- function(ss, cores = getOption("mc.cores", 1), ...) {
  if (length(ss[["taus"]]) > 1) {
    fits <- parallel::mclapply(ss[["taus"]], function(tau) {
      fit <- do.call("rq", args = list(
        formula = ss[["formula"]],
        data = quote(ss[["df"]]),
        weights = ss[["weights"]],
        tau = tau,
        ...
      ))
      ocall <- as.character(fit$call)
      ocall[5] <- "weights"
      fit$call <- as.call(str2expression(ocall))
      return(fit)
    }, mc.cores = cores)
    names(fits) <- ss[["taus"]]
  } else {
    fits <- do.call("rq", args = list(
      formula = ss[["formula"]],
      data = quote(ss[["df"]]),
      weights = ss[["weights"]],
      tau = ss[["taus"]],
      ...
    ))
    ocall <- as.character(fits$call)
    ocall[5] <- "weights"
    fits$call <- as.call(str2expression(ocall))
  }
  return(fits)
}

#' Ease of use mgcv wrapper function for fitting gams specified by \code{growthSS}
#'
#' Helper function generally called from \link{fitGrowth}.
#'
#' @param ss A list generated by \code{growthSS}.
#' @param ... Additional arguments passed to \code{mgcv::gam}.
#' @keywords mgcv gam
#' @return An \code{gam} object, see \code{?gam} for details.
#' @importFrom mgcv gam s
#' @export

fitGrowthmgcvgam <- function(ss, ...) {
  fit <- do.call("gam", args = list(
    formula = ss[["formula"]],
    data = ss[["df"]],
    ...
  ))
  return(fit)
}


#' Ease of use wrapper function for fitting growth models specified by \code{growthSS}
#'
#' Helper function generally called from \link{fitGrowth}.
#'
#' @param ss A list generated by \code{growthSS}.
#' @param ... Additional arguments passed to \code{survival::survreg}.
#' @keywords nlme
#' @return A \code{survreg} object.
#' @importFrom survival survreg Surv
#' @export

fitGrowthsurvreg <- function(ss, ...) {
  fit <- do.call("survreg", args = list(
    formula = ss[["formula"]],
    data = quote(ss[["df"]]),
    dist = ss[["distribution"]],
    ...
  ))
  return(fit)
}

#' Ease of use wrapper function for fitting growth models specified by \code{growthSS}
#'
#' Helper function generally called from \link{fitGrowth}.
#'
#' @param ss A list generated by \code{growthSS}.
#' @param ... Additional arguments passed to \code{flexsurv::flexsurvreg}.
#' @keywords flexsurv
#' @return A \code{survreg} object.
#' @export

fitGrowthflexsurv <- function(ss, ...) {
  fit <- do.call(eval(parse(text = "flexsurv::flexsurvreg")), args = list(
    formula = ss[["formula"]][["f1"]],
    anc = ss[["formula"]][["f2"]],
    data = quote(ss[["df"]]),
    dist = ss[["distribution"]],
    ...
  ))
  return(fit)
}