R/fit-SurvData.R

In morseTKTD: Bayesian Inference of TKTD Models

#' @name SurvFit
#' 
#' @title Fits a TKTD model for survival analysis using Bayesian inference
#'
#' @description
#' This function estimates the parameters of a TKTD model ('SD' or 'IT')
#' for survival analysis using Bayesian inference. In this model,
#' the survival rate of individuals is modeled as a function of the chemical compound
#' concentration with a mechanistic description of the effects on survival over
#' time.
#' 
#' The function returns the parameter estimates of 
#' Toxicokinetic-toxicodynamic (TKTD) models
#' \code{SD} for 'Stochastic Death' or \code{IT} fo 'Individual Tolerance'.
#' TKTD models, and particularly the General Unified Threshold model of
#' Survival (GUTS), provide a consistent process-based
#' framework to analyse both time and concentration dependent datasets.
#' In GUTS-SD, all organisms are assumed to have the same internal concentration 
#' threshold (denoted \eqn{z}), and, once exceeded, the instantaneous probability
#' to die increases linearly with the internal concentration.
#' In GUTS-IT, the threshold concentration is distributed among all the organisms, and once 
#' exceeded in one individual, this individual dies immediately.
#'
#' @param data An object of class \code{survDataCstExp} or \code{survDataVarExp}.
#' @param model_type Can be \code{"SD"} or \code{"IT"} to choose
#'   between "Stochastic Death" or "Individual Tolerance" models
#'   (resp.). See the modeling vignette for details.
#' @param hb_value If \code{TRUE}, the background mortality \code{hb} is taken into account.
#' If \code{FALSE}, parameter \code{hb} is set to 0. The default is \code{TRUE}.
#' @param \dots Further arguments to be passed to generic methods
#' using argument of \link[rstan]{sampling} function.
#' 
#' @references Jager, T., Albert, C., Preuss, T. G. and Ashauer, R. (2011) 
#' General unified threshold model of survival-a toxicokinetic-toxicodynamic
#'  framework for ecotoxicology, \emph{Environmental Science and Technology}, 45, 2529-2540.
#' 303-314.
#' 
#' @keywords estimation 
#' 
#' @return An object of class `stanfit` returned by `rstan::sampling`
#' 
#' @export
#' 
fit <- function(data,
                model_type,
                hb_value,
                ...){
  UseMethod("fit")
}

#' @name SurvFit
#' @export
fit.SurvDataCstExp <- function(data, model_type, hb_value = NULL, ...) {
  standata = modelData(data, model_type, hb_value)
  if (model_type == "SD") {
    outfit <- rstan::sampling(stanmodels$cstSD, data = standata, ...)
  }
  if (model_type == "IT") {
    outfit <- rstan::sampling(stanmodels$cstIT, data = standata, ...)
  }
  out <- list(standata = standata, stanfit = outfit, model_type = model_type)
  class(out) <- append(c("SurvFitCstExp", "SurvFit"), class(out))
  return(out)
}

#' @name SurvFit
#' @export
fit.SurvDataVarExp <- function(data, model_type,  hb_value = NULL, ...) {
  standata = modelData(data, model_type, hb_value)
  if (model_type == "SD") {
    outfit <- rstan::sampling(stanmodels$varSD, data = standata, ...)
  }
  if (model_type == "IT") {
    outfit <- rstan::sampling(stanmodels$varIT, data = standata, ...)
  }
  out <- list(standata = standata, stanfit = outfit, model_type = model_type)
  class(out) <- append(c("SurvFitVarExp", "SurvFit"), class(out))
  return(out)
}


#############################################################################
#
#    survFit_TKTD_params
#
#############################################################################
#' @title Table of posterior estimated parameters
#' 
#' @description
#' create the table of posterior estimated parameters for the survival analyses
#' 
#' @param mcmc list of estimated parameters for the model with each item
#'  representing a chain
#' @param model_type model type `SD` or `IT` 
#' @param hb_value TRUE or FALSE, conservning the use of hb parameter in the model.
#' 
#' @return a `data.frame` with 3 columns (values, CIinf, CIsup) and
#'  3-4rows (the estimated parameters)
survFit_TKTD_params <- function(mcmc, model_type, hb_value = TRUE) {

    # Retrieving parameters of the model
    res.M <- summary(mcmc)
    
    kd <- 10^res.M$quantiles["kd_log10", "50%"]
    kd_inf95 <- 10^res.M$quantiles["kd_log10", "2.5%"]
    kd_sup95 <- 10^res.M$quantiles["kd_log10", "97.5%"]
    
    if (hb_value == TRUE) {
      hb <- 10^res.M$quantiles["hb_log10", "50%"]
      hb_inf95 <- 10^res.M$quantiles["hb_log10", "2.5%"]
      hb_sup95 <- 10^res.M$quantiles["hb_log10", "97.5%"]
    }
    
    if (model_type == "SD") {
      kk <- 10^res.M$quantiles["kk_log10", "50%"]
      kk_inf95 <- 10^res.M$quantiles["kk_log10", "2.5%"]
      kk_sup95 <- 10^res.M$quantiles["kk_log10", "97.5%"]
      
      z <- 10^res.M$quantiles["z_log10", "50%"]
      z_inf95 <- 10^res.M$quantiles["z_log10", "2.5%"]
      z_sup95 <- 10^res.M$quantiles["z_log10", "97.5%"]
      
      if (hb_value == TRUE) {
        res <- data.frame(parameters = c("kd", "hb", "z", "kk"),
                          median = c(kd, hb, z, kk),
                          Q2.5 = c(kd_inf95, hb_inf95, z_inf95, kk_inf95),
                          Q97.5 = c(kd_sup95, hb_sup95, z_sup95, kk_sup95))
      } else{
        res <- data.frame(parameters = c("kd", "z", "kk"),
                          median = c(kd, z, kk),
                          Q2.5 = c(kd_inf95, z_inf95, kk_inf95),
                          Q97.5 = c(kd_sup95, z_sup95, kk_sup95))
      }
      
    } else if (model_type == "IT") {
      alpha <- 10^res.M$quantiles["alpha_log10", "50%"]
      alpha_inf95 <- 10^res.M$quantiles["alpha_log10", "2.5%"]
      alpha_sup95 <- 10^res.M$quantiles["alpha_log10", "97.5%"]
      
      beta <- 10^res.M$quantiles["beta_log10", "50%"]
      beta_inf95 <- 10^res.M$quantiles["beta_log10", "2.5%"]
      beta_sup95 <- 10^res.M$quantiles["beta_log10", "97.5%"]
      
      if (hb_value == TRUE) {
        res <- data.frame(parameters = c("kd", "hb", "alpha", "beta"),
                          median = c(kd, hb, alpha, beta),
                          Q2.5 = c(kd_inf95, hb_inf95, alpha_inf95, beta_inf95),
                          Q97.5 = c(kd_sup95, hb_sup95, alpha_sup95, beta_sup95))
      } else{
        res <- data.frame(parameters = c("kd", "alpha", "beta"),
                          median = c(kd, alpha, beta),
                          Q2.5 = c(kd_inf95, alpha_inf95, beta_inf95),
                          Q97.5 = c(kd_sup95, alpha_sup95, beta_sup95))
      }
    } else {
      stop("please, provide the 'model_type': 'IT' or 'SD'")
    }
    
    return(res)
}