---
title: "Getting Started with BayesianDEB"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Getting Started with BayesianDEB}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
HAS_CMDSTAN <- requireNamespace("cmdstanr", quietly = TRUE) &&
  isTRUE(nzchar(tryCatch(cmdstanr::cmdstan_path(),
                         error = function(e) "")))
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  fig.width  = 6,
  fig.height = 4,
  eval = HAS_CMDSTAN
)
library(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

```{r install, eval = FALSE}
# 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

```{r 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

```{r 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

```{r fit}
fit <- bdeb_fit(mod,
                chains = 2, iter_warmup = 300, iter_sampling = 300,
                seed = 123, refresh = 100)
fit
```

### 4. Diagnostics

```{r diagnose}
bdeb_diagnose(fit)
plot(fit, type = "trace")
```

```{r diagnose-pairs, eval = HAS_CMDSTAN && requireNamespace("gridExtra", quietly = TRUE)}
# `bayesplot::mcmc_pairs` requires gridExtra (a Suggests of bayesplot).
plot(fit, type = "pairs", pars = c("p_Am", "p_M", "kappa"))
```

### 5. Posterior Predictive Checks

```{r ppc}
ppc <- bdeb_ppc(fit, type = "growth")
plot(ppc)
```

### 6. Derived Quantities

```{r derived}
bdeb_derived(fit, quantities = c("L_inf", "growth_rate"))
```

### 7. Trajectory Plot

```{r trajectory}
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.

```{r hierarchical}
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

```{r debtox-data}
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.

```{r debtox-fit}
fit_tox <- bdeb_fit(mod_tox, algorithm = "variational",
                    seed = 123, refresh = 0)
```

```{r debtox-ec50}
bdeb_ec50(fit_tox)
```

```{r debtox-plot}
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

```{r 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`:

```{r obs}
# 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()))
```

## Further Reading

- Kooijman, S.A.L.M. (2010). *Dynamic Energy Budget Theory for Metabolic
  Organisation*. 3rd edition. Cambridge University Press.
- 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.
- Carpenter, B. et al. (2017). Stan: A probabilistic programming
  language. *Journal of Statistical Software*, 76(1).
- Gelman, A. et al. (2013). *Bayesian Data Analysis*. 3rd edition.
  Chapman & Hall/CRC.