In virgile-baudrot/rstanTKTD: Bayesian fitting of the General Unified Threshold model of Survival (GUTS) for data in ecotoxicology

knitr::opts_chunk$set(
  fig.width = 7,
  fig.height = 4,
  cache = TRUE,
  collapse = TRUE,
  comment = "#>"
)

library(gutsRstan)

The package gutsRstan is devoted to the analysis of data from standard toxicity tests. It provides a simple workflow to calibrate GUTS models (Jager et al., 2011 ; Jager and Ashauer, 2018). This document illustrates a typical use of gutsRstan on survival data, which can be followed step-by-step to analyze new data.

Analysis of results by visualization is possible with the package morse.

Bayesian inference of GUTS models

There is only 4 functions in the gutsRstan package. We want to keep it simple.

• One function for the Bayesian inference of GUTS models. In this initial version we only impleted the GUTS-RED versions. We hope to implement all other verions as exposed in the GUTS book in the future.
• One function to quickly plot the fit, compare your data with the fit.
• Two functions to transform a stanguts object into: (i) a survFit object to be used with the morse package ; or (ii) a stanfit object to use the set of functions of the rstan package and all other functions from packages developped around the rstan environment.

Loading data

Here is a typical session to analyse concentration-dependent time-course data using the General Unified Threshold model of Survival (GUTS) (see Jager et al., 2011 or Jager and Ashauer, 2018).

# (1) load dataset
data("data_Diazinon")

The data set can be checked using the morse package where a function for GUTS model is also implemented using the JAGS language for Bayesian inference.

## ---- OPTIONAL - using 'morse' package 
library(morse)
# (2) : check structure and integrity of the dataset, with the package morse
survDataCheck(data_Diazinon)
# (3) OPTIONAL - using package 'morse': represent the number of survivors as a function of time
plot(survData(data_Diazinon), pool.replicate = FALSE)
# (4) check information on the experimental design
summary(survData(data_Diazinon))

Fitting models

• To fit the Stochastic Death model, GUTS-SD, we have to specify the model_type as "SD".
• To fit the Individual Tolerance model, GUTS-IT, we have to specify the model_type as "IT".
• To fit the general model, GUTS-PROPER, we have to specify the model_type as "PROPER".

Both model IT and PROPER may be fitted with either the loglogistic or the lognormal distribution. The default distribution is loglogistic. The lognormal distribution can be specified with the argument distribution.

Before fitting models, you can set the number of cores to use when executing the Markov chains in parallel, which defaults to 1 but we recommend setting the mc.cores option to be as many processors as the hardware and RAM allow (up to the number of chains).

For instance, within gutsRstan the default number of chains is 3, so we recommend setting mc.cores = 3, if your machine has at least 3 cores which should be the case in a recent computer.

# (5) OPTION for the number of cores:
options(mc.cores = 3)

Fitting a data set is easy with the function stan_guts().

The argument model_type is required with "SD", "IT" or "PROPER", there is no default value.

The argument distribution matters for "IT" or "PROPER" models. The default is "loglogistic", the other distribution is "lognormal". If you would like to test some other distributions not implemented yet, we would be happy to add any new ones: please let us know by opening a new github issue.

Here is the single line for fitting an "SD" model on data set "data_Diazinon". We put all the outputs in an object named fit_SD_diaz.

# (6-SD) fit the TK-TD model SD
fit_SD_diaz <- stan_guts(data_Diazinon, model_type = "SD")

Here is the line for a fit of an "IT" model on data set "data_Diazinon". The fit object is named fit_IT_diaz.

# (6-IT_ll) fit the GUTS model IT (default: with distribution 'loglogistic')
fit_IT_diaz <- stan_guts(data_Diazinon, model_type = "IT")

The IT model can also be fitted with the "lognormal" distribution.

# (6-IT_ln) fit the GUTS model IT with distribution 'lognormal'
fit_IT_lN_diaz <- stan_guts(data_Diazinon, model_type = "IT", distribution = "lognormal")

Here is the same line for a fit with thz "PROPER" model and the default "loglogistic" distribution.

# (6-PROPER_ll) fit the GUTS model PROPER (default with distribution 'loglogistic')
fit_PROPER_ll_diaz <- stan_guts(data_Diazinon, model_type = "PROPER")

And finally, here is the line for a fit with "PROPER" model using the "lognormal" distribution.

# (6-PROPER_ln) fit the GUTS model PROPER with distribution 'lognormal'
fit_PROPER_lN_diaz <- stan_guts(data_Diazinon, model_type = "PROPER", distribution = "lognormal")

Plots of fitting results based on internal functions

Plot of survival rate

# Survival rate as a function of time
plot_stanguts(fit_SD_diaz)

plot_stanguts(fit_IT_diaz)
plot_stanguts(fit_IT_lN_diaz)
plot_stanguts(fit_PROPER_lN_diaz)
plot_stanguts(fit_PROPER_ll_diaz)

Plot of number of survivors

# Number of survivors as a function of time
plot_stanguts(fit_SD_diaz, data_type = "Number")

plot_stanguts(fit_IT_diaz, data_type = "Number")
plot_stanguts(fit_IT_lN_diaz, data_type = "Number")
plot_stanguts(fit_PROPER_lN_diaz, data_type = "Number")
plot_stanguts(fit_PROPER_ll_diaz, data_type = "Number")

Plots of fitting results based on functions from the rstan package.

The function stanguts_to_stanfit() allows to use the output of the function rstan_guts() with packages developped around the R interface to the Stan C++ library for Bayesian inference rstan, shinystan, bayesplot, loo packages

# Built of stanfit object
stanfit_SD_diaz <- stanguts_to_stanfit(fit_SD_diaz)
stanfit_IT_diaz <- stanguts_to_stanfit(fit_IT_diaz)
stanfit_IT_lN_diaz <- stanguts_to_stanfit(fit_IT_lN_diaz)
stanfit_PROPER_lN_diaz <- stanguts_to_stanfit(fit_PROPER_lN_diaz)
stanfit_PROPER_ll_diaz <- stanguts_to_stanfit(fit_PROPER_ll_diaz)

library(rstan)

Explore the MCMC

From the R-package rstan, there are many functions to scrutinize the inference process. A simple one is the pairs plot:

## Note: to reduce the size of the vignettes, we did not saved previous objects from the 'stanguts_to_stanfit()' function, so we recall this function in pairs plots.
# ---
pairs(stanguts_to_stanfit(fit_SD_diaz),
      pars = c("hb_log10", "kd_log10", "z_log10", "kk_log10"))

pairs(stanguts_to_stanfit(fit_IT_diaz),
      pars = c("hb_log10", "kd_log10", "alpha_log10", "beta_log10"))
pairs(stanguts_to_stanfit(fit_IT_lN_diaz),
      pars = c("hb_log10", "kd_log10", "alpha_log10", "beta_log10"))
pairs(stanguts_to_stanfit(fit_PROPER_lN_diaz),
      pars = c("hb_log10", "kd_log10", "kk_log10", "alpha_log10", "beta_log10"))
pairs(stanguts_to_stanfit(fit_PROPER_ll_diaz),
      pars = c("hb_log10", "kd_log10", "kk_log10", "alpha_log10", "beta_log10"))

Explore outputs through Shinystan webpages

You can also explore the fitting results with the R-package shinystan :

library(shinystan)
launch_shinystan(stanfit_SD_diaz)
launch_shinystan(stanfit_IT_diaz)
launch_shinystan(stanfit_PROPER_lN_diaz)
launch_shinystan(stanfit_PROPER_ll_diaz)

Plots of fitting results based on functions from package morse

At that time, PROPER models are not yet handled by the stable CRAN version of the morse package. So we only return output for SD and IT models.

First of all, we have to convert the object of class stanguts into an object of class survFit:

survFit_SD_diaz <- stanguts_to_survFit(fit_SD_diaz)
survFit_IT_diaz <- stanguts_to_survFit(fit_IT_diaz)

Summary

library(morse)

The summary function provides parameter estimates as medians and 95\% credible intervals.

## Note: to reduce the size of the vignettes, we did not saved previous object from 'stanguts_to_survFit()' function, so we recall this function
# ---
summary(stanguts_to_survFit(fit_SD_diaz))
# summary(stanguts_to_survFit(fit_IT_diaz))

Plot

The plot function provides a representation of the fitting results for each replicate

# plot(stanguts_to_survFit(fit_SD_diaz))
# plot(stanguts_to_survFit(fit_IT_diaz))

PPC

The ppc function allows to check posterior predictions

ppc(stanguts_to_survFit(fit_SD_diaz))
# ppc(stanguts_to_survFit(fit_IT_diaz))

Lethal concentration prediction

Compared to the target time analysis, TKTD modelling allows to compute and plot the lethal concentration for any x percentage and at any time-point. The chosen time-point can be specified with time_LCx ; by default the maximal time-point in the data set is used.

# LC50 at the final time-point:
LC50_SD_diaz <- LCx(stanguts_to_survFit(fit_SD_diaz),
                    X = 50)
plot(LC50_SD_diaz)
# LC30 at the time 4:
LC30t4_SD_diaz <- LCx(stanguts_to_survFit(fit_SD_diaz),
                      X = 30, time_LCx = 4)
plot(LC30t4_SD_diaz)

Multiplication Factor

data_4predict <- data.frame(time = c(1, 2, 2.1, 5.9, 6, 7, 7.1, 10),
                           conc = c(10, 10, 0, 0, 10, 10, 0, 0))
plot(x = data_4predict$time, y = data_4predict$conc, type = "l",
     las = 1, xlab = "Time", ylab = "Concentration")

# MF50 at the maximum time-point:
MF50_SD_diaz <- MFx(stanguts_to_survFit(fit_SD_diaz),
                    data_predict = data_4predict, X = 50)
plot(MF50_SD_diaz)
# MF30 at the maximum time 20:
MF30t20_SD_diaz <- MFx(stanguts_to_survFit(fit_SD_diaz),
                      data_predict = data_4predict,
                      X = 30, time_MFx = 7)
plot(MF30t20_SD_diaz)

References

Jager, T., Albert, C., Preuss, T. G. and Ashauer, R. (2011) General unified threshold model of survival - a toxicokinetic-toxicodynamic framework for ecotoxicology, Environmental Science & Technology, 45, 2529-2540.

Jager, T. and Ashauer, R. (2018) Modelling survival under chemical stress. A comprehensive guide to the GUTS framework. Version 1.0., Leanpub.

virgile-baudrot/rstanTKTD documentation built on May 29, 2019, 9:57 a.m.