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}\]
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 represents a predator’s prey and size preference, together with factors such as spatial overlap between predator and prey distributions.
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)\]
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\]
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)\]
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 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 rationPvalue Tuning parameter \(\varphi_p\) scaling the maximum consumption
used for rationfday 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 gramsCB \(\beta^{\delta}_p\), the slope of that
functionQc, 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
belowration_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.
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)\]
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}
\]
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 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.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:
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.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.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.
Input: Selectivity = 0 or "Fixed", with
selectivity-at-age supplied in emp_sel.
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}}\]
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} ) } )\]
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}) \]
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} ) } ))\]
\[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.
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}) \]
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} ) } )\]
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.
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}\)
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 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
catchabilityCatchability_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)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)\]
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\()\)