Delay Difference (DD_TMB) Model
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
Growth in weight-at-age \(w_a\) follows the recursive Ford-Brody equation: \[w_a = \tilde\alpha + \rho 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 \(a \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. The age corresponding to the length of 50% maturity is chosen for \(k\).
Dynamics equations
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 later) at age \(k\) and \(w_k\) is the weight of recruits. Survival \(s_t\) is defined as \[ s_t = \exp(-Z_t), \] where the total mortality rate \(Z_t = F_t + M\) is the sum of fishing mortality \(F_t\) and natural mortality \(M\).
Conditioning on catch
If the model is conditioned on catch (\(C_t\)), then \(F_t\) is solved iteratively such that \[ C_t = \dfrac{F_t}{Z_t} [1 - \exp(-Z_t)]B_t. \]
The predicted biomass-based index for survey \(s\) is \[\hat{I}_{s,t} = \hat{q}^I_s \hat{B}_t ,\] where \(q^I_s\) is the catchability coefficient scaling the stock size to the index and the circumflex \(^\) denotes the model estimate. Otherwise, the predicted abundance-based index, is \[\hat{I}_{s,t} = \hat{q}^I_s \hat{N}_t .\]
The likelihood \(L^I\) of the observed indices \(I_{s,t}\), using a lognormal distribution, is \[L^I = \sum_s \lambda_s \sum_t \left(-\log(\sigma_{s,t}) - 0.5\left[\dfrac{\log(I_{s,t}) - \log(\hat{I}_{s,t})}{\sigma_{s,t}}\right]^2\right),\] where \(\sigma_{s,t}\) is the standard deviation of the index.
Conditioning on effort
Otherwise, if the model is conditioned on effort (as the ratio of the catch and index), then fishing mortality is set proportional to effort \[ F_t = q^f f_t, \] where \(q^f\) is the coefficient that scales effort and \(f_t\) is the effort in year \(t\).
The predicted catch is \[\hat{C}_t = \frac{\hat{F}_t}{\hat{Z}_t} [1 - \exp(-\hat{Z}_t)] \hat{B}_t,\] where the circumflex \(^\) denotes the model estimate.
The likelihood \(L^C\) of the observed catch \(C_t\), following a lognormal distribution, is \[L^C = \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.
Stock-recruit parameters
Beverton-Holt relationship
With a Beverton-Holt stock recruit relationship, and spawning occurring before fishing in each annual time step, the 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\) is an estimated parameter and steepness \(h\) can either be fixed or estimated. Unfished biomass \(B_0\) is 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))\).
Ricker equation
With 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 18.
Continuous Delay-Differential model (cDD)
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
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\).
Dynamics equations
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 29 and 30 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\), and \(\Delta\) is the length of the time step, i.e., \(\Delta = 1\) for an annual time step.
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 using Newton’s method.
The predicted index for survey \(s\) is \[\hat{I}_{s,t}=\hat{q}_s \hat{B}_{s,t}\] or \[\hat{I}_{s,t}=\hat{q}_s \hat{N}_{s,t}\] for a biomass and abundance-based index, respectively.
The likelihood of the observed index \(I_t\), using a lognormal distribution, is \[L^I = \sum_s \lambda_s \sum_t \left(-\log(\sigma_{s,t}) - 0.5\left[\dfrac{\log(I_{s,t}) - \log(\hat{I}_{s,t})}{\sigma_{s,t}}\right]^2\right),\] where \(\lambda_s\) is a weighting coefficient and \(\sigma_{s,t}\) is the standard deviation of the index.
Stock-recruit parameters
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 &= R_0 \phi_0 = R_0 \dfrac{w_k + \dfrac{\kappa W_{\infty}}{M}}{M + \kappa}\\ N_0 &= \dfrac{R_0}{M}. \end{align}\]
State-space version (DD_SS and cDD_SS)
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).\]
Mean weight
If either model is fitted to mean weight data \(\bar{w}_t\), then the corresponding likelihood is
\[L^W = \sum_t \left(-\log(\sigma^W) - 0.5\left[\dfrac{\log(\bar{w}_t) - \log(\hat{B_t}/\hat{N_t})}{\sigma^W}\right]^2\right)\] where \(\sigma_W\) is the mean weight standard deviation and \(\hat{B_t}/\hat{N_t}\) is the model predicted mean weight.
References
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.