2a. Mathematical description of the RCM

Quang Huynh (quang@bluematterscience.com)

2021-03-09


1 Selectivity and mortality

For fleets, selectivity is defined by blocks (indexed by \(b\)) which can then be assigned to any fleet \(f\) and year \(y\) to allow for time-varying selectivity. By default, each fleet is assigned to its own block for all years (no time-varying selectivity).

For flat-topped selectivity in block \(b\), two parameters are used and expressed in terms of length units: the length of 5% selectivity (\(L^5_b\)) and the length of full selectivity \(L^{\textrm{FS}}_b\). For dome selectivity, a third parameter, the selectivity at \(L_{\infty}\), \(V^{L_{\infty}}_b\) is also used. Length-based selectivity is converted to age-based selectivity in the age-structured model as:

\[ v_{y,a,b} = \begin{cases} 2^{-[(L_{y,a} - L^{\textrm{FS}}_b)/(\sigma^{\textrm{asc}}_b)]^2} & \textrm{if } L_{y,a} < L^{\textrm{FS}}_b\\ 1 & \textrm{if logistic and } L_{y,a} \ge L^{\textrm{FS}}_b,\\ 2^{-[(L_{y,a} - L^{\textrm{FS}}_b)/(\sigma^{\textrm{des}}_b)]^2} & \textrm{if dome and } L_{y,a} \ge L^{\textrm{FS}}_b \end{cases} \]

where \(L_{y,a}\) is the mean length-at-age and \(\sigma^{\textrm{asc}}_b = (L^5_b - L^{\textrm{FS}}_b)/\sqrt{-\log_2(0.05)}\) and \(\sigma^{\textrm{des}}_b = (L_{\infty} - L^{\textrm{FS}}_b)/\sqrt{-\log_2(V^{L_{\infty}}_b)}\) control the shape of the ascending and descending limbs, respectively, of the selectivity function. In this parameterization, length-based selectivity is constant over time. The corresponding age-based selectivity within each block is constant so long as growth is not time-varying.

Selectivity can also be parameterized where \(v_{y,a,b}\) are free independent parameters. In this case, selectivity does not vary among years.

Total mortality \(Z\) in year \(y\) and for age \(a\) is the sum of fishing mortality \(F\) from all fleets and natural mortality \(M\),

\[ Z_{y,a} = M_{y,a} + \Sigma_f v_{y,a,f} F_{y,f},\]

where \(v_{y,a,f}\) is the fleet selectivity after assigning blocks to fleets.

2 Survey selectivity

Survey selectivity is constant over time and is denoted as \(v_{a,s}\) with either logistic, dome, or free parameterizations.

3 Initial population distribution

The population age distribution in the first year of the model \(y=1\) is in equilibrium where \[ N_{1,a} = \begin{cases} R^{\textrm{eq}} \exp(-\Sigma_{i=1}^{a-1}Z^{\textrm{eq}}_i) & a = 1, \ldots, A-1\\ \dfrac{R^{\textrm{eq}} \exp(-\Sigma_{i=1}^{a-1}Z^{\textrm{eq}}_i)}{1 - \exp(-Z^{\textrm{eq}}_A)} & a = A, \end{cases} \] where the \(R^{\textrm{eq}}\) is the equilibrium recruitment and \(Z^{\textrm{eq}}_a = M_{1,a} + \Sigma_f v_{1,a,f} F^{\textrm{eq}}_f\) is the equilibrium total mortality rate. Unfished conditions are modeled by setting \(F^{\textrm{eq}}_f = 0\). To estimate \(F^{\textrm{eq}}_f\), the corresponding equilibrium catch in weight \(Y^{\textrm{eq}}_f\) prior to the first year of the model should be provided. In the equilibrium yield curve, \(F^{\textrm{eq}}_f\) would be the fishing mortality corresponding to fishing at \(F^{\textrm{eq}}_f\). Once \(Z^{\textrm{eq}}_a\) is obtained, then the equilibrium recruitment is calculated as:

\[ R^{\textrm{eq}} = \begin{cases} \dfrac{\alpha^{\textrm{B}}\phi^{\textrm{eq}} - 1}{\beta^{\textrm{B}}\phi^{\textrm{eq}}} & \textrm{if Beverton-Holt stock-recruit relationship}\\ \dfrac{\log(\alpha^{\textrm{R}}\phi^{\textrm{eq}})}{\beta^{\textrm{R}}\phi^{\textrm{eq}}} & \textrm{if Ricker stock-recruit relationship} \end{cases}, \] where \(\phi^{\textrm{eq}}\) is the spawners-per-recruit when the mortality is \(Z^{\textrm{eq}}_a\). From steepness \(h\), \(\alpha^{\textrm{B}} = \frac{4h}{(1-h)\phi_0}\), \(\beta^{\textrm{B}} = \frac{5h-1}{(1-h)B^S_0}\), \(\alpha^{\textrm{R}} = \frac{(5h)^{1.25}}{\phi_0}\), \(\beta^{\textrm{R}} = \frac{\log(5h)}{B^S_0}\), where \(\phi_0\) and \(B^S_0\) are unfished spawners-per-recruit and unfished spawning biomass, respectively.

4 Dynamics equations

After setting the equilibrium population age distribution in the first year of the model, the population abundance \(N_{y,a}\) in subsequent years is \[ N_{y,a} = \begin{cases} R_y & a = 1\\ N_{y-1,a-1} \exp(-Z_{y-1,a-1}) & a = 2, \ldots, A - 1,\\ N_{y-1,a-1} \exp(-Z_{y-1,a-1}) + N_{y-1,a} \exp(-Z_{y-1,a}) & a = A \end{cases} \] where \(R_y\) is the recruitment and \(A\) is the maximum-age as the plus-group. Recruitment is modelled as \[ R_y = \begin{cases} \dfrac{\alpha^{\textrm{BH}} B^S_{y-1}}{1 + \beta^{\textrm{BH}}B^S_{y-1}} \exp(\delta_y - 0.5 \tau^2) & \textrm{if Beverton-Holt stock-recruit relationship}\\ \alpha^{\textrm{Ricker}} B^S_{y-1} \exp(-\beta^{\textrm{Ricker}} B^S_{y-1})\exp(\delta_y - 0.5 \tau^2) & \textrm{if Ricker stock-recruit relationship} \end{cases}, \] where \(\delta_y\) are recruitment deviates and \(\tau\) is the standard deviation of the deviates.

The spawning biomass is \(B^S_y\) is \[B^S_y = \sum_a w_{y,a} m_{y,a} N_{y,a},\] where \(m_{y,a}\) and \(w_{y,a}\) are the maturity at age and weight at age, respectively.

5 Catch at age

The catch (in numbers) \(C\) at age for fleet \(f\) is \[ C_{y,a,f} = \dfrac{v_{y,a,f} F_{y,f}}{Z_{y,a}} N_{y,a} (1 - \exp[-Z_{y,a}]).\]

If the model is conditioned on catch, then \(F_{y,f}\) can be estimated as parameters or solved iteratively to match the observed catch. If the model is conditioned on effort, then \[ F_{y,f} = q_f E_{y,f},\] where \(E_{y,f}\) is the observed effort and \(q_f\) is the scaling coefficient.

The proportion of the catch-at-age is \[ p_{y,a,f} = \dfrac{C_{y,a,f}}{\sum_a C_{y,a,f}}.\]

6 Catch-at-length

The catch at length is calculated assuming a normally distributed length-at-age \(P(\ell,a)\), where \[ C_{y,\ell,f} = \sum_a C_{y,a,f} P(\ell|a) \] and

\[ P(\ell|a) = \begin{cases} \phi(L'_{\ell+1}) & \ell = 1\\ \phi(L'_{\ell+1}) - \phi(L'_\ell) & \ell = 2, \ldots, L - 1,\\ 1 -\phi(L'_\ell) & \ell = L \end{cases} \] with \(L'_{\ell}\) as the length at the lower boundary of length bin \(\ell\) and \(\phi(L'_{\ell})\) as the cumulative distribution function of a normal variable with mean \(\tilde{L}_{y,a}\) (the expected mean length at age \(a\)) and standard deviation \(\tilde{L}_{y,a} \times CV^L\) (\(CV^L\) is the coefficient of variation in mean length at age).

The proportion of the catch-at-length is \[ p_{y,\ell,f} = \dfrac{C_{y,\ell,f}}{\sum_{\ell}C_{y,\ell,f}}.\]

7 Catch in weight

The catch in weight \(Y\) is \[ Y_{y,f} = \sum_a C_{y,a,f} w_{y,a}.\]

8 Mean size

The mean length of the catch \(\bar{L}_{y,f}\) is \[ \bar{L}_{y,f} = \dfrac{\sum_{\ell} L_{\ell} C_{y,\ell,f}}{\sum_{\ell} C_{y,\ell,f}},\] where \(L_\ell\) is the midpoint of the length bin \(\ell\).

The mean weight of the catch \(\bar{w}_{y,f}\) is \[ \bar{w}_{y,f} = \dfrac{\sum_a C_{y,a,f}w_{y,a}}{\sum_a C_{y,a,f}},\]

9 Survey

If the \(s^{\textrm{th}}\) survey is biomass-based, then the survey value \(I_{y,s}\) is calculated as \[ I_{y,s} = q_s \sum_a v_{y,a,s} N_{y,a} w_{y,a}, \] where \(q\) is the scaling coefficient and \(s\) indexes survey.

If the survey is abundance-based, then \[ I_{y,s} = q_s \sum_a v_{y,a,s} N_{y,a} . \]

The proportion-at-age vulnerable to the survey is \[ p_{y,a,s} = \dfrac{v_{a,s}N_{y,a}}{\sum_a v_{a,s}N_{y,a}}.\]

The proportion-at-length vulnerable to the survey is \[ p_{y,\ell,s} = \dfrac{\sum_a v_{a,s} N_{y,a} P(\ell|a)}{\sum_{\ell} \sum_a v_{a,s} N_{y,a} P(\ell|a)}.\]

10 Estimated parameters

The estimated parameters, denoted in this section as \(x\), are unconstrained over all real numbers and then transformed in order to constrain the corresponding model parameters. For optimization, the transformation is also designed to reduce the scale of all estimated parameters to within an order of magnitude.

10.1 Selectivity

For a fleet block \(b\) for which selectivity is estimated, then parameters \(x^{LFS}_b\) and \(x^{L5}_b\) are estimated over all real numbers, where

\[ \begin{align} L^{\textrm{FS}}_b &= 0.99 \times L_{\infty} \times \textrm{logit}^{-1}(x^{LFS}_b)\\ L^5_b &= L^{\textrm{FS}}_b - \exp(x^{L5}_b) \end{align}\]

If a third parameter \(x^{V}_b\)is estimated for dome selectivity, then \[ V^{L_{\infty}}_b = \textrm{logit}^{-1}(x^V_b)\]

If selectivity is parameterized as free parameters, then \[ v_{y,a,f} = \textrm{logit}^{-1}(x^v_{y,a,f}).\]

For surveys, parameterizations are identical except with indexing for survey \(s\).

10.2 Fishing mortality

If \(F_{y,f}\) are estimated parameters (condition = "catch"), then one parameter \(x^F_f\) is the estimated \(F\) in log-space in the middle of the time series is estimated and all others are subsequent deviations, represented as \(x^{F_{dev}}_{y,f}\):

\[ F_{y,f} = \begin{cases} \exp(x^F_f) & y \textrm{ is midpoint of the time series}\\ \exp(x^F_f) \times \exp(x^{F_{dev}}_{y,f}) & \textrm{otherwise}\\ \end{cases} \]

If condition = "effort", then \(q_f\) is estimated in log space, where \[F_{y,f} = q_f E_{y,f} = \exp(x^q_f) \times E_{y,f}\]

10.3 Index catchability

To scale biomass to index values, the index catchability \(q_s\) is solved analytically in the model:

\[ q_s = \exp\left(\dfrac{\sum_y \log(I^{\textrm{obs}}_{y,s}) - \sum_y \log(\sum_a v_{y,a,s}N_{y,a,s})}{n_s}\right),\] or \[ q_s = \exp\left(\dfrac{\sum_y \log(I^{\textrm{obs}}_{y,s}) - \sum_y \log(\sum_a v_{y,a,s}N_{y,a,s}w_{y,a})}{n_s}\right),\] for an abundance or biomass based index, respectively, where \(n_s\) is the number of years with index values and the summation is over those \(n_s\) years.

10.4 Other parameters

Unfished recruitment is estimated in log-space, \(R_0 = \dfrac{1}{z}\exp(x^{R_0})\) where \(z\) is an optional rescaler, e.g. mean historical catch, to reduce the magnitude of the \(x^{R_0}\) estimate.

Recruitment deviations \(\delta_y\) are estimated in log space.

The support of the steepness parameter \(x^h\) is over all real numbers and is transformed. With the Beverton-Holt SR function: \[ h = 0.8 \times \dfrac{1}{1 + \exp(-x^h)} + 0.2. \] With a Ricker SR function: \[ h = \exp(x^h) + 0.2.\]

Steepness is fixed unless a prior is used.

Age- and time-constant natural mortality \(M\) can be estimated with a prior. Otherwise, it is fixed to values specified in the OM object. The parameter \(x^M\) is estimated in log space.

11 Objective function

11.1 Total objective function

The total objective function \(\textrm{LL}\) to be maximized is \[\textrm{LL} = \sum_{i=1}^9\Lambda_i + \sum \textrm{pr} + \sum_y\sum_f\textrm{pen}_{y,f},\]

including the log likelihoods \(\Lambda_i\), logarithm of the parameter priors \(pr\), and penalty functions, described below.

11.2 Likelihoods

If the model is conditioned on catch and fishing mortality rates are estimated parameters, then the log-likelihood component \(\Lambda_1\) of the catch is \[\Lambda_1 = \sum_f \left[\lambda^{Y}_f \sum_y \left(-\log(\omega_{y,f}) - \dfrac{[\log(Y^{\textrm{obs}}_{y,f}) - \log(Y^{\textrm{pred}}_{y,f})]^2}{2 \omega _{y,f}^2}\right)\right],\]

where \(\textrm{obs}\) and \(\textrm{pred}\) indicate observed and predicted quantities, respectively, \(\lambda\) are likelihood weights, and \(\omega\) is the catch standard deviation. With a small standard deviation for the catch likelihood relative to the variance in other likelihood components, i.e., \(\omega = 0.01\), the predicted catch should match the observed catch in the model.

The log-likelihood component \(\Lambda_2\) of survey data is \[\Lambda_2 = \sum_s \left[ \lambda^I_s \sum_y \left(-\log(\sigma_{y,s}) - \dfrac{[\log(I^{\textrm{obs}}_{y,s}) - \log(I^{\textrm{pred}}_{y,s})]^2}{2\sigma_{y,s}^2}\right) \right].\]

The log-likelihood component \(\Lambda_3\) of catch-at-age data is \[\Lambda_3 = \sum_f \lambda^A_f \left[\sum_y O^A_{y,f} \sum_a p^{\textrm{obs}}_{y,a,f} \log(p^{\textrm{pred}}_{y,a,f})\right],\] where \(O^A\) is the annual sample sizes for the age compositions.

The log-likelihood component \(\Lambda_4\) of catch-at-length data is \[\Lambda_4 = \sum_f \lambda^L_f \left[ \sum_y O^L_{y,f} \sum_{\ell} p^{\textrm{obs}}_{y,\ell,f} \log(p^{\textrm{pred}}_{y,\ell,f})\right]\] where \(O^L\) is the annual sample sizes for the length compositions.

The log-likelihood component \(\Lambda_5\) for the observed mean size in this catch is:

\[\Lambda_5 = \sum_f \lambda^{\bar{L}}_f\left[ \sum_y \left(-\log(\eta_{y,f}) - \dfrac{[\bar{L}^{\textrm{obs}}_{y,f} - \bar{L}^{\textrm{pred}}_{y,f}]^2}{2 \eta^2_{y,f}}\right)\right],\] for mean lengths, or

\[\Lambda_5 = \sum_f \lambda^{\bar{w}}_f\left[ \sum_y \left(-\log(\eta_{y,f}) - \dfrac{[\bar{w}^{\textrm{obs}}_{y,f} - \bar{w}^{\textrm{pred}}_{y,f}]^2}{2 \eta^2_{y,f}}\right)\right],\] for mean weights, where \(\eta_{y,f}\) is the standard deviation of the mean size. With constant coefficient of variation (CV), \(\eta_{y,f} = \bar{w}^{\textrm{obs}}_{y,f} CV^{\eta}_f\).

The log-likelihood component \(\Lambda_6\) of annual estimated recruitment deviates \(\delta_y\) in log space is \[\Lambda_6 = \Sigma_y\left(-\log(\tau) - \dfrac{\delta_y^2}{2 \tau^2}\right),\] where \(\tau\) is the standard deviation of recruitment deviates.

The log-likelihood component \(\Lambda_7\) of the equilibrium catch is \[\Lambda_7 = \sum_f \lambda^{Y}_f \left(-\log(\omega^{\textrm{eq}}_f) - \dfrac{[\log(Y^{\textrm{eq,obs}}_f) - \log(Y^{\textrm{eq,pred}}_f)]^2}{2 \times (\omega^{\textrm{eq}}_f)^2}\right),\]

The log-likelihood component \(\Lambda_8\) of the survey proportion-at-age data is \[\Lambda_8 = \sum_s \lambda^{IA}_s \left[\sum_y O^{IA}_{y,s} \sum_a p^{\textrm{obs}}_{y,a,s} \log(p^{\textrm{pred}}_{y,a,s})\right],\] where \(O^{IA}\) is the annual sample sizes for the survey age compositions.

The log-likelihood component \(\Lambda_9\) of the survey proportion-at-length data is \[\Lambda_9 = \sum_s \lambda^{IL}_s \left[ \sum_y O^{IL}_{y,s} \sum_{\ell} p^{\textrm{obs}}_{y,\ell,s} \log(p^{\textrm{pred}}_{y,\ell,s})\right]\] where \(O^{IL}\) is the annual sample sizes for the survey length compositions.

11.3 Priors

Vague priors are added to selectivity parameters to aid in convergence, with \[ \begin{align} x^{\textrm{LFS}}_b &\sim N(0,3)\\ x^{\textrm{L5}}_b &\sim N(0,3)\\ V^{L_{\infty}}_b &\sim \textrm{Beta}(1.01, 1.01) \end{align}\]

If free selectivity parameters are estimated, then \[v_{y,a,f} \sim \textrm{Beta}(1.01, 1.01)\]

11.4 Optional priors

A lognormal prior on \(R_0\) can be specified. The penalty term \(\textrm{pr}\) to the objective function is: \[\textrm{pr}^{R_0} = -\log(\sigma^{R_0}) - 0.5\left[\dfrac{\log(R_0/\mu^{R_0})}{\sigma^{R_0}}\right]^2,\] where \(\mu^{R_0}\) and \(\sigma^{R_0}\) are the mean and standard deviations for the prior provided by the user.

Recall that the parameter corresponding to Beverton-Holt steepness \(x^h\) is estimable over all real numbers and is transformed such that \(h = 0.8 \times \textrm{logit}^{-1}(x^h) + 0.2\). The prior uses a beta distribution on \(y = (h - 0.2)/0.8\). The penalty to the objective function is: \[\textrm{pr}^{h} = \log(y) (\alpha-1) + \log(1-y) * (\beta-1) - \log(y - y^2),\] where \(\alpha = \mu^h \left(\frac{\mu^h (1 - \mu^h)}{(\sigma^h)^2} - 1\right)\), \(\beta = (1-\mu^h)\left(\frac{\mu^h (1 - \mu^h)}{(\sigma^h)^2} - 1\right)\), and the last term is the log Jacobian transform of the logit function.

With a Ricker SR function, a normal distribution is used: \[\textrm{pr}^{h} = -\log(\sigma^{h}) - 0.5\left(\dfrac{h - \mu^{h}}{\sigma^{h}}\right)^2.\]

The prior for \(M\) is lognormal: \[\textrm{pr}^{M} = -\log(\sigma^{M}) - 0.5\left[\dfrac{\log(M/\mu^{M})}{\sigma^{M}}\right]^2.\]

Priors for survey catchability \(q_s\) are normally distributed: \[\textrm{pr}^{q} = \sum_s\left(-\log(\sigma^{q}_s) - 0.5\left(\dfrac{q_s - \mu^{q}_s}{\sigma^{q}_s}\right)^2\right).\]

11.5 Penalties

A penalty to the likelihood is added when \(F_{y,f} \ge F_{\textrm{max}}\) for any year and fleet. The penalty is

\[\textrm{pen}_{y,f} = \begin{cases} 0.01 (F_{y,f} - F_{\textrm{max}})^2 & F_{y,f} \ge F_{\textrm{max}}\\ 0 & \textrm{otherwise} \end{cases}.\]

12 Additional resources

The help file for RCM will provide more information on the possible inputs for the model. The help file for the RCModel class (obtained by typing class?RCModel into the R console) describes the outputs from the function. An additional vignette is available to describe set up of fleet and survey selectivity in the function call.