tutorial.R
In morseTKTD: Bayesian Inference of TKTD Models

## ----include=FALSE, echo=FALSE------------------------------------------------
knitr::opts_chunk$set(fig.width = 7,
                      fig.height = 4,
                      # cache = TRUE,
                      warning = FALSE)

## ----load_package, echo=FALSE, results='hide', message=FALSE------------------
library(morseTKTD)

## ----TKTDcst------------------------------------------------------------------
# (1) load data set
data(propiconazole)

# (2) check structure and integrity of the data set
survDataCheck(propiconazole)

# (3) create a `survData` object
PRZcst <- survData(propiconazole)

# (4) represent the number of survivors as a function of time
plot(PRZcst)

## ----fitcstSD, echo=TRUE, warning=FALSE---------------------------------------
# (6) fit the TKTD model SD
fit_cstSD <- fit(PRZcst, model_type = "SD", refresh = 0)

## ----priorsParams-------------------------------------------------------------
in_param <- priors_distribution(fit_cstSD)
head(in_param)

## ----fitcst_pp1---------------------------------------------------------------
# prior/posterior estimates
pp_cstSD = priorPosterior(fit_cstSD)
head(pp_cstSD)

## ----plot_pp_cstSD------------------------------------------------------------
# plot of parameters estimates
plot(pp_cstSD)

## ----plot_fit_cstSD-----------------------------------------------------------
plot(fit_cstSD)

## ----plot_fit_cstSD_withData--------------------------------------------------
plot(fit_cstSD, add_data = FALSE)

## ----plot_ppc_cstSD-----------------------------------------------------------
# PPC_fit_cstSD <- ppc(fit_cstSD)
# plot(PPC_fit_cstSD)

## ----extractNsurvPPC----------------------------------------------------------
extNsurvPPC_cstSD <- extract_Nsurv_ppc(fit_cstSD)
head(extNsurvPPC_cstSD)

## ----extractNsurvSIM----------------------------------------------------------
extNsurvSIM_cstSD <- extract_Nsurv_sim(fit_cstSD)
head(extNsurvSIM_cstSD)

## ----fitcstSDFIXhb, warning=FALSE---------------------------------------------
# fit the TKTD model SD with fixed hb value
fit_cstSDFIXhb <- fit(PRZcst, model_type = "SD", hb_value = 0.2, refresh = 0)

## ----plotcstSDFIXhb_PrintVal--------------------------------------------------
plot(fit_cstSDFIXhb)

## ----summarycstSDFIXhbSummary-------------------------------------------------
pp_cstSDFIXhb = priorPosterior(fit_cstSDFIXhb)
head(pp_cstSDFIXhb[pp_cstSDFIXhb$pp == "posterior", ])

## ----fitcstIT-----------------------------------------------------------------
fit_cstIT <- fit(PRZcst, model_type = "IT", refresh = 0)

## ----PLOTpriors_post_cstIT----------------------------------------------------
PP_cstIT <- priorPosterior(fit_cstIT)
plot(PP_cstIT)

## ----plot_fit_cstIT-----------------------------------------------------------
plot(fit_cstIT)

## ----ppc_fit_cstIT------------------------------------------------------------
ppc_cstIT = ppc(fit_cstIT)
head(ppc_cstIT)

## ----plot_ppc_fit_cstIT-------------------------------------------------------
# plot(ppc_cstIT)

## ----predict------------------------------------------------------------------
# (1) upload or build a data frame with the exposure profile
# argument `replicate` is used to provide several profiles of exposure
data_4prediction <- data.frame(time = c(1:10, 1:10),
                               conc = c(c(0,0,40,0,0,0,40,0,0,0),
                                        c(21,19,18,23,20,14,25,8,13,5)),
                               replicate = c(rep("pulse", 10), rep("random", 10)))

# (2) Use the fit on constant exposure propiconazole with model SD (see previously)
predict_PRZ_cstSD_4pred <- predict(fit_cstSD, data_4prediction)

## ----predictPlot--------------------------------------------------------------
# (3) Plot the predicted survival rate under the new exposure profiles.
plot(predict_PRZ_cstSD_4pred)

## ----predictPlotBackground----------------------------------------------------
# (3) Plot the predicted survival rate under the new exposure profiles.
plot(predict_PRZ_cstSD_4pred, background_concentration = TRUE)

## ----predictPlotBackgroundLegend----------------------------------------------
# (3) Plot the predicted survival rate under the new exposure profiles.
plot(predict_PRZ_cstSD_4pred, background_concentration = TRUE, add_legend = TRUE)

## ----MCMCPredicExternal, warning=FALSE----------------------------------------
Nsurv_predict_PRZ_cst4pred = compute_Nsurv(predict_PRZ_cstSD_4pred, Ninit = 10)
Nsurv_predict_PRZ_cst4pred$df_quantile

## ----hb_value-----------------------------------------------------------------
# Use the same data set profile to predict without 'hb'
predict_PRZ_cstSD_4pred_hbOUT <- predict(
  fit = fit_cstSD,
  display.exposure = data_4prediction,
  hb_value = 0)
# Plot the prediction:
plot(predict_PRZ_cstSD_4pred_hbOUT)

## ----hb_valueFIX2-------------------------------------------------------------
# Use the same data set profile to predict without 'hb'
predict_PRZ_cstSD_4pred_hbFIX2 <- predict(
  fit = fit_cstSD,
  display.exposure = data_4prediction,
  hb_value = 2)
# Plot the prediction:
plot(predict_PRZ_cstSD_4pred_hbFIX2)

## ----cstTOcst, warning=FALSE--------------------------------------------------
predict_psurv_PRZ_SD_cstTOcst <- predict(fit_cstSD, propiconazole)
compute_Nsurv_PRZ_SD_cstTOcst <- compute_Nsurv(
  predict_psurv_PRZ_SD_cstTOcst,
  Ninit = compute_Ninit(propiconazole)
)
head(predict_psurv_PRZ_SD_cstTOcst
     $df_quantile)

## ----cstTOvar, warning=FALSE--------------------------------------------------
data(propiconazole_pulse_exposure)
predict_psurv_PRZ_SD_cstTOvar <- predict(fit_cstSD, propiconazole_pulse_exposure)
predict_Nsurv_PRZ_SD_cstTOvar <- compute_Nsurv(
  predict_psurv_PRZ_SD_cstTOvar,
  Ninit = compute_Ninit(propiconazole)
)
head(predict_Nsurv_PRZ_SD_cstTOvar$df_quantile)

## ----checkNsurvPRED_1---------------------------------------------------------
# compute_check(compute_Nsurv_PRZ_SD_cstTOvar)

## ----plotPredict_Nsurv--------------------------------------------------------
# plot(predict_Nsurv_PRZ_SD_cstTOvar)

## ----ppcPredict_Nsurv---------------------------------------------------------
# ppc(predict_Nsurv_PRZ_SD_cstTOvar)

## ----nameEFSA_SD--------------------------------------------------------------
# head(priors_distribution(fit_cstSD, EFSA_name = TRUE))
# plot_prior_post(fit_cstSD, EFSA_name = TRUE)

## ----nameEFSA_IT--------------------------------------------------------------
# head(priors_distribution(fit_cstIT, EFSA_name = TRUE))
# plot_prior_post(fit_cstIT, EFSA_name = TRUE)

## ----cstSDLCx-----------------------------------------------------------------
# LC50 at the maximum time-point:
LCx_cstSD <- lcxt(fit_cstSD, x = 0.5)
plot(LCx_cstSD)

# LC50 at time = 2
lcxt(fit_cstSD, x = 0.5, t = 2) |> plot()
## Note the use of the pipe operator, `|>`, which is a powerful tool for clearly expressing a sequence of multiple operations.
## For more information on pipes, see: http://r4ds.had.co.nz/pipes.html

## ----cstSDLCx_3015------------------------------------------------------------
# LC50 at time = 15
lcxt(fit_cstSD, x = 0.5, t = 15) |> plot()

## ----cstITLCx-----------------------------------------------------------------
# LC50 at the maximum time-point:
LCx_cstIT <- lcxt(fit_cstIT, x = 0.5)
plot(LCx_cstIT)

# LC50 at time = 2
lcxt(fit_cstIT, x = 0.5, t = 2) |> plot()

# LC30 at time = 15
lcxt(fit_cstIT, x = 0.3, t = 15) |> plot()

## ----MFx_compt----------------------------------------------------------------
# (1) upload or build a data frame with the exposure profile
data_4MFx <- data.frame(time = 1:10,
                        conc = c(0,0.5,8,3,0,0,0.5,8,3.5,0),
                        replicate = "A")

# (2) Use the fit on constant exposure propiconazole with model SD (see previously)
MFx_PRZ_cstSD_4MFx <- lpxt(fit = fit_cstSD, display.exposure = data_4MFx)

## ----MFx_plot-----------------------------------------------------------------
# (3) Plot the survival rate as function of the multiplication factors.
plot(MFx_PRZ_cstSD_4MFx)

## ----MFx_plotMfX--------------------------------------------------------------
# (3 bis) Plot the  binary-search of MFx
plot(MFx_PRZ_cstSD_4MFx, plot  = "MFx")

## ----print_df_MFx-------------------------------------------------------------
MFx_PRZ_cstSD_4MFx$df_q50

## ----MFx_x10------------------------------------------------------------------
# (2 bis) fit on constant exposure propiconazole with model SD (see previously)
MFx_PRZ_cstSD_4MFx_x10 <- lpxt(fit = fit_cstSD, x = 0.1, display.exposure = data_4MFx)
plot(MFx_PRZ_cstSD_4MFx_x10)

## ----MFx_x10plot2-------------------------------------------------------------
plot(MFx_PRZ_cstSD_4MFx_x10, plot = "MFx")

## ----MFx_IT-------------------------------------------------------------------
# (2) Use the fit on constant exposure propiconazole with model IT. No print of run messages.
MFx_PRZ_cstIT_4pred <- lpxt(
  fit = fit_cstIT, x = 0.5, t = 4,
  display.exposure = data_4MFx)

# (3) Plot the survival rate versus multiplication factors.
plot(MFx_PRZ_cstIT_4pred)

# (4) Plot the serching
plot(MFx_PRZ_cstIT_4pred, plot = "MFx")

## ----MFx_ITplot---------------------------------------------------------------
# (2) Use the fit on constant exposure propiconazole with model IT. No print of run messages.
MFx_PRZ_cstIT_4pred <- lpxt(fit = fit_cstIT, x = 0.1, t = 4, display.exposure = data_4MFx)
# 
plot(MFx_PRZ_cstIT_4pred)
# 
plot(MFx_PRZ_cstIT_4pred, plot = "MFx")