## ----setup, include = FALSE--------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5, fig.align = "center", eval = FALSE ) ## ----prerequisites------------------------------------------------------------ # library(BayesianDEB) # library(ggplot2) # library(posterior) # for summarise_draws() ## ----check-stan--------------------------------------------------------------- # check_cmdstanr() # informative error if missing ## ----eisenia-explore---------------------------------------------------------- # 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) ## ----eisenia-plot, fig.cap="Growth trajectories of 21 *E. fetida* individuals. Structural length $L = V^{1/3}$ measured weekly over 12 weeks."---- # 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") ## ----ind-data----------------------------------------------------------------- # df1 <- eisenia_growth[eisenia_growth$id == 5, ] # dat1 <- bdeb_data(growth = df1, f_food = 1.0) # dat1 ## ----ind-model---------------------------------------------------------------- # 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) # )) # mod1 ## ----ind-fit------------------------------------------------------------------ # fit1 <- bdeb_fit(mod1, # chains = 4, # iter_warmup = 1000, # iter_sampling = 2000, # adapt_delta = 0.9, # seed = 42 # ) # fit1 ## ----ind-diag----------------------------------------------------------------- # diag1 <- bdeb_diagnose(fit1) ## ----ind-trace, fig.cap="MCMC trace plots for core DEB parameters. Well-mixed chains should appear as overlapping 'hairy caterpillars'."---- # plot(fit1, type = "trace", # pars = c("p_Am", "p_M", "kappa", "sigma_L")) ## ----ind-pairs, fig.cap="Bivariate posterior scatter. Strong correlation between $\\{p_{Am}\\}$ and $[p_M]$ is expected: both control ultimate size $L_\\infty = \\kappa \\{p_{Am}\\} / [p_M]$."---- # plot(fit1, type = "pairs", # pars = c("p_Am", "p_M", "kappa", "E_G")) ## ----ind-summary-------------------------------------------------------------- # bdeb_summary(fit1, # pars = c("p_Am", "p_M", "kappa", "v", "E_G", "sigma_L"), # prob = 0.95) ## ----ind-ppc, fig.cap="Posterior predictive check: grey lines are replicated growth trajectories, red points are observed data."---- # ppc1 <- bdeb_ppc(fit1, type = "growth") # plot(ppc1, n_draws = 200) ## ----ind-traj, fig.cap="Posterior predicted trajectories (blue) with observed data (black points). The spread reflects parameter uncertainty."---- # plot(fit1, type = "trajectory", n_draws = 200) ## ----ind-derived-------------------------------------------------------------- # 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)) ## ----ind-food, fig.cap="Posterior distributions of $L_\\infty$ at $f = 1.0$ (blue) and $f = 0.7$ (orange)."---- # 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") ## ----hier-data---------------------------------------------------------------- # dat_all <- bdeb_data(growth = eisenia_growth, f_food = 1.0) # dat_all # 21 individuals, 273 observations ## ----hier-model--------------------------------------------------------------- # mod_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) # )) ## ----hier-fit----------------------------------------------------------------- # 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 # ) ## ----hier-diag---------------------------------------------------------------- # bdeb_diagnose(fit_h) ## ----hier-trace, fig.cap="Trace plots for population-level hyperparameters $\\mu_{\\log p_{Am}}$ and $\\sigma_{\\log p_{Am}}$."---- # plot(fit_h, type = "trace", # pars = c("mu_log_p_Am", "sigma_log_p_Am")) ## ----hier-post, fig.cap="Marginal posterior densities for shared parameters."---- # plot(fit_h, type = "posterior", # pars = c("mu_log_p_Am", "sigma_log_p_Am", "p_M", "kappa")) ## ----hier-pop----------------------------------------------------------------- # 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) ## ----hier-shrinkage, fig.cap="Individual-level $\\{p_{Am}\\}$ estimates (points: posterior means; bars: 90% CI) compared to the population mean (dashed red line). Shrinkage toward the mean is visible for individuals with noisier data."---- # 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") ## ----hier-new----------------------------------------------------------------- # bdeb_summary(fit_h, pars = "p_Am_new", prob = 0.95) ## ----compare------------------------------------------------------------------ # 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) ## ----debtox-explore, fig.cap="Growth trajectories under 4 toxicant concentrations. Higher concentrations suppress growth through reduced assimilation."---- # 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") ## ----debtox-prep-------------------------------------------------------------- # 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_tox ## ----debtox-model------------------------------------------------------------- # mod_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_tox ## ----debtox-fit--------------------------------------------------------------- # 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 # ) ## ----debtox-diag-------------------------------------------------------------- # bdeb_diagnose(fit_tox) ## ----debtox-trace-tox, fig.cap="Trace plots for the three toxicological parameters. Good mixing is essential for reliable EC$_{50}$ and NEC estimates."---- # plot(fit_tox, type = "trace", pars = c("k_d", "z_w", "b_w")) ## ----debtox-post-tox, fig.cap="Marginal posterior densities for toxicological parameters."---- # plot(fit_tox, type = "posterior", pars = c("k_d", "z_w", "b_w")) ## ----debtox-pairs, fig.cap="Posterior pairs for toxicological parameters. A correlation between $z_w$ and $b_w$ is expected since both determine the shape of the dose-response curve."---- # plot(fit_tox, type = "pairs", pars = c("z_w", "b_w", "k_d")) ## ----debtox-summary----------------------------------------------------------- # 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) ## ----debtox-ec50, fig.cap="Posterior distribution of EC$_{50}$ (blue histogram) with the posterior median (red dashed line). The full distribution — not just a point estimate — is available for regulatory risk assessment."---- # 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") ## ----debtox-dr, fig.cap="Dose-response curve with posterior uncertainty bands (blue lines: individual posterior draws). The dashed horizontal line marks 50% effect; vertical dashed lines mark the NEC (green) and EC$_{50}$ (red)."---- # plot_dose_response(fit_tox, n_draws = 200) ## ----debtox-sens-------------------------------------------------------------- # 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) ## ----convert-length----------------------------------------------------------- # # 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.720 ## ----prior-pred--------------------------------------------------------------- # 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) ## ----obs-models--------------------------------------------------------------- # # 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())) ## ----arrhenius---------------------------------------------------------------- # # 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 / cT ## ----fluxes------------------------------------------------------------------- # 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)) ## ----repro-convert------------------------------------------------------------ # 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 40 ## ----session------------------------------------------------------------------ # sessionInfo()