R/survFit.survDataVarExp.R

In morse: Modelling Reproduction and Survival Data in Ecotoxicology

#' @rdname survFit
#'
#' @return The function returns an object of class \code{survFitVarExp}, which is
#' a list with the following information:
#' \item{estim.par}{a table of the estimated parameters as medians and 95\%
#' credible intervals}
#' \item{mcmc}{an object of class \code{mcmc.list} with the posterior
#' distribution}
#' \item{model}{a JAGS model object}
#' \item{dic}{return the Deviance Information Criterion (DIC) if \code{dic.compute} is \code{TRUE}}
#' \item{warnings}{a table with warning messages}
#' \item{parameters}{a list of parameter names used in the model}
#' \item{n.chains}{an integer value corresponding to the number of chains used
#' for the MCMC computation}
#' \item{mcmcInfo}{a table with the number of iterations, chains, adaptation, warmup and the thinning interval} 
#' \item{jags.data}{a list of the data passed to the JAGS model}
#' \item{model_type}{the type of TKTD model used: \code{SD} or \code{IT}}
#' 
#' 
#' @examples
#'
#' # Example with time-variable exposure profile
#' # (1) Load the survival data
#' data("propiconazole_pulse_exposure")
#' # (2) Create an object of class "survData"
#' dataset <- survData(propiconazole_pulse_exposure)
#' \donttest{
#' # (3) Run the survFit function with TKTD model 'SD' or 'IT' 
#' out <- survFit(dataset , model_type = "SD")
#' # (4) Summarize look the estimated parameters
#' summary(out)
#' # (5) Plot the fitted curve
#' plot(out, adddata = FALSE)
#' # (6) Plot the fitted curve with ggplot style and CI as spaghetti
#' plot(out, spaghetti = TRUE)
#' }
#' 
#' @import rjags
#' 
#' @export
survFit.survDataVarExp <- function(data,
                                 model_type = NULL,
                                 quiet = FALSE,
                                 n.chains = 3,
                                 n.adapt = 1000,
                                 n.iter = NULL,
                                 n.warmup = NULL,
                                 thin.interval = NULL,
                                 limit.sampling = TRUE,
                                 dic.compute = FALSE,
                                 dic.type = "pD",
                                 hb_value = TRUE,
                                 hb_valueFIXED = NA,
                                 extend_time = 100,
                                 ...){
  
  ##
  ## Pre modelling measure and tests
  ##

  ### ensures model_type is one of "SD" and "IT"
  if(is.null(model_type) || ! (model_type %in% c("SD","IT"))) {
    stop("You need to specify a 'model_type' among 'SD' or 'IT'")
  }
  ### check number of sample for the diagnostic procedure
  if (n.chains < 2) {
    stop('2 or more parallel chains required')
  }
  ### warning message when hb_value = NULL
  if(hb_value==FALSE){
    warning("This is not an error message: the parameter 'hb' is fixed. This means that the correlation between
            'hb' and other parameters is ignored.")
    ### Set default hb_valueFIXED
    if(is.na(hb_valueFIXED)){
      hb_valueFIXED = 0
    }
  }

  ##
  ## Data and Priors for model
  ##
  
  globalData <- modelData(x = data, model_type = model_type, extend_time = extend_time)
  
  ### Remove the information of replicate since this is not used in JAGS, and so a warning message would be show
  
  jags.data <- globalData$modelData
  
  jags.data_fit <- jags.data
  
  jags.data_fit$replicate <- NULL
  jags.data_fit$conc <- NULL
  jags.data_fit$replicate_long <- NULL
  
  priorsData = globalData$priorsMinMax
  
  ##
  ## Define model
  ##
  
  if(model_type == "SD"){
    
    jags.data_fit$time = NULL # remove jags.data_fit$time for varSD model
    
    if(hb_value == TRUE){
      ### Determine sampling parameters
      jags.data_fit$hb_value = 1
      jags.data_fit$hb_valueFIXED = -1 # just to have it in JAGS
      parameters_sampling <- c("kd_log10", "hb_log10", "z_log10", "kk_log10")
      parameters <- c("kd_log10", "hb_log10", "z_log10", "kk_log10", "hb", "psurv", "Nsurv_ppc", "Nsurv_sim")
    } else{
      ### Determine sampling parameters
      jags.data_fit$hb_value = 0
      jags.data_fit$hb_valueFIXED = hb_valueFIXED
      parameters_sampling <- c("kd_log10", "z_log10", "kk_log10")
      parameters <- c("kd_log10", "z_log10", "kk_log10", "hb", "psurv", "Nsurv_ppc", "Nsurv_sim")
    }
    file_to_use <- jags_TKTD_varSD
    
      
  } else if(model_type == "IT"){
    ### Determine sampling parameters
    
    if(hb_value == TRUE){
      jags.data_fit$hb_value = 1
      jags.data_fit$hb_valueFIXED = -1 # just to have it in JAGS
      parameters_sampling <- c("kd_log10", "hb_log10","alpha_log10", "beta_log10")
      parameters <- c("kd_log10", "hb_log10","alpha_log10", "beta_log10", "hb", "psurv", "Nsurv_ppc", "Nsurv_sim")
    } else{
      jags.data_fit$hb_value = 0
      jags.data_fit$hb_valueFIXED = hb_valueFIXED
      parameters_sampling <- c("kd_log10","alpha_log10", "beta_log10")
      parameters <- c("kd_log10","alpha_log10", "beta_log10", "psurv", "hb", "Nsurv_ppc", "Nsurv_sim")
    }
    file_to_use <- jags_TKTD_varIT

  }

  model <- survLoadModel(model.program = file_to_use,
                         data = jags.data_fit,
                         n.chains = n.chains,
                         Nadapt = n.adapt,
                         quiet = quiet)

  
  ##
  ## estimate the number of iteration required for convergency of chains
  ## by using the raftery.diag
  ##
  
  if(is.null(n.warmup) | is.null(thin.interval) | is.null(n.iter)){
    
    sampling.parameters <- modelSamplingParameters(model,
                                                   parameters_sampling,
                                                   n.chains = n.chains, quiet = quiet)
    if (sampling.parameters$niter > 5e5)
      stop("The model needs too many iterations to provide reliable parameter estimates !")
    
    n.warmup = sampling.parameters$burnin
    thin.interval = sampling.parameters$thin
    n.iter = sampling.parameters$niter
    
  }

  ### model to check priors with the model
  update(model, n.warmup)
  
  if(dic.compute == TRUE){ # Deviance Information Criterion
    dic <- dic.samples(model,
                       n.iter = n.iter,
                       thin = thin.interval,
                       type = dic.type) 
  } else dic = NULL
  
  mcmc =  coda.samples(model,
                       variable.names = parameters,
                       n.iter = n.iter,
                       thin = thin.interval)
  
  ##
  ## Cheking posterior range with data from experimental design:
  ##
  
  estim.par <- survFit_TKTD_params(mcmc, model_type = model_type, hb_value = hb_value)
  
  warnings <- msgTableCreate()
  
  if (filter(estim.par, parameters == "kd")$Q97.5 > priorsData$kd_max){
    ## store warning in warnings table
    msg <- "The estimation of the dominant rate constant (model parameter kd) lies 
    outside the range used to define its prior distribution which indicates that this
    rate is very high and difficult to estimate from this experiment !"
    warnings <- msgTableAdd(warnings, "kd_outRange", msg)
    ## print the message
    warning(msg, call. = FALSE)
  }
  
  if(hb_value == TRUE){
    if (filter(estim.par, parameters == "hb")$Q2.5 < priorsData$hb_min){
      ## store warning in warnings table
      msg <- "The estimation of the natural instantaneous mortality rate (model 
    parameter hb) lies outside the range used to define its prior distribution 
    which indicates that this rate is very low and so difficult to estimate 
    from this experiment !"
      warnings <- msgTableAdd(warnings, "hb_outRange", msg)
      ## print the message
      warning(msg, call. = FALSE)
    }
  }

  ### for SD model
  if(model_type == "SD"){
    if (filter(estim.par, parameters == "kk")$Q97.5 > priorsData$kk_max){
      ## store warning in warnings table
      msg <- "The estimation of the killing rate (model parameter k) lies
      outside the range used to define its prior distribution which indicates
      that this rate is very high and difficult to estimate from this experiment !"
      warnings <- msgTableAdd(warnings, "kk_outRange", msg)
      ## print the message
      warning(msg, call. = FALSE)
    }
    
    if (filter(estim.par, parameters == "z")$Q2.5 < priorsData$conc_min ||
        filter(estim.par, parameters == "z")$Q97.5 > priorsData$conc_max){
      ## store warning in warnings table
      msg <- "The estimation of Non Effect Concentration threshold (NEC) 
      (model parameter z) lies outside the range of tested concentration 
      and may be unreliable as the prior distribution on this parameter is
      defined from this range !"
      warnings <- msgTableAdd(warnings, "z_outRange", msg)
      ## print the message
      warning(msg, call. = FALSE)
    }
  }
  
  ### for IT model
  if(model_type == "IT"){
    
    if (filter(estim.par, parameters == "alpha")$Q2.5 < priorsData$conc_min ||
        filter(estim.par, parameters == "alpha")$Q97.5 > priorsData$conc_max){
      ## store warning in warnings table
      msg <- "The estimation of log-logistic median (model parameter alpha) 
      lies outside the range of tested concentration and may be unreliable as 
      the prior distribution on this parameter is defined from this range !"
      warnings <- msgTableAdd(warnings, "alpha_outRange", msg)
      ## print the message
      warning(msg, call. = FALSE)
    }
  }
  
  ##
  ## MCMC information
  ## 
  mcmcInfo = data.frame(n.iter = n.iter,
                        n.chains = n.chains,
                        n.adapt = n.adapt,
                        thin.interval = thin.interval,
                        n.warmup = n.warmup)
  

  ##
  ##
  ##
  transformed.data <- data.frame(
    replicate = jags.data$replicate,
    time = jags.data$time,
    conc = jags.data$conc,
    Nsurv = jags.data$Nsurv
  ) %>%
    group_by(replicate) %>%
    mutate(Ninit = max(Nsurv, na.rm = TRUE))
  
  ##
  ## OUTPUT
  ##
  
  OUT <- list(estim.par = estim.par,
              mcmc = mcmc,
              model = model,
              dic = dic,
              parameters = parameters,
              mcmcInfo = mcmcInfo,
              jags.data = jags.data,
              warnings = warnings,
              model_type = model_type,
              transformed.data = transformed.data,
              original.data = data,
              hb_valueFIXED = hb_valueFIXED)
  
  class(OUT) <- c("survFitVarExp","survFit")
  return(OUT)
}