lbsprHakeTutorial.Rmd

---
title: "Length Based Spawning Potential Ratio (LB-SPR) with Southern Hake"
author: "S. Cerviño"
date: "`r Sys.Date()`"
output:
  html_document:
    number_sections: yes
    toc: yes
    toc_float: yes
  pdf_document:
    number_sections: yes
    toc: yes
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


# Introduction
## What is Spawning Potential ratio (SPR)?
- Spawning Potential ratio (SPR) is defined as the proportion of Spawning Biomass per recruit (SBPR) in an exploited stock with regards to SBPR in an unfished (virgin) stock. 
- The rationale behind is that the abundance at length in the population decreases with ageing (length) because of the mortality (M and F). 
- A virgin population will have a larger amount of large mature individuals than an exploited population. 
- The SPR ranges between 1 (virgin population) and 0. 
  - A SPR in the Range of 0.35-0.4 are ususaly considered a population at MSY level although this is a quite variable parameter. 
  - A population with SPR below 0.1-0.15 are considered collapsed.

## LBSPR Modes
LBSPR package contains functions to run the Length-Based Spawning Potential Ratio (LBSPR) method. The LBSPR package can be used in two ways: 

- **Simulating** the expected length composition, growth curve, and SPR and yield curves using the LBSPR model and 

- **Fitting** to empirical length data to provide an estimate of the spawning potential ratio (SPR).

## LBSPR data needed
The LBSPR method has been developed for data-limited fisheries, Data requested are: 

- A representative sample of the **size structure of the catch**.

- An understanding of the **life history** of the species. 

The LBSPR method does not require knowledge of the natural mortality rate (*M*), but instead uses the ratio of natural mortality and the von Bertalanffy growth coefficient (*K*) (*M*/*K*), which is believed to vary less across stocks and species than *M* (Prince et al. 2015).  

## LBSPR Assumptions
LBSPR model relies on a number of simplifying **assumptions**:

- The LBSPR models are equilibrium based.

- The length composition data is representative of the exploited population at steady state. Whith logistic selection.

See the publications listed in the reference list for full details of the assumptions of the model, including simulation testing to evauate the effect of violations of these assumptions.  

## SPR versions

### Age-Structured Length-Based Model
The LBSPR model described by Hordyk et al. (2015a, b), and tested in a MSE framework (Hordyk et al. 2015c), use a conventional age-structured equilibrium population model.  An important assumption of this model structure is that selectivity is age-based not length-based.

### Length-Structured Growth-Type-Group Model
Hordyk et al. (2016) describe a length-structured version of the LBSPR model that uses growth-type-groups (GTG) to account for size-based selectivity. The GTG-LBSPR model also has the ability to include variable *M* at size (by default *M* is assumed to be constant).  The GTG-LBSPR model typically estimates a lower fishing mortality rate for a given size structure compared to the earlier age-structured model.  This is because the age-structured model has a 'regeneration assumption', where, because of the age-based selectivity assumption, large individuals are expected even at high fishing mortality (large, young fish).

The default setting for the LBSPR package is to use the GTG-LBSPR model for all simulation and estimation. Control options in the simulation and estimation functions can be used to switch to the age-structured LBSPR model. 

Please make sure you understand the data and the biological parameters (and how the model treats these) and critically evaluate any output of the LBSPR model. 

## Installing and loading the Package
The LBSPR package is now available on CRAN:
```{r, eval=FALSE}
install.packages("LBSPR")
```
You can install the development version of the package from GitHub using the `devtools` package:
```{r, eval=FALSE}
install.packages("devtools")
devtools::install_github("AdrianHordyk/LBSPR")
```
Load the Package
```{r}
library(LBSPR) 
```

# Simulation Mode Tutorial
The LBSPR package can be used to generate the expected size composition, the SPR, and relative yield for a given set of biological and exploitation pattern parameters. 

## `LB_pars` Object
The first thing to do is to create a `LB_pars` object that contains all of the required parameters for the simulation model. `LB_pars` is an S4 class object.  

The S4 system is different to the S3 system that is commonly used in R, and that R users are familiar with.  Don't worry if you've never used S4 objects before. The main thing to know is that elements in a S4 object are called **slots** and you access them using the `@` symbol (rather than the `$` symbol that is used with data.frames). 

You can read more about S4 objects [here](http://adv-r.had.co.nz/S4.html).

### Create a new `LB_pars` Object
To create a new `LB_pars` object you use the `new` function:
```{r}
MyPars <- new("LB_pars")
```

You can see the elements or **slots** of the `LB_pars` object using the `slotNames` function:
```{r}
slotNames(MyPars)
```

`MyPars` is an object of class `LB_pars`.  You can access the help file for classes by using the `?` symbol (similar to how you find the help file for functions):
```{r, eval=FALSE}
class?CLASSNAME
class?LB_pars
```

### Populate the `LB_pars` Object
The `LB_pars` object has `r length(slotNames(MyPars))` slots. However, not all parameters need to be specified for the simulation model. 

Some parameters are essential, and a warning message should appear if you attempt to progress without values (please let me know if there are issues). 

Default values will be used for some of the other parameters if no value is specified. For example, the first slot (`r slotNames(MyPars)[1]`) is a character object that can be used for the species name. If this slot is left empty, the simulation model will populate it with a default value.

A message should alert you any time a default value is being used. 

The minimum parameters that are needed for the simulation model are:

**Biology**

 - von Bertalanffy asymptotic length (`Linf`)
 
 - M/K ratio (natural mortality divided by von Bertalanffy K coefficient) (`MK`)
 
 - Length at 50% maturity (`L50`)
 
 - Length at 95% maturity (`L95`)
 
**Exploitation**

 - Length at 50% selectivity (`SL50`)
 
 - Length at 95% selectivity (`SL95`)
 
 - F/M ratio (`FM`) **or** SPR (`SPR`). If you specify both, the F/M value will be ignored. 
 
**Size Classes**

 - Width of the length classes (`BinWidth`)
 
Type `class?LB_pars` in the console for help.
 
To create an example parameter object:

```{r}
MyPars@Linf <- 130 
MyPars@L50 <- 66 
MyPars@L95 <- 70
MyPars@MK <- 1.5 
MyPars@SL50 <- 50 
MyPars@SL95 <- 65
MyPars@SPR <- 0.4
MyPars@BinWidth <- 5
```

## Running the Simulation Model 
Now we are ready to run the LBSPR simulation model. To do this we use the `LBSPRsim` function:
```{r}
MySim <- LBSPRsim(MyPars)
```
You will notice some messages in the console alerting you that default values have been used.  You can change these by specifying values in `MyPars` and re-running the `LBSPRsim` function.

We'll manually set those values here so we don't keep seeing the messages throughout the vignette.

```{r}
MyPars@BinMax <- 150
MyPars@BinMin <- 0
```

We can also choose to set the units for the length parameters:
```{r}
MyPars@L_units <- "mm"
```

### The `LB_obj` Object
The output of the `LBSPRsim` function is an object of class `LB_obj`.  This is another S4 object, and contains all of the information from the `LB_pars` object and the output of the `LBSPRsim` function. 

Many of the functions in the LBSPR package return an object of class `LB_obj`.  You should not modify the `LB_obj` object directly. Rather, make changes to the `LB_pars` object (`MyPars` in this case), and re-run the simulation model (or other functions, covered later in the vignette).

### Simulation Output
Let's take a look at some of the simulated output.

```{r}
MySim@SPR 
```
The simulated SPR is the same as our input value (`MyPars@SPR`).

What is the ratio of fishing mortality to natural mortality in this scenario?

```{r}
MySim@FM 
```

It is important to note that the F/M ratio reported in the LBSPR model refers to the *apical F* over the adult natural mortality rate. That is, the value for fishing mortality refers to the highest level of *F* experienced by any single size class. 

If the selectivity pattern excludes all but the largest individuals from being exploited, it is possible to have a **very** high F/M ratio in a sustainable fishery (high SPR).

A couple of more simulations with alternative values:

Specify F/M instead of SPR:
```{r}
MyPars@SPR <- numeric() # remove value for SPR 
MyPars@FM <- 1 # set value for FM
MySim <- LBSPRsim(MyPars)
round(MySim@SPR, 2) # SPR at F/M = 1 
```

Change the life history parameters:
```{r}
MyPars@MK <- 2.0 
MySim <- LBSPRsim(MyPars)
round(MySim@SPR, 2) # SPR 
MyPars@MK <- 0.5
MySim <- LBSPRsim(MyPars)
round(MySim@SPR, 2) # SPR 
MyPars@Linf <- 120
MySim <- LBSPRsim(MyPars)
round(MySim@SPR, 2) # SPR 
```

Change selectivity parameters:
````{r}
MyPars@MK <- 1.5 
MyPars@SL50 <- 10
MyPars@SL95 <- 15 
MyPars@FM <- 1 
MySim <- LBSPRsim(MyPars)
round(MySim@SPR, 2) # SPR 
MyPars@SL50 <- 80
MyPars@SL95 <- 85 
MySim <- LBSPRsim(MyPars)
round(MySim@SPR, 2) # SPR 
```
### Control Options
There are a number of additional parameters that can be modified to control other aspects of the simulation model. 
For example, by default the LBSPR model using the Growth-Type-Group model (Hordyk et at. 2016).  The `Control` argument can be used to switch to the Age-Structured model (Hordyk et al. 2015a, b):
```{r}
MyPars@Linf <- 100 
MyPars@SL50 <- 50
MyPars@SL95 <- 55 
MyPars@FM <- numeric()
MyPars@SPR <- 0.4
MySim <- LBSPRsim(MyPars, Control=list(modtype="absel"))
MySim@FM 
MySim <- LBSPRsim(MyPars, Control=list(modtype="GTG"))
MySim@FM  # lower F/M for the GTG model
```
See the help file for the `LBSPRsim` function for additional parameters for the `Control` argument.  

### Plotting the Simulation
The `plotSim` function can be used to plot `MySim`:
```{r, fig.width = 6, fig.height=6}
plotSim(MySim) 
```
(Azul tenemos la población virgen, mientras que en rojo tenemos lo que queda después de pescar, no lo que pescamos)
 
By default the function plots: a) the expected (equilibrium) size structure of the catch and the expected unfished size structure of the vulnerable population, b) the maturity and selectivity-at-length curves, c) the von Bertalanffy growth curve with relative age, and d) the SPR and relative yield curves as a function of relative fishing mortality (see note above on the *F/M* ratio).
The `plotSim` function can be controlled in a number of ways.  For example, you can plot the expected unfished and fished size structure of the population by changing the `lf.type` argument:
```{r, fig.width = 6, fig.height=6}
plotSim(MySim, lf.type="pop")
```

Santi: Why fished and unfished populations are different at size below ~50 cm, which is not fishable populationW Same doubt with relative F/M plot with yield and SSB/SSB0 equal to zero at F/M~=2 (Fcrash reference point). Probably there are any kind of compensatory recruitment. There is a default steepness parameter set as `r MyPars@Steepness`. This means that recruitmet is 0.7 times R0 (e.g. asintotic recruitment in Bev-Holt model) when SSB is 0.2 times virging biomass (SSB0).

To check impact of recruitment just compare with a Per Rec model (See ?MyPlot)
```{r, fig.width = 6, fig.height=6}
plotSim(MySim, lf.type="pop", perRec = TRUE)
```

Santi: How to get F/M corresponding with maximun yield? How to access to data outputs?

Individual plots can be created using the `type` argument:
```{r, fig.width = 4, fig.height=4}
plotSim(MySim, type="len.freq")
```
See `?plotSim` for more options for plotting the output of the LBSPR simulation model. 

# Tutorial fitting Empirical Length Data 
Apart of simulate a exploited population, LB-SPR can be used to estimate the exploitation parameters (selection, F/M and SPR) if catch-at-length data (observations or frequency; yearly or time series) are available. Two objects are required to fit the LBSPR model to length data: `LB_pars` which contains the life-history parameters (described above) and `LB_lengths`, which contains the length frequency data.

## Creating a `LB_lengths` object
A `LB_lengths` object can be created in two ways.  The `new` function can be used to create an empty object which can be manually populated:
```{r}
MyLengths <- new("LB_lengths")
slotNames(MyLengths)
```
However, it is probably easier to create the `LB_lengths` object by directly reading in a CSV file.
A number of CSV files containing example data have been included in the LBSPR package.  To find the location of the data files on your machine, use the `DataDir` function:
```{r, eval=FALSE}
datdir <- DataDir()
```

The available example files can be printed out using `list.files`:
```{r}
datdir <- DataDir()
list.files(datdir, pattern=".csv")
```
The length data can be either **raw** data, that is, individual length measurements, or **length frequency** data, where the first column of the data set must contain the mid-points of the length bins, and the remaining columns contain the counts for each length class.
The data type (`freq` or `raw`) must be specified in the call to the `new` function.  
A valid `LB_pars` object must also be supplied to the `new` function.
Following the normal convention in R, if a header row is present in the data file, the argument `header` in the call to `new` must be set to `TRUE`. 
Finally, the full file path must be supplied to the `new` function in order to read in the CSV file.  If the file is not found, or the `file` argument is left empty, a blank `LB_lengths` object with be created.

## Reading in Example CSV 
A valid `LB_pars` object must be first created (see sections above):
```{r}
MyPars <- new("LB_pars")
MyPars@Species <- "MySpecies"
MyPars@Linf <- 100 
MyPars@L50 <- 66 
MyPars@L95 <- 70
MyPars@MK <- 1.5 
MyPars@L_units <- "mm"
```
Note that only the life history parameters need to be specified for the estimation model. The exploitation parameters will be estimated.  
A length frequency data set with multiple years:
```{r}
Len1 <- new("LB_lengths", LB_pars=MyPars, file=paste0(datdir, "/LFreq_MultiYr.csv"), 
            dataType="freq")
```
A length frequency data set with multiple years and a header row (identical to Len1 data, but with a header row):
```{r}
Len2 <- new("LB_lengths", LB_pars=MyPars, file=paste0(datdir, "/LFreq_MultiYrHead.csv"), 
            dataType="freq", header=TRUE)
```
A raw length data set with multiple years:
```{r}
Len3 <- new("LB_lengths", LB_pars=MyPars, file=paste0(datdir, "/LRaw_MultiYr.csv"), 
            dataType="raw")
```
Notice that for raw length measurements you must specify the parameters for the length bins (maximum, minimum, and width of length classes) in the `LB_pars` object.  If these are left blank, default values are used.

## Plotting Length Data
The `plotSize` function can be used to plot the imported length data. This is usually a good idea to do before proceeding with fitting the model, to confirm that everything has been read in correctly:
```{r, fig.width=7, fig.height=4}
plotSize(Len1)
plotSize(Len2)
plotSize(Len3)
```

## Fit the Model
The LBSPR model is fitted using the `LBSPRfit` function:
```{r, echo=FALSE, message=FALSE}
myFit1 <- LBSPRfit(MyPars, Len1)
myFit2 <- LBSPRfit(MyPars, Len2)
```
```{r, eval=FALSE}
myFit1 <- LBSPRfit(MyPars, Len1)
myFit2 <- LBSPRfit(MyPars, Len2)
```
Note that the `Control` argument can be used to modify the additional parameters or LBSPR model type (see description in earlier section).

## Examine and Plot Results
The LBSPR package uses a Kalman filter and the Rauch-Tung-Striebel smoother function (see `FilterSmooth`) to smooth out the multi-year estimates of SPR, F/M, and selectivity parameters. 
The smoother parameter estimates can be accessed from the `myFit` object (which is an object of class `LB_obj` [see earlier section for details]):
```{r}
myFit1@Ests
```
Note that by default the smoothed estimates are used in the plotting routines.
The individual point estimates for each year can be accessed from the `LB_obj` object:
```{r}
data.frame(rawSL50=myFit1@SL50, rawSL95=myFit1@SL95, rawFM=myFit1@FM, rawSPR=myFit1@SPR)
```
The `plotSize` function can also be used to show the model fit to the data:
```{r, fig.width=7, fig.height=4}
plotSize(myFit1)
```
Similarly, the `plotMat` function can be used to show the specified maturity-at-length curve, and the estimated selectivity-at-length curve:
```{r, fig.width=5, fig.height=5}
plotMat(myFit1)
```
Finally, the `plotEsts` function can be used to visually display the estimated parameters. Note that this works for all data sets, but only makes sense when there are several years of data:
```{r, fig.width=8, fig.height=5}
plotEsts(myFit1)
```
By default the plotting function adds the smoother line to the estimated points. 

# Fitting Southern hake data
Southern hake life history. 2 options:

1. ICES model LH: MK=2.3; M = 0.4; Linf=129; L50=33 and L95=50 

2. LH Invariants analysis: MK=1.72; M = 0.28; Linf=100; L50=38 and L95=55 

Populate LB_pars with Sothern hake ICES LH option
```{r}
Hke1Pars <- new("LB_pars")
Hke1Pars@Linf <- 129 
Hke1Pars@L50 <- 33 
Hke1Pars@L95 <- 50
Hke1Pars@MK <- 2.3 
Hke1Pars@M <- 0.4 
Hke1Pars@L_units <- "cm"
```
Populate LB_pars with Sothern hake Invariants LH option
```{r}
Hke2Pars <- new("LB_pars")
Hke2Pars@Linf <- 100 
Hke2Pars@L50 <- 38 
Hke2Pars@L95 <- 55
Hke2Pars@MK <- 1.72 
Hke2Pars@M <- 0.28
Hke2Pars@L_units <- "cm"
```

LB_lengths object read from Southerh hake length frequency cvs file and plot
```{r}
HkeLenFreq <- new("LB_lengths", LB_pars=Hke1Pars, file="LFreqSHake.csv", dataType="freq", header=TRUE)
HkeLenFreq@L_units <- Hke1Pars@L_units
plotSize(HkeLenFreq)
```

Fitting the hake model
```{r, echo=FALSE, message=FALSE}
hke1Fit <- LBSPRfit(Hke1Pars, HkeLenFreq)
```

Plots

Plotting length distribution fit. 
```{r}
plotSize(hke1Fit)
```
Plotting yearly selectivity and maturity
```{r}
plotMat(hke1Fit)
```
Plotting yearly selectivity, F/M and SPR
```{r}
plotEsts(hke1Fit)
```

Fit 1 conclusion

- Good length distribution fit 

- Selectivity parameters are not far away than those estimeted with GADGET at ICES

- However

  - F/M ratio quite high. Around 3 in recent years. F high = 3 * 0.4 = 1.2
  
  - SPR extremely low (around 0.03) when optimum is between 0.35-0.4 and collapsed below 0.1
  
## Lets test LH based on invariants
LB_lengths object read from Southerh hake length frequency cvs file and plot

```{r, message=FALSE}
Hke2LenFreq <- new("LB_lengths", LB_pars=Hke2Pars, file="LFreqSHake.csv", dataType="freq", header=TRUE)
Hke2LenFreq@L_units <- Hke2Pars@L_units
hke2Fit <- LBSPRfit(Hke2Pars, Hke2LenFreq)
```

Fitting length distribution
```{r}
plotSize(hke2Fit)
```

Comparing Model Fit (NegLogLik). Model 1 fits best 2010-11 and worst 2012-18
```{r}
hke1Fit@NLL / hke2Fit@NLL
```

Plotting yearly selectivity, F/M and SPR
```{r}
plotEsts(hke2Fit)
```

## Conclusions:  

- Even with quite different LH the fishery perception and stock status is prety much the same.

- Estimated selection is similar than the ICES estimated selection

- What are the sources of such a hugh differences between LBSPR and GADGET?

  - GADGET considers much more aditional information (e.g. abundance indices trends)
  
  - Check whether length distributions provide different signal than abundance indices.
  
  - Check whether alternative selection (dome shaped) can provide less conflicting signals  


# References

Hordyk, A.R., Ono, K., Sainsbury, K.J., Loneragan, N., and Prince, J.D. 2015a. Some explorations of the life history ratios to describe length composition, spawning-per-recruit, and the spawning potential ratio. ICES J. Mar. Sci. 72: 204 - 216.

Hordyk, A.R., Ono, K., Valencia, S.R., Loneragan, N.R., and Prince, J.D. 2015b. A novel length-based empirical estimation method of spawning potential ratio (SPR), and tests of its performance, for small-scale, data-poor fisheries. ICES J. Mar. Sci. 72: 217 – 231. 

Hordyk, A.R., Loneragan, N.R., and Prince, J.D. 2015c. An evaluation of an iterative harvest strategy for data-poor fisheries using the length-based spawning potential ratio assessment methodology. Fish. Res. 171: 20– 32.

Hordyk, A., Ono, K., Prince, J.D., and Walters, C.J. 2016. A simple length-structured model based on life history ratios and incorporating size-dependent selectivity: application to spawning potential ratios for data-poor stocks. Can. J. Fish. Aquat. Sci. 13: 1– 13. doi: 10.1139/cjfas-2015-0422.

Prince, J.D., Hordyk, A.R., Valencia, S.R., Loneragan, N.R., and Sainsbury, K.J. 2015. Revisiting the concept of Beverton–Holt life-history invariants with the aim of informing data-poor fisheries assessment. ICES J. Mar. Sci. 72: 194 - 203.

# Web resources
https://raw.githubusercontent.com/ices-tools-dev/ICES_MSY/master/R/LBSPR.R

https://raw.githubusercontent.com/ices-tools-dev/ICES_MSY/master/R/LBI.R

https://github.com/AdrianHordyk/LBSPR

https://github.com/merrillrudd/LIME/wiki

http://barefootecologist.com.au/