---
title: "Guidelines for the stochastic production model in continuous time (SPiCT)"
author: "T.K. Mildenberger, A. Kokkalis, C.W. Berg"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
  rmarkdown::html_vignette:
    number_sections: false
    toc: false
  pdf_document:
    number_sections: false
    toc: false
    keep_tex: false
    citation_package: natbib    # or "biblatex", or just omit for default pandoc cites
bibliography: bibliography.bib
# optional: choose a bibliography style for PDF if you use natbib
# biblio-style: apalike

header-includes:
  - \usepackage[utf8]{inputenc}

vignette: >
  %\VignetteIndexEntry{SPiCT Guidelines}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, echo = FALSE, warning = FALSE, message = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "  ", fig.align = 'center', cache=FALSE)
library(knitr)
```
This is a living document, be sure to check for the latest update
([link](https://github.com/DTUAqua/spict/commits/master)). The SPiCT package
is actively developed, check for the most recent package version and
report problems at ([link](https://github.com/DTUAqua/spict/releases)).

| Document version | SPiCT version | Comments                              |
| :--------------- | :-----------: | ---------------------------:          |
|  v 1.0           | v 1.2.7       | Initial version (2019-09-23)          |
|  v 1.1           | v 1.2.8       | Management section added              |
|  v 1.2           | v 1.3.0       | Adjusted to package update            |
|  v 1.2           | v 1.3.4       | Adjusted after WKMSYSPICT and WKLIFEX |
|  v 1.3           | v 1.3.8.1     | Updated management section            |

This document provides specific guidelines for the use of the stochastic
production model in continuous time [SPiCT, @spict] and is divided into three
parts, containing: (i) the main assumptions and data requirements of SPiCT, (ii)
a checklist for the acceptance of a SPiCT assessment, (iii) options for
assessment tuning, and (iv) harvest control rules for SPiCT.

Accompanying sample code that assumes that

- input data is contained in a list called `inp`,
- the results are contained in a list called `fit` after fitting SPiCT to the
input data with `fit <- fit.spict(inp)`, and
- calculating the one-step-ahead residuals with `fit <-
calc.osa.resid(fit)`.


## Main assumptions and input data for SPiCT

- Catch data should be representative of both landings and bycatch. It is also
possible to use landings only, but then the interpretation of the results
changes. If available, seasonal catches should be used as input. Catches are
assumed to be taken over a time interval (e.g. years or quarters), thus the
associated time vector in SPiCT `inp$timeC` should reflect the beginning of each
catch interval (e.g. 2002.25 and 2002.75 for the second and fourth quarter
catches, respectively). Additionally, the vector `inp$dtc` should reflect the
length of each time interval (e.g. 1 for annual and 0.25 for quarterly catches,
respectively).

- Stock size indices should be in terms of biomass (not numbers) and
representative of the part of the stock vulnerable to the commercial fleets, the
so called exploitable stock biomass (ESB). In many cases, the gear selectivity
of the commercial and scientific fleets do not coincide and thus the stock size
indices have to be corrected to exclude individuals that are not represented in
the commercial fleets.

- Biomass indices are assumed to be snapshots at given points in time.
Therefore, the timing of survey indices `inp$timeI` has to be given as decimal
years reflecting the timing of the survey (e.g. 1995.5 for the middle of the
year). The timing of the survey will be matched to the closest model time which
is dependent on `inp$dteuler` (see below). Commercial CPUE indices should be
associated with the midpoint of the interval of the corresponding catches, i.e.
middle of the year if they are based on yearly aggregated catches and effort.

- If input data spans long periods and there is ecological evidence
for potential changes in productivity, it is possible to fit
productivity regime shifts or long-term gradual changes in
productivity with SPiCT [@tvp].

- The Euler discretisation has to be finer than the data, i.e. for yearly data
the Euler time step `inp$dteuler` has to be smaller than 1. A value of 1 changes
the model to a discrete time production model with different model assumptions.
The default value is 1/16 year, which is sufficient in most cases.


## Checklist for the acceptance of a SPiCT assessment

ICES category 3 stocks can be managed using the official advice rule 3.1.1 and
3.1.2 in @msycat34, which requires an accepted SPiCT assessment. Following
points should be considered for the acceptance of a SPiCT assessment.

1. The assessment **converged** (`fit$opt$convergence` equals $0$).

2. All **variance parameters** of the model parameters are **finite**
(`all(is.finite(fit$sd)`) should be `TRUE`).

3. **No violation of model assumptions** based on one-step-ahead
residuals (bias, auto-correlation, normality). This means, that
p-values are insignificant ($> 0.05$), indicated by green titles in
the graphs of `plotspict.diagnostic(fit)`. Slight violations of these
assumptions do not necessarily invalidate model results.

4. **Consistent patterns in the retrospective analysis** (`fit <-
retro(fit)`). This means that there is no tendency of consistent
under- or overestimation of the relative fishing mortality
($\frac{F}{F_{MSY}}$) and relative biomass ($\frac{B}{B_{MSY}}$) in
successive assessment. The retrospective trajectories of those two
quantities should be inside the credible intervals of the base run.

5. **Realistic production curve**. The shape of the production curve
should not be too skewed ($\frac{B_{MSY}}{K}$ should be between 0.1
and 0.9). Low values of $\frac{B_{MSY}}{K}$ allow for an infinite
population growth rate (`calc.bmsyk(fit)`).

6. **High assessment uncertainty** can indicate a lack of contrast in
the input data or violation of the ecological model assumptions. The
main variance parameters (`logsdb`, `logsdc`, `logsdi`, `logsdf`)
should not be unrealistically high. Credible intervals for
$\frac{B}{B_{MSY}}$ and $\frac{F}{F_{MSY}}$ should not span more than
1 order of magnitude (`calc.om(fit)`).

7. **Initial values do not influence the parameter estimates** (`fit
<- check.ini(fit)`). The estimates should be the same for all initial
values (`fit$check.ini$resmat`). Runs which did not converge should
not be considered in this regard.

## Optional model and assessment tuning

**Disclaimer**: An uncertain assessment is okay; the model should not
be tweaked and parameters should not be fixed unless there is
sufficient information and evidence to make such assumptions. Rather
than avoiding uncertainty, it should be accounted for by means of
stochastic harvest control rules [@wklifeix].

* Increase **iterations of optimisation**. The special error code `false
convergence (8)` could indicate that the optimisation exceeded the maximum
number of iterations `max.iter`. `inp$optimiser.control = list(iter.max = 1e5,
eval.max = 1e5)`

* If the catch time series is longer than the survey(s), **shortening
the catch time series** to cover only the period where there is an
available biomass index may help model convergence. `inp =
shorten.inp(inp, 2005, 2018)`

* **Adjust** the prior for the parameter `logn` determining the **shape of the
  production curve**.

    - Alternative priors for `logn` could be based on meta studies, e.g. for all
      species pooled based on @thorson2012: `inp$priors$logn <- c(log(1.478),
      0.6, 1)`

    - Tighter Schaefer prior for `logn` (or any other prior e.g. from meta
    studies): `inp$priors$logn <- c(log(2), 0.5, 1)`

    - Fixing n to resemble the Schaefer production model (or the meta study,
    alternatively): `inp$ini$logn <- log(2); inp$phases$logn <- -1`

Be aware that these modifications might cause smaller credible bands.

* Use a **prior for the initial depletion level** based on available
information, FOR EXAMPLE:

    - If evidence or expert knowledge allows to infer that there was **low or no**
    exploitation before the beginning of the available data: initial depletion
    level could be assumed to be close to the carrying capacity (**e.g.**
    `inp$priors$logbkfrac <- c(log(0.8),0.5,1)`)

    - If evidence or expert knowledge allows to infer that there was **high**
    exploitation before the beginning of the available data: initial depletion
    level could be assumed to be a certain fraction of the carrying capacity
    (**e.g.** `inp$priors$logbkfrac <- c(log(0.2),0.5,1)`)

* If information on the level of uncertainty in the biomass index and/or the
uncertainty of the catch is available prior distribution could be used for the
observation error term of the indices `logsdi` and catches `logsdc`,
respectively. This requires to remove the priors for the ratios of process to
observation errors `logalpha` and `logbeta`.

## Management based on an accepted SPiCT assessment

SPiCT estimates MSY based reference levels, which can be used to calculate
quantities relevant for fisheries management, such as the Total Allowable Catch
(TAC). The TAC can be calculated based on fishing at $F_{MSY}$ or based more
elaborate harvest control rules (HCR) which include additional components, e.g.
a linearly reduction in F if the biomass is below a certain reference level
(e.g., $\text{MSY} B_{\text{trigger}}$) as the hockey-stick MSY rule recommended
by @msycat34.

@mildenberger2022implementing demonstrated the importance of the
quantification and consideration of assessment uncertainty into the calculation
of the catch advice. The authors introduced a modification of the hockey-stick
MSY rule that reduces catch advice as a function of uncertainty by using a
percentile lower than 50 (i.e., median) for relative biomass
$\frac{B}{B_{\text{MSY}}}$ and catch $C_{y+1}$ and higher than 50 for relative
fishing mortality $\frac{F}{F_{\text{MSY}}}$.

Based on this work, ICES currently recommends calculating TAC based on an MSY
hockey-stick rule with a biomass trigger reference point at
$\text{MSY} B_{\text{trigger}} = 0.5 B_{\text{MSY}}$, a biomass limit reference
point at $B_{\text{lim}} = 0.3 B_{\text{MSY}}$, and the 35th percentile of the
predicted catch distribution [@ices2025advice]. Under this rule, the target
fishing mortality $F^{\tau}$ equals $F_{\text{MSY}}$ when $B$ is above
$\text{MSY}B_{\text{trigger}}$; when $B$ is between 0 and
$\text{MSY}B_{\text{trigger}}$, $F^{\tau}$ is reduced linearly with $B$; and
when $B$ falls below $B_{\text{lim}}$, $F^{\tau}$ is set to zero. This can be
expressed as

\begin{equation*}
F^{\tau}_y =
\begin{cases}
F_{\text{MSY}} \min\!\left(1, \dfrac{B_y}{\text{MSY}B_{\text{trigger}}}\right),
& \text{if } B_y \ge B_{\text{lim}},\\[6pt]
0, & \text{if } B_y < B_{\text{lim}}.
\end{cases}
\end{equation*}

Given the target fishing mortality, the advised catch $A_y$ for year $y$ is
taken as the 35th percentile of the predicted catch distribution:

\begin{equation*}
A_y = \Phi^{-1}_{C_{\text{pred},y} \mid F^{\tau}_y}(0.35),
\end{equation*}

where $\Phi^{-1}$ denotes the inverse cumulative distribution function, so that
$\Phi^{-1}_{C_{\text{pred},y} \mid F^{\tau}_y}(0.35)$ is the catch
corresponding to the 35th percentile of the predicted catch distribution.

This HCR is one of the 8 default management scenarios of the `manage()` function
in spict but can als be run exclusively with:

`fit <- manage(fit, scenarios = "ices")`.

This function is actually only a wrapper function and calls internally the
`add.man.scenario()` function that makes the three key components of the ICES
(2025) HCR explicit:

`fit <- add.man.scenario(fit, "ICES (2025)", fractiles = list(catch = 0.35),
breakpointB = 0.5, limitB = 0.3)`.

Another way to calculate the TAC based on the ICES (2025) HCR directly is using
the `get.TAC()` function:

`get.TAC(fit, fractiles = list(catch = 0.35), breakpointB = 0.5, limitB = 0.3)`.

It is important to note that informative priors (i.e.\ priors with smaller
standard deviations) can influence parameter estimates and credible intervals in
SPiCT assessments, and thereby affect the resulting TAC.

## References
