Choosing Growth Options

Fish of the same age are rarely the same length, and an assessment that fits length data has to state how much they vary and how that spread changes with age. SPoRC holds this in the size-age transition matrix, which gives the probability that a fish of age $a$ falls in length bin $l$ . It is used in four places: length compositions are compared against it, length-based selectivity is applied through it, conditional age-at-length rows are read from it, and when weight at age is derived rather than supplied it converts numbers to biomass, and so enters spawning biomass and the reference points.

All of the options that build it are arguments to Setup_Mod_Biologicals(). In the following, we describe what each option does and the circumstances under which we would use it, working from the choice between fixed and estimated growth, through the curve itself, the spread around it, the plus group and the treatment of seasons, and then on to growth that changes over time. The algebra is in vignette("c_model_equations") and the full argument list is in vignette("t_model_options"). Two case studies run the module end to end: vignette("ad_goa_rex_sole_case_study") for constant parametric growth with conditional age-at-length, and vignette("ae_ebs_pacific_cod_case_study") for time-varying growth kept cohort by cohort.

The examples use the packaged EBS Pacific cod inputs, which have length bins, a weight-length relationship and length compositions, so every option below can be switched on against them.

library(SPoRC)
data("sgl_rg_ebs_pcod_data")

dat <- sgl_rg_ebs_pcod_data
yrs <- dat$years
n_yrs <- length(yrs)
ages <- dat$ages
n_ages <- length(ages)
n_reg <- 1
n_sex <- 1

# lens are the bin midpoints the population is kept on; the transition is
# built on the lower edges, which is what growth_len_lower supplies
input_list <- Setup_Mod_Dim(
  years = yrs, ages = ages, lens = dat$lens,
  n_regions = n_reg, n_sexes = n_sex, n_fish_fleets = 1, n_srv_fleets = 1,
  n_seas = 1, n_pop = 1, natal_region = 1, verbose = FALSE
)

# Setup_Mod_Rec() has to run before the biologicals, because the maturity check
# reads rec_lag
input_list <- Setup_Mod_Rec(
  input_list = input_list, rec_model = "mean_rec", rec_lag = 0, t_spawn = 0,
  sigmaR_spec = "fix", ln_sigmaR = array(log(dat$rec$sigmaR), dim = c(2, 1, n_reg)),
  sigmaR_switch = 1, do_rec_bias_ramp = 0,
  init_age_strc = 2, equil_init_age_strc = 2, ln_global_R0 = dat$mle$ln_R0
)

Notation

Every symbol used below. Growth parameters additionally have population, region and sex subscripts $p, r, s$ , which are dropped throughout because sharing across them is a separate choice (growth_spec).

Symbol	Meaning	Kind
$a$	integer age class, $a = a_{1}, \dots, a_{+}$	index
$a_{+}$	accumulator age, holding every fish that old and older	index
$a^{\ast}$	first age past $A_{1}$ , where cohort propagation starts	index
$y$	model year	index
$\tau$	season, $\tau = 1, \dots, n_{\tau}$	index
$l$	length bin, $l = 1, \dots, n_{l}$	index
$d_{\tau}$	duration of season $\tau$ as a fraction of the year, $\sum_{\tau} d_{\tau} = 1$ (`seasdur`)	data
$t$	fraction of a season elapsed at a reading point (`t_fish`, `t_srv`, `t_spawn`)	data
$e$	fraction of the year elapsed at that reading point	derived
$x$	real age, $x = a + e$	derived
$\ell_{l}$	lower edge of length bin $l$ (`growth_len_lower`)	data
$\tilde{\ell}_{l}$	midpoint of length bin $l$ (`lens`)	data
$A_{1}, A_{2}$	reference ages for $L_{1}$ and $L_{2}$ (`growth_A1`, `growth_A2`)	data
$A_{2}^{cv}$	age at and above which $CV_{2}$ applies	data
$L_{0}$	length at age zero, anchoring the linear phase (`growth_L0`)	data
$L_{1}, L_{2}$	mean length at $A_{1}$ and at $A_{2}$	estimated
$K$	rate of approach to the asymptote	estimated
$\rho$	Richards coefficient; $\rho = 1$ under `"vb_schnute"`	estimated
$CV_{1}, CV_{2}$	spread parameters at the two reference ages	estimated
$L_{\infty}$	asymptotic length, solved from $L_{1}$ , $L_{2}$ , $K$	derived
$\bar{L}(x)$	mean length at real age $x$	derived
$\bar{L}_{+}$	mean length of the plus group	derived
$c(x)$	interpolated spread parameter at $x$	derived
$\sigma_{L}(x)$	standard deviation of length at $x$	derived
$g_{y}(L)$	growth increment over one year under year $y$ ’s parameters	derived
$A_{l,a}$	size-age transition, $P(\ell_{l} \mid a)$	derived
$N_{y,a}$	numbers at age at the start of year $y$	derived
$W_{a}$	weight at age	derived
$\delta_{y}$	time-varying deviation on a growth parameter in year $y$	estimated
$\varepsilon_{y,a}$	semi-parametric deviation on log mean length at age	estimated
$P, P^{lo}, P^{hi}$	a growth parameter and its bounds (`growth_par_bounds`)	see above
$\alpha, \beta$	weight-length parameters in $W = \alpha L^{\beta}$ (`wt_len_pars`)	data
$\Phi$	standard normal distribution function	function

Data or parameters

Under growth_model = "none", the default, SizeAgeTrans is data: an array [p x r x y x seas x l x a x s] of $P(\ell \mid a)$ built outside the model, column-stochastic, one slice per year and season. Weight at age is data too. Nothing about size is estimated.

Under growth_model = "vb_schnute" or "richards", the size-age transition is a function of five or six estimable parameters. SizeAgeTrans is then ignored and may be NA. One is built per fleet at that fleet’s timing (t_fish, t_srv), plus one at spawning time.

We would generally keep "none" when growth has been determined externally and is not itself under study, when bridging an assessment that fixes it, or when there are no length data to inform a curve.

An estimated curve is preferable when length compositions or conditional age-at-length are fit, when size at age may have changed over the series, when the uncertainty in growth should propagate into the biomass, or when selectivity is length-based and would otherwise absorb a misspecified transition, although this comes at the cost of additional parameters and a longer tape build.

Prerequisites

The setup function will stop on each of these:

Requirement	Where
Length bin midpoints	`lens` in `Setup_Mod_Dim()`
`fit_lengths = 1`	`Setup_Mod_Biologicals()`
`growth_len_lower`	lower edges of those same bins, one per bin
`growth_A1`, `growth_A2`	the two reference ages the curve is anchored at

Which data you fit matters more than any switch below. Marginal length compositions identify the curve weakly, because a length distribution is a mixture over ages and abundance at age has to be separated from size at age. Conditional age-at-length constrains growth much more directly: a CAAL row is the age composition of the otoliths taken from one length bin, so it is informative about size at age and has little information about abundance. If you have aged otoliths with their lengths, fit them as CAAL (ObsSrv_caal, do_caal = 1) rather than as marginal age compositions.

The curve and its anchors

MatAA <- array(0, dim = c(1, n_reg, n_yrs, 1, n_ages, n_sex))
MatAA[1, 1, , 1, , 1] <- dat$MatAA

input_list <- Setup_Mod_Biologicals(
  input_list = input_list,
  WAA = NULL, MatAA = MatAA, AgeingError = dat$AgeingError,
  M_spec = "fix", Fixed_natmort = array(0.3866, dim = c(1, n_reg, n_yrs, n_ages, n_sex)),
  fit_lengths = 1, SizeAgeTrans = NA,
  growth_model = "vb_schnute",
  ln_growth_pars = array(log(c(L1 = 14, L2 = 100, K = 0.12, CV1 = 0.15, CV2 = 0.08)),
                         dim = c(1, n_reg, n_sex, 5)),
  growth_A1 = 1.5, growth_A2 = 20,
  growth_len_lower = dat$lens_lower, growth_L0 = dat$lens_lower[1],
  waa_model = "wt_len", wt_len_pars = dat$wtlen
)

WAA = NULL and SizeAgeTrans = NA because both come out of the growth module. The later sections change one argument of this call at a time; ... in them stands for the rest of it unchanged.

The Schnute curve and three of the switches that shape it. Top left: mean length at age, with the linear phase below growth_A1 shaded and a band one standard deviation either side of the mean. Top right: the Richards coefficient changes the shape of the approach to the asymptote, and equals von Bertalanffy at one. Bottom left: the same $L_{2}$ of 70 read as the length at age 8, from which the asymptote is solved, or read as the asymptote itself. Bottom right: the plus-group adjustment adds 4.7 cm at an accumulator age of 8.

Growth is parameterized in Schnute form rather than as $L_{\infty}$ , $K$ , $t_{0}$ . The parameters are lengths at named ages, so starting values can be read off a length-at-age scatter and bounds are lengths. Writing $x$ for real age (integer age plus the fraction of the year elapsed), and dropping the population, region and sex subscripts on every parameter:

$\bar{L}(x) = \begin{cases} L_{0} + \dfrac{x}{A_{1}}\left( L_{1} - L_{0} \right) & x < A_{1} \\[2ex] L_{\infty} + \left( L_{1} - L_{\infty} \right)e^{-K(x - A_{1})} & x \geq A_{1} \end{cases}$

with the asymptote solved from the two anchors,

$L_{\infty} = L_{1} + \frac{L_{2} - L_{1}}{1 - e^{-K(A_{2} - A_{1})}}$

and reported as Linf. The $L_{\infty}$ , $K$ , $t_{0}$ form is strongly correlated between $L_{\infty}$ and $K$ ; the Schnute form is less so, because each parameter is tied to an age with data behind it.

growth_model = "richards" applies the same curve to lengths raised to a power $\rho$ :

$\bar{L}(x)^{\rho} = L_{\infty}^{\rho} + \left( L_{1}^{\rho} - L_{\infty}^{\rho} \right)e^{-K(x - A_{1})}$

with $L_{\infty}^{\rho} = L_{1}^{\rho} + \left( L_{2}^{\rho} - L_{1}^{\rho} \right) / \left( 1 - e^{-K(A_{2} - A_{1})} \right)$ . At $\rho = 1$ this is the von Bertalanffy curve exactly, so the Richards form nests it. The linear phase below $A_{1}$ is unchanged.

growth_A1 and growth_A2 are data rather than parameters. We would set growth_A1 at or just above the youngest age present in the length data, since below it the curve is a straight line from growth_L0 and is therefore an assumption rather than a fit, and set growth_A2 old enough that fish are near the asymptote while young enough that the sample at that age is not small.

growth_A2 = "Linf" reads $L_{2}$ as the asymptote itself, with no second reference age from which to solve it, and we would use it when matching a source parameterized in that way. The bottom-left panel shows that the same $L_{2}$ gives a different curve under each reading.

growth_L0 anchors the linear phase and defaults to the lower edge of the first length bin. It is rarely estimable, and matters only when recruits appear in the length compositions.

$\rho$ gives a more flexible approach to the asymptote, at the cost of one parameter that is often poorly identified. We would start at $\rho = 1$ , and if the estimate barely moves from it while its standard error remains large, revert to "vb_schnute".

The spread around the mean

The curve described above gives only the mean length at age, whereas the size-age transition also requires the spread of lengths around that mean, which is set by the three switches described in this section.

Top left: the CV runs from CV1 to CV2 between the reference ages, interpolated on mean length or on age. Top right: the same two parameters read as CVs, which scale the mean, or as standard deviations, which do not. Bottom left: the two dimensionally coherent pairings at a spread of 0.35, normal with growth_sd_type = "cv" and lognormal with "sd"; note the mass the normal puts in the smallest bin at age 1. Bottom right: the resulting size-age transition, on a square root color scale.

Two parameters hold the spread, $CV_{1}$ and $CV_{2}$ at the two reference ages, interpolated between them. growth_cv_type sets what the interpolation runs on:

$c(x) = \begin{cases} CV_{1} & x < A_{1} \\[1ex] CV_{1} + \left( CV_{2} - CV_{1} \right)\dfrac{\bar{L}(x) - L_{1}}{L_{2} - L_{1}} & A_{1} \leq x < A_{2}^{cv}, \quad \texttt{"len"} \\[2ex] CV_{1} + \left( CV_{2} - CV_{1} \right)\dfrac{x - A_{1}}{A_{2}^{cv} - A_{1}} & A_{1} \leq x < A_{2}^{cv}, \quad \texttt{"age"} \\[2ex] CV_{2} & x \geq A_{2}^{cv} \end{cases}$

where $A_{2}^{cv} = A_{2}$ , or the oldest age evaluated when growth_A2 = "Linf". The branch is taken on age either way; only the middle branch’s value depends on $\bar{L}(x)$ . The two agree at the reference ages and differ most where growth is fast, because interpolating on length compresses the change into the young ages. "len" is the Stock Synthesis convention; "age" gives a CV that declines steadily with age.

growth_sd_type then turns $c(x)$ into the standard deviation the transition is built from:

$\sigma_{L}(x) = \begin{cases} c(x)\,\bar{L}(x) & \texttt{"cv"} \\[1ex] c(x) & \texttt{"sd"} \end{cases}$

Under "cv" the spread scales with size, so that a 20 cm and a 60 cm fish have standard deviations in the ratio 1:3, whereas under "sd" the absolute spread is the same at every age. The "cv" form gives an implausibly tight distribution at the youngest ages when the linear phase runs down to a very small growth_L0, and we would use "sd" in that case. A CV that declines with size ( $CV_{2} < CV_{1}$ ) is the usual estimate, since the lengths of young fish vary with birth date and early growth while old fish have converged on the asymptote.

growth_dist sets the distribution that is binned. The two are close at the spreads most assessments estimate, and the panel above uses 0.35 in order to make the difference visible. We would use "lognormal" when the spread is large enough that a normal distribution places appreciable mass below zero length.

Note that these two switches have to be paired carefully, because the transition is always built from the $\sigma_{L}(x)$ that growth_sd_type produced and the lognormal transform then reads that value as a log-scale standard deviation. For this reason "lognormal" belongs with growth_sd_type = "sd", since combining it with "cv" would pass a length-scale spread into a log-scale transform.

The plus group

The accumulator age holds every fish that old and older, so its mean size is not the curve’s value there.

growth_plus_group = "mixture", the default, treats it as a mixture of the ages it contains. The share of fish $k$ years past the accumulator age $a_{+}$ is taken as proportional to $e^{-0.2k}$ , the survivorship under a total mortality of 0.2 per year, and mean length is taken to rise linearly from $\bar{L}(a_{+})$ to $L_{\infty}$ over a second lifetime of $a_{+}$ years:

$\bar{L}_{+} = \frac{\sum_{k=0}^{a_{+}} e^{-0.2k}\left[ \bar{L}(a_{+}) + \dfrac{k}{a_{+}}\left( L_{\infty} - \bar{L}(a_{+}) \right) \right]}{\sum_{k=0}^{a_{+}} e^{-0.2k}}$

The 0.2 is a fixed assumption about how fast numbers decline with age, not an estimated rate, and it sets how much weight the older and larger fish have. $\bar{L}_{+}$ lies between $\bar{L}(a_{+})$ and $L_{\infty}$ . growth_plus_group = "curve" uses $\bar{L}(a_{+})$ directly.

The two agree closely when the accumulator age is old enough that little growth remains, and diverge as the plus group moves younger. With an accumulator age well short of the asymptote, "curve" undersizes the plus group and shifts its length composition left.

Binning

The mean and the spread are turned into the size-age transition by integrating the distribution over each length bin. With $\ell_{l}$ the lower edge of bin $l$ , each edge is standardized,

$z_{l,a} = \frac{\ell_{l} - \bar{L}(a)}{\sigma_{L}(a)}$

and each bin takes the probability between its edges, with the two tails accumulated into the end bins:

$A_{l,a} = \begin{cases} \Phi(z_{2,a}) & l = 1 \\[1ex] \Phi(z_{l+1,a}) - \Phi(z_{l,a}) & 1 < l < n_{l} \\[1ex] 1 - \Phi(z_{n_{l},a}) & l = n_{l} \end{cases}$

The column sums to one for any $\bar{L}(a)$ and $\sigma_{L}(a)$ , including means outside the binned range entirely, which is what makes each column a proper $P(\ell \mid a)$ . Under growth_dist = "lognormal" the same differences are taken on the log scale, about a median-corrected mean so that $E[L] = \bar{L}(a)$ :

$z_{l,a} = \frac{\log \ell_{l} - \left( \log \bar{L}(a) - \sigma_{L}^{2}(a)/2 \right)}{\sigma_{L}(a)}$

The value of $\sigma_{L}(a)$ therefore sets how many bins an age’s mass spreads over, in that large values flatten the column while small values concentrate it near $\bar{L}(a)$ .

Timing within the year

Growth parameters are annual. growth_pars_y has no season dimension, and none of $L_{1}, L_{2}, K, CV_{1}, CV_{2}, \rho$ changes within a year. What changes within the year is the age the curve is evaluated at.

Where each reading point falls

Each fleet reads its own size-age transition and weight at age at the point its observations are taken: fishery fleet $f$ at $t^{fish}_{\tau,f}$ , survey $g$ at $t^{srv}_{\tau,g}$ , and the spawning weight at $t^{spawn}$ in the spawning season. Each $t$ is a fraction of its own season, in $[0, 1]$ , not a fraction of the year.

Converting one to a fraction of the year gives the elapsed time $e$ . For a reading at fraction $t$ of season $\tau$ :

$e(\tau, t) \;=\; \underbrace{\sum_{\tau' = 1}^{\tau - 1} d_{\tau'}}_{\text{seasons already completed}} \;+\; \underbrace{t\, d_{\tau}}_{\text{part of season } \tau \text{ elapsed}}$

where $d_{\tau'}$ is the duration of season $\tau'$ as a fraction of the year, from seasdur, and the seasons sum to a year, $\sum_{\tau = 1}^{n_{\tau}} d_{\tau} = 1$ . The empty sum is zero, so $e(1, 0) = 0$ is the start of the year.

The curve is then evaluated at the real age

$x_{a} = a + e(\tau, t)$

for every age class $a$ , giving $\bar{L}(x_{a})$ and $\sigma_{L}(x_{a})$ from the same two equations as before, and nothing is interpolated between integer ages.

A worked example, four equal quarters so $d_{\tau} = 0.25$ for all $\tau$ , with a survey read at the midpoint of the third quarter, $\tau = 3$ and $t = 0.5$ :

$e = (0.25 + 0.25) + 0.5 \times 0.25 = 0.625$

so an age 4 fish is placed at $x = 4.625$ on the curve. Note that $t = 0.5$ puts the reading halfway through that quarter, not halfway through the year. With a single season, $d_{1} = 1$ , the two coincide and $t = 0.5$ does give $x = a + 0.5$ .

Why seasons add no approximation

Write the one-year increment under year $y$ ’s parameters as

$g_{y}(L) = \left[ L_{\infty,y}^{\rho_{y}} + \left( L^{\rho_{y}} - L_{\infty,y}^{\rho_{y}} \right)e^{-K_{y}} \right]^{1/\rho_{y}}$

and the same map over an arbitrary elapsed time $e$ as $g_{y}(L; e)$ , which replaces $e^{-K_{y}}$ with $e^{-K_{y}e}$ . Applied to a length already on the curve it returns the curve’s own value $e$ later:

$g_{y}\!\left( \bar{L}(x);\, e \right) = \bar{L}(x + e)$

The identity is exact for the von Bertalanffy form and holds to machine precision for the Richards form (a maximum absolute difference of $3 \times 10^{-14}$ cm at $\rho = 0.4$ and $\rho = 2.5$ ). Because the module evaluates $\bar{L}(x_{a})$ directly rather than stepping season by season, the season structure cannot introduce error: with growth constant over years, the same mean length at age comes out whatever $n_{\tau}$ and $d_{\tau}$ are.

Only distinct values of $e$ within a season are evaluated, so that fleets which share a timing cost nothing extra.

The two exceptions

Two quantities are advanced by $g_{y}(\cdot\,; e)$ rather than read off the curve, because neither sits on it:

$\bar{L}_{\tau,a} = \begin{cases} g_{y}\!\left( \bar{L}_{+};\, e \right) & a = a_{+} \\[1.5ex] g_{y}\!\left( \bar{L}_{y,a};\, e \right) & a \geq a^{\ast}, \quad \texttt{growth\_tv\_type = "cohort"} \\[1.5ex] \bar{L}(x_{a}) & \text{otherwise} \end{cases}$

The plus group starts from its mixture-adjusted size $\bar{L}_{+}$ , which is a weighted average of ages rather than a point on the curve. Under cohort growth, every age at or past $a^{\ast}$ starts from $\bar{L}_{y,a}$ , the size that cohort actually reached.

# one curve per sex, shared across regions, with the second CV kept fixed
input_list <- Setup_Mod_Biologicals(
  input_list = input_list, ...,
  growth_model = "vb_schnute",
  growth_spec = "est_shared_r",
  growth_fix = c(L1 = FALSE, L2 = FALSE, K = FALSE, CV1 = FALSE, CV2 = TRUE)
)

growth_spec uses the same vocabulary as every other _spec argument: "est_all" (one set per population, region and sex), "est_shared_r", "est_shared_s", "est_shared_r_s", or "fix". growth_fix overrides it for individual parameters, holding named ones at their starting values regardless of the sharing.

Sex: if the sexes differ in size at age and the compositions are sex-split, we would estimate growth separately for each, because sharing it across sexes that differ leaves the difference to be absorbed by selectivity instead.

Region: we would generally share growth across regions unless regional length data are available to identify a difference between them.

CV2 is usually the first parameter we would fix, since it is informed only by the oldest ages, where samples tend to be small and the plus group confounds it, and rho is generally best held at one unless the estimate moves away from it.

Starting values matter more here than in most parts of the model. ln_growth_pars defaults to the ends of the length bins with a rate of 0.15 and CVs of 0.1, which is a placeholder rather than a prior. Supply your own on the log scale, dimensioned [n_pop x n_regions x n_sexes x n_growth_pars] in the order L1, L2, K, CV1, CV2 and then rho. An external growth fit, even a rough one, is a reasonable starting point.

Growth that changes over time

Growth can be allowed to change over time through two mechanisms, described in turn below, and the two may be used together.

Varying the parameters

input_list <- Setup_Mod_Biologicals(
  input_list = input_list, ...,
  growth_model = "richards",
  growth_tv_model = c(L1 = "iid", K = "iid"),
  growth_tv_years = list(L1 = 2000:2024, K = 2000:2024),
  growth_tv_link = "logit",
  growth_par_bounds = dat$growth$bounds,
  growth_tv_sigma_spec = "fix",
  growth_tv_type = "cohort"
)

growth_tv_model names a structure per parameter, "none", "iid" or "rw", supplied either as a vector in parameter order or named with the rest left constant. Each varying parameter gains a deviation series in ln_growth_devs and a process error standard deviation. Growth parameters trade off against each other, so that letting all five vary at once tends to give a well-fitting series which is difficult to interpret, and we would name only those parameters there is a specific reason to move.

growth_tv_years restricts each parameter’s deviations to given calendar years, supplied either as a vector applied to all of them or as a list named per parameter, with deviations outside that range kept at zero. We would generally use the years over which the length data are informative.

growth_tv_link sets the scale the deviation enters on. For a base parameter $P$ with deviation $\delta_{y}$ , the value in year $y$ is

$P_{y} = \begin{cases} P\exp(\delta_{y}) & \texttt{"log"} \\[1ex] P^{lo} + \left( P^{hi} - P^{lo} \right)\,\mathrm{logit}^{-1}\left[ \mathrm{logit}\!\left( \dfrac{P - P^{lo}}{P^{hi} - P^{lo}} \right) + \delta_{y} \right] & \texttt{"logit"} \end{cases}$

The log link is unbounded and can drive a CV negative or a $K$ to an implausible value. The logit link keeps the parameter strictly inside $(P^{lo}, P^{hi})$ from growth_par_bounds however large the deviation. Growth parameters are all bounded by construction, so "logit" with sensible bounds is the better default. The realized $P_{y}$ are reported as growth_pars_y.

growth_tv_sigma_spec holds the process error standard deviations at their starting values ("fix") or estimates them ("est"). The starting values are taken from the first slot of growth_pe_pars, one entry per growth parameter. Fixing them makes the deviations a penalized-likelihood problem in which you set how much variation is allowed, whereas estimating them requires the deviations to be integrated out, as described below.

growth_tv_spec shares the deviation series across regions and sexes with the same vocabulary as growth_spec. A deviation series per stratum adds parameters quickly, so this is often worth sharing even when the base parameters are not. growth_rw_init_sigma gives the first year of a random walk its own standard deviation, defaulting to 5, which leaves the starting point close to free.

Curve or cohort

Top: the two growth parameters that vary in the EBS Pacific cod model, as the assessment’s estimate realizes them year by year. Bottom: what the same deviations do to mean length at age under each reading. The assessment’s deviations are small enough that the two nearly coincide, so the bottom panel exaggerates them into one sustained drop in $K$ from 2000 to make the mechanism visible. It is an illustration, not the assessment’s estimate.

growth_tv_type = "curve", the default, reads every year’s sizes off that year’s curve, so that nothing is kept between years: if $K$ drops, every age is immediately shorter, including fish that grew through twenty prior years. We would use this option when the deviations are thought to represent something acting on all ages at once, or when the curve is better understood as describing the population sampled in a given year rather than a growth process.

"cohort" propagates size at age forward instead. Writing $g_{y}$ for the growth increment over one year under year $y$ ’s parameters,

$g_{y}(L) = \left[ L_{\infty,y}^{\rho_{y}} + \left( L^{\rho_{y}} - L_{\infty,y}^{\rho_{y}} \right)e^{-K_{y}} \right]^{1/\rho_{y}}$

the start-of-year mean length advances as

$\bar{L}_{y+1,a} = \begin{cases} L_{0} + \dfrac{a}{A_{1}}\left( L_{1,\,y+1-a} - L_{0} \right) & a < A_{1} \\[2ex] \bar{L}_{y+1}(a) & a = a^{\ast} \\[2ex] g_{y}\left( \bar{L}_{y,a-1} \right) & a^{\ast} < a < a_{+} \\[2ex] \dfrac{\left( N_{y,a_{+}-1} + 0.01 \right)g_{y}\left( \bar{L}_{y,a_{+}-1} \right) + \left( N_{y,a_{+}} + 0.01 \right)g_{y}\left( \bar{L}_{y,a_{+}} \right)}{N_{y,a_{+}-1} + N_{y,a_{+}} + 0.02} & a = a_{+} \end{cases}$

with $a^{\ast}$ the first age past $A_{1}$ . Three consequences follow from this. Ages still in the linear phase take the length at $A_{1}$ that their birth year’s parameters gave them, so that a cohort born in a poor year keeps that smaller start for the rest of its life. The plus group blends the cohort entering it with the fish already there by their numbers at age, which is the only place in the model where growth depends on abundance. And a poor growth year tracks the affected cohorts up through the ages for the remainder of the series, which is why we would use this option when the deviations are thought to represent conditions a fish lived through.

Two further consequences follow. Because the plus group reads abundance, growth has to be evaluated inside the population dynamics year loop, and the model therefore costs more to tape. In addition, $c(x)$ is evaluated once from the first year’s curve and then kept while the mean moves, so that the deviations change mean size without also changing the spread.

Semi-parametric growth

growth_semipar adds a year-by-age surface of deviations multiplying the parametric mean length at age,

$\bar{L}_{y,a}^{\ast} = \bar{L}_{y,a}\exp\left( \varepsilon_{y,a} \right)$

applied after the curve and after any cohort propagation. The curve remains the parametric part and $\varepsilon$ holds the deviations from it. This represents changes in the shape of size at age that no parameter of the curve can express: a curve with five parameters cannot fit a year in which only the four-year-olds were small, and an unconstrained transition per year is not identified.

input_list <- Setup_Mod_Biologicals(
  input_list = input_list, ...,
  growth_model = "vb_schnute",
  growth_semipar = "2dar1",
  growth_semipar_spec = "est",
  growth_semipar_ages = 2:12,
  growth_semipar_years = 1990:2024
)

Top: a known year-by-age surface used to simulate length compositions and conditional age-at-length. Bottom: mean length at age recovered from those data by a model with a 2D AR(1) surface and by one with the parametric curve alone, against the truth. Over the whole surface the semi-parametric model is within 8.9% of the simulated mean length at age; the parametric curve alone is within 15.8%.

The available forms are the process errors semi-parametric selectivity uses: "iid", "rw" (a random walk over years within an age), "2dar1" (a separable first-order autoregression over ages and years), and "3dmarg" / "3dcond" (a three-dimensional Gaussian Markov random field over age, year and cohort, on the marginal or conditional variance).

The correlated forms are what make the surface estimable, because a year-by-age deviation is not identified from length data one cell at a time and the process error therefore has to borrow strength across neighboring cells. "iid" provides none of that correlation and is mainly useful as a null model against which the others can be compared, whereas "2dar1" is generally a sensible default. The three-dimensional forms add a hyperparameter and correlate along the cohort diagonal as well as over ages and years, which is likely to be worthwhile when a cohort signal is expected.

growth_semipar_ages and growth_semipar_years hold deviations outside them at zero, and we would restrict the surface to the ages and years that the length data actually cover, since elsewhere it adds parameters that nothing constrains.

The spread follows the deviated mean, so that under growth_cv_type = "len" a deviation which lengthens a fish also moves it along the CV ramp, whereas under "age" the spread at age is unaffected.

The surface and the parametric time variation compose, and we would generally use the parameters for changes that can be attributed to a named mechanism and the surface for what remains, while avoiding asking both to explain the same variation.

Why the two sets of deviations are separate

They are indexed over different dimensions and applied at different points:

ln_growth_devs         # [n_pop, n_regions, n_yrs, n_growth_pars, n_sexes]  4th dim: parameter
ln_growth_semipar_devs # [n_pop, n_regions, n_yrs, n_ages,  n_sexes]  4th dim: age

ln_growth_devs is applied before the curve is evaluated, perturbing L1, L2, K, CV1, CV2 and rho. ln_growth_semipar_devs multiplies mean length at age after the curve and after any cohort propagation.

Three consequences follow from this. The correlated process errors correlate over age, year and cohort, which is only meaningful when one of the indices is age, so that "2dar1" and the three-dimensional forms are available to the surface but not to the parametric deviations, for which "iid" and "rw" are the only sensible structures. Cohort propagation works off the growth increment that the year’s parameters imply and therefore acts on the parametric deviations, whereas the surface multiplies in afterwards and is not propagated forward. Finally, parametric deviations keep size at age within the shape the curve can express, whereas the surface allows it to leave that shape entirely.

The two sets of deviations do share one array, growth_pe_pars, in which slot [,,,,1] holds the time-varying standard deviations and slot [,,,,2] the surface’s hyperparameters, with the slots that a given form does not read mapped off.

Deviations as random effects

Both sets of deviations can be integrated out with the Laplace approximation rather than penalized, as the selectivity deviations are in vignette("n_single_region_ebs_pollock_randomeff_case_study"). Pass the parameter name to fit_model():

# time-varying growth parameters, with their process error estimated
fit <- fit_model(input_list$data, input_list$par, input_list$map,
                 random = "ln_growth_devs", newton_loops = 3, silent = TRUE)

# a semi-parametric surface, or both together
fit <- fit_model(input_list$data, input_list$par, input_list$map,
                 random = c("ln_growth_devs", "ln_growth_semipar_devs"),
                 newton_loops = 3, silent = TRUE)

This is what growth_tv_sigma_spec = "est" and growth_semipar_spec = "est" require, because maximizing over a process error standard deviation and the deviations it penalizes at the same time, with both treated as fixed effects, drives that standard deviation toward zero. Integrating the deviations out instead makes it a variance parameter that the data can inform. When the process error is held fixed, the fixed-effect treatment is adequate and is cheaper to fit.

fit_model() uses the tape Hessian for its Newton refinement only when there are no random effects, so a random-effects fit takes the slower path.

Weight at age

input_list <- Setup_Mod_Biologicals(
  input_list = input_list,
  WAA = NULL, MatAA = MatAA, ...,
  growth_model = "richards",
  waa_model = "wt_len", wt_len_pars = dat$wtlen  # a, b in W = a L^b
)

waa_model = "data", the default, reads WAA, WAA_fish and WAA_srv from the arguments of the same name.

waa_model = "wt_len" builds them from the size-age transition and the weight-length relationship $W = \alpha L^{\beta}$ applied at the bin midpoints $\tilde{\ell}_{l}$ :

$W_{a} = \sum_{l} A_{l,a}\,\alpha\,\tilde{\ell}_{l}^{\beta}$

Weight at age then reflects the spread of length at age rather than being the weight of the mean length, which differ when $\beta \neq 1$ . The spawning weight uses the transition at spawning time and each fleet’s weight the transition at that fleet’s timing, so weight at age becomes fleet-specific, and time-varying whenever growth is. WAA may be NULL.

We would use "wt_len" when growth is estimated and the biomass should inherit its uncertainty and time variation, and "data" when empirical weight at age is available, which is generally preferable to a weight-length relationship, or when matching an assessment that supplies it.

Under "wt_len" the weight arrays are model output, but reference point and projection code read data$WAA. Copy the reported arrays into the data list before calling Get_Reference_Points() or Do_Population_Projection():

rep <- fit$rep
input_list$data$WAA <- rep$WAA
input_list$data$WAA_fish <- rep$WAA_fish
input_list$data$WAA_srv <- rep$WAA_srv

A related switch sits with the fleets rather than here: fish_waa_selected and srv_waa_selected put catch and index biomass on the selection-weighted weight at age, which applies only with length-based selectivity and a derived weight.

Compositions on coarser bins

LenBinMap maps the model’s length bins onto the bins the compositions are recorded on, a matrix [n_lens x n_obs_lens] whose rows sum to one. Growing fish on fine bins and fitting them on coarse ones gives a smoother transition and better-behaved tails while keeping the compositions on the bins they were collected on. The expected compositions are mapped through it inside the likelihood, as the ageing error matrix maps ages.

Checking the setup

Setup_Mod_Biologicals() prints what it built when verbose = TRUE. The EBS Pacific cod setup, with a semi-parametric surface added on top, reports:

Growth is estimated (Richards); SizeAgeTrans is built inside the model
Growth parameters varying over time: L1 (iid), K (iid); link logit; size at age
  read from cohort propagation from 2000
Weight at age (spawning, fishery, survey) is derived from growth and the
  weight-length relationship
Semi-parametric growth: 2dar1 deviations on mean length at age over 48 years
  and 11 ages, process error fix
Length compositions are recorded on 24 bins, mapped from the model's 121 bins
  inside the likelihood

This reports which parameters vary and under what structure, the link they enter on, the year from which cohort propagation starts, how many ages the surface covers, and where weight at age comes from.

After fitting, the report contains:

Field	What it holds
`growth_pars_y`	the parameters in effect, `[p x r x y x 6 x s]`
`Linf`	the derived asymptote, `[p x r x y x s]`
`L_beg`	mean length at age at the start of the year
`mean_LAA_fish`, `mean_LAA_srv`, `mean_LAA_spawn`	mean length at age at each fleet’s timing and at spawning
`sd_LAA_*`	the matching spreads
`SizeAgeTrans_fish`, `SizeAgeTrans_srv`, `SizeAgeTrans_spawn`	the size-age transitions
`growth_tv_nLL`, `growth_semipar_nLL`	the deviation penalties

Plot growth_pars_y against year before interpreting a time-varying fit, and mean_LAA_* in a few years against any external length-at-age data. A growth series that swings without an identifiable mechanism is usually absorbing something else, most often selectivity or a change in sampling.

Summary

Question	Argument	Default and when to move off it
Estimate growth at all?	`growth_model`	`"none"`; move when fitting length or CAAL data
Which curve?	`growth_model`	`"vb_schnute"`; `"richards"` only if $\rho$ moves
Where to anchor?	`growth_A1`, `growth_A2`	no default; ages the data cover
Is $L_{2}$ the asymptote?	`growth_A2 = "Linf"`	no; use when matching a source that says so
How does the CV interpolate?	`growth_cv_type`	`"len"`; `"age"` for a steadier decline
CVs or SDs?	`growth_sd_type`	`"cv"`; `"sd"` if young ages get implausibly tight
Length distribution?	`growth_dist`	`"normal"`; `"lognormal"` at large spread, paired with `growth_sd_type = "sd"`
Plus group size?	`growth_plus_group`	`"mixture"`; `"curve"` only near the asymptote
Share across strata?	`growth_spec`, `growth_fix`	`"est_all"`; share regions, split sexes, fix `CV2`
Does growth move?	`growth_tv_model`	`NULL`; name only defensible parameters
Bounded deviations?	`growth_tv_link`	`"log"`; prefer `"logit"` with bounds
Memory of past years?	`growth_tv_type`	`"curve"`; `"cohort"` for lived conditions
Shape changes too?	`growth_semipar`	`"none"`; `"2dar1"` when the curve cannot express it
Penalized or integrated?	`random` in `fit_model()`	penalized; integrate to estimate a process error
Where does weight come from?	`waa_model`	`"data"`; `"wt_len"` to inherit growth’s uncertainty