TagIntegrated model equations

Process equations

Population dynamics

The order of processes in an annual cycle follow

1. recruitment and release of tagged fish (we apply initial tag induced mortality here),
2. total mortality and ageing,
3. markovian movement
4. annual tag shedding (applied as a mortality process),

Before movement, the partition is updated following \[\begin{align*} N_{a,r,y,s} = \begin{cases} R_{r,y} p^R_{s,y}, & 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*}\]

\(p^R_{s,y}\) is proportion recruits for sex \(s\) in year \(y\) and movement is then applied

\[\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.

** An important input is whether movement occurs to new recruits, which is specified by the input do_recruits_move. If do_recruits_move == 0 i.e., they don’t get movement applied. In the code this applied by applying movement and then resetting the numbers at age by the recruited values.

Initialisation

An equilibrium age structure is derived by iterating the annual cycle is run \(n_a - 1\) times (i.e., 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 global or region specific deviations to allow the model to start with a non-equilibrium age-structure

\[ N_{a, r} = N_{a, r} e^{\psi_a} \ \ \forall \ r \ \& \forall a \in (a_{min} + 1, \dots, a_{max} - 1) \] These estimatble intiial age-deviations are not applied on the first or last age cohort. This is what is done in the current Sablefish assessment, I think it has something to do with helping estimation of \(R_0\) and such. There is the ability to estimate less age-deviations than age using data$n_init_rec_devs. If there are 30 age-cohorts in the model you can set up initial age-deviations so only 15 are estimated with ages 16-29 getting applied the same last \(e^{\psi_a}\) value.

To help with estimation, there is a penalty on \(\psi_a\) that assumes a central tendancy of zero with an estimable variance parameter (\(\sigma_{\psi}^2\)). \[ \psi_a \sim \mathcal{N}(0, \sigma_{\psi}) \]

In general these parameters can be highly uncertain. The variance on these deviation parameters is often fixed at a low value which requires a strong signal in the data to deviate away from the initial equilibrium age-structure.

You can also apply/estimate additional fishing mortality during the initial equilibrium calculation denoted by \(F^{init}_{a,s}\) for age \(a\) and sex \(s\). If F_method = 0 then \[ F^{init}_{a,s} = F^{fixed} S_{a,s}^{fixed} p_{hist} \] where, \(F^{fixed}\) is the average annual fishing mortality for the fixed gear (ln_fixed_F_avg), \(S_{a,s}^{fixed}\) is the fixed gear selectivity in the first time-block and \(p_{hist}\) is an input value that specifies the proportion of average longline fishing mortality to apply during initialisation. If F_method = 1 then \[ F^{init} = \widehat{F}^{init}S_{a,s}^{fixed} p_{hist} \] where, \(\widehat{F}^{init}\) is the estimable parameter ln_init_F_avg. It is advised to set \(p_{hist} = 1.0\) so that the estimated parameter has a natural interpretation.

There are two initial biomass quantities that are reported, these are Bzero and Binit. Bzero is the spawning biomass from a population that has been exposed to natural mortality, recruitment (\(R_0\)), and movement only. Binit is the spawning biomass from a population that has been exposed to natural mortality and initial fishing mortality (\(F^{init}\)), recruitment (\(R_0\)), and movement.

Growth

Empirical length at age matrices are supplied for all years where sufficient age-length data was available and growth was assumed the same across all regions. Mean weight at age is also an user input. Growth cannot vary across regions in the current models, so users must supply length at age matrices and weight at age vectors for each year, age and sex in the model.

Recruitment

How to spatially apportion recruits in a way that is consistent with the data and isn’t confounded with movement of young fish? There are two considerations, firstly looking at the AF and LFs there are regions that seem to have more young fish than others (link to figure and describe). However, there is very little information on spawning grounds and behavior. This coupled with a complex early life history make it difficult to a priori set regional apportionment of recruits, which would be ideal, an alternative is to estimate these as free parameters. This may introduce confounding when we start estimating movement which is believed to be ontogenetic.

Terminology regarding recruitment and spatial population structures can be ambiguous and confusing. There is a fair bit of literature describing and comparing metapopulations, panmictic populations and spatially heterogeneous populations. I believe this model can represent the later two A metapopulation implies there is some natal identifiability when fish are mixing which this model doesn’t do. However, this would only be advantageous when there is some sort of natal dynamic e.g. natal homing or natal productivity difference such as growth. Neither of these are really evident in the observed data for sablefish and so having a natal attribute in the partition doesn’t make much sense.

Spawning biomass

Spawning biomass is based on mature female weight. SSB in a given year and region is calculated as,

\[ SSB_{y,r} \sum\limits_a N_{a,s = 2, r, y} exp(-Z_{a,s = 2, r, y})^{prop_Z} S^{mat}_{a,y} \bar{w}_{a,s = 2, y} \]

where, \(s = 2\) denotes the female sex, \(S^{mat}_{a,y}\) is the maturity at age, \(\bar{w}_{a,s = 2, y}\) is female mean weight at age and \(prop_Z\) indicates the proportion of total mortality applied before SSB is calculated.

Fishing mortality

When the hybrid \(F\) method is assumed, tagged fish were not included when internally solving the fishing mortality nuisance parameters. 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. Often tagging data has only length at release information, we have used age-length keys to convert numbers at length of release to numbers at age and sex at release that are then input into the model (data$male_tagged_cohorts_by_age and data$female_tagged_cohorts_by_age). Tagged fish from release event \(k\) are denoted in the partition by \(N^k_{a,r,y,s}\), and are tracked for \(n_{T}\) (data$n_years_to_retain_tagged_cohorts_for) years before migrating into an accumulation tag group, at which point we loose release-year information but do maintain release region. At present, tagged fish are assumed to take on the exact same population processes as the untagged elements of the partition i.e., no mixing periods.

• Something to consider, given there is uncertainty in the age-length key method, when converting numbers at length to numbers at age and sex and the numbers of tag-releases are assumed known without error. Consider rerunning the model with bootstrapped tag releases from bootstrapped age-length keys. This will give you an indication of the uncertainty the model has to this input assumption.

Selectivities

Describe the functional forms available in the models. The selectivity component is quite modular and you can easily add a new selectivity function in the source code here. There are three selectivities in the model one for the survey and two for the fisheries. Selectivities types are set by the model data objects ending with sel_type e.g., fixed_sel_type, trwl_sel_type and srv_dom_ll_sel_type. These objects are vectors each element defines a selectivity type that is applied in a time-block. The time-block is defined by the corresponding sel_by_year_indicator objects. The following is an example of how to specify and interpret these containers. Say data$years = c(1990, 1991, 1992, 1993) and we were not interested in doing projections (require more elements in sel_by_year_indicator objects). Say we wanted to have 1990 and 1991 have one logistic selectivity and then 1992 and 1993 to have seperate logistic selectivity for the fixed gear fishery, you would specify the containers as

• data$fixed_sel_type = c(0, 0)
• data$fixed_sel_by_year_indicator = c(0, 0, 1, 1) this points each model year to a selectivity in data$fixed_sel_type

Selectivity types are

• 0 Logistic \[ S(a | \theta_1, \theta_2) = \frac{1}{1 + \exp(-\theta_2 \times (a - \theta_1)} \]

• 1 Double normal - Punt et. al 1996 gamma parameterization (that is a comment from the ADMB assessment code) \[ S(a | \theta_1, \theta_2) = \frac{a}{\theta_1}^{\theta_1/0.5\sqrt{\theta_1^2 + 4\theta_2^2}} exp\left(\frac{\theta_1 - a}{0.5 \left( \sqrt{\theta_1^2 + 4\theta_2^2} - \theta_1 \right)} \right) \]

• 2 power function (Be careful with this selectivity, if your min age is not = 1. Then I don’t think you will end up with a max selectivity value of 1) \[ S(a | \theta_1) = \frac{1}{a ^{\theta_1}} \]

• 3 Alternative logistic \[ S(a | \theta_1, \theta_2) = \frac{1}{1 + 19^{(\theta_1 - a) / \theta_2}} \]

• 4 Exponential decay \[ S(a | \theta_1) = exp(-a * \theta_1); \]

• 5 Double normal 3 parameterisation, where \(\boldsymbol{\theta} = (\theta_1, \theta_2, \theta_3) = (\mu, \sigma_r, \sigma_l)\), where \(\mu\) is the age that the selectivity = 1, \(\sigma_r\) is the standard deviation for the right hand side of the selectivity curve and \(\sigma_l\) is the standard deviation for the left hand side of the selectivity curve \[ S(a | \mu, \sigma_r, \sigma_l) = j_a \exp\left(\log(0.5) \left(\frac{a - \mu}{\sigma_r^2}\right)^2\right) + \left( 1 - j_a \right)\exp\left(\log(0.5) \left(\frac{a - \mu}{\sigma_l^2}\right)^2\right) \] where, \[ j_a = 1 / \left(1 + \exp\{-5(a - \mu)\}\right) \]

Observation equations

There are five observation types available in the TagIntegrated model

• Relative indices of abundance from a longline survey (its called longline survey but could be any survey)
• Age composition disaggregated by sex for the fixed gear fishery and longline survey
• Length composition disaggregated by sex for the trawl and fixed gear fishery
• Annual observed catch by fishery
• Tag recoveries from the fixed gear fishery

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}}, \]

if the multinomial distribution is 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 using data$srv_dom_ll_obs_is_abundance = 1. There are currently two likelihoods available for this observation

• if data$srv_dom_ll_bio_comp_likelihood == 0 \[ ll_{r,y} = \frac{\left( \log {I}^s_{r,y} + 0.0001 - \log \widehat{I}^s_{r,y} + 0.0001 \right)^2}{2 \left( SE[{I}^s_{r,y}] / {I}^s_{r,y} \right)} \]

• if data$srv_dom_ll_bio_comp_likelihood == 1 \[ ll_{r,y} = \frac{1}{{I}^s_{r,y} SE[{I}^s_{r,y}] \sqrt{ 2\pi}} \exp \left( - \frac{ \left(\log {I}^s_{r,y} - \log \widehat{I}^s_{r,y} - 0.5SE[{I}^s_{r,y}]^2\right)^2}{2 SE[{I}^s_{r,y}]^2}\right) \]

The catchability coefficient can also be calculated as a nuisance parameter instead of a free estimable parameter using the input data$q_is_nuisance = 1. This will solve the value for q in each region conditional on the other input free parameters. If the abundance or biomass observation is assumed to be lognormally distributed, then we can algebraically solve for the catchability coefficient.

For simplicity in the following equations we use \(E_i\) to denote the model expected value \(O_i\) to denote the observed value and \(\sigma_i\) to denote the standard error for element \(i\) of a series. In this case a time-series is within each region. So we will calculate a nuisance q for each region.

Technically, this derivation is only correct when data$srv_dom_ll_bio_comp_likelihood == 1 however, it will be considered an approximation for the other likelihood. Assuming a specific element contributes to the negative log likelihood in the following manor,

\[ -ll = \sum_i log (\sigma_i) + 0.5\sum_i \bigg(\frac{log(O_i) - log(qE_i) + 0.5\sigma_i^2}{\sigma_i} \bigg)^2 \ , \]

then, if you differentiate this term with respect to \(q\) you get \[ \frac{\partial }{\partial q}(-ll) = \frac{-1}{q} \sum_i\bigg( \frac{log(O_i/E_i) - log(q) + 0.5\sigma_i^2}{\sigma_i^2}\bigg) \, \]

if you solve for when the derivative is zero (MLE estimate) then you obtain the following result

\[ \hat q = exp\frac{0.5n + S_3}{S_4} \, \] where \(S_3 = \sum_i (log(O_i /E_I)/\sigma_i^2)\) and \(S_4 = \sum_i(1/\sigma_i^2)\). \(E_i\) in the initial calculation is derived assuming a \(q = 1\), before \(\hat q\) is calculated.

Tag recovery observations

When tagged fish are released in to the model partition they are tracked by the tag release event index denoted by \(k\) for \(n_{T}\) years after release. Due to tags being predominately recovered by the fixed gear fishery, all recoveries are assumed from the fixed gear fishery. A tag recovery observation will depend on if it is release conditioned or recapture conditioned. If tag-observations are release conditioned (the default) then for each release event there is a maximum of \(n_{T} \times n_r + 1\) recovery possibilities, assuming the fixed gear fishery operates in all regions and all years we tracked the release event for. The plus one is for the not-recovered event, although it is not clear to me how this is different from a zero recovery event? (ask Dan). Let \(m\) index a recovery event for a given release event and let \(n_{m|k} = n_{T} \times n_r + 1\) where \(n_{m|k}\) is the number of recovery events for release event \(k\). Let \(N^k_{m}\) denote the number of tag-recoveries for a given release and recovery combination, and \(\boldsymbol{N^k} = (N^k_{1}, N^k_{2}, \dots, N^k_{n_{m|k}})^T\) represents a vector of all possible observed recoveries for release event \(k\). Model fitted values for \(N^k_{m}\) are calculated as

\[ \widehat{N}^k_{m} = \sum_a \sum_s N^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_y \ \ y,r \in m|k \] where, \(N^k_{a,r,y,s}\) is the number of tagged fish alive from tag-release event \(k\), and \(\delta_y\) is the reporting rate in year \(y\).

There are three likelihoods used in the literature for tag-recoveries,

• The Poisson (tag_likelihood = 0) \[ N^k_{m} \sim \mathcal{Poisson}(\widehat{N}^k_{m}) \]
• The Negative Binomial (tag_likelihood = 1) \[ N^k_{m} \sim \mathcal{NB}(\widehat{N}^k_{m}, \hat{\phi}) \] where, \(\hat{\phi}\) is the estimable dispersion parameter
• Multinomial (tag_likelihood = 2) \[ \boldsymbol{N^k} \sim \text{Multinomial}\left(\boldsymbol{\hat{\theta}}^k, \text{Neff}^k\right) \] All the recovery states include all regions, and years of recoveries including the non recaptured, where, \(\boldsymbol{\hat{\theta}}^k = (\hat{\theta}^k_1, \hat{\theta}^k_2\dots, \hat{\theta}^k_{n_{m|k}})\) is the proportions of all recovery states (Vincent, Brenden, and Bence 2020; Goethel, Legault, and Cadrin 2014), with

\[\begin{equation} \hat{\theta}^{k}_{m} = \frac{\widehat{N}^k_{m}}{\sum_m \widehat{N}^k_{m}} \end{equation}\]

For recapture conditioned tag-recovery observations, we calculate the probability of capture a fish at time \(t\) being released in region \(r\) and recovered in region \(r'\). This was explored by McGarvey and Feenstra (2002). In this study they show how many of the terms can be dropped when describing the probability of that recovery event. The term \(S_i[a_t, a_m]\) represents the probability of a tagged individual tagged at \(a_t\) surviving until \(a_m\). In our model this surviorship will be regional specific, i.e., fish in certain regions will be exposed to higher fishing mortality rates than other regions over time.

Here, model predicted recaptures are recapture conditioned, which means the predicted proportions are year specific, and represent the proportion of recaptures in each region, \(\boldsymbol{\hat{\theta}}^k_y = (\hat{\theta}^{k}_{y,1},\dots, \hat{\theta}^k_{y,n_r})'\) \[\begin{equation} \hat{\theta}^{k}_{y,r} = \frac{\widehat{N}^k_{r,y}}{\sum_r\widehat{N}^k_{r,y}} \end{equation}\]

• How to deal with the high number of zeros in \(\boldsymbol{N^k}\) zero inflated?

These are derived for each recovery event (also year and region specific) denoted by the index \(m\) and tag-release event by using the 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 of going backwards and forwards through the age-length transition matrix, which appears to be more problematic.

• Multinomial (Not implemented) \[ {N}^k_{r,y} \sim \text{Multinomial}\left(\boldsymbol{\hat{\theta}}^k_y, \text{Neff}^k_y\right) \] Here, model predicted recaptures are recapture conditioned, which means the predicted proportions are year specific, and represent the proportion of recaptures in each region, \(\boldsymbol{\hat{\theta}}^k_y = (\hat{\theta}^{k}_{y,1},\dots, \hat{\theta}^k_{y,n_r})'\) \[\begin{equation} \hat{\theta}^{k}_{y,r} = \frac{\widehat{N}^k_{r,y}}{\sum_r\widehat{N}^k_{r,y}} \end{equation}\]

Release conditioned are defined as the proportion of recaptures over all states denoted as \(i\), for each release event. All the states, include all regions, and years of recoveries including the non recaptured, where, \(\boldsymbol{\hat{\theta}}^k = (\hat{\theta}^k_i, \hat{\theta}^k_2\dots, \hat{\theta}^k_{n_i}) = (\hat{\theta}^k_{1,1}, \hat{\theta}^k_{1,2}, \hat{\theta}^k_{2,2}, \dots, \hat{\theta}^k_{n_y,n_r}, \hat{\theta}^k_{NR})\) is the proportions of all states (Vincent, Brenden, and Bence 2020; Goethel, Legault, and Cadrin 2014), with

\[\begin{equation} \hat{\theta}^{k}_{i} = \frac{\widehat{N}^k_{i}}{\sum_i \widehat{N}^k_{i}} \end{equation}\]

Likelihood formulations

Multinomial

This likelihood has no estimable parameters

\[ \boldsymbol{X} \sim \text{Multinomial}\left(\boldsymbol{\widehat{P}}\right) \]

To avoid \(\widehat{P}_i = 0\), in which case the multinomial can be undefined if \(X_i = 0\) all values of \(\widehat{P}_i\) are run through a posfun function which adds a small constant and penalty to the likelihood. See this TMB issue here or Ben Bolker’s write up found here.

Dirichlet-Multinomial

The Dirichlet-Multinomial distribution used for composition data, with input sample size denoted by \(n\), observed proportion for age \(a\) denoted by \(p_a\) and expected proportions denoted by \(\widehat{p}_a\). There is an estimable parameter denoted by \(\theta\). This follows the linear-parameterisation from Thorson et al. (2017), with the following density function

\[\begin{equation} f(x) = \frac{\Gamma\left(n + 1\right)}{\sum_a n{p}_a + 1}\frac{\Gamma\left(\theta n\right) }{\Gamma\left(n + \theta n\right)}\prod_{a} \frac{\Gamma\left(n {p}_a + n \theta \widehat{p}_a\right)}{\Gamma\left(n \theta \widehat{p}_a\right)} \end{equation}\]

The effectictive sample size is calcualted as,

\[ n_{eff} = \frac{1 + \theta n}{1 + \theta} \]

Simulating from this distribution is done by first simulating a dirichlet variable, using independent gamma draws (normalised), with shape parameter set as the expected proportion. Then draw from the multinomial with expected proportions based on the dirichlet draw.

Lognormal

As mentioned in the abundance observation section there are two alternative lognormal forms in the model. Pearon residuals are calculated as

\[ r_i = \frac{log(O_i/E_i) + 0.5\sigma_i^2}{\sigma_i} \]

Projecting

Recruitment

There are two methods for future recruitment, parametric and empirical. Parametric will simulate from the the lognormal distribution with \(\sigma_R\). Empirical will resample input recruitment deviations between certain years

• Parametric future_recruitment_type == 0 \[ \epsilon_y \sim \mathcal{N}\left(0, \sigma_R^2 \right) \] resulting recruitment multipliers, also termed year class strengths (\(YCS_y\)) are bias adjusted \(YCS_y = \exp(\epsilon_y - 0.5\sigma_R^2)\)

• Empirical future_recruitment_type == 1 this will sample with replacement input recruitment deviations between the values specified in year_ndx_for_empirical_resampling. Note elements of year_ndx_for_empirical_resampling are C++ vector dimensions. If year_ndx_for_empirical_resampling = c(0,n_years - 1) this will resample from all input recruitment deviations. if year_ndx_for_empirical_resampling = c(0,9) this will resample from the first 10 input recruitment deviations. if year_ndx_for_empirical_resampling = c((n_years - 9):(n_years - 1)) this will resample from the last 10 input recruitment deviations.

• Deterministic recruitment future_recruitment_type == 2. This will assume \(YCS_y = 1\) which result in the model applying mean recruitment for all future years

Fishing

Fishing selectivity that was assumed in the last year of the model is used during the projection period. There are two methods users can use for future fishing, assume fishing mortality rates or catches.

• User specifies F values during projection (unlikely useage) future_fishing_type = 0 populate future_fishing_inputs_trwl and future_fishing_inputs_fixed with F values

• User specifies Catch values during projection (more likely) future_fishing_type = 1 populate future_fishing_inputs_trwl and future_fishing_inputs_fixed with future catches. An F is calculated based on the hybrid fishing mortality method, which should result in predicted catches close to these but not exact.

Example code snippets

## once you have done MLE and evaluated fits to a model
## then you can running the following code to find
## deterministic reference points based on %B0
n_proj_years = 100
regional_spr = find_regional_Fspr(data = data, 
                                  MLE_report = mle_report, n_years_for_fleet_ratio = 5, 
                                  percent_Bzero = 40, n_future_years = n_proj_years, verbose = T)

# setup projection data
proj_data = setup_proj_data(mle_obj = mle_obj, n_proj_years  = 100)
## use the Fs from Fspr function
proj_data$future_fishing_inputs_trwl = matrix(regional_spr$Fspr * (1 - regional_spr$fixed_gear_F_proportion), 
                                              byrow = F, ncol = n_proj_years, nrow = n_regions)
proj_data$future_fishing_inputs_fixed = matrix(regional_spr$Fspr * regional_spr$fixed_gear_F_proportion, 
                                               byrow = F, ncol = n_proj_years, nrow = n_regions)

# check projection data
validate_input_data_and_parameters(data = proj_data, parameters = parameters)
# build proj object
proj_obj <- MakeADFun(proj_data, parameters, map = map_fixed_pars,
                      DLL="SpatialSablefishAssessment_TMBExports", hessian = T, silent = T)
# run the projection model with future F's and MLE parameter values
proj_rep = proj_obj$report(mle_spatial$par)
# Check the depletion levels are as expected
plot_SSB(MLE_report = proj_rep, region_key = region_key, depletion = T) +
  geom_hline(yintercept = 40, col = "gray", linetype = "dashed", linewidth = 1.1) +
  ylim(0, NA)

## once this checks out, you should run projections with stochasticity in recruitment and parameters estimates
## to get projected quantities with uncertainty.

Symbol Notation

Symbol	Description
\(y\)	Year index
\(T\)	Terminal year of the model
\(s\)	Sex index \(s \in \{1,2\}\)
\(a\)	Model age cohort, i.e., \(a = a_0, a_0 + 1, \dots\)
\(a_{1}\)	Recruitment age to the model = 2
\(a_+\)	Plus-group age class (oldest age considered plus all older ages)
\(n_a\)	Number of age classes modeled \(a_+ \ - a_1\)
\(l\)	length class
\(n_l\)	Number of length classes
\(g\)	gear type index, i.e. longline survey, fixed gear fishery, trawl fishery
\(x\)	log-likelihoos index
\(\bar{w}_{a,y, s}\)	Average weight at age \(a\), year \(y\) and sex \(s\)
\(\phi_{a,y}\)	Proportion of female mature by age and year
\(p^s_{y}\)	Proportion of recruits for sex \(s\). Often assumed = 0.5
\(\ln \mu_{r}\)	Average log-recruitment
\(p^R_y\)	proportion recruitment male
\(\ln \mu_{f}\)	Average log-fishing mortality
\(\phi_{y,g}\)	annual fishing mortality deviation by gear (log space)
\(\tau_{y}\)	annual recruitment deviation \(\sim LogNormal\left(0,\sigma_r\right)\)
\(\sigma_r\)	Recruitment standard deviation
\(N_{a,y,s}\)	Numbers of fish at age \(a\) in year \(y\) of sex \(s\)
\(M\)	Natural mortality
\(F^g_{a,y}\)	Fishing mortality for year \(y\), age \(a\) and gear \(g\)
\(F_{hist}\)	Historical proportion of Fishing mortality
\(Z_{a,y}\)	Total mortality for year \(y\), age \(a\) \(=\sum\limits_g F^g_{a,y} + M\)
\(R_{y}\)	Annual recruitment
\(B_{y}\)	Spawning biomass in year \(y\)
\(S^g_{a,y,s}\)	Selectivity at age \(a\) for gear type \(g\) and sex \(s\)
\(a_50\)	age at 50% selection for ascending limb
\(d_50\)	age at 50% selection for descending limb
\(\delta\)	slope/shape parameters for different logistic curves
\(\boldsymbol{A}\)	ageing-error matrix dimensions \(n_a \ \times \ n_a\)
\(\boldsymbol{A}^l_s\)	age to length conversion matrix by sex. dimensions \(n_a \ \times \ n_l\)
\(q_g\)	abundance index catchability coeffecient by gear
\(\lambda_x\)	Statistical weight (penalty) for component \(x\)
\(P^g_{l,y,s}\)	Observed proportions at length for gear \(g\) in year \(y\) and sex \(s\)
\(P^g_{a,y,s}\)	Observed proportions at age for gear \(g\) in year \(y\) and sex \(s\)

References

Goethel, Daniel R, Christopher M Legault, and Steven X Cadrin. 2014. “Demonstration of a Spatially Explicit, Tag-Integrated Stock Assessment Model with Application to Three Interconnected Stocks of Yellowtail Flounder Off of New England.” ICES Journal of Marine Science 72 (1): 164–77.

McGarvey, Richard, and John E Feenstra. 2002. “Estimating Rates of Fish Movement from Tag Recoveries: Conditioning by Recapture.” Canadian Journal of Fisheries and Aquatic Sciences 59 (6): 1054–64.

Thorson, James T, Kelli F Johnson, Richard D Methot, and Ian G Taylor. 2017. “Model-Based Estimates of Effective Sample Size in Stock Assessment Models Using the Dirichlet-Multinomial Distribution.” Fisheries Research 192: 84–93.

Vincent, Matthew T, Travis O Brenden, and James R Bence. 2020. “Parameter Estimation Performance of a Recapture-Conditioned Integrated Tagging Catch-at-Age Analysis Model.” Fisheries Research 224: 105451.