---
title: "Growth estimation with linkages and priors"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{Growth estimation with linkages and priors}
  %\VignetteEncoding{UTF-8}
---
  
```{r, include = FALSE}
  knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5,
  eval = FALSE
  )
```
  
## Overview
  
This vignette demonstrates how to fit growth using Rceattle's built-in
von Bertalanffy and Richards growth functions, and how to use the
formula-driven `linkage_spec()` API to attach priors to growth
linkages.

We use the internal `whamGrowthData` dataset, which includes conditional
age-at-length data (`caal_data`) and the stock assessment data needed to
estimate growth parameters.

## Load data and inspect the growth dataset
  
```{r}
library(Rceattle)

data(whamGrowthData)
plot_data(whamGrowthData)

# The dataset includes CAAL data for the survey.
head(whamGrowthData$caal_data)

# Use logistic selectivity for all fleets in this example.
whamGrowthData$fleet_control$Selectivity <- "Logistic"
whamGrowthData$fleet_control$Selectivity_dimension <- "Length"
```
  
## Build a growth model with linkages
  
Rceattle supports growth linkages through `build_growth()` and
`linkage_spec()`. The linkage table can alter any von Bertalanffy growth
parameter such as `K`, `Linf`, or `L1` (natural-scale names; the default
`link = "log"` keeps the underlying TMB parameter on the log scale), and
it can also carry priors on the linked coefficients.
  
```{r}
# Create a simple environmental covariate for the example data.
set.seed(123)
yrs <- whamGrowthData$styr:whamGrowthData$endyr
whamGrowthData$env_data <- data.frame(
  Year = yrs,
  temp = rnorm(length(yrs), 0, 0.5)
)


growth_spec <- build_growth(
  fun = "vonBertalanffy",
  linkages = list(
    Linf = linkage_spec(
      formula = ~ temp,             # 2 columns: (Intercept) and `temp`
      by      = ~ species,          # one beta per (species, column)
      priors  = list(
        # Intercept prior is on natural-scale Linf (log link back-transforms).
        # Use lognormal() if you want a Normal on log(Linf) instead.
        "(Intercept)" = normal(85, 5),
        temp = normal(0, 0.1)) # Tight prior on the slope
    )
  )
)

```

### What this linkage does

- `formula = ~ temp` fits an intercept and a temperature slope.
- `by = ~ species + sex` gives each species/sex combination its own
intercept and slope.
- `link = "log"` (the default) keeps the underlying TMB parameter on the
log scale, so the slope rows act multiplicatively on natural `Linf`.
- `priors` attaches a Normal prior to both the intercept and the
temperature slope.

The prior on `(Intercept)` is effectively a prior on the base
`Linf` value (natural scale, because the default log link is back-transformed
inside the prior block), and intercept-bearing linkages replace the base
parameter and carry its level.

## Fit the model

```{r}
vbgf_prior_fit <- fit_mod(
  data_list   = whamGrowthData,
  growthFun   = growth_spec,
  estimateMode = 0,
  msmMode     = 0,
  fit_control = fit_control(verbose = 0, phase = FALSE)
)
summary(vbgf_prior_fit)
```

## Inspect the linkage table and prior contributions

```{r}
# The linkage table is stored on the fitted model.
vbgf_prior_fit$data_list$linkage_table

# The design matrix used by the linkages.
head(vbgf_prior_fit$data_list$linkage_X)

# Prior contribution appears in the joint NLL components.
vbgf_prior_fit$quantities$jnll_comp["Linkage-table priors", ]
```

## Compare with an unlinked baseline model and Richards model

```{r}
baseline_fit <- Rceattle::fit_mod(
  data_list = whamGrowthData,
  inits = NULL,       # Initial parameters at default
  estimateMode = 0,   # Estimate
  growthFun = build_growth(fun = "vonBertalanffy"), # Von Bert
  random_rec = FALSE, # No random recruitment
  msmMode = 0,        # Single species mode
  initMode = "NonEquilibrium",       # Unfished disequilibrium
  fit_control = fit_control(phase = FALSE, verbose = 1))


richards_model <- Rceattle::fit_mod(
  data_list = whamGrowthData,
  inits = NULL, # Initial parameters from above
  estimateMode = 0,   # Estimate
  growthFun = build_growth(fun = "Richards"), # Richards
  random_rec = FALSE, # No random recruitment
  msmMode = 0,        # Single species mode
  initMode = "NonEquilibrium",       # Unfished disequilibrium
  fit_control = fit_control(phase = FALSE, verbose = 1))



plot_biomass(
  list(baseline_fit, richards_model, vbgf_prior_fit),
  species = 1,
  model_names = c("Baseline VBGF", "Richards", "VBGF with K linkage and priors"),
  incl_proj = TRUE
)

baseline_fit$opt$AIC
vbgf_prior_fit$opt$AIC
```

## Notes on priors through linkages

- A prior on `"(Intercept)"` in a growth linkage is a prior on the
base growth parameter value itself.
- A prior on a covariate column such as `temp` constrains the effect of
that driver on the parameter.
- The linkage table API supports all R-style formulas and per-species
or per-sex stratification.
- The SD endpoints `sd_L1` and `sd_Linf` (the length-at-age SDs anchored
at `L1` and `Linf`, stored internally on the log scale via `growth_log_sd`)
are linkage targets too -- they thread `init` / `bounds` / `priors` onto
`growth_log_sd`. Only intercept-bearing formulas (typically `~ 1`) are
honored; the SDs do not vary by year in `growth.hpp`, so slope rows have
no effect and slope-only formulas error. See the "Growth SD endpoints"
section of `vignette("environmental-linkages-and-priors")` for an example.