---
title: "8. Model parameterizations"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{8. Model parameterizations}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

# In progress!

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```
# Population dynamics
CEATTLE is an age-structured model that follows cohorts of individuals of a species $i$, age $a$, and sex $s$ through time $y$. Annual recruitment and the initial age-structure are estimated and are treated as random effects or as penalized deviates, with the variance in recruitment estimated or fixed, respectively. Note that because the initial age-structure is estimated (assuming the population is depleted prior to the first year of the model), the mean recruitment parameter ($\bar{R}_{i}$ will be biased low and projected recruitment will have to account for that (e.g. resemble recruitment deviates from hindcast). Similarly, biomass reference points should use the mean of annual recruitment. Populations can be modelled as 1- or 2-sexes via `nsex` in the data. The sex-ratio at birth $\rho_{i1}$ is set in the `sex_ratio` data for the first age. Natural mortality is partitioned into predation mortality (M2) and non-predation related mortality (M1), the latter of which is either input or estimated (`estM1`), and accounts for mortality related to causes other than predation by the species included in the model. Predation can be turned off by setting `msmMode = 0` (see below for more detail).

$$N_{is1,y} = R_{isy} = \bar{R}_{i }e^{\tau_{i,y}}*\rho_{i1}$$
$$N_{isa,1} = \left \{ \begin{array}{} \bar{R}_ie^{ \left( -j M1_{isa} \right)}*\rho_{i1} \\ \frac{\bar{R}_ie^{ \left( -\sum^{a-1} M1_{isa} \right)}} {\left(1 - e^{\left( -\sum^{a-1} M1_{is A_i} \right)} \right)} *\rho_{i1} \end{array} \right.$$
$$N_{is a+1, y+1} = N_{isa,y} e^{-Z_{isa,y}}$$
$$N_{is A_i, y+1} = N_{is A_{i}-1,y} e^{-Z_{is A_{i}-1,y}} + N_{is A_{i},y} e^{-Z_{is A_{i},y}}$$
$$B_{isa,y} = N_{isa,y} W_{isa,y}$$
$$SSB_{isa,y} = B_{isa,y} \rho_{ia}$$
$$Z_{isa,y} = M1_{isa} + M2_{isa,y} + F_{isa,y}$$

# Predation sub-model
Multiple age-based predation parameterizations are supported by Rceattle and are controlled by `suitMode` defining the predator-prey preference (suitability/predator selectivity) and `msmMode` defining the functional response of predation. The data used to parameterize the predation model include the proportion of prey-at-age in the diet of a predator-at-age (`UobsAgeWt` in the data) and the ration/consumption data (defined below).

## Suitability
Suitability represents a predator’s prey and size preference and factors such as spatial overlap in the distribution of predators and prey. 

### MSVPA Type 2 based suitability
Suitability is derived empirically from the ratio of mean proportion (by weight) of prey-at-age in the stomach of predators-at-age $\bar{U}_{pbjisa}$ (`UobsAgeWt` in the data). Magnusson (1995) provides details on the derivation of suitability coefficients. Magnusson (1995) provides details on the derivation of suitability coefficients.
`suitMode = 0`
Predator-prey suitability for predators of species $p$, sex $b$, age $j$, and prey of species $i$, sex $s$ and age $a$
$$\hat{S}_{pbjisa} = \frac{1}{n_y} \sum_y \left( \frac{ \frac{\bar{U}_{pbjisa}} {B_{isa,y}} } {\sum_{isa} \left( \frac{\bar{U}_{pbjisa}} {B_{isa,y}} \right) + \frac{1 + \sum_{isa} \bar{U}_{pbjisa}} {B_p^{other}} } \right)$$

## Predation mortality

### Single-species mode
This is essentially multiple single-species age structured models where all mortality is due to input or estimated natural mortality ($M = M1$) and fishing.
`msmMode = 0`
$$M2_{isa,y} = 0$$

### MSVPA Type 2 based predation mortality
`msmMode = 1`
Predation mortality can be modelled following previous MSCAA models (Curti et al., 2013; Holsman et al., 2016; Tsehaye et al., 2014) based on MSVPA. Predation mortality due to a predator-at-age in the model is derived from estimated or empirically derived suitability coefficients $\hat{S}_{pbjisa}$, input estimates of annual ration/consumption (see below), estimated or input biomass, and input biomass of other prey in the system for predator $p$ ($B_p^{other}$). Estimated average biomass is used to derive predation mortality and because of the interdependence between predation and biomass, predation mortality and the underlying age-structured dynamics are iteratively estimated until convergence. Preliminary analysis suggested that three to five iterations are sufficient for convergence. Models are generally insensitive to a range of $B_p^{other}$, but historically they have been derived from ecosystem-wide biomass estimates. Predation mortality equations for predators of species $p$, sex $b$, age $j$, and prey of species $i$, sex $s$ and age $a$
$$M2_{isa,y} = \sum_{pbj} \left( \frac{N_{pbj,y} \delta_{pbj,y} {S}_{pbjisa}} {\sum_{isa} \left({S}_{pbjisa} B_{isa,y} \right) + B_p^{other} \left(1 - \sum_{isa} \left({S}_{pbjisa} \right) \right)} \right)$$

### MSVPA Type 3 based predation mortality
`msmMode = 2`
A Type III functional response variant of the MSVPA predation model. The available-food denominator replaces prey biomass $B_{isa,y} = N_{isa,y} W_{isa,y}$ with $N^2_{isa,y} W_{isa,y}$, creating a sigmoidal (accelerating-then-saturating) relationship between prey density and predation rate. This gives low-density prey a refuge from predation and can stabilise multi-species dynamics at low stock sizes.
$$M2_{isa,y} = \sum_{pbj} \left( \frac{N_{pbj,y} \delta_{pbj,y} {S}_{pbjisa}} {\sum_{i'a's'} \left({S}_{pbji'a's'} N^2_{i'a's',y} W_{i'a's',y} \right) + B_p^{other} \left(1 - \sum_{i'a's'} \hat{S}_{pbji'a's'} \right)} \right)$$


## Consumption
Multiple age-based consumption parameterizations are supported by Rceattle. They include a environmentally driven bioenergetics consumption and direction input of consumption-at-age. The consumption model is controlled by the following objects in the data (`bioenergetics_control` sheet of the excel data input):

* `Ceq`	Integer: switch for which bioenergetics equation to use for each species for ft to scale max consumption: 1 = Exponential (Stewart et al 1983), 2 = Temperature-dependence for warm-water species (Kitchell et al 1977; sensu Holsman et al 2015), 3 = temperature dependence for cool and cold-water species (Thornton and Lessem 1979); 4 = temperature function held constant at 1.0 — bypasses bioenergetics and uses `ration_data` as direct empirical annual ration input (set `CA = 1`, `fday = 1`, `CB = 0`)
* `Cindex`	Integer: which environmental index in env_data to use to use for $T_y$ for bioenergetics based ration
* `Pvalue`	Tuning parameter $\varphi_p$ that scales the maximum consumption used for ration for each species
* `fday` Data $Fday_p$ number of foraging days per year for each species
* `CA`	Parameter $\alpha^{\delta}_p$ for slope of allometric mass function for calculating maximum consumption: CA * Weight ^ CB
* `CB`	Parameter $\beta^{\delta}_p$ for slope of allometric mass function for calculating maximum consumption: CA * Weight ^ 1+CB
* `Qc`	Parameter $Q^c_p$ for temperature scaling function of maximum consumption specified by `Ceq`
* `Tco`	Parameter $T^{co}_p$ for temperature scaling function of maximum consumption specified by `Ceq`
* `Tcm`	Parameter $T^{cm}_p$ for temperature scaling function of maximum consumption specified by `Ceq`
* `Tcl`	Parameter $T^{cl}_p$ for temperature scaling function of maximum consumption specified by `Ceq`
* `CK1`	Parameter $CK1_p$ for temperature scaling function of maximum consumption specified by `Ceq`
* `CK4`	Parameter $CK4_p$ for temperature scaling function of maximum consumption specified by `Ceq`

`ration_data` $RFR_{pbj,y}$ is the relative foraging rate data used to modify ration. Essentially relative foraging rate and, subsequently, ration can be tuned via `fday` and `pvalue`.

The weight-at-age used for the consumption equations is specified by `pop_wt_index`

### Bioenergetics based ration
`Ceq %in% 1:3`
Individual specific ration ($\delta_{pbj,y}$ $(kg~kg^{-1}~yr^{-1})$) of predator $p$ of sex $b$ and age $j$ in year $y$:
$$\delta_{pbj,y} = RFR_{pbj,y}*Fday_p *\alpha^{\delta}_p W_{pbj,y}^{\left(1+\beta^{\delta}_p\right)} *f\left(T_y \right)_p$$
Ration can be tuned by adjusting `Pvalue` ($\varphi_p$) to help convergence:
$$\delta_{pbj,y} = \varphi_p * \delta_{pbj,y}$$

*Temperature scaling algorithms. *

`Ceq = 1`: Exponential function from Stewart et al. 1983:
$$f \left(T_y \right)_p = e^{Q^c_p *T_y}$$

`Ceq = 2`: Temperature dependence for warm-water-species from Kitchell et al. 1977:
$$f \left(T_y \right)_p = V^X e^{\left(X \left(1-V \right) \right)}$$
$$V = \left( T^{cm}_p - T_y \right)/ \left( T^{cm}_p - T^{co}_p \right)$$
$$X = \left(Z^2 \left(1+ \left( 1 + 40 /Y \right)^{0.5} \right)^2 \right)/400$$
$$Z = ln\left(Q^c_p \right) \left( T^{cm}_p - T^{co}_p \right)$$
$$Y = ln\left(Q^c_p \right) \left( T^{cm}_p - T^{co}_p + 2\right)$$

`Ceq = 3`: Temperature dependence for cool and cold-water species from Thornton and Lessem 1979:
$$f \left(T_y \right)_p = Ka * Kb$$
$$Ka = C^{K1}_p * L_1 / (1+C^{K1}_p * (L_1 - 1 ))$$
$$L_1 = e^{G_1 * (T_y - Q^c_p)}$$
$$G_1 = \frac{1}{T^{co}_p-Q^c_p}*\log\left(\frac{0.98*(1-CK1_p)}{CK1_p*0.02}\right)$$
$$Kb = \frac{CK4_p *L_2} {1+CK4_p*(L_2-1)}$$
$$L_2 = e^{G_2 (T^{cl}_p - T_y)}$$
$$G_2 = \frac{1}{T^{cl}_p-T^{cm}_p}*\log\left( \frac{0.98*(1-CK4_p)}{CK4_p*0.02}\right)$$

### Input ration-at-age
`Ceq = 4` and for directly input of consumption-at-sex/age in the `ration_data` data sheet for $RFR_{pbj,y}$. Set `CA = 1`, `fday = 1`, and `CB = 0`.
$$\delta_{pbj,y} = RFR_{pbj,y}* W_{pbj,y} $$


# Fishery and survey observation model
CEATTLE is fit to time series of fishery and survey biomass, sex-ratio, and length- or age-composition data. Log-fishery catch, log-survey biomass and sex-ratio are assumed to be normally distributed, while age- and length-composition data are fit assuming multinomial distributions. Variance parameters for the lognormal distributions and initial input sample size for the multinomial distributions can be assumed known or for surveys also estimated as a free parameter or estimated analytically following Walters and Ludwig (1994). Separate selectivity and catchability functions can be estimated or input for each survey or fishery.

Catch is estimated in the model following the Baranov catch equation:
$$\hat{C}_{isa,y} = \frac{F_{isa,y}} {Z_{isa,y}} \left( 1 - e^{-Z_{isa,y}} \right) N_{isa,y}$$
$$\hat{C}_{i,y} = \sum_s\sum_a^{A_i} \frac{F_isa,y} {Z_{isa,y}} \left( 1 - e^{-Z_{isa,y}} \right) N_{isa,y} W_{isa,y}$$
$$F_{isa,y} = \bar{F}_{f_i}e^{\epsilon_{{f_i},y} s_{{f_i}sa}}$$
CPUE data is estimated given catchability $q_{f_iy}$, selectivity $S_{f_isa,y}$, and fleet month $Month_{f_i}$ and can either be in biomass $\hat{CPUE:B}_{{f_i},y}$ or numbers $\hat{CPUE:I}_{{f_i},y}$:
$$\hat{CPUE:B}_{{f_i},y} = \sum_s\sum_a \left( N_{isa,y} e^{-Month_{f_i} Z_{isa,y}} W_{{f_i}sa,y} S_{{f_i}isa,y} * q_{{f_i}y} \right)$$
$$\hat{CPUE:I}_{{f_i},y} = \sum_s\sum_a \left( N_{isa,y} e^{-Month_{f_i} Z_{isa,y}} S_{{f_i}isa,y} * q_{{f_i}y} \right)$$


## Selectivity parameterizations
Multiple age-based selectivity functions are supported by Rceattle. They are controlled by `Selectivity_index`, `Selectivity`, `N_sel_bins`, `Time_varying_sel`, `Time_varying_sel_sd_prior`, and `Bin_first_selected` in fleet_control sheet of an Rceattle data file. They are defined as follows:

* `Selectivity_index`:	index to use if selectivities of different surveys are to be the same
* `Selectivity`:	Selectivity parameterization to use for the species: 0 = empirical selectivity provided in `srv_emp_sel`; 1 = logistic selectivity; 2 = non-parametric selecitivty sensu Ianelli et al 2018; 3 = double logistic; 4 = descending logistic; 5 = non-parametric selectivity sensu Taylor et al 2014 (Hake), 6 or "2DAR1" for age by year, 7 or "3DAR1" sensu Cheng et al (2025).
* `N_sel_bins`:	Number of ages to estimate non-parametric selectivity. For example, if `minage = 1` and selectivity parameters are estimated up till age-6, `N_sel_bins = 6` and if `minage = 0` and selectivity parameters are estimated up till age-6, `N_sel_bins = 7`.  Not used otherwise
* `Time_varying_sel`:	Whether a time-varying selectivity should be estimated for logistic, double logistic selectivity, descending logistic, or non-parametric selectivity (`Selectivity = 5 or "Hake"`). 0 = no, 1 = penalized deviates given sel_sd_prior, 2 = penalized deviates and estimate sel_sd_prior, 3 = time blocks with no penalty, 4 = random walk following Dorn, 5 = random walk on ascending portion of double logistic only. 
* `Time_varying_sel_sd_prior`:	The sd to use for the random walk of time varying selectivity if set to 1. 
* `Bin_first_selected`:	Age/length bin at which selectivity is non-zero. Selectivity before this age will be set to 0.

### Fixed selectivity
Input:
`Selectivity = 0 or "Fixed"` and input selectivity-at-age in the `emp_sel` data.

### Non-parametric (Type 1)
Input:
`Selectivity = 2 or "NonParametric"`, `N_sel_bins = 6`, `Sel_curve_pen1 = 12.5`, `Sel_curve_pen2 = 200`, and `Bin_first_selected = 1`

Equation:
$$sel_{f_i sa}=e^{\phi_{f_isa}}$$

### Logistic
Input:
`Selectivity = 1 or "Logistic"`, `N_sel_bins = NA`, `Time_varying_sel = 0`, and `Time_varying_sel_sd_prior = NA`

Equation:
$$sel_{f_i sa}=1/(1+e^{-Slp_{f_i} (a-Inf_{f_i} ) } )$$

### Time-varying logistic
Input:
`Selectivity = 1 or "Logistic"`, `N_sel_bins = NA`, and `Time_varying_sel_sd_prior = 0.1` or other desired value. 

Equation:
$$sel_{f_i sa}=1/(1+e^{-(Slp_{f_i}+e^{\phi_{f_i sy}^{slp} } ) (a-(Inf_{f_i}+\phi_{f_i sy}^{Inf} ) ) } )$$
Penalized likelihood: `Time_varying_sel = 1`
$$\phi_{f_i sy} \sim N(0,\sigma^{\phi}_{f_i}) $$
Random effect: `Time_varying_sel = 2` $\hat\sigma^{\phi}_{f_i}$ is estimated
$$\phi_{f_i sy} \sim N(0,\hat\sigma^{\phi}_{f_i}) $$
Block: `Time_varying_sel = 3` main parameters are set to 0 and the following is specificed:
$$\phi_{f_i sy} = \phi_{f_i sy_{block}}$$
Penalized likelihood random walk `Time_varying_sel = 4`
$$\phi_{f_i sy}-\phi_{f_i sy-1} \sim N(0,\sigma^{\phi}_{f_i}) $$

### Double logistic
Input:
`Selectivity = 3 or "DoubleLogistic"`, `N_sel_bins = NA`, `Time_varying_sel = 0`, and `Time_varying_sel_sd_prior = NA`

Equation:
$$sel_{f_i sa}=1/(1+e^{-Slp1_{f_i} (a-Inf1_{f_i} ) } )(1-1/(1+e^{-Slp2_{f_i} (a-Inf2_{f_i} ) } ))$$

### Time-varying double logistic 1	(ascending and descending time-varying parameters)
$$sel_{f_i sa}=1/(1+e^{-(Slp1_{f_i}+e^{\phi_{f_i sy}^{slp1} } ) (a-(Inf1_{f_i}+\phi_{f_i sy}^{Inf1} ) ) } )(1-1/(1+e^{-(Slp2_{f_i}+e^{\phi_{f_i sy}^{slp2} } ) (a-(Inf2_{f_i}+\phi_{f_i sy}^{Inf2} ) ) } ))$$

Penalized likelihood: `Time_varying_sel = 1`
$$\phi_{f_i sy} \sim N(0,\sigma^{\phi}_{f_i}) $$
Random effect: `Time_varying_sel = 2` $\hat\sigma^{\phi}_{f_i}$ is estimated
$$\phi_{f_i sy} \sim N(0,\hat\sigma^{\phi}_{f_i}) $$
Block: `Time_varying_sel = 3` main parameters are set to 0 and the following is specificed:
$$\phi_{f_i sy} = \phi_{f_i sy_{block}}$$
Penalized likelihood random walk `Time_varying_sel = 4`
$$\phi_{f_i sy}-\phi_{f_i sy-1} \sim N(0,\sigma^{\phi}_{f_i}) $$

### Time-varying double logistic 2	(random walk ascending-time varying parameters)
Input:
`Selectivity = 3 or "DoubleLogistic"`, `N_sel_bins = NA`, `Time_varying_sel = 5`, and `Time_varying_sel_sd_prior = NA`

Equation:
$$sel_{f_i sa}=1/(1+e^{-(Slp1_{f_i}+e^{\phi_{f_i sy}^{slp1} } ) (a-(Inf1_{f_i}+\phi_{f_i sy}^{Inf1} ) ) } )(1-1/(1+e^{-Slp2_{f_i} (a-Inf2_{f_i} ) } ))$$

$$\phi_{f_i sy} - \phi_{f_i sy-1} \sim N(0,\sigma^{\phi}_{f_i}) $$

### Descending logistic
Input:
`Selectivity = 4 or "DescendingLogistic"`, `N_sel_bins = NA`, `Time_varying_sel = 0`, and `Time_varying_sel_sd_prior = NA`

Equation:
$$sel_{f_i sa}=1-1/(1+e^{-Slp2_{f_i} (a-Inf2_{f_i} ) } )$$

### Time-varying descending logistic	
Input:
`Selectivity = 4 or "DescendingLogistic"`, `N_sel_bins = NA`, `Time_varying_sel = 1` or `2` , and `Time_varying_sel_sd_prior = NA`

Equation:
$$sel_{f_i sa}=1-1/(1+e^{-(Slp2_{f_i}+e^{\phi_{f_i sy}^{slp2} } ) (a-(Inf2_{f_i}+\phi_{f_i sy}^{Inf2} ) ) } )$$

Penalized likelihood: `Time_varying_sel = 1`
$$\phi_{f_i sy} \sim N(0,\sigma^{\phi}_{f_i}) $$
Random effect: `Time_varying_sel = 2` $\hat\sigma^{\phi}_{f_i}$ is estimated
$$\phi_{f_i sy} \sim N(0,\hat\sigma^{\phi}_{f_i}) $$
Block: `Time_varying_sel = 3` main parameters are set to 0 and the following is specificed:
$$\phi_{f_i sy} = \phi_{f_i sy_{block}}$$
Penalized likelihood random walk `Time_varying_sel = 4`
$$\phi_{f_i sy}-\phi_{f_i sy-1} \sim N(0,\sigma^{\phi}_{f_i}) $$

### Non-parametric (Type 2)
For each age $a \geq A_{min}$ there is an incremental selectivity parameter $p_{f_i say}$ for the fleet $f_i$. Selectivity at age is calculated as follows:

$$sel_{f_i sa}=e^{sel^{`}_{f_i sa} - max(sel^{`}_{f_i sa})}$$
where
$$ sel^{`}_{f_i sa} = \sum_{i=A_{min}}^a p_{f_i sa}$$
Selectivity is fixed at $sel_{f_i sa} = 0$ for $a < A_{min}$ and $p_{f_i sa} = 0$ for ages above `minage + N_sel_bins`, giving constant selectivity beyond the last estimated value. `minage` is specified in the data and is the minimum age of the modeled population.

Input:
`Selectivity = 5 or "Hake"`, `N_sel_bins = 6`, `Time_varying_sel = NA`,  `Time_varying_sel_sd_prior = NA`, and `Bin_first_selected = `$A_{min}$

### Time-varying non-parametric (Type 2 for hake)
As above, but a set of deviations on $p_a$ is used  to control annual changes in selectivity given a fixed or estimated standard deviation ($\sigma^{\phi}_{f_i}$):
$$p_{f_i say}=p_{f_i sa}+\phi_{f_i say}$$

Input:
`Selectivity = 5 or "Hake"`, `N_sel_bins = 6`, `Time_varying_sel = NA`,  `Time_varying_sel_sd_prior = NA`, and `Bin_first_selected = `$A_{min}$

Penalized likelihood: `Time_varying_sel = 1`
$$\phi_{f_i sy} \sim N(0,\sigma^{\phi}_{f_i}) $$
Random effect: `Time_varying_sel = 2` $\hat\sigma^{\phi}_{f_i}$ is estimated
$$\phi_{f_i sy} \sim N(0,\hat\sigma^{\phi}_{f_i}) $$


## Catchability parameterizations
Separate selectivity and catchability functions can be estimated for each survey or fishery. Additionally, catchability can be set the same for multiple surveys by setting `Q_index` to the same value. Catchability is controlled by the following parameters in the `fleet_control` sheet of the data

* `Q_index`	index to use if catchability coefficients are to be set the same
* `Catchability`	Text or integer switch specifying catchability. 0 or "Fixed" = fixed at prior; 1 or "Estimated" = Estimate single parameter; 2 or "Estimated-with-prior" = Estimate single parameter with prior; 3 or "Analytical" = Estimate analytical q  from Ludwig and Walters 1994; 4 or "PowerEquation" = Estimate power equation; - 5 or "Environmental" = Linear equation log(q_y) = q_mu + beta * index_y); 6 or "AR1" = annual AR1 catchability deviates are fit to environmental index sensu Rogers et al 2025.
* `Q_prior`	Starting value or fixed value for catchability
* `Q_sd_prior`	Variance of q prior: dnorm (log_q, log_q_prior, q_sd_prior)
* `Time_varying_q`	Whether a time-varying q should be estimated. 0 = no, 1 = penalized deviate, 2 = random effect, 3 = time blocks with no penalty; 4 = random walk from mean following Dorn 2018 (dnorm(q_y - q_y-1, 0, sigma). If Catchability = 5, this determines the environmental index to be used in the equation log(q_y) = q_mu + beta * index_y
* `Time_varying_q_sd_prior`	The sd to use for the random walk of time varying q if set to 1

### Linear time-invariant catchability formulations
`Time_varying_q = 0` for all formulations
`Catchability = 0`

$q_{f_i} =$ `Q_prior` from data

`Catchability = 1` A freely estimated catchability parameter $\hat{q}_{f_i}$ is estimated for survey/index fleet $f_i$

$$q_{f_i} = \hat{q}_{f_i}$$

`Catchability = 2` A freely estimated catchability parameter $\hat{q}_{f_i}$ is estimated for survey/index fleet $f_i$ assuming a lognormal prior

$$ \hat{q}_{f_i} \sim lognormal(log(Q\ prior), Q\ sd\ prior) $$

`Catchability = 3` Catchability $\hat{q}_{f_i}$ is analytically derived for survey/index fleet $f_i$ following Walters and Ludwig (1994). 
* For CPUE in weight with time-invariant standard error of the survey:
$$\hat{q}_{f_i} = \exp\left( \sum_y log \left[ \frac{CPUE:B_{{f_i},y}}  {\sum_s\sum_a \left( N_{isa,y} e^{-Month_{f_i} Z_{isa,y}} W_{{f_i}sa,y} S_{{f_i}isa,y} \right)} \right]/n \right)$$
* For CPUE in numbers with time-invariant standard error of the survey
$$\hat{q}_{f_i} = \exp\left( \sum_y log \left[ \frac{CPUE:I_{{f_i},y}}  {\sum_s\sum_a \left( N_{isa,y} e^{-Month_{f_i} Z_{isa,y}} S_{{f_i}isa,y} \right)} \right]/n \right)$$
* For CPUE in weight with time-varying standard error of the survey
$$\hat{q}_{f_i} = \exp\left( \sum_y \left( log \left[ \frac{CPUE:B_{{f_i},y}}  {\sum_s\sum_a N_{isa,y} e^{-Month_{f_i} Z_{isa,y}} W_{{f_i}sa,y} S_{{f_i}isa,y}} \right]/\sigma^{2}_{f_iy} \right)/ \sum_y{\frac{1}{\sigma^{2}_{f_iy}}} \right)$$
* For CPUE in numbers with time-varying standard error of the survey
$$\hat{q}_{f_i} = \exp\left( \sum_y \left( log \left[ \frac{CPUE:I_{{f_i},y}}  {\sum_s\sum_a N_{isa,y} e^{-Month_{f_i} Z_{isa,y}} S_{{f_i}isa,y}} \right]/\sigma^{2}_{f_iy} \right)/ \sum_y{\frac{1}{\sigma^{2}_{f_iy}}} \right)$$

### Time-varying catchability formulations
*Annual catchability estimated via penalized likelihood or as random effects* 

NOTE: If `Catchability = 2` a prior will be put on mean catchability as above.
`Catchability %in% c(1:2)` and `Time_varying_q %in% c(1:2)`
$$\hat{q}_{f_iy} = e^{\bar{q}_{f_i}+ \omega_{f_iy}}$$
$\omega_{f_iy} \sim N(0,$`Time_varying_q_sd_prior`$)$ where `Time_varying_q_sd_prior` is either fixed `Time_varying_q = 1` or estimated `Time_varying_q = 2` & `random_q = TRUE` in \function{fit_mod}

*Catchability time-blocks*

`Catchability %in% c(1:2)`, `Time_varying_q = 3` where the blocks are based on the blocks in `index_data` data sheet.
$$\hat{q}_{f_iy} = e^{\omega_{f_i block_y}}$$

*Annual catchability via random walk*

`Catchability %in% c(1:2)`, `Time_varying_q = 4`
$$\hat{q}_{f_iy} = e^{\bar{q}_{f_i}+ \omega_{f_iy}}$$
$\omega_{f_iy} \sim N(\omega_{f_iy-1},$`Time_varying_q_sd_prior`$)$

