Growth estimation with linkages and priors • Rceattle

Overview

This vignette shows how to estimate growth with Rceattle’s built-in von Bertalanffy and Richards growth functions, and how to use the formula-driven linkage_spec() API to add covariates and priors to growth parameters.

The examples use the bundled whamGrowthData dataset, which supplies conditional age-at-length data (caal_data) plus the stock assessment inputs needed to estimate growth.

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

Growth linkages are built with build_growth() and linkage_spec(). A linkage can act on any von Bertalanffy parameter, K, Linf, or L1 (natural-scale names; the default link = "log" keeps the underlying TMB parameter on the log scale), and can hold 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 gives each species its own intercept and slope (add + sex to further split by sex).
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 (Intercept) prior is a prior on the base Linf value (natural scale, because the default log link is back-transformed inside the prior block). An intercept-bearing linkage replaces the base parameter and holds its level.

Time-varying growth

The same grammar holds time variation, with a random-effect term on the right-hand side instead of (or alongside) a covariate:

# A random walk on the von Bertalanffy k. Use ar1(1 | Year) for a stationary
# autoregression, or (1 | Year) for independent annual deviations; under ar1()
# the `sigma` init is the MARGINAL standard deviation.
build_growth(
  fun = "vonBertalanffy",
  linkages = list(K = linkage_spec(~ rw(1 | Year)))
)

This is the only route to time-varying growth. The growth_re switch that earlier versions held in the data_list was never implemented, nothing consumed it, the deviation array was mapped off in every configuration, and the model gave it no density, so it is removed as of 5.9.0 and setting it now raises a deprecation message. Unlike the switch, a random linkage holds a real density, so it is estimated, reported in the linkage table, and redrawn by sim_mod(process = "growth").

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 Linf linkage and priors"),
  incl_proj = TRUE
)

baseline_fit$opt$AIC
vbgf_prior_fit$opt$AIC

Matching a Stock Synthesis growth specification

These options let an estimated-growth model reproduce how Stock Synthesis 3.30 turns a growth curve into an age-length key, weight-at-age and spawning output. Each defaults to Rceattle’s usual behaviour, so nothing changes until you set one.

SS3 control input	Rceattle
Population length bins (data file)	`build_growth(pop_lengths = <lower edges, cm>)`
`CV_Growth_Pattern` 0 (CV a function of length-at-age)	`build_growth(sd_form = "CV")`; `sd_L1` / `sd_Linf` are then CV_young / CV_old
`CV_Growth_Pattern` 2 (SD a function of length-at-age)	`sd_form = "SD"` (default)
`Growth_Age_for_L2 = 999`	`sd_plus_group = "WHAM"` (plus group at CV_old / SD_old)
`Linf_decay` -998 / -999 / d > 0	`plus_group_length = "none"` / `"SS3.24"` / `"decay"` with `plus_group_decay = d`
Maturity option 1 (length logistic), fecundity = weight	control columns `L50_mat_len`, `slope_mat_len` (slope positive: SS3’s sign flipped)
`Growth_Age_for_L1`	`build_growth(growth_age_L1 =)`

For example, a Richards model with SS3 1-cm population bins from 0.5 to 142.5 cm, CVs, and no plus-group adjustment:

whamGrowthData$L50_mat_len   <- 55.43   # cm
whamGrowthData$slope_mat_len <- 0.155   # per cm
growth_ss3 <- build_growth(
  fun = "Richards",
  growth_age_L1 = 1,
  sd_form = "CV",
  plus_group_length = "none",
  sd_plus_group = "WHAM",
  pop_lengths = seq(0.5, 142.5, by = 1)
)

With maturity-at-length set, spawning biomass integrates maturity x weight over the length distribution at spawning; the fitted object reports the per-fish result as quantities$spawn_output.

Notes on priors through linkages

A prior on "(Intercept)" in a growth linkage constrains the base growth parameter itself.
A prior on a covariate column such as temp constrains that driver’s effect on the parameter.
The linkage API accepts any R-style formula and supports per-species or per-sex stratification.
The length-at-age SD endpoints sd_L1 and sd_Linf are also linkage targets, so you can set init / bounds / priors on them. Only intercept formulas (e.g. ~ 1) are honored, the SDs do not vary by year, so slope terms have no effect. See the “Growth SD endpoints” section of vignette("environmental-linkages-and-priors") for an example.