Chapter 4 Spatial stock assessment model for Alaskan sablefish

The spatial model used for this research is available as an R package and is best described by the R package documentation found here. I have described the general model here for completness, but the package documentation is the best resource for specific equations.

Process equations

Population dynamics

The order of processes in an annual cycle follow

Recruitment and release of tagged fish (we apply initial tag induced mortality at this point)
Total mortality and ageing
Markovian movement
Annual tag shedding (applied as a mortality process)

Before applying movement, the partition is updated following \[\begin{align*} N_{a,r,y,s} = \begin{cases} R_{r,y} 0.5, & a = a_1\\ N_{a - 1,r,y - 1,s} \exp\bigg( -Z_{a - 1,r, y - 1,s} \bigg), & a_1 < a < a_+\\ N_{a - 1,r,y - 1,s} \exp\bigg( -Z_{a - 1,r, y - 1,s} \bigg) + N_{a,r,y - 1,s} \exp\bigg( -Z_{a,r, y - 1,s} \bigg), & a = a_+ \end{cases} \end{align*}\] where, \(N_{a,r,y,s}\) is the numbers at age \(a\) in region \(r\), year \(y\) for sex \(s\), \(Z_{a,r,y,s} = M + \sum_g F_{a,r,y,s}^g\) is total mortality and \(R_{r,y}\) is annual recruitment for region \(r\).

Once ageing and mortality have taken place movement is then applied as

\[\begin{equation*} \boldsymbol{N}'_{a,y,s} = \boldsymbol{N}_{a,y,s} \boldsymbol{M} \ \ \forall \ a \end{equation*}\]

where, \(\boldsymbol{N}'_{a,y,s} = (N'_{a,1,y,s}, \dots, N'_{a,n_r,y,s})\) denotes the numbers for age \(a\) across all regions after movement and \(\boldsymbol{M}\) is an \(n_r \times n_r\) movement matrix, which will move age cohort \(a\) among the regions based on the movement matrix.

Initialisation

An equilibrium age structure is derived following the method described in Chapter 18, but for completeness we will briefly describe it here. The annual cycle is run \(n_a - 1\) times to populate all age cohorts prior to the plus group. Then, iterate the annual cycle one more time and calculate the number of individuals that moved into each regions plus age cohort, denoted by \(c^r_{a+}\). This will be the result of ageing, mortality and movement. The equilibrium plus group for region \(r\) is then calculated as

\[ N_{a+, r} = N_{a+ - 1, r} \frac{1}{1 - c^r_{a+}} \ . \]

After the equilibrium age-structure is calculated, there is an option to estimate age specific deviations to allow the model to start with a non-equilibrium age-structure denoted by \(e^{\epsilon_a}\)

\[ N_{a, r} = N_{a, r} e^{\epsilon_a} \ \ \forall \ r \ a \in [a_2..(a_+ - 1)] . \] To help with estimation there is a penalty on \(\epsilon_a\) that assumes a central tendancy of zero with an estimable variance parameter (\(\sigma_{\epsilon}^2\)). \[ \epsilon_a \sim \mathcal{N}(0, \sigma_{\epsilon}) \]

Future model generalizations should consider making this inital age-deviations regional specific.

Growth

Empirical age-length matrices are supplied for all years by sex, these are the same matricies used in the current single area sex. We did not consider spatially varying growth because when the data was visually inspected, there did not seem to be obvious differences. The other reason was when spatially varying growth is included in the model, we would need to track length in the partition (add an extra dimension) to avoid fish drastically changing length as they move between areas and different growth curves.

Mean weight at age was calculated using an allometric length weight relationship with time and space invariant parameters \(\alpha\) and \(\beta\)

\[ \bar{w}_{a,y} = \alpha \bar{l}_{a,y}^{\beta} \]

Recruitment

There are two options in the model for recruitment, regional mean recruitment with global recruitment deviations and regional mean recruitment with region recruitment deviations.

TODO - Add regional stock recruitment relationships as a model option - Add a global stock recruitment relationship

Fishing mortality

When the hybrid \(F\) method is assumed (this is my default), tagged fish were not included when internally solving the fishing mortality nuisance parameters (Chapter 17). The ratio of tagged to untagged numbers of fish in the partition in any year was assumed to be small enough not to effect \(F\) estimates. However, tagged fish were included when the model calculates predicted catch-at-age/length and catch for a fleet. This decision was made to reduce the computational overhead this would require to implement. When \(F\) parameters are estimated as free parameters, then this is not a problem and \(F\) will be applied to both tagged and untagged fish.

Tag release events

Tags release event denoted by the index \(k\) have an implied region \(r\) and year \(y\) dimension. Each tag release event has known sex and age frequency at release. This is derived using the survey age-length key. This is reasonable given the survey is responsible for releasing most tags. A downside of using the age-length key approach to convert unsexed lengths at release into tag releases by age and sex is a tagged fish that would have an exact age, sex and length will be represented in the model as a fraction of a fish across multiple ages and sexes. For release years where there were no age length keys, we used age-length keys pooled over years 1981, 1985 and 1987 for earlier perios (prior to 1981) and the closest age-length key for later periods. (How can we explore uncertainty in this input? re-run the model with numbers at age and length from bootstrapping age-length keys?)

Tagged fish from release event \(k\) are denoted in the partition as \(T^k_{a,r,y,s}\), and are tracked for \(n_{T}\) years before migrating into a pooled tag group, at which point we loose release-year information but do maintain its region of release. At present, tagged fish are assumed to take on the exact same population processes as the untagged elements of the partition (instantaneous mixing).

Other considerations for tag-releases

Do we want to include the inshore/State (Chatham Strait and Clarence Strait) tag-releases? there are alot of fish tagged in these regions.
Do we care about the amount of tags that leave the stock boundaries? See Figure 7.3

Observation equations

There are four observation types in the model

Relative indices of abundance
Age composition disaggregated by sex for the fixed gear fishery and longline survey
Length composition disaggregated by sex for the trawl and early period of the fixed gear fishery
Tag-recovery observations

Catch at age

Fishery dependent catch at age observations are available for the fixed gear fishery, but are also needed to calculate catch at length observations for the trawl fishery. Catch at age for fishery \(g\) is denoted by \({C}^g_{a,r,y,s}\) and model fitted values are calculated following

\[\begin{equation} {C}^g_{a,r,y,s} = \frac{F^g_{a,r,y,s}}{Z_{a,r,y,s}} N_{a,r,y,s} \left(1 - S_{a,r,y,s} \right) \end{equation}\]

Observed values were proportions with respect to age and sex, final model fitted proportions were \[ {P}^g_{a,r,y,s} = \frac{{C}^g_{a,r,y,s}}{\sum_a \sum_s {C}^g_{a,r,y,s}}, \] and initially the multinomial likelihood was assumed

\[ \boldsymbol{X}^g_{r,y} \sim \text{Multinomial}\left(\boldsymbol{\widehat{P}}^g_{r,y}\right) \] where, \(\boldsymbol{X}^g_{r,y} = \boldsymbol{P}^g_{r,y}N^{eff}_{r,y}\) and \(\boldsymbol{P}^g_{r,y}\) is the observed proportions, \(N^{eff}_{r,y}\) is the effective sample size and \(\boldsymbol{P}^g_{r,y} = (P^g_{1,r,y,1}, \dots, P^g_{a_+,r,y,1}, P^g_{1,r,y,2}, \dots, P^g_{a_+,r,y,2})\) is the vector of observed proportions across all ages and sexs in year \(y\) and region \(r\), and \(\boldsymbol{\widehat{P}}^g_{r,y}\) is the model fitted values which have the same dimension (\(\sum_a \sum_s \widehat{P}^g_{a,r,y,s} = 1\)).

Catch at length

Catch at length observations are available for the trawl fishery. Model fitted values are derived by multiplying the catch at age (see above) by an age-length transition matrix denoted by \(\boldsymbol{A}^l_{y,s}\) (dimensions of \(\boldsymbol{A}^l_{y,s}\) are \(n_a \ \times \ n_l\) and its rows must sum to 1),

\[\begin{equation} \widehat{\boldsymbol{Cl}}^g_{r,y,s} = \left(\boldsymbol{A}^l_{y,s} \right)^T \ \times \ \widehat{\boldsymbol{C}}^g_{r,y,s} \end{equation}\] where, \(\widehat{\boldsymbol{C}}^g_{r,y,s}\) is a column vector of catch at age (dimension \(n_a \ \times \ 1\)) at the beginning of year \(y\) in region \(r\) for sex \(s\), and \(\widehat{\boldsymbol{Cl}}^g_{r,y,s}\) is a column vector of catch at length (dimension \(n_l \ \times \ 1\)).

Proportions at age

Survey age composition data denoted by \({N}^s_{a,r,y,s}\) is available for the longline survey where model fitted numbers are derived as,

\[\begin{equation} \widehat{N}^s_{a,r,y,s} = N_{a,r,y,s} \left(1 - \exp^{\delta_y Z_{a,r,y,s}}\right)S^s_{y,r,a,s} \end{equation}\] where, \(\delta_y \in (0,1)\) is the proportion of time in the year that the observation occurs during year \(y\) and \(S^s_{y,r,a,s}\) is the survey selectivity.

Relative abundance indices

\[\begin{equation} \widehat{I}^s_{r,y} = \sum_s\sum_a \widehat{N}^s_{a,r,y,s} \bar{w}_{a,y,s} \end{equation}\] where, \(\bar{w}_{a,y,s}\) is mean weight at age, this can be omitted if the observation is in numbers i.e., abundance instead of biomass.

Tag recovery observations

When tagged fish are released in to the model partition they are tracked by the tag release event \(k\) for \(n_{T}\) years. The model expects tag-recovery observations to be known by age and sex. These are derived for each recovery event (also year and region specific) and tag-release event by using the same age-length key that was used to release each recovered fish to obtain an sex specific age-frequency. This age-frequency is then aged by the number of years this tag-recovery event was at liberty to derive observed tag-recoveries by age and sex.

This assumes the age-length conversion between releases and recoveries is consistent and allows us to use recovered fish with no length or sex recorded, but has the down side as mentioned earlier in the tag release section of smearing a single recovered fish across multiple age bins and sexes. However, it does mitigate the problem highlighted in chapter 7 (tag release section) of going backwards and forwards through the age-length transition matrix, which appears to be more problematic.

All tag recoveries are from the fixed gear (longline and pot gear) fishery, due to the recent trend in recoveries by gear method (Figure 7.2), we only considered tags recovered before 2018. The switch from longline to pot fishing in the fixed gear fishery has not seen anywhere near the same recoveries to be reported from the pot gear compared to the longling gear. To avoid this complication we will only consider recoveries from a period when the longline method was the dominant method of fixed gear catch

Model fitted values for tag event \(k\) recovered in region \(r\) and year \(y\) denoted by the set index \(m = \{r,y\}\) are denoted by \(\widehat{N}^k_{m}\), and are derived as

\[ \widehat{N}^k_{m} = \sum_a \sum_s T^k_{a,r,y,s} \frac{F^{LL}_{a,r,y,s}}{Z_{a,r,y,s}} (1 - e^{-Z_{a,r,y,s}})\delta_{r,y}, \quad m = \{r,y\} \] where, \(T^k_{a,r,y,s}\) are the number of tagged fish from tag-release event \(k\), and \(\delta_{r,y}\) is the reporting rate in year \(y\). Initial model runs used the Poisson

\[ N^k_{m} \sim \mathcal{Poisson}(\widehat{N}^k_{m}) \]

However, the negative binomial and multinomial are in development. See Chapter 13 for simple exploration of different tag-likelihoods for estimating movement parameters.