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
.
There is only 4 functions in the gutsRstan
package. We want to keep it simple.
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.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))
model_type
as "SD"
.model_type
as "IT"
.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")
# 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)
# 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")
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)
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"))
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)
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)
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))
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))
The ppc
function allows to check posterior predictions
ppc(stanguts_to_survFit(fit_SD_diaz)) # ppc(stanguts_to_survFit(fit_IT_diaz))
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)
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)
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.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.