Examples

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)
#> 
#> SAMPLING FOR MODEL 'PBK_AD' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 8.6e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.86 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.548 seconds (Warm-up)
#> Chain 1:                0.344 seconds (Sampling)
#> Chain 1:                0.892 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'PBK_AD' NOW (CHAIN 2).
#> Chain 2: 
#> Chain 2: Gradient evaluation took 4.1e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.41 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2: 
#> Chain 2: 
#> Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 2: 
#> Chain 2:  Elapsed Time: 0.994 seconds (Warm-up)
#> Chain 2:                3.9 seconds (Sampling)
#> Chain 2:                4.894 seconds (Total)
#> Chain 2: 
#> 
#> SAMPLING FOR MODEL 'PBK_AD' NOW (CHAIN 3).
#> Chain 3: 
#> Chain 3: Gradient evaluation took 3.6e-05 seconds
#> Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.36 seconds.
#> Chain 3: Adjust your expectations accordingly!
#> Chain 3: 
#> Chain 3: 
#> Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 3: 
#> Chain 3:  Elapsed Time: 0.528 seconds (Warm-up)
#> Chain 3:                0.572 seconds (Sampling)
#> Chain 3:                1.1 seconds (Total)
#> Chain 3: 
#> 
#> SAMPLING FOR MODEL 'PBK_AD' NOW (CHAIN 4).
#> Chain 4: 
#> Chain 4: Gradient evaluation took 4.5e-05 seconds
#> Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.45 seconds.
#> Chain 4: Adjust your expectations accordingly!
#> Chain 4: 
#> Chain 4: 
#> Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 4: 
#> Chain 4:  Elapsed Time: 0.933 seconds (Warm-up)
#> Chain 4:                0.898 seconds (Sampling)
#> Chain 4:                1.831 seconds (Total)
#> Chain 4:
#> Warning: There were 2088 divergent transitions after warmup. See
#> https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> to find out why this is a problem and how to eliminate them.
#> Warning: Examine the pairs() plot to diagnose sampling problems
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#bulk-ess
#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#tail-ess
plot(fitPBK_MGS)
#> Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
#> ℹ Please use tidy evaluation idioms with `aes()`.
#> ℹ See also `vignette("ggplot2-in-packages")` for more information.
#> ℹ The deprecated feature was likely used in the rPBK package.
#>   Please report the issue at
#>   <https://gitlab.in2p3.fr/mosaic-software/rPBK/-/issues>.
#> This warning is displayed once per session.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.

library(loo)
#> This is loo version 2.10.1
#> - Online documentation and vignettes at mc-stan.org/loo
#> - As of v2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session.
log_lik_MGS <- loo::extract_log_lik(fitPBK_MGS$stanfit, merge_chains = FALSE)
WAIC_MGS <- waic(log_lik_MGS)
#> Warning: 
#> 1 (12.5%) p_waic estimates greater than 0.4. We recommend trying loo instead.

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)
#> $ku_nest
#> uptake intestin    uptake reste   uptake caecum uptake cephalon 
#>               1               1               1               1 
#> 
#> $ke_nest
#> excretion intestin    excretion reste   excretion caecum excretion cephalon 
#>                  1                  1                  1                  1 
#> 
#> $k_nest
#>          intestin reste caecum cephalon
#> intestin        0     1      1        1
#> reste           1     0      1        1
#> caecum          1     1      0        1
#> cephalon        1     1      1        0
fitPBK_C4 <- fitPBK(modelData_C4, chains = 1, iter = 1000)
#> 
#> SAMPLING FOR MODEL 'PBK_AD' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.000514 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 5.14 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
#> Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
#> Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
#> Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
#> Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
#> Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
#> Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
#> Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 5.044 seconds (Warm-up)
#> Chain 1:                10.349 seconds (Sampling)
#> Chain 1:                15.393 seconds (Total)
#> Chain 1:
#> Warning: There were 500 divergent transitions after warmup. See
#> https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> to find out why this is a problem and how to eliminate them.
#> Warning: There were 1 chains where the estimated Bayesian Fraction of Missing Information was low. See
#> https://mc-stan.org/misc/warnings.html#bfmi-low
#> Warning: Examine the pairs() plot to diagnose sampling problems
#> Warning: The largest R-hat is NA, indicating chains have not mixed.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#r-hat
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#bulk-ess
#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#tail-ess
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)
#> Warning: 
#> 6 (7.1%) p_waic estimates greater than 0.4. We recommend trying loo instead.
print(WAIC_C4)
#> 
#> Computed from 500 by 84 log-likelihood matrix.
#> 
#>           Estimate   SE
#> elpd_waic   -324.7 12.3
#> p_waic        13.5  2.9
#> waic         649.4 24.6
#> 
#> 6 (7.1%) p_waic estimates greater than 0.4. We recommend trying loo instead.

Compute LOO:

r_eff_C4 <- relative_eff(exp(log_lik_C4))
LOO_C4 <- loo(log_lik_C4, r_eff = r_eff_C4, cores = 2)
#> Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
print(LOO_C4)
#> 
#> Computed from 500 by 84 log-likelihood matrix.
#> 
#>          Estimate   SE
#> elpd_loo   -324.1 12.2
#> p_loo        12.9  2.5
#> looic       648.2 24.3
#> ------
#> MCSE of elpd_loo is NA.
#> MCSE and ESS estimates assume MCMC draws (r_eff in [0.0, 0.1]).
#> 
#> Pareto k diagnostic values:
#>                           Count Pct.    Min. ESS
#> (-Inf, 0.63]   (good)     79    94.0%   1       
#>    (0.63, 1]   (bad)       4     4.8%   <NA>    
#>     (1, Inf)   (very bad)  1     1.2%   <NA>    
#> See help('pareto-k-diagnostic') for details.

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
#>          intestin reste caecum cephalon
#> intestin        0     1      1        1
#> reste           1     0      1        1
#> caecum          1     1      0        1
#> cephalon        1     1      1        0

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)
#>      [,1] [,2] [,3] [,4]
#> [1,]    0    1    1    1
#> [2,]    1    0    1    1
#> [3,]    1    1    0    0
#> [4,]    1    1    1    0

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)
#>      [,1] [,2] [,3] [,4]
#> [1,]    0    1    1    1
#> [2,]    0    0    1    1
#> [3,]    0    0    0    0
#> [4,]    0    0    0    0
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)
#> $ku_nest
#> uptake intestin    uptake reste   uptake caecum uptake cephalon 
#>               1               1               1               1 
#> 
#> $ke_nest
#> excretion intestin    excretion reste   excretion caecum excretion cephalon 
#>                  1                  1                  1                  1 
#> 
#> $k_nest
#>          intestin reste caecum cephalon
#> intestin        0     1      1        1
#> reste           0     0      1        1
#> caecum          0     0      0        0
#> cephalon        0     0      0        0
fitPBK_C42 <- fitPBK(modelData_C42, chains = 1, iter = 1000)
#> 
#> SAMPLING FOR MODEL 'PBK_AD' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.000727 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 7.27 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
#> Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
#> Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
#> Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
#> Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
#> Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
#> Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
#> Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 7.515 seconds (Warm-up)
#> Chain 1:                118.598 seconds (Sampling)
#> Chain 1:                126.113 seconds (Total)
#> Chain 1:
#> Warning: There were 453 divergent transitions after warmup. See
#> https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> to find out why this is a problem and how to eliminate them.
#> Warning: There were 47 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
#> https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
#> Warning: Examine the pairs() plot to diagnose sampling problems
#> Warning: The largest R-hat is NA, indicating chains have not mixed.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#r-hat
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#bulk-ess
#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#tail-ess
plot(fitPBK_C42)

log_lik_C42 <- loo::extract_log_lik(fitPBK_C42$stanfit, merge_chains = FALSE)
WAIC_C42 <- waic(log_lik_C42)
#> Warning: 
#> 8 (9.5%) p_waic estimates greater than 0.4. We recommend trying loo instead.
print(WAIC_C42)
#> 
#> Computed from 500 by 84 log-likelihood matrix.
#> 
#>           Estimate   SE
#> elpd_waic   -264.1 18.2
#> p_waic        11.3  2.1
#> waic         528.3 36.5
#> 
#> 8 (9.5%) p_waic estimates greater than 0.4. We recommend trying loo instead.

Compare WAIC with previous model

comp_C4_C42 <- loo_compare(WAIC_C4, WAIC_C42)
print(comp_C4_C42)
#>   model elpd_diff se_diff p_worse diag_diff diag_elpd
#>  model2       0.0     0.0      NA                    
#>  model1     -60.6     9.2    1.00   N < 100
#> 
#> Diagnostic flags present.
#> See ?`loo-glossary` (sections `diag_diff` and `diag_elpd`)
#> or https://mc-stan.org/loo/reference/loo-glossary.html.

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).