In MSEtool, assessment models are of class Assess. This appendix provides a brief description and references for the Assess objects. Further details regarding parameterization, e.g., fixing parameters, and tuning, e.g., adjusting start parameters, are provided in the function documentation.
For LaTeX equation rendering, it is recommended that this vignette be viewed in a HTML browser. This can be done with the browseVignettes function in R:
The statistical catch-at-age model uses a time series of total catch (in weight), index, and catch-at-age observations, as well as information on weight, maturity, natural mortality at age.
The population abundance \(N_{a,t}\) of age \(a\) in year \(t\) is \[ N_{a,t} = \begin{cases} R_t & a = 1\\ N_{a-1,t-1} \exp(- v_{a-1} F_{t-1} - M_{a-1}) & a = 2, \ldots, A-1\\ N_{a-1,t-1} \exp(- v_{a-1} F_{t-1} - M_{a-1}) + N_{a,t-1} \exp(- v_a F_{t-1} - M_a) & a = A \end{cases}, \] where \(R_t\) is the recruitment (age-1), \(v_a\) is the vulnerability at age \(a\), \(F_t\) is the apical fishing mortality rate, \(M_a\) is the instantaneous natural mortality rate at age \(a\), and \(A\) is the maximum age in the model as a plus-group accumulator age.
Assuming logistic vulnerability, the vulnerability is: \[v_a = \left[1 + \exp\left(-\log(19) \dfrac{a - a_{50}}{a_{95} - a_{50}}\right)\right]^{-1}, \] where \(a_{50}\) and \(a_{95}\) are the estimated ages of 50% and 95% vulnerability, respectively.
Assuming dome vulnerability, a double Gaussian formulation is used: \[ v_a = \begin{cases} f(a; a_{asc}, \sigma_{asc}) & a \le a_{asc}\\ 1 & a_{asc} \lt a \le a_{des}\\ f(a; a_{des}, \sigma_{des}) & a \gt a_{des} \end{cases}, \] where \(f(a; \mu, \sigma) = \exp(-0.5 (a - \mu)^2/\sigma^2)\) is the normal probability density function scaled to one at \(\mu\). Four parameters are estimated: \(a_{50}\) the age of 50% vulnerability (ascending limb), \(a_{asc}\) the first age of full vulnerability, \(a_{des}\) the last age of full vulnerability, and \(v_A\) the vulnerability at the maximum age in the model. The \(\mu\) and \(\sigma\) for both the ascending and descending limbs of the double-normal equation are estimated parameters. From these four parameters, \(\sigma_{asc} = \sqrt{(a_{50} - \mu_{asc})^2/\log(4)}\) and \(\sigma_{des} = \sqrt{-0.5(A - \mu_{des})^2/\log(v_A)}\) can be derived.
The vulnerable biomass \(VB_t\) at the beginning of year \(t\) is \[VB_t = \sum_{a=1}^A v_a w_a N_{a,t},\] where weight-at-age \(w_a\) is given by \[w_a = W_{\infty}(1 - \exp[K\{a-a_0\}])^b.\]
The mature spawning biomass \(E_t\) is given by \[E_t = \sum_{a=1}^A m_a w_a N_{a,t},\] where maturity at age \(m_a\) is \[m_a = \left[1 + \exp\left(-\log(19) \dfrac{a - \tilde{a}_{50}}{\tilde{a}_{95} - \tilde{a}_{50}}\right)\right]^{-1}, \] where \(\tilde{a}_{50}\) and \(\tilde{a}_{95}\) are the ages of 50% and 95% maturity, respectively.
The estimated catch-at-age \(\hat{C}_{a,t}\) is obtained from the Baranov equation: \[\hat{C}_{a,t} = \dfrac{\hat{v}_a \hat{F}_t}{\hat{v}_a \hat{F}_t + M_a} [1 - \exp(- \hat{v}_a \hat{F}_t - M_a)] \hat{N}_{a,t}.\] The predicted total catch in weight \(\hat{Y}_t\) is \(\hat{Y}_t = \sum_a w_a \hat{C}_{a,t}.\)
The estimated index \(\hat{I}_t\), assuming that it is an index for total biomass, is \[ \hat{I}_t = \hat{q} \hat{B}_t,\] where \(B_t = \sum_{a=1}^A w_a N_{a,t}\). A function argument allows for the user to specify that the observed index is for vulnerable or spawning biomass.
The likelihood of the observed catch at age \(C_{a,t}\), assuming a multinomial distribution, is \[ C_{a,t} \sim \textrm{Multinomial}(O_t, \hat{p}_{a,t}), \] where \(O_t\) is the assumed sample size of catch-at-age observations in year \(t\) and \(\hat{p}_{a,t} = \hat{C}_{a,t}/\sum_a\hat{C}_{a,t}\) is annual predicted proportions of catch-at-age.
If a lognormal distribution for the observed proportions at age is assumed, then the likelihood is \[\log(p_{a,t}) \sim N(\log[\hat{p}_{a,t}], 0.01/\hat{p}_{a,t}).\]
The likelihood of the observed catch \(Y_t\), assuming a lognormal distribution, is \[\log(Y_t) \sim N(\log[\hat{Y}_t], \omega^2).\]
The likelihood of the observed index \(I_t\), assuming a lognormal distribution, is \[\log(I_t) \sim N(\log[\hat{I}_t], \sigma^2)\]
A variant of the SCA is the SCA_Pope function, which fixes the predicted catches to the observed catches and uses Pope’s approximation to calculate the annual harvest rate in the midpoint of the year.
The abundance at age is
\[ N_{a,t} = \begin{cases} R_t & a = 1\\ N_{a-1,t-1} (1 - v_{a-1} u_{t-1}) \exp(-M_{a-1}) & a = 2, \ldots, A-1\\ N_{a-1,t-1} (1 - v_{a-1} u_{t-1}) \exp(-M_{a-1}) + N_{a,t-1} (1 - v_a u_{t-1}) \exp(-M_a) & a = A \end{cases}, \] where \(u_t\) is the harvest rate.
The vulnerable biomass in the midpoint of the year is \[VB^{mid}_t = \sum_{a=1}^A v_a w_a N_{a,t} \exp(-0.5 M_a).\]
By conditioning the model on catch in weight \(Y_t\), the estimated annual harvest rate \(\hat{u}_t\) is \[\hat{u}_t = Y_t / \widehat{VB^{mid}}_t .\] The estimated catch at age \(\hat{C}_{a,t}\) is \[\hat{C}_{a,t} = \hat{v}_a \hat{u}_t \hat{N}_{a,t} \exp(-0.5 M_a).\]
There are two variants of the statistical catch-at-age model for estimation of recruitment and reference points. In function SCA, productivity parameters \(R_0\) and \(h\) are estimated in the assessment model. Annual recruitment is estimated as deviations from the resulting stock-recruitment relationship within the model. MSY reference points are derived from the estimates of \(R_0\) and \(h\).
In SCA2, no stock-recruit relationship is assumed in the assessment model, i.e., annual recruitment is estimated as deviations from the mean recruitment over the observed time series, similar to Cadigan (2016). After the assessment, a stock-recruit function can be fitted post-hoc to the recruitment and spawning stock biomass estimates from the assessment model to obtain MSY reference points.
Recruitment \(R_t\) in year \(t\) is \[ R_t = \dfrac{\alpha E_{t-1}}{1 + \beta E_{t-1}} \exp(\delta_t - 0.5 \tau^2),\] where \(\delta_t \sim N(0, \tau^2)\) are recruitment deviations from the stock-recruit relationship in lognormal space and \(\tau\) is the standard deviation of the recruitment deviations. Parameters \(\alpha\) and \(\beta\) are defined as \[ \begin{align} \alpha &= \dfrac{4hR_0}{(1-h)B_0}\\ \beta &= \dfrac{5h-1}{(1-h)B_0}, \end{align}\] where \(B_0 = R_0 \phi_0\). The biomass per recruit \(\phi_0\) is calculated as \(\phi_0 = \sum_{a=1}^A m_a w_a l_a\), where \[ l_a = \begin{cases} 1 & a = 1\\ l_{a-1} \exp(-M_{a-1}) & a = 2, \ldots, A-1\\ \dfrac{l_{a-1} \exp(-M_{a-1})}{1 - \exp(-M_a)} & a = A \\ \end{cases}. \]
Recruitment \(R_t\) in year \(t\) is \[ R_t = \alpha E_{t-1} \exp(-\beta E_{t-1}) \exp(\delta_t - 0.5 \tau^2),\] where \[ \begin{align} \alpha &= \dfrac{(5h)^{1.25} R_0}{B_0}\\ \beta &= \dfrac{5}{4B_0}\log(5h). \end{align}\]
Recruitment \(R_t\) in year \(t\) is \[R_t = \bar{R} \exp(\delta_t - 0.5 \tau^2), \] where \(\delta_t \sim N(0, \tau^2)\) are recruitment deviations in lognormal space from the estimated mean recruitment \(\bar{R}\) and \(\tau\) is the standard deviation of the recruitment deviations. Typically, \(\tau\) is set to 1 so that recruitment is estimated almost as free parameters (Cadigan, 2016).
The likelihood of the estimated log-recruitment deviations \(\hat{\delta}_t\) is \[\hat{\delta}_t \sim N(0, \tau^2).\]
Cadigan, N.G. 2016. A state-space stock assessment model for northern cod, including under-reported catches and variable natural mortality rates. Canadian Journal of Fisheries and Aquatic Science 72:296-308.