---
title: "cat3advice"
author: "Simon H. Fischer"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output: 
  rmarkdown::github_document:
      toc: true
      toc_depth: 2
bibliography: references.bib
vignette: >
  %\VignetteIndexEntry{cat3advice}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: console
---

```{r, include = FALSE, warning=FALSE, message=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

# `cat3advice`

## Vignette change log

+---------+----------------------+------------+------------------------------------------------------------------------------+
| Version | `cat3advice` version | Date       | Changes                                                                      |
+=========+======================+============+==============================================================================+
| 1       | 0.0.7                | 2023       |                                                                              |
+---------+----------------------+------------+------------------------------------------------------------------------------+
| 2       | 0.1                  | 2025-04-01 | Updated vignette to latest changes in packages;                              |
|         |                      |            |                                                                              |
|         |                      |            | Added example of chr rule with tuned control parameters and discard survival |
+---------+----------------------+------------+------------------------------------------------------------------------------+

# cat3advice R package

`cat3advice` is an R package that allows the application of the ICES category 3 data-limited harvest control rules (rfb/rb/chr rules) and follows the ICES Technical Guidelines [@ICES2022_cat23_tech_guidelines] (<https://doi.org/10.17895/ices.pub.28506179>).

# Documentation

The package documentation contains help files for its functions which also include code examples. See `help(package = "cat3advice")` for available functions.

The main functions are `rfb()`, `rb()`, and `chr()`. Each of these functions has a help file with code examples (see `?rfb`, `?rb` `?chr`).

# Installation

The easiest way to install the `cat3advice` package is to install it as a binary package from ICES r-universe:

```{R, eval = FALSE}
install.packages("cat3advice", repos = c("https://ices-tools-prod.r-universe.dev", "https://cran.r-project.org"))
```

It is also possible to install it directly from the GitHub repository. However, this means building the package locally and requires the necessary build tools (e.g. on Windows RTools <https://cran.r-project.org/bin/windows/Rtools/> and the `devtools` R package).

```{R, eval = FALSE}
library(remotes)
install_github("shfischer/cat3advice")
```

# Tutorial

This tutorial uses data from the ICES Western English Channel plaice stock (ple.27.7e) to illustrate the application of the rfb/rb/chr rules. The data are included in the `cat3advice` R package.

Before reading this vignette, please first read the ICES Technical Guidelines [@ICES2022_cat23_tech_guidelines].

For more details on the rfb rule, please refer to Fischer et al. [-@Fischer2020_catch_rule; -@Fischer2021_GA; @Fischer2021_GA_PA; -@Fischer2022_risk_equivalence] and for the chr rule, please refer to Fischer et al. [-@Fischer2022_hr; -@Fischer2022_risk_equivalence].

```{r, message=FALSE, warning=FALSE}
### load package
library(cat3advice)
```

# The rfb rule

The rfb rule is an index adjusted harvest control rule that uses a biomass index and catch length data. The method is defined as Method 2.1 in the ICES Technical Guidelines [@ICES2022_cat23_tech_guidelines] as

$$
A_{y+1} = A_y \times r \times f \times b \times m
$$

where $A_{y+1}$ is the new catch advice, $A_y$ the previous advice, $r$ the biomass ratio from a biomass index, $f$ the fishing pressure proxy from catch length data, $b$ a biomass safeguard and $m$ a precautionary multiplier. Furthermore, the change in the advice is restricted by a stability clause that only allows changes of between $+20\%$ and $-30\%$ relative to the previous advice, but the application of the stability clause is conditional on $b=1$ and turned off when $b<1$.

The rfb rule should be applied biennially, i.e. the catch advice is valid for two years.

Please note that any change from the default configuration should be supported by case-specific simulations.

## Reference catch $A_y$

The reference catch $A_y$ is usually the most recently advised catch. In a typical ICES setting, an assessment is conducted in an assessment (intermediate) year $y$ to give advice for the following advice year $y+1$, this is the advice for year $y$:

```{r}
### load plaice catch and advice data
data("ple7e_catch")
tail(ple7e_catch)
### get reference catch
A <- A(ple7e_catch, units = "tonnes")
A
```

The ICES Technical Guidelines [@ICES2022_cat23_tech_guidelines] specify that if the realised catch is very different from the advised catch, the reference catch could be replaced by an average of recent catches:

```{r}
### use 3-year average
A(ple7e_catch, units = "tonnes", basis = "average catch", avg_years = 3)
```

The reference catch can also be defined manually:

```{r}
### define manually
A(2000, units = "tonnes")
```

## Biomass index trend (ratio) $r$

The biomass index trend $r$ calculates the trend in the biomass index over last the five years, by dividing the mean of the last two years by the mean of the three preceding years:

$$
r = \Sigma_{i=y-2}^{y-1}(I_i/2)\ / \ \Sigma_{i=y-5}^{y-3}(I_i/3)
$$

The ratio is calculated with the function `r`. Index data should be provided as a `data.frame` with columns `year` and `index`.

```{r, fig.align='center', fig.width=3, fig.height=2, dpi=300, out.width=500}
### load plaice data
data("ple7e_idx")
tail(ple7e_idx)

### calculate biomass trend
r <- r(ple7e_idx, units = "kg/hr")
r

### ICES advice style table
advice(r)

### plot index
### horizontal orange lines correspond to the the 2/3-year averages
plot(r)

### when the value of r is known
r(1)
```

Biomass index data are usually available until the year before the advice year. More recent data can be used and the function automatically picks the most recent data provided to it.

## Biomass safeguard $b$

The biomass safeguard reduces the advice when the biomass index $I$ falls below a threshold value $I_{\text{trigger}}$:

$$
b = \text{min}\{1,\ I_{y-1}/I_{\text{trigger}}\}
$$ By default, the trigger value is defined based on the lowest observed index value $I_{\text{loss}}$ as $I_{\text{trigger}} = 1.4I_{\text{loss}}$.

The biomass safeguard is calculated with the function `b`:

```{r, fig.align='center', fig.width=3, fig.height=2, dpi=300, out.width=500}
### use same plaice data as before
### application in first year with new calculation of Itrigger
b <- b(ple7e_idx, units = "kg/hr")
b

### ICES advice style table
advice(b)

### plot
plot(b)
### plot b and r in one figure
plot(r, b)
```

**Please note that** $I_{\text{trigger}}$ should only be defined once in the first year of the application of the rfb rule. In the following years, the same value should be used. For this purpose, `b` allows the direct definition of $I_{\text{trigger}}$, or, more conveniently, $I_{\text{trigger}}$ can be based on the year in which $I_{\text{loss}}$ is defined:

```{r}
### in following years, Itrigger should NOT be updated
### i.e. provide value for Itrigger
b(ple7e_idx, units = "kg/hr", Itrigger = 0.3924721)
### alternatively, the reference year for Iloss can be used
b(ple7e_idx, units = "kg/hr", yr_ref = 2007)
```

## Fishing pressure proxy $f$

Catch length data are used to approximate the fishing pressure. The mean length of fish in the catch compared to a reference length is used as the indicator.

### Length data

The fishing pressure proxy requires length data from the catch. Ideally, length data for several years are provided in a `data.frame` with columns `year`, `length` and `numbers`. An additional column `catch_category` specifying the catch category, such as discards and landings, is optional.

```{r}
data("ple7e_length")
head(ple7e_length)
```

### Length at first capture $L_c$

Only length data above the length at first capture $L_c$ are used to avoid noisy data from fish that are not fully selected. $L_c$ is defined as the first length class where the abundance is at or above 50% of the mode of the length distribution and can be calculated with the function `Lc()`:

```{r fig.align='center', fig.height=4, fig.width=5, warning=FALSE, dpi=300, out.width=500}
lc <- Lc(ple7e_length)
lc
plot(lc)
```

$L_c$ can change from year to year. Therefore, it is recommended to pool data from several (e.g. 5) years:

```{r, fig.align='center', fig.width=4, fig.height=3, dpi=300, out.width=500, warning=FALSE}
lc <- Lc(ple7e_length, pool = 2017:2021)
lc
plot(lc)
```

If length data are noisy, the size of the length classes can be increased:

```{r, fig.align='center', fig.width=4, fig.height=3, dpi=300, out.width=500, warning=FALSE}
### use 20mm length classes
plot(Lc(ple7e_length, pool = 2017:2021, lstep = 20))
```

Once defined, $L_c$ should be kept constant and the same value used for all data years. $L_c$ should only be changed if there are strong changes in the fishery or data.

### Mean length

After defining $L_c$, the mean (annual) catch length $L_{\text{mean}}$ above $L_c$ can be calculated:

```{r, fig.align='center', fig.width=5, fig.height=3, dpi=300, out.width=700}
### calculate annual mean length
lmean <- Lmean(data = ple7e_length, Lc = lc, units = "mm")
lmean
plot(lmean)
```

If length data are noisy, the size of the length classes can be increased:

```{r, fig.align='center', fig.width=5, fig.height=3, dpi=300, out.width=700}
### use 20mm length classes
plot(Lmean(data = ple7e_length, Lc = lc, units = "mm", lstep = 20))
```

### Reference length

The reference length follows the concepts of @Beverton1957_dynamics_populations and is calculated as derived by @Jardim2015_HCR:

$$
L_{F=M} = 0.75L_c + 0.25L_\infty
$$ where $L_{F=M}$ is the MSY reference length, $L_c$ the length at first capture as defined above, and $L_\infty$ the von Bertalanffy asymptotic length. This simple equation assumes that fishing at $F=M$ can be used as a proxy for MSY and that $M/k=1.5$ (where $M$ is the natural mortality and $k$ the von Bertalanffy individual growth rate). The reference length can be calculated with

```{r}
lref <- Lref(Lc = 264, Linf = 585, units = "mm")
lref
```

The reference length $L_{F=M}$ should only be defined once in the first year of the application of the rfb rule. In the following years, the same value should be used.

Deviations from the assumptions of $F=M$ and $M/k=1.5$ are possible following Equation A.3 of @Jardim2015_HCR $L_{F=\gamma M, k = \theta M}=(\theta L_\infty + L_c(\gamma + 1))/(\theta + \gamma +1)$ and can be used by providing arguments `gamma` and `theta` to `Lref()`.

### Indicator

The length indicator $f$ is defined as

$$
f = L_{\text{mean}}/L_{F=M}
$$

and can be calculated with `f()`:

```{r, fig.align='center', fig.width=4, fig.height=3, dpi=300, out.width=500}
f <- f(Lmean = lmean, Lref = lref, units = "mm")
f
### ICES advice style table
advice(f)
### plot
plot(f)
```

In this case, the mean catch length (orange curve) is always below the MSY proxy reference length (blue horizontal line), indicating that the fishing pressure was above $HR_{\text{MSY proxy}}$.

## Multiplier $m$

The multiplier $m$ is a tuning parameter and ensures that the catch advice is precautionary in the long term.

The value of $m$ is set depending on the von Bertalanffy parameter $k$ (individual growth rate), with $m=0.95$ for stocks with $k<0.2/\text{year}$ and $m=0.90$ for stocks with $0.2\text{year}^{-1} \leq k < 0.32\text{year}^{-1}$. The multiplier can be calculated with the function `m()`:

```{r}
# for k=0.1/year
m <- m(hcr = "rfb", k = 0.1)
m

### alias for rfb rule
rfb_m(k = 0.1)

### ICES advice style table
advice(m)
```

Please note that the multiplier $m$ does not lead to a continuous reduction in the catch advice. The components of the rfb rule are multiplicative, this means that $m$ could be considered as part of component $f$ and essentially adjusts the reference length $L_{F=M}$ to $L'_{F=M}$:

$$
A_{y+1} = A_y\ r\ f\ b\ m = A_y\ r\ \frac{L_{\text{mean}}}{L_{F=M}}\ b\ m = A_y\ r\ \frac{L_{\text{mean}}}{L_{F=M}/m}\ b = A_y\ r\ \frac{L_{\text{mean}}}{L'_{F=M}}\ b
$$

## Application of rfb rule

Now we have all the components of the rfb rule and can apply it:

```{r, fig.align='center', fig.width=4, fig.height=3, dpi=300, out.width=500}
advice <- rfb(A = A, r = r, f = f, b = b, m = m, discard_rate = 27)
advice
```

A discard rate in % can be provided to the argument `discard_rate` and this means the advice is provided for the catch and landings.

The rfb rule includes a stability clause that restricts changes to $+20\%$ and $-30\%$. This stability clause is conditional on the biomass safeguard and is only applied if $b=1$ but turned off when $b<1$.

`cat3advice` can print a table similar to the table presented in ICES advice sheets in which all numbers are rounded following ICES rounding rules:

```{r}
advice(advice)
```

# The rb rule

The rb rule is essentially a simplified version of the rfb rule without component f, i.e. applicable to stocks without reliable catch length data. This method is meant as a last resort and should be avoid if possible.

The rb rule is an index adjusted harvest control rule that adjusts the catch advice based on a biomass index but does not have a target. The method is defined as Method 2.3 in the ICES Technical Guidelines [@ICES2022_cat23_tech_guidelines] as

$$
A_{y+1} = A_y \times r \times b \times m
$$

where $A_{y+1}$ is the new catch advice, $A_y$ the previous advice, $r$ the biomass ratio from a biomass index, $b$ a biomass safeguard and $m$ a precautionary multiplier. Furthermore, the change in the advice is restricted by a stability clause that only allows changes of between $+20\%$ and $-30\%$ relative to the previous advice, but the application of the stability clause is conditional on $b=1$ and turned off when $b<1$.

The rb rule should be applied biennially, i.e. the catch advice is valid for two years.

Please note that any change from the default configuration should be supported by case-specific simulations.

## Reference catch $A_y$

The reference catch $A_y$ is usually the most recently advised catch. In a typical ICES setting, an assessment is conducted in an assessment (intermediate) year $y$ to give advice for the following advice year $y+1$, this is the advice for year $y$:

```{r}
### load plaice catch and advice data
data("ple7e_catch")
tail(ple7e_catch)
### get reference catch
A <- A(ple7e_catch, units = "tonnes")
A
```

The ICES Technical Guidelines [@ICES2022_cat23_tech_guidelines] specify that if the realised catch is very different from the advised catch, the reference catch could be replaced by an average of recent catches:

```{r}
### use 3-year average
A(ple7e_catch, units = "tonnes", basis = "average catch", avg_years = 3)
```

The reference catch can also be defined manually:

```{r}
### use 3-year average
A(2000, units = "tonnes")
```

## Biomass index trend (ratio) $r$

The biomass index trend $r$ calculates the trend in the biomass index over last the five years, by dividing the mean of the last two years by the mean of the three preceding years:

$$
r = \Sigma_{i=y-2}^{y-1}(I_i/2)\ / \ \Sigma_{i=y-5}^{y-3}(I_i/3)
$$

The ratio is calculated with the function `r`. Index data should be provided as a `data.frame` with columns `year` and `index`.

```{r, fig.align='center', fig.width=3, fig.height=2, dpi=300, out.width=500}
### load plaice data
data("ple7e_idx")
tail(ple7e_idx)

### calculate biomass trend
r <- r(ple7e_idx, units = "kg/hr")
r

### ICES advice style table
advice(r)

### plot index
### horizontal orange lines correspond to the the 2/3-year averages
plot(r)

### when the value of r is known
r(1)
```

Biomass index data are usually available until the year before the advice year. More recent data can be used and the function automatically picks the most recent data provided to it.

## Biomass safeguard $b$

The biomass safeguard reduces the advice when the biomass index $I$ falls below a threshold value $I_{\text{trigger}}$:

$$
b = \text{min}\{1,\ I_{y-1}/I_{\text{trigger}}\}
$$ By default, the trigger value is defined based on the lowest observed index value $I_{\text{loss}}$ as $I_{\text{trigger}} = 1.4I_{\text{loss}}$.

The biomass safeguard is calculated with the function `b`:

```{r, fig.align='center', fig.width=3, fig.height=2, dpi=300, out.width=500}
### use same plaice data as before
### application in first year with new calculation of Itrigger
b <- b(ple7e_idx, units = "kg/hr")
b

### ICES advice style table
advice(b)

### plot
plot(b)
### plot b and r in one figure
plot(r, b)
```

**Please note that** $I_{\text{trigger}}$ should only be defined once in the first year of the application of the rb rule. In the following years, the same value should be used. For this purpose, `b` allows the direct definition of $I_{\text{trigger}}$, or, more conveniently, $I_{\text{trigger}}$ can be based on the year in which $I_{\text{loss}}$ is defined:

```{r}
### in following years, Itrigger should NOT be updated
### i.e. provide value for Itrigger
b(ple7e_idx, units = "kg/hr", Itrigger = 0.3924721)
### alternatively, the reference year for Iloss can be used
b(ple7e_idx, units = "kg/hr", yr_ref = 2007)
```

## Multiplier $m$

The multiplier $m$ is a tuning parameter and ensures that the catch advice is precautionary in the long term.

The value of $m$ is set to $m=0.5$ for all stocks. The multiplier can be calculated with the function `m()`:

```{r}
m <- m(hcr = "rb")
m

### alias for rb rule
rb_m()

### ICES advice style table
advice(m)
```

Please note that for the rb rule, the multiplier $m$ does lead to a continuous reduction in the catch advice because the rb rule does not include a target.

## Application of rb rule

Now we have all the components of the rb rule and can apply it:

```{r, fig.align='center', fig.width=4, fig.height=3, dpi=300, out.width=500}
advice <- rb(A = A, r = r, b = b, m = m, discard_rate = 27)
advice
```

A discard rate in % can be provided to the argument `discard_rate` and this means the advice is provided for the catch and landings.

The rb rule includes a stability clause that restricts changes to $+20%$ and $-30%$. This stability clause is conditional on the biomass safeguard and is only applied if $b=1$ but turned off when $b<1$.

`cat3advice` can print a table similar to the table presented in ICES advice sheets in which all numbers are rounded following ICES rounding rules:

```{r}
advice(advice)
```

# The chr rule

The chr rule is a (relative) harvest rate-based harvest control rule. The relative harvest rate is defined by dividing the catch by the values from a biomass index. It is not an absolute harvest rate because the absolute stock size is unknown. The method is defined as Method 2.2 in the ICES Technical Guidelines [@ICES2022_cat23_tech_guidelines] as

$$
A_{y+1} = I_{y-1} \times HR_{\text{MSY proxy}} \times b \times m
$$

where $A_{y+1}$ is the new catch advice, $I_{y-1}$ the last biomass index value, $HR_{\text{MSY proxy}}$ the (relative) harvest rate target, $b$ a biomass safeguard and $m$ a precautionary multiplier. Furthermore, the change in the advice is restricted by a stability clause that only allows changes of between $+20\%$ and $-30\%$ relative to the previous advice, but the application of the stability clause is conditional on $b=1$ and turned off when $b<1$.

The chr rule should be applied annually, i.e. the catch advice is valid for one year.

Please note that any change from the default configuration should be supported by case-specific simulations.

## Reference catch $A_y$

A reference catch $A_y$ is needed for the application of the stability clause. The reference catch $A_y$ is usually the most recently advised catch. In a typical ICES setting, an assessment is conducted in an assessment (intermediate) year $y$ to give advice for the following advice year $y+1$, this is the advice for year $y$:

```{r}
### load plaice catch and advice data
data("ple7e_catch")
tail(ple7e_catch)
### get reference catch
A <- A(ple7e_catch, units = "tonnes")
A
```

The ICES Technical Guidelines [@ICES2022_cat23_tech_guidelines] specify that if the realised catch is very different from the advised catch, the reference catch could be replaced by an average of recent catches:

```{r}
### use 3-year average
A(ple7e_catch, units = "tonnes", basis = "average catch", avg_years = 3)
```

The reference catch can also be defined manually:

```{r}
### manual definition
A(2000, units = "tonnes")
```

## Biomass index value $I_{y-1}$

$I_{y-1}$ is the most recent value from the biomass index. This is usually a value from the year before the assessment year ($y$) but a more recent value can be used if available.

Index data should be provided as a `data.frame` with columns `year` and `index`.

```{r, fig.align='center', fig.width=3, fig.height=2, dpi=300, out.width=500}
### load plaice data
data("ple7e_idx")
tail(ple7e_idx)

### get most recent value
i <- I(ple7e_idx, units = "kg/hr")
i

### ICES advice style table
advice(i)

### plot index
plot(i)

### when the value of I is known
I(1)
```

## Biomass safeguard $b$

The biomass safeguard reduces the advice when the biomass index $I$ falls below a threshold value $I_{\text{trigger}}$:

$$
b = \text{min}\{1,\ I_{y-1}/I_{\text{trigger}}\}
$$ By default, the trigger value is defined based on the lowest observed index value $I_{\text{loss}}$ as $I_{\text{trigger}} = 1.4I_{\text{loss}}$.

The biomass safeguard is calculated with the function `b`:

```{r, fig.align='center', fig.width=3, fig.height=2, dpi=300, out.width=500}
### use same plaice data as before
### application in first year with new calculation of Itrigger
b <- b(ple7e_idx, units = "kg/hr")
b

### ICES advice style table
advice(b)

### plot
plot(b)
```

**Please note that** $I_{\text{trigger}}$ should only be defined once in the first year of the application of the chr rule. In the following years, the same value should be used. For this purpose, `b` allows the direct definition of $I_{\text{trigger}}$, or, more conveniently, $I_{\text{trigger}}$ can be based on the year in which $I_{\text{loss}}$ is defined:

```{r}
### in following years, Itrigger should NOT be updated
### i.e. provide value for Itrigger
b(ple7e_idx, units = "kg/hr", Itrigger = 0.3924721)
### alternatively, the reference year for Iloss can be used
b(ple7e_idx, units = "kg/hr", yr_ref = 2007)
```

## Target harvest rate $HR_{\text{MSYproxy}}$

The target harvest rate $HR_{\text{MSYproxy}}$ defines the target for the chr rule and is a proxy for MSY. The standard approach to define $HR_{\text{MSYproxy}}$ is to use catch length data, find years in which the mean catch length $L_{\text{mean}}$ is above a reference length ($L_{F=M}$), calculate the harvest rate for these years, and use their average. The approach is the same as the one used for component $f$ of the rfb rule described above. This needs to be done only once in the first year of the application of the chr rule. In subsequent years, no length data are required.

### Length data

Ideally, length data for several years are provided in a `data.frame` with columns `year`, `length` and `numbers`. An additional column `catch_category` specifying the catch category, such as discards and landings, is optional.

```{r}
data("ple7e_length")
head(ple7e_length)
```

### Length at first capture $L_c$

Only length data above the length at first capture $L_c$ are used to avoid noisy data from fish that are not fully selected. $L_c$ is defined as the first length class where the abundance is at or above 50% of the mode of the length distribution and can be calculated with the function `Lc()`:

```{r, fig.align='center', fig.width=5, fig.height=4, dpi=300, out.width=500, warning=FALSE}
lc <- Lc(ple7e_length)
lc
plot(lc)
```

$L_c$ can change from year to year. Therefore, it is recommended to pool data from several (e.g. 5) years:

```{r, fig.align='center', fig.width=4, fig.height=3, dpi=300, out.width=500, warning=FALSE}
lc <- Lc(ple7e_length, pool = 2017:2021)
lc
plot(lc)
```

If length data are noisy, the size of the length classes can be increased:

```{r, fig.align='center', fig.width=4, fig.height=3, dpi=300, out.width=500, warning=FALSE}
### use 20mm length classes
plot(Lc(ple7e_length, pool = 2017:2021, lstep = 20))
```

Once defined, $L_c$ should be kept constant and the same value used for all data years. $L_c$ should only be changed if there are strong changes in the fishery or data.

### Mean length

After defining $L_c$, the mean (annual) catch length $L_{\text{mean}}$ above $L_c$ can be calculated:

```{r, fig.align='center', fig.width=5, fig.height=3, dpi=300, out.width=700}
### calculate annual mean length
lmean <- Lmean(data = ple7e_length, Lc = lc, units = "mm")
lmean
plot(lmean)
```

If length data are noisy, the size of the length classes can be increased:

```{r, fig.align='center', fig.width=5, fig.height=3, dpi=300, out.width=700}
### use 20mm length classes
plot(Lmean(data = ple7e_length, Lc = lc, units = "mm", lstep = 20))
```

### Reference length

The reference length follows the concepts of @Beverton1957_dynamics_populations and is calculated as derived by @Jardim2015_HCR:

$$
L_{F=M} = 0.75L_c + 0.25L_\infty
$$ where $L_{F=M}$ is the MSY reference length, $L_c$ the length at first capture as defined above, and $L_\infty$ the von Bertalanffy asymptotic length. This simple equation assumes that fishing at $F=M$ can be used as a proxy for MSY and that $M/k=1.5$ (where $M$ is the natural mortality and $k$ the von Bertalanffy individual growth rate). The reference length can be calculated with

```{r}
lref <- Lref(Lc = 264, Linf = 585, units = "mm")
lref
### use a dummy value here for illustrative purposes of the chr rule
lref <- Lref(330, units = "mm")
```

Deviations from the assumptions of $F=M$ and $M/k=1.5$ are possible following Equation A.3 of @Jardim2015_HCR $L_{F=\gamma M, k = \theta M}=(\theta L_\infty + L_c(\gamma + 1))/(\theta + \gamma +1)$ and can be used by providing arguments `gamma` and `theta` to `Lref()`.

### Indicator

The mean catch length relative to the reference $f = L_{\text{mean}}/L_{F=M}$ is used as an indicator. The same function (`f()`) as used in the rfb rule (described above) can be used to calculate the indicator time series:

```{r, fig.align='center', fig.width=4, fig.height=3, dpi=300, out.width=500}
f <- f(Lmean = lmean, Lref = lref, units = "mm")
plot(f)
```

In this case, the mean catch length (orange curve) is above the reference length in two years, indicating that the fishing pressure was below $HR_{\text{MSY proxy}}$ in these two years.

ICES advice sheets typically show the inverse indicator ratio. This can be plotted by adding an `inverse = TRUE` argument to `plot()`:

```{r, fig.align='center', fig.width=4, fig.height=3, dpi=300, out.width=500}
plot(f, inverse = TRUE)
```

The time series of the indicator values can be printed with

```{r, fig.align='center', fig.width=4, fig.height=3, dpi=300, out.width=500}
indicator(f)
```

and the inverse values also with

```{r, fig.align='center', fig.width=4, fig.height=3, dpi=300, out.width=500}
inverse_indicator(f)
```

### Harvest rate

For the estimation of the target harvest rate, the harvest rate needs to be calculated with `HR()`.

```{r, fig.align='center', fig.width=4, fig.height=3, dpi=300, out.width=500}
### catch data
data("ple7e_catch")
### index data
data("ple7e_idx")

### combine catch and index data into single data.frame
df <- merge(ple7e_catch, ple7e_idx, all = TRUE) # combine catch & index data

### calculate harvest rate
hr <- HR(df, units_catch = "tonnes", units_index = "kg/hr")
hr
```

The harvest rate can only be calculated for years in which both catch and index data are available. The harvest rate and its input data can be plotted automatically:

```{r, fig.align='center', fig.width=4, fig.height=3, dpi=300, out.width=700}
plot(hr)
```

### Harvest rate target $HR_{\text{MSYproxy}}$

Now we can use the indicator and (relative) harvest rate time series to calculate the target harvest rate:

```{r, fig.align='center', fig.width=4, fig.height=3, dpi=300, out.width=700}
### calculate (relative) target harvest rate
HR <- F(hr, f) 
HR
plot(HR)
```

The years selected for the target harvest rate are indicated by orange points.

## Multiplier $m$

The multiplier $m$ is a tuning parameter and ensures that the catch advice is precautionary in the long term.

By default, the value of $m$ is set to $m=0.5$. The multiplier can be calculated with the function `m()`:

```{r}
m <- m(hcr = "chr")
m

### alias for chr rule
chr_m()

### ICES advice style table
advice(m)
```

Please note that the multiplier $m$ does not lead to a continuous reduction in the catch advice. The components of the chr rule are multiplicative, this means that $m$ could be considered as adjusting the target harvest rate $HR_{\text{MSY proxy}}$ to $HR'_{\text{MSY proxy}}$:

$$
A_{y+1} = I \times HR_{\text{MSY proxy}} \times b \times m = I \times \frac{HR_{\text{MSY proxy}}}{1/m} \times b = I \times HR'_{\text{MSY proxy}} \times b
$$

## Application of chr rule

Now we have all the components of the chr rule and can apply it:

```{r}
advice <- chr(A = A, I = i, F = HR, b = b, m = m, discard_rate = 27)
advice
```

A discard rate in % can be provided to the argument `discard_rate` and this means the advice is provided for the catch and landings.

The chr rule includes a stability clause that restricts changes to $+20\%$ and $-30\%$. This stability clause is conditional on the biomass safeguard and is only applied if $b=1$ but turned off when $b<1$.

`cat3advice` can print a table similar to the table presented in ICES advice sheets in which all numbers are rounded following ICES rounding rules:

```{r}
advice(advice)
```

# chr rule with custom parameters and discard survival

The ICES plaice stock in the western English Channel was included in the WKBPLAICE [@ICES2025_WKBPLAICE] benchmark in 2024. During this benchmark, a stock-specific management strategy evaluation (MSE) was developed to tune the chr rule. This means the control parameters of the chr rule for this stock differ from the default values. Furthermore, discard survival is also considered. This section of the vignette illustrates how the advice can be calculated with the `cat3advice` R package for this situation.

## Data

Example data is available in the `ple7e_WKBPLAICE` object:

```{r}
data("ple7e_WKBPLAICE")
tail(ple7e_WKBPLAICE)
```

## Reference catch $A_y$

The reference catch $A_y$ is based on the previous catch advice. However, there is an assumption of 50% discard survival, and the reference catch only considers the dead part of the catch advice (landings + dead discards). To calculate this, the data object passed to `A_chr()` needs to include the columns `advice_landings` and `advice_discards`:

```{r}
tail(ple7e_WKBPLAICE[, c("year", "advice_landings", "advice_discards")])

### define discard survival
discard_survival <- 50 ### 50%, from WKBPLAICE 2024

A <- chr_A(ple7e_WKBPLAICE, units = "tonnes", 
           basis = "advice", advice_metric = "catch", 
           discard_survival = discard_survival)
A
```

`A_chr()` calculated the reference catch by adding the landings corresponding to the last advice and 50% of the discards corresponding to the last advice.

## Biomass index value

In this specific case, the biomass index value `I` is set to the average of the values of the last two years (not just the last value):

```{r}
I <- chr_I(ple7e_WKBPLAICE, n_yrs = 2, lag = 1, 
           units = "kg/(hr m beam)")
I
```

The data set has data up to 2024 but the last biomass index value is from 2023. The argument `lag = 1` defines that 2023 is the last year to use. `n_yrs = 2` specifies that the average over two years should be used (2022 and 2023).

## Target harvest rate

Usually, the target harvest rate $HR_\text{MSYproxy}$ is defined by looking at historical catch length data. However, WKBPLAICE conducted an MSE and the values were defined through simulations.

First, we need to calculate the historical harvest rates. Because of the assumed discard survival, the harvest rate only considers the dead catch, i.e. the dead catch divided by the biomass index. The argument `split_discards = TRUE` specifies that the discards should be split into dead and surviving discards, and the argument `discard_survival` defines the survival:

```{r}
### 1st: calculate historical harvest rates
hr <- HR(ple7e_WKBPLAICE, units_catch = "tonnes", 
         units_index = "kg/(hr m beam)", split_discards = TRUE,
         discard_survival = discard_survival)
hr

plot(hr)
```

WKBPLAICE defined $HR_\text{MSYproxy}$ as the average of all historical harvest rates (`yr_ref = 2003:2023`), multiplied by 0.66 (`multiplier = 0.66`):

```{r, warning=FALSE}
### 2nd: calculate harvest rate target
HR <- F(hr, yr_ref = 2003:2023, MSE = TRUE, multiplier = 0.66)
HR

plot(HR)
```

## Biomass safeguard

The biomass safeguard reduces the catch advice when the biomass index $I$ falls below $I_\text{trigger}$. This values was defined by WKBPLAICE as the lowest observed biomass index value (observed in 2007) multiplied by 3.7:

```{r}
b <- chr_b(I, ple7e_WKBPLAICE, units = "kg/(hr m beam)", 
           yr_ref = 2007, w = 3.7)
b

plot(b)
```

## Multiplier

The multiplier is meant to adjust the target harvest rate. However, for this tuned version of the chr rule, the target harvest rate is already a tuned value and the multiplier is set to 1:

```{r}
m <- chr_m(1, MSE = TRUE)
m
```

## Application of the (tuned) chr rule

First, we need to define a discard rate to split the catch advice into landings and discards:

```{r}
### discard rate in % for advice:
### average of 2012 until last year
discard_rate <- ple7e_WKBPLAICE %>%
  dplyr::mutate(discard_rate = discards/catch * 100) %>%
  dplyr::filter(year >= 2012) %>%
  dplyr::summarise(discard_rate = mean(discard_rate, na.rm = TRUE)) %>% 
  as.numeric()
discard_rate
```

The chr rule can then be applied (including discard survival):

```{r}
advice <- chr(A = A, I = I, F = HR, b = b, m = m,
              frequency = "biennial",
              discard_rate = discard_rate,
              discard_survival = discard_survival,
              units = "tonnes", advice_metric = "catch")
advice

advice(advice)
```

The target harvest rate refers to the dead catch (landings + dead part of the discards). Consequently, the advice calculations are based on the dead catch but are topped up to derive the total catch (landings, dead discards, surviving discards). This is done automatically in `chr()` and the advice table illustrates the calculations. The stability clause (limiting changes in the catch advice to +20% and -30%) is also applied to the dead part of the catch advice.

For further details, see the ICES WKBPLAICE report [@ICES2025_WKBPLAICE].

# References
