Environmental linkages and priors

Overview

Environmental drivers (temperature, prey density, climate indices) enter stock assessments through recruitment, natural mortality, growth, catchability, and selectivity. Rceattle holds them in a long-format linkage table: one row per estimated coefficient tying a design column to a process parameter. Describe a linkage with linkage_spec(), a formula-driven helper, and pass it to the relevant build_*() function; the process applies the resulting offset to its parameter, year by year.

The linkage path is wired for growth (von Bertalanffy and Richards), natural mortality (M1), recruitment (R0, alpha, beta, R_init), catchability (q), selectivity (each form’s own parameters), and, prior-only, the Dirichlet-multinomial composition weights.

A first example: temperature on K

Von Bertalanffy growth with K responding to a temperature index, one slope per species and sex (whamGrowthData has one of each, so one slope here):

library(Rceattle)
data(whamGrowthData)
plot_data(whamGrowthData)

# Use logistic selectivity for all fleets in this example.
whamGrowthData$fleet_control$Selectivity <- "Logistic"
whamGrowthData$fleet_control$Selectivity_dimension <- "Length"

# 1. Attach an env_data table to the data list. Must include `Year`
#    plus whatever covariates show up in your formulas. Here `temp` is
#    a placeholder; in real use it would be a standardized SST series,
#    PDO index, etc.
yrs <- whamGrowthData$styr:whamGrowthData$endyr
whamGrowthData$env_data <- data.frame(
  Year = yrs,
  temp = rnorm(length(yrs), 0, 0.1)
)

# 2. Build a growth specification with one linkage per parameter that
#    you want to vary. linkage_spec() captures one process parameter's
#    formula, strata (by / species / sex / fleet), link, init, bounds,
#    priors, random-effect handling (formula terms, integrate) and any
#    state-space covariate (observe, obs_sd, obs_sd_est). ?linkage_spec
#    documents each.
growth_spec <- build_growth(
  fun = "vonBertalanffy",
  linkages = list(
    K = linkage_spec(
      formula = ~ temp,             # 2 columns: (Intercept) and `temp`
      by      = ~ species + sex,    # one beta per (species, sex, column)
      priors  = list(temp = normal(0, 0.1)) # Tight prior
    )
  )
)

# 3. Fit as usual.
fit <- fit_mod(
  data_list   = whamGrowthData,
  growthFun   = growth_spec,
  estimateMode = 0, 
  msmMode     = 0,
  fit_control = fit_control(verbose = 0)
)
summary(fit)

This fits \(K_{s,\mathrm{sex},y} = K_{s,\mathrm{sex}}\,e^{\beta_{s,\mathrm{sex}}\,\mathrm{temp}_y}\) for species \(s\); How the pieces fit together below gives the general form.

The linkage_spec() API

linkage_spec(
  formula,
  param     = NULL,
  by        = ~ species,
  species   = NULL,                 # ids or spnames
  sex       = NULL,
  fleet     = NULL,                 # Fleet_codes or Fleet_names
  link      = "log",                # or "identity", "exponential"
  init      = NULL,
  bounds    = NULL,
  priors    = NULL,
  est_phase = 1L
)

Key arguments (see ?linkage_spec for the full list):

• formula: one-sided R formula describing the linear predictor (anything model.matrix() understands: ~ 1, ~ temp, ~ temp + PDO, ~ poly(temp, 2), factor predictors, etc.). A (1 | Year) (IID), rw(1 | Year) (random walk) or ar1(1 | Year) (AR1) term adds a random effect with one deviation per level of the variable after the bar (Random effects below).

• param: target parameter name (e.g. "K"). Filled in automatically when the spec is registered under a linkages = list(K = ...) list key.

• by: stratifying factors (default ~species); use ~species + sex for per-(species, sex) coefficients, or NULL to share a single coefficient set.

• species / sex / fleet: restrict a spec to specific strata. Use with the multi-spec form (Per-species / Per-sex formulas below) to give different species, sexes, or fleets different formulas. species and fleet accept either 1-based ids or the model’s own names (data_list$spnames / fleet_control$Fleet_name); sex accepts integers (1L, 2L) or strings ("Females"/"Males", case-insensitive). Prefer names for fleets: an id that is wrong but in range attaches the linkage to a different fleet and the model still fits, while a misspelled name errors and lists the fleets you do have.

• link: relates the linear predictor to the natural-scale target. "identity" adds to the natural-scale parameter, "log" (the default) multiplies it, and "exponential" raises it to a covariate-dependent power, q^exp(beta * x) — the power is of the parameter, not of the covariate (equations under How the pieces fit together).

  This is Stock Synthesis’s environmental link type 1, and "exponential" is SS3’s own name for it (case 1: // exponential env link). It is not SS3’s q power function, which is a separate option that raises the vulnerable biomass to the power 1 + p (pow(vbio, 1 + p)) – that is what Catchability = "PowerEquation" is reserved for.

  SS3 multiplies the parameter on whatever scale it stores, so that scale decides where the link belongs:

  process	SS3 scale	SS3 type 1 is	Rceattle
  catchability, lognormal index	log (Svy_log_q)	q^exp(beta x)	"exponential"
  catchability, natural-scale index (MVN, MVNORM, Normal, TruncatedNormal)	natural	q * exp(beta x)	"log"; "exponential" warns here
  natural mortality, growth	natural	parm * exp(beta x)	"log"
  recruitment	log (SR_LN(R0))	R0^exp(beta x)	not yet wired

  So a NatM or growth line carrying env-var 1xx bridges with "log", not this, and Rceattle refuses "exponential" there rather than let it through: on a natural-scale parameter it would raise the value to a power instead of scaling it, and on M – where log M is always negative, M being below 1 – that inverts the sign of the covariate effect. Recruitment is the one case where the form is right and simply has no accumulator yet, so the refusal says that instead of offering "log" (which is SS3’s type 2).

  Two further conditions on catchability. The link multiplies log q, so the base has to be estimated: a formula with no intercept, or an intercept fixed at est_phase = 0, masks index_log_q and freezes the value being multiplied, and is refused. A free base is not second-guessed at build time: fit_mod(inits = ) and an intercept init both override Catchability_init, so where log q really starts is only knowable once the model is built. It also cannot share a fleet with a random-effect q linkage, because those deviations accumulate into q_linkage_offset, which sits inside the multiply – they would be scaled by exp(beta * x) while their density still scored them at a constant sigma. For an environmental effect plus q deviations, express the deviations as Time_varying_q, whose index_q_dev the template keeps outside the multiply, as SS3 does its dev blocks.

  Because beta multiplies a log, the size of the effect depends on how far log q sits from zero: at q = 1 the covariate does nothing whatever beta is, and below q = 1 the sign of the effect inverts. Both are SS3’s properties, not Rceattle’s.

  Report beta with its base, never alone. To first order log q_y = log q * (1 + beta * x_y), so what the index informs is the product log q * beta, not beta. With an estimated base the fit is free to cross q = 1 and re-express the same curve with the opposite sign of beta: in the simulate-and-refit check in test-linkage-exponential-link.R, a truth of q = 0.3, beta = +0.4 came back as q = 467.0, beta = -0.0775, preserving log q * beta to 1.1% with the fitted q series correlating 0.984 with the truth. The test asserts that invariant rather than those values, which are from one platform. Interpret and plot the fitted catchability (fit$quantities$index_q), and pin the level with a q prior (Catchability = "Estimated-with-prior") if you need beta itself to be identified. convergence_diagnostics() flags a fit sitting at q = 1 as exponential_q_near_one, a WARN when the fitted catchability is within 0.1% of 1.

  On catchability the effect stops at the last hindcast year. The q linkage tensors are built over hindcast years only (the recruitment, M and growth tensors do span the projection), and a projection-year index is predicted at the last hindcast year’s catchability, so a covariate value for a projection year does not move q there. SS3 carries its environmental link across the forecast (for (y1 = styr - 1; y1 <= YrMax; y1++)), so a bridge that compares forecast index predictions will diverge.

  Set a bound. beta enters a double exponential here, so the usable range is far narrower than for "log", and it depends on the base: the exponent is log q * exp(beta * x), so at q = 0.0012 (log q = -6.7) a beta * x of only 5 gives -994 and index_q underflows to exactly 0 – where both the likelihood and its gradient are dead and no line search can recover. Under "log" the same beta is harmless, because beta * x enters the exponent linearly and overflow starts near 709. Pass bounds on any "exponential" slope, and scale them to the base: underflow starts once |beta * x| > log(709 / |log q|), which at log q = -6.7 is beta * x of about 4.7 – so a bound of c(-2, 2) is not safe there, since a standardized covariate reaches x = 2.4 over a long series.

• init / bounds: named lists keyed by design-matrix column (init = list(temp = 0.05), bounds = list(temp = c(-2, 2))). A slope is on the link scale, an (Intercept) on the parameter’s natural scale — init = list(intercept = 0.3) means an M1 of 0.3/yr, not a log. Mind that split: it is the one place a linkage value can be read on the wrong scale and still fit.

• priors: named list keyed by design-matrix column. Each entry is either an Rceattle_prior or a list keyed by species id; see Priors below.

• est_phase: the estimation phase ordinal (0 fixes, 1+ estimates).

Priors

Four families are available, plus the implicit “none”:

prior_normal(mean, sd)
prior_lognormal(meanlog, sdlog)
prior_gamma(shape, rate)
prior_beta(shape1, shape2)

Inside linkage_spec(priors = ...) the unprefixed normal(), lognormal(), gamma(), beta() resolve to those constructors through a private data mask scoped to the priors argument, so base::gamma() / base::beta() are untouched elsewhere:

linkage_spec(
  formula = ~ temp + PDO,
  by      = ~ species,
  priors  = list(
    temp = normal(0, 1),     # informative, ~95% in (-2, 2)
    PDO  = normal(0, 5)      # weak
  )
)

Prior semantics. On (Intercept) rows the prior is evaluated on the natural scale (e.g. normal(0.3, 0.1) on K means dnorm(K_natural, 0.3, 0.1)). Use lognormal(meanlog, sdlog) to put a normal directly on the log scale; under the default fit_control(bias_adjust_proc = TRUE) it is centred at meanlog - sdlog^2/2, so exp(meanlog) is the prior mean (the median with FALSE). On covariate rows the prior applies to the coefficient as-is, regardless of link function.

For programmatic prior assembly (e.g. species-specific priors built in a loop) call the prior_* constructors directly:

# One prior SD per species, however you derive them.
species_specific_sds <- c(0.3, 0.7, 0.5)
priors <- lapply(species_specific_sds, function(s) prior_normal(0, s))
names(priors) <- c("temp", "PDO", "salinity")
linkage_spec(formula = ~ temp + PDO + salinity, priors = priors)

Each prior’s negative log-density contributes to the Linkage-table priors row of the joint NLL:

fit$quantities$jnll_comp["Linkage-table priors", ]

Putting a prior on the parameter itself

A prior attached to an (Intercept) row is re-targeted to the base parameter, so the same grammar puts an informative prior straight onto a growth parameter (K, L1, Linf, m), onto natural mortality (M1), or onto a recruitment parameter (R0, alpha, beta).

# Example: prior on K for von Bertalanffy growth
growth_spec <- build_growth(
  fun = "vonBertalanffy",
  linkages = list(
    K = linkage_spec(
      formula = ~ 1,                              # intercept only, no env effects
      by      = ~ species + sex,
      init    = list(`(Intercept)` = 0.3),   # start at K = 0.3
      # Use `lognormal(log(0.3), 0.1)` to get a normal prior on log(K).
      priors  = list(`(Intercept)` = lognormal(log(0.3), 0.1))
    )
  )
)

An M1 intercept prior replaces build_M1(M1_use_prior = TRUE) for that species, and fit_mod() refuses the pair unless M2_use_prior = TRUE (a prior on total M).

# Prior on the SD at L1 anchored at exp(log(5)) = 5 length units,
# with a tight normal prior on the log scale.
growth_spec_sd <- build_growth(
  fun = "vonBertalanffy",
  linkages = list(
    sd_L1   = linkage_spec(
      formula = ~ 1,
      by      = ~ species + sex,
      init    = list(`(Intercept)` = 5),
      # Prior on log-scale SD; switch to `normal(5, 0.3)`
      # for a Normal on sd_L1.
      priors  = list(`(Intercept)` = lognormal(log(5), 0.3))
    ),
    sd_Linf = linkage_spec(
      formula = ~ 1,
      by      = ~ species + sex,
      init    = list(`(Intercept)` = 15)
    )
  )
)

# After fit_mod(), the init values set the SD at L1 and the SD at Linf
# (on the log scale), and the priors evaluate against those SD
# parameters via the intercept re-targeting path described above.

rec_spec <- build_srr(
  srr_fun = 0,
  linkages = list(
    R0 = linkage_spec(
      formula = ~ 1,
      by      = ~ species,
      init    = list(`(Intercept)` = 1000),
      # Species-keyed priors on log-scale R0.
      priors  = list(`(Intercept)` = list(
        `1` = lognormal(log(1000), 0.2),
        `2` = lognormal(log(1500), 0.3),
        `3` = lognormal(log(800),  0.25)
      ))
    )
  )
)

fit$quantities$jnll_comp["Linkage-table priors", ]

fit$data_list$linkage_table[fit$data_list$linkage_table$param == "K", ]
fit$estimated_params$beta_linkage          # (Intercept) rows are 0; slopes are estimated
fit$estimated_params$log_growth_pars       # base K / L1 / Linf / m (on log scale internally)

data("GOApollock") # Using GOA pollock as an example

# One-species model, lognormal(log(0.34), 0.01) prior on M1
# applied via the linkage intercept (equivalent to a Normal on log_M1).
# The structural M1_model is kept at "sex_age_invariant" (= 1).
M1_spec <- build_M1(
  M1_model = "sex_age_invariant",
  linkages = list(
    M1 = linkage_spec(
      formula = ~ 1,
      by      = ~ species,
      init    = list(`(Intercept)` = 0.34),
      # Prior is on the natural-scale for M1.
      # `lognormal(log(0.34), 0.01)` evaluates Normal on log(M1).
      # Use `normal(0.34, 0.01)` for a Normal directly on natural M1.
      priors  = list(`(Intercept)` = lognormal(log(0.34), 0.01))
    )
  )
)

# Fit the model - no manual mapping needed
fit <- fit_mod(
  data_list = GOApollock,
  M1Fun     = M1_spec,
  estimateMode = 0,
  fit_control = fit_control(phase = TRUE)
)

Species- and (species, sex)-specific priors

Each entry under priors = list(...) may be:

• a single prior, which applies to every (species, sex) row using that design column;
• a list keyed by species id, one prior per named species; or
• a two-level list keyed by species id and then sex id, one prior per named (species, sex).

# Species-only keying
linkage_spec(
  formula = ~ temp,
  by      = ~ species,
  priors  = list(
    temp = list(
      `1` = normal(0, 0.1),    # tight prior for species 1
      #`2` = normal(0, 0.5),   # No prior species 2
      `3` = normal(0, 1.0)     # weakest for species 3
      # any species not listed gets no prior
    )
  )
)

For sexually dimorphic stocks the same temp column can take different priors per (species, sex). At each level the value is either an Rceattle_prior, which applies to every sub-stratum, or a list that drills down further:

linkage_spec(
  formula = ~ temp,
  by      = ~ species + sex,
  priors  = list(
    temp = list(
      # Species 1: different tightness per sex.
      `1` = list(`1` = normal(0, 0.1),    # sp 1, sex 1 (e.g. females)
                 `2` = normal(0, 0.2)),   # sp 1, sex 2 (e.g. males)
      # Species 2: scalar applies to every sex.
      `2` = normal(0, 0.5)
      # Any (species, sex) combination not listed gets no prior.
    )
  )
)

The same syntax works for any prior family. Reach for it when prior strength differs across stocks — one species has published temperature sensitivity and another has none — or across the sexes within a stock. Built programmatically:

sd_by_species <- c(`1` = 0.1, `2` = 0.5, `3` = 1.0)
priors_per_sp <- lapply(sd_by_species, function(s) prior_normal(0, s))

linkage_spec(
  formula = ~ temp,
  by      = ~ species,
  priors  = list(temp = priors_per_sp)
)

Per-species formulas

Where species have qualitatively different relationships — species 1 responds to temperature only, species 2 also to PDO — register multiple specs against the same parameter and scope each with its species argument:

build_growth(
  fun = "vonBertalanffy",
  linkages = list(
    K = list(
      linkage_spec(
        formula = ~ temp,
        by      = ~ species,
        species = 1L,                                    # only sp 1
        priors  = list(temp = normal(0, 0.5))
      ),
      linkage_spec(
        formula = ~ temp + PDO,
        by      = ~ species,
        species = 2L,                                    # only sp 2
        priors  = list(temp = normal(0, 1),
                       PDO  = normal(0, 1))
      )
    )
  )
)

Each parameter key accepts both forms:

linkages = list(K = single_spec)             # one formula for all sp
linkages = list(K = list(spec_a, spec_b))    # different formulas per sp

A build_*() output’s linkages$K is either an Rceattle_linkage_spec or a list of them; after fitting, fit$data_list$linkage_table shows the materialized rows with species and X_col.

Each spec holds its own priors, which can themselves be species-keyed:

linkage_spec(
  formula = ~ temp,
  by      = ~ species + sex,
  species = c(1L, 2L),
  priors  = list(
    temp = list(
      `1` = normal(0, 0.3),
      `2` = normal(0, 0.7)
    )
  )
)

Per-sex formulas

The sex argument mirrors species: it scopes a spec to a subset of sex ids, so the multi-spec form can register separate specs per sex against one parameter — sex-dimorphic M priors, or different env responses for females and males:

# Female and male M1 each get their own intercept-only prior.
M1_spec <- build_M1(
  M1_model = "sex_specific",
  linkages = list(
    M1 = list(
      linkage_spec(
        formula = ~ 1,
        by      = ~ species + sex,
        sex     = "Females",                              # or sex = 1L
        init    = list(`(Intercept)` = 0.12),
        priors  = list(`(Intercept)` = lognormal(log(0.12), 0.0001))
      ),
      linkage_spec(
        formula = ~ 1,
        by      = ~ species + sex,
        sex     = "Males",                                # or sex = 2L
        init    = list(`(Intercept)` = 0.20),
        priors  = list(`(Intercept)` = lognormal(log(0.20), 0.0001))
      )
    )
  )
)

The sex filter applies only when by includes sex; otherwise fit_mod() warns that it has no effect, and a filter matching none of the model’s levels warns that the spec is dropped. Reach for the multi-spec form when the sexes need different RHS terms, init/bounds, or est_phase. For different prior values on one shared formula, the nested-priors form above is simpler: one spec, one set of design columns, priors keyed by species then sex.

Polynomial and basis-expansion formulas

The formula passed to linkage_spec() is fed straight to stats::model.matrix(), so any expression it understands works: polynomial bases, splines, interactions, and arithmetic transformations.

Polynomials. poly() gives an orthogonal basis and poly(..., raw = TRUE) raw monomials. Each order is its own design column with its own coefficient:

# 4th-order orthogonal polynomial in temp, by species.
linkage_spec(
  formula = ~ poly(temp, 4),
  by      = ~ species
)
# design columns produced:
#   (Intercept), poly(temp, 4)1, poly(temp, 4)2,
#   poly(temp, 4)3, poly(temp, 4)4

poly() orthogonalises against every row of env_data supplied, including any beyond the projection horizon, so the columns are zero-mean over that span and the linear slope, the quadratic curvature and so on are estimated independently. Prefer it unless you have a scientific reason to want raw monomials.

Arbitrary arithmetic. Use I() to insert a literal expression (otherwise temp^2 is interpreted as a formula crossing operator and reduces to temp):

# Linear + raw quadratic in temp.
linkage_spec(
  formula = ~ temp + I(temp^2),
  by      = ~ species,
  priors  = list(`I(temp^2)` = normal(0, 0.1))   # tight prior on curvature
)

The design column name is the deparsed call (I(temp^2)), so init, bounds and priors must key that same string — backticks let you write it as a name.

Splines. splines::ns() gives natural cubic splines: flexibility without a high-order polynomial’s tail behaviour:

linkage_spec(
  formula = ~ splines::ns(temp, df = 4),
  by      = ~ species
)

Interactions. temp:PDO adds a single product column; temp * PDO is shorthand for temp + PDO + temp:PDO:

linkage_spec(
  formula = ~ temp * PDO,
  by      = ~ species,
  priors  = list(`temp:PDO` = normal(0, 0.5))
)

Mixing forms across species. With the per-species formula machinery one species can take a polynomial response and another a linear one, without duplicating columns:

build_growth(
  fun = "vonBertalanffy",
  linkages = list(
    K = list(
      linkage_spec(formula = ~ poly(temp, 3),
                   by = ~ species, species = 1L),
      linkage_spec(formula = ~ temp,
                   by = ~ species, species = 2L)
    )
  )
)

Same-named columns must be numerically identical across the specs that use them. Deterministic formulas evaluated against one env_data give that automatically; only specs built against different env_data frames can break it.

Time blocks

A time block is an ordinary factor on the right-hand side, and cut() makes one from Year. Three blocks running through 1990, 1991 to 2010, and 2011 on:

# Blocks on the ascending inflection of a logistic fishery selectivity. The
# legacy Block deviate is added to the inflection, so match it with identity.
build_selectivity(linkages = list(
  inf_asc = linkage_spec(~ cut(Year, breaks = c(-Inf, 1990, 2010, Inf)),
                         by = ~ fleet, fleet = "GOA_pollock_fishery",
                         link = "identity")))

cut() closes each interval on the right, so a break year joins the block before it: 1990 is in the first block, not the second.

model.matrix() folds the first block into (Intercept), so three blocks give two offset columns and each coefficient is the step away from the first block’s level. ~ 0 + cut(Year, ...) gives one coefficient per block instead, and with it a base parameter held fixed (see Behind the scenes below). The column name adds the level label to the deparsed call, as in cut(Year, breaks = c(-Inf, 1990, 2010, Inf))(1.99e+03,2.01e+03], which init, bounds and priors key exactly. cut()’s default dig.lab = 3 is what rounds 1990 to 1.99e+03; pass dig.lab = 4 for (1990,2010].

Time_varying_sel = "Block" and Time_varying_q = "Block" fit blocks this way too, taking them from the Selectivity_block column of catch_data or index_data: one free deviate per block, with no penalty on it. The selectivity switch blocks a fleet’s slope and inflection at once, with separate deviates per sex, and fixes both base parameters, so a close rewrite names each parameter, passes by = ~ fleet + sex, and uses the ~ 0 + form. A fixed-effect coefficient is free either way, so a block needs no integrate = FALSE. An IID or RandomWalk rewrite does.

Any factor works. A Selectivity_block id reused across non-adjacent year spans is not an interval, so give env_data a per-year block column and write ~ factor(block).

A block needs rows in env_data. The design column is built over the rows supplied, so a block with no rows gets an all-zero column and a coefficient nothing identifies, and a hindcast year past the last row takes a zero offset, leaving the base parameter unshifted. Selectivity and catchability act in the hindcast years only, so rows past endyr move neither; M, growth and recruitment blocks do run through the projection.

There is no block() helper. A block term would be a second grammar for what cut() already expresses in base R.

The right-hand side is evaluated against env_data. The fixed part goes to model.matrix(data = env_data), so ~ cut(Year, ...) needs a Year column, which is otherwise optional. species, sex and fleet are by arguments rather than design columns, so ~ temp * fleet fails with object 'fleet' not found. model.matrix() falls back to the calling environment, so a stray object of that name becomes a design column instead.

How the pieces fit together

fit_mod() pools the linkages from every build_*() call into one linkage table and one shared design matrix (both stored on data_list), then adds the per-process offset to the base parameter each year; priors and bounds travel in the same table and add to the joint negative log-likelihood. With no linkages every offset is zero, so an unlinked model is unaffected.

Each row of the table is one coefficient \(\beta_i\) on one covariate column \(x_i(y)\) of the design matrix, for one parameter and stratum \(s\) (species, sex, age bin, fleet). The (Intercept) row is held at 0, so the base parameter \(p_s\) keeps the level, and the intercept’s init, bounds and prior go to \(p_s\) (the exception is R_init; see Recruitment). For year \(y\), write \(S_{\log}\), \(S_{\mathrm{id}}\) and \(S_{\exp}\) for \(\sum_i \beta_i x_i(y)\) over the stratum’s log-, identity- and exponential-link rows:

A parameter’s own deviations (those M1_re, Time_varying_q or Time_varying_sel turn on) still apply. With both a log and an identity link on one parameter, the identity sum is added after the multiplication, \(p_{s,y} = p_s e^{S_{\log}} + S_{\mathrm{id}}\). On the selectivity inflection slots and the DoubleNormalSS3 parameters it is added before, \(p_{s,y} = (p_s + S_{\mathrm{id}})\, e^{S_{\log}}\), on the scale each is stored on: a length (cm), or on an age-based fleet a bin number (age - minage + 1), for an inflection or peak; a logit for DoubleNormalSS3 P2, P5, P6 and DoubleNormal’s floor; a log for DoubleNormalSS3’s widths (P3, P4) and LogisticPM’s first-bin selectivity. There a log link multiplies the logit or the log.

In the first example, \(K_{s,\mathrm{sex},y} = K_{s,\mathrm{sex}}\,e^{\beta_{s,\mathrm{sex}}\,\mathrm{temp}_y}\): a log-scale slope, so \(\beta = 0.1\) raises \(K\) by about 10.5% (\(e^{0.1} - 1\)) per unit of temperature.

A random effect such as (1 | Year) adds one row per year in each stratum, each with an indicator column, whose \(\beta\) is that year’s deviation; each stratum’s deviations share one SD and a structure (IID, random walk or AR1), and their density is part of the objective.

Natural mortality

build_M1() takes the same linkages argument with one key, M1. The offset enters additively on the log scale, so a coefficient of 0.3 on a normalized temperature index multiplies M1 by exp(0.3 * temp[yr]) each year.

Why M1 and not M? CEATTLE decomposes total natural mortality as M = M1 + M2. Only M1 (residual / non-predation) is a parameter; M2 (predation) is derived from predator abundance, suitability, and ration each iteration. An environmental effect on predation therefore propagates through upstream parameters: a temp effect on predator K, on recruitment, or on bioenergetics ration inputs.

yrs <- GOApollock$styr:GOApollock$endyr
GOApollock$env_data <- data.frame(
  Year = yrs,
  temp = rnorm(length(yrs), 0, 0.1)
)

M1_spec <- build_M1(
  M1_model = "sex_age_invariant",            # or M1_model = 1 (back-compat)
  linkages = list(
    M1 = linkage_spec(
      formula = ~ temp,
      by      = ~ species,
      priors  = list(temp = normal(0, 0.1))
    )
  )
)

fit <- fit_mod(
  data_list   = GOApollock,
  M1Fun       = M1_spec,
  estimateMode = 3,   # build the objective; see ?fit_mod
  msmMode     = 0
)

The structural choice (M1_model) and the random-effect choice (M1_re) accept integer codes or readable strings interchangeably:

build_M1(M1_model = "sex_age_invariant", M1_re = "ar1_age")   # = 1, 4
build_M1(M1_model = 1,                   M1_re = 4)           # equivalent

Random effects on M1 need a level that leaves them identifiable, so M1_re > 0 is available for M1_model 0 ("fixed"), 1 and 2, and for 3 only with a by-year family (2 or 5). Under "fixed" the input schedule is held and the deviations carry M’s variation on their own. Under 3 the level is already age-specific, so the age-specific log_M1 and an age-varying deviation enter the likelihood only as their sum and are not separately identifiable. For an environmentally driven M with deviations, use the linkage grammar: linkage_spec(~ temp + (1 | Year)).

linkages on M is age-aware: by default a row’s age_bin is NA and the offset broadcasts across every age. To pin an effect to one age slice, set by = ~ species + age_bin and supply strata$age_bin, a bin index (age - minage + 1) rather than an age; fit_mod() supplies it automatically, so only a direct materialize call needs it.

The M1_indices column-index argument for M1_model = 4 / 5 adds the same log-scale environmental effect and is soft-deprecated. Prefer linkages: it composes with priors, bounds, and per-species formulas. Supplying both for one coefficient double-counts the effect.

Inspecting the result:

fit$data_list$linkage_table         # rows have process == "M"
fit$quantities$M_linkage_offset     # [nspp, max_sex, max_age, nyrs]
fit$quantities$jnll_comp["Linkage-table priors", ]

Recruitment

build_srr() accepts the same linkages argument with four parameter keys:

• R0: acts under mean recruitment (srr_fun = "mean", penalty form included); the right target if recruitment should respond to env without fitting an SRR. Under a curve fitted in the hindcast a single-species R0 is derived from alpha and beta, so an R0 linkage is refused. In a multispecies model R0 sets only the initial level (R_init), so only an intercept-only linkage (~ 1, a prior or fixed value) is accepted there.
• alpha (productivity) and beta (density-dependence) do work only where the model has a curve: srr_fun (or, for the penalty form, srr_pred_fun) is "BevertonHolt" or "Ricker".
• R_init: the initial recruitment level, a log-scale multiplier on the initial age-structure, read at year 0 because the initial state is one year. Use it where the series begins after a regime shift, so the stock starts at a recruitment other than \(R_0\): without it the level has to go into init_dev, which is penalised per age, and the fit pays a recruitment-deviate penalty for sitting where the data say it sits. Unlike the three above it has no base parameter in rec_pars, so an intercept-only ~ 1 keeps the level on its own estimated coefficient — and init / bounds on that intercept are the natural-scale multiplier on \(R_0\) (init = 0.25 is a quarter), logged onto the coefficient in place of a base parameter. One design column per species, since only year 0 is read. It needs link = "log", and is refused under initMode = "FreeParams" (which never reads \(R_{init}\)) and "OffsetEquilibrium" (which already scales the same ages by rec_dev[, 1]).

The level is informed only by the first year’s observations, so it is weakly identified at the low end: once the initial cohorts are effectively annihilated the data cannot tell one very small level from another, and a fit can settle there with a clean gradient. Bound it on a real assessment — the bound is the natural-scale multiplier — and read the estimated multiplier rather than the convergence flag alone. tools/verify/verify-sim-recovery-r-init.R measures both the recovery and the rate at which a replicate runs away.

# A free initial level, Rceattle's form of Stock Synthesis's SR_regime block.
# SS3 shrinks that block at sigma_R / ave_age; this is free, so add a
# lognormal(0, sigma_R / ave_age) prior on the intercept to match SS3.
# It must travel through build_srr(): fit_mod() overwrites data_list$srr_linkages
# from recFun, so assigning the spec onto the data_list drops it silently.
init_spec <- build_srr(linkages = list(R_init = linkage_spec(
  formula = ~ 1,
  init    = list(`(Intercept)` = 0.25),   # start at a quarter of R0
  bounds  = list(`(Intercept)` = c(0.05, 20)))))

# Mean-recruitment with a temperature linkage on R0.
rec_spec <- build_srr(
  srr_fun  = "mean",                       # = 0
  linkages = list(
    R0 = linkage_spec(
      formula = ~ temp,
      priors  = list(temp = normal(0, 0.5))
    )
  )
)

# A three-species data_list for the fits below, still fitted in single-species
# mode (msmMode = 0): the species-keyed priors need more than one species, and
# the linkages need `temp` in env_data.
data_list <- BS2017SS
yrs <- data_list$styr:data_list$endyr
data_list$env_data <- data.frame(
  Year      = yrs,
  temp      = as.numeric(scale(sin(seq_along(yrs) / 4))),
  cold_pool = as.numeric(scale(cos(seq_along(yrs) / 5))))

# Beverton-Holt with env effects on alpha (and species-keyed
# priors so each stock has its own sensitivity).
rec_spec_bh <- build_srr(
  srr_fun  = "BevertonHolt",               # = 2
  linkages = list(
    alpha = linkage_spec(
      formula = ~ temp,
      priors  = list(
        temp = list(`1` = normal(0, 0.3),
                    `2` = normal(0, 0.7))
      )
    )
  )
)

fit <- fit_mod(
  data_list   = data_list,
  recFun      = rec_spec,
  estimateMode = 3,   # build the objective; see ?fit_mod
  msmMode     = 0
)

# Lognormal prior on alpha (mean 5, log-scale SD 0.5), species 1 only.
rec_prior <- build_srr(
  srr_fun  = "BevertonHolt",
  linkages = list(
    alpha = linkage_spec(~ 1, species = 1,
                         priors = list(`(Intercept)` = prior_lognormal(log(5), 0.5)))
  )
)

# Alpha fixed at 20: an intercept held at its `init`.
rec_fixed <- build_srr(
  srr_fun  = "BevertonHolt",
  linkages = list(
    alpha = linkage_spec(~ 1, init = list(`(Intercept)` = 20), est_phase = 0)
  )
)

build_srr(srr_fun = "mean")           # = 0
build_srr(srr_fun = "BevertonHolt")   # = 2
build_srr(srr_fun = "Ricker")         # = 4

fit$data_list$linkage_table             # rows have process == "recruitment"
fit$quantities$recruitment_linkage_offset  # [nspp, n_rec_params, nyrs]
                                           # n_rec_params = 3 (R0, alpha, beta)

Catchability

Survey/index catchability is indexed by fleet, not by species/sex/age, so its linkages use by = ~ fleet. build_catchability() takes the same linkages argument with a single parameter key, q:

# A q linkage shifts an ESTIMATED log q, so restrict it to fleets that have
# one: a fishery carries no survey index, and a Fixed or Analytical q is
# refused because the linkage would either have no effect or make a fixed q
# time-varying.
q_data <- data_list
srv <- q_data$fleet_control$Fleet_name[q_data$fleet_control$Fleet_type == 2]
q_data$fleet_control$Catchability[q_data$fleet_control$Fleet_type == 2] <- "Estimated"

fit <- fit_mod(
  data_list = q_data,
  qFun      = build_catchability(
    linkages = list(
      q = linkage_spec(~ temp, by = ~ fleet, fleet = srv)  # one per survey fleet
    )
  ),
  estimateMode = 0,
  msmMode      = 0
)

# One coefficient set per fleet; the offset lands here:
fit$quantities$q_linkage_offset      # [n_flt, nyrs_hind]
fit$quantities$q_linkage_log_mult    # [n_flt, nyrs_hind]; exponential link only

Restrict a spec to particular fleets with the fleet argument, and attach priors as for any other process. Name the fleets rather than numbering them: names are checked against fleet_control$Fleet_name when the model is assembled, whereas a Fleet_code typo that stays in range quietly links the wrong survey:

build_catchability(linkages = list(
  q = linkage_spec(~ temp, by = ~ fleet, fleet = c("Pollock", "EIT_Pollock"),
                   priors = list(temp = normal(0, 1)))))

The Catchability = "Environmental" / Time_varying_q column pair does the same thing, naming covariates by their position in a comma-separated list of env_data column indices; such a model still fits, with a deprecation warning pointing at build_catchability(). As a linkage the effect is a formula, holds priors/bounds/phase, and names the fleet explicitly.

An intercept-bearing formula (~ temp) leaves the base catchability estimable to hold the level and fixes the intercept row at 0, as elsewhere.

build_catchability(linkages = list(
  # Survey q for fleet 2 believed near 1, with moderate uncertainty.
  q = linkage_spec(~ 1, by = ~ fleet, fleet = 2,
                   priors = list(`(Intercept)` = lognormal(log(1), 0.2)))))

Selectivity

Selectivity is also fleet-indexed (by = ~ fleet) and, unlike the other processes, has several parameters.

Because ~ fleet is the default, fleets sharing a Selectivity_index or a Catchability_index estimate one coefficient for the whole group and keep identical selectivity — provided the linkage reaches every member (5.57.0). Before that each member got its own coefficient and the copies diverged, which is the defect the shared index exists to prevent.

The proviso is the thing to remember. The base parameter is shared, but the offset a linkage adds is applied per fleet, so a fleet = filter naming some members moves only those: the group shares a coefficient and still fits different curves. Name every member, or pass by = NULL for a single row that reaches all of them. Four things are refused on a mirrored block: naming a strict subset; a random-effect linkage (each member would get its own deviation series and SD); two rows for one design column on one fleet (the offset is added once per row, so that fleet moves twice as far); and two specs disagreeing on link, bounds or init. est_phase may differ — a fixed donor holds the group, which is how a reference model that fixes some time blocks is reproduced.

A linkage on the lead fleet alone does not propagate, and that includes the apical height. The one thing re-targeted onto the shared base is a prior, which is why it belongs on a single member: named on two it is counted twice on the coefficient they share, halving its stated variance. The usual spelling is a prior on one member and a plain row for the others. An intercept is also re-targeted, so an intercept init or bounds goes on the donor.

build_selectivity() accepts these keys, which name the slots shared across the parametric forms under whichever alias your fleet’s form actually has:

A per-sex apical height. Fishing mortality is one log_F per fleet and year shared by the sexes, so a sex difference in F can only come from selectivity, and no form has a height parameter of its own. apical multiplies one sex’s whole curve by exp(log_sel_apical), bin by bin, after the form and before normalization, so it works for every estimated form. Name the fleet and the sex that holds it; the other sex is the reference, as in Stock Synthesis’s male-offset option:

build_selectivity(linkages = list(
  apical = linkage_spec(~ 1, by = ~ fleet + sex, fleet = 3, sex = "male",
                        priors = list(intercept = lognormal(0, 0.5)))))

It adds to the form rather than replacing part of it: the form’s own parameters are estimated separately for each sex already, so a two-sex Logistic fleet estimates a slope and an inflection per sex plus this multiplier, and the named sex’s curve has three parameters to the reference sex’s two. An SS male offset usually ties the male shape to the female’s; this does not.

The multiplier equals the ratio of the sexes’ peak heights only where their shapes peak equally, which the logistic family does on an age axis; for a dome with sex-specific shape it is the pointwise multiplier, and the peak ratio is read from fit$quantities$sel_at_age. Only the contrast between the sexes is identified, so the fit is refused when both sexes or no sex hold it, when no fleet is named, on a one-sex species, on a Fixed, AR1 or Off fleet, on a fleet that takes another’s block through a shared Selectivity_index (name every fleet in the group, or pass by = NULL), with link = "identity", and where Sel_norm_scope = "WithinSex" normalization would divide the offset out (use "AcrossSexes", under which the more-selected sex peaks at 1, or turn Sel_norm_bin off). The contrast is informed only by joint compositions (comp_data$Sex = 3); with single-sex compositions it rests on its prior. On GOAatf’s fishery (Logistic, across-sex normalization) the male multiplier fits at 2.15 (log offset 0.77, SE 0.21), which is also the peak ratio there because both logistic curves saturate.

# Only the parametric forms are wired for selectivity linkages, so name the
# fleets that use one -- here the Logistic surveys.
logistic_flts <- data_list$fleet_control$Fleet_name[
  data_list$fleet_control$Selectivity %in% c(1, "Logistic")]

fit <- fit_mod(
  data_list = data_list,
  selFun    = build_selectivity(
    linkages = list(
      inf_asc = linkage_spec(~ cold_pool, by = ~ fleet, fleet = logistic_flts)
    )
  ),
  estimateMode = 0,
  msmMode      = 0
)

Each spec carries its own link, so one fit can multiply the log-space slope (link = "log") and add to the natural-space inflection (link = "identity").

Supported forms: the parametric selectivities — Logistic, DoubleLogistic, DescendingLogistic, DoubleNormal, LogisticPM, and DoubleNormalSS3, which reads only its own six parameters (dn_peak, top_logit, ascend_se, descend_se, start_logit, end_logit). A linkage on any other form errors at fit time naming the fleet, rather than being estimated to no effect. apical is the exception: it multiplies the finished curve, so every estimated form takes it. The non-parametric forms (coff) are excluded on principle: they mean-centre their coefficients every year, so a per-year offset applied across all bins cancels exactly, and the per-bin covariate that would be needed is not something the formula grammar expresses. Time_varying_sel process-error deviates are a separate mechanism and compose additively with a selectivity linkage.

Data weighting (Dirichlet-multinomial)

The Dirichlet-multinomial (DM) likelihood (Thorson et al. 2017) is a self-weighting alternative to Francis (2011) or McAllister-Ianelli (1997) tuning: it estimates each composition dataset’s effective sample size inside the fit rather than iterating an external weight. Where a fleet has only a handful of comp years that weight is poorly identified and can be pushed to an implausible extreme; a prior on it keeps it in a believable range without returning to hand-tuning:

# Use the DM likelihood on the fleets you want it to self-weight ...
data_list$fleet_control$Comp_distribution[c(1, 2)] <- "DirichletMultinomial"

# ... and (optionally) constrain its weight with a prior.
fit <- fit_mod(
  data_list = data_list,
  compFun   = build_composition(
    linkages = list(
      theta_comp = linkage_spec(
        ~ 1, by = ~ fleet, fleet = c(1, 2),
        # Gamma(2, 0.5): mean 4, mildly informative. Keeps the weight off 0
        # (comps ignored) and off very large (comps fit as if unweighted).
        priors = list(`(Intercept)` = gamma(2, 0.5)))
    )
  ),
  estimateMode = 0,
  msmMode      = 0
)

Three keys, one per kind of composition data:

The prior is on the weight itself (a larger weight = a larger effective sample size = the composition data trusted more). It only makes sense where the DM is actually estimating that weight, so a theta_* linkage on a fleet or predator whose likelihood is not DirichletMultinomial errors up front rather than sitting inert. And because the weight is a single number per fleet, not a time series, this is a prior only: a covariate (~ temp), an init or an est_phase on a composition spec is rejected. Start or fix the weight through the model’s usual inits / map if you need to.

A multi-process example

The same pattern composes when more than one parameter is linked:

growth_spec <- build_growth(
  fun = "Richards",
  linkages = list(
    K    = linkage_spec(
      formula = ~ temp,
      by      = ~ species + sex,
      priors  = list(temp = normal(0, 1))
    ),
    Linf = linkage_spec(
      formula = ~ temp + PDO,
      by      = ~ species,
      priors  = list(temp = normal(0, 0.5),
                     PDO  = normal(0, 0.5))
    ),
    m    = linkage_spec(
      formula = ~ 1,             # intercept only, no env effect
      by      = ~ species
    )
  )
)

At fit time these specs are pooled. The result:

• one shared design matrix X = model.matrix(~ temp + PDO, env_data) of three columns,
• 16 linkage-table rows: 4 (sp × sex) + 3 (Linf cols) × 3 (sp) + 3 (m intercepts),
• same-named columns ((Intercept), temp) matched across specs, so nothing is duplicated in X.

Per-process logic — only Richards has m — is enforced at build time:

build_growth(
  fun = "vonBertalanffy",
  linkages = list(m = linkage_spec(~ temp))
)

Growth, M and recruitment can all be linked in the same fit:

# whamGrowthData, not data_list: estimating growth needs caal_data (the age
# composition within a length bin is what informs the growth curve), and
# BS2017SS carries none.
fit <- fit_mod(
  data_list   = whamGrowthData,
  growthFun   = build_growth(
    fun      = "vonBertalanffy",
    linkages = list(
      K = linkage_spec(formula = ~ temp)
    )
  ),
  M1Fun       = build_M1(
    M1_model = "sex_age_invariant",
    linkages = list(
      M1 = linkage_spec(formula = ~ temp)
    )
  ),
  recFun      = build_srr(
    srr_fun  = "mean",
    linkages = list(
      R0 = linkage_spec(formula = ~ temp)
    )
  ),
  estimateMode = 3,   # build the objective; see ?fit_mod
  msmMode     = 0
)

# Linkage table contains rows for all three processes:
table(fit$data_list$linkage_table$process)
#>      growth           M  recruitment
#>           4           4            4

Inspecting and debugging

Diagnostics worth knowing:

# The materialized table, one row per estimated coefficient.
fit$data_list$linkage_table

# The shared design matrix, one column per unique RHS term.
fit$data_list$linkage_X

# Per-(species, sex, growth_param, year) offsets actually applied to
# growth_parameters this iteration. Useful for verifying the linkage
# fired as you intended.
str(fit$quantities$growth_linkage_offset)
str(fit$quantities$M_linkage_offset)
str(fit$quantities$recruitment_linkage_offset)

# Per-row coefficient values in the same order as `linkage_table`.
fit$estimated_params$beta_linkage

# The prior contribution from the "Linkage-table priors" row of jnll_comp.
fit$quantities$jnll_comp["Linkage-table priors", ]

Behind the scenes: automatic base-parameter handling

For three species × two sexes × ~ temp (intercept + slope) the pooled table has 12 rows, an (Intercept) and a temp slope per (species, sex). Priors and inits supplied for (Intercept) flow to the base parameter (K, M1, R0, …); the slope rows hold the year-by-year offset that combines with the base. What happens to the base parameter depends on whether the formula holds an intercept:

Intercept-bearing formulas (~ 1, ~ temp): the base parameter stays estimable — R0 for species s lives on rec_pars[s, 1] exactly as it does without any linkages. The (Intercept) row is a hook for init and priors; its beta_linkage slot is held at 0 and not estimated, so the recruitment, mortality and growth processes read base + offset(yr) with the offset being exactly the slope contribution.

Slope-only formulas (~ 0 + temp): there is no (Intercept) row, so the base parameter is held fixed, at its default initial value, at the value supplied by fit_mod(inits = ), or at the spec’s own init for the intercept.

# Start M1 at 0.06 without estimating an extra intercept beta.
build_M1(
  M1_model = "sex_age_invariant",
  linkages = list(
    M1 = linkage_spec(
      formula = ~ 1,
      init    = list(`(Intercept)` = 0.06)
    )
  )
)

Random effects

A linkage formula may hold a random effect, so a parameter varies year-to-year (or over any grouping) as a set of deviations damped by an estimated density rather than fixed covariate slopes:

• ~ (1 | Year), IID deviations, N(0, sigma);
• rw(1 | Year): a random walk, N(0, sigma) on first differences (first deviate fixed so the base parameter holds the level, fixed at its inits value, 0 unless you supply another);
• ar1(1 | Year): a stationary AR1, \(u_t = \rho\, u_{t-1} + e_t\), with sigma the marginal SD of \(u_t\) and rho the correlation.

These structures are the general way to write a time-varying parameter, and cover what the Time_varying_q, Time_varying_sel and M1_re switches do. A survey catchability that used Time_varying_q = "RandomWalk" becomes

build_catchability(linkages = list(
  q = linkage_spec(~ rw(1 | Year), by = ~ fleet)))

time-varying selectivity becomes a linkage on the relevant selectivity parameter (e.g. inf_asc = linkage_spec(~ rw(1 | Year), by = ~ fleet)), and a time-varying natural-mortality M1_re model becomes build_M1(linkages = list(M1 = linkage_spec(~ (1 | Year)))). The Time_varying_* / M1_re switches are soft-deprecated: they still fit with their exact numerics and fit_mod() emits a nudge when it sees them, while the grammar adds a prior on the deviation SD, or a free estimate of it, in place of a fixed input value.

sigma means different things across structures. Under rw() it is the innovation (step) SD, the SD of each year-to-year change, a random walk being non-stationary with no overall spread. Under ar1() it is the marginal (stationary) SD, the overall spread of the deviations (sigma_marginal = sigma_innovation / sqrt(1 - rho^2)), following the glmmTMB convention. The two agree at rho = 0, where ar1 reduces to IID, so a sigma estimated under rw() is not comparable with one under ar1().

Each grouping level gets one deviation; each random term in each stratum estimates its own SD. The SD and (for ar1) the correlation route through the same init / priors arguments as every other parameter: init = list(sigma = 0.1) fixes the SD (the same as fixing Time_varying_*_sd), init = list(rho = 0.7) fixes the correlation, and priors = list(sigma = lognormal(...)) / priors = list(rho = normal(...)) place priors. rw() / ar1() need a numeric grouping variable (a real elapsed-time lag).

build_selectivity(linkages = list(
  slp_asc = linkage_spec(~ rw(1 | Year), fleet = "GOA_pollock_fishery",
                         init = list(sigma = 0.05), integrate = FALSE)))

Current limitations

• Time-variation switches are not auto-translated. The soft-deprecated Time_varying_q / Time_varying_sel / M1_re switches (above) are not rewritten into the grammar for you; rewrite the spec by hand when you want the grammar’s features. The switches estimate their deviations as penalized fixed effects with a fixed SD, so a faithful rewrite pairs rw() with integrate = FALSE and init = list(sigma = ) set to the switch’s Time_varying_*_sd; an integrated rw() is a different model. Two traps when rewriting a selectivity switch: the legacy penalty weights the inflection deviate at 4x the slope deviate (and swaps which limb gets the 4x between RandomWalk and RandomWalkAscending), so read the sigmas off ceattle.cpp’s penalty block rather than assuming; and an inflection linkage needs link = "identity", because the legacy deviate is added to the inflection whereas the default log link multiplies it. The nudge is silent on Catchability = "Environmental" and Catchability = "AR1", where Time_varying_q holds an env_data column index rather than a mode: the first carries its own deprecation warning, and the second (the Rogers QAR1 form) is refused outright because it never fitted. It is silent too on the non-parametric selectivity forms, which have no additive slot, and on the separable M1_re = 6 (age × year), which has no 1-D grammar structure; both keep the switch-based path.
• Only the grammar’s random effects are simulated. sim_mod(process = ) and self_test(process = ) redraw catchability and selectivity where they are written as a random linkage; the legacy Time_varying_q / Time_varying_sel deviations hold densities in the model but no draw, so they keep their fitted values and sim_mod() warns that they did. Recruitment and M are drawn from their own switches (M1_re) as well as from linkages. Observed AR1 (Rogers QAR1) linkages are the one grammar form that is deliberately not redrawn: the latent state is measured by an observed covariate series, and redrawing the state alone would leave the two describing different histories. A switch-based time-varying process therefore cannot have that process redrawn in a self-test.
• sigma means different things in the two systems. In a linkage, init = list(sigma = ) under ar1() is the marginal standard deviation. The M1_re and 2DAR1/3DAR1 selectivity switches instead take the AR(1) innovation SD, and scale it to a marginal of sigma / sqrt(1 - rho^2), the WHAM convention. The two coincide only at rho = 0.
• Age-dependent linkages on processes other than M: only the M accumulator consumes the age_bin stratum (a row’s age_bin == NA broadcasts across ages; a specific value pins the offset to that age slice). The schema is reserved for other processes that may eventually vary by age.
• Non-parametric selectivity: a coff linkage needs a per-bin covariate, since a per-year scalar cancels under the annual mean-centring; the formula grammar does not yet express that.
• Legacy retirement: srr_fun %in% c(1, 3, 5) and srr_indices are errors from 5.32.0: express the environmental effect as a recruitment linkage. M1_indices and M1_model %in% c(4, 5) remain deprecated. See the migration table in NEWS.md.

See also

• ?linkage_spec, ?build_growth, ?build_M1, ?build_srr, ?build_catchability, ?build_selectivity, ?build_composition, ?prior_normal
• vignette("model-parameterizations", package = "Rceattle") for the full set of process options.