8. Model parameterizations

In progress!

Population dynamics

CEATTLE is an age-structured model that follows cohorts of a species $i$, age $a$, and sex $s$ through time $y$. Annual recruitment and the initial age structure are estimated, either as random effects with the recruitment variance estimated, or as penalized deviates with it fixed. Because the initial age structure is estimated assuming the population was already depleted before the first model year, the mean recruitment parameter $\bar{R}_{i}$ is biased low: projected recruitment should instead resemble the hindcast recruitment deviations, and biomass reference points should use the mean of annual recruitment. Populations are 1- or 2-sex via nsex, with the sex ratio at birth $\rho_{i1}$ taken from the first age of sex_ratio. Natural mortality is partitioned into predation mortality (M2) and residual mortality from all other causes (M1), which is either input on the M1_base sheet or estimated with M1Fun = build_M1(M1_model = ...) in fit_mod(). Setting msmMode = 0 turns predation off (below).

$$N_{is1,y} = R_{isy} = \bar{R}_{i }e^{\tau_{i,y}}*\rho_{i1}$$$$N_{isa,1} = \left \{ \begin{array}{} \bar{R}_ie^{ \left( -j M1_{isa} \right)}*\rho_{i1} \\ \frac{\bar{R}_ie^{ \left( -\sum^{a-1} M1_{isa} \right)}} {\left(1 - e^{\left( -\sum^{a-1} M1_{is A_i} \right)} \right)} *\rho_{i1} \end{array} \right.$$$$N_{is a+1, y+1} = N_{isa,y} e^{-Z_{isa,y}}$$$$N_{is A_i, y+1} = N_{is A_{i}-1,y} e^{-Z_{is A_{i}-1,y}} + N_{is A_{i},y} e^{-Z_{is A_{i},y}}$$$$B_{isa,y} = N_{isa,y} W_{isa,y}$$$$SSB_{isa,y} = B_{isa,y} \rho_{ia}$$$$Z_{isa,y} = M1_{isa} + M2_{isa,y} + F_{isa,y}$$

Predation sub-model

suitMode sets predator-prey preference (suitability, that is predator selectivity) and msmMode the functional response. The predation model is parameterized from the proportion by weight of prey-at-age in a predator-at-age’s diet (diet_data) and from ration/consumption (below).

Suitability

Suitability represents a predator’s prey and size preference, together with factors such as spatial overlap between predator and prey distributions.

MSVPA Type 2 based suitability

Suitability is derived empirically from the mean proportion by weight of prey-at-age in the stomachs of predators-at-age, $\bar{U}_{pbjisa}$ (diet_data). Magnusson (1995) derives the suitability coefficients. suitMode = 0 Predator-prey suitability for predators of species $p$, sex $b$, age $j$, and prey of species $i$, sex $s$ and age $a$$$\hat{S}_{pbjisa} = \frac{1}{n_y} \sum_y \left( \frac{ \frac{\bar{U}_{pbjisa}} {B_{isa,y}} } {\sum_{isa} \left( \frac{\bar{U}_{pbjisa}} {B_{isa,y}} \right) + \frac{1 + \sum_{isa} \bar{U}_{pbjisa}} {B_p^{other}} } \right)$$

Predation mortality

Single-species mode

Multiple independent single-species age-structured models: all mortality is fishing plus input or estimated natural mortality ($M = M1$). msmMode = 0 $$M2_{isa,y} = 0$$

MSVPA Type 2 based predation mortality

msmMode = 1 An MSVPA-based formulation following previous MSCAA models (Curti et al., 2013; Holsman et al., 2016; Tsehaye et al., 2014). Predation mortality is built from suitability coefficients $\hat{S}_{pbjisa}$ (estimated or empirical), input annual ration/consumption (below), mean annual prey biomass whether estimated or input, and the input biomass of other prey available to predator $p$ ($B_p^{other}$). Predation and biomass are interdependent, so predation mortality and the age-structured dynamics are iterated to convergence; three to five iterations suffice. Fits are generally insensitive to $B_p^{other}$ over a range of values, which has historically been derived from ecosystem-wide biomass estimates. For predators of species $p$, sex $b$, age $j$ and prey of species $i$, sex $s$, age $a$$$M2_{isa,y} = \sum_{pbj} \left( \frac{N_{pbj,y} \delta_{pbj,y} {S}_{pbjisa}} {\sum_{isa} \left({S}_{pbjisa} B_{isa,y} \right) + B_p^{other} \left(1 - \sum_{isa} \left({S}_{pbjisa} \right) \right)} \right)$$

MSVPA Type 3 based predation mortality

msmMode = 2 A Type III functional response variant of the MSVPA predation model. The available-food denominator replaces prey biomass $B_{isa,y} = N_{isa,y} W_{isa,y}$ with $N^2_{isa,y} W_{isa,y}$, creating a sigmoidal (accelerating-then-saturating) relationship between prey density and predation rate. This gives low-density prey a refuge from predation and can stabilize multi-species dynamics at low stock sizes. $$M2_{isa,y} = \sum_{pbj} \left( \frac{N_{pbj,y} \delta_{pbj,y} {S}_{pbjisa}} {\sum_{i'a's'} \left({S}_{pbji'a's'} N^2_{i'a's',y} W_{i'a's',y} \right) + B_p^{other} \left(1 - \sum_{i'a's'} \hat{S}_{pbji'a's'} \right)} \right)$$

Consumption

Consumption is either an environmentally driven bioenergetics model or direct input of consumption-at-age, controlled per species on the bioenergetics_control sheet:

• Ceq Which bioenergetics equation scales maximum consumption with temperature: 1 = exponential (Stewart et al 1983); 2 = warm-water temperature dependence (Kitchell et al 1977; sensu Holsman et al 2015); 3 = cool- and cold-water temperature dependence (Thornton and Lessem 1979); 4 = temperature function held at 1.0, which bypasses bioenergetics and reads ration_data as empirical annual ration (set CA = 1, fday = 1, CB = 0)
• Cindex Which env_data index supplies $T_y$ for bioenergetics-based ration
• Pvalue Tuning parameter $\varphi_p$ scaling the maximum consumption used for ration
• fday Foraging days per year, $Fday_p$
• CA $\alpha^{\delta}_p$, the intercept of the allometric mass function for maximum consumption, Cmax = CA * Weight ^ (1 + CB), with Weight in grams
• CB $\beta^{\delta}_p$, the slope of that function
• Qc, Tco, Tcm, Tcl, CK1, CK4 Parameters $Q^c_p$, $T^{co}_p$, $T^{cm}_p$, $T^{cl}_p$, $CK1_p$ and $CK4_p$ of the temperature scaling function; which are read depends on Ceq, per the algorithms below

ration_data holds the relative foraging rate $RFR_{pbj,y}$; it, and so ration, can be tuned through fday and Pvalue. The weight-at-age used in the consumption equations is set by pop_wt_index.

Bioenergetics based ration

Ceq %in% 1:3 Individual specific ration ($\delta_{pbj,y}$$(kg~kg^{-1}~yr^{-1})$) of predator $p$ of sex $b$ and age $j$ in year $y$: $$\delta_{pbj,y} = RFR_{pbj,y}*Fday_p *\alpha^{\delta}_p W_{pbj,y}^{\left(1+\beta^{\delta}_p\right)} *f\left(T_y \right)_p$$

where $W$ is weight-at-age in grams (the model stores it in kg and multiplies by 1000 here) and $\delta$ comes back in kg per predator per year, so $\alpha^{\delta}_p$ is the coefficient for a gram-scale weight. Ration can be tuned by adjusting Pvalue ($\varphi_p$) to help convergence: $$\delta_{pbj,y} = \varphi_p * \delta_{pbj,y}$$

Temperature scaling algorithms.

Ceq = 1: Exponential function from Stewart et al. 1983: $$f \left(T_y \right)_p = e^{Q^c_p *T_y}$$

Ceq = 2: Temperature dependence for warm-water-species from Kitchell et al. 1977: $$f \left(T_y \right)_p = V^X e^{\left(X \left(1-V \right) \right)}$$$$V = \left( T^{cm}_p - T_y \right)/ \left( T^{cm}_p - T^{co}_p \right)$$$$X = \left(Z^2 \left(1+ \left( 1 + 40 /Y \right)^{0.5} \right)^2 \right)/400$$$$Z = ln\left(Q^c_p \right) \left( T^{cm}_p - T^{co}_p \right)$$$$Y = ln\left(Q^c_p \right) \left( T^{cm}_p - T^{co}_p + 2\right)$$

Ceq = 3: Temperature dependence for cool and cold-water species from Thornton and Lessem 1979: $$f \left(T_y \right)_p = Ka * Kb$$$$Ka = C^{K1}_p * L_1 / (1+C^{K1}_p * (L_1 - 1 ))$$$$L_1 = e^{G_1 * (T_y - Q^c_p)}$$$$G_1 = \frac{1}{T^{co}_p-Q^c_p}*\log\left(\frac{0.98*(1-CK1_p)}{CK1_p*0.02}\right)$$$$Kb = \frac{CK4_p *L_2} {1+CK4_p*(L_2-1)}$$$$L_2 = e^{G_2 (T^{cl}_p - T_y)}$$$$G_2 = \frac{1}{T^{cl}_p-T^{cm}_p}*\log\left( \frac{0.98*(1-CK4_p)}{CK4_p*0.02}\right)$$

Input ration-at-age

Ceq = 4 reads consumption-at-sex/age directly from ration_data as $RFR_{pbj,y}$. Set CA = 1, fday = 1 and CB = 0. $$\delta_{pbj,y} = RFR_{pbj,y}* W_{pbj,y}$$

Fishery and survey observation model

CEATTLE is fit to time series of fishery catch, survey biomass, and length- or age-composition data. Log-catch and log-index observations are normally distributed; age and length compositions are multinomial, with their sample sizes input. There is no separate sex-ratio likelihood component: sex-ratio information enters through joint (comp_data$Sex = 3) composition data, whose female and male cells are normalized together. The log-scale index SD is input, estimated as a free parameter, or derived analytically following Walters and Ludwig (1994), per Estimate_index_sd. Selectivity and catchability are estimated or input separately for each fishery and survey.

Catch is estimated in the model following the Baranov catch equation: $$\hat{C}_{isa,y} = \frac{F_{isa,y}} {Z_{isa,y}} \left( 1 - e^{-Z_{isa,y}} \right) N_{isa,y}$$$$\hat{C}_{i,y} = \sum_s\sum_a^{A_i} \frac{F_isa,y} {Z_{isa,y}} \left( 1 - e^{-Z_{isa,y}} \right) N_{isa,y} W_{isa,y}$$$$F_{isa,y} = \bar{F}_{f_i}e^{\epsilon_{{f_i},y} s_{{f_i}sa}}$$ An index is predicted from catchability $q_{f_iy}$ and selectivity $S_{f_isa,y}$, in biomass $\hat{I:B}_{{f_i},y}$ (mt) or in numbers $\hat{I:N}_{{f_i},y}$ (thousands of fish), per Observation_units:

$$\hat{I:B}_{{f_i},y} = \sum_s\sum_a \left( \tilde{N}_{isa,y}\, W_{{f_i}sa,y}\, S_{{f_i}sa,y}\, q_{{f_i}y} \right)$$$$\hat{I:N}_{{f_i},y} = \sum_s\sum_a \left( \tilde{N}_{isa,y}\, S_{{f_i}sa,y}\, q_{{f_i}y} \right)$$

What differs between the two fleet types is the numbers the index sees, $\tilde{N}_{isa,y}$. A survey index is a snapshot: the survey is on the grounds at one point in the year, $Month_{{f_i},y}$, and sees the numbers that have survived to it. That month is read from the index_data row, not from the fleet, so a fleet whose survey timing has moved carries a different value per year.

$$\tilde{N}_{isa,y} = N_{isa,y}\, e^{-(Month_{{f_i},y}/12) Z_{isa,y}} \qquad \text{(Fleet\_type = Survey)}$$

A fishery index is CPUE, which accumulates over the whole year alongside the catch rather than at an instant. Baranov gives $\hat{C}_{isa,y} = F_{isa,y} \bar{N}_{isa,y}$, and effort cancels the $F$ (since $F_{isa,y} = q^{E}_{f_i} E_{f_i,y} s_{{f_i}sa}$), leaving the same year-average numbers the catch equation uses:

$$\tilde{N}_{isa,y} = \bar{N}_{isa,y} = N_{isa,y}\, \frac{1 - e^{-Z_{isa,y}}}{Z_{isa,y}} \qquad \text{(Fleet\_type = Fishery)}$$

$\bar{N}$ is also what a fishery’s age composition is built on, so a fleet’s index, catch and comps are consistent by construction, and $Month_{{f_i},y}$ is not read for a fishery’s index rows. The model keys this on Fleet_type; there is no column to set. The distinction moves trend as well as scale, so a constant $q$ cannot absorb it: at Month = 6, $e^{-Z/2}$ approximates $\bar{N}/N$ to within 1% at $Z = 0.5$ but is 9.6% low at $Z = 1.5$. Stock Synthesis splits the same way on its per-fleet survey timing (SS_expval.tpl: $timing \ge 0 \rightarrow e^{-Z \cdot timing}$; $timing < 0$, required of every fishing fleet $\rightarrow (1 - e^{-Z})/Z$). See vignette("model-options-and-functionality") for the seasonal-fishery case and the measured trend errors.

Selectivity parameterizations

Selectivity is controlled by Selectivity_index, Selectivity, N_sel_bins, Time_varying_sel, Time_varying_sel_sd and Bin_first_selected on the fleet_control sheet:

• Selectivity_index: Fleets given the same value share one selectivity parameter block.
• Selectivity: Selectivity parameterization for the fleet: 0 or "Fixed" = empirical selectivity provided in emp_sel; 1 or "Logistic"; 2 or "NonParametric" sensu Ianelli et al 2018; 3 or "DoubleLogistic"; 4 or "DescendingLogistic"; 5 or "Hake", non-parametric sensu Taylor et al 2014; 6 or "2DAR1" for bin by year; 7 or "3DAR1" sensu Cheng et al (2024); 8 or "DoubleNormal", a four-parameter double normal (peak, ascending and descending widths, and a right-tail floor); this is the Stock Synthesis pattern 24 shape without its initial-selectivity and plateau-width parameters; 9 or "NonParametricPM", Ianelli non-parametric with the ADMB AMAK ("pm") penalty; 11 or "LogisticPM", the AMAK logistic with a free age-1 selectivity; 13 or "NonParametricIntegrable", the Ianelli base curve whose deviations have a proper density and so integrate under random_sel (Time_varying_sel picks the structure: "Off", "IID" or "RandomWalk"). There is no 10, 12 or 14; the integer codes are frozen, because renumbering one would silently reinterpret every saved config and every fitted .rds that holds it. Whether a form is age- or length-based is set by Selectivity_dimension, not by the code.
• N_sel_bins: Number of age or length bins (per Selectivity_dimension) estimated for the non-parametric and AR1 forms; not read by the others. It is a count of bins, not an age: estimating up to age 6 means N_sel_bins = 6 when minage = 1 and 7 when minage = 0.
• Time_varying_sel: Whether time-varying selectivity is estimated, for the logistic, double-logistic, descending-logistic and non-parametric forms. 0 or "Off" = no; 1 or "IID" = penalized deviates given Time_varying_sel_sd, or random effects with that sd estimated under fit_mod(random_sel = TRUE); 3 or "Block" = time blocks with no penalty; 4 or "RandomWalk" = random walk following Dorn; 5 or "RandomWalkAscending" = random walk on the ascending portion of the double logistic only. Value 2 ("AR1") was removed in 5.16.0: it was never an AR1, and an autocorrelated deviation is now a selectivity linkage, linkage_spec(~ ar1(1 | Year)). The switch as a whole is soft-deprecated in favour of that grammar.
• Time_varying_sel_sd: The sd of the time-varying selectivity deviates. Under fit_mod(random_sel = FALSE) it is fixed at this value, per fleet. Under random_sel = TRUE the deviates are integrated out and the sd is estimated, so this is only a starting value, and fleets sharing a Selectivity_index share one estimated sd started at the geometric mean of theirs. The 2DAR1 and 3DAR1 forms read it under every Time_varying_sel, as the conditional sd of their deviation field, so they always need a positive value. Time_varying_q_sd is the catchability equivalent, under random_q.
• Bin_first_selected: First age/length bin with non-zero selectivity; the curve is 0 below it. A 1-based bin ordinal on the fleet’s Selectivity_dimension.

Selectivity penalty weights

Sel_curve_pen1/2/3 are penalty weights. Each penalty is a Gaussian sum-of-squares ($\text{weight}\cdot x^2 = x^2 / (2\sigma^2)$), so you may instead supply a standard deviation and let switch_check() convert it with $\text{weight} = 1/(2\sigma^2)$. Supplying an SD sets the matching Sel_curve_pen column; a Sel_curve_pen you set directly is never overwritten.

Which form reads which slot, and what it charges:

Slot	SD column	Reads the sign?	NonParametric (2)	NonParametricPM (9)	NonParametricIntegrable (13)	LogisticPM (11)
Sel_curve_pen1	Sel_shape_sd	only on 9 under Sel_shape_mode = "Directional"	shape, on the realized curve	shape, on the realized curve	shape, on the base curve	random walk on realized log-selectivity
Sel_curve_pen2	Sel_curvature_sd	no	2nd-difference curvature, on the realized curve	2nd-difference curvature, on the realized curve	2nd-difference curvature, on the base curve	not read
Sel_curve_pen3	Sel_devmag_sd	no	not read	deviation magnitude	not read	random walk on the age-1 deviates

Realized or base. Forms 2 and 9 charge the shape and curvature penalties on the realized curve, accumulated over sexes and years (form 2 under "IID" or "RandomWalk", form 9 from Sel_start_year under "RandomWalk"; 9 refuses "IID"). Form 13 charges them once per sex on the base coefficients. The same weight therefore bites harder on a time-varying 2 or 9 – retune rather than copy one across forms.

Three points that decide whether a value you set does anything:

• The direction. Sel_shape_dir = "Increasing" is implemented only by NonParametricPM (9) under Sel_shape_mode = "Directional", where the template switches the one-sided penalty from the decreasing side to the increasing one. Everywhere else – including 9 under "Smooth", which is two-sided – the weight is applied as written, so "Increasing" is refused rather than silently negating it.
• The sign. A weight multiplies a squared deviation, so a negative one rewards that deviation. data_check() refuses a negative Sel_curve_pen wherever the template reads it as a weight – which leaves four cases it does not refuse: the directional case above, a Fleet_type = "Off" fleet (charged no selectivity penalty at all), a fleet that follows another’s Selectivity_index and shares its selectivity form (the group’s penalty is charged once, on the lead), and a slot charged on deviates the fleet does not estimate, per the next point.
• Whether the deviates exist. The three cells above marked as random walks or deviation magnitude are charged on the time-varying deviates, so under Time_varying_sel = "Off" they are identically zero and the weight is never read. Forms 2 and 13 do estimate deviates, but score them with a Gaussian density on Time_varying_sel_sd, which is why they never read slot 3.

The 2D/3D-AR1 forms (6/7) reuse Sel_curve_pen1/2/3 for logit-scale correlations, not weights, and keep those columns. The Hake form (5) reads no Sel_curve_pen slot.

Fixed selectivity

Input: Selectivity = 0 or "Fixed", with selectivity-at-age supplied in emp_sel.

Non-parametric (Type 1)

Input: Selectivity = 2 or "NonParametric", N_sel_bins = 6, Sel_curve_pen1 = 12.5, Sel_curve_pen2 = 200, and Bin_first_selected = 1. See Selectivity penalty weights above for the Sel_curve_pen slots and the standard-deviation columns that set them.

Equation: $$sel_{f_i sa}=e^{\phi_{f_isa}}$$

Logistic

Input: Selectivity = 1 or "Logistic", N_sel_bins = NA, Time_varying_sel = 0, and Time_varying_sel_sd = NA

Equation: $$sel_{f_i sa}=1/(1+e^{-Slp_{f_i} (a-Inf_{f_i} ) } )$$

Time-varying logistic

Input: Selectivity = 1 or "Logistic", N_sel_bins = NA, and Time_varying_sel_sd = 0.1 or another starting value.

Equation: $$sel_{f_i sa}=1/(1+e^{-(Slp_{f_i}+e^{\phi_{f_i sy}^{slp} } ) (a-(Inf_{f_i}+\phi_{f_i sy}^{Inf} ) ) } )$$ Penalized likelihood: Time_varying_sel = 1 $$\phi_{f_i sy} \sim N(0,\sigma^{\phi}_{f_i})$$ Random effect: Time_varying_sel = 1 with fit_mod(random_sel = TRUE), where the deviations are integrated out by the Laplace approximation and $\hat\sigma^{\phi}_{f_i}$ is estimated rather than fixed at Time_varying_sel_sd $$\phi_{f_i sy} \sim N(0,\hat\sigma^{\phi}_{f_i})$$ Block: Time_varying_sel = 3, where the main parameters are set to 0 and $$\phi_{f_i sy} = \phi_{f_i sy_{block}}$$ Penalized likelihood random walk: Time_varying_sel = 4 $$\phi_{f_i sy}-\phi_{f_i sy-1} \sim N(0,\sigma^{\phi}_{f_i})$$

Double logistic

Input: Selectivity = 3 or "DoubleLogistic", N_sel_bins = NA, Time_varying_sel = 0, and Time_varying_sel_sd = NA

Equation: $$sel_{f_i sa}=1/(1+e^{-Slp1_{f_i} (a-Inf1_{f_i} ) } )(1-1/(1+e^{-Slp2_{f_i} (a-Inf2_{f_i} ) } ))$$

Time-varying double logistic 1 (ascending and descending time-varying parameters)

$$sel_{f_i sa}=1/(1+e^{-(Slp1_{f_i}+e^{\phi_{f_i sy}^{slp1} } ) (a-(Inf1_{f_i}+\phi_{f_i sy}^{Inf1} ) ) } )(1-1/(1+e^{-(Slp2_{f_i}+e^{\phi_{f_i sy}^{slp2} } ) (a-(Inf2_{f_i}+\phi_{f_i sy}^{Inf2} ) ) } ))$$

The four deviate structures – Time_varying_sel 1 penalized, 1 as a random effect under fit_mod(random_sel = TRUE), 3 block, and 4 random walk – are as under Time-varying logistic above.

Time-varying double logistic 2 (random walk ascending-time varying parameters)

Input: Selectivity = 3 or "DoubleLogistic", N_sel_bins = NA, Time_varying_sel = 5, and Time_varying_sel_sd = 0.05 or another starting value

Equation: $$sel_{f_i sa}=1/(1+e^{-(Slp1_{f_i}+e^{\phi_{f_i sy}^{slp1} } ) (a-(Inf1_{f_i}+\phi_{f_i sy}^{Inf1} ) ) } )(1-1/(1+e^{-Slp2_{f_i} (a-Inf2_{f_i} ) } ))$$

$$\phi_{f_i sy} - \phi_{f_i sy-1} \sim N(0,\sigma^{\phi}_{f_i})$$

Descending logistic

Input: Selectivity = 4 or "DescendingLogistic", N_sel_bins = NA, Time_varying_sel = 0, and Time_varying_sel_sd = NA

Equation: $$sel_{f_i sa}=1-1/(1+e^{-Slp2_{f_i} (a-Inf2_{f_i} ) } )$$

Time-varying descending logistic

Input: Selectivity = 4 or "DescendingLogistic", N_sel_bins = NA, Time_varying_sel = 1, and Time_varying_sel_sd = 0.05 or another starting value

Equation: $$sel_{f_i sa}=1-1/(1+e^{-(Slp2_{f_i}+e^{\phi_{f_i sy}^{slp2} } ) (a-(Inf2_{f_i}+\phi_{f_i sy}^{Inf2} ) ) } )$$

The four deviate structures are as under Time-varying logistic above.

Non-parametric (Type 2)

For each age $a \geq A_{min}$ there is an incremental selectivity parameter $p_{f_i say}$ for the fleet $f_i$. Selectivity at age is calculated as follows:

$$sel_{f_i sa}=e^{sel^{`}_{f_i sa} - max(sel^{`}_{f_i sa})}$$ where $$ sel^{`}_{f_i sa} = \sum_{i=A_{min}}^a p_{f_i sa}$$ Selectivity is 0 for $a < A_{min}$, and $p_{f_i sa} = 0$ above minage + N_sel_bins, so selectivity is constant beyond the last estimated bin. minage, set in the data, is the youngest age modeled.

Input: Selectivity = 5 or "Hake", N_sel_bins = 6, Time_varying_sel = NA, Time_varying_sel_sd = NA, and Bin_first_selected =$A_{min}$

Time-varying non-parametric (Type 2 for hake)

As above, but deviations on $p_a$ carry annual changes in selectivity, given a fixed or estimated standard deviation ($\sigma^{\phi}_{f_i}$): $$p_{f_i say}=p_{f_i sa}+\phi_{f_i say}$$

Penalized likelihood: Time_varying_sel = 1 $$\phi_{f_i sy} \sim N(0,\sigma^{\phi}_{f_i})$$ Random effect: Time_varying_sel = 1 with fit_mod(random_sel = TRUE), where the deviations are integrated out by the Laplace approximation and $\hat\sigma^{\phi}_{f_i}$ is estimated rather than fixed at Time_varying_sel_sd $$\phi_{f_i sy} \sim N(0,\hat\sigma^{\phi}_{f_i})$$

Catchability parameterizations

Catchability is controlled by the following fleet_control columns:

• Catchability_index Fleets given the same value share one catchability parameter.
• Catchability Text or integer switch specifying catchability. 0 or “Fixed” = fixed at Catchability_init; 1 or “Estimated” = a single free parameter; 2 or “Estimated-with-prior” = a single free parameter with a prior; 3 or “Analytical” = analytical q from Walters and Ludwig 1994; 4 or “PowerEquation” = power equation, declared but not yet implemented; 5 or “Environmental” = linear equation log(q_y) = q_mu + beta * index_y; 6 or “AR1” = REMOVED in 5.12.0 (it never estimated the deviates, so q was constant) – express the Rogers et al. form as a q linkage with ar1(1 | Year) and observe =; 7 or “AnalyticalArith” = arithmetic-mean analytical q, the AMAK/ebswp form, paired with the MVN survey likelihood.
• Catchability_init Starting or fixed value for catchability
• Catchability_prior_sd Log-scale SD of the q prior: dnorm(log_q, log(Catchability_init) - bias_adjust_proc * sd^2/2, sd)
• Time_varying_q Whether a time-varying q is estimated. 0 = no; 1 = penalized deviate; 2 = AR1 (REMOVED in 5.16.0; data_check() errors); 3 = time blocks with no penalty; 4 = random walk from the mean following Dorn 2018 (dnorm(q_y - q_y-1, 0, sigma)). Under Catchability = 5 it instead holds the env_data index for log(q_y) = q_mu + beta * index_y.
• Time_varying_q_sd Fixed or starting sd of the time-varying q deviations, read under Time_varying_q 1 and 4 and estimated under fit_mod(random_q = TRUE)

Linear time-invariant catchability formulations

All take Time_varying_q = 0.

Catchability = 0 fixes $q_{f_i}$ at Catchability_init.

Catchability = 1 A free parameter $\hat{q}_{f_i}$ per index fleet $f_i$

$$q_{f_i} = \hat{q}_{f_i}$$

Catchability = 2 As above, with a lognormal prior on $q$

$$log(\hat{q}_{f_i}) \sim N\left(log(Q\ prior) - b \cdot Q\ sd\ prior^2/2,\ Q\ sd\ prior\right), \quad b = \texttt{bias\_adjust\_proc}$$

Q prior (Catchability_init) is the prior mean of $q$ under fit_control(bias_adjust_proc = TRUE), the default, and its median otherwise; Q sd prior is on the log scale.

Catchability = 3 $\hat{q}_{f_i}$ is derived analytically following Walters and Ludwig (1994): Walters, C. and Ludwig, D. (1994) ‘Calculation of Bayes posterior probability distributions for key population parameters’, Canadian Journal of Fisheries and Aquatic Sciences, 51(3), pp. 713-722, doi:10.1139/f94-071. * CPUE in weight, time-invariant survey SD: $$\hat{q}_{f_i} = \exp\left( \sum_y log \left[ \frac{CPUE:B_{{f_i},y}} {\sum_s\sum_a \left( N_{isa,y} e^{-(Month_{{f_i},y}/12) Z_{isa,y}} W_{{f_i}sa,y} S_{{f_i}isa,y} \right)} \right]/n \right)$$ * CPUE in numbers, time-invariant survey SD: $$\hat{q}_{f_i} = \exp\left( \sum_y log \left[ \frac{CPUE:I_{{f_i},y}} {\sum_s\sum_a \left( N_{isa,y} e^{-(Month_{{f_i},y}/12) Z_{isa,y}} S_{{f_i}isa,y} \right)} \right]/n \right)$$ * CPUE in weight, time-varying survey SD: $$\hat{q}_{f_i} = \exp\left( \sum_y \left( log \left[ \frac{CPUE:B_{{f_i},y}} {\sum_s\sum_a N_{isa,y} e^{-(Month_{{f_i},y}/12) Z_{isa,y}} W_{{f_i}sa,y} S_{{f_i}isa,y}} \right]/\sigma^{2}_{f_iy} \right)/ \sum_y{\frac{1}{\sigma^{2}_{f_iy}}} \right)$$ * CPUE in numbers, time-varying survey SD: $$\hat{q}_{f_i} = \exp\left( \sum_y \left( log \left[ \frac{CPUE:I_{{f_i},y}} {\sum_s\sum_a N_{isa,y} e^{-(Month_{{f_i},y}/12) Z_{isa,y}} S_{{f_i}isa,y}} \right]/\sigma^{2}_{f_iy} \right)/ \sum_y{\frac{1}{\sigma^{2}_{f_iy}}} \right)$$

Time-varying catchability formulations

Under Catchability = 2 the prior above still applies to mean catchability.

Annual catchability estimated via penalized likelihood or as random effects

Catchability %in% c(1:2) and Time_varying_q %in% c(1, 4) $$\hat{q}_{f_iy} = e^{\bar{q}_{f_i}+ \omega_{f_iy}}$$$\omega_{f_iy} \sim N(0,$Time_varying_q_sd$)$, with Time_varying_q_sd either fixed or estimated under fit_mod(random_q = TRUE)

Catchability time-blocks

Catchability %in% c(1:2), Time_varying_q = 3, with the blocks taken from index_data. $$\hat{q}_{f_iy} = e^{\omega_{f_i block_y}}$$

Annual catchability via random walk

Catchability %in% c(1:2), Time_varying_q = 4 $$\hat{q}_{f_iy} = e^{\bar{q}_{f_i}+ \omega_{f_iy}}$$$\omega_{f_iy} \sim N(\omega_{f_iy-1},$Time_varying_q_sd$)$