Dynamic Energy Budget (DEB) theory provides a mechanistic, thermodynamically consistent framework for describing how organisms acquire and utilise energy for maintenance, growth, development, and reproduction (Kooijman 2010). Classical DEB calibration relies on deterministic optimisation, yielding a single best-fit parameter vector without formal uncertainty quantification.
BayesianDEB embeds the DEB ordinary differential equation (ODE) system within a Bayesian state-space model, using Stan’s Hamiltonian Monte Carlo (HMC) sampler (Carpenter et al. 2017) to explore the full joint posterior distribution of all parameters. This approach naturally provides:
This vignette demonstrates the complete workflow on two standard ecotoxicological test organisms:
library(BayesianDEB)
library(ggplot2)
library(posterior) # for summarise_draws()BayesianDEB requires cmdstanr and a working CmdStan installation for model fitting. Data preparation, prior specification, and utility functions work without Stan.
check_cmdstanr() # informative error if missingThe eisenia_growth dataset contains simulated weekly
length measurements for 21 Eisenia fetida individuals over 84
days (12 weeks). The simulation used the standard 2-state DEB model
(reserve \(E\), structure \(V\)) with parameters representative of
E. fetida from the AmP collection: \(\{p_{Am}\} = 5.0\) J d\(^{-1}\) cm\(^{-2}\), \([p_M]
= 0.5\) J d\(^{-1}\) cm\(^{-3}\), \(\kappa
= 0.75\), \(v = 0.2\) cm d\(^{-1}\), \([E_G]
= 400\) J cm\(^{-3}\).
Individual variation in \(\{p_{Am}\}\)
(CV \(\approx\) 10%) and Gaussian
observation error (\(\sigma_L = 0.015\)
cm) were added.
data(eisenia_growth)
# Structure: 273 obs, 3 variables (id, time, length)
str(eisenia_growth)
length(unique(eisenia_growth$id)) # 21 individuals
length(unique(eisenia_growth$time)) # 13 time points (days 0–84)ggplot(eisenia_growth, aes(time, length, group = id)) +
geom_line(alpha = 0.3, colour = "steelblue") +
geom_point(size = 0.8, alpha = 0.4) +
theme_bw(base_size = 12) +
labs(x = "Time (days)", y = expression(paste("Structural length ", L, " (cm)")),
title = "Eisenia fetida: 21 individuals, 12 weeks")Key features visible in the data:
We start with a single individual to validate the approach.
df1 <- eisenia_growth[eisenia_growth$id == 5, ]
dat1 <- bdeb_data(growth = df1, f_food = 1.0)
dat1The f_food = 1.0 argument specifies ad libitum feeding
(\(f = 1\), the ratio of actual to
maximum ingestion rate; Kooijman 2010, Eq. 2.3).
The individual model tracks two state variables: reserve energy \(E\) (J) and structural volume \(V\) (cm\(^3\)), governed by:
\[ \frac{dE}{dt} = f\{p_{Am}\}L^2 - \frac{EvL}{E + [E_G]V}, \qquad \frac{dV}{dt} = \frac{\kappa \dot{p}_C - [p_M]V}{[E_G]} \]
where \(L = V^{1/3}\) is structural length. Observed lengths are assumed to follow \(L_\text{obs} \sim \mathcal{N}(\hat{L}, \sigma_L)\).
We set biologically informed priors based on published AmP values for earthworms:
| Parameter | Prior | Rationale |
|---|---|---|
| \(\{p_{Am}\}\) | LogNormal(1.5, 0.5) | Median \(e^{1.5} \approx 4.5\); AmP range 3–8 |
| \([p_M]\) | LogNormal(−1.0, 0.5) | Median \(e^{-1} \approx 0.37\); typical 0.15–0.8 |
| \(\kappa\) | Beta(3, 2) | Mode 0.67; earthworms allocate ~60–80% to soma |
| \(v\) | LogNormal(−1.5, 0.5) | Median \(e^{-1.5} \approx 0.22\) cm/d |
| \([E_G]\) | LogNormal(6.0, 0.5) | Median \(e^6 \approx 403\); typical 200–800 |
| \(\sigma_L\) | HalfNormal(0.05) | Measurement precision ~0.01–0.03 cm |
mod1 <- bdeb_model(dat1, 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),
v = prior_lognormal(mu = -1.5, sigma = 0.5),
E_G = prior_lognormal(mu = 6.0, sigma = 0.5),
sigma_L = prior_halfnormal(sigma = 0.05)
))
mod1Unspecified priors (E0, L0) are filled
automatically from prior_default("individual").
The model is compiled to C++ and sampled using the No-U-Turn Sampler (NUTS; Hoffman and Gelman 2014) with the stiff BDF ODE solver.
fit1 <- bdeb_fit(mod1,
chains = 4,
iter_warmup = 1000,
iter_sampling = 2000,
adapt_delta = 0.9,
seed = 42
)
fit1Diagnostics follow the recommendations of Vehtari et al. (2021):
diag1 <- bdeb_diagnose(fit1)What to check:
adapt_delta (e.g., 0.95).plot(fit1, type = "trace",
pars = c("p_Am", "p_M", "kappa", "sigma_L"))plot(fit1, type = "pairs",
pars = c("p_Am", "p_M", "kappa", "E_G"))bdeb_summary(fit1,
pars = c("p_Am", "p_M", "kappa", "v", "E_G", "sigma_L"),
prob = 0.95)Posterior predictive checks (PPCs) compare data replicated from the fitted model (\(L^\text{rep}\)) with the observed data (\(L^\text{obs}\)). If the model fits well, the observed points should fall within the envelope of replicated trajectories (Gelman et al. 2013, Ch. 6).
ppc1 <- bdeb_ppc(fit1, type = "growth")
plot(ppc1, n_draws = 200)plot(fit1, type = "trajectory", n_draws = 200)BayesianDEB computes biologically meaningful quantities directly from the posterior, automatically propagating uncertainty:
| Quantity | Formula | Interpretation |
|---|---|---|
| \(L_m\) | \(\kappa \{p_{Am}\} / [p_M]\) | Maximum structural length (\(f = 1\)) |
| \(L_\infty\) | \(f \cdot L_m\) | Ultimate structural length at food \(f\) |
| \(k_M\) | \([p_M] / [E_G]\) | Somatic maintenance rate constant |
| \(g\) | \([E_G] \, v / (\kappa \{p_{Am}\})\) | Energy investment ratio |
| \(\dot{r}_B\) | \(k_M \, g \,/\, 3(f + g)\) | von Bertalanffy growth rate (Eq. 3.23) |
Note: all lengths are structural (\(L = V^{1/3}\)), not physical. Physical length \(L_w = L / \delta_M\) where \(\delta_M\) is the species-specific shape coefficient (not estimated by this package).
der1 <- bdeb_derived(fit1,
quantities = c("L_m", "L_inf", "k_M", "g", "growth_rate"), f = 1.0)
summarise_draws(der1,
"mean", "sd",
"q2.5" = ~quantile(.x, 0.025),
"q97.5" = ~quantile(.x, 0.975))How does ultimate size change if food availability drops to 70%? Since \(L_\infty \propto f\), we can compute this directly:
d_f10 <- bdeb_derived(fit1, quantities = "L_inf", f = 1.0)
d_f07 <- bdeb_derived(fit1, quantities = "L_inf", f = 0.7)
df_compare <- data.frame(
L_inf = c(d_f10$L_inf, d_f07$L_inf),
food = rep(c("f = 1.0", "f = 0.7"), each = nrow(d_f10))
)
ggplot(df_compare, aes(x = L_inf, fill = food)) +
geom_density(alpha = 0.4) +
theme_bw(base_size = 12) +
labs(x = expression(L[infinity] ~ "(cm)"),
y = "Posterior density",
fill = "Food level")When multiple individuals are available, a hierarchical model is preferred. It estimates population-level distributions for parameters that vary across individuals, while sharing parameters that are species-level constants.
The hierarchical model places a lognormal random effect on the assimilation rate:
\[ \log\{p_{Am}\}_j = \mu_{\log p_{Am}} + \sigma_{\log p_{Am}} \cdot z_j, \qquad z_j \sim \mathcal{N}(0, 1) \]
This non-centred parameterisation (Betancourt and Girolami 2015) avoids the pathological funnel geometry that arises when \(\sigma\) is small. The parameters \([p_M]\), \(\kappa\), \(v\), \([E_G]\) are shared across all individuals.
dat_all <- bdeb_data(growth = eisenia_growth, f_food = 1.0)
dat_all # 21 individuals, 273 observationsmod_h <- bdeb_model(dat_all, type = "hierarchical",
priors = list(
mu_log_p_Am = prior_normal(mu = 1.5, sigma = 0.5),
sigma_log_p_Am = prior_exponential(rate = 2),
p_M = prior_lognormal(mu = -1.0, sigma = 0.5),
kappa = prior_beta(a = 3, b = 2),
v = prior_lognormal(mu = -1.5, sigma = 0.5),
E_G = prior_lognormal(mu = 6.0, sigma = 0.5),
sigma_L = prior_halfnormal(sigma = 0.05)
))The prior on \(\sigma_{\log p_{Am}}\) uses an Exponential(2) distribution, following the guidance of Gelman (2006) for hierarchical variance parameters. This places most prior mass near zero while allowing substantial variation if the data support it.
With 21 individuals, each requiring an independent ODE solve per
iteration, the hierarchical model benefits from within-chain
parallelism via Stan’s reduce_sum. Setting
threads_per_chain = 4 distributes the 21 ODE solves across
4 threads per chain.
fit_h <- bdeb_fit(mod_h,
chains = 4,
iter_warmup = 1000,
iter_sampling = 2000,
adapt_delta = 0.95,
max_treedepth = 12,
threads_per_chain = 4,
seed = 123
)bdeb_diagnose(fit_h)plot(fit_h, type = "trace",
pars = c("mu_log_p_Am", "sigma_log_p_Am"))plot(fit_h, type = "posterior",
pars = c("mu_log_p_Am", "sigma_log_p_Am", "p_M", "kappa"))bdeb_summary(fit_h,
pars = c("mu_log_p_Am", "sigma_log_p_Am",
"p_M", "kappa", "v", "E_G", "sigma_L"),
prob = 0.95)A key feature of hierarchical models is shrinkage: individuals with sparse or noisy data are pulled toward the population mean. This is partial pooling — a principled compromise between complete pooling (ignoring individual variation) and no pooling (fitting each individual independently).
ind_summary <- bdeb_summary(fit_h,
pars = paste0("p_Am_ind[", 1:21, "]"),
prob = 0.90)
pop_summary <- bdeb_summary(fit_h, pars = "mu_log_p_Am")
pop_mean_pAm <- exp(as.data.frame(pop_summary)$mean)
ind_df <- as.data.frame(ind_summary)
ind_df$individual <- 1:21
ggplot(ind_df, aes(x = individual, y = mean)) +
geom_pointrange(aes(ymin = `5%`, ymax = `95%`),
colour = "steelblue", size = 0.4) +
geom_hline(yintercept = pop_mean_pAm, linetype = "dashed",
colour = "red", linewidth = 0.8) +
theme_bw(base_size = 12) +
labs(x = "Individual", y = expression({p[Am]} ~ "(J/d/cm"^2*")"),
title = "Individual assimilation rates with 90% CI")The generated quantities block draws
p_Am_new from the population distribution — useful for
predicting the performance of an unobserved individual from the same
population:
bdeb_summary(fit_h, pars = "p_Am_new", prob = 0.95)s_ind <- bdeb_summary(fit1, pars = c("p_Am", "p_M", "kappa"), prob = 0.90)
s_hier <- bdeb_summary(fit_h, pars = c("p_M", "kappa"), prob = 0.90)
cat("=== Individual model (id = 5, n = 1) ===\n")
print(as.data.frame(s_ind), digits = 3, row.names = FALSE)
cat("\n=== Hierarchical model (n = 21) ===\n")
print(as.data.frame(s_hier), digits = 3, row.names = FALSE)The hierarchical model yields narrower credible intervals for shared parameters (\([p_M]\), \(\kappa\)) because it pools information across 21 individuals.
The DEBtox framework (Jager et al. 2006; Jager and Zimmer 2012) extends DEB with toxicokinetic-toxicodynamic (TKTD) components. A scaled internal damage variable \(D_w\) tracks the toxicant’s effect:
\[ \frac{dD_w}{dt} = k_d\bigl(\max(C_w - z_w, 0) - D_w\bigr) \]
where \(k_d\) is the damage recovery rate, \(C_w\) is the external concentration, and \(z_w\) is the no-effect concentration (NEC). The damage causes a stress factor \(s = b_w \cdot D_w\) that reduces assimilation:
\[ \dot{p}_A = f\{p_{Am}\}L^2 \cdot \max(1 - s, \; 0) \]
At steady state (\(D_w = C_w - z_w\) for \(C_w > z_w\)), the EC50 for 50% assimilation reduction is:
\[ \text{EC}_{50} = z_w + \frac{0.5}{b_w} \]
The debtox_growth dataset simulates growth under 4
toxicant concentrations (0, 20, 80, 200 arbitrary units), 10 individuals
per group, measured weekly over 6 weeks. True parameters: NEC = 15,
\(b_w = 0.003\).
data(debtox_growth)
ggplot(debtox_growth,
aes(time, length, colour = factor(concentration), group = id)) +
geom_line(alpha = 0.3) +
geom_point(size = 0.8, alpha = 0.4) +
facet_wrap(~concentration, labeller = label_both) +
theme_bw(base_size = 11) +
scale_colour_brewer(palette = "RdYlBu", direction = -1) +
labs(x = "Time (days)", y = "Structural length (cm)",
colour = "Concentration") +
theme(legend.position = "none")conc_levels <- unique(debtox_growth$concentration)
conc_map <- setNames(conc_levels, as.character(conc_levels))
dat_tox <- bdeb_data(
growth = debtox_growth,
concentration = conc_map,
f_food = 1.0
)
dat_toxmod_tox <- bdeb_tox(dat_tox, stress = "assimilation",
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),
v = prior_lognormal(mu = -1.5, sigma = 0.5),
E_G = prior_lognormal(mu = 6.0, sigma = 0.5),
sigma_L = prior_halfnormal(sigma = 0.05),
k_d = prior_lognormal(mu = -1.0, sigma = 1.0),
z_w = prior_lognormal(mu = 2.5, sigma = 1.0),
b_w = prior_lognormal(mu = -5.0, sigma = 2.0)
))
mod_toxPrior rationale for toxicological parameters:
| Parameter | Prior | Median | 95% prior range | Rationale |
|---|---|---|---|---|
| \(k_d\) | LogNormal(−1, 1) | 0.37 d\(^{-1}\) | 0.05–2.7 | Damage recovery: hours to days |
| \(z_w\) | LogNormal(2.5, 1) | 12.2 | 1.6–89 | NEC within tested range 0–200 |
| \(b_w\) | LogNormal(−5, 2) | 0.0067 | 0.00009–0.5 | Weakly informative on effect intensity |
fit_tox <- bdeb_fit(mod_tox,
chains = 4,
iter_warmup = 1000,
iter_sampling = 2000,
adapt_delta = 0.95,
max_treedepth = 12,
threads_per_chain = 2,
seed = 77
)bdeb_diagnose(fit_tox)plot(fit_tox, type = "trace", pars = c("k_d", "z_w", "b_w"))plot(fit_tox, type = "posterior", pars = c("k_d", "z_w", "b_w"))plot(fit_tox, type = "pairs", pars = c("z_w", "b_w", "k_d"))bdeb_summary(fit_tox,
pars = c("p_Am", "p_M", "kappa", "v", "E_G",
"k_d", "z_w", "b_w", "sigma_L"),
prob = 0.95)The EC50 is computed analytically in the Stan
generated quantities block, giving the full posterior
distribution without post-hoc root-finding.
ec <- bdeb_ec50(fit_tox, prob = 0.95)
print(ec$summary, digits = 3)
hist(ec$draws, breaks = 50, col = "steelblue", border = "white",
main = expression("Posterior distribution of EC"[50]),
xlab = "Concentration", freq = FALSE)
abline(v = ec$summary$median[1], col = "red", lwd = 2, lty = 2)
legend("topright", "Posterior median",
col = "red", lty = 2, lwd = 2, bty = "n")Interpretation for risk assessment:
plot_dose_response(fit_tox, n_draws = 200)A Bayesian analysis should always report the sensitivity of key conclusions to prior choices. We refit with a tighter prior on \(z_w\):
mod_tox2 <- bdeb_tox(dat_tox, stress = "assimilation",
priors = list(
z_w = prior_lognormal(mu = 3.0, sigma = 0.3), # tighter
b_w = prior_lognormal(mu = -5.0, sigma = 2.0)
))
fit_tox2 <- bdeb_fit(mod_tox2, chains = 4, adapt_delta = 0.95,
threads_per_chain = 2, seed = 78)
cat("=== Original: z_w ~ LogNormal(2.5, 1.0) ===\n")
bdeb_summary(fit_tox, pars = c("z_w", "b_w"), prob = 0.95)
cat("\n=== Tighter: z_w ~ LogNormal(3.0, 0.3) ===\n")
bdeb_summary(fit_tox2, pars = c("z_w", "b_w"), prob = 0.95)If the posteriors agree despite different priors, the data are informative and the inference is robust. If they diverge, the parameter is prior-dominated and should be reported as weakly identified.
BayesianDEB works with structural length \(L = V^{1/3}\) (cube root of structural volume), not the physical body length \(L_w\) that you measure with a ruler. The two are related by the species-specific shape coefficient \(\delta_M\):
\[L = \delta_M \times L_w\]
| Species | \(\delta_M\) | Source |
|---|---|---|
| Eisenia fetida | 0.24 | AmP |
| Folsomia candida | 0.19 | AmP |
| Daphnia magna | 0.37 | AmP |
Before fitting, convert your measured lengths:
# Example: measured body lengths in mm for E. fetida
L_physical_mm <- c(12, 18, 25, 30)
delta_M <- 0.24
# Convert to structural length in cm
L_structural_cm <- delta_M * L_physical_mm / 10
L_structural_cm
# [1] 0.288 0.432 0.600 0.720If you pass physical lengths directly, the estimated DEB parameters will absorb the shape coefficient and will not be comparable with AmP values. BayesianDEB warns if maximum length exceeds 10 cm, which is unusually large for structural length.
Before fitting, it is good practice to verify that priors produce biologically plausible predictions (Gabry et al., 2019). Sample directly from the prior distributions and compute derived quantities:
set.seed(42)
n_sim <- 4000
# Sample from priors
p_Am_sim <- rlnorm(n_sim, 1.5, 0.5)
p_M_sim <- rlnorm(n_sim, -1.0, 0.5)
kappa_sim <- rbeta(n_sim, 3, 2)
v_sim <- rlnorm(n_sim, -1.5, 0.5)
E_G_sim <- rlnorm(n_sim, 6.0, 0.5)
# Prior predictive for L_inf
L_inf_prior <- kappa_sim * p_Am_sim / p_M_sim
hist(L_inf_prior, breaks = 50, col = "steelblue", border = "white",
main = "Prior predictive: ultimate structural length",
xlab = expression(L[infinity] ~ "(cm)"), xlim = c(0, 50))
# Should cover plausible range for earthworms (~2-20 cm structural)If the prior predictive distribution covers unreasonable values (e.g., \(L_\infty > 100\) cm for an earthworm), tighten the priors.
By default, growth observations use a Gaussian likelihood and
reproduction uses negative binomial. You can switch the observation
model via the observation argument — the Stan likelihood is
controlled by integer flags, so no recompilation is
needed:
# Robust to outliers: Student-t with 5 df
mod_robust <- bdeb_model(dat1, type = "individual",
observation = list(growth = obs_student_t(nu = 5)))
# Multiplicative error (constant CV)
mod_logn <- bdeb_model(dat1, type = "individual",
observation = list(growth = obs_lognormal()))
# For reproduction: Poisson instead of NegBin
# (appropriate when overdispersion is negligible)
mod_pois <- bdeb_model(dat_gr, type = "growth_repro",
observation = list(growth = obs_normal(),
reproduction = obs_poisson()))Available observation families:
| Endpoint | Family | Function | When to use |
|---|---|---|---|
| Growth | Gaussian | obs_normal() |
Default; additive error |
| Growth | Log-normal | obs_lognormal() |
Multiplicative error (constant CV) |
| Growth | Student-t | obs_student_t(nu) |
Outlier-robust |
| Reproduction | Neg. binomial | obs_negbinom() |
Default; overdispersed counts |
| Reproduction | Poisson | obs_poisson() |
Equidispersed counts |
DEB rate parameters scale with temperature via the Arrhenius relationship (Kooijman 2010, Eq. 1.2):
\[ c_T = \exp\!\left(\frac{T_A}{T_\text{ref}} - \frac{T_A}{T}\right) \]
# Experiment at 22 C, reference 20 C, typical T_A for ectotherms
cT <- arrhenius(temp = 273.15 + 22, T_ref = 273.15 + 20, T_A = 8000)
cat("Temperature correction factor:", round(cT, 3), "\n")
# Rate at reference temperature: p_Am_ref = p_Am_obs / cTInspect the energetics at a specific state:
fl <- deb_fluxes(E = 10, V = 0.5, f = 1.0,
p_Am = 5, p_M = 0.5, kappa = 0.75,
v = 0.2, E_G = 400)
cat(sprintf("Assimilation (p_A): %.3f J/d\n", fl$p_A))
cat(sprintf("Mobilisation (p_C): %.3f J/d\n", fl$p_C))
cat(sprintf("Maintenance (p_M): %.3f J/d\n", fl$p_M))
cat(sprintf("Growth (p_G): %.3f J/d\n", fl$p_G))
cat(sprintf("Struct. length (L) : %.3f cm\n", fl$L))
cat(sprintf("Scaled reserve (e) : %.3f\n", fl$e))Many protocols report cumulative offspring. The
repro_to_intervals() function converts these to the
interval format required by BayesianDEB:
cumul <- data.frame(
id = rep(1, 5),
time = c(0, 7, 14, 21, 28),
cumulative = c(0, 10, 30, 60, 100)
)
repro_to_intervals(cumul)
# id t_start t_end count
# 1 1 0 7 10
# 2 1 7 14 20
# 3 1 14 21 30
# 4 1 21 28 40The eisenia_growth dataset was simulated using DEB
parameters from the Add-my-Pet (AmP) collection entry for Eisenia
fetida (Marques et al., 2018). The table below compares the
simulation truth with published AmP estimates and the expected posterior
recovery from BayesianDEB:
| Parameter | Symbol | True (simulation) | AmP estimate | Units |
|---|---|---|---|---|
| Assimilation rate | \(\{p_{Am}\}\) | 5.0 | 3.9–6.2 | J d\(^{-1}\) cm\(^{-2}\) |
| Maintenance rate | \([p_M]\) | 0.5 | 0.3–0.8 | J d\(^{-1}\) cm\(^{-3}\) |
| Allocation fraction | \(\kappa\) | 0.75 | 0.6–0.85 | — |
| Energy conductance | \(v\) | 0.2 | 0.1–0.3 | cm d\(^{-1}\) |
| Cost of structure | \([E_G]\) | 400 | 200–600 | J cm\(^{-3}\) |
| Max. structural length | \(L_m\) | 7.5 | 5–12 | cm |
The simulation truth falls within the published AmP ranges for all parameters. When BayesianDEB is fitted to these data (see Section @ref(individual)), the posterior medians should recover values close to the simulation truth, providing a closed-loop validation: known parameters → simulated data → Bayesian recovery → comparison with truth.
This is not a substitute for fitting real experimental data, but it demonstrates that:
For real-data applications, we recommend comparing estimated parameters with the AmP entry for the species of interest as a sanity check.
This vignette demonstrated three analysis types:
| Analysis | Model type | Individuals | Key output |
|---|---|---|---|
| Single growth | "individual" |
1 | DEB posteriors, PPC, \(L_\infty\) |
| Population growth | "hierarchical" |
21 | \(\mu/\sigma\) of \(\{p_{Am}\}\), shrinkage, prediction for new individual |
| Toxicant effect | "debtox" |
40 (4 groups) | EC50, NEC with full uncertainty |
What the Bayesian framework provides:
reduce_sum
accelerates hierarchical and DEBtox models on multi-core machines.Betancourt, M. and Girolami, M. (2015). Hamiltonian Monte Carlo for hierarchical models. In: Upadhyay, S.K. et al. (eds) Current Trends in Bayesian Methodology with Applications. CRC Press, pp. 79–101.
Carpenter, B., Gelman, A., Hoffman, M.D. et al. (2017). Stan: A probabilistic programming language. Journal of Statistical Software, 76(1), 1–32. doi: 10.18637/jss.v076.i01
ECHA (2017). Guidance on Information Requirements and Chemical Safety Assessment, Chapter R.10: Characterisation of dose [concentration]- response for environment. European Chemicals Agency.
Gelman, A. (2006). Prior distributions for variance parameters in hierarchical models. Bayesian Analysis, 1(3), 515–534. doi: 10.1214/06-BA117A
Gelman, A. and Hill, J. (2006). Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press.
Gelman, A., Carlin, J.B., Stern, H.S. et al. (2013). Bayesian Data Analysis. 3rd edition. Chapman & Hall/CRC.
Hoffman, M.D. and Gelman, A. (2014). The No-U-Turn Sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15(47), 1593–1623.
Jager, T., Heugens, E.H.W. and Kooijman, S.A.L.M. (2006). Making sense of ecotoxicological test results: towards application of process-based models. Ecotoxicology, 15(3), 305–314. doi: 10.1007/s10646-006-0060-x
Jager, T. and Zimmer, E.I. (2012). Simplified Dynamic Energy Budget model for analysing ecotoxicity data. Ecological Modelling, 225, 74–81. doi: 10.1016/j.ecolmodel.2011.11.012
Kooijman, S.A.L.M. (2010). Dynamic Energy Budget Theory for Metabolic Organisation. 3rd edition. Cambridge University Press. doi: 10.1017/CBO9780511805400
Marques, G.M., Augustine, S., Lika, K. et al. (2018). The AmP project: comparing species on the basis of dynamic energy budget parameters. PLOS Computational Biology, 14(5), e1006100. doi: 10.1371/journal.pcbi.1006100
Vehtari, A., Gelman, A., Simpson, D. et al. (2021). Rank-normalization, folding, and localization: an improved \(\hat{R}\) for assessing convergence of MCMC. Bayesian Analysis, 16(2), 667–718. doi: 10.1214/20-BA1221
sessionInfo()