In SAMtool, 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:
browseVignettes("SAMtool")There has been a rich history of development for the delay difference model for catch and index data. For the formulation used in SAMtool, the most relevant citations are Chapter 9 of Hilborn and Walters (1992) and Carruthers et al. (2012).
Growth in weight-at-age \(w_a\) follows the recursive Ford-Brody equation: \[w_a = \rho + \tilde\alpha w_{a-1}.\] We can obtain \(\tilde\alpha\) and \(\rho\) for the delay difference model if weight is also described by the equation \[w_a = W_{\infty}(1 - \exp[K\{a-a_0\}])^b.\] Parameter \(\tilde\alpha\) is solved in the limiting case where \(w_a = w_{a-1} = W_{\infty}\) as \(t \rightarrow \infty\), \[\tilde\alpha = W_{\infty}(1 - \rho). \] Substitution of equation 3 into equation 1 solves for \(\rho\), \[\rho = \dfrac{w_a - W_{\infty}}{w_{a-1} - W_{\infty}}.\] In SAMtool, \(a = k+2\) is arbitrarily chosen to calculate \(\rho\), where \(k\) is the age of knife-edge selectivity. From catch and effort data alone, the age corresponding to the length of 50% maturity is chosen for \(k\).
The population biomass \(B_t\) and abundance \(N_t\) in year \(t\) is given by \[ \begin{align} B_t &= s_{t-1}(\tilde{\alpha} N_{t-1} + \rho B_{t-1}) + w_k R_t\\ N_t &= s_{t-1} N_{t-1} + R_t, \end{align}\] where \(R_t\) is the recruitment (defined in the next section) at age \(k\) and \(w_k\) is the weight of recruits. Survival \(s_t\) is defined as \[ s_t = (1 - u_t)\exp(-M), \] where \(u_t\) is the annual harvest rate and \(M\) is the instantaneous natural mortality rate.
If the model is conditioned on catch (\(C_t\)), the harvest rate \(u_t\) is \[ u_t = C_t/B_t. \]
The predicted biomass-based index is \[\hat{I}_t = \hat{q}^I \hat{B}_t ,\] where the circumflex \(^\) denotes the model estimate. Otherwise, the predicted abundance-based index, is \[\hat{I}_t = \hat{q}^I \hat{N}_t .\]
The likelihood \(L\) of the observed index \(I_t\), assuming a lognormal distribution, is \[L^I = \sum_t \left(-\log(\sigma) - 0.5\left[\dfrac{\log(I_t) - \log(\hat{I}_t)}{\sigma}\right]^2\right),\] where \(\sigma\) is the standard deviation of the index.
Otherwise, if the model is conditioned on effort (as the ratio of the catch and index), the harvest rate is \[ u_t = 1 - \exp(-q^f f_t), \] where \(q\) is the estimated coefficient that scales effort and \(f_t\) is the effort in year \(t\).
The predicted catch is \[\hat{C}_t = \hat{u}_t \hat{B}_t,\] where the circumflex \(^\) denotes the model estimate.
The likelihood \(L\) of the observed catch \(C_t\), assuming a lognormal distribution, is \[L^I = \sum_t \left(-\log(\omega) - 0.5\left[\dfrac{\log(C_t) - log(\hat{C}_t)}{\omega}\right]^2\right).\] where \(\omega\) is the standard deviation of the catch.
Assuming a Beverton-Holt stock recruit relationship and spawning occurring after fishing in each annual time step, then recruitment (at age \(k\)) in year \(t\) is: \[ R_t = \dfrac{\alpha B_{t-k}}{1 + \beta B_{t-k}},\] where \[ \begin{align} \alpha &= \dfrac{4hR_0}{(1-h)B_0}\\ \beta &= \dfrac{5h-1}{(1-h)B_0}, \end{align}\]
Unfished recruitment \(R_0\) and steepness \(h\) are estimated parameters, with unfished biomass \(B_0\) calculated as \[B_0 = R_0 \phi_0.\] The unfished biomass per recruit \(\phi_0\) is \[\phi_0 = \dfrac{\tilde{\alpha} \exp(-M) + w_k (1 - \exp(-M))}{1 - \rho \exp(-M)}\] and is obtained by solving the equilibrium equation for biomass, \(B_0 = \exp(-M)(\tilde{\alpha}N_0 + \rho B_0) + w_k R_0\), is solved for \(B_0/R_0\), with \(N_0 = R_0/(1−\exp(-M))\).
Assuming a Ricker stock-recruit relationship, the recruitment is \[ R_t = \alpha B_{t-k}\exp[-\beta B_{t-k}],\] where \[ \begin{align} \alpha &= \dfrac{(5h)^{1.25} R_0}{B_0}\\ \beta &= \dfrac{5}{4B_0}\log(5h), \end{align}\]
and \(B_0\) is calculated as in equation 14.
Compared to the discrete delay-difference (annual time-step in production and fishing), the delay-differential model (cDD) is based on continuous recruitment and fishing mortality within a time-step. The continuous model works much better for populations with high turnover (e.g. high F or M, continuous reproduction).
Growth in weight is modeled as a von Bertalanffy equation: \[ \dfrac{dw_{a,t}}{da} = \kappa (W_{\infty} - w_{a,t}).\] A solution to Equation 19 is \[w_{a+1,t} = W_{\infty} + (w_{a,t} - W_{\infty})\exp(-\kappa).\] From a mean weight-at-age schedule for ages \(a = k, k+1, \ldots\), a non-linear regression can be used to obtain \(\kappa\).
The governing equations for the pooled biomass \(B_t\) and abundance \(N_t\) over time \(t\) can be written as \[ \begin{align} \dfrac{dN_t}{dt} &= \dfrac{d}{dt} \int N_{a,t}da\\ \dfrac{dB_t}{dt} &= \dfrac{d}{dt} \int w_{a,t}N_{a,t}da, \end{align}\] where the integration is over ages \(k\) to \(\infty\).
To evaluate the integral, we make substitutions based on the following: \[ \begin{align} \dfrac{dN_{a,t}}{dt} &= \dfrac{dN_{a,t}}{da}\dfrac{da}{dt} = -Z_t N_{a,t}\\ \dfrac{dw_{a,t}}{dt} &= \dfrac{dw_{a,t}}{da}\dfrac{da}{dt} = \kappa (W_{\infty} - w_{a,t}). \end{align}\] After substitution and evaluation of the integrals, the governing equations are \[ \begin{align} \dfrac{dN_t}{dt} &= R_t - Z_t N_t\\ \dfrac{dB_t}{dt} &= w_k R_t + \kappa W_{\infty} N_t - (Z_t + \kappa) B_t, \end{align}\] where \(R_t\) is the abundance of recruits and \(w_k R_t\) is the weight of recruits, both of which serve as the constants of integration describing the input of abundance and biomass, respectively, to the population.
Solving the differential equations 25 and 26 leads to the dynamics equations: \[ \begin{align} N_{t+\delta} &= N_{\infty,t} + (N_t - N_{\infty,t})\exp(-Z_t\delta)\\ B_{t+\delta} &= B_{\infty,t} + (N_t - N_{\infty,t}) \dfrac{\kappa W_{\infty}}{Z_t + \kappa} + \left[B_t - B_{\infty,t} - (N_t - N_{\infty,t}) \dfrac{\kappa W_{\infty}}{Z_t + \kappa}\right]\exp(-[Z_t+\kappa]\delta), \end{align}\] where \(Z_t = F_t + M\), and \(N_{\infty,t} = R_t/Z_t\) and \(B_{\infty,t} = \dfrac{w_k + \dfrac{\kappa W_{\infty}}{Z_t}}{Z_t + \kappa}R_t\) are the equilibrium abundance and biomass respectively, conditional on \(R_t\) and \(Z_t\).
With a constant and continuous fishing mortality rate \(F_t\) over time step \(s = t\) to \(s = t + \delta\), the accumulated catch in weight \(C_t\) is \[\begin{align} C_t &= \int F_t B_s ds\\ &= F_t\left[B_{\infty,t}\delta + (N_t - N_{\infty,t}) \dfrac{\kappa W_{\infty}}{Z_t + \kappa} \delta + \dfrac{B_t - B_{\infty,t} - (N_t - N_{\infty,t}) \dfrac{\kappa W_{\infty}}{Z_t + \kappa}[1 - \exp(-[Z_t+\kappa]\delta)]}{Z_t + \kappa} \right] \end{align}\] To match the predicted catch to the observed catch in the model, \(F_t\) is solved iteratively.
The predicted index is \[\hat{I}_t=\hat{q}\hat{B}_t\] or \[\hat{I}_t=\hat{q}\hat{N}_t\] for a biomass and abundance-based index, respectively.
The likelihood of the observed index \(I_t\), assuming a lognormal distribution, is \[L^I = \sum_t \left(-\log(\sigma) - 0.5\left[\dfrac{\log(I_t) - \log(\hat{I}_t)}{\sigma}\right]^2\right),\] where \(\sigma\) is the standard deviation of the index.
The stock-recruit parameters are estimated in the same way as the delay difference model, except unfished abundance \(N_0\) and biomass \(B_0\) are calculated as \[\begin{align} B_0 &= \dfrac{w_k + \dfrac{\kappa W_{\infty}}{M}}{M + \kappa}R_0\\ N_0 &= \dfrac{R_0}{M}. \end{align}\]
In the state-space version, annual recruitment deviates from the stock-recruit relationship are estimated. The recruitment in year \(t\) is \[ R_t = \dfrac{\alpha B_{t-k}}{1 + \beta B_{t-k}} \exp(\delta_t - 0.5 \tau^2)\] or \[ R_t = \alpha B_{t-k}\exp(-\beta B_{t-k})\exp(\delta_t - 0.5 \tau^2),\] where \(\delta_t\) are recruitment deviations in lognormal space and \(\tau\) is the standard deviation of the recruitment deviations.
Log-recruitment deviations \(\hat{\delta}_t\) are typically estimated as penalized parameters towards the likelihood, with the penalty: \[L^{\delta} = \sum_t \left(-\log(\tau) - 0.5\left[\dfrac{\hat{\delta}_t}{\tau}\right]^2\right).\]
Carruthers, T., Walters, C.J., and McAllister, M.K. 2012. Evaluating methods that classify fisheries stock status using only fisheries catch data. Fisheries Research 119-120:66-79.
Hilborn, R., and Walters, C. 1992. Quantitative Fisheries Stock Assessment: Choice, Dynamics and Uncertainty. Chapman and Hall, New York.