---
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}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = FALSE
)
```

## 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}
# 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")
```

## Example: Individual Growth Model

### 1. Prepare Data

```{r data}
library(BayesianDEB)
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 = 4, iter_sampling = 1000, seed = 123)
fit
```

### 4. Diagnostics

```{r diagnose}
bdeb_diagnose(fit)
plot(fit, type = "trace")
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)

```{r hierarchical}
dat_all <- bdeb_data(growth = eisenia_growth, f_food = 1.0)

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

fit_hier <- bdeb_fit(mod_hier, chains = 4, adapt_delta = 0.9,
                     threads_per_chain = 4)  # within-chain parallelism

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

## Example: DEBtox Model

```{r debtox}
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")
fit_tox <- bdeb_fit(mod_tox, chains = 4, adapt_delta = 0.95)

# EC50 and NEC
bdeb_ec50(fit_tox)

# Dose-response plot
plot_dose_response(fit_tox)
```

## 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.