Growth estimation with linkages and priors

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

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.

# 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

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

# 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

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.