In tcarruth/MSEtool: Management Strategy Evaluation Toolkit

knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
knitr::opts_chunk$set(dpi=85)
options(width = 650)

Introduction

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:

browseVignettes("MSEtool")

Virtual Population Analysis (VPA)

The VPA model reconstructs the historical population based on an expanded catch at age matrix and an index of abundance. From an estimate of the fishing mortality in the terminal year, historical abundance and fishing mortality are back-calculated. Natural mortality is a required input parameter. For a biomass-based index, weight at age is also required. Maturity at age is also required if spawning biomass is to be calculated. The dynamics equations are based on VPA-2BOX (Porch 2018), although 2-stock mixing and diffusion processes are not modeled here. The maximum age in the model is a plus-group accumulator age.

Dynamics equations

In this model, fishing mortality $F_{a,t}$ of age $a$ in year $t$ is generally solved from the observed catch of the corresponding age and year $C_{a,t}$ and the estimated abundance of the same cohort in the following year $N_{a+1,t+1}$: $$F_{a,t} = \dfrac{Z_{a,t}C_{a,t}}{[\exp(Z_{a,t}) - 1]N_{a+1,t+1}},$$ where $Z_{a,t} = F_{a,t} + M_a$, and $M$ is natural mortality. Equation 1 is Baranov's equation with the substitution: $N_{a,t}=N_{a+1,t+1}\exp(Z_{a,t})$. Newton's method is used to numerically solve Equation 1.

With $F_{a,t}$, the abundance of the corresponding age and year is solved also from Baranov's equation: $$N_{a,t} = \dfrac{C_{a,t}Z_{a,t}}{F_{a,t}[1 - \exp(-Z_{a,t})]}.$$

There are two exceptions to this algorithm. First, for the terminal year $t=T$, F is estimated in the model, followed by the calculations for the terminal year abundances with Equation 2. This starts the back-calculations for $F$ and $N$ for the preceding years with equations 1 and 2, respectively. The terminal F-at-age vector $F_{a,T}$ can be free parameters with age or constrained as a logistic or double-normal (dome) function with age.

Second, the fishing mortality rates of the plus-group age $a=A$ and the preceding age $a=A-1$ are subject to additional constraints and are solved separately. First, note that the abundance of the plus-group in year $t+1$ must satisfy the equation $$ N_{A,t+1} = N_{A,t}\exp(-Z_{A,t}) + N_{A-1,t}\exp(-Z_{A-1,t}),$$ with $$\begin{align} N_{A,t} &= \dfrac{C_{A,t}Z_{A,t}}{F_{A,t}[1 - \exp(-Z_{A,t})]}\ N_{A-1,t} &= \dfrac{C_{A-1,t}Z_{A-1,t}}{F_{A-1,t}[1 - \exp(-Z_{A-1,t})]} \end{align}.$$

An additional constraint is needed for the plus-group $F$, where $$F_{A,t} = \alpha F_{A-1,t},$$ where $\alpha$ is a numeric constant either fixed (often to 1) or estimated in the model.

From $N_{A,t+1}$, $C_{A,t}$, $C_{A-1,t}$, $M_A$, and $M_{A-1}$, we can solve $F_{A-1,t}$ and $F_{A,t}$ as the values that satisfy equations 3-6. The corresponding abundances $N_{A-1,t}$ and $N_{A,t}$ are then obtained with equation 2.

The VPA is tuned to an index via maximum likelihood. The predicted biomass index is $$ \hat{I}t = \hat{q}\sum_a N{a,t} w_a.$$ The likelihood of the observed index is $$\log(I_t) \sim N(\log(\hat{I}_t), \sigma^2).$$

From the VPA output, reference points can be calculated. Estimates of recruitment (the youngest age class $\tilde{a}$) and spawning biomass are used to estimate a stock-recruitment relationship to obtain MSY and unfished reference points.

Additional constraints

In any VPA model, estimates of F and N of the youngest age classes are highly uncertain. To stabilize estimates, random walk penalties can be applied to the recruitment and vulnerabilty estimates in the most recent years: $$ \begin{align} \log(N_{\tilde{a},t}) &\sim N(\log(N_{\tilde{a}, t-1}), \sigma_{\tilde{a}}^2)\ \log(v_{a,t}) &\sim N(\log(v_{a,t-1}, \sigma_v^2) \end{align},$$ where $\tilde{a}$ is the smallest age in the VPA and $v_{a,t} = F_{a,t}/\max_a F_{a,t}$.

Note that these penalties can significantly alter model results, and generate spurious trends in F and N especially if there are significantly shifts in recent vulnerability and recruitment. The implications of incorporating these random walk penalties need to be evaluated on a case-by-case basis.