Conditioning Model
The Southern Bluefin Tuna (SBT) operating model is used for management strategy evaluation (Punt et al. 2016) and has been developed and refined through the Commission for the Conservation of Southern Bluefin Tuna (CCSBT) process (Anon. 2024). This article documents the current RTMB implementation and keeps links between mathematical notation and the R code used by the package.
Model setup using RTMB
Implemented in: R/sbt-fit.R (sbt_fit, sbt_obj), R/fit-lifecycle.R (sbt_build_object, sbt_optimise), and R/model.R (sbt_model).
The supported workflow keeps processed data, configuration, optimisation, and posterior state in one portable sbt_fit. Missing default configuration is resolved when the objective is first built or optimised:
fit <- sbt_fit(data_in = data_in)
fit <- sbt_optimise(
fit,
control = list(eval.max = 10000L, iter.max = 10000L)
)
obj <- sbt_obj(fit)For specialised model-development work, sbt_build_object() exposes the transient RTMB objective before optimisation. Parameters can be fixed by supplying a map through sbt_add_map(); entries mapped to factor(NA) are fixed.
Model Structure
Implemented in: R/get-data.R (get_data) and consumed by R/model.R (sbt_model).
The model is a single age-structured stock with ages 0 to 30 and two seasons per year. Six fisheries contribute catch dynamics, and a seventh selectivity block is used for CPUE.
| Fishery | Description | Model pulse/season | Catch-year accounting |
|---|---|---|---|
| LL1 | Japanese LL areas 4-9 and other longline catches | 2 | Jan-Dec |
| LL2 | Taiwanese albacore LL and gillnet | 2 | Jan-Dec |
| LL3 | Japanese LL in Area 2 | 1 | Jul-Jun |
| LL4 | Japanese spawning (Area 1) | 1 | Jul-Jun |
| Indonesian | Indonesian spawning | 1 | Jul-Jun |
| Australian Surface | Surface fishery | 1 | Jul-Jun |
| CPUE block | Selectivity-only block used in the CPUE likelihood | 2 | Index year; not catch |
The two model seasons locate fishing pulses within the annual population dynamics. They do not redefine the fishery-specific catch-year convention: LL1 and LL2 catches are accounted January-December, while LL3, LL4, Indonesian, and Australian surface catches are accounted July-June.
Population dynamics
Implemented in: R/dynamics.R (do_dynamics, get_harvest_rate) and called from R/model.R (sbt_model).
Let N_{y,s,a} be numbers-at-age at year y, season s \in \{1,2\}, age a \in \{0,\dots,A\}, and let S_a = \exp(-0.5M_a).
Season 1 to season 2:
N_{y,2,a} = N_{y,1,a}(1-h_{y,1,a})S_a
Season 2 to next year:
N_{y+1,1,a+1} = N_{y,2,a}(1-h_{y,2,a})S_a, \quad a=0,\dots,A-1
Plus-group update:
N_{y+1,1,A} = N_{y,2,A-1}(1-h_{y,2,A-1})S_{A-1} + N_{y,2,A}(1-h_{y,2,A})S_A
Recruitment:
N_{y+1,1,0} = R_{y+1}
where h_{y,s,a} = \sum_{f=1}^{6} H_{y,s,f,a}.
For fisheries with removal_switch_f = 0, harvest is:
F_{y,s,f} = \frac{C_{y,s,f}}{\sum_a N_{y,s,a}s_{f,y,a}w_{f,y,a}}, \quad H_{y,s,f,a} = F_{y,s,f}s_{f,y,a}
For fisheries with removal_switch_f = 1, age-specific removals are:
D_{y,s,f,a} = \frac{C_{y,s,f}A_{y,s,f,a}} {\sum_j A_{y,s,f,j}w_{f,y,j}}, \qquad H_{y,s,f,a} = \frac{D_{y,s,f,a}}{N_{y,s,a}}
with A_{y,s,f,a} from sliced age composition (af_sliced_ysfa).
Catch predictions:
\hat{C}^{\text{numbers}}_{f,y,a} = \sum_{s=1}^{2} D_{y,s,f,a}
\hat{C}^{\text{weight}}_{y,s,f} = \sum_a D_{y,s,f,a}w_{f,y,a}
Raw combined survival is
S^{\rm raw}_{y,s,a}=1-\sum_f H_{y,s,f,a}.
Accepted fits and posterior draws require S^{\rm raw}_{y,s,a}\geq0.1 for every populated year, season, and age. Trial parameter vectors below this limit retain their exact conditioned catches, receive the fixed penalty 10^6\,0.01(S^{\rm raw}_{y,s,a}-0.1)^2, and use a positive posfun continuation only for population propagation. Raw and continued harvest are reported separately; a non-zero continuation penalty is never accepted as a production fit or draw.
An assessment can additionally configure a smooth preventive wall on the raw combined harvest rate. With strength A, onset H_0, biological ceiling H_c, and scale s, define
g(H)=s\log\left(1+\exp\left(\frac{H-H_0}{s}\right)\right), \qquad p_{y,s,a}^{\rm wall}=A\left\{\frac{g(h_{y,s,a})}{g(H_c)}\right\}^2.
The normalization gives p^{\rm wall}=A for one cell at H=H_c. The model adds lp_{\rm wall}=\sum_{y,s,a}p_{y,s,a}^{\rm wall} to the objective and reports it separately from the feasibility continuation. Setting A=0 disables the wall. The wall never rescales conditioned catches or changes the population update; the raw H\leq0.9, positive-abundance, and zero-continuation acceptance gates remain in force.
Spawning biomass is:
\text{SSB}_y = \sum_a N_{y,1,a} \phi_{y,a}
Stock-Recruitment
Implemented in: R/recruitment.R (get_recruitment, get_recruitment_prior, get_rho) and R/dynamics.R (get_initial_numbers).
Recruitment follows Beverton-Holt with a depensation term and lognormal deviation:
R_y = \frac{\alpha \,\text{SSB}_y}{\beta + \text{SSB}_y} \left(1 - \exp\left(\frac{\log(0.5)\,\text{SSB}_y}{\text{sr\_dep}\,B_0}\right)\right) \exp(\varepsilon_y - 0.5\sigma_r^2)
In R/recruitment.R, sr_dep defaults to 1e-10.
Initial equilibrium quantities from get_initial_numbers():
\text{rel}N_0=1,\quad \text{rel}N_a=\text{rel}N_{a-1}\exp(-M_{a-1}),\quad \text{rel}N_A=\frac{\text{rel}N_A}{1-\exp(-M_A)}
R_0 = \frac{B_0}{\sum_a \phi_{1,a}\text{rel}N_a},\quad \alpha=\frac{4hR_0}{5h-1},\quad \beta=\frac{B_0(1-h)}{5h-1},\quad N^{\text{init}}_a = R_0\text{rel}N_a
Recruitment-deviation prior in get_recruitment_prior():
\varepsilon_{1:(Y-3)} \sim N(0,\sigma_r^2), \quad \varepsilon_{(Y-2):Y} \sim \text{AR1}(\phi, \sigma_r)
where \phi is estimated by get_rho() as the empirical correlation of \epsilon_y and \epsilon_{y+1} over years 1965 to the terminal year - 5.
The unfished spawning biomass B_0 is parameterized through the initial equilibrium age structure. For a reference recruitment R_0, the initial numbers are:
N^{\text{init}}_a = R_0\text{rel}N_a
and the Beverton-Holt parameters can be expressed in terms of steepness h:
\alpha = \frac{4hR_0}{5h - 1}, \quad \beta = \frac{B_0(1-h)}{5h - 1}
The historical ADMB operating model allowed a regime shift in carrying capacity from 1978 onward. In the current RTMB package this is controlled through the parameter list and get_map() rather than through ADMB control-file flags.
Selectivity
Implemented in: R/selectivity.R (get_selectivity, get_selectivity_prior) and called from R/model.R (sbt_model).
Selectivity is estimated for seven blocks (six fisheries plus CPUE). For fishery f, let a^{\min}_f and a^{\max}_f be the estimated age range, and let c_f be years where sel_change_year_fy[f,y] = 1.
At each change year y \in c_f, unnormalized selectivity is:
s'_{f,y,a} = \exp(\lambda_{f,r,a}), \quad a \in [a^{\min}_f, a^{\max}_f]
and is normalized by the mean over the estimated age range:
s_{f,y,a} = \frac{s'_{f,y,a}}{\text{mean}_{a^{\min}_f:a^{\max}_f}(s'_{f,y,a})}
For ages above a^{\max}_f, if sel_end_f[f] = 1 then s_{f,y,a}=s_{f,y,a^{\max}_f}, otherwise these ages remain zero. For years without a change point, selectivity is copied forward from the previous year.
The current implementation normalizes each selectivity vector by the mean over the estimated age range. This corresponds to the preferred ADMB jim_select = 1 parameterization; older notes that normalize by a reference age are retained only as historical context.
The prior for each fishery’s log-selectivity matrix \Lambda_f is separable AR1 across change-year and age dimensions:
\log p(\Lambda_f) = \log \left[\text{dseparable}\left(\text{AR1}(\rho^y_f), \text{AR1}(\rho^a_f)\right)\right]
with scale
\sigma_f = \frac{\exp(\eta_f)}{\sqrt{1-(\rho^y_f)^2}\sqrt{1-(\rho^a_f)^2}}
where \eta_f = \text{par\_log\_sel\_sigma}[f].
A more complete explanation perhaps is detailed as follows:
Selectivity Prior (Separable AR1)
For fishery f, let the log-selectivity parameters be
\Lambda_f = \{ \lambda_{r,a} \},
where r indexes selectivity change-years and a indexes age.
Vectorizing the matrix,
\mathrm{vec}(\Lambda_f) \sim \mathrm{MVN}(0, \Sigma_f).
We assume a separable covariance structure:
\Sigma_f = \sigma_f^2 \left( \Sigma^{(y)}_f \otimes \Sigma^{(a)}_f \right),
where \otimes denotes the Kronecker product.
AR1 Structure Across Change-Years
For change-years,
\text{Cov}(\lambda_{r,a}, \lambda_{r',a}) = \sigma_f^2 \, (\rho^y_f)^{|r-r'|}.
AR1 Structure Across Ages
For ages,
\text{Cov}(\lambda_{r,a}, \lambda_{r,a'}) = \sigma_f^2 \, (\rho^a_f)^{|a-a'|}.
Combined Separable Covariance
Under separability,
\text{Cov}(\lambda_{r,a}, \lambda_{r',a'}) = \sigma_f^2 \, (\rho^y_f)^{|r-r'|} (\rho^a_f)^{|a-a'|}.
Scale Parameter
The marginal standard deviation is defined as
\sigma_f = \frac{\exp(\eta_f)} {\sqrt{1-(\rho^y_f)^2} \sqrt{1-(\rho^a_f)^2}},
where
\eta_f = par \log sel \sigma[f].
Growth
Implemented in: R/length-weight.R (get_length_at_age, get_dl) and giR/dynamics.R (get_phi).
Mean length-at-age is time varying by year and season (\text{length\_mu}_{y,s,a}), and SD-at-age is fixed (\text{length\_sd}_a). Age-length keys (alk_ysal) are computed from these inputs and used in composition likelihoods. The dl_l distribution weights used for the spawning-output integral are common across years and ages because the integration grid is standardized from -1.98 to 1.98 length standard deviations.
Spawning output-at-age in year y (get_phi) is:
\phi_{y,a} = \frac{\sum_l d_l\,l^{3\psi}\,\text{mat}(l)}{\phi^*_{y,A}}
with maturity-at-length
\text{mat}(l)=\frac{1}{1 + 19^{(L_{50}-l)/(L_{95}-L_{50})}}
and normalization by the oldest age value \phi^*_{y,A}.
Natural Mortality
Implemented in: R/natural-mortality.R (get_M) and called from R/model.R (sbt_model).
Natural mortality-at-age is generated by get_M(min_age, max_age, age_increase_M, m0, m4, m10, m30).
Define:
p = \frac{\log\left(\frac{m_0 - m_4}{m_0 - m_{10}}\right)}{\log(1/3)}, \quad \Delta_M = \frac{m_{30} - m_{10}}{A - a_{\text{inc}}}
where A is the maximum age and a_{\text{inc}} is age_increase_M.
Then:
M_a = \begin{cases} m_0 & a \in \{0,1\} \\ m_0 - (m_0 - m_{10})\left(\frac{a-1}{9}\right)^p & 2 \le a \le 9 \\ m_{10} & 10 \le a \le a_{\text{inc}} \\ M_{a-1} + \Delta_M & a_{\text{inc}} < a < A \\ m_{30} & a = A \end{cases}
In the current RTMB implementation, all four mortality parameters are estimable unless mapped out in get_map().
Tagging Model
Implemented in: R/likelihoods.R (get_tag_like) and called from R/model.R (sbt_model).
Tag recapture data are modeled by cohort index k, tagger group t, release age i, and recapture age j (get_tag_like).
Adjusted releases (same-year recaptures removed and expanded for reporting):
N^{\text{adj}}_{k,t,i} = N_{k,t,i} - \operatorname{round}\left(\frac{R_{k,t,i,i}}{\nu_{k+i-2,i}}\right)
where \nu_{y,a} is the age- and year-specific reporting rate matrix (tag_rep_rates_ya).
To account for incomplete mixing in season 1, the tagged-fish harvest component is scaled by:
h^+_{k,j} = H_{\text{tag}}\,h_{y,1,j+1}, \quad y = \text{minK} + k - 1 + j
with H_{\text{tag}} = \exp(\text{par\_log\_tag\_H\_factor}).
Define survival and recapture terms:
S^{(1)}_{k,t,j} = (1-h^+_{k,j})(1-h_{y,2,j+1})\exp(-M_{j+1}-\Omega_t) S^{(2)}_{k,t,j} = (1-h^+_{k,j})(1-h_{y,2,j+1})\exp(-M_{j+1}-2\Omega_t) f^{(1)}_{k,t,j} = h^+_{k,j} + (1-h^+_{k,j})\exp\{-0.5(M_{j+1}+\Omega_t)\}h_{y,2,j+1} f^{(2)}_{k,t,j} = h^+_{k,j} + (1-h^+_{k,j})\exp\{-0.5(M_{j+1}+2\Omega_t)\}h_{y,2,j+1}
For the year immediately after release:
S^{\ast(1)}_{k,t,i} = \exp(-M_{i+1}-\Omega_t), \quad S^{\ast(2)}_{k,t,i} = \exp(-M_{i+1}-2\Omega_t)
Recapture probability (at least one tag returned):
p_{k,t,i,j} = \begin{cases} \left(2\xi_t S^{\ast(1)}_{k,t,i}f^{(1)}_{k,t,j} - \xi_t^2 S^{\ast(2)}_{k,t,i}f^{(2)}_{k,t,j}\right)\nu_{k+j-2,j}, & j=i+1 \\ \left(2\xi_t S^{\ast(1)}_{k,t,i}\prod_{s=i+1}^{j-1}S^{(1)}_{k,t,s}f^{(1)}_{k,t,j} - \xi_t^2 S^{\ast(2)}_{k,t,i}\prod_{s=i+1}^{j-1}S^{(2)}_{k,t,s}f^{(2)}_{k,t,j}\right)\nu_{k+j-2,j}, & j>i+1 \end{cases}
where \xi_t and \Omega_t correspond to immediate and continuous tag shedding (tag_shed_immediate, tag_shed_continuous).
Predicted recaptures are:
\hat{R}_{k,t,i,j} = N^{\text{adj}}_{k,t,i} p_{k,t,i,j}
This is a Brownie-style tag return model. Tagging is assumed to occur at the start of the calendar year, newly tagged fish are allowed to be incompletely mixed in the first season after release, reporting rates vary by year and age, and tag retention depends on tagger group through immediate retention and continuous shedding parameters. The incomplete-mixing sensitivity from the ADMB implementation is represented in RTMB by the multiplicative harvest-rate factor H_{\text{tag}}.
Catch underreporting
Implemented in data preprocessing: R/get-data.R (get_data).
In the current RTMB workflow, catch adjustments are handled in get_data() using scenario multipliers (catch_LL1_case, catch_surf_case) plus additional historical undercatch components (catch_UA). The retired catch_UR_on interface is rejected: there is no standalone I42-style underreporting switch in R/model.R.
Predicted quantities
Implemented across R/dynamics.R and R/likelihoods.R; section-specific mappings are listed below.
Catch-at-age and catch-at-length
Implemented in: R/dynamics.R (do_dynamics) for catch predictions, and R/likelihoods.R (get_age_like, get_length_like) for observation models.
Predicted catch-at-age used by age and length likelihoods is accumulated across seasons:
\hat{C}^{\text{num}}_{f,y,a} = \sum_{s=1}^{2} H_{y,s,f,a}N_{y,s,a}
For length data, predicted proportions are generated using ALKs:
\hat{L}_{i,l} \propto \sum_a \hat{C}^{\text{num}}_{f_i,y_i,a}\,\text{ALK}_{y_i,s_i,a,l}
with lower bins optionally pooled via lf_minbin.
CPUE
Implemented in: R/likelihoods.R (get_cpue_like).
The CPUE likelihood uses fishery-7 selectivity and season-2 abundance (get_cpue_like). Define a creep adjustment:
a_1 = a_{\text{init}}, \quad a_i = a_{i-1} + \delta_{\text{cpue}} \; (i>1)
and effective abundance index:
\tilde{N}_i = \frac{\sum_{a=4}^{A} s_{7,y_i,a}N_{y_i,2,a}}{\text{mean}(s_{7,y_i,a_1:a_2})}
Unscaled log prediction:
\eta_i = \log(a_i) + \omega_{\text{cpue}}\log(\tilde{N}_i)
The series is centered and shifted by q:
\log \hat{I}_i = \eta_i - \log\left(\frac{1}{n}\sum_{i=1}^n e^{\eta_i}\right) + \log q_{\text{cpue}}
where \omega_{\text{cpue}} = \exp(\text{par\_log\_cpue\_omega}).
Tag returns
Implemented in: R/likelihoods.R (get_tag_like).
Predicted tag returns are:
\hat{R}_{k,t,i,j} = N^{\text{adj}}_{k,t,i} p_{k,t,i,j}
using the recapture probabilities defined in the Tagging section.
Aerial survey
Implemented in: R/likelihoods.R (get_aerial_survey_like).
Legacy data note: the aerial survey is retained because it is still implemented in the codebase, but it is treated as a legacy data type.
The unscaled aerial index for observation year y_i is:
\tilde{I}_i = \sum_{a=2}^{4} s^{\text{aerial}}_a\,w_{6,y_i,a}\,N_{y_i,1,a}
Selectivity vector s^{\text{aerial}} depends on aerial_switch:
-
0: (0.5,1,1), likelihood off -
1: (\exp(\theta_1),1,\exp(\theta_2)) -
2: (1,1,1) -
3: (0.33,1,0.33) -
4: (0.5,1,1)
A scaling parameter is estimated analytically each evaluation and applied to \tilde{I}_i.
Trolling survey
Implemented in: R/likelihoods.R (get_troll_like).
The trolling index is based on age-1 abundance in season 1:
\tilde{T}_i = N_{y_i,1,1}
A weighted least-squares scale factor is estimated:
q_{\text{troll}} = \frac{\sum_i \tilde{T}_i\,O_i/\sigma_i^2}{\sum_i \tilde{T}_i^2/\sigma_i^2}, \quad \sigma_i^2 = \text{SD}_{i}^2 + \tau_{\text{troll}}^2
and predictions are \hat{T}_i = q_{\text{troll}}\tilde{T}_i.
Objective Function
Implemented in: R/model.R (sbt_model), summing likelihood and prior components from helper functions.
The implemented objective in sbt_model() is:
\begin{aligned} \text{nll} = {}& lp_{\text{prior}} + \sum_f lp_{\text{sel},f} + lp_{\text{rec}} + lp_{\text{penalty}} + lp_{\text{wall}} \\ &+ \sum_i lp_{\text{af},i} + \sum_i lp_{\text{lf},i} + \sum_i lp_{\text{cpue\_lf},i} \\ &+ \sum_i lp_{\text{cpue},i} + \sum_i lp_{\text{aerial},i} + lp_{\text{aerial\_tau}} \\ &+ \sum_i lp_{\text{troll},i} + \sum_i lp_{\text{tag},i} + \sum_i lp_{\text{pop},i} \\ &+ \sum_i lp_{\text{hsp},i} + \sum_i lp_{\text{gt},i} \end{aligned}
CPUE Data
Implemented in: R/likelihoods.R (get_cpue_like).
Observed CPUE is modeled on log scale with year-specific SD:
\log O_i \sim N\left(\log \hat{I}_i,\; \sigma_i\right), \quad \sigma_i = \sqrt{\text{cpue}_\sigma_i^2 + \sigma_{\text{cpue}}^2}
where \sigma_{\text{cpue}} = \exp(\text{par\_log\_cpue\_sigma}).
Catch-at-age and Catch-at-length
Implemented in: R/likelihoods.R (get_age_like, get_length_like, get_cpue_length_like).
For fisheries with removal_switch_f = 0, both age- and length-composition likelihoods support:
-
switch = 1: multinomial (dmultinom) -
switch = 2: Dirichlet (ddirichlet) -
switch = 3: Dirichlet-multinomial (ddirmult) -
switch = 9: multinomial/KL approximation
Age compositions additionally support af_switch = 4, the uncorrelated logistic-normal LN1 likelihood recommended by Francis (2014), building on the compositional formulation of Schnute and Haigh (2007). Observed and predicted proportions are normalized within observation, with end bins pooled as defined in get_age_like() and get_length_like(). For LN1, af_logistic_normal_min_prop (default 10^{-6}) identifies the first and last observed age bins retained in each year. All younger bins are pooled into the first endpoint and all older bins into the last endpoint, for both observations and predictions. Every bin between those endpoints is retained, including small positive bins below the threshold. Exact internal observed zeros are replaced by af_logistic_normal_epsilon (default 10^{-6}), after which the composition is renormalized. Thus definitive tail zeros are compressed while happenstance internal zeros receive a small-value replacement, as recommended by Francis (2014). Compression is defined only by the observations, so its grouping does not change with model predictions. A row whose retained span contains fewer than two bins has no likelihood contribution.
For age-composition row y in fishery f, let O_{yb} and E_{yb} be the compressed observed and predicted proportions across B_y retained bins. Using the final retained bin as the additive-log-ratio reference, define
w_{yb} = \log\left(\frac{O_{yb}}{O_{yB_y}}\right) - \log\left(\frac{E_{yb}}{E_{yB_y}}\right), \qquad b=1,\ldots,B_y-1.
Under LN1 the latent log-abundance errors are independent with common variance. Therefore, with d_y=B_y-1,
q_y = \mathbf{w}_y^\mathsf{T} \left(\mathbf{I}-\frac{\mathbf{J}}{B_y}\right) \mathbf{w}_y = \sum_b w_{yb}^2-\frac{\left(\sum_b w_{yb}\right)^2}{B_y}.
The Francis sample-size adjustment and the conditionally profiled variance are
W_y^2=\frac{\bar N_f}{N_y}, \qquad \hat{\sigma}_f^2 = \frac{\sum_{y \in f}q_y/W_y^2}{\sum_{y \in f}d_y},
where \bar N_f is the mean sample size over positive-sample rows for fishery f. One dispersion is profiled separately for each age-composition fishery. The additive-log-ratio vector is evaluated with RTMB::dmvnorm() using covariance \mathbf I+\mathbf J and scale \hat{\sigma}_fW_y. The compositional Jacobian is then added explicitly. The resulting row-level negative log likelihood is
\begin{aligned} lp_{\mathrm{af},y} = {}& \frac{d_y}{2}\log(2\pi) +\sum_{b=1}^{B_y}\log O_{yb} +d_y\log\hat{\sigma}_f +d_y\log W_y\\ &+\frac{1}{2}\log B_y +\frac{q_y}{2\hat{\sigma}_f^2W_y^2}. \end{aligned}
The fitted \hat{\sigma}_f values are available in the RTMB report as af_logistic_normal_sigma_f. The Dirichlet concentration parameter par_log_af_alpha is not used by this option. The reported af_pred remains the standard model-age composition after observation-specific end pooling and before LN1 tail compression or zero replacement, so existing plots and posterior summaries retain their original age support.
Tag returns likelihood
Implemented in: R/likelihoods.R (get_tag_like).
For each non-empty release set (k,t,i) with recapture ages j=i+1,\dots,J_k, the model uses a Dirichlet-multinomial form.
Define:
\omega_{k,t,i} = \max\left(\frac{N_{k,t,i} - \phi_{\text{tag}}}{\phi_{\text{tag}} - 1},\;10^{-3}\right)
with \phi_{\text{tag}} = \text{tag\_var\_factor}, and
R_{k,t,i,\neg} = N^{\text{adj}}_{k,t,i} - \sum_{j=i+1}^{J_k} R_{k,t,i,j}, \quad p_{k,t,i,\neg} = 1 - \sum_{j=i+1}^{J_k} p_{k,t,i,j}
Then (up to constants):
\begin{aligned} \log L_{k,t,i} = {}& \log\Gamma(\omega_{k,t,i}) - \log\Gamma(N^{\text{adj}}_{k,t,i}+\omega_{k,t,i}) \\ &+ \sum_{j=i+1}^{J_k}\left[\log\Gamma(R_{k,t,i,j}+\omega_{k,t,i}p_{k,t,i,j}) - \log\Gamma(\omega_{k,t,i}p_{k,t,i,j})\right] \\ &+ \log\Gamma(R_{k,t,i,\neg}+\omega_{k,t,i}p_{k,t,i,\neg}) - \log\Gamma(\omega_{k,t,i}p_{k,t,i,\neg}) \end{aligned}
and lp_{\text{tag}} = -\log L_{k,t,i}.
Aerial Survey Likelihood
Implemented in: R/likelihoods.R (get_aerial_survey_like).
Legacy data note: this likelihood remains documented for code fidelity, but the aerial stream is considered legacy in current use.
Let \Sigma = \text{aerial\_cov} + \tau_{\text{aerial}}^2 I, with \tau_{\text{aerial}} = \exp(\text{par\_log\_aerial\_tau}).
Using unscaled predictions \tilde{I}, define:
\log q_{\text{aerial}} = \frac{\mathbf{1}^T\Sigma^{-1}(\log O - \log \tilde{I})}{\mathbf{1}^T\Sigma^{-1}\mathbf{1}}
\log \hat{I} = \log \tilde{I} + \log q_{\text{aerial}}
and
\log O \sim \text{MVN}(\log \hat{I}, \Sigma)
For diagnostics, the implementation also reports the log-determinant component:
lp_{\text{aerial\_tau}} = 0.5\log|\Sigma|
This component is already included in the complete MVN likelihood and is not added to the model objective a second time. When aerial_switch = 0, both the aerial likelihood and this diagnostic contribution are zero.
Trolling Survey Likelihood
Implemented in: R/likelihoods.R (get_troll_like).
With \tau_{\text{troll}} = \exp(\text{par\_log\_troll\_tau}) and
\sigma_i = \sqrt{\sigma_{i,\text{troll}}^2 + \tau_{\text{troll}}^2}
the model fits
\log O_i \sim N(\log \hat{T}_i, \sigma_i), \quad \hat{T}_i = q_{\text{troll}}N_{y_i,1,1}
where q_{\text{troll}} is profiled on the log scale at each evaluation:
\log q_{\text{troll}} = \frac{\sum_i \left(\log O_i - \log N_{y_i,1,1}\right) / \sigma_i^2} {\sum_i 1 / \sigma_i^2}.
Close-kin Data
Implemented in data preparation and likelihood helpers: R/get-data.R (get_data), R/likelihoods.R (get_POP_like, get_HSP_like, get_GT_like).
The close-kin data are now an agreed dataset for use within the CCSBT Operating Model (OM). While the full dataset contains many details, for modeling purposes the simplified structure is as follows:
-
Juveniles: There are individuals i \in \Theta with a year of capture y_i, length-at-capture l_i, and either:
- A distribution of ages conditional on length, p(a_i \mid l_i) (via a growth relationship), or
- Direct age a_i (in the case of a known POP).
The cohort year of birth is calculated as c_i = y_i - a_i.
-
Adults: There are individuals j \in \Xi with a year of capture y_j, length-at-capture l_j, and either:
- A distribution p(a_j \mid l_j), or
- Direct age a_j (as in a POP match).
Defining the Parental Probability
Implemented in: R/likelihoods.R (get_POP_like).
The genetic data allow us to determine whether juvenile i and adult j form a Parent-Offspring Pair (POP) with high certainty. In a simple cartoon version, if there are N adults, the probability that any one is the parent of a juvenile is 2/N (either a mother or a father).
However, in a more realistic model with multiple cohorts and age structure, the following two factors must be considered:
- Relative reproductive output varies with age (e.g., maturity, batch fecundity, spawning history).
- Relative abundance decreases with age due to mortality and recruitment dynamics.
To compute the probability that adult j is the parent of juvenile i, we balance these effects. Let:
\quad c_i = y_i - a_i: birth year of juvenile i
\quad a_j - a_i: age of adult j in the birth year of juvenile i
\quad \varphi_{a_j - a_i}: reproductive output of the adult at that age
Then the parental probability \pi_{ij} is:
\pi_{ij} = \frac{2 \varphi_{a_j - a_i}}{ \sum_a N_{y_i - a, a} \varphi_a }
where N_{y,a} is the abundance of fish of age a in year y, and \sum_a N_{c_i,a} \varphi_a is the total reproductive output in the birth year of juvenile i.
Close-kin Likelihood
Implemented in: R/likelihoods.R (get_POP_like).
At the individual level, the likelihood of observing the parent-offspring pairs (POPs), denoted \Psi, is a product of independent Bernoulli distributions:
L(\Psi \mid \circ) = \prod_{i \in \Theta} \prod_{j \in \Xi} \pi_{ij}^{k_{ij}} (1 - \pi_{ij})^{1 - k_{ij}}
where k_{ij} = 1 if the \{i, j\} pair is a POP and 0 otherwise.
Given that the full comparison set includes millions of \{i, j\} combinations, this base-level likelihood is computationally infeasible. To simplify:
- Group comparisons by:
- Adult capture age a
- Adult capture year y
- Juvenile cohort year c
Let:
\quad \Omega_{c,y,a} be the number of comparisons in group (c, y, a)
\quad \Psi_{c,y,a} be the number of POPs observed in that group
Then the grouped likelihood becomes a binomial likelihood:
L(\Psi_{c,y,a} \mid \cdot) \propto \prod_{c} \prod_{y} \prod_{a} \pi_{c,y,a}^{\Psi_{c,y,a}} (1 - \pi_{c,y,a})^{\Omega_{c,y,a} - \Psi_{c,y,a}}
Standardised Residuals
Implemented status: conceptual diagnostic description; no dedicated close-kin residual helper in current RTMB source.
As in catch-at-age and tagging models, residuals help evaluate model fit.
Under the binomial assumption, the standardised residual for close-kin group (c, y, a) is:
\hat{\varepsilon}_{c,y,a} = \frac{ \Psi_{c,y,a} - \Omega_{c,y,a} \pi_{c,y,a} }{ \sqrt{ \Omega_{c,y,a} \pi_{c,y,a} (1 - \pi_{c,y,a}) } }
If the binomial assumption holds, \hat{\varepsilon}_{c,y,a} should follow a standard normal distribution.
Aggregated Residuals: Higher-Level Approximation
Implemented status: conceptual diagnostic description; no dedicated aggregated-residual helper in current RTMB source.
A key challenge with close-kin data is deciding the appropriate level to evaluate residuals. At the base level, individual parental probabilities \pi_{ij} are often extremely low, with many comparisons and very few POP detections (often only one or two). This structure is uncommon in ecological models, making standard residual analysis (e.g., expecting variance ≈ 1) potentially misleading.
To address this, one more nuanced approach is to aggregate residuals to a higher level — such as cohort c and adult capture year y — by summing across adult ages.
Using the normal approximation to the binomial, the expected number of POPs at this aggregated level is:
\hat{\mu}_{c,y} = \sum_a \Omega_{c,y,a} \pi_{c,y,a}
with corresponding variance:
\hat{\sigma}^2_{c,y} = \sum_a \Omega_{c,y,a} \pi_{c,y,a} (1 - \pi_{c,y,a})
This allows a higher-level standardized residual:
\hat{\varepsilon}_{c,y} = \frac{ \Psi_{c,y} - \hat{\mu}_{c,y} }{ \sqrt{ \hat{\sigma}^2_{c,y} } }
While this normal approximation is useful for visual diagnostics, it is not ideal with very small probabilities or uneven sample sizes. In such cases, it is advisable to compute the exact probability from the sum of binomial variables, especially if assessing how well the OM explains the close-kin data.
The POP likelihood is binomial (pop_switch = 1) or beta-binomial (pop_switch = 2). The beta-binomial concentration pop_od is a model parameter fixed by get_map().
Half-Sibling Pairs (HSP)
Implemented in: R/likelihoods.R (get_HSP_like).
For each HSP data record, the model computes a probability from age-specific reproductive output and survival between cohorts (get_HSP_like):
\gamma_{y,a} = \frac{N_{y,2,a}\phi_{y,a}}{\text{SSB}_y}
and
p^{\text{HSP}}_i = \frac{4q_{\text{hsp}}}{\text{SSB}_{c_{\max}}} \sum_a \gamma_{c_{\min},a}\,\text{Surv}(a,c_{\max}-c_{\min})\,\phi_{c_{\max},a'}
with a binomial (hsp_switch = 1) or beta-binomial (hsp_switch = 2) likelihood. The beta-binomial concentration hsp_od is a model parameter fixed by get_map(). For the binomial case,
K_i \sim \text{Binomial}(N_i, p^{\text{HSP}}_i \times \text{hsp\_false\_negative})
where q_{\text{hsp}} = \exp(\text{par\_log\_hsp\_q}).
Gene Tagging (GT)
Implemented in: R/likelihoods.R (get_GT_like).
For each GT observation (get_GT_like):
p^{\text{GT}}_i = \frac{N^{\text{rel}}_i}{q_{\text{gt}}N_{y^{\text{rel}}_i,1,a^{\text{rel}}_i}}, \quad q_{\text{gt}} = \exp(\text{par\_log\_gt\_q})
Likelihood is binomial (gt_switch = 1) or beta-binomial (gt_switch = 2). The beta-binomial concentration gt_od is a model parameter fixed by get_map().
Prior distributions
Implemented in: R/priors.R (get_priors, evaluate_priors), R/recruitment.R (get_recruitment_prior), R/selectivity.R (get_selectivity_prior), and parameter mapping in R/parameters.R (get_map).
The prior term in the objective is:
lp_{\text{prior}} = -\sum_{m} \log p(\theta_m)
where active priors are defined in get_priors() (and may be changed in code). Current defaults are:
- \text{par\_log\_psi} \sim N(\log(1.75), 0.122^2)
- \text{par\_log\_m10} \sim N(\log(0.1), 0.06^2)
- \text{par\_log\_h} \sim N(\log(1.0), 1.5^2)
- \text{par\_log\_cpue\_omega} \sim N(0.875, 0.1^2)
In addition to lp_prior, two structured priors are always included:
- Recruitment-deviation prior
lp_recfromget_recruitment_prior() - Selectivity prior \sum_f lp_{\text{sel},f} from
get_selectivity_prior()
Parameters can be fixed by mapping in get_map(), which removes them from estimation (and effectively from prior contribution).
Historical ADMB control files also included priors or penalties for steepness, recruitment residual variance, carrying-capacity shifts, selectivity shape, time-varying selectivity changes, and natural mortality. In RTMB these enter through the explicit prior table returned by get_priors(), the recruitment prior in get_recruitment_prior(), and the separable AR1 selectivity prior in get_selectivity_prior(). This means prior assumptions are now controlled by R objects rather than by positional ADMB flags.
MSY Calculation
Implemented in: R/msy.R (msy_calc, find_msy, run_msy).
The current implementation solves an equilibrium two-season, multi-fishery MSY problem using fully-selected fishery mortalities F^{\text{full}}_f.
run_msy() applies this calculation to every retained posterior draw and every fitted catch year by default. For each draw-year cell it uses the year-specific model-predicted catch allocation corresponding to the observed catch equation, the first six normalized fishery selectivities, catch and spawning weights, reproductive output, natural mortality, and stock–recruitment parameters. Chain, within-chain iteration, and calendar year remain explicit in the returned summary; fishery- and age-specific quantities are stored in separate arrays.
MSY State Equations
Implemented in: R/msy.R (msy_calc).
Fishery-at-age mortality is:
F_{f,a} = F^{\text{full}}_f s_{f,a}
Seasonal exploitation rates are accumulated as:
ER_{2,a} = \sum_{f=1}^{2} F_{f,a}, \qquad ER_{1,a} = \sum_{f=3}^{n_f} F_{f,a}
To avoid invalid exploitation rates near 1, each seasonal maximum is smoothed with a posfun-style adjustment. Let x_s = 1 - \max_a(ER_{s,a}) and \epsilon = 0.05. Then:
\text{sr\_adj}_s = \begin{cases} x_s, & x_s \ge \epsilon \\ \dfrac{\epsilon}{2 - x_s/\epsilon}, & x_s < \epsilon \end{cases}
\text{adj}_s = \frac{1-\text{sr\_adj}_s}{\max_a(ER_{s,a})}
and ER_{s,a} and F_{f,a} are rescaled by adj. The penalty contribution is:
\text{pen} = 10^6 \sum_s 0.01\,(x_s-\epsilon)^2 \, \mathbf{1}(x_s < \epsilon)
Relative equilibrium survivorship-per-recruit is:
Rel_0 = 1
Rel_a = Rel_{a-1}(1-ER_{1,a-1})(1-ER_{2,a-1})e^{-M_{a-1}}, \quad a=1,\dots,A
Plus-group correction:
Rel_A \leftarrow \frac{Rel_A}{1-(1-ER_{1,A})(1-ER_{2,A})e^{-M_A}}
Spawning biomass per recruit and equilibrium recruitment:
SBR = \sum_{a=0}^{A} Rel_a \phi_a
Rec = \max\left(1,\; \alpha - \frac{\beta}{SBR}\right)
Biomass metrics:
SSB^{\text{msy}} = Spbio = Rec \cdot SBR
Tbio = Rec \sum_{a=2}^{A} Rel_a \, spwt_a
Season-2-start abundance proxy for fisheries 1–4:
Rel2_a = Rel_a(1-ER_{1,a})e^{-0.5M_a}
Catch by fishery:
C_f = \begin{cases} \sum_a Rec \cdot Rel_a \cdot F_{f,a}\cdot cwt_{f,a}, & f \in \{5,6\} \\ \sum_a Rec \cdot Rel2_a \cdot F_{f,a}\cdot cwt_{f,a}, & f \in \{1,2,3,4\} \end{cases}
C^{\text{msy}} = Ctot = \sum_{f=1}^{n_f} C_f
Estimated fishery catch split:
\widehat{p}_f = \begin{cases} C_f/C^{\text{msy}}, & C^{\text{msy}} > 0\\ 0, & C^{\text{msy}} = 0 \end{cases}
Annual total mortality-at-age reported by the routine:
F^{\text{tot}}_a = -\log\left[(1-ER_{1,a})(1-ER_{2,a})\right]
Unfished spawning biomass per recruit:
Rel^0_0 = 1,\qquad Rel^0_a = Rel^0_{a-1}e^{-M_{a-1}}
Rel^0_A \leftarrow \frac{Rel^0_A}{1-e^{-M_A}},\qquad SBR_0 = \sum_{a=0}^{A} Rel^0_a \phi_a
MSY Optimization Objective
Implemented in: R/msy.R (find_msy).
find_msy() estimates F^{\text{full}} using bounded nlminb:
\min_{F^{\text{full}}_f \in [10^{-5},\,0.5]} \left[ -\log(C^{\text{MSY}}+1) + \lambda \sum_f (\widehat{p}_f - p_f^{\text{target}})^2 + \text{pen} \right]
The historical two-stage solve is retained. Fisheries with annual catch share no greater than 0.002 are fixed at zero. Active fisheries first use \lambda=100 from the year-specific F_{\text{ini}}, followed by the final \lambda=10^5 solve. The second phase uses \log(F^{\text{full}}) as the optimizer coordinate while retaining the same [10^{-5},0.5] bounds and objective. This reparameterization avoids the severe numerical stiffness of the final allocation penalty without changing the MSY solution.
Projection Model
Implemented in: R/projections.R (run_projections, project_selectivity, project_rec_devs, project_ctp_schedule) with core dynamics from R/dynamics.R (do_dynamics).
Model Structure
Implemented in: R/projections.R (run_projections) through projected catch/selectivity inputs to the dynamics engine.
The projection dynamics retain all six conditioning-model fisheries: LL1, LL2, LL3, LL4, Indonesia, and Australia. run_projections() constructs a year-by-season-by-fishery catch array. Observed catches are retained for years covered by the conditioning data; later fixed or CTP-derived TACs are allocated across all six fisheries using projection_tac_split_yf (or a repeated projection_tac_split_f). When no split is supplied, fishery proportions are calculated from the most recent five observed catch years. Each fishery’s TAC is then divided between seasons using its proportions in the final observed catch year. Fishery allocation rows are required to sum to one.
Population Model
Implemented in: R/projections.R (run_projections) via calls to R/dynamics.R (do_dynamics).
Identical to the conditioning model.
Initial Abundances
Implemented in: R/projections.R (run_projections) using terminal reconstruction state and projected recruitment deviates from project_rec_devs.
Initial age-specific abundances at the start of the projection are estimated from the conditioning model. To represent process error in recruitment, lognormal autocorrelated error is added to initial abundances for ages 0–2.
Let Y_1 be the first projection year. Then:
N_{Y_1,4} = \hat{N}_{Y_1,4} \cdot \exp(0.4z - 0.08)
N_{Y_1,3} = \hat{N}_{Y_1,3} \cdot \exp(0.4z - 0.08)
N_{Y_1,2} = \hat{N}_{Y_1,2} \cdot \exp(\epsilon_{Y_1})
N_{Y_1,1} = \hat{N}_{Y_1,1} \cdot \exp(\hat{\rho} \epsilon_{Y_1-2} + \epsilon_{Y_1-1})
N_{Y_1,0} = \hat{N}_{Y_1,0} \cdot \exp(\hat{\rho}^2 \epsilon_{Y_1-2} + \hat{\rho} \epsilon_{Y_1-1} + \epsilon_{Y_1})
Where:
-
z \sim {N}(0,1)
-
\epsilon_y \sim {N}\left(0, (1 - \hat{\rho}^2) \sigma_R^2\right)
- \sigma_R = 0.6 (and an extra lognormal SD = 0.4 is added for \hat{N}_{Y_1-4,0} and \hat{N}_{Y_1-3,0})
These imply:
\tau_{Y_1} = \hat{\tau}_{Y_1} + \hat{\rho}^2 \epsilon_{Y_1-2} + \hat{\rho} \epsilon_{Y_1-1} + \epsilon_{Y_1}
and subsequent years follow:
\tau_{Y_1+1} = \hat{\rho} \tau_{Y_1} + \epsilon_{Y_1+1}, \quad \text{etc.}
This formulation assumes autocorrelated recruitment starting in Y_1-2.
However, it may over-propagate uncertainty from Y_1-3 due to grid-based point estimates.
A sensitivity run with uncorrelated \tau_{Y_1-2} is used to address this.
Selectivity in Projections
Implemented in: R/projections.R (project_selectivity).
Unlike the conditioning model, the projection model does not use random-walk selectivity due to the risk of unrealistic drift. The default projection option sets future selectivity to the mean of the last 10 historical years before the first projection year:
For fishery f and age a:
s^{\text{proj}}_{f,y,a} = \frac{1}{10} \sum_{j=1}^{10} s_{f,Y_0-j,a}, \quad y = Y_0, \ldots, Y_P
where Y_0 is the first projection year and Y_P is the final projection year. The same mean selectivity vector is repeated for every projected year.
An older lognormal option is still available in project_selectivity() via option = "lognormal", where each fishery-age selectivity is drawn from the historical log-scale mean and standard deviation over the selected sample years.
Australian Surface Fishery (Special Case)
Implemented status: legacy/ADMB-style narrative; no dedicated age-3 boost rule in current RTMB projection helpers (R/projections.R).
Let:
P_{y,3} = \frac{N_{y,3}}{\sum_{a=1}^5 N_{y,a}}, \quad \bar{P}_3 = \frac{1}{10} \sum_{y=y_{n2}-9}^{y_{n2}} P_{y,3}
If P_{y,3} \ge \bar{P}_3:
s_{6,y,a} = s_{6,y-1,a} \cdot \exp(\varepsilon_{6,y,a}) \quad \text{for } a = 1, \ldots, 5, \quad \varepsilon_{6,y,a} \sim {N}(0, 0.1^2)
Otherwise, boost selectivity at age 3:
s_{6,y,3} = s_{6,y-1,3} \cdot \left(1 + 0.5 \cdot \frac{\bar{P}_3 - P_{y,3}}{P_{y,3}} \right), \quad s_{6,y,a} = s_{6,y-1,a} \cdot \exp(\varepsilon_{6,y,a}) \text{ for } a = 1,2,4,5
Fishing Mortality in Projections
Implemented in: R/projections.R (run_projections) via R/dynamics.R (get_harvest_rate, do_dynamics).
Year subscripts are omitted for simplicity. Catch is defined as:
C = \sum_a \sum_f s_{f,a} F_f N_a
For each fishery f:
C = \sum_f C_f, \quad \text{where } C_f = \left( \sum_a s_{f,a} N_a \right) F_f
Solving for F_f:
F_f = \frac{C_f}{\sum_a s_{f,a} N_a}
Additionally:
- Catch-at-age: C_{f,a} = s_{f,a} F_f N_a
- Catch at age: C_a = \left( \sum_f s_{f,a} F_f \right) N_a
Note: A problem occurs if \sum_f s_{f,a} F_f > 1, which implies C_a > N_a.
Harvest Rate Bounding and Selectivity Adjustment
Implemented status: partial. Current RTMB uses posfun penalty in R/dynamics.R (get_harvest_rate); explicit A1-A6 rescaling routine is legacy narrative.
In early MP trials, the model enforced a hard upper bound on age-specific exploitation rates:
\sum_f s_{f,a} F_f \le 0.99
When this bound was exceeded, the model reduced the catch at that age, but not the selectivity or fishing mortality for other ages, which led to underutilization of TAC.
An improved method adjusts selectivities at the age(s) where the bound is exceeded (from MP Workshop II, Queenstown, April 2003).
Case of One Fleet (or Non-Overlapping Selectivities)
Implemented status: legacy/ADMB analytical description; no dedicated standalone helper in current RTMB source.
Assume a single fleet. If F \le 0.9, no adjustment is needed. If F > 0.9, apply:
C = \sum_a g(s_a F) N_a \tag{A1}
where the adjusted selectivity is:
s_a^* = \frac{g(s_a F)}{F} \tag{A2}
The proposed function g(x) is:
g(x) = \begin{cases} x & \text{if } x \le 0.9 \\ 0.9 + 0.08 \left[ 1 - \exp(-10(x - 0.9)) \right] & \text{if } x > 0.9 \end{cases} \tag{A3}
Notes: - g(x) < 1 ensures C_a = g(s_a F) N_a < N_a - g(x) is continuous and differentiable at x = 0.9 - Newton-Raphson is used to solve equation (A1) for F
Extension to Multiple Fleets
Implemented status: partially represented by multi-fleet harvest in R/dynamics.R (get_harvest_rate); the explicit A4-A6 derivation is legacy documentation.
If:
\sum_f s_{f,a} F_f \le 0.9 \quad \text{for all } a
then proceed as normal. Otherwise, for any age a:
C = \sum_a g\left(\sum_f s_{f,a} F_f\right) N_a \tag{A4}
To proportionally reduce selectivity for each fleet at age a, define:
s^*_{f,a} = s_{f,a} \cdot \left[ \frac{g\left( \sum_f s_{f,a} F_f \right)}{\sum_f s_{f,a} F_f} \right] \tag{A5}
Then:
C_{f,a} = s^*_{f,a} F_f N_a = s_{f,a} F_f N_a \cdot \left[ \frac{g\left( \sum_f s_{f,a} F_f \right)}{\sum_f s_{f,a} F_f} \right] \tag{A6}
This guarantees:
C_a = \sum_f C_{f,a} = g\left( \sum_f s_{f,a} F_f \right) N_a < N_a
Solving for F_f When Bounds Are Exceeded
Implemented status: current RTMB computes F_f algebraically by fishery in R/dynamics.R (get_harvest_rate) and applies posfun penalty rather than an explicit multivariate root-solver.
Once selectivity is adjusted via:
C_a = \sum_f C_{f,a} = \sum_f s^*_{f,a} F_f N_a = g\left(\sum_f s_{f,a} F_f \right) N_a
Then the coupled non-linear equations must be solved for each F_f:
C_f = \sum_a C_{f,a} = \sum_a s^*_{f,a} F_f N_a
Substituting s^*_{f,a} gives:
C_f = \sum_a s_{f,a} F_f \left[ \frac{g\left( \sum_{f'} s_{f',a} F_{f'} \right)}{\sum_{f'} s_{f',a} F_{f'}} \right] N_a
This system is solved via a multivariate root-finding method, e.g., extended Newton-Raphson.
CPUE and Aerial Survey Deviations
Implemented status: partial. CPUE projection observations are simulated in R/projections.R (run_projections) using R/likelihoods.R (get_cpue_like); aerial-deviation text below is legacy narrative.
CPUE Deviations
Implemented in: R/projections.R (run_projections) plus CPUE predictor from R/likelihoods.R (get_cpue_like).
Expected CPUE is calculated from projected abundance, projected CPUE selectivity, and the fitted CPUE parameters in get_cpue_like(): catchability block, abundance-power parameter, and fitted creep rate. The projection is scaled consistently with the historical unscaled CPUE series. Future observed CPUE is then generated as
\log I_y^{\mathrm{obs}} = \log I_y^{\mathrm{pred}} + \epsilon_y,
with the fitted-history AR(1)
\epsilon_y = \rho_I \epsilon_{y-1} + \sqrt{1-\rho_I^2}\,\sigma_{\mathrm{sim},y} z_y, \qquad z_y \sim N(0,1), \qquad \sigma_{\mathrm{sim},y} = \sqrt{\sigma_{I,y}^2 + \sigma_{\mathrm{cpue}}^2},
where \sigma_{I,y} is the projected sampling-error SD derived from the CPUE input and \sigma_{\mathrm{cpue}} is the fitted additional-error SD. For each posterior draw, \rho_I is the empirical lag-one correlation of the historical fitted log residuals, \log I_t^{\mathrm{obs}}-\log I_t^{\mathrm{pred}}. The recursion is initialized with the final historical fitted residual, so the first future CPUE observation is conditioned directly on the last fitted CPUE year. The \sqrt{1-\rho_I^2} term follows the legacy ADMB projection and treats \sigma_{\mathrm{sim},y} as the marginal residual SD. The process is advanced through every intervening calendar year if projection monitoring begins after a gap. Observations outside the configured monitoring years are stored as missing.
This is autocorrelation in the projected CPUE observation residual, not a separate latent catchability process. There is no separate fixed 0.5% increase; any trend in expected CPUE comes from the fitted par_cpue_creep value.
Aerial Survey Deviations
Implemented status: legacy narrative; no dedicated aerial-deviation simulation helper in current R/projections.R.
Aerial indices are simulated by adding lognormal deviations. Initially, empirical variance and autocorrelation were used, but:
- When \tau_{\text{aerial}} was estimated as free, it produced high SDs (0.50–0.60)
- Residual autocorrelation was very low
To stabilize, \tau_{\text{aerial}} was fixed to 0.18, which with sampling error gives SD \approx 0.30 in log space.
Simulation Grid for Deviation Parameters
Implemented status: scenario documentation; no dedicated grid helper for this table in current RTMB code.
| Scenario | SD (CPUE) | \rho_{\text{cpue}} |
|---|---|---|
| Base | 0.20 | Empirical, correlated to last residual |
| highCpueCV | 0.30 | Empirical, correlated to last residual |
For the aerial survey, autocorrelation is set to zero due to short residual history and poor reliability.
Lags in Data Availability and Schedules of TAC Changes
Implemented status: policy/process documentation; not encoded as executable scheduling logic in current RTMB source.
When a TAC for year y is determined, the projection code assumes the following data are available:
- Catch data up to year y - 3
- TAC up to year y - 1
- CPUE up to year y - 3
- Age-composition data up to year y - 3
- Gene tagging index (age-2 abundance) in year y - 4
- Pop-HSP index of SSB for year y - 7
Options for TAC Update Frequency
Implemented status: policy documentation; not represented by a dedicated function in R/*.R.
OMMP8 recommended calculating the first TAC under the new MP in 2020 (for implementation in 2021), to allow more time for stakeholder consultation. This skips the usual lag for only the first decision.
The control file includes four TAC update schedule options:
-
Option (a): First TAC in 2021, then every year — no lag
-
Option (b): First TAC in 2021, then every other year — no lag
-
Option (c): First TAC in 2021, then every 3 years — no lag
- Option (d): First TAC in 2021, then every 3 years — with one-year lag
Example Schedule: Options (c) and (d)
Implemented status: illustrative schedule table; not generated by a dedicated function in current RTMB source.
In all cases, the first TAC is calculated in 2019, using:
- CPUE: up to 2018
- Gene tagging: 2016, 2017
- POP/HSP SSB: up to 2014
| Decision Year | Catch Data | CPUE Data | Gene Tagging | POP/HSP SSB | TAC Year | TAC Change? |
|---|---|---|---|---|---|---|
| 2017 | 2016 | 2016 | — | 2012 | 2019 | hardwired |
| 2018 | 2017 | 2017 | 2016 | 2013 | 2020 | hardwired |
| 2019 | 2018 | 2018 | 2017 | 2014 | 2021 | Yes |
| 2020 | 2019 | 2019 | 2018 | 2015 | 2020 | No |
| 2021 | 2020 | 2020 | 2019 | 2016 | 2021 | No |
| 2022 | 2021 | 2021 | 2020 | 2017 | 2022 | Yes |
| 2023 | 2022 | 2022 | 2021 | 2018 | 2023 | No |
Note: Shaded years (2016–2018) represent real data available in 2019 (not OM-simulated).
Tuning Levels
Implemented status: policy/tuning documentation; not encoded in a dedicated tuning module in current RTMB source.
The SFMWG proposed six tuning options based on the probability of rebuilding the stock to 20% of SSB_0 in 25 or 30 years. Short-term checkpoints were also set at 12 or 15 years from present, to check the probability of rebuilding to 10% of SSB_0 or to 2 \times SSB_{2009}.
Tuning Options and Checkpoint Years
Implemented status: documentation table; not generated by a dedicated function in current RTMB source.
| Tuning Level | Tuning Year | P[SSB ≥ 0.2 SSB_0] | Short-term Checkpoint Year |
|---|---|---|---|
| 1 | 2035 | 0.60 | 2022 |
| 2 | 2035 | 0.70 | 2022 |
| 3 | 2035 | 0.90 | 2022 |
| 4 | 2040 | 0.60 | 2025 |
| 5 | 2040 | 0.70 | 2025 |
| 6 | 2040 | 0.90 | 2025 |
Maximum–Minimum TAC Changes
Implemented status: policy constraint discussion; bounds can be explored by user-specified projection inputs, not by a dedicated hardcoded helper.
The SFMWG recommended examining two options for TAC change bounds:
- Maximum change: 3000 t or 5000 t
- Minimum change: 100 t
These are not hardcoded, so users can explore alternative values.
However, prior analyses indicated that these bounds act as binding constraints, especially at high tuning levels where the TAC was often pushed against the maximum limit.
Appendix: Transition from ADMB to RTMB
Implemented in: R/parameters.R (get_map), R/model.R (sbt_model), and the staged fit lifecycle in R/sbt-fit.R and R/fit-lifecycle.R.
This model was originally implemented in ADMB using control file flags. These have been replaced in RTMB with mapped parameters and script-based control. Key transitions include:
| ADMB Flag/Control Concept | RTMB Equivalent |
|---|---|
flag_rec_devs = 1 |
Estimate par_rdev_y; use sbt_add_map() to fix selected entries |
flag_est_q = 0 |
Map the relevant catchability parameter, such as par_log_cpue_q, to factor(NA)
|
| Estimation phase definitions | Configure the fit in R, then call sbt_optimise()
|
| Profile and Bayesian runs | Use sbtprofile() and sbt_mcmc() respectively |
This makes model setup portable and scriptable, and separates scientific configuration from optimisation and sampling controls.
Appendix: Mathematical notation to R code cross-reference
This table maps the mathematical symbols used in the equations above to their corresponding variable names in the R implementation.
Population dynamics
| Mathematical symbol | Description | R code variable | Location |
|---|---|---|---|
| N_{y,s,a} | Numbers-at-age | number_ysa[y, s, a] |
R/dynamics.R |
| N^*_{y,a} | Mid-year numbers | number_ysa[y, 2, a] |
R/dynamics.R |
| M_a | Natural mortality at age | M_a |
R/natural-mortality.R |
| R_y | Recruitment in year y | recruitment_y[y] |
R/dynamics.R |
| S_y | Spawning biomass | spawning_biomass_y[y] |
R/dynamics.R |
| H_{y,s,f,a} | Harvest rate | hrate_ysfa[y, s, f, a] |
R/dynamics.R |
| F_{y,s,f} | Fishing mortality | F_ysf[y, s, f] |
R/dynamics.R |
| s_{f,y,a} | Selectivity | sel_fya[f, y, a] |
R/selectivity.R |
| w_{f,y,a} | Weight-at-age (tonnes per fish) | weight_fya[f, y, a] |
data input |
| C_{y,s,f} | Catch | catch_obs_ysf[y, s, f] |
data input |
Stock-recruitment
| Mathematical symbol | Description | R code variable | Location |
|---|---|---|---|
| \alpha | Beverton-Holt alpha parameter | alpha |
R/dynamics.R |
| \beta | Beverton-Holt beta parameter | beta |
R/dynamics.R |
| B_0 | Unfished spawning biomass | B0 |
R/model.R |
| R_0 | Unfished recruitment | R0 |
R/dynamics.R |
| h | Steepness | par_h |
R/model.R |
| \sigma_R | Recruitment SD | sigma_r |
R/model.R |
| \epsilon_y | Recruitment deviation | rdev_y[y] |
R/model.R |
| \rho | Autocorrelation | rho |
R/recruitment.R |
| \nu | Depensation parameter | sr_dep |
R/recruitment.R |
| \psi | Fecundity exponent parameter |
par_log_psi, psi
|
R/model.R |
| \phi_{y,a} | Spawning output-at-age | phi_ya[y, a] |
R/dynamics.R |
Natural mortality
| Mathematical symbol | Description | R code variable | Location |
|---|---|---|---|
| m^1 | M at young ages |
par_log_m0 / M_a
|
R/model.R, R/natural-mortality.R
|
| m^4 | M at age 4 |
par_log_m4 / M_a
|
R/model.R, R/natural-mortality.R
|
| m^{10} | M at age 10 |
par_log_m10 / M_a
|
R/model.R, R/natural-mortality.R
|
| m^{30} | M at age 30+ |
par_log_m30 / M_a
|
R/model.R, R/natural-mortality.R
|
| p | M slope parameter | M_slope |
R/natural-mortality.R |
Tagging
| Mathematical symbol | Description | R code variable | Location |
|---|---|---|---|
| N_{k,t,i} | Tag releases | tag_release_cta[k, t, i] |
data input |
| R_{k,t,i,j} | Tag recaptures | tag_recap_ctaa[k, t, i, j] |
data input |
| \xi_t | Immediate tag retention | shed1[t] |
data input |
| \Omega_t | Continuous shedding rate | shed2[t] |
data input |
| \nu_{y,a} | Tag reporting rate | tag_rep_rates_ya[y, a] |
data input |
| \phi_{\text{tag}} | Tag overdispersion | tag_var_factor |
data input |
| H_{\text{tag}} | Tagged-fish harvest multiplier | tag_H_factor |
parameter input |
Close-kin genetics
| Mathematical symbol | Description | R code variable | Location |
|---|---|---|---|
| \pi_i | Parental probability | pop_prob[i] |
R/likelihoods.R |
| \phi_{y,a} | Reproductive output | phi_ya[y, a] |
R/dynamics.R |
| POP counts | Parent-offspring pairs | pop_obs |
data input |
| HSP probability | Half-sibling probability | hsp_prob[i] |
R/likelihoods.R |
CPUE and surveys
| Mathematical symbol | Description | R code variable | Location |
|---|---|---|---|
| q | Catchability coefficient | cpue_q |
R/likelihoods.R |
| CPUE creep | Annual q increase | cpue_creep |
R/likelihoods.R |
| \omega | CPUE power parameter | cpue_omega |
R/likelihoods.R |
| \sigma_{CPUE} | CPUE observation error | cpue_sigma |
R/likelihoods.R |
| Aerial survey q | Aerial catchability | exp(par_log_aerial_tau) |
R/likelihoods.R |
Composition data
| Mathematical symbol | Description | R code variable | Location |
|---|---|---|---|
| Age comp | Observed age frequencies | af_obs |
data input |
| Length comp | Observed length frequencies | lf_obs |
data input |
| Sliced age comp | LF-derived age frequencies | af_sliced_ysfa |
data input |
| ALK | Age-length key | alk_ysal |
data input |
| \alpha_{AF} | Age comp overdispersion | exp(par_log_af_alpha) |
R/model.R |
| \alpha_{LF} | Length comp overdispersion | exp(par_log_lf_alpha) |
R/model.R |
Many parameters are estimated on the log scale (for example par_log_B0) and transformed inside sbt_model() using exp() to ensure positivity.