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

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

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

4. Diagnostics

bdeb_diagnose(fit)
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"))

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)

bdeb_diagnose(fit_hier)
summary(fit_hier, pars = c("mu_log_p_Am", "sigma_log_p_Am", "p_M", "kappa"))

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

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)

bdeb_ec50(fit_tox)

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

# 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()))