Modelling_stock_recruitment_with_FLSR.Rmd

---
title: Modelling Stock-Recruitment with FLSR
date: "`r format(Sys.time(), '%d %B, %Y')`"
tags: [FLR]
output:
github_document:
  mathjax: TRUE
pdf_document:
  fig_width: 6 
  fig_height: 4 
license: Creative Commons Attribution-ShareAlike 4.0 International Public License
---

```{r, ini, echo=FALSE, results='hide', message=FALSE, warnings=FALSE, cache=FALSE}
library(knitr)
source("R/ini.R")
```
 
`FLSR` is an [S4](https://stat.ethz.ch/R-manual/R-devel/library/methods/html/Classes_Details.html) class for Stock-Recruitment (SR) models, an extension of `FLModel` , and part of the `FLCore` package.  Commonly used or custom-tailored SR models can be fitted directly to `FLStock` objects, providing estimates of uncertainty. `FLSR` class objects can then be used to visualize the fitted models, calculate biological reference points using `FLBPR`, and perform stock projections. 

## Required packages

To follow this tutorial you should have installed the following packages:

- CRAN: [ggplot2](https://cran.r-project.org/web/packages/ggplot2/index.html)
- FLR: [FLCore](http://www.flr-project.org/FLCore/), [ggplotFL](http://www.flr-project.org/ggplotFL/)

You can do so as follows,

```{r, eval=FALSE}
install.packages(c("ggplot2"))
install.packages(c("FLCore"), repos="http://flr-project.org/R")
install.packages(c("ggplotFL"), repos="http://flr-project.org/R")
```

Initially, the libraries need to be called.

```{r, pkgs}
library(FLCore)
library(ggplotFL)
```

The user can load and visualize the results of an assessment (VPA) already performed and stored in the ple4 `FLStock` object.

```{r, data}
# Load the ple4 FLStock object 
data(ple4)
```

```{r figA, fig.margin=FALSE}
# Plot the assesment output
plot(ple4)
```

# The Stock-Recruitment (SR) relationship

Given that recruitment and spawning stock biomass (SSB) are provided as an output of the assessment, their relationship can be visualized simply by ploting the recruits against the SSB.

```{r figB}
# Plot the SSB-Recruits graph
ggplot(aes(ssb, rec), data=model.frame(FLQuants(ple4, "ssb", "rec"))) +
  geom_point() + geom_smooth(method="loess")
```

## Working with `FLSR` objects
An empty `FLSR` object can be directly created simply by:

```{r, FLSRobject1}
sr1 <- FLSR()
```

An `FLSR` object can be also be created by directly converting an `FLStock` object:

```{r, FLSRobject2}
p4sr <- as.FLSR(ple4)
```

The contents of the `FLSR` object are the following:

```{r, FLSRobject2_contents}
summary(p4sr)
```

In the case of the `ple4` `FLStock` object, recruits are fish of age=1. Hence, the lag between ssb and rec is 1 year. The starting year for SSB is 1957, whereas for recruits it is 1958. 

```{r, rec_ssb_lag}
# Outputs the contents of the first year of the rec and ssb slots of the FLSR object
ssb(p4sr)[,1]
rec(p4sr)[,1]
```

The user can change the recruitment age by triming the `FLStock` object while converting it into an `FLSR` object:

```{r, set_rec_age1}
# You can set a different recruitment age, e.g. age 2, 
# by trimming the FLStock object as follows:
p4sr2 <-as.FLSR(ple4[-1])
```

In this case, the lag between SSB and recruitment is 2 years. The starting year for SSB is 1957, whereas for recruits it is 1959.

```{r, set_rec_age2}
ssb(p4sr2)[,1]
rec(p4sr2)[,1]
```

# Fitting an SR model 

To fit an SR model, a series of commonly-used stock-recruit models are already available, including their corresponding likelihood functions and calculation of initial values. See [SRModels](https://github.com/flr/FLCore/blob/master/R/SRmodels.R) for more details and the exact formulation implemented for each of them. Each method is defined as a function returning a list with one or more elements as follows:

* `modelFormula` for the model, using the slot names (`rec` and `ssb`) to refer to the usual inputs

* `loglFunction` to calculate the loglikelihood of the given model when estimated through Maximum Likelihood Estimation (MLE, see `fmle`)

* `initialFunction` to provide initial values for all parameters in the minimisation algorithms called by `fmle` or `nls`. If required, this function also has two attributes, `lower` and `upper`, that give lower and upper limits for the parameter values, respectively. This is used by some of the methods defined in `optim`, like `"L-BFGS-B"`.

The _model()_ <-- method for `FLModel` can then be called with the value being a list thus described, the name of the function returning such a list, or the function itself. 

The available SR models are: bevholt(), bevholt.ar1(), bevholt.c.a(), bevholt.c.b(), bevholt.d(), bevholt.ndc(), bevholt.sv(), geomean(), logl.ar1(rho, sigma2, obs, hat), ricker(), ricker.ar1(), ricker.c.a(), ricker.c.b(), ricker.d(), ricker.sv(), segreg(), shepherd(), shepherd.ar1(), shepherd.d(), shepherd.d.ar1(), shepherd.ndc(), shepherd.ndc.ar1(), sv2ab(steepness, vbiomass, spr0, model). 

The user can assign e.g. a Ricker SR model to the `FLStock` object. The user can also obtain the model formula of the fitted model, as well as the log-likelihood. The `fmle` method fits the model specified in an `FLModel` object using R's `optim` function to minimise the negative of the log-likelihood function, in the `logl` slot, through calls to the minimisaton routine. The default algorithm for optim is Nelder-Mead; however other options are available (e.g. `"L-BFGS-B"`, see `?optim`).

```{r, fit_SR_model, results="hide"}
# Assign a Ricker SR model and fit it with fmle (which uses logl and R's optim model fitting through MLE)
model(p4sr) <- ricker()
p4sr<-fmle(p4sr)
## model formula
# model(p4sr)
## log-likelihood
# logl(p4sr)
```

The user can extract the initial parameters used by the optimiser, as well as the lower and upper limits of these parameters.

```{r, SR_init_params_lmts}
# initial values for the optimiser
initial(p4sr)
# lower and upper limits for the parameters
lower(p4sr)
upper(p4sr)
```

Diagnostic plots can be produced by simply calling the `plot` function on the `FLSR` object.

```{r, figC}
plot(p4sr)
```

## NS Herring stock-recruitment dataset example

The user can experiment with North Sea herring data where a Ricker model has already been fitted.

```{r, nsher_ricker}
data(nsher)
summary(nsher)
```

The user can change the fitted SR model if so desired. Below bevholt() and cushing() models are used.

```{r, nsher_bh_cs, results="hide"}
# Assign nsher with ricker model to a new object
nsher_ri <- nsher
# change model to bevholt
model(nsher) <- bevholt()
# fit through MLE
nsher_bh <- fmle(nsher)
# change model to cushing
model(nsher) <- cushing()
# fit through MLE
nsher_cs <- fmle(nsher)
```

The three fits can then be inspected visually.

```{r, ri_bh_cs_plots}
plot(nsher_ri)
plot(nsher_bh)
plot(nsher_cs)
```

They can also be compared by using the AIC,

```{r, ri_bh_cs_AIC}
print(paste0('Ricker: ',round(AIC(nsher_ri),4),' ',
             'Beverton-Holt: ',round(AIC(nsher_bh),4),' ',
             'Cushing: ',round(AIC(nsher_cs),4)))
```

or Schwarz's Bayesian Information Criterion.

```{r, ri_bh_cs_BIC}
# this chunk plots the fits from the 3 different SR models
print(paste0('Ricker: ',round(BIC(nsher_ri),4),' ',
             'Beverton-Holt: ',round(BIC(nsher_bh),4),' ',
             'Cushing: ',round(BIC(nsher_cs),4)))
```

Additionally, a profiling of the model parameters can be visualised for each fitted model.

```{r, figD, fig.width=6, fig.height=3}
# Profile the likelihood to check the fit
par(mfrow=c(1,3))
profile(nsher_ri, main="Ricker")
profile(nsher_bh, main="Beverton-Holt")
profile(nsher_cs, main="Cushing")
```

# Advanced topics

Please note: some of the code below is provided for demonstration purposes only, as the used datasets are not necessarily adequate for estimating more than 2 parameters of an SR model. 

SR model parameters can also be fixed. In this case, _steepness_ is fixed to a value of 0.8. Details on the model parameterization can be found in  [SRmodels](https://github.com/flr/FLCore/blob/master/R/SRmodels.R).

```{r, figE1, results="hide"}
# Fit a bevholtSV model with fixed steepness at 0.8
model(p4sr) <- bevholtSV
p4sr <- fmle(p4sr, fixed = list(s = 0.8))
```
```{r, figE2}
# Plot the SR model and show parameters
par(mfrow=c(1,1))
plot(p4sr)
params(p4sr)
```

Custom SR models can be implemented. To define a new model requires the specification of its

1. functional form,
2. likelihood,
3. bounds, and
4. starting values.

For example, the user can fit the Deriso-Schnute model below.

```{r, SR_custom1}
# Define a custom SR model (Deriso Schnute)
dersch<-function(){
  ## log-likelihood
  logl <- function(a,b,c,rec,ssb) {
          res<-loglAR1(log(rec), log(a*ssb*(1-b*c*ssb)^(1/c)))
          return(res)
          }
  ## initial parameter values
  initial <- structure(function(rec, ssb){
      slopeAt0 <- max(quantile(c(rec)/c(ssb), 0.9, na.rm = TRUE))
      maxRec   <- max(quantile(c(rec), 0.75, na.rm = TRUE))
      ### Bevholt by default c=-1
      return(FLPar(a=slopeAt0, b=1/maxRec, c=-1))},
    lower=rep(-Inf, 3),
	  upper=rep( Inf, 3))
  ## model to be fitted
  model  <- rec~a*ssb*(1-b*c*ssb)^(1/c)
  return(list(logl = logl, model = model, initial = initial))}
```
```{r, SR_custom2, results="hide"}
# Fit the custom SR model
model(nsher)<-dersch()
nsher_dersch<-fmle(nsher,fixed=list(c=-1))
```
```{r, SR_custom3}
# Plot the custom SR model
plot(nsher_dersch)
```

An SR model with AR1 autocorrelation can be also be fitted. 
```{r, SR_custom_AR1}
# Define a custom SR AR1 model
rickerAR1 <- function()
  {
  ## log-likelihood
  logl <- function(a, b, rho, rec, ssb)
          loglAR1(log(rec), log(a*ssb*exp(-b*ssb)), rho=rho)
  ## initial parameter values
  initial <- structure(function(rec, ssb) {
      res  <-coefficients(lm(c(log(rec/ssb))~c(ssb)))
      return(FLPar(a=max(exp(res[1])), b=-max(res[2]), rho=0))},
	  lower=rep(-Inf, 3),
	  upper=rep( Inf, 3))
  ## model to be fitted
	model  <- rec~a*ssb*exp(-b*ssb)
	return(list(logl=logl, model=model, initial=initial))}
```
```{r, SR_custom_AR2, results="hide"}
# Fit the custom SR AR1 model
model(nsher)<-rickerAR1()
nsherAR1 <-fmle(nsher)
```
```{r, SR_custom_AR3}
# Plot the custom SR AR1 model
plot(nsherAR1)
```
Finally, an SR model with covariates (e.g. the NAO index) can be used to model environmental effects on the stock recruitment relationship.

```{r, SR_custom_covars1}
# Read in the data to represent the covariate
nao  <-read.table(url("https://www.esrl.noaa.gov/psd/data/correlation/nao.data"), 
       skip=1, nrow=62, na.strings="-99.90")
dnms <-list(quant="nao", year=1948:2009, unit="unique", season=1:12, area="unique")
nao  <-FLQuant(unlist(nao[,-1]), dimnames=dnms, units="nao")

# Include NAO as the covariate (covar) and adjust the model.
# (Note that covar must be an FLQuants with a single component called `covar`
# that matches the year span of the data.)
nsherCovA <- nsher
nsherCovA <- transform(nsherCovA,ssb=ssb/1000,rec=rec/1000)

# Define the custom SR model with covariate
# (modified so temperature affects larval survival)
rickerCovA <- function(){
  ## log likelihood
  logl <- function(a, b, c, rec, ssb, covar){
          loglAR1(log(rec), log(a*(1+c*covar[[1]])*ssb*exp(-b*ssb)))}
  ## initial parameter values
  initial <- structure(function(rec, ssb, covar) {
      res  <-coefficients(lm(c(log(rec/ssb))~c(ssb)))
      return(FLPar(a=max(exp(res[1])), b=-max(res[2]), c=0.0))},
	  lower=rep(-Inf, 3),
	  upper=rep( Inf, 3))
  ## model to be fitted
	model  <- rec~a*(1+c*covar[[1]])*ssb*exp(-b*ssb)
	return(list(logl=logl, model=model, initial=initial))}
```
```{r, SR_custom_covars2, results="hide"}
# Fit the custom SR model with covariate
model(nsherCovA)<-rickerCovA()
covar(nsherCovA)<-FLQuants(covar=seasonMeans(trim(nao, year=dimnames(ssb(nsherCovA))$year)))
nsherCovA       <-fmle(nsherCovA,fixed=list(c=0))
```
```{r, SR_custom_covars3}
# Plot the custom SR model with covariate
plot(nsherCovA)
```

# References
Beverton, R.J.H. and Holt, S.J. (1957) On the dynamics of exploited fish populations. MAFF Fish. Invest., Ser: II 19, 533. ISBN: [1930665946](http://www.worldcat.org/title/on-the-dynamics-of-exploited-fish-populations/oclc/469978238)

Needle, C.L. Recruitment models: diagnosis and prognosis. Reviews in Fish Biology and Fisheries 11: 95-111, 2002. doi: [10.1023/A:1015208017674](https://doi.org/10.1023/A:1015208017674)

Ricker, W.E. (1954) Stock and recruitment. J. Fish. Res. Bd Can. 11, 559-623. doi: [10.1139/f54-039](https://doi.org/10.1139/f54-039)

Shepherd, J.G. (1982) A versatile new stock-recruitment relationship for fisheries and the construction of sustainable yield curves. J. Cons. Int. Explor. Mer 40, 67-75. doi: [10.1093/icesjms/40.1.67](https://doi.org/10.1093/icesjms/40.1.67)


# 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('ggplotFL')`
* ggplot2: `r packageVersion('ggplot2')`
* **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

**Nikolaos NIKOLIOUDAKIS**. Institute of Marine Research (IMR), Pelagic Fish Group, Nordnesgaten 33, P.O. Box 1870, 5817 Bergen, Norway. <http://www.imr.no/>