Getting Started with BayesianDEB

Overview

BayesianDEB provides a complete Bayesian framework for Dynamic Energy Budget (DEB) modelling. It wraps pre-written Stan models with a clean R interface for data preparation, model specification, fitting, diagnostics, and visualisation.

The package implements four model types:

Type	Stan model	Description
`individual`	2-state (E, V)	Single individual growth
`growth_repro`	3-state (E, V, R)	Growth + reproduction
`hierarchical`	2-state + random effects	Multi-individual with partial pooling
`debtox`	4-state (E, V, R, D)	Toxicokinetic-toxicodynamic

Installation

# Install cmdstanr (required backend)
install.packages("cmdstanr",
  repos = c("https://stan-dev.r-universe.dev", getOption("repos")))
cmdstanr::install_cmdstan()

# Install BayesianDEB
# remotes::install_github("sciom/BayesianDEB")

The chunks below are executed only when cmdstanr and a working CmdStan installation are available — otherwise they are skipped. To keep build time within JSS limits we use 2 chains and 300 + 300 iterations, which is sufficient for illustration but not for publication-grade inference. Re-run with chains = 4, iter_warmup = iter_sampling = 1000 for production fits.

Example: Individual Growth Model

1. Prepare Data

data(eisenia_growth)

# Use first individual
df1 <- eisenia_growth[eisenia_growth$id == 1, ]
dat <- bdeb_data(growth = df1, f_food = 1.0)
dat
#> 
#> ── BDEB Data ──
#> 
#> ℹ Individuals: 1
#> ℹ Endpoints: growth
#> ℹ Functional response (f): 1
#> → Growth: 13 observations, t = [0, 84]

2. Specify Model

mod <- bdeb_model(
  data   = dat,
  type   = "individual",
  priors = list(
    p_Am    = prior_lognormal(mu = 1.5, sigma = 0.5),
    p_M     = prior_lognormal(mu = -1.0, sigma = 0.5),
    kappa   = prior_beta(a = 3, b = 2),
    sigma_L = prior_halfnormal(sigma = 0.05)
  )
)
mod
#> 
#> ── BDEB Model Specification ──
#> 
#> ℹ Type: individual
#> ℹ Stan model: bdeb_individual_growth
#> ℹ Individuals: 1
#> ℹ Endpoints: growth
#> 
#> ── Priors
#> →   p_Am: LogNormal(1.5, 0.5)
#> →   p_M: LogNormal(-1.0, 0.5)
#> →   kappa: Beta(3.0, 2.0)
#> →   sigma_L: HalfNormal(0.05)
#> →   v: LogNormal(-1.5, 1.0)
#> →   E_G: LogNormal(6.0, 1.0)
#> →   E0: LogNormal(0.0, 1.0)
#> →   L0: LogNormal(-2.0, 1.0)

Unspecified priors are filled from prior_default("individual").

3. Fit Model

fit <- bdeb_fit(mod,
                chains = 2, iter_warmup = 300, iter_sampling = 300,
                seed = 123, refresh = 100)
#> ℹ Compiling Stan model: 'bdeb_individual_growth'
#> ℹ Running MCMC (2 chains, 300 iterations each)
#> Running MCMC with 2 parallel chains...
#> 
#> Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
#> Chain 2 Iteration:   1 / 600 [  0%]  (Warmup)
#> Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 2 Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmp1KA7pd/model-c47577a84097.stan', line 111, column 4 to column 35)
#> Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 2
#> Chain 1 Iteration: 100 / 600 [ 16%]  (Warmup) 
#> Chain 2 Iteration: 100 / 600 [ 16%]  (Warmup) 
#> Chain 1 Iteration: 200 / 600 [ 33%]  (Warmup) 
#> Chain 2 Iteration: 200 / 600 [ 33%]  (Warmup) 
#> Chain 2 Iteration: 300 / 600 [ 50%]  (Warmup) 
#> Chain 2 Iteration: 301 / 600 [ 50%]  (Sampling) 
#> Chain 1 Iteration: 300 / 600 [ 50%]  (Warmup) 
#> Chain 1 Iteration: 301 / 600 [ 50%]  (Sampling) 
#> Chain 2 Iteration: 400 / 600 [ 66%]  (Sampling) 
#> Chain 1 Iteration: 400 / 600 [ 66%]  (Sampling) 
#> Chain 2 Iteration: 500 / 600 [ 83%]  (Sampling) 
#> Chain 1 Iteration: 500 / 600 [ 83%]  (Sampling) 
#> Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
#> Chain 2 finished in 4.5 seconds.
#> Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
#> Chain 1 finished in 5.0 seconds.
#> 
#> Both chains finished successfully.
#> Mean chain execution time: 4.7 seconds.
#> Total execution time: 5.1 seconds.
fit
#> 
#> ── BDEB Fit ──
#> 
#> ℹ Model type: individual
#> ℹ Algorithm: sampling (NUTS)
#> ℹ Chains: 2, Warmup: 300, Sampling: 300
#> ✔ No divergent transitions

4. Diagnostics

bdeb_diagnose(fit)
#> 
#> ── BDEB Diagnostics (individual) ──
#> 
#> ✔ No divergent transitions.
#> ✔ Treedepth OK.
#> ✔ E-BFMI OK (all > 0.3).
#> 
#> ── Parameter Summary
#> ✔ All R-hat < 1.01.
#> ! Low bulk ESS (<400) for: E_G, L0, sigma_L, x_sol[1,2], x_sol[2,2], x_sol[3,2], x_sol[9,2], x_sol[10,2], x_sol[11,2], x_sol[12,2], x_sol[13,2], L_hat[1], L_hat[2], L_hat[3], L_hat[9], L_hat[10], L_hat[11], L_hat[12], L_hat[13]
#>  variable     mean       sd   median       5%      95% rhat ess_bulk ess_tail
#>      p_Am   5.2470 2.86e+00   4.5464   2.1325 1.10e+01    1      708      351
#>       p_M   0.4251 2.30e-01   0.3677   0.1679 8.67e-01    1      745      421
#>     kappa   0.5942 1.97e-01   0.6011   0.2413 9.14e-01    1      639      254
#>         v   0.3442 4.15e-01   0.2309   0.0647 9.52e-01    1      450      415
#>       E_G 465.0713 3.42e+02 387.3190 110.9284 1.08e+03    1      336      284
#>        E0   1.5403 1.58e+00   1.0254   0.2406 4.93e+00    1      692      441
#>        L0   0.1097 1.09e-02   0.1093   0.0930 1.28e-01    1      327      136
#>   sigma_L   0.0185 4.75e-03   0.0173   0.0125 2.77e-02    1      362      373
#> ℹ 39 latent-state rows hidden; use `print(x, full = TRUE)` or `summary(x)$table` to see all.
plot(fit, type = "trace")

# `bayesplot::mcmc_pairs` requires gridExtra (a Suggests of bayesplot).
plot(fit, type = "pairs", pars = c("p_Am", "p_M", "kappa"))

5. Posterior Predictive Checks

ppc <- bdeb_ppc(fit, type = "growth")
plot(ppc)

6. Derived Quantities

bdeb_derived(fit, quantities = c("L_inf", "growth_rate"))
#> # A draws_df: 600 iterations, 1 chains, and 2 variables
#>    L_inf growth_rate
#> 1   13.6     0.00018
#> 2    4.6     0.00065
#> 3   11.5     0.00043
#> 4    6.7     0.00021
#> 5   16.3     0.00068
#> 6    3.1     0.00011
#> 7    8.3     0.00092
#> 8    7.9     0.00030
#> 9   13.5     0.00038
#> 10   3.4     0.00022
#> # ... with 590 more draws
#> # ... hidden reserved variables {'.chain', '.iteration', '.draw'}

7. Trajectory Plot

plot(fit, type = "trajectory", n_draws = 200)

Example: Hierarchical Model (Multiple Individuals)

For illustration we use a 5-individual subset; the full multi-individual analysis with all 21 individuals is in the replication archive.

dat_all <- bdeb_data(
  growth = eisenia_growth[eisenia_growth$id %in% 1:5, ],
  f_food = 1.0
)

mod_hier <- bdeb_model(dat_all, type = "hierarchical")

fit_hier <- bdeb_fit(mod_hier,
                     chains = 2, iter_warmup = 300, iter_sampling = 300,
                     seed = 123, refresh = 100)
#> ℹ Compiling Stan model: 'bdeb_hierarchical_growth'
#> ℹ Running MCMC (2 chains, 300 iterations each)
#> Running MCMC with 2 parallel chains...
#> 
#> Chain 1 Iteration:   1 / 600 [  0%]  (Warmup)
#> Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 1 Exception: Exception: ode_bdf_tol: ode parameters and data is inf, but must be finite! (in '/tmp/Rtmp1KA7pd/model-c47538a8e01f.stan', line 51, column 6 to line 55, column 8) (in '/tmp/Rtmp1KA7pd/model-c47538a8e01f.stan', line 149, column 2 to line 153, column 43)
#> Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 1
#> Chain 2 Iteration:   1 / 600 [  0%]  (Warmup)
#> Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 2 Exception: Exception: ode_bdf_tol: ode parameters and data is inf, but must be finite! (in '/tmp/Rtmp1KA7pd/model-c47538a8e01f.stan', line 51, column 6 to line 55, column 8) (in '/tmp/Rtmp1KA7pd/model-c47538a8e01f.stan', line 149, column 2 to line 153, column 43)
#> Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 2
#> Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
#> Chain 1 Exception: Exception: CVode(cvodes_mem, t_final, nv_state_, &t_init, CV_NORMAL) failed with error flag -1:
#> Chain 1 The solver took mxstep internal steps but could not reach tout. (in '/tmp/Rtmp1KA7pd/model-c47538a8e01f.stan', line 51, column 6 to line 55, column 8) (in '/tmp/Rtmp1KA7pd/model-c47538a8e01f.stan', line 149, column 2 to line 153, column 43)
#> Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
#> Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
#> Chain 1
#> Chain 2 Iteration: 100 / 600 [ 16%]  (Warmup) 
#> Chain 1 Iteration: 100 / 600 [ 16%]  (Warmup) 
#> Chain 2 Iteration: 200 / 600 [ 33%]  (Warmup) 
#> Chain 1 Iteration: 200 / 600 [ 33%]  (Warmup) 
#> Chain 2 Iteration: 300 / 600 [ 50%]  (Warmup) 
#> Chain 2 Iteration: 301 / 600 [ 50%]  (Sampling) 
#> Chain 1 Iteration: 300 / 600 [ 50%]  (Warmup) 
#> Chain 1 Iteration: 301 / 600 [ 50%]  (Sampling) 
#> Chain 2 Iteration: 400 / 600 [ 66%]  (Sampling) 
#> Chain 2 Iteration: 500 / 600 [ 83%]  (Sampling) 
#> Chain 1 Iteration: 400 / 600 [ 66%]  (Sampling) 
#> Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
#> Chain 2 finished in 64.2 seconds.
#> Chain 1 Iteration: 500 / 600 [ 83%]  (Sampling) 
#> Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
#> Chain 1 finished in 90.9 seconds.
#> 
#> Both chains finished successfully.
#> Mean chain execution time: 77.6 seconds.
#> Total execution time: 91.0 seconds.
#> Warning: 20 of 600 (3.0%) transitions ended with a divergence.
#> See https://mc-stan.org/misc/warnings for details.

bdeb_diagnose(fit_hier)
#> 
#> ── BDEB Diagnostics (hierarchical) ──
#> 
#> ✖ Divergent transitions: 20
#> ℹ Consider: increase adapt_delta, reparameterise, or tighten priors.
#> ✔ Treedepth OK.
#> ✔ E-BFMI OK (all > 0.3).
#> 
#> ── Parameter Summary
#> ✖ R-hat > 1.01 for: sigma_log_p_Am, z_log_p_Am[2], z_log_p_Am[4], kappa, v, L0[4], p_Am_ind[2], p_Am_ind[5]
#> ! Low bulk ESS (<400) for: mu_log_p_Am, sigma_log_p_Am, z_log_p_Am[1], z_log_p_Am[2], z_log_p_Am[5], p_M, kappa, v, E_G, L0[1], L0[2], L0[3], L0[4], sigma_L, p_Am_ind[2], p_Am_ind[3], p_Am_ind[4]
#>        variable     mean       sd   median      5%      95% rhat ess_bulk
#>     mu_log_p_Am   1.5556 9.04e-01   1.5269  0.0902   3.1094 1.01      391
#>  sigma_log_p_Am   0.5418 5.40e-01   0.3924  0.0339   1.5244 1.02      142
#>             p_M   1.9723 2.81e+00   1.1812  0.1339   5.9670 1.01      108
#>           kappa   0.5146 2.06e-01   0.5175  0.1760   0.8573 1.01      250
#>               v   0.2090 1.95e-01   0.1514  0.0498   0.5986 1.01      167
#>             E_G 192.7620 2.52e+02 112.8990 33.3285 624.7077 1.01      245
#>              E0   1.5791 2.15e+00   0.9718  0.1878   4.8291 1.00      640
#>         sigma_L   0.0163 1.65e-03   0.0161  0.0141   0.0191 1.00      320
#>  ess_tail
#>     395.0
#>     286.1
#>      35.9
#>     187.9
#>      76.0
#>     380.7
#>     382.5
#>     245.7
#> ℹ 15 latent-state rows hidden; use `print(x, full = TRUE)` or `summary(x)$table` to see all.
summary(fit_hier, pars = c("mu_log_p_Am", "sigma_log_p_Am", "p_M", "kappa"))
#> # A tibble: 4 × 9
#>   variable        mean    sd median   `5%` `95%`  rhat ess_bulk ess_tail
#>   <chr>          <dbl> <dbl>  <dbl>  <dbl> <dbl> <dbl>    <dbl>    <dbl>
#> 1 mu_log_p_Am    1.56  0.904  1.53  0.0902 3.11   1.01     391.    395. 
#> 2 sigma_log_p_Am 0.542 0.540  0.392 0.0339 1.52   1.02     142.    286. 
#> 3 p_M            1.97  2.81   1.18  0.134  5.97   1.01     108.     35.9
#> 4 kappa          0.515 0.206  0.517 0.176  0.857  1.01     250.    188.

Example: DEBtox Model

data(debtox_growth)

# Concentration mapping
conc_map <- setNames(
  c(0, 20, 80, 200),
  c("1", "2", "3", "4")
)

dat_tox <- bdeb_data(
  growth = debtox_growth,
  concentration = conc_map,
  f_food = 1.0
)

mod_tox <- bdeb_tox(dat_tox, stress = "assimilation")
#> Warning: ! 4 concentration group(s) contain multiple individuals.
#> ℹ DEBtox fits one ODE per concentration group, not per individual.
#> ℹ Aggregating to group means per time point. For individual-level TKTD, a
#>   hierarchical DEBtox extension is needed (not yet implemented).

For this demo only, we use variational inference (ADVI) which gives an approximation of the posterior in seconds rather than minutes. The replication archive uses full HMC (NUTS) which is the publication-grade method.

fit_tox <- bdeb_fit(mod_tox, algorithm = "variational",
                    seed = 123, refresh = 0)
#> ℹ Compiling Stan model: 'bdeb_debtox'
#> ℹ Running variational inference (ADVI, mean-field; approximation, NOT exact MCMC).
#> ------------------------------------------------------------ 
#> EXPERIMENTAL ALGORITHM: 
#>   This procedure has not been thoroughly tested and may be unstable 
#>   or buggy. The interface is subject to change. 
#> ------------------------------------------------------------ 
#> Rejecting initial value: 
#>   Error evaluating the log probability at the initial value. 
#> Exception: Exception: CVode(cvodes_mem, t_final, nv_state_, &t_init, CV_NORMAL) failed with error flag -4:  
#> Convergence test failures occurred too many times during one internal time step or minimum step size was reached. (in '/tmp/RtmpyWTkBO/model-6de365c93a1a.stan', line 72, column 6 to line 78, column 8) (in '/tmp/RtmpyWTkBO/model-6de365c93a1a.stan', line 206, column 2 to line 213, column 59) 
#> Gradient evaluation took 0.005141 seconds 
#> 1000 transitions using 10 leapfrog steps per transition would take 51.41 seconds. 
#> Adjust your expectations accordingly! 
#> Begin eta adaptation. 
#> Iteration:   1 / 250 [  0%]  (Adaptation) 
#> Iteration:  50 / 250 [ 20%]  (Adaptation) 
#> Iteration: 100 / 250 [ 40%]  (Adaptation) 
#> Iteration: 150 / 250 [ 60%]  (Adaptation) 
#> Iteration: 200 / 250 [ 80%]  (Adaptation) 
#> Success! Found best value [eta = 1] earlier than expected. 
#> Begin stochastic gradient ascent. 
#>   iter             ELBO   delta_ELBO_mean   delta_ELBO_med   notes  
#>    100           64.232             1.000            1.000 
#>    200           77.982             0.588            1.000 
#>    300           77.786             0.393            0.176 
#>    400           77.399             0.296            0.176 
#>    500           79.454             0.242            0.026 
#>    600           72.759             0.217            0.092 
#>    700           74.871             0.190            0.028 
#>    800           78.083             0.171            0.041 
#>    900           77.548             0.153            0.028 
#>   1000           78.590             0.139            0.028 
#>   1100           78.587             0.039            0.026 
#>   1200           79.835             0.023            0.016 
#>   1300           78.441             0.025            0.018 
#>   1400           78.605             0.024            0.018 
#>   1500           74.630             0.027            0.018 
#>   1600           78.496             0.023            0.018 
#>   1700           77.040             0.022            0.018 
#>   1800           79.432             0.021            0.018 
#>   1900           76.063             0.024            0.019 
#>   2000           77.425             0.025            0.019 
#>   2100           79.030             0.027            0.020 
#>   2200           76.741             0.028            0.030 
#>   2300           78.230             0.028            0.030 
#>   2400           77.773             0.029            0.030 
#>   2500           78.263             0.024            0.020 
#>   2600           78.778             0.020            0.019 
#>   2700           79.487             0.019            0.019 
#>   2800           78.240             0.017            0.018 
#>   2900           79.690             0.015            0.018 
#>   3000           78.667             0.014            0.016 
#>   3100           79.524             0.013            0.013 
#>   3200           78.824             0.011            0.011 
#>   3300           78.383             0.010            0.009   MEDIAN ELBO CONVERGED 
#> Drawing a sample of size 1000 from the approximate posterior...  
#> COMPLETED. 
#> Finished in  11.1 seconds.

bdeb_ec50(fit_tox)
#> 
#> ── DEBtox Effect Concentrations
#>  parameter  mean median   sd  lower upper
#>       EC50 28.01   6.40 93.3 0.4941 127.3
#>        NEC  8.59   1.15 28.9 0.0338  39.1

plot_dose_response(fit_tox, n_draws = 20, n_conc = 25, dt = 1.0)

Prior Specification

BayesianDEB follows the prior recommendations of Kooijman (2010) and the AmP collection (Marques et al., 2018):

Positive rate parameters (p_Am, p_M, v, E_G): Log-normal priors
Allocation fraction (kappa): Beta prior
Scale parameters (sigma_L): Half-normal priors
Hierarchical SDs: Exponential priors

# View defaults
prior_default("individual")
#> $p_Am
#> → lognormal(mu = 1, sigma = 1)
#> 
#> $p_M
#> → lognormal(mu = -1, sigma = 1)
#> 
#> $kappa
#> → beta(a = 2, b = 2)
#> 
#> $v
#> → lognormal(mu = -1.5, sigma = 1)
#> 
#> $E_G
#> → lognormal(mu = 6, sigma = 1)
#> 
#> $E0
#> → lognormal(mu = 0, sigma = 1)
#> 
#> $L0
#> → lognormal(mu = -2, sigma = 1)
#> 
#> $sigma_L
#> → halfnormal(sigma = 0.1)

# Override specific priors
my_priors <- list(
  p_Am  = prior_lognormal(mu = 2.0, sigma = 0.3),
  kappa = prior_beta(a = 5, b = 2)
)

Observation Models

The default likelihood is Gaussian for growth and negative binomial for reproduction. Switch via observation:

# Robust to outliers
mod <- bdeb_model(dat, type = "individual",
  observation = list(growth = obs_student_t(nu = 5)))

# Multiplicative error
mod <- bdeb_model(dat, type = "individual",
  observation = list(growth = obs_lognormal()))