Examples
In rPBK: Inference and Prediction of Generic Physiologically-Based Kinetic Models

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(rPBK)

Single Compartment and Single Exposure : Male Gammarus Single

data("dataMaleGammarusSingle")
# work only when replicate have the same length !!!
data_MGS <- dataMaleGammarusSingle[dataMaleGammarusSingle$replicate == 1,]

modelData_MGS <- dataPBK(
  object = data_MGS,
  col_time = "time",
  col_replicate = "replicate",
  col_exposure = "expw",
  col_compartment = "conc",
  time_accumulation = 4,
  nested_model = NA)

fitPBK_MGS <- fitPBK(modelData_MGS)

plot(fitPBK_MGS)

library(loo)
log_lik_MGS <- loo::extract_log_lik(fitPBK_MGS$stanfit, merge_chains = FALSE)
WAIC_MGS <- waic(log_lik_MGS)

Multiple Compartiment, Single Exposure - Default interaction

data("dataCompartment4")
data_C4 <- dataCompartment4

modelData_C4 <- dataPBK(
  object = data_C4,
  col_time = "temps",
  col_replicate = "replicat",
  col_exposure = "condition",
  col_compartment = c("intestin", "reste", "caecum", "cephalon"),
  time_accumulation = 7)

You can have a look at the assumption on the interaction

nested_model(modelData_C4)

fitPBK_C4 <- fitPBK(modelData_C4, chains = 1, iter = 1000)

plot(fitPBK_C4)

Compute WAIC with loo library:

library(loo)
log_lik_C4 <- loo::extract_log_lik(fitPBK_C4$stanfit, merge_chains = FALSE)
WAIC_C4 <- waic(log_lik_C4)
print(WAIC_C4)

Compute LOO:

r_eff_C4 <- relative_eff(exp(log_lik_C4))
LOO_C4 <- loo(log_lik_C4, r_eff = r_eff_C4, cores = 2)
print(LOO_C4)

Multiple Compartiment, Single Exposure : Change nesting

You can have a look at the assumption on the interaction

nm_C4 = nested_model(modelData_C4)

We want to change the interaction between organs. For now, all organs interact with each other but not with themselve, the the interaction matrix look like:

nm_C4$k_nest

which can be written like:

matrix(c(
  c(0,1,1,1),
  c(1,0,1,1),
  c(1,1,0,0),
  c(1,1,1,0)),
  ncol=4,nrow=4,byrow=TRUE)

Let assume interaction are only one way, so a triangular matrix:

matrix(c(
  c(0,1,1,1),
  c(0,0,1,1),
  c(0,0,0,0),
  c(0,0,0,0)),
  ncol=4,nrow=4,byrow=TRUE)

modelData_C42 <- dataPBK(
  object = data_C4,
  col_time = "temps",
  col_replicate = "replicat",
  col_exposure = "condition",
  col_compartment = c("intestin", "reste", "caecum", "cephalon"),
  time_accumulation = 7,
  ku_nest = c(1,1,1,1), # No Change here
  ke_nest = c(1,1,1,1), # No Change here
  k_nest = matrix(c(
            c(0,1,1,1),
            c(0,0,1,1),
            c(0,0,0,0),
            c(0,0,0,0)),
            ncol=4,nrow=4,byrow=TRUE) # Remove 
  )

nested_model(modelData_C42)

fitPBK_C42 <- fitPBK(modelData_C42, chains = 1, iter = 1000)

plot(fitPBK_C42)

log_lik_C42 <- loo::extract_log_lik(fitPBK_C42$stanfit, merge_chains = FALSE)
WAIC_C42 <- waic(log_lik_C42)
print(WAIC_C42)

Compare WAIC with previous model

comp_C4_C42 <- loo_compare(WAIC_C4, WAIC_C42)
print(comp_C4_C42)

The first column shows the difference in ELPD relative to the model with the largest ELPD. In this case, the difference in elpd and its scale relative to the approximate standard error of the difference) indicates a preference for the second model (model2).