Description of Model Equations

Process Models and Equations

SPoRC assumes an annual time-step, where the following processes are applied in order:

1. Recruitment and tag releases initially occur (tag releases occur if tagging data are used),
2. Markovian movement of all individuals then occurs,
3. Total mortality, chronic tag loss (if applicable), and ageing processes then take place,

The aforementioned processes are applied to four main population partitions, which include region ($r$), year ($y$), age ($a$ and $a_+$, where the latter indicates the plus group), and sex ($s$). Although SPoRC can accommodate spatial process, all of the equations described in this vignette can collapse and be generalized to a single region (i.e., $r = 1$).

Recruitment Processes

Recruitment processes occur at the beginning of the year and can be described by either mean recruitment or with a Beverton-Holt stock recruitment relationship. Annual recruitment generated by a mean recruitment relationship is governed by the following equation:

$$Rec_{r,y,s} = \mu^{Rec}\exp\left(\epsilon_{r,y}^{Rec}-\frac{\sigma^2_{Rec}}{2}b_y\right)\psi_s\zeta_r$$

where $Rec_{r,y,s}$ represents recruitment across regions, years, and sexes. $\mu^{Rec}$ represents the mean recruitment parameter, $\epsilon_{r,y}^{Rec}$ are the recruitment deviations constrained by a normal distribution. $b_y$ represents values determined from a recruitment bias correction ramp, $\psi_s$ is the recruitment sex ratio, and $\zeta_r$ are the recruitment apportionment parameters. Recruitment parameterized as a Beverton-Holt relationship assuming localized density-dependence can be expressed as:

$$Rec_{r,y,s} = \frac{4R0\zeta_rh_rSSB_{r,y-\text{RecLag}}}{(1-h_r)SSB0_r+5(h_r-1)SSB_{r,y-\text{RecLag}}}\exp\left(\epsilon_{r,y}^{Rec}-\frac{\sigma^2_{Rec}}{2}b_y\right)\psi_s$$

where $R0$ is the virgin unfished recruitment, $h_r$ is the steepness parameter representing the fraction of $R0\zeta_r$ that would be produced when at 20% of $SSB0_r$. The steepness parameter is constrained to be between values of 0.2 and 1 and are estimated in bounded logit space:

$$h_r = 0.2 + (1 - 0.2) \frac{\exp(h^{'}_r)}{1 + \exp(h^{'}_r)}$$

where $h^{'}_r$ is the steepness parameter in bounded logit space. $SSB0_r$ represents the unfished spawning stock biomass computed as the spawning biomass per recruit multiplied by $R0\zeta_r$ and $RecLag$ represents the delay between spawning and when recruits enter the population. If $RecLag > y$, the model utilizes $SSB0_r$ instead of $SSB_(r,y-RecLag)$ to compute deterministic recruitment. $SSB_(r,y)$ is the spawning stock biomass, which is computed as:

$$SSB_{r,y} = \sum_a^{a+} N_{r,y,a,s=1}W_{r,y,a,s=1}Mat_{r,y,a,s=1}$$

$N_{r,y,a,s=1}$ here represents the numbers at age across regions and years for females ($s = 1$), $W_{r,y,a,s=1}$ are values of weight at age for females, and $Mat_{r,y,a,s=1}$ are values of maturity at age. In a single-sex model, $SSB_{r,y}$ is simply computed as: $SSB_{r,y} = \sum_a^{a+} N_{r,y,a,s=1}W_{r,y,a,s=1} \cdot 0.5$.

In a spatial model, regional recruitment apportionment parameters $\zeta_r$ are derived using a multinomial logit transformation. There are also additional options that can be specified for recruitment processes in a spatial model. In particular, users can also specify whether recruitment occurs globally instead of locally. In such a case, the Beverton-Holt relatinoship is rewritten as:

$$Rec_{r,y,s} = \frac{4R0h\sum_rSSB_{r,y-\text{RecLag}}}{(1-h)SSB0+5(h-1)\sum_rSSB_{r,y-\text{RecLag}}}\exp\left(\epsilon_{y}^{Rec}-\frac{\sigma^2_{Rec}}{2}b_y\right)\psi_s\zeta_r$$

where $SSB$ is summed across regions, recruitment apportionment parameters are applied to the overall deterministic recruitment relationship (as opposed to $R0$ in the numerator), and recruitment deviations are estimated globally.

Population Projection

Following recruitment processes, the population can then be projected forward. Movement processes first occur: $$\boldsymbol{N}{_{y,a,s}} = (\boldsymbol{N}{_{y,a,s}})^T \boldsymbol{M}{_{y,a,s}}$$

where $\boldsymbol{M}{_{y,a,s}}$ is a first-order Markov matrix representing movement. In a single-region case, no movement occurs (i.e., $\boldsymbol{M}{_{y,a,s}}$ is an identity matrix). Movement for a given year, age, and sex combination is constrained by a multinomial logit transformation:

$$M_{k,j,y,a,s} = \frac{\exp(\omega_{k,j,y,a,s})}{\sum_{k=1}^J\exp(\omega_{k,j,y,a,s})}, \omega_{j,1,y,a,s} = 0$$

where $\omega_{k,j,y,a,s}$ are movement parameters in multinomial logit space. Given that movement parameters must sum to 1 across rows, only $n_r - 1$ movement parameters are estiamted for a given year, age, and sex combination. After the population undergoes movement, individuals then experience mortality and ageing processes following an exponential mortality model:

$$\begin{aligned} N_{r,y+1,a+1,s} = N_{r,y,a,s}\exp(-Z_{r,y,a,s}), \quad \text{for } 1 \leq a < a_+ \\ N_{r,y+1,a,s} = N_{r,y,a-1,s}\exp(-Z_{r,y,a-1,s}) + N_{r,y,a,s}\exp(-Z_{r,y,a,s}), \quad \text{for } a = a_+ \end{aligned}$$

Here, $Z_{r,y,a,s}$ denotes the total instantaneous mortality:

$$Z_{r,y,a,s} = Natmort_{r,y,a,s} + \sum_{f} Fmort_{r,y,f} Sel_{r,y,a,s,f}^{Fsh}$$

where $Natmort_{r,y,a,s}$ is the instantaneous natural mortality rate, $Fmort_{r,y,f}$ is the annual instantaneous fishing mortality rate and for fishery fleet $f$, and $Sel_{r,y,a,s,f}^{Fsh}$represents fishery selectivity processes. $Fmort_{r,y,f}$ is determined by:

$$Fmort_{r,y,f} = \mu^{Fmort}_{r,f} \cdot exp(\epsilon_{r,y,f}^{FishMort})$$

$\mu^{Fmort}_{r,f}$ are mean fishing mortality parameters for a given region and fleet, and $\epsilon_{r,y,f}^{FishMort}$ are deviations from the mean, constrained by a normal distribution.

Population Initialization

Several options are available for initializing the population age structure. In the first approach, the initial population age structure is derived by iterating the population to equilibrium:

$$\begin{aligned} Z_{r,a,s}^{\text{Init}} &= \text{Natmort}_{r,1,a,s} + F_{\text{mort,InitProp}} \cdot \mu_{f=1} \cdot \text{Sel}^{\text{Fsh}}_{r,1,a,s,f=1}, \quad \text{for } a < a_+ \\ \text{cumsum}_{Z_{r,a,s}^{\text{Init}}} &= \sum_{a'=1}^{a} Z_{r,a',s}^{\text{Init}} \end{aligned}$$

$Z_{r,a,s}^{Init}$ is the initial instantaneous total mortality rate determined by the sum of natural mortality and an inital fishing mortality rate. The initial instantaneous fishing mortality rate is a user-specified value for $Fmort_{InitProp}$, which represents the proportion of $\mu_{f=1}$ to apply to the initial age structure. In SPoRC, it is generally recommended to structure associated data inptus for fishery fleet 1 to represent the dominant fishery fleet. $cumsum_{Z_{r,a,s}^{Init}}$ is the cumulative sum of the initial total mortality rate which is initially applied to the equilibrium age structure following an exponential motality model:

$$N_{r,a,s}^{Init} = \begin{cases} R0 \psi_s \zeta_r, & \text{if } a = 1 \\ R0 \exp\left(-cumsum_{Z_{r,a,s}^{Init}}\right)\psi_s\zeta_r, & \text{if } a > 1 \end{cases}$$

If mean recruitment dynamics are assumed, $R0$ in the above equation can simply be replaced with $\mu^{Rec}_{r}$. Following determination of the equilibrium age structure, the initial population is then iterated $n_a \cdot 5$ times where movement is first applied, followed mortality and ageing processes. In general, users might want to iterate the equilibrium age structure if movement processes occur given that there is no simple closed form solution to derive the plus-group of the population when movement occurs.

Conversely, users can also specify an initial age structure determined by a geometric series solution:

$$N_{r,a,s}^{Init} = \begin{cases} R0 \exp\left(-(a-1)(Natmort_{r,1,a,s} + Fmort_{init}\mu_{f=1}Sel^{Fsh}_{r,1,a,s,f=1})\right)\psi_s\zeta_r, & \text{if } 2 \leq a < a_+ \\ \frac{R0 \exp\left(-(a_+ - 1)(Natmort_{r,1,a,s} + Fmort_{InitProp} \cdot \mu_{f=1}Sel^{Fsh}_{r,1,a,s,f=1})\right)}{1 - \exp\left(-(Natmort_{r,1,a,s} + Fmort_{init}\mu_{f=1}Sel^{Fsh}_{r,1,a,s,f=1})\right)}\psi_s\zeta_r, & \text{if } a = a_+ \end{cases}$$

This approach is generally more appropriate when movement does not occur among regions, given that a closed form solution can easily be derived for the plus-group. Initial age deviations can then be applied to the equilibrium age structure:

$$N_{r,a,s}^{Init} = N_{r,a,s}^{Init} \exp(\epsilon_{r,a}^{Init})$$

where $\epsilon_{r,y}^{Init}$ are initial age deviations constrained by a normal distribution. Note that there are several options for initializing the age structure. In particular, users can specify equilibrium age structure, where deviations are assumed to be 0, stochastic age structure with the plus group following the equilibrium solution, and stochastic age structure with the plus group following the equilibrium solution but with an additional deviation applied to the plus group.

Selectivity Processes (Fishery and Survey)

Currently, several options are available for parameterizing fishery and survey selectivity. Firstly, selectivity can be specified as either being age-based or length-based. In the context of age-based selectivity, the usual calculations occur where an age vector is utilized to compute selectivity, given a specified functional form. In the case of length-based selectivity, a length vector is utilized to compute selectivity, after which the dot product of selectivity-at-length and a size-age transition matrix is computeed to derive selectivity-at-age:

$$\boldsymbol{Sel^a} = (\boldsymbol{A}^l_{r,y,s})^T \boldsymbol{Sel^l}$$

This derived selectivity-at-age is then used for all subsequent calculations, given that the model is an age-based model.

Various functional forms can also be specified for selectivity processes. In particular, two forms of logistic selectivity can be specified. Throughout this section, we will use $b$ to denote a genearlized bin. The first form is:

$$Sel_{b} = \frac{1}{1 + \exp[-k(b - b^{50})]}$$

where $k$ regulates the slope/steepness of the logistic curve and $b^{50}$ is the bin-at-50% selection. The second logistic form is expressed as:

$$Sel_{b} = \frac{1}{1 + 19^{\left(\frac{b^{50} - b}{b^{95}}\right)}}$$

Here, $ab^{50}$ is the bin-at-50% selection and $b^{95}$ is the bin-at-95% selection. Selectivity can also be parameterized as a gamma dome-shaped function:

$$\begin{aligned} p = 0.5 \left[ \sqrt{(b_{s,f}^{\text{max}})^2 + 4\gamma_{s,f}^2} - b_{s,f}^{\text{max}} \right] \\ Sel_{b} = \left( \frac{b}{b_{s,f}^{\text{max}}} \right)^{\left(\frac{b_{s,f}^{\text{max}}}{p}\right)} e^{\left( \frac{b_{s,f}^{\text{max}} - b}{p} \right)} \end{aligned}$$

where $p$ is a derived power parameter, $\gamma$ is the shape parameter the describes the steepness of the descending limb, and $b^{max}$ is an estimated parameter describing the bin-at-maximum selection. Other dome-shaped functions are also available. In particular, this includes the power function:

$$Sel_{b} = \frac{1}{b^{\phi}}$$

$\phi$ here is a power parameter determining the descending limb (larger values are steeper).

The other function that is available includes the 6 parameter double normal selectivity form. We define a parametric selectivity function using a flexible compound form. Parameters are transformed from input vector $\hat{p}_1$ through $\hat{p}_6$ to define the shape.

$$\begin{align*} p_1 &= \text{min}(b) + [\text{max}(b) - \text{min}(b)] \cdot \frac{1}{1 + \text{exp}(-\hat{p}_1)} && \text{Peak bin — beginning of the plateau} \\ p_2 &= p_1 + 1 + \frac{0.99 + \text{max}_x - p_1 - 1}{1 + \exp(-\hat{p}_2)} && \text{Width of plateau, smoothly approaches maximum bin} \\ p_3 &= \exp(\hat{p}_3), \quad p_4 = \exp(\hat{p}_4) && \text{Ascending and descending widths (on log scale)} \\ p_5 &= \frac{1}{1 + \exp(-\hat{p}_5)}, \quad p_6 = \frac{1}{1 + \exp(-\hat{p}_6)} && \text{Selectivity at first and last bins (logistic scale)} \end{align*}$$

Here, at low bins, the selectivity follows the ascending curve; at intermediate bins, it plateaus around 1; at higher bins, it transitions smoothly to the descending curve.

$$\begin{align*} \text{asc}(b) &= \exp\left( -\frac{(a - p_1)^2}{p_3} \right) && \text{Gaussian rise centered at } p_1 \\ \text{asc}_{\text{scaled}}(b) &= p_5 + (1 - p_5) \cdot \text{asc}(b) && \text{Rescales minimum to } p_5 \end{align*}$$

$$\begin{align*} \text{desc}(b) &= \exp\left( -\frac{(a - p_2)^2}{p_4} \right) && \text{Gaussian fall centered at } p_2 \\ \text{stj} &= \exp\left( -\frac{(40 - p_2)^2}{p_4} \right) && \text{Reference value for rescaling } \text{desc} \\ \text{desc}_{\text{scaled}}(b) &= 1 + (p_6 - 1) \cdot \frac{\text{desc}(b) - 1}{\text{stj} - 1} && \text{Descending scaled to reach } p_6 \text{ at } a = 40 \end{align*}$$

$$\begin{align*} \text{join}_1(b) &= \frac{1}{1 + \exp\left( -20 \cdot \frac{a - p_1}{1 + |a - p_1|} \right)} && \text{Logistic transition from asc to plateau} \\ \text{join}_2(b) &= \frac{1}{1 + \exp\left( -20 \cdot \frac{a - p_2}{1 + |a - p_2|} \right)} && \text{Logistic transition from plateau to desc} \end{align*}$$

$$\begin{align*} \text{Sel}(b) &= \text{asc}_{\text{scaled}}(b) \cdot (1 - \text{join}_1(b)) \\ &\quad + \text{join}_1(b) \cdot \left[1 \cdot (1 - \text{join}_2(b)) + \text{desc}_{\text{scaled}}(b) \cdot \text{join}_2(b)\right] \end{align*}$$

In addition to the various functional forms that can be specified to describe selectivity processes, several parameterizations of time-varying processes can be specified. In the simplest form, time-varying selectivity can be specified as time blocks, where new selectivity parameters are estimated for an arbitrary number of user-specified periods. Several options are also available for specifying continuous time-varying selectivity, where deviations are restricted to vary parametrically (i.e., deviations applied to the parameters of a given selectivity form) or semi-parametrically (i.e., deviations applied about the values of a selectivity form). For example, if logistic selectivity is specified and deviations are parametric, the following expression is invoked:

$$\begin{aligned} Sel_{y,b} = \frac{1}{1 + \exp[-k_{y}(b - b_{y}^{50})]} \\ k_{y} = k \exp(\epsilon_{y,1}^{Sel}) \\ b_{y}^{50} = b^{50} \exp(\epsilon_{y,2}^{Sel}) \end{aligned}$$

where the parameters of the logistic form are allowed to vary over time. By contrast, in the case of semi-parametric deviations for logistic selectivity, the following expression is used:

$$\begin{aligned} Sel_{y,b} = \frac{1}{1 + \exp[-k(b - b^{50})]} \exp(\epsilon_{y,b}^{Sel}) \\ Sel_{y,b} = \frac{Sel_{y,b}}{mean(\boldsymbol{Sel})} \end{aligned}$$

where deviations are placed about the parametric form and selectivity values are log mean standardized to aid with interpretability. Note that if age-based selectivity is specified and continuous deviations are estimated, these deviations are placed on the ages themselves. By contrast, if length-based selectivity is specified and continuous deviations are estimated, these deviations are placed on the lengths themselves. For further details on options available, see the selectivity process error section.

Observation Models and Equations

Fishery Observation Model

The fishery observation model describes the expected catch-at-age, catch-at-length, catch (in units of biomass), and fishery indices. Expected catch-at-age ($C^a_{r,y,a,s,f}$) is calculated using Baranov’s catch equation:

$$C^a_{r,y,a,s,f} = \frac{Fmort_{r,y,f}Sel^{Fsh}_{r,y,a,s,f}}{Z_{r,y,a,s}} N_{r,y,a,s} \left[1-\exp(-Z_{r,y,a,s})\right]$$

Catch-at-length ($C^l_{r,y,l,s,f}$) can then be derived using catch-at-age ($C^a_{r,y,a,s,f}$) and a size-age transition matrix ($\boldsymbol{A}^l_{r,y,s}$):

$$C^l_{r,y,l,s,f} = (\boldsymbol{A}^l_{r,y,s})^T \boldsymbol{C^a}_{r,y,s,f}$$

$\boldsymbol{A}^l_{r,y,s}$ is a data input derived external to the assessment model and is required if users are fitting to length composition data, where the matrix represents the probability of being in a length bin, given some age.

Expected catch ($\text{Catch}_{r,y,f}$) is then derived by simply taking the summation of the expected catch-at-age across ages and sexes, multiplied by their respective weights-at-age:

$$\text{Catch}_{r,y,f} = \sum_a^{a_+} \sum_s^{n_s} C^a_{r,y,a,s,f} W_{r,y,a,s}$$

Fishery indices ($\text{FshIdx}_{r,y,f}$) can be specified as either abundance-based (first equation below) or biomass-based (second equation below):

$$\begin{aligned} \text{FshIdx}_{r,y,f} = q^{Fsh}_{r,y,f} \sum_a^{a_+} \sum_s^{n_s} N_{r,y,a,s} \exp(-Z_{r,y,a,s}0.5) Sel_{r,y,a,s,f}^{Fsh} \\ \text{FshIdx}_{r,y,f} = q^{Fsh}_{r,y,f} \sum_a^{a_+} \sum_s^{n_s} N_{r,y,a,s} \exp(-Z_{r,y,a,s}0.5) Sel_{r,y,a,s,f}^{Fsh} W_{r,y,a,s} \end{aligned}$$

where $q^{Fsh}_{r,y,f}$ is the catchability coefficient for a given fishery and the value 0.5 represents the fraction of year in which the fishery index should be computed at. Currently, fishery catchability coefficients can only vary over time as time blocks.

Survey Observation Model

The survey observation model describes the expected survey catch-at-age, survey catch-at-length, and survey indices. Expected survey catch-at-age ($I^a_{r,y,a,s,sf}$) is calculated using the following equation:

$$I^a_{r,y,a,s,sf} = N_{r,y,a,s} \exp(-Z_{r,y,a,s}0.5) Sel_{r,y,a,s,sf}^{Srv}$$

where $Sel_{r,y,a,s,sf}^{Srv}$ here is the survey selectivity pattern. Expected survey catch-at-length ($I^l_{r,y,l,s,sf}$) can then be derived as:

$$I^l_{r,y,l,s,sf} = (\boldsymbol{A}^l_{r,y,s})^T \boldsymbol{I^a}_{r,y,s,f}$$

Again, survey indices ($\text{SrvIdx}_{r,y,sf}$) can be computed either as abundance-based (first equation below), or biomass-based (second equation below):

$$\begin{aligned} \text{SrvIdx}_{r,y,sf} = q^{Srv}_{r,y,sf} \sum_a^{a_+} \sum_s^{n_s} I^a_{r,y,a,s,sf} \\ \text{SrvIdx}_{r,y,sf} = q^{Srv}_{r,y,sf} \sum_a^{a_+} \sum_s^{n_s} I^a_{r,y,a,s,sf} W_{r,y,a,s} \end{aligned}$$

where $q^{Srv}_{r,y,sf}$ is the survey catchability coefficient. Currently, survey catchability coefficients can only vary over time as time blocks.

Tagging Observation Model

Tag cohorts ($T^k_{r,y,a,s}$) are dimensioned by regions, years, ages, and sexes, where index $k$ denotes a unique tag cohort (release region and release year combination). Tag cohorts are tracked for a user-specific maximum duration (maximum tag liberty; $n_L$), after which calculations for the tag cohort are no longer computed. This is done to help with computation time. Additionally, tag cohorts also have the option to undergo a fraction of mortality in the year of release, such that if a given cohort is released during the middle of the year, mortality processes are discounted. If mortality discounting is specified, movement does not occur in the year of release. In general, the process dynamics for the tagged cohort mimics those specified for the overall population. Immediately following release, tag cohorts are decremented by an initial tag-induced mortality rate:

$$T^k_{r,y,a,s} = T^k_{r,y,a,s} \exp(-\tau)$$

where $\tau$ is the initial tag-induced mortality rate, which can either be estimated freely or fixed at a user-specified value. In practice, this parameter is often cofounded with other mortality processes and needs to be fixed. If tagged cohorts are released at the beginning of a calendar year, Markovian movement occurs (movement does not occur otherwise):

$$\boldsymbol{T}^k_{y,a,s} = (\boldsymbol{T}^k_{y,a,s})^T \boldsymbol{M}_{y,a,s}$$

Mortality and ageing of the tagged cohort then occurs:

$$\begin{aligned} T^k_{r,y+1,a+1,s} = T^k_{r,y,a,s}\exp(-Z^{Tag}_{r,y,a,s}), \quad \text{for } 1 \leq a < a_+ \\ T^k_{r,y+1,a,s} =T^k_{r,y,a-1,s}\exp(-Z^{Tag}_{r,y,a-1,s}) + T^k_{r,y,a,s}\exp(-Z^{Tag}_{r,y,a,s}), \quad \text{for } a = a_+ \end{aligned}$$

where tagged cohorts follow an exponential mortality model, with accumulation of individuals in the plus-group. Total mortality for the tagged cohort can be specified as:

$$Z^{Tag}_{r,y,a,s} = (\kappa + NatMort^{Tag}_{r,y,a,s} + FishMort^{Tag}_{r,y,a,s})$$

$\kappa$ is a parameter describing chronic tag loss, which can be estimated in the model. However, this parameter should be fixed in practice, given confounding with mortality processes. Natural mortality for tagged individuals can be derived in several ways:

1. Natural mortality applied to tagged cohorts is averaged across ages and sexes,
2. Natural mortality applied to tagged cohorts is averaged across sexes, but is age-specific,
3. Natural mortality applied to tagged cohorts is averaged across ages, but is sex-specific, and
4. Natural mortality is both age-and sex-specific:

$$\begin{aligned} NatMort^{Tag}_{r,y,a,s} = \frac{\sum_a\sum_s NatMort_{r,y,a,s}}{n_an_s} \\ NatMort^{Tag}_{r,y,a,s} = \frac{\sum_s NatMort_{r,y,a,s}}{n_s} \\ NatMort^{Tag}_{r,y,a,s} = \frac{\sum_a NatMort_{r,y,a,s}}{n_a} \\ NatMort^{Tag}_{r,y,a,s} = NatMort_{r,y,a,s} \end{aligned}$$

The different calculations for natural mortality applied to tagged cohorts might be used if the observed tagging data are not age-or sex-specific such that calculations between observed and expected values are consistent with the nature of tag recapture data (e.g., if tagging data are not age-or sex-specific, one might use option 1). Similarly, fishing mortality experienced by tagged cohorts can also be specified in a multitude of ways:

1. Fishing mortality at age and sex from fleet 1 is applied to all ages and sexes with uniform selectivity,
2. Fishing mortality at age and sex from fleet 1 is multiplied by sex-averaged selectivity values,
3. Fishing mortality at age and sex from fleet 1 is multiplied by sex-specific selectivity values,
4. Fishing mortality at age and sex is derived by taking the weighted sum of annual fishing mortality across fleets with sex-averaged selectivity values, and
5. Fishing mortality at age and sex is derived by taking the weighted sum of annual fishing mortality across fleets with sex-specific selectivity values:

$$\begin{aligned} FishMort^{Tag}_{r,y,a,s} = Fmort_{r,y,f=1} \\ FishMort^{Tag}_{r,y,a,s} = Fmort_{r,y,f=1} \cdot \frac{\sum_s Sel_{r,y,a,s,f=1}^{Fsh}}{n_s} \\ FishMort^{Tag}_{r,y,a,s} = Fmort_{r,y,f=1} \cdot Sel_{r,y,a,s,f=1}^{Fsh} \\ FishMort^{Tag}_{r,y,a,s} = \sum_f^{n_f} Fmort_{r,y,f} \\ FishMort^{Tag}_{r,y,a,s} = \sum_{f=1} \left[Fmort_{r,y,f} \cdot \frac{\sum_s Sel_{r,y,a,s,f}^{Fsh}}{n_s} \right] \\ FishMort^{Tag}_{r,y,a,s} = \sum_{f=1} \left[Fmort_{r,y,f} \cdot Sel_{r,y,a,s,f}^{Fsh} \right] \end{aligned}$$

Again, users might specify different options depending on the nature of observed recapture data. For instance, if sex-specific information was not available, one might specify option 2 or 4 (sex-averaged selectivity values).

To derive the expected number of recaptures, a modified version of Baranov’s catch equation is used:

$$\text{Recap}^k_{r,y,a,s} = \beta_{r,y} \frac{FishMort^{Tag}_{r,y,a,s}}{Z^{Tag}_{r,y,a,s}} T^k_{r,y,a,s} \left[1-\exp(-Z^{Tag}_{r,y,a,s})\right]$$

where $\beta_{r,y}$ represents a tag reporting rate parameter that can vary across space and time. Currently, time-variation in $\beta_{r,y}$ can only be specified as a time-block. Additionally, $\beta_{r,y}$ is estimated as a logit variable:

$$\beta_{r,y} = \frac{\exp(\beta^{'}_{r,y})}{1 + \exp(\beta^{'}_{r,y})}$$

and $\beta^{'}_{r,y}$ represents the tag reporting parameter in logit space.

Observation Likelihoods

Currently, the generalized assessment model has the capability to fit to a variety of data sources, which include:

1. fishery catches,
2. fishery indices,
3. fishery age compositions,
4. fishery length compositions,
5. survey indices,
6. survey age compositions,
7. survey length compositions, and
8. tagging data

Each of these data sources are fit to using a statistical likelihood, where each likelihood component summed together to compute the joint negative log likelihood. Therefore, the joint negative log likelihood represents the function to be minimized.

Fishery Catches

Fishery catches can fit to using a lognormal likelihood, given by the following:

$$L(log(\text{ObsCatch}_{r,y,f})) = \lambda_{\text{ObsCatch}} \frac{1}{\sqrt{2\pi} \sigma^{\text{ObsCatch}}} \exp(-\frac{[log(\text{ObsCatch}_{r,y,f}) - log(\text{Catch}_{r,y,f})]^2}{2\sigma^{2^{\text{ObsCatch}}}})$$

where $\text{ObsCatch}_{r,y,f}$ represent the observed catches and $\sigma^{2^{\text{ObsCatch}}}$ is the standard deviation of the catch. In general, $\sigma^{2^{\text{ObsCatch}}}$ is not estimable and should be fixed at a small value (e.g., 1e-03) such that catch is assumed to be known, which is the default in SPoRC. Users should be able to easily modify this value in the parameter list to something other than 1e-03 if desired. $\lambda_{\text{ObsCatch}}$ in this case represents a likelihood weight to apply to catch data (generally not recommended to use and should be set to 1). In a case where only spatially-aggregated catches are available for a period of time, the model has the option to switch to an aggregated likelihood for those defined periods:

$$L(log(\text{ObsCatch}_{y,f})) = \lambda_{\text{ObsCatch}} \frac{1}{\sqrt{2\pi} \sigma^{\text{ObsCatch}}} \exp(-\frac{[log(\text{ObsCatch}_{y,f}) - log(\sum_r \text{Catch}_{r,y,f})]^2}{2\sigma^{2^{\text{ObsCatch}}}})$$

where a lognormal likelihood is still assumed, but the expected catches are summed across regions to derive spatially-aggregated catch.

Fishery and Survey Indices

Fishery indices can be fit to assuming a lognormal likelihood:

$$L(log(\text{ObsFshIdx}_{r,y,f})) = \lambda_{\text{ObsFshIdx}} \frac{1}{\sqrt{2\pi} \sigma^{\text{ObsFshIdx}}} \exp(-\frac{[log(\text{ObsFshIdx}_{r,y,f}) - log(\text{FshIdx}_{r,y,f})]^2}{2\sigma^{2^{\text{ObsFshIdx}}}})$$

where $\text{ObsFshIdx}_{r,y,f}$ are the observed fishery indices and $\sigma^{2^{\text{ObsFshIdx}}}$ describes the variance of the fishery index (can be estimated but should theoretically be known). $\lambda_{\text{ObsFshIdx}}$ represents a likelihood weight applied to fishery indices (should ideally not be used and set at a value of 1 to allow the model to objectively weight data sources using the specified variance).

Similarly, survey indices are fit to assuming a lognormal likelihod:

$$L(log(\text{ObsSrvIdx}_{r,y,sf})) = \lambda_{\text{ObsSrvIdx}} \frac{1}{\sqrt{2\pi} \sigma^{\text{ObsSrvIdx}}} \exp(-\frac{[log(\text{ObsSrvIdx}_{r,y,sf}) - log(\text{SrvIdx}_{r,y,sf})]^2}{2\sigma^{2^{\text{ObsSrvIdx}}}})$$

where $\text{ObsSrvIdx}_{r,y,sf}$ are the observed survey indices and $\sigma^{2^{\text{ObsSrvIdx}}}$ describes the variance of the survey index (can be estimated but should theoretically be known). $\lambda_{\text{ObsSrvIdx}}$ represents a likelihood weight applied to survey indices (should ideally not be used and set at a value of 1 to allow the model to objectively weight data sources using the specified variance).

Fishery and Survey Compositions

At present, a total of 3 multivariate likelihoods are available for use to fit to composition data. These include the multinomial, the Dirichlet-multinomial, and logistic-normal likelihoods.

In the case of a multinomial likelihood, the following expression is invoked: $$L((\text{ObsCompositionData}_{r,y,f}) = \lambda_{\text{ObsCompositionData}} \cdot ISS \cdot \prod_{b = 1}^{n_b} E_b^{O_b}$$ where $ISS$ is the input sample size that determines the initial weight to apply to a given composition dataset, and $\lambda_{\text{ObsCompositionData}}$ is the weight to apply to the input sample size after some form of iterative re-weighting process (e.g., Francis). $E_b$ describes the expected values for a given bin $b$, while $O_b$ describes the observed values.

If a Dirichlet-multinomial likelihood (linear parameterization) is assumed, the following expression is used: $$L(\text{ObsCompositionData}_{r,y,f}) = \lambda_{\text{ObsCompositionData}} \cdot \frac{\Gamma(ISS + 1)}{\prod_b ISS \cdot O_b + 1} \frac{\Gamma(\theta \cdot ISS)}{\Gamma(ISS + \theta \cdot ISS)} \prod_{b = 1}^{n_b} \frac{\Gamma(ISS \cdot O_b + \theta \cdot ISS \cdot E_b)}{\theta \cdot ISS \cdot E_b}$$

where $\theta$ here represents the overdispersion parameter such that the input sample size can be adjusted to obtain an effective sample size ($ESS$):

$$ESS = \frac{1}{1 + \theta} + ISS \frac{\theta}{1 + \theta}$$

Thus, if $\theta$ was large such that ($\theta >> ISS$), then $ESS -> ISS$, while if $\theta$ was small such that ($\theta << ISS$) then $\theta$ would approximate the ratio of $ESS$ and $ISS$.

Lastly, a logistic-normal likelihood can also be invoked, where several covariance structures can be parameterized. In it’s simplest form, an iid covariance can be specified:

$$\boldsymbol{\Sigma} = \boldsymbol{I} \cdot \theta^2$$

where $\boldsymbol{\Sigma}$ is the covariance structure scaled by $\theta^2$ to allow for overdispersion and inter-annual variation in weighting. $\boldsymbol{I}$ is an identity matrix dimensioned by $n_b - 1 \cdot n_b - 1$. A first-order autoregressive process can also be specified across ages:

$$\boldsymbol{\Sigma} = \boldsymbol{R_a} \cdot \theta^2$$

where $\boldsymbol{R_a}$ is a matrix representing a first-order autoregressive process. Two-dimensional autoregression can also be specified across ages and sexes:

$$\boldsymbol{\Sigma} = \boldsymbol{R_a} \otimes \boldsymbol{R_s} \cdot \theta^2$$

where $\boldsymbol{R_s}$ in this case is a constant correlation matrix by sex. Following the specification of a covariance matrix, a logistic transformation is then conducted on the observed and expected compositions, and a multivariate normal likelihood is invoked on logistic random variables:

$$\begin{aligned} O^{'}_b = log(O_{-n_b}) - log(O_{n_b}) \\ E^{'}_b = log(E_{-n_b}) - log(E_{n_b}) \\ L(\text{ObsCompositionData}_{r,y,f}) = \lambda_{\text{ObsCompositionData}} \cdot 2\pi^{-k/2} \text{det}(\boldsymbol{\Sigma})^{-1/2} \exp\left[-\frac{1}{2} (O^{'}_b - E^{'}_b)^T \boldsymbol{\Sigma}^{-1}(O^{'}_b - E^{'}_b)\right] \end{aligned}$$

where $O^{'}_b$ and $E^{'}_b$ are the observed and expected compositions after a logistic additive transformation.

Structuring Observed and Expected Compositions

Four options can be used for structuring age and length compositions. These are:

1. Aggregated age or length compositions across regions and sexes,
2. ‘Split’ age or length compositions for each region and sex (i.e., no implicit information about region and sex-ratios),
3. ‘Split’ age or length compositions for each region, coupled with ‘Joint’ compositions across sexes (i.e., implicit information about sex-ratios, but no information about regional-ratios), and

In the first option (termed ‘Aggregated’), both observed and expected age or length compositions within a given year and fleet (survey or fishery) are aggregated across both sexes and regions. For simplicity, observed age and length compositions will be denoted as $O_{r,b,s}$ and expected age and length compositions will be denoted as $E_{r,b,s}$. Here, $b$ represents either age or length bins. Note that observed compositions will need to be aggregated prior to the modelling process (but need not be normalized, as it is normalized within the modeling process). Expected compositions are then first normalized by dividing the sum of a given region and/or sex such that proportions sum to 1 within a given region and a given sex. Following that, the average of $E_{r,b,s}$ across regions and sexes is taken. If age compositions have ageing error associated with them, the expected age compositions are then multiplied by an ageing-error matrix to derive the compositions to fit to. Note that a non-square (e.g., rectangular) ageing-error matrix can be provided to reduce the number of modelled ages to align with the number of observed composition ages. This approach is common for many rockfish species in Alaska, where the model includes a greater number of ages than are represented in the observed compositions. This allows for a more accurate representation of population dynamics (by modelling more ages), even though age estimates for older fish are often unreliable. This process is described by the equations below:

$$\begin{aligned} O_b = \frac{O_b}{\sum_b O_b} \\ E_b = \frac{1}{n_s \cdot n_r} \sum_{i = 1}^{n_s \cdot n_r} \frac{E_{r,b,s}}{\sum_b E_{r,b,s}} \\ \boldsymbol{E} = \boldsymbol{E} \textbf{AgeError} \\ \boldsymbol{O} = \left\{ O_b \;\middle|\; b = 1, \ldots, n_b \right\} \\ \boldsymbol{E} = \left\{ E_b \;\middle|\; b = 1, \ldots, n_b \right\} \end{aligned}$$

In the second option (termed ‘SplitRegionSplitSex’), both observed and expected compositions for a given year and fleet are normalized to sum to 1 within a given sex and region. For age compositions, an ageing-error matrix can then applied to a given region and sex. This process is mathematically described in the following equation:

$$\begin{aligned} O_{r,b,s} = \frac{O_{r,b,s}}{\sum_b O_{r,b,s}} \\ E_{r,b,s} = \frac{E_{r,b,s}}{\sum_b E_{r,b,s}} \\ \boldsymbol{E}_{r,s} = \boldsymbol{E}_{r,s} \textbf{AgeError} \\ \boldsymbol{O}_{r,s} = \left\{ O_{r,b,s} \;\middle|\; b = 1, \ldots, n_b \right\} \\ \boldsymbol{E}_{r,s} = \left\{ E_{r,b,s} \;\middle|\; b = 1, \ldots, n_b \right\} \end{aligned}$$

In the third option (termed ‘SplitRegionJointSex’), observed and expected compositions for a given year and fleet are normalized to sum to 1 within a given region, but jointly across sexes. For age compositions, the ageing-error matrix can then extended to accommodate such that it is dimensioned by $n_bn_s \cdot n_bn_s$, and is given by the following equations:

$$\begin{aligned} O_{r,b,s} = \frac{O_{r,b,s}}{\sum_{b,s} O_{r,b,s}} \\ E_{r,b,s} = \frac{E_{r,b,s}}{\sum_{b,s} E_{r,b,s}} \\ \boldsymbol{E}_{r} = (\boldsymbol{E}_r)^T (\textbf{I} \otimes \textbf{AgeError}) \\ \boldsymbol{O}_r = \left\{ O_{r,b,s} \;\middle|\; b = 1,\ldots,n_b;\; s = 1,2 \right\} \\ \boldsymbol{E}_r = \left\{ E_{r,b,s} \;\middle|\; b = 1,\ldots,n_b;\; s = 1,2 \right\} \end{aligned}$$

Tagging

The generalized assessment model has the option to fit to tagging data, which tracks a tag cohort across regions, ages, and sexes. However, in many cases, age or sex information may not be available for a given species. Thus, when tagging data does not include age or sex information, the user must specify all tagged individuals to occupy the first age class and first sex category as a starting point. Users should then also specify tagging data to be fit via pooling across all ages and sexes using the move_age_tag_pool and move_sex_tag_pool data inputs under the scenario that no information is available to distinguish age-and sex-specific tag recaptures.

Given that tagging data is highly multi-dimensional, use of these data can greatly increase the computational demand of the model. Therefore, to partially reduce these demands, the model allows for pooling of data across age and sex partitions. This process groups observations and expected values into their respective age blocks. It is recommended that if age-varying movement is estimated in age-blocks, pooling of tag data should align with these age-blocks. Moreover, it is important to note that when tag recoveries are fit to by pooling observations and expected values, likelihood values for alternative pooling parameterizations are not directly comparable.

Additionally, a central tenant of using tagging data is that tagged individuals are representative of the overall population such that these individuals need to be well-mixed. To address the well-mixed assumption, the model allows for a specification of mixing period, which specifies the time-at-liberty for which tags are fit to. For example, if all individuals are well-mixed followed release, then the mixing period would be specified at a value of 1. By contrast, if individuals are only well-mixed after 2 years, then the mixing period would be specified at a value of 2. Therefore, the model requires the entire history of a given tag cohort to be input (starting from the initial release and recapture year).

Currently, several tag likelihoods can be utilized to describe tag recapture data. These options are: 1) Poisson likelihood, 2) Negative Binomial likelihood, 3) Multinomial likelihood, assuming release conditioned dynamics, and 4) Multinomial likelihood, assuming recapture conditioned dynamics

The Poisson likelihood for tag recaptures is given by the following expression: $$L(\text{ObsRecap}_{r,y,a,s}) = \frac{\exp(-\text{Recap}^k_{r,y,a,s})(\text{Recap}^k_{r,y,a,s})^{\text{ObsRecap}_{r,y,a,s}}}{\text{ObsRecap}_{r,y,a,s}!}$$

where $\text{Recap}^k_{r,y,a,s}$ are observed tag recaptures. If pooling is specified, then observed recoveries are fit to but summing the observed and expected tag recoveries across the respective age and sex groups. In the model, the likelihood is specified in log space, such that the factorial term in the denominator is specified as a gamma function, allowing for recoveries to be fit to as non-integer terms (as would be the case if an age-length key were used to convert observed lengths to ages).

The negative binomial likelihood for tag recaptures is given by the expression: $$L(\text{ObsRecap}^k_{r,y,a,s}) = \frac{\Gamma(\text{ObsRecap}_{r,y,a,s} + \eta)}{\Gamma(\eta)\Gamma(\text{ObsRecap}_{r,y,a,s} + 1)} \left(\frac{\text{Recap}_{r,y,a,s}}{\text{Recap}_{r,y,a,s} + \eta}\right)^{\text{ObsRecap}_{r,y,a,s}} \left(\frac{\eta}{\text{Recap}_{r,y,a,s}^k + \eta}\right)^\eta$$

where $\eta$ represents the overdispersion parameter for tagging data.

The multinomial release conditioned likelihood for tag recoveries fits to both recaptured and non-recaptured states. Proportions of observed ($\text{PObsRecap}^k_{r,y,a,s}$) and expected recaptured ($\text{PRecap}^k_{r,y,a,s}$) individuals are given by:

$$\begin{aligned} \text{PRecap}^k_{r,y,a,s} = \frac{\text{Recap}^k_{r,y,a,s}}{\text{InitTag}^k} \\ \text{PObsRecap}^k_{r,y,a,s} = \frac{\text{ObsRecap}^k_{r,y,a,s}}{\text{InitTag}^K} \end{aligned}$$

$\text{InitTag}^k$ are the total tags released for a given tag cohort for a given release region and year. Non-recaptured states can then be written as:

$$\begin{aligned} \text{PNonRecap}^k = 1 - \sum_r \sum_y \sum_a \sum_s \text{PRecap}^k_{r,y,a,s} \\ \text{PObsNonRecap}^k = 1 - \sum_r \sum_y \sum_a \sum_s \text{PObsNonRecap}^k_{r,y,a,s} \end{aligned}$$

where $\text{PNonRecap}^k$ and $\text{PObsNonRecap}^k$ are the expected and observed non-recapture proportions, respectively. Rewriting both recaptured and non-recaptured states for observed and expected values gives:

$$\begin{aligned} \boldsymbol{E}^k_{\text{Tagging}} = \{\textbf{PRecap}^k, \text{PNonRecap}^k \} \\ \boldsymbol{O}^k_{\text{Tagging}} = \{\textbf{PObsRecap}^k, \text{PObsNonRecap}^k \} \end{aligned}$$

Following that, a standard multinomial likelihood with the aforementioned vectors can be written as:

$$L(\text{ObsRecap}^k_{r,y,a,s}) = \text{InitTag}^k \prod_i (E^k_{i,\text{Tagging}})^{O^k_{i,\text{Tagging}}}$$

where $i$ indexes a given element in $E^k_{i,\text{Tagging}}$ and $O^k_{i,\text{Tagging}}$.

The multinomial recapture conditioned likelihood for tag recoveries assumes that tag shedding, tag induced mortality, and tag reporting rates are spatially-invariant, given that these terms cancel out in the denominator. Thus, these terms need not be specified or estimated when assuming this likelihood. The following expression is used to transform observed and expected recaptured individuals to proportions:

$$\begin{aligned} \text{PRecap}^k_{r,y,a,s} = \frac{\text{Recap}^k_{r,y,a,s}}{\sum_r\sum_a\sum_s \text{Recap}^k_{r,y,a,s}} \\ \text{PObsRecap}^k_{r,y,a,s} = \frac{\text{ObsRecap}^k_{r,y,a,s}}{\sum_r\sum_a\sum_s \text{ObsRecap}^k_{r,y,a,s}} \\ \text{TotRecap}_y^k = \sum_r\sum_a\sum_s \text{ObsRecap}^k_{r,y,a,s} \end{aligned}$$

where recapture probabilities are normalized by the sum of all recaptures in a given year and $\text{TotRecap}_y^k$ are the total number of observed recaptures within a given year. The multinomial likelihood is then expressed as:

$$L(\text{ObsRecap}^k_{r,y,a,s}) = \text{TotRecap}_y^k \prod_i (\text{PRecap}^k_{r,y,a,s})^{\text{PObsRecap}^k_{r,y,a,s}}$$

Priors

Several priors are available to help inform the estimation of parameters. In particular, priors are available for natural mortality, catchability, steepness, movement, and tag reporting rates.

Natural Mortality

For natural mortality, a lognormal prior can be invoked:

$$P(log(NatMort_{s = 1}) = \frac{1}{\sqrt{2\pi} \sigma^{\text{NatMortPrior}}} \exp(-\frac{[log(NatMort_{s = 1}) - log(Natmort_{s = 1}^{\text{Prior}})]^2}{2\sigma^{2^{\text{NatMortPrior}}}})$$

where the lognormal prior is only invoked on the first sex (female if sex-structured), $\sigma^{2^{\text{NatMortPrior}}}$ represents the variance of the prior, and $Natmort_{s = 1}^{\text{Prior}}$ represents the prior value itself.

Fishery and Survey Catchability

For fishery or survey catchability, a lognormal prior can similarly be used:

$$P(log(q_{r,y,f})) = \frac{1}{\sqrt{2\pi} \sigma_{r,y,f}^{\text{qPrior}}} \exp(-\frac{[log(q_{r,y,f}) - log(q_{r,y,f}^{\text{Prior}})]^2}{2\sigma^{2^{\text{qPrior}}}_{r,y,f}})$$

where $\sigma^{2^{\text{qPrior}}}$ represents the variance of the prior, and $q_{r,y,f}^{\text{qPrior}}$ represents the prior value itself.

Steepness

If a Beverton-Holt stock recruitment relationship is assumed, priors for steepness can be imposed. Currently, a scaled beta prior (bounded between 0.2 and 1) can be invoked:

$$\begin{aligned} a = \left(\frac{h^{\text{PriorMu}}_r - 0.2}{1 - 0.2}\right) \left(\frac{h^{\text{PriorSD}}_r }{1 - 0.2}\right)^2 \\ b = \left[1 - \left(\frac{h^{\text{PriorMu}}_r - 0.2}{1 - 0.2}\right)\right] \left[\left(\frac{h^{\text{PriorSD}}_r }{1 - 0.2}\right)^2\right] \\ P(h_r) = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}h^{a-1}_r(1-h_r)^{b-1} \end{aligned}$$

where $h^{\text{PriorMu}}_r$ represents the prior mean steepness value in a given region (bounded between 0.2 and 1) and $h^{\text{PriorSD}}_r$ is the standard deviation for those priors. $a$ and $b$ are the shape parameters for a beta distribution.

Movement

Priors on movement parameters can also be invoked, where these priors are imposed on the rows of the movement matrix itself (as opposed to the multinomial logit parameters). Here a Dirichlet prior can be invoked:

$$P(M_{k,y,a,s}) = \frac{\Gamma(\sum_{j=1}^{n_r} \alpha_{j,y,a,s})}{\prod_{j=1}^{n_r}\Gamma(\alpha_{j,y,a,s})} \prod_{j=1}^{n_r}M_{k,j,y,a,s}^{\alpha_{j,y,a,s}-1}$$

where $j$ indexes movement towards a given region and $\alpha_{j,y,a,s}$ are the Dirichlet prior parameters. These priors have the flexibility to be specified for movement of individuals towards region j for any year, age, or sex partition. If a user were to specify a uniform prior, all values in $\alpha_{j,y,a,s}$ would have the same values. Increasing the values of $\alpha_{j,y,a,s}$ increases the concentration (i.e., the weight of the prior).

Tag Reporting Rates

Two types of priors can be specified for facilitating the estimation of tag reporting rates. In particular, the first type of prior includes a symmetric beta parameterization, which is given by the following:

$$P(\beta_{r,y}) = (\beta_{r,y} + 1e-4)^{\beta^{\text{PriorScale}} } (1 - \beta_{r,y} + 1e-4)^{\beta^{\text{PriorScale}}}$$

where $\beta^{\text{PriorScale}}$ determines the scale of the tag reporting parameter and controls how strongly to penalize estimates when they approach the bounds of $(0,1)$. Larger values of $\beta^{\text{PriorScale}}$ result in larger penalties. By contrast, a standard beta distribution can also be specified, with shape parameters $a$ and $b$. This is given by the following expression:

$$\begin{aligned} a = \frac{\beta^{\text{PriorMu}}}{\left(\beta^{\text{PriorSD}}\right)^2} \\ b = \left[\frac{1 - \beta^{\text{PriorMu}}}{\left(\beta^{\text{PriorSD}}\right)^2}\right] \\ P(\beta_{r,y}) = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)} \, \beta_{r,y}^{a-1} \left(1 - \beta_{r,y}\right)^{b-1} \end{aligned}$$

where $\beta^{\text{PriorMu}}$ defines the mean of the distribution, while $\beta^{\text{PriorSD}}$ defines the standard deviation of the beta distribution. Currently, only a tag reporting rate priors are not region or year-specific (i.e., the same value is applied across all regions and years).

Penalties (Process Error)

Initial Age Deviations

Estimation of initial age deviations is constrained by a normal distribution:

$$\epsilon_{r,y}^{Init} \sim N\left(0,\sigma^2_{Init}\right)$$

with a mean of 0 and a variance of $\sigma^2_{Init}$. In SPoRC, these can either be estimated in a penalized likelihood framework, where the variance parameter is fixed at a known value and $\epsilon_{r,y}^{Init}$ are estimated as fixed-effects. Conversely, a random effects framework can also be used where $\epsilon_{r,y}^{Init}$ are treated as random effects and are integrated out using Laplace approximation.

Recruitment

Similarly, estimation of recruitment deviations is constrained by a normal distribution:

$$\epsilon_{r,y}^{Rec} \sim N\left(0,\sigma^2_{Rec}\right)$$

with a mean of 0 and a variance of $\sigma^2_{Rec}$. Again, recruitment deviations can either be estimated in a penalized likelihood framework, where the variance parameter is fixed at a known value and $\epsilon_{r,y}^{Rec}$ are estimated as fixed-effects. Conversely, $\epsilon_{r,y}^{Rec}$ can also be treated as random effects and integrated out using Laplace approximation.

Fishing Mortality

Fishing mortality deviations are also constrained by a normal distribution:

$$\epsilon_{r,y,f}^{Fmort} \sim N\left(0,\sigma^2_{Fmort}\right)$$

where $\sigma^2_{Fmort}$ is fixed at an arbitrary value to ensure that observed catches are fit to adequately. Given that values of observed catch are generally known with negligible error, the value specified for $\sigma^2_{Fmort}$ should be inconsequential. It is generally necessary to constrain fishing mortality deviations to ensure that all parameters are identifiable. Currently, $\sigma^2_{Fmort}$ can be specified as either random effects or penalized likelihood. Under a penalized likelihood framework, this value is fixed at 1, although this can be easily modified by the user.

Selectivity (Fishery and Survey)

Several options can be specified to parameterize continuous time-varying selectivity. For all forms, users can specify deviations to occur in a penalized likelihood framework, or in a random effects framework. The simplest form represents iid deviations about selectivity parameter $i$ for a given functional form: $$\epsilon_{r,y,i,s,f}^{Sel} \sim N\left(0,\sigma^2_{i,Sel}\right)$$ where $\sigma^2_{i,Sel}$ is the variance of selectivity parameter $i$ (e.g., $a^{50}$). Using this form, selectivity is assumed to vary across time but in a parametric manner.

Parametric selectivity deviations can also be specified assuming a random walk process: $$\begin{aligned} \epsilon_{r,y,i,s,f}^{Sel} \begin{cases} \sim N\left(0,25\right), & \text{if } y = 1 \\ \sim N\left(\epsilon_{r,y-1,i,s,f}^{Sel},\sigma^2_{i,Sel}\right) & \text{if } y > 1 \end{cases} \end{aligned}$$

where the first parameter deviation is initialized with a broad distribution (essentially a fixed effect parameter), and subsequent deviations are dependent on the previous time point, constrained by an estimable variance for a given selectivity parameter ($\sigma^2_{i,Sel}$). For both iid and random walk deviations, a smoothness penalty is applied onto selectivity values to induce temporal stability:

$$P_{\text{year smoothness}} = \sum_{s=1}^S\sum_{y=2}^{Y} \sum_{b=1}^{B} (log(Sel_{y,b,s}) - log(Sel_{y-1,b,s}))^2$$

Semi-parametric deviations can also be specified, where deviations can occur across years, bins (age or length), and cohorts. For all semi-parametric deviations, a curvature penalty is applied on selectivity values across bins and across years for regularity:

$$\begin{aligned} P_{\text{bin curvature}} = \sum_{s=1}^S\sum_{y=1}^{Y} \sum_{b=2}^{b-1} (log(Sel_{y,b+1,s}) - 2log(Sel_{y,b,s}) - log(Sel_{y,b-1,s}))^2 \\ P_{\text{year curvature}} = \sum_{s=1}^S\sum_{y=2}^{Y-1} \sum_{b=1}^{B} (log(Sel_{y+1,b,s}) - 2log(Sel_{y,b,s}) - log(Sel_{y-1,b,s}))^2 \\ \end{aligned}$$

Two options can be specified to allow for deviations across the dimensions of ages, years, and cohorts (and note that this is only possible when age-based selectivity is specified), which include a marginal stationary and a conditional non-stationary variance option. The primary difference between the two is that the marginal variance version does not have a closed form solution for $\boldsymbol{\Omega}$ (see equations below) and needs to be iteratively solved, while the conditional variance version has a closed form solution. The following equations describe the conditional variance version:

$$\boldsymbol{\epsilon}_{r,s,f}^{Sel} \sim \text{MVN}(0,\textbf{Q})$$

where process error arises from a multivariate normal distribution with precision matrix ($Q$; inverse of the covariance). Here, $Q$ is constructed as:

$$\begin{aligned} \textbf{Q} = (\textbf{I} - (\textbf{B})^T) \boldsymbol{\Omega} (\textbf{I} - \textbf{B}) \\ \text{diag}(\boldsymbol{\Omega}) = \sigma_{Sel}^{-2} \end{aligned}$$

$\textbf{I}$ is an identity matrix and $\boldsymbol{\Omega}$ is a diagonal matrix that determines the variance of the multivariate normal process. $\textbf{B}$ is a square matrix representing the partial effect of $\boldsymbol{\epsilon}_{r,s,f}^{Sel}$ on preceeding ages and/or years, governing by partial correlation coefficients for ages, years, and cohorts. To demonstrate the formulation of $\textbf{B}$, a simplified example is provided with rows representing ages $a \epsilon \{1,2\}$ and columns representing years $t \epsilon \{1,2\}$. In this example, $\textbf{B}$ is a $4 x 4$ matrix, where both the rows and columns represent combinations of age and year. For instance, $B_{1,2}$ captures the correlation within the same year between age 1 in year 1, and age 2 in year 1, whereas $B_{4,1}$ constructs the correlation within the same cohort between age 1 in year 1, and age 2 in year 2:

$$\begin{aligned} \mathbf{B} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ \rho_y & 0 & 0 & 0 \\ \rho_a & 0 & 0 & 0 \\ \rho_c & \rho_a & \rho_y & 0 \end{bmatrix} \end{aligned}$$

Here, $\rho_y$, $\rho_a$, and $\rho_c$ are parameters describing the partial autocorrelation among years within a given age, among ages within a given year, and years within a cohort, respectively.

Lastly, semi-parametric deviations can also be specified to vary across years and bins (i.e., in a two-dimensional manner; this option can be utilized for both age-and length-based selectivity). Similar to the process described above, selectivity deviations are constrained by a precision matrix $Q$, following a multivarite normal process:

$$\begin{aligned} \boldsymbol{\epsilon}_{r,s,f}^{Sel} \sim \text{MVN}(0,\textbf{Q}^{-1}) \\ \textbf{Q}^{-1} = \frac{\sigma^2_{Sel}}{(1-\rho_y)^2(1-\rho_b)^2}\textbf{R}_y \otimes \textbf{R}_b \\ \end{aligned}$$

where is determined by the Kronecker product of a first-order autoregressive process across bins $\textbf{R}_b$ and across years $\textbf{R}_y$.

Movement Deviations

Continuous movement deviations can also be estimated under a random effects or penalized likelihood framework. Currently, only iid deviations are allowed, where deviations are placed upon fixed-effects movement parameters in logit space:

$$\begin{aligned} \omega_{k,-1,y,a,s} = \omega_{k,-1,s} + \epsilon^{\text{Move}}_{k,j,y,a,s} \\ \epsilon^{\text{Move}}_{k,-1,y,a,s} \sim N(0,\sigma^2_{k,a,s}) \end{aligned}$$

where $\epsilon^{\text{Move}}_{k,-1,y,a,s}$ are movement deviatons in logit space. Because of the sum to 1 constraint, movement deviations are only estimated for $n_r - 1$. Deviations then arise from a normal distribution with process variation constrained by $\sigma^2_{k,a,s}$. Thus, process variance can be distinct by the origin region and age (although other parameterizations that share all parameters across origin region, ages, and sexes are also allowed). Similarly, process deviations can be invoked only across years, ages, years and ages, years, ages, and sexes, etc. This flexibility is maintained given the computational demand for estimating random effects for movement.