Simulating-stocks-using-FLR.Rmd

---
title: "Simulating Fish Stocks using FLR"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
  html_document:
    theme: united
    highlight: tango
github_document:
  mathjax: TRUE
pdf_document:
  fig_width: 6 
  fig_height: 4
header-includes:
  \usepackage{underscore}
tags: [FLR]
license: Creative Commons Attribution-ShareAlike 4.0 International Public License
---


```{r, ini, echo=FALSE, results='hide', message=FALSE, warnings=FALSE, cache=FALSE}
library(knitr)
library(ggplotFL)
source("R/ini.R")
options(tinytex.verbose = TRUE)
```

Simulating fish stocks using FLR
================================

Introduction
------------

We will show you how to simulate stocks from scratch using the libraries [FLife](http://www.flr-project.org/FLife/) and [FLasher](http://www.flr-project.org/FLasher/) from the [Fisheries Library in R](http://www.flr-project.org/). The combination of those two libraries allows the user to initiate and project a fishstock based on the understanding of life-history theory *(quote the big cahuna)*. This tutorial is structured in two parts; the first part shows the quick way and the second part highlights the flexibility of the method.


We will use `FLife` to setup the stock and `FLasher` to project the stock. The process can be summarised as follows:

1. Setting life-history parameters for stock.
    * at least one parameter needs to be set: $L_\infty$. 
    * rest of parameters is calculated using *big cahuna*'s relationship.
    * creation of a simulation-setup object.
2. Calculation of equilibrium states of fish stocks using life-history parameters.
    * the simulation-setup object is used to calculate stock parameters at a range of fisheries mortalities.
3. Forward projections, using a single stock at equilibrium. 
    * fisheries trajectories are defined
    * Stock-recruitment functions are defined.
    * Deviances are added.

# Simple simulation of a model with default parameters

First we need to load both libraries:

```{r, echo=TRUE, results='hide', message=FALSE }
library(FLife)
library(FLasher)
```

## Setting life-history parameters

To initialise a stock is to use the function `FLPar()` specifying the argument `linf` to create a *FLPar*-object for a given $L_\infty$ (asymptotic length parameter of the von Bertalanffy growth curve). The function `lhPar()` will extend the required attributes using the *big cahuna*'s life history relationships. The resulting *FLPar*-object contains the life history parameters required to create a *FLBRP*-object using `lhEql()`.

```{r, echo = TRUE, results = 'show'}
(par <- FLPar(linf = 100))

(qpar <- lhPar(par))
```

Now the simulation-setup object needs to be created using `lhEql()`, which upon closer inspection contains all the "at-age" information required to create the stock in the next step. In this particular case the stock model will use 41 age classes. The fisheries mortalities (`F_bar =  seq(0, F_crash, length.out = 101)`) are stored as "years" in the object, which is confusing but should not bother the user. The next step will use each fisheries mortality separately and calculate a stock-size at equilibrium for that specific *F*. The overall selectivity is a combination of two quantities: `landings.sel` and `discards.sel`. They are fixed at this stage (to default if not specified) and will be used for the entire projection. Means to change the levels of `F_bar` and the selectivity during the simulation will be shown further down.

```{r, echo = TRUE}
eql <- lhEql(qpar)

summary(qpar)

```

## Calculation of equilibrium states of fish stocks using life-history parameters

This step is very simple. All it it requires is the following little code:

```{r, echo = TRUE}
stk <- as(eql, "FLStock")
```

At this point the reader is invited to browse the *FLStock_*-object at his leisure using `str(stk)`. We will show the output of `summary(stk)` as it is all that is required to highlight what has happened:

```{r summarystk}
summary(stk)
```

In contrast to the `eql` object, the new stock object has 101 "yearly" entries (read stock levels at 101 different `Fbar` levels). The quantities, that are calculated for each age class, have a "41" at the first position (*e.g.* `catch.n`) since we specified 41 age classes. The compounded quantities contain a "1" at the first position (*e.g.* `catch`) as they are applicable to or the result of all classes. The plot of the stock object is a good summary:

```{r, plotstk, fig.cap="Line plots from our stock object."}
plot(stk)
```

The x-axis is an index of the *`Fbar`* sequence specified in `eql` and **not** of time. As such the figure should be read in columns: for a given *`Fbar`* the catch, SSB and recruitment will be in equilibrium at a the level indicated on the plot. The creation of the stock-object also changed *`Fbar`* to `Fage`* via a selection-at-age curve. 
  
  
This selectivity curve is not available separately in this particular object (*FLStock*) and is implied in the distribution of `Fbar`*. This is important as we see later: in order to create an unfished stock, one needs to use a stock with a tiny amount of fisheries mortality (*i.e. $F_{_{bar} \le 0.1}$) to maintain the selectivity pattern within the stock-object. As we will see later, we **can** forward project using zero fisheries mortality. However, if we choose the unfished parameter space generated within `stk` we will not be able to introduce fisheries mortality during the projection stage as the selectivity for all age groups will be 0.

## Forward projections, using a single stock at equilibrium.

In this section we will proiject/simulate the stock for 100 years using the previously generated *FLStock*-object (`stk`). The first step is to create an output object into which we will project the stock. As mentioned before, we must use a stock at equilibrium with a minimum of fisheries mortality in order to maintain the selection structure. As the index of `stk` are linked to fisheries mortalities, we will use the second realisation `stk[,2]` (as the first is `Fbar` = 0, which can be checked by `fbar(stk[,1])`). 

```{r}
omStore <- fwdWindow(stk[,2], eql, end = 100)
summary(omStore)

```

In this storage object are also the relevant parameters for each year and age class. Some of the parameters can be changed prior to the projections as we will see later. In order to project we need to specify two more functions/paramter: the stock-recruitment relationship and the trajectory of *`Fbar`*. As this is the simplest of projections, we will use the default Beverton-Holt stock-recruitment function, parameterised according to the *big cahuna*'s relationships:

```{r, fig.cap = "The Beverton-Holt stock-recruitment model used to simulate the stocks."}
bhm <- as(eql, "predictModel")
nssb <- FLQuant(seq(1, 2000, length=100))

ggplot(model.frame(FLQuants(ssb=nssb, bevholt=predict(bhm, ssb=nssb))),
  aes(x=ssb, y=bevholt)) + geom_line()
```

As for the fishing mortality, we shall fix it at *$F_{MSY}$* for the entire run. 

```{r}
FMSY <- c(fmsy(eql))
hfc <- fwdControl(year=3:100, quant="f", value=FMSY)
```

Note we start at year 3 rather than year 1 or 2. This is due to the internal implementation; when we used the crutch replacing years with levels of `Fbar` and then picked the second `Fbar`, internally the first two slots are filled. This is not a big issue as long we remember to make the time series long enough and cut off the unwanted early years. 

Now we are ready to project a stock, in a deterministic fashion as no noise has been added.

```{r, fig.cap = "Fish-stock behaviour over time. Note the deterministic nature as well as the artefact at the first data point. The artefact is a result of the way the simulation is initialised. See further down how to remove it."}
OpMod <- fwd(omStore, sr = bhm, control = hfc)

plot(OpMod)
```


# Doing stuff to stocks

## Changing life history parameters

We will show how to change and implement different life history parameters during the setup (*i.e.* during the first step). As before we use `FLPar` to set the initial life history parameters. 

The parameters available in a `FLPar` object are:

|Parameter|Description      |
|:--------|:---------------:|
|linf     |$L_\infty$ von Bertalanffy|
|k        |k von Bertalanffy         |
|t0       |size at $t_0$ von Bertalanffy   |
|a        |length-weight relationship|
|b        |length-weight relationship|
|ato96    |age @ 95% maturity \*   |
|a50      |age @ 50% maturity    |
|asym     |no idea          |
|bg       |no idea          |
|m1       |natural mortality parameter for function|
|m2       |same (but not equal) as m1|
|a1       |no idea          |
|sl       |selectivity \@ age \**         |
|sr       |Stock Recruitment Function    |
|s        |Steepnes Bev-Holt |
|v        |Virgin Biomass   |

\* age at 95% maturity is used as an offset to the age at which 50% maturity occurs.

\** selectivity at age parameter, &sigma;<sup>2</sup> of lefthand limb of double normal distribution

```{r}
par <- FLPar(linf=90, a=0.00001, sl=1, sr=2000, a1=4, s=0.65, v=1000)
parg <- lhPar(par)

range <- c(min=1, max=10, minfbar =3, maxfbar = 6, plusgroup=10)
eql <- lhEql(parg, range=range)
```

In this case we changed some life history parameters with FLPar and estimated the other required parameters using lhPar. In the following step, we restricted our age groups to 10 in total, starting at age 1 not 0. The plusgroup is in the 10th age group. If we do change the number of age groups, we also need to specify the minimum and maximum fully selected age classes using minfbar and maxfbar. Once this is all set, we use lhEql to create the simulation-setup object. For the simulation and projection, the syntax will be exactly the same as before:

```{r, fig.cap = "Stock with different life history parameters and the first-year artefact removed."}
stk <- as(eql, "FLStock")
omStore <- fwdWindow(stk[,2], eql, end=100)
bhm <- as(eql, "predictModel")

FMSY <- c(fmsy(eql))
hfc <- fwdControl(year=3:100, quant="f", value = FMSY)
OpMod <- fwd(omStore, sr = bhm, control = hfc)
plot(OpMod[,-1,,,,]) #removing the data artefact
```

```{r, echo = FALSE, results = 'hide', message = FALSE }
rm(list=ls())
library(FLife)
library(FLasher)
library(knitr)
source("R/ini.R")
```

The following chunk of code needs to be run before the code below will work.

```{r}
par <- FLPar(linf=90, a=0.00001, sl=1, sr=2000, a1=4, s=0.65, v=1000)
parg <- lhPar(par)

range <- c(min=1, max=10, minfbar =3, maxfbar = 6, plusgroup=10)
eql <- lhEql(parg, range=range)
stk <- as(eql, "FLStock")
omStore <- fwdWindow(stk[,2], eql, end=100)
bhm <- as(eql, "predictModel")

FMSY <- c(fmsy(eql))
hfc <- fwdControl(year=3:100, quant="f", value=FMSY)

```

## Changing natural mortality (*M*)

Natural mortality (M) is not that easily specified. Due to the way the FLBRP-object is coded, we must write a function that specifies natural mortality by either (a) age, (b) length, (c) weight. There are two default funtions *gislason* (length) and *lorenzen* (weight). Creating an `FLBRP`-object using the default `lorenzen` method:

```
eql <- lhEql(parg, range=range, m = lorenzen)
```
We are not obligated to use the coded functions, as long as we respect a few simple rules. The function must have two arguments, of which one must be `params`. The second argument must be one of the following three: *age*, *length*, *weight*. Although `params` **must** be included, there is no necessity to use it within the function. Here are two examples:
```{r, eval = FALSE}
# length based function with the use of params
mJensen <- function(length, params){
  length[]=params["a50"]
  1.45 / length
}
eql <- lhEql(parg, range=range, m = mJensen)

# age based equation without use of params
mFun <- function(age, params){
 (0.2+ 1.64* exp(-age))
}
eql <- lhEql(parg, range=range, m = mFun)
```

## Changing fisheries mortality (*F*)

The fisheries mortality is changed by specifying different values in the `fwdControl()` function (it is assumed `eql` has already been generated):

```{r, fig.cap = "Fish stock with increasing fisheries mortality to a maximum of 2 x FMSY."}
ifc <- fwdControl(year=3:100, quant="f", value=c(((2*FMSY)/100) * c(3:100)))

OpMod <- fwd(omStore, sr = bhm, control=ifc)
plot(OpMod[,-1,,,,])
```


The value that is changed is *`Fbar`*, the selectivity is not affected by this. Your imagination is the limit when it comes to set fisheries mortality levels. A pattern often used is the "rollercoaster" - up, stable, down:

```{r, fig.cap = "Fish stock simulated using a \"rollercoaster\" fishing behviour. The fishing mortality increases to two times FMSY, remains at that level and then reduces to 1.1 x FMSY."}
rfc <- fwdControl(year=3:100, quant="f", value=c(rep(FMSY,50), (FMSY + ((FMSY/10)*c(1:10))), rep(2*FMSY, 20), (2*FMSY - (((1.1 * FMSY)/10)*c(1:10))), rep(FMSY, 8)))

OpMod <- fwd(omStore, sr = bhm, control=rfc)
plot(OpMod[,-1,,,,])
```

## Selectivity

Selectivity can be set by altering the `harvest` attribute of the output object known as `omStore` in our case:

```{r, fig.cap = "Shift in selectivity at age."}
omStore <- fwdWindow(stk[,2], eql, end=100)
harvest(omStore) <-  c(0.1, 0.5, 1, 0.8, 0.65, 0.55, rep(0.5,4))
harvest(omStore)[,c(75:99)] <- c(0,0, 0.1, 0.5, 1, 0.8, 0.65, 0.55, rep(0.5,2))

moo <- as.data.frame(harvest(omStore))
plot(data~age, moo, type = "n", ylab = "selectivity")
for (i in c(4,77)) lines(data~age, moo[moo$year == i,], col = i, lwd = 2)
legend("bottomright", legend = c("selectivity before shift", "selectivity after shift"), col = c(4, 77), lwd = 2)
```

```{r, fig.cap = "Fish stock with shift in selectivity at year 75."}
OpMod <- fwd(omStore, sr = bhm, control=rfc)
plot(OpMod[,-1,,,,])
```



# Adding stochasticity to the runs

We are showing the methodology to add uncertainty in the recruitment process as well as the fishing mortality. To add uncertainty in natural mortality, I assume that you would could do this by modifying the function above, but it has not been tried **THERE BE DRAGONS** (and let us know how you got on).

In both cases the deviances have to be sampled outside of the FLxxx-process and then added to the parameters. For the recruitment process there is a handy argument within `fwd()`, called `deviances`. There are two ways of introducing deviations from the recruitment curve (we are assuming that we created a omStore object already):

```{r, fig.cap = "Added stochasticity in recruitment."}
lndev03 <- rlnorm(1, FLQuant(0, dimnames=list(year=3:100)), 0.3)
OpMod <- fwd(omStore, sr=bhm, control=hfc, deviances=lndev03)

## or ##

lndev03 <- rlnorm(100, 0, 0.7)

OpMod <- fwd(omStore, sr=bhm, control=hfc, deviances=FLQuant(lndev03, dimnames=list(year=3:100)))

plot(OpMod[,-1,,,,])
```



This can be used to produce any kind of deviations from the stock-recruitment function, e.g. pulse recruitment:

```{r, fig.cap = "Stochastic pulse recruitment events."}
nrPu <- 12
pDev <- rep(1, (100 -2))

#setting the years of the pulses around a mean number of years
s1 <-cumsum(rlnorm(nrPu, 0, .3)*(100/(nrPu+3)))
s1a <- round(s1)[s1 <100]

# setting the pulses
pDev[s1a] <- 6*rlnorm(length(s1a),0,0.6)
puMod <- fwd(omStore, sr=bhm, control=hfc, deviances=FLQuant(pDev, dimnames = list(year = 3:100)))

plot(puMod[,-1,,,,])
```

Adding noise to fisheries mortality is done by altering the control object that controls the fishing pattern. We use a similar method as described above, except we add noise to the average fishing patter:

```{r, fig.cap = "Stochasticity in fishing mortality"}
FMSY <- c(fmsy(eql))

ranfc <- fwdControl(year=3:100, quant="f", value= 1.5 * FMSY *rlnorm(98,0,0.2))

OpMod <- fwd(omStore, sr = bhm, control=ranfc)

plot(OpMod[,-1,,,,])
```

## Bringing it all together

A stock with a lot of stuff going on, see if you can identify it from the code:

```{r}

set.seed(5432)

par <- FLPar(linf=90, a=0.00001, sl=1, sr=2000, a1=4, s=0.65, v=1000)
parg <- lhPar(par)

range <- c(min=1, max=10, minfbar =3, maxfbar = 6, plusgroup=10)
mFun <- function(age, params){
 (0.2+ 1.64* exp(-age))
}
eql <- lhEql(parg, range=range, m = mFun)

stk <- as(eql, "FLStock")
bhm <- as(eql, "predictModel")

omStore <- fwdWindow(stk[,2], eql, end=100)


harvest(omStore) <-  c(0.1, 0.5, 1, 0.8, 0.65, 0.55, rep(0.5,4))
harvest(omStore)[,c(75:99)] <- c(0,0, 0.1, 0.5, 1, 0.8, 0.65, 0.55, rep(0.5,2))

nrPu <- 12
pDev <- rep(1, (100 -2))

#setting the years of the pulses around a mean number of years
s1 <-cumsum(rlnorm(nrPu, 0, .3)*(100/(nrPu+3)))
s1a <- round(s1)[s1 <100]

# setting the pulses
pDev[s1a] <- 6*rlnorm(length(s1a),0,0.6)

moo <- c(rep(FMSY,50), (FMSY + ((FMSY/10)*c(1:10))), rep(2*FMSY, 20), (2*FMSY - (((1.1 * FMSY)/10)*c(1:10))), rep(FMSY, 8))
ranfc <- fwdControl(year=3:100, quant="f", value= moo * rlnorm(98,0,0.2))


OpMod <- fwd(omStore, sr = bhm, deviances=FLQuant(pDev, dimnames = list(year = 3:100)), control=ranfc)

#removing the first year
plot(OpMod[,-1,,,,])
```

# More information

* You can submit bug reports, questions or suggestions on this tutorial at <https://github.com/flr/doc/issues>.
* Or send a pull request to <https://github.com/flr/doc/>
* For more information on the FLR Project for Quantitative Fisheries Science in R, visit the FLR webpage, <http://flr-project.org>.

## Software Versions

* `r version$version.string`
* FLCore: `r packageVersion('FLCore')`
* ggplotFL: `r packageVersion('FLasher')`
* ggplot2: `r packageVersion('FLife')`
* **Compiled**: `r date()`

## License

This document is licensed under the [Creative Commons Attribution-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-sa/4.0) license.

## Author information

**Christoph KONRAD**. European Commission Joint Research Centre (JRC), Ispra VA, Italy. <https://ec.europa.eu/jrc/>

**Iago MOSQUEIRA**. Wageningen Marine Research, Ijmuiden, Netherlands.<https://www.wur.nl/en/Research-Results/Research-Institutes/marine-research.htm>