---
title: "2a. Mathematical description of the SRA model"
author: "Quang Huynh (<q.huynh@oceans.ubc.ca>)"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    toc: true
    number_sections: true
vignette: >
  %\VignetteIndexEntry{2a. Mathematical description of the SRA model}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
---
<script type="text/x-mathjax-config">
  MathJax.Hub.Config({ TeX: { equationNumbers: {autoNumber: "all"} } });
</script>

<style type="text/css">

body{ /* Normal  */
   font-size: 12px;
}
td {  /* Table  */
   font-size: 8px;
}
h1 { /* Header 1 */
 font-size: 18px;
 color: DarkBlue;
}
h2 { /* Header 2 */
 font-size: 15px;
 color: DarkBlue;
}
h3 { /* Header 3 */
 font-size: 14px;
 color: DarkBlue;
}
code.r{ /* Code block */
  font-size: 10px;
}
pre { /* Code block */
  font-size: 10px
}
</style>


```{r set options, echo = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
knitr::opts_chunk$set(dpi=85)
options(width = 650)
```

<br>

# Selectivity and mortality
For fleets, selectivity is defined by blocks (indexed by $b$) which can then be assigned to any fleet $f$ and year $y$ to allow for time-varying selectivity. By default, each fleet is assigned to its own block for all years (no time-varying selectivity).

For flat-topped selectivity in block $b$, two parameters are used and expressed in terms of length units: the length of 5% selectivity ($L^5_b$) and the length of full selectivity $L^{\textrm{FS}}_b$. For dome selectivity, a third parameter, the selectivity at $L_{\infty}$, $V^{L_{\infty}}_b$ is also used. Length-based selectivity is converted to age-based selectivity in the age-structured model as:

$$ v_{y,a,b} = 
\begin{cases}
2^{-[(L_{y,a} - L^{\textrm{FS}}_b)/(\sigma^{\textrm{asc}}_b)]^2} & \textrm{if } L_{y,a} < L^{\textrm{FS}}_b\\ 
1 & \textrm{if logistic and } L_{y,a} \ge L^{\textrm{FS}}_b,\\
2^{-[(L_{y,a} - L^{\textrm{FS}}_b)/(\sigma^{\textrm{des}}_b)]^2} & \textrm{if dome and } L_{y,a} \ge L^{\textrm{FS}}_b
\end{cases}
$$


where $L_{y,a}$ is the mean length-at-age and $\sigma^{\textrm{asc}}_b = (L^5_b - L^{\textrm{FS}}_b)/\sqrt{-\log_2(0.05)}$ and $\sigma^{\textrm{des}}_b = (L_{\infty} - L^{\textrm{FS}}_b)/\sqrt{-\log_2(V^{L_{\infty}}_b)}$ control the shape of the ascending and descending limbs, respectively, of the selectivity function. In this parameterization, length-based selectivity is constant over time. The corresponding age-based selectivity within each block is constant so long as growth is not time-varying.

Selectivity can also be parameterized where $v_{y,a,b}$ are free independent parameters. In this case, selectivity does not vary among years.

Total mortality $Z$ in year $y$ and for age $a$ is the sum of fishing mortality $F$ from all fleets and natural mortality $M$,

$$ Z_{y,a} = M_{y,a} + \Sigma_f v_{y,a,f} F_{y,f},$$

where $v_{y,a,f}$ is the fleet selectivity after assigning blocks to fleets.

# Survey selectivity
Survey selectivity is constant over time and is denoted as $v_{a,s}$ with either logistic, dome, or free parameterizations.

# Initial population distribution
The population age distribution in the first year of the model $y=1$ is in equilibrium where
$$ N_{1,a} = 
\begin{cases}
      R^{\textrm{eq}} \exp(-\Sigma_{i=1}^{a-1}Z^{\textrm{eq}}_i) & a = 1, \ldots, A-1\\
      \dfrac{R^{\textrm{eq}} \exp(-\Sigma_{i=1}^{a-1}Z^{\textrm{eq}}_i)}{1 - \exp(-Z^{\textrm{eq}}_A)} & a = A,
\end{cases}
$$
where the $R^{\textrm{eq}}$ is the equilibrium recruitment and $Z^{\textrm{eq}}_a = M_{1,a} + \Sigma_f v_{1,a,f} F^{\textrm{eq}}_f$ is the equilibrium total mortality rate. Unfished conditions are modeled by setting $F^{\textrm{eq}}_f = 0$. To estimate $F^{\textrm{eq}}_f$, the corresponding equilibrium catch in weight $\tilde{C}^{\textrm{eq}}_f$ prior to the first year of the model should be provided. In the equilibrium yield curve,  $F^{\textrm{eq}}_f$ would be the fishing mortality corresponding to fishing at $F^{\textrm{eq}}_f$. Once $Z^{\textrm{eq}}_a$ is obtained, then the equilibrium recruitment is calculated as:

$$ R^{\textrm{eq}} = 
\begin{cases}
      \dfrac{\alpha^{\textrm{B}}\phi^{\textrm{eq}} - 1}{\beta^{\textrm{B}}\phi^{\textrm{eq}}} & \textrm{if Beverton-Holt stock-recruit relationship}\\
      \dfrac{\log(\alpha^{\textrm{R}}\phi^{\textrm{eq}})}{\beta^{\textrm{R}}\phi^{\textrm{eq}}} & \textrm{if Ricker stock-recruit relationship}
\end{cases},
$$
where $\phi^{\textrm{eq}}$ is the spawners-per-recruit when the mortality is $Z^{\textrm{eq}}_a$. From steepness $h$, $\alpha^{\textrm{B}} = \frac{4h}{(1-h)\phi_0}$, $\beta^{\textrm{B}} = \frac{5h-1}{(1-h)B^S_0}$, $\alpha^{\textrm{R}} = \frac{(5h)^{1.25}}{\phi_0}$, $\beta^{\textrm{R}} = \frac{\log(5h)}{B^S_0}$, where $\phi_0$ and $B^S_0$ are unfished spawners-per-recruit and unfished spawning biomass, respectively.

# Dynamics equations
After setting the equilibrium population age distribution in the first year of the model, the population abundance $N_{y,a}$ in subsequent years is
$$ N_{y,a} =
\begin{cases}
      R_y & a = 1\\
      N_{y-1,a-1} \exp(-Z_{y-1,a-1}) & a = 2, \ldots, A - 1,\\
      N_{y-1,a-1} \exp(-Z_{y-1,a-1}) + N_{y-1,a} \exp(-Z_{y-1,a}) & a = A
\end{cases}
$$
where $R_y$ is the recruitment and $A$ is the maximum-age as the plus-group. Recruitment is modelled as
$$ R_y = 
\begin{cases}
      \dfrac{\alpha^{\textrm{BH}} B^S_{y-1}}{1 + \beta^{\textrm{BH}}B^S_{y-1}} \exp(\delta_y - 0.5 \tau^2) & \textrm{if Beverton-Holt stock-recruit relationship}\\
      \alpha^{\textrm{Ricker}} B^S_{y-1} \exp(-\beta^{\textrm{Ricker}} B^S_{y-1})\exp(\delta_y - 0.5 \tau^2) & \textrm{if Ricker stock-recruit relationship}
\end{cases},
$$
where $\delta_y$ are recruitment deviates and $\tau$ is the standard deviation of the deviates.

The spawning biomass is $B^S_y$ is
$$B^S_y = \sum_a w_{y,a} m_{y,a} N_{y,a},$$
where $m_{y,a}$ and $w_{y,a}$ are the maturity at age and weight at age, respectively.

# Catch at age
The catch (in numbers) $C^N$ at age for fleet $f$ is
$$ C^N_{y,a,f} = \dfrac{v_{y,a,f} F_{y,f}}{Z_{y,a}} N_{y,a} (1 - \exp[-Z_{y,a}]).$$

If the model is conditioned on catch, then $F_{y,f}$ can be estimated as parameters or solved iteratively to match the observed catch. If the model is conditioned on effort, then
$$ F_{y,f} = q_f E_{y,f},$$ where $E_{y,f}$ is the observed effort and $q_f$ is the scaling coefficient.

# Catch-at-length 
The catch at length is calculated assuming a normally distributed length-at-age $P(\ell,a)$, where
$$ C^N_{y,\ell,f} = \sum_a C^N_{y,a,f} P(\ell|a) $$
and

$$ P(\ell|a) = 
\begin{cases}
      \phi(L'_{\ell+1}) & \ell = 1\\
      \phi(L'_{\ell+1}) - \phi(L'_\ell) & \ell = 2, \ldots, L - 1,\\
      1 -\phi(L'_\ell) & \ell = L
\end{cases}
$$
with $L'_{\ell}$ as the length at the lower boundary of length bin $\ell$ and $\phi(L'_{\ell})$ as the cumulative distribution function of a normal variable with mean $\tilde{L}_{y,a}$ (the expected mean length at age $a$) and standard deviation $\tilde{L}_{y,a} \times CV^L$ ($CV^L$ is the coefficient of variation in mean length at age). 

The catch in weight $\tilde{C}$ is
$$ \tilde{C}_{y,f} = \sum_a C^N_{y,a,f} w_{y,a}.$$

The mean length of the catch $\bar{L}_{y,f}$ is
$$ \bar{L}_{y,f} = \dfrac{\sum_{\ell} L_{\ell} C^N_{y,\ell,f}}{\sum_{\ell} C^N_{y,\ell,f}},$$
where $L_\ell$ is the midpoint of the length bin $\ell$.

The proportion of the catch-at-age is
$$ p_{y,a,f} = \dfrac{C^N_{y,a,f}}{\sum_a C^N_{y,a,f}}.$$

The proportion of the catch-at-length is
$$ p_{y,\ell,f} = \dfrac{C^N_{y,\ell,f}}{\sum_{\ell}C^N_{y,\ell,f}}.$$

# Survey 
If the $s^{\textrm{th}}$ survey is biomass-based, then the survey value $I_{y,s}$ is calculated as
$$ I_{y,s} = q_s \sum_a v_{y,a,s} N_{y,a} w_{y,a}, $$
where $q$ is the scaling coefficient and $s$ indexes survey.

If the survey is abundance-based, then 
$$ I_{y,s} = q_s \sum_a v_{y,a,s} N_{y,a} . $$

The proportion-at-age vulnerable to the survey is
$$ p_{y,a,s} = \dfrac{v_{a,s}N_{y,a}}{\sum_a v_{a,s}N_{y,a}}.$$

The proportion-at-length vulnerable to the survey is
$$ p_{y,\ell,s} = \dfrac{\sum_a v_{a,s} N_{y,a} P(\ell|a)}{\sum_{\ell} \sum_a v_{a,s} N_{y,a} P(\ell|a)}.$$

# Likelihoods

If the model is conditioned on catch and fishing mortality rates are estimated parameters, then the log-likelihood component $\Lambda_1$ of the catch is
$$\Lambda_1 = \sum_f \left[\lambda^{\tilde{C}}_f \sum_y \left(-\log(0.01) - \dfrac{[\log(\tilde{C}^{\textrm{obs}}_{y,f}) - \log(\tilde{C}^{\textrm{pred}}_{y,f})]^2}{2 \times 0.01^2}\right)\right],$$

where $\textrm{obs}$ and $\textrm{pred}$ indicate observed and predicted quantities, respectively, and $\lambda$ are likelihood weights. With a small standard deviation for the catch likelihood relative to the variance in other likelihood components, the predicted catch should match the observed catch.

The log-likelihood component $\Lambda_2$ of survey data is
$$\Lambda_2 = \sum_s \left[ \lambda^I_s \sum_y \left(-\log(\sigma_{y,s}) - \dfrac{[\log(I^{\textrm{obs}}_{y,s}) - \log(I^{\textrm{pred}}_{y,s})]^2}{2\sigma_{y,s}^2}\right) \right].$$

The log-likelihood component $\Lambda_3$ of catch-at-age data is
$$\Lambda_3 = \sum_f \lambda^A_f \left[\sum_y O^A_{y,f} \sum_a p^{\textrm{obs}}_{y,a,f} \log(p^{\textrm{pred}}_{y,a,f})\right],$$
where $O^A$ is the annual sample sizes for the age compositions.

The log-likelihood component $\Lambda_4$ of catch-at-length data is
$$\Lambda_4 = \sum_f \lambda^L_f \left[ \sum_y O^L_{y,f} \sum_{\ell} p^{\textrm{obs}}_{y,\ell,f} \log(p^{\textrm{pred}}_{y,\ell,f})\right]$$
where $O^L$ is the annual sample sizes for the length compositions.

The log-likelihood component $\Lambda_5$ of observed mean lengths in the catch is
$$\Lambda_5 = \sum_f \lambda^{\bar{L}}_f\left[ \sum_y \left(-\log(\omega_f) - \dfrac{[\bar{L}^{\textrm{obs}}_{y,f} - \bar{L}^{\textrm{pred}}_{y,f}]^2}{2 \omega^2_f}\right)\right],$$
where $\omega_f$ is the standard deviation of mean lengths.

The log-likelihood component $\Lambda_6$ of annual estimated recruitment deviates $\delta_y$ in log space is
$$\Lambda_6 = \Sigma_y\left(-\log(\tau) - \dfrac{\delta_y^2}{2 \tau^2}\right),$$
where $\tau$ is the standard deviation of recruitment deviates.

The log-likelihood component $\Lambda_7$ of the equilibrium catch is
$$\Lambda_7 = \sum_f \lambda^{\tilde{C}}_f \left(-\log(0.01) - \dfrac{[\log(\tilde{C}^{\textrm{eq,obs}}_f) - \log(\tilde{C}^{\textrm{eq,pred}}_f)]^2}{2 \times 0.01^2}\right),$$

The log-likelihood component $\Lambda_8$ of the survey proportion-at-age data is
$$\Lambda_8 = \sum_s \lambda^{IA}_s \left[\sum_y O^{IA}_{y,s} \sum_a p^{\textrm{obs}}_{y,a,s} \log(p^{\textrm{pred}}_{y,a,s})\right],$$
where $O^{IA}$ is the annual sample sizes for the survey age compositions.

The log-likelihood component $\Lambda_9$ of the surey proportion-at-length data is
$$\Lambda_9 = \sum_s \lambda^{IL}_s \left[ \sum_y O^{IL}_{y,s} \sum_{\ell} p^{\textrm{obs}}_{y,\ell,s} \log(p^{\textrm{pred}}_{y,\ell,s})\right]$$
where $O^{IL}$ is the annual sample sizes for the survey length compositions.

The total log-likelihood $\textrm{LL}$ to be maximized is
$$\textrm{LL} = \sum_{i=1}^9\Lambda_i.$$

# Estimated parameters

The estimated parameters, denoted in this section as $x$, are unconstrained over all real numbers and then transformed in order to constrain the corresponding model parameters. For optimization, the transformation is also designed to reduce the scale of all estimated parameters to within an order of magnitude.

## Selectivity

For a fleet block $b$ for which selectivity is estimated, then parameters $x^{LFS}_b$ and $x^{L5}_b$ are estimated over all real numbers, where

$$ \begin{align}
L^{\textrm{FS}}_b &= 0.99 \times L_{\infty} \times \textrm{logit}^{-1}(x^{LFS}_b)\\
L^5_b &= L^{\textrm{FS}}_b - \exp(x^{L5}_b) 
\end{align}$$

If a third parameter $x^{V}_b$is estimated for dome selectivity, then
$$ V^{L_{\infty}}_b = \textrm{logit}^{-1}(x^V_b)$$

If selectivity is parameterized as free parameters, then
$$ v_{y,a,f} = \textrm{logit}^{-1}(x^v_{y,a,f}).$$

For surveys, parameterizations are identical except with indexing for survey $s$.

## Fishing mortality

If $F_{y,f}$ are estimated parameters (`condition = "catch"`), then one parameter $x^F_f$ is the estimated $F$ in log-space in the middle of the time series is estimated and all others are subsequent deviations, represented as $x^{F_{dev}}_{y,f}$:

$$ F_{y,f} = 
\begin{cases}
\exp(x^F_f) & y \textrm{ is midpoint of the time series}\\ 
\exp(x^F_f) \times \exp(x^{F_{dev}}_{y,f}) & \textrm{otherwise}\\
\end{cases}
$$

If `condition = "effort"`, then $q_f$ is estimated in log space, where $$F_{y,f} = q_f E_{y,f} = \exp(x^q_f) \times E_{y,f}$$

## Index catchability

To scale biomass to index values, the index catchability $q_s$ is solved analytically in the model:

$$ q_s = \exp\left(\dfrac{\sum_y \log(I^{\textrm{obs}}_{y,s}) - \sum_y \log(\sum_a v_{y,a,s}N_{y,a,s})}{n_s}\right),$$
or
$$ q_s = \exp\left(\dfrac{\sum_y \log(I^{\textrm{obs}}_{y,s}) - \sum_y \log(\sum_a v_{y,a,s}N_{y,a,s}w_{y,a})}{n_s}\right),$$
for an abundance or biomass based index, respectively, where $n_s$ is the number of years with index values and the summation is over those $n_s$ years.

## Other parameters

Unfished recruitment is estimated in log-space, $R_0 = \dfrac{1}{z}\exp(x^{R_0})$ where $z$ is an optional rescaler, e.g. mean historical catch, to reduce the magnitude of the $x^{R_0}$ estimate. Recruitment deviations $\delta_y$ are directly estimated.

# Additional resources

The help file for `SRA_scope` will provide more information on the possible inputs for the model. The help file for the SRA class (obtained by typing `class?SRA` into the R console) describes the outputs from the function. An additional <a href="SRA_scope_sel.html">vignette</a> is available to describe set up of fleet and survey selectivity in the function call.

