Tutorial_4.Rmd

---
title: "Rfishpop (Fourth tutorial)"
author: "Marta Cousido Rocha, Santiago Cerviño López, Maria Grazia Pennino"
date: "29/7/2020"
output: html_document
---

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

## Population dynamic tools in support of fisheries managment

As we mentioned in the previous tutorials the aim of the package is to implement a completed MSE (management strategy evaluation) cycle.

The firsts steps of this cycle have been explained in the previous tutorials:

- The operating model  used to generate "true" ecosystem dynamics including the
natural variations in the system.

- The sampling procedure from the operating model to mimic collection of fishery dependent data and research surveys (and their inherent variability).

After these steps the next point is to pass these data to the **assessment model**. Based on this assessment and the HCR (Haverst Control Rule), a
management action is determined (e.g., a change in the TAC).

<span style="color:red">NOTE: Our package does not implement any assessment models, the idea is to use available implementations
of the assessment models. The package contains specific functions to change the format of the data into the required format
of the assessment model function. Now, the package contains such functions for the data-poor methods, LBI (Length Based Indicators) and LB-SPR (Length Based Spawning Potential Ratio). The list of functions will be expanded when exploring the application of other assessment models to the data reported by our package. THE PACKAGE IS AN OPEN PROJECT AND NEW FUNCTIONS ARE INCLUDED CONSTANTLY. </span> 

<span style="color:red">WARNING: Due to the open ISSUES than the user can see in the corresponding page we recommend to be careful with the interpretation of the results corresponding with the application of the assessment models.</span> 

In this tutorial we focus on explaining the use of the <span style="color:blue">functions for passsing the data to the LBI and LB-SPR methods.</span> 
Furthermore, the package also contains a <span style="color:blue">function to define and fill a FLStock object</span>  using the information provided by the main function Population.Modeling in order to use any of the tools developed in the FLR project. Therefore the tutorial finishes explaining this function.


### <span style="color:blue"> Data for Length Based Indicators (LBI)</span> 


The function provides required information for computing Length Based Indicators: Length distributions of catches and the corresponding average weight per length.


The arguments of this function are described above:

- *Pop.Mod*	 A list containing the components returned by Population.Modeling function (main function).

- *CV*	The coefficient of variation associated to the log-normal distribution used in Distribution.length function (see Details of such function).

```{r, eval=FALSE}
??Distribution.length
```

- *RF.value* The number of values generated for each age (given a year and an iteration) from the log-normal distribution used in Distribution.length function (see Details of such function). By default RF.value=1000.

The function **reports** the length distributions of catches for each year and iteration in our Pop.Mod object. Furthermore, the corresponding average weight per length (for each year and iteration) is also provided.

The first point to illustrate the use of this function is to create a population. In this case as we can see below logistic selectivity and Beverton-Holt recruitment model are used.

```{r}
library(Rfishpop)
ctrPop<-list(years=seq(1980,2020,1),niter=1,N0=15000,ages=0:15,minFage=2,
maxFage=5,tc=0.5,seed=NULL)
number_ages<-length(ctrPop$ages);number_years<-length(ctrPop$years)
Mvec=c(1,0.6,0.5,0.4,0.35,0.35,0.3,rep(0.3,9))
M<-matrix(rep(Mvec,number_years),ncol = number_years)
colnames(M)<-ctrPop$years
rownames(M)<-ctrPop$ages
ctrBio<-list(M=M,CV_M=0, L_inf=20, t0=0, k=0.3, CV_L=0, CV_LC=0, a=6*10^(-6), b=3,
             a50_Mat=4, ad_Mat=-0.2,CV_Mat=0)
ctrSEL<-list(type="Logistic", par=list(a50_Sel=2.3, ad_Sel=-0.2),CV_SEL=0)
f=matrix(rep(0.5,number_years),ncol=number_years,nrow=1,byrow=TRUE)
ctrFish<-list(f=f,ctrSEL=ctrSEL)
a_BH=15000; b_BH=50; CV_REC_BH=0
SR<-list(type="BH",par=c(a_BH,b_BH,CV_REC_BH))
Pop.Mod<-Population.Modeling(ctrPop=ctrPop,ctrBio=ctrBio,ctrFish=ctrFish,SR=SR)
```

Then *Data.to.LBI* function is used to obtain the length distributions and average weight per length. NOTE that in the next lines of code we access to the first element of the list which corresponds to the first iteration (in this case we only have one since niter=1, in other examples the next iterations are reported in the next elements of the list).

```{r}
resul=Data.to.LBI(Pop.Mod,CV=0.2)
freq=resul$length[[1]];head(freq)
wal=resul$weight[[1]];head(wal)
```
Furthermore than the data provided by Data.to.LBI function the LBI method also needs values for some life history parameters. Below, we explain how to access to some of them.


- Parameters of Von Bertalanffy Growth Model.
```{r}
L_inf=Pop.Mod$Info$ctrBio$L_inf 
k=Pop.Mod$Info$ctrBio$k
t0=Pop.Mod$Info$ctrBio$t0
```

See details of the parameters using the following code.

```{r,eval=FALSE}
?Length_VB
```

- Now, we can use the Von Bertalanffy Growth Model to obtain the length at 50% maturity (L50=$L_{50}$).

```{r}
x50=Pop.Mod$Info$ctrBio$a50_Mat
L50=Length_VB(L_inf, k, x50, t0) 
```

Since M is constant through the year but age dependent, the natural mortality divided by von Bertalanffy k coefficient is computed as follows, using the mean of the age vector.

```{r}
M.vec=Pop.Mod$Matrices$M[,1,1]
MK <-mean(M.vec)/k
```

Finally, running the following line of code, after load the required code and libraries, LBI are computed.

```{r, include=FALSE,message=FALSE, warning=FALSE}
library(LBSPR) # se usan algunas funciones gráficas
library(reshape2)
library(ggplot2) 
library(tidyr)
#library(ReporteRs) # para generar tablas y Documento-Resumen (no disponible en CRAN) 

# require (rJava)
# .jinit()
# .jcall('java.lang.System','S','getProperty','java.version')
# [1] "1.8.0_211"
#devtools::install_github('davidgohel/ReporteRsjars')
#devtools::install_github('davidgohel/ReporteRs')

source("https://raw.githubusercontent.com/ices-tools-dev/LBIndicator_shiny/master/utilities.R") # incluye la opción de m_k
source("utilities_vpaz.R") 
```

```{r}
LBI=lb_ind(data=freq,binwidth=3,linf=L_inf,lmat=L50,mk_ratio=MK,weight=wal)
LBI
```


```{r}

lb_plot(data=freq,
       binwidth=3,
       linf=L_inf,
       lmat=L50,
       mk_ratio=MK,
       weight=wal,l_units="cm")
```



### <span style="color:blue"> Length Based Spawning Potential Ratio (LB-SPR)</span> 

The function provides required information for computing Length Based Spawning Potential Ratio (LB-SPR): Length distributions of catches.


The **arguments** of this function are described above:

- *Pop.Mod*	A list containing the components returned by Population.Modeling function (main function).

- *CV*	The coefficient of variation associated to the log-normal distribution used in Distribution.length function (see Details of such function).

Remember:

```{r, eval=FALSE}
??Distribution.length
```

- *RF.value* The number of values generated for each age (given a year and an iteration) from the log-normal distribution used in Distribution.length function (see Details of such function). By default RF.value=1000.


The function **reports** the length distributions of catches for each year and iteration in our Pop.Mod object.


<span style="color:red"> THE FUNCTION PROVIDES THE LENGTH DISTRIBUTION OF THE CATCHES AS THE PREVIOUS FUNCTION Data.to.LBI. THE DIFFERENCE IS THAT Data.to.LB.SPR DOES NOT PROVIDE THE AVERAGE WEIGHT PER LENGTH (not necessary here). ON THE OTHER HAND IN THE FOLLOWING EXAMPLE WE PROVIDE THE LINES OF CODE TO OBTAIN THE FORMAT REQUIRED BY THE LB-SPR IMPLEMENTATION. THESE LINES OF CODE ARE NOT INCLUDED IN THE FUNCTION BECAUSE TO CREATE THE REQUIRED OBJECT WE NEED TO FIX THE LIFE PARAMETERS AS WELL, AND THERE ARE A HUGE LIST OF PARAMETERS WHICH CAN BE FIXED AND INFINITY POSSIBILITIES OF DEFINITION. THEN IT IS EASIER AND MORE FLEXIBLE IF WE PROVIDE THE REQUIRED LINES OF CODE OUTSIDE THE FUNCTION. THE USER ONLY NEEDS TO FIX THE VALUES OF THE PARAMETERS AND RUN THE LINES.</span>


To illustrate the practical use of the function we can use the population simulated previously.

```{r}
resul=Data.to.LB.SPR(Pop.Mod,CV=0.2)
```

Some life history parameters have been defined previously in the LBI section. The length at 95% maturity (L95=$L_{95}$) is fixed below. Note that first we use the logistic equation to obtain the age at 95% maturity, and then $L_{95}$.

```{r}
xd=Pop.Mod$Info$ctrBio$ad_Mat
x95=(-xd)*log(19)+x50
L95=Length_VB(L_inf, k, x95, t0)
```

All the information collected above must be used through the following lines of code to apply the LB-SPR model.

First of all we create the object LB_pars.

```{r}
MyPars <- new("LB_pars")
MyPars@Species <- "MySpecies"
MyPars@Linf <- L_inf # von Bertalanffy asymptotic length
MyPars@L50 <-L50     # Length at 50% maturity (L50)
MyPars@L95 <-L95     # Length at 95% maturity (L95)
MyPars@MK <-MK       # The natural mortality divided by von Bertalanffy k #coefficient
```

We need to add such information to the length distribution as follows.

```{r}
freq=resul[[1]]
write.csv(freq, file="len.csv")
Len <- new("LB_lengths", LB_pars=MyPars, file=paste0("len.csv"),
dataType="freq", header=TRUE)
```

Now, the distribution length has the required format to be
introduced in the LBSPRfit().

```{r}
myFit<- LBSPRfit(MyPars,Len)
plotEsts(myFit)
```



### <span style="color:blue"> A FLStock class object</span> 

The function defines and fills a FLStock object using the information provided by the main function Population.Modeling.


The unique **argument** of this function is *Pop.Mod*.

The function **returns** the FLStock object (see R package FLCore) corresponding to the simulated population generated by Population.Modeling function.


If we use the population simulated before the next line allows us to apply the function.

```{r}
resul_FLStock=FLStock.from.Rfishpop(Pop.Mod)
resul_FLStock
```

<span style="color:red"> 
The advantage of being able to transfer the information from Pop.Mod object to a FLStock object is that now we can use tools developed in the FLR project.</span> 


For example, we can apply code/tools of the tutorial "Modelling stock recruitment with FLSR" of FLR project to our population created using Rfishpop package.

Initially, the libraries need to be called.

```{r,message=FALSE,warning=FALSE}
library(FLCore)
library(ggplotFL)
```

The user can load and visualize the results.

```{r}
plot(resul_FLStock)
```


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}
ggplot(aes(ssb, rec), data=model.frame(FLQuants(resul_FLStock, "ssb", "rec"))) +
  geom_point() + geom_smooth(method="loess")
```


An empty FLSR object can be directly created simply by:

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

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

```{r}
p4sr <- as.FLSR(resul_FLStock)
```

The contents of the FLSR object are the following:

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


Fitting an SR model.


```{r,message=FALSE,warning=FALSE}
model(p4sr) <-bevholt()
p4sr<-fmle(p4sr)
```

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

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



Of course this is only one example but any tools of FLR project can be used. Below we apply also some code of "Short Term Forecasting for advice using FLash" tutorial of FLR project.

Load the required packages.

```{r}
library(FLAssess)
library(FLash)
library(ggplotFL)
library(FLBRP)
```

Extending the stock object for the projections. A simple 3-year forecast for the weights, natural mortality, etc., assuming these are equal to their averages over the last 3 years.

```{r}
maxyr_stk <- range(resul_FLStock)[["maxyear"]]
ple4_stf <- stf(resul_FLStock,nyears=3,wts.nyears=3, na.rm=TRUE)
maxyr_stf <- range(ple4_stf)[["maxyear"]]
range(ple4_stf)
stock.wt(ple4_stf)[,ac((maxyr_stf-5):maxyr_stf)]
```

Notice that the future fishing mortality has also been set (average of the last 3 years, by default).

```{r}
ggplot(harvest(ple4_stf)[,ac((maxyr_stf-5):maxyr_stf)]) + geom_line(aes(x=age, y=data)) + facet_wrap(~year)
```

The stock numbers at age and the catch numbers at age are not forecast yet - this is what the fwd() function will perform later.

```{r}
stock.n(ple4_stf)[,ac((maxyr_stf-5):maxyr_stf)]
```


The stock-recruitment relationship (SRR). A short term forecast does not use an SRR (in the traditional sense). Instead, it generally assumes that recruitment in the future is some mean (e.g. geometric mean) of the historic recruitments. However, we still need to have an SRR that contains this mean value, which is what we mimic for this example. First, we estimate the geometic mean recruitment, that we then add to an SRR object.

```{r}
mean_rec <- exp(mean(log(rec(resul_FLStock))))
ple4_sr <- as.FLSR(resul_FLStock, model="geomean")
params(ple4_sr)['a',] <- mean_rec
params(ple4_sr)
```

The control object. The final thing we need to set up is the control object. This tells the projection what to do, i.e. what level of fishing mortality to use.

```{r}
fbar_SQ <- mean(fbar(resul_FLStock)[,as.character(maxyr_stk)])
ctrl_target <- data.frame(year = 2021:2023, quantity = "f", val = fbar_SQ)
ctrl_f <- fwdControl(ctrl_target)
ctrl_f
```

Running the STF.

```{r}
ple4_sq <- fwd(ple4_stf, ctrl = ctrl_f, sr = ple4_sr)
```

The stock numbers are calculated using the recruitment and future mortality assumptions.

```{r}
stock.n(ple4_sq)[,ac((maxyr_stf-5):maxyr_stf)]
```