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 surplus production model uses the Fletcher (1978) formulation. The biomass \(B_t\) in year \(t\) is \[B_t = B_{t-1} + P_{t-1} - C_{t-1},\] where \(C_t\) is the observed catch and \(P_t\) is the surplus production given by: \[P_t = \gamma \times MSY \times \left(\dfrac{B_t}{K}-\left[\dfrac{B_t}{K}\right]^n\right), \] where \(K\) is the carrying capacity, \(MSY\) is the estimated maximum sustainable yield, and \(n\) is the parameter that controls shape of the production curve, and \(\gamma\) is \[\gamma = \dfrac{1}{n-1}n^{n/(n-1)}.\]
By conditioning the model on observed catch, the predicted index \(\hat{I}_t\) is \[\hat{I}_t = \hat{q} \hat{B}_t \] and the harvest rate is \[\hat{F}_t = \dfrac{C_t}{\hat{B}_t}.\] The dynamics equations above use an annual time step. Optionally, smaller time steps are used in the model to approximate continuous production and fishing. Given the biomass in the start of the year and assuming a constant fishing mortality over the time steps within a year, the fishing mortality that produces the observed annual catch is solved iteratively.
The likelihood of the observed index \(I_t\), assuming a lognormal distribution, is \[\log(I_t) \sim N(\log[\hat{I}_t], \sigma^2).\]
From estimates of leading parameters \(F_{MSY}\) and \(MSY\), the biomass \(B_{MSY}\) at \(MSY\) is \[B_{MSY} = \dfrac{MSY}{F_{MSY}},\] the carrying capacity \(K\) is \[K = n^{1/(n-1)} \times B_{MSY} ,\] and the intrinsic rate of population increase \(r\) is \[ r = n \times F_{MSY}.\] The production parameter \(n\) is typically fixed and the model has a symmetric productive curve (\(B_{MSY}/K = 0.5\)) when \(n = 2\).
The production parameter \(n\) is typically fixed and the Schaefer model with a symmetric productive curve (\(B_{MSY}/K = 0.5\)) is produced when \(n = 2\).
The Fox model is the limiting case of the Fletcher parameterization as \(n \rightarrow 1\), where
\[K = e \times B_{MSY}\] \[r = F_{MSY}\] \[P_t = -e \times MSY \times \dfrac{B_t}{K} \times \log\left(\dfrac{B_t}{K}\right) \]
In the state-state version, annual biomass deviates are estimated as random effects. Similar to Meyer and Millar (1999), the biomass \(B_t\) in year \(t\) is \[B_t = (B_{t-1} + P_{t-1} - C_{t-1})\exp(\delta_t - 0.5 \tau^2),\] where \(\delta_t \sim N(0, \tau^2)\) are biomass deviations in lognormal space and \(\tau\) is the standard deviation of the biomass deviations.
The log-likelihood of the estimated deviations \(\hat{\delta}_t\) is \[\hat{\delta}_t \sim N(0, \tau^2).\]
To generate the prior for the intrinsic rate of increase, natural mortality \(M_a\) and steepness \(h\) are sampled from a distribution. Natural mortality is modelled to be age-invariant and is sampled from a lognormal distribution. Assuming either a Beverton-Holt or Ricker stock-recruit relationship, steepness is sampled from a transformed beta or transformed lognormal distribution, respectively.
For each pair of sampled M and h values, the corresponding value of \(r\) is obtained by solving a modified Euler-Lotka equation: \[\Sigma_{a=1}^A l_a m_a \exp(-r \times a) = 1.\]
Equation 15 is modified to include the \(\alpha\) term from the stock-recruit relationship (Stanley et al. 2009). In this way, the recruits-per-spawner at low stock sizes, i.e., as spawning biomass approaches zero, is considered for calculating \(r\).
The numbers-per-recruit at age \(a\) is \[ 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}. \]
Fecundity at age \(m_a\) is \[m_a = \dfrac{\alpha w_a}{\left[1 + \exp\left(-\log(19) \dfrac{a - \tilde{a}_{50}}{\tilde{a}_{95} - \tilde{a}_{50}}\right)\right]}, \] where \(\tilde{a}_{50}\) and \(\tilde{a}_{95}\) are the ages of 50% and 95% maturity, respectively.
Weight-at-age \(w_a\) is \[w_a = W_{\infty}(1 - \exp[K\{a-a_0\}])^b.\]
The recruits per spawner at the origin of the stock-recruit relationship \(\alpha\) is \[\alpha = \dfrac{4h}{(1-h)\phi_0}, \] or \[\alpha = \dfrac{(5h)^{1.25}}{\phi_0},\] for a Beverton-Holt and Ricker stock-recruit relationship, respectively, where unfished recruits-per-spawner \(\phi_0\) is \[\phi_0 = \Sigma_{a=1}^A l_a w_a m_a.\]
A normal distribution is assumed for the prior with the mean and standard deviation calculated from the values of \(r\) calculated using the above procedure.
Fletcher, R.I. 1978. On the restructuring of the Pella-Tomlinson system. Fishery Bulletin 76:515-521.
Meyer, R., and Millar, R.B. 1999. BUGS in Bayesian stock assessments. Canadian Journal of Fisheries and Aquatic Science 56:1078-1086.
Stanley, R.D., M. McAllister, P. Starr and N. Olsen. 2009. Stock assessment for bocaccio (Sebastes paucispinis) in British Columbia waters. DFO Can. Sci. Advis. Sec. Res. Doc. 2009/055. xiv + 200 p.