R/fitBeeGUTS.R

In BeeGUTS: General Unified Threshold Model of Survival for Bees using Bayesian Inference

#' Fit a GUTS model for bees survival analysis using Bayesian Inference (stan)
#'
#' @description The function \code{fitBeeGUTS} estimates the parameters of a GUTS model
#' for the stochastic death (SD) or individual tolerance (IT) death mechanisms for
#' survival analysis using Bayesian inference.
#'
#' @param data An object of class \code{beeSurvData}
#' @param modelType A model type between \code{"SD"} for Stochastic Death and
#' \code{"IT"} for Individual Tolerance.
#' @param distribution A distribution for the IT death mechanism. To be chosen between
#' \code{"loglogistic"} and \code{"lognormal"}. Default is \code{"loglogistic"}
#' @param priorsList A list containing the prior distribution for the parameter considered.
#' By default, when no priors are provided (default is \code{NULL}), priors are set automatically
#' based on the experimental design (adapted from Delignette-Muller et al 2017)
#' @param parallel Logical indicating whether parallel computing should be used or not.  Default is \code{TRUE}
#' @param nCores A positive integer specifying the number of cores to use.
#' Default is one core less than maximum number of cores available
#' @param nChains A positive integer specifying the number of MCMC chains to run. Default is 3.
#' @param nIter A positive integer specifying the number of iteration to monitor for each MCMC chain. Default is 2000
#' @param nWarmup A positive integer specifying the number of warmup iteration per chain. Default is half the number of iteration
#' @param thin A positive integer specifying the interval between the iterations to monitor. Default is 1 (all iterations are monitored)
#' @param adaptDelta A double, bounded between 0 and 1 and controlling part of the sampling algorithms.
#' See the \code{control} in the function \code{stan} [rstan::stan()] of the package \code{rstan}. The default is 0.95.
#' @param odeIntegrator A string specifying the integrator used to solve the system of
#' differential equations (ODE) in the \code{stan} module. To be chosen between
#' \code{"rk45"} and \code{"bdf"}. Default is \code{"rk45"}.
#' @param relTol A double, bounded between 0 and 1 and controlling the relative tolerance
#' of the accuracy of the solutions generated by the integrator. A smaller tolerance produces
#' more accurate solution at the expanse of the computing time. Default is 1e-8
#' @param absTol A double, bounded between 0 and 1 and controlling the absolute tolerance
#' of the accuracy of the solutions generated by the integrator. A smaller tolerance produces
#' more accurate solution at the expanse of the computing time. Default is 1e-8
#' @param maxSteps A double controlling the maximum number of steps that can be
#' taken before stopping a runaway simulation. Default is 1000
#' @param ... Additional parameters to be passed to \code{sampling} from \code{stan}
#'
#' @details The automated prior determination is modified from Delignette-Muller et al.
#' by considering that the minimal concentration for the prior can be close to 0 (1e-6)
#' whereas the original paper considered the lowest non-zero concentration.
#' Similarly, the minimal kd considered for the prior calculation was reduced to allow
#' more chance to capture slow kinetics.
#'
#' @return The function \code{fitBeeGUTS} returns the parameter estimates
#' of the General Unified Threshold model of Survival (GUTS) in an object
#' of class \code{beeSurvFit}. This object is a list composed of the following:
#' \item{stanFit}{An object of S4 class \code{stanfit}. More information is available
#' in the package \code{rstan}. }
#' \item{data}{The data object provided as argument of the function}
#' \item{dataFit}{A list of data passed to the Stan model object}
#' \item{setupMCMC}{A list containing the setup used for the MCMC chains}
#' \item{modelType}{A character vector specifying the type of GUTS model used between
#' \code{SD} and \code{IT}}
#' \item{distribution}{A character vector specifying the type of distribution used in case \code{IT} was used;
#' \code{NA} otherwise}
#' \item{messages}{A character vector containing warning messages}
#'
#' @references
#' Delignette-Muller, M.L., Ruiz P. and Veber P. (2017).
#' Robust fit of toxicokinetic-toxicodynamic models using prior knowledge contained in the design of survival toxicity tests.
#' \doi{10.1021/acs.est.6b05326}
#'
#' @export
#'
#' @examples
#' \donttest{
#' data(betacyfluthrinChronic)
#' fit <- fitBeeGUTS(betacyfluthrinChronic, modelType = "SD", nIter = 1000, nCores = 2)
#' }
fitBeeGUTS <- function(data, # CHECK CORRECT DATA OBJECT IS USED !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                       modelType = NULL,
                       distribution = "loglogistic",
                       priorsList = NULL,
                       parallel = TRUE,
                       nCores = parallel::detectCores()-1L,
                       nChains = 3,
                       nIter = 2000,
                       nWarmup = floor(nIter / 2),
                       thin = 1,
                       adaptDelta = 0.95,
                       odeIntegrator = "rk45",
                       relTol = 1e-8,
                       absTol = 1e-8,
                       maxSteps = 1000,
                       ...) {
  # Check correct user inputs
  if (!is(data,"beeSurvData")) {
    stop("fitBeeGUTS: a 'data' object of class 'beeSurvData' is expected")
  }
  if (is.null(modelType) || !(modelType %in% c("SD", "IT"))) {
    stop("You need to specifiy a correct 'modelType' amongst 'SD' and 'IT'.
         'PROPER' is not yet implemented. When selecting 'IT' please also
         provide a 'distribution' parameter amongst 'loglogistic' and 'lognormal'.")
  }
  if (!(distribution %in% c("loglogistic", "lognormal"))) {
    stop("You need to specifiy a correct 'distribution' amongst 'loglogistic' and 'lognormal'.")
  }
  if (!(odeIntegrator %in% c("rk45"))) {
    stop("You need to specifiy a correct 'odeIntegrator' amongst 'rk45'.
         'bdf' is not yet implemented")
  }

  # Regroup control for the ode solver
  odeControl <- list(relTol = relTol, absTol = absTol, maxSteps = maxSteps)

  # Prepare data for inference with stan
  lsFullData <- dataFitStan(data, modelType, odeControl, priorsList)
  lsStanData <- lsFullData
  lsStanData$replicateConc <- NULL # NECESSARY?????????????????????????????????????????????????????????
  lsStanData$replicateNsurv <- NULL # NECESSARY?????????????????????????????????????????????????????????
  lsStanData$Ninit <- NULL # NECESSARY?????????????????????????????????????????????????????????

  if (modelType == "SD") {
    modelObject <- stanmodels$GUTS_SD
  }
  if (modelType == "IT") {
    modelObject <- stanmodels$GUTS_IT
    lsStanData$distribution <- switch(distribution, loglogistic = 1, lognormal = 2)
  }

  # Set options for parallel computing
  if (parallel == TRUE) {
    op <- options()
    options(mc.cores = as.integer(nCores))
    on.exit(options(op))

  }

  # Sample MCMC chains
  fit <- rstan::sampling(
    object = modelObject,
    data = lsStanData,
    chains = nChains,
    iter = nIter,
    warmup = nWarmup,
    thin = thin,
    control = list(adapt_delta = adaptDelta),
    ...)

  # cleanup parallel computing options
  if (parallel == TRUE) {
    options(op)
  }

  # Infos on MCMC chains
  setupMCMC <- data.frame(nIter = nIter,
                          nChains = nChains,
                          thinInterval = thin,
                          nWarmup = nWarmup)

  ## Warnings on fit quality

  outRhat <- rstan::summary(fit)$summary[, "Rhat"]

  if (!all(outRhat < 1.1, na.rm = TRUE)){
    msg <- "
    *** Markov chains did not converge! Do not analyze results! ***.
    Plot MCMC chains and try the following options:
    (1) if one or more chain are a simple stable line, increase 'adapt_delta' (default is 0.95).
    (2) if the variabbility between chain is great, you can increase the number of iteration (default is 2000 iteration).
    (3) if 'Conditional_Psurv_hat' is greater than 1, the ODE integration is wrong. So you can reduce the tolerance of the ODE integrator."
    warning(msg, call. = FALSE)
    print(outRhat)
  } else {
    msg <- "NA"
  }

  # Return
  lsOut <- list(stanFit = fit,
                data = data,
                dataFit = lsFullData,
                setupMCMC = setupMCMC,
                modelType = modelType,
                distribution = ifelse(modelType == "IT", distribution, "NA"),
                messages = msg)
  class(lsOut) <- "beeSurvFit"

  return(lsOut)
}