A_quick_introduction_to_FLR.Rmd

---
title: "A quick introduction to FLR"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
github_document:
  mathjax: TRUE
pdf_document:
  fig_width: 6 
  fig_height: 4 
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)
source("R/ini.R")
```

The Fisheries Library in R (FLR) is a collection of tools for quantitative fisheries science, developed in the R language, that facilitates the construction of bio-economic simulation models of fisheries systems as well as the application of a wide range of quantitative analysis.

FLR builds on the powerful R environment and syntax to create a [domain-specific language](https://en.wikipedia.org/wiki/Domain-specific_language) for the quantitative analysis of the expected risks and impacts of fisheries management decisions. The classes and methods in FLR consider uncertainty an integral part of our knowledge of fisheries systems.

## Required packages

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

- CRAN: latticeExtra, gridExtra, ggplot2
- FLR: [FLCore](http://www.flr-project.org/FLCore/), [FLa4a](http://www.flr-project.org/FLa4a/), [ggplotFL](http://www.flr-project.org/ggplotFL/), [FLBRP](http://www.flr-project.org/FLBRP/), [FLash](http://www.flr-project.org/FLash/)

You can do so as follows:

```{r, eval=FALSE}
install.packages(c("latticeExtra", "gridExtra", "ggplot2", "triangle", "copula", "coda", "mgcv"))
install.packages(c("FLCore", "ggplotFL", "FLa4a", "FLBRP", "FLash"), repos="http://flr-project.org/R")
```

# Classes and methods

The *FLCore* package defines a series of data structures that represent different elements of the fishery system. These structures are defined using R's S4 class system, while functions to operate on those classes are defined as methods. In this way function call (e.g. `plot()`) can be used on different classes and their behaviour is consistent with the class structure and content. For further information on the S4 class system, please look at the [introductory text on the **methods** package](https://stat.ethz.ch/R-manual/R-patched/library/methods/html/Introduction.html).

## FLQuant

The [FLQuant](http://www.flr-project.org/FLCore/reference/FLQuant.html) class is the basic data structure for **FLR**, and it is in essence a 6D array on which data and outputs can be stored, structured along dimensions related to time, space, age or length, etc. It is also able to keep track of the units of measurement employed, as a simple character vector, although various operations on these objects are aware of the units of measurements contained in them.

We will be creating and altering *FLQuant* objects as we go along, but we can look briefly at some basic operations, as the syntax will then be used in other classes in *FLR*.

Let's first load the necessary packages in our session. The **FLCore** package

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

As it is the case for R classes, for example *data.frame*, classes in *FLR* can be constructed using a constructor method with the same name as the class, in this case `FLQuant()`.

```{r, flquant}
FLQuant(1:10)
```

The *FLQuant* constructor can take a variety of inputs, like a vector, a matrix or an array, and also allows you to specify a number of arguments. We will only look now at three of them, `units` for the units of measurement, `quant` for the name of the first dimension, and `dimnames` for the dimension names of the object. For example, to construct an object holding (random) data in tonnes on catches for ages 1 to 4, over a period of 6 years, one could call

```{r, flquantage}
flq <- FLQuant(rlnorm(60), dimnames=list(age=1:4, year=2012:2017), units="t")

flq
```

The object has six dimensions, although some of them (`unit`, `season`, `area` and `iter`) are not used in this case. We can inspect the object in various ways

```{r, flquantinspect}
# A summary of structure and data
summary(flq)

# dimnames
dimnames(flq)

# dims
dim(flq)

# units
units(flq)
```

It can be also subset and modified

```{r, flquantmodify}
# Extract first year
flq[, 1]

# Extract year 2013
flq[, "2013"]

# Set catches on age 1 to zero
flq[1,] <- 0
flq
```

as it can be used in many operations

```{r, flquantarith}
# Product with scalar
flq * 10

# Addition with another FLQuant
flq + (flq * 0.20)

# Sum along years
yearSums(flq)
```

### More information

To learn more about the class, check the [FLQuant help page](http://www.flr-project.org/FLCore/reference/FLQuant.html).

# Loading your data

The first step in an analysis is to get the relevant data into the system. We will need to construct objects of certain **FLR** classes, but first data is loaded into the R session using any of the tools available in the language: `read.csv` for CSV files, `readVPA` and others for fisheries legacy file formats, ...

In this example we will create an object of class *FLStock*, a representation of both inputs and outputs involved in fitting a stock assessment model, from a CSV file containing a table with five columns:

- *slot*, the name of the input quantity, in this case either landings at age, *landings.n*, weight-at-age in the landings, *landings.wt*, and the same quantities for discards, *discards.n* and *discards.wt*.
- *age* and *year*
- *data*, the numeric values themselves.
- *units*, a character string for the units of measurement employed: *kg* and *1000* in this case for kilograms and thousands.

The file is downloaded into a temporary folder, and uncompressed. Simply change the value of `dir` to save the file in another folder.

```{r, getfiles, message=FALSE}
dir <- tempdir()
download.file("http://www.flr-project.org/doc/src/ple4.csv.zip", file.path(dir, "ple4.csv.zip"))
unzip(file.path(dir, "ple4.csv.zip"), exdir=dir)
```

The CSV file can now be loaded as a `data.frame` using `read.csv` and inspected

```{r, loadple4}
dat <- read.csv(file.path(dir, "ple4.csv"))
head(dat)
```

This data.frame contains the time series of landings and discards at age, in thousands, and the corresponding mean weights-at-age, in kg, for [North Sea plaice](http://standardgraphs.ices.dk/ViewCharts.aspx?key=8980) (*Pleuronectes platessa*, ICES ple.27.420).

We can create an object to store the landings-at-age data by subsetting it from the *data.frame*

```{r, subsetlandingsn}
landn <- subset(dat, slot=="landings.n", select=-slot)
```

and then convert into an *FLQuant* using

```{r, convertlandingsn}
landsn <- as.FLQuant(landn)
```

The object can now be inspected and plotted using the FLR-defined methods

```{r, plotlandings.n}
summary(landsn)

plot(landsn)
```

In a similar way, we can now convert the input *data.frame*, containing data for four data elements, as specified in the `slots` column, into an *FLStock* object, `nple4`

```{r, convertple4}
nple4 <- as.FLStock(dat)

summary(nple4)
```

To complete this *FLStock* object we will need to specify the natural mortality, *m*, in this case as a constant value of 0.1 for all ages and years,

```{r, nple4m}
m(nple4) <- 0.1
```

the proportion of natural and fishing mortality that takes place before spawning, assumed to be zero in both cases,

```{r, nple4spwn}
m.spwn(nple4) <- harvest.spwn(nple4) <- 0
```

and the maturity at age, as the proportion mature, as a vector repeated for all years.

```{r, nple4mat}
mat(nple4) <- c(0, 0.5, 0.5, rep(1, 7))
```

We now must compute the overall landings and discards, in biomass, from the age-disaggregated values,

```{r, nple4compute}
landings(nple4) <- computeLandings(nple4)
discards(nple4) <- computeDiscards(nple4)
```

and then the catch slots from both landings and discards

```{r, nple4catch}
catch(nple4) <- computeCatch(nple4, slot="all")
```

The mean weight-at-age in the stock needs to be provided. In this case we will assume it is the same as the one computed from the catch sampling programme

```{r, nple4stockwt}
stock.wt(nple4) <- catch.wt(nple4)
```

We finalize by specifying the fully selected age range, used in the calculation of an overall index of fishing mortality, for example as the mean, using  `fbar()`, or as the maximum value, using `fapex()` across those ages. This information is part of the object's *range*.

```{r, nple4range}
range(nple4, c("minfbar", "maxfbar")) <- c(2, 6)
```

If we now inspect the resulting object we can see that all calculated slots have been assigned the corresponding *units* of measurement, except for those that will hold the estimates coming from a stock assessment model: *stock*, *stock.n* and *stock.wt* for the estimates of abundance, and *harvest* for the estimates of fishing mortality at age.

```{r, nple4}
summary(nple4)

plot(metrics(nple4, Catch=catch, Landings=landings))
data(ple4)
```

### More information

Please see the [Loading your data into FLR](Loading_your_data_into_FLR.html) tutorial for more examples on how to get your data into into **FLR**.

# Visualizing and plotting

As already shown, **FLR** classes always have a `plot` method defined that provides a basic visual exploration of the object contents. But you are not limited to those plots, and other methods are available to allow you to build the plots you need.

The *FLCore* package provides a set of those plot methods based on R's [lattice package](https://cran.r-project.org/web/packages/lattice/index.html), as well as versions of lattice's methods that work directly on **FLR** objects, for example `xyplot`

```{r, xyplot}
xyplot(data~year|age, harvest(ple4), xlab="", ylab="", type="b", cex=0.5, pch=19)
```

But the *ggplotFL* package, which we have already loaded, provides a new set of plot methods based on [ggplot2](http://ggplot2.tidyverse.org/), from which all the plots shown so far in this tutorial have been generated.

One important advantage of *ggplot2* is that plots are returned as objects that can be modified and extended. For example, the basic plot for the abundances at age *FLQuant* using *ggplotFL* would be this one

```{r, plotstockn}
plot(stock.n(ple4))
```

in which each age is shown in a different horizontal panel (*facet*). Now we can modify this plot in various ways by adding extra *ggplot2* commands to it, for example by adding a label to the `y` axis, shading a certain period of interest

```{r, plotstocknplus}
plot(stock.n(ple4)) +
  # Add y label
  ylab("Biomass (t)") +
  # Draw rectangle between years 1990 and 2000
  annotate("rect", xmin = 1990, xmax = 2000, ymin = 0, ymax = Inf,
    # in semi-transparent red
    alpha = .2, fill='red')
```

The basic method in *ggplot2* taking a data input, `ggplot`, has also been defined for some *FLCore* class, which simplifies constructing plots based on objects of those classes from scratch. For example, we can plot our catch-at-age matrix from *ple4* to show the signal of strong cohorts in the data using

```{r, plotcatchn}
ggplot(data=catch.n(ple4), aes(x=year, y=data, group=age)) +
  geom_line(aes(colour=as.factor(age))) +
  ylab("Total catch (t)") + xlab("") + theme(legend.position="none")
```

*ggplotFL* deals with multiple iterations in **FLR** objects by calculating quantiles along the *iter* dimension and plotting them as areas and lines. For example, an object with random lognormal noise around the catch time series would look like this

```{r, plotcatch}
plot(rlnorm(250, log(catch(ple4)), 0.5))
```

where the 50% and 95% intervals are shown as a red region and a dotted line respectively, while the median value is shown as a black line.

Note that plot methods on both *FLCore* and *ggplotFL* work by converting various **FLR** objects to a class that *lattice* and **ggplot2** can understand, *data.frame*s. Coercing a 6D array like an *FLQuant* to a 2D *data.frame* requires generating one column per dimension, named as those, plus an extra column for the values stored in the object. This column is named `data`, thus the use of `data` in various plot formulas above. For classes with multiple slots, like *FLStock*, an extra column is necessary, which will be called `slots`. Please check [this help page](http://www.flr-project.org/FLCore/reference/coerce-methods.html) to know more about the coercion methods for **FLR** objects.

### More information

For further explanations on the use of the two plotting platforms available for **FLR** classes, please refer to the [ggplotFL plotting FLR objects with ggplot2](http://www.flr-project.org/doc/ggplotFL_plotting_FLR_objects_with_ggplot2.html) and [Plotting FLR objects using lattice](http://www.flr-project.org/doc/Plotting_FLR_objects_using_lattice.html) tutorials.

# Running a stock assessment

An important step in providing management advice for a fish stock is the estimation of current status, past history and stock productivity through a stock assessment (SA) model. The *FLStock* object we have created contains some of the inputs required for most SA models: catches, natural mortality, and maturity. What is now needed is some indication of changes in abundance over the period of exploitation, an index of abundance, derived from either research surveys or catch-per-unit-effort of commercial fleets.

## Index of abundance

An index of abundance for the North Sea plaice stock, as an object of class [FLIndex](http://www.flr-project.org/FLCore/reference/FLIndex.html), is available as an example dataset in *FLCore*, and can be loaded using the *data* command

```{r, loadple4index}
data(ple4.index)
```

This time series of abundances at age, covering from `r  dims(ple4.index)$minyear` to `r  dims(ple4.index)$maxyear`, was obtained during the [ICES International Bottom Trawl Survey](http://ocean.ices.dk/Project/IBTS/), as is one of the indices used by the Working Group on the Assessment of Demersal Stocks in the North Sea and Skagerrak (WGNSSK).

The *FLIndex* class is able to contain not only the index of abundance, in this case as an age-structured *FLQuant*, but also other information on variance, associated catches, effort, selectivity and catchability, which might be used by certain methods.

```{r, summaryple4index}
summary(ple4.index)

plot(ple4.index)
```

The *range* slot in this class contains information about the time period over which the survey takes place, as proportions of the year, in this case from the beginning of August to the end of September.

```{r, ple4indexrange}
range(ple4.index)[c("startf", "endf")]
```

## a4a statistical catch-at-age

The [FLa4a](http://flr-project.org/FLa4a) package provides an **FLR** implementation of the *a4a* model, a statistical catch-at-age model of medium complexity that has been designed to make use of the power and flexibility of R's formula-based model notation (as used for example in [lm](https://stat.ethz.ch/R-manual/R-patched/library/stats/html/lm.html)), while giving access to a solid and fast minimization algorithm implemented using [AD Model Builder](http://www.admb-project.org/).

```{r, fla4apkg}
library(FLa4a)
```

A first go at fitting this model to the *ple4* and *ple4.index* datasets requires a single line of code

```{r, sca}
fit <- sca(ple4, FLIndices(BTS=ple4.index))
```

This returns an object of class `r c(class(fit))` that contains the results of the model fit, including the estimated catch-at-age, and the derived abundances and fishing mortality at age

```{r, summarya4afit}
summary(fit)
```

which we can use to update the *FLStock* object

```{r, ple4fit}
stk <- ple4 + fit
plot(stk)
```

Among other methods, *FLa4a* provides a range of tools for inspecting the quality of the model fit. 
<!-- This should be uncommented once a tutorial exists to link to
For example, a set of residuals, for the fits to both the index olf abundance and 


### More information

Please see the [](.html) tutorial for more examples on how to .
 -->

# Exploring the stock-recruitment relationship

Fitting a model for the relationship between stock abundance and expected recruitment is necessary for forecasting the future dynamics of the stock under various levels of fishing mortality, the central element of fisheries management advice. Certain stock assessment models, for example most statistical catch-at-age models, include the stock-recruitment model in the estimated parameters, but for others, like Virtual Population Analysis, this is a setup carried out once the stock and population abundance have been obtained.

Both inputs and outputs for stock-recruitment models are stored in an object of the [FLSR](http://www.flr-project.org/FLCore/reference/FLSR.html) class which we can create directly from the SA results by coercion using

```{r, flsrcreate}
plsr <- as.FLSR(stk)
```

The resulting *FLSR* object contains the series of recruitment `rec`, the abundance of the first age, and the stock reproductive potential, in this case the spawning stock biomass `ssb`. Both series are already displaced by the age of recruitment, 1 in this case.

```{r, flsrsummary}
summary(plsr)
```

We now need to choose a model to fit, and *FLCore* provides a number of them already defined in terms of:

- the model formula, kept in the `model` slot
- a function returning the log-likelihood, to be used for minimization, in the `logl` slot
- a function to obtain initial values for the minimization procedure, in the `initial` slot

See the [SR Models](http://www.flr-project.org/FLCore/reference/FLSR.html) help page for further details of the available models.

We can now assign the output of the chosen model function, in this case `ricker()`, to the existing object using

```{r, flsrmodel}
model(plsr) <- ricker()
```

The model object is ready to fit, through Maximum Likelihood Estimation (MLE), using the `fmle` method defined in *FLCore*. This makes use of R's `optim`, but defaults to the Melder-Mead algorithm.

```{r, flsrfit, results="hide"}
plsr <- fmle(plsr)
```

Another look at the object will show us that a number of slots now contain the results of the model fit, among them:

- `fitted` for the estimated recruitment, as an *FLQuant*
- `residuals` for the residuals in log space, *FLQuant*
- `params` for the estimated parameters, as an object of class *FLPar*
- `vcov` for the variance-covariance matrix, as an *array*

The standard plot for the class shows the model fit and a series of useful diagnostics that will help us evaluating its quality and usefulness for future projections. This can be further explored by likelihood profiling over a range of values of the estimated parameters by simply calling the `profile` method of the object

```{r, flsrprofile}
profile(plsr)
```

Finally, the model can be used to predict future recruitments if provided with a set of new values of SSB, as in

```{r, flsrpredict}
predict(plsr, ssb=FLQuant(rnorm(10, 25e4, sd(ssb(plsr))), dimnames=list(age=1, year=2008:2017)))
```

### More information

Please see the [Modelling stock recruitment with FLSR](Modelling_stock_recruitment_with_FLSR.html) tutorial for more examples and information on the stock-recruitment model fitting capabilities of **FLR**.

# Estimating reference points

To assess current and future status of a stock we need reference points which we hold to represent desirable states for the stock, usually either in terms of stock abundance, like $B_{MSY}$, the biomass expected to produce, on average and under the current exploitation pattern, the Maximum Sustainable Yield ($MSY$), or of fishing mortality, such as the fishing mortality expected to produce $MSY$, $F_{MSY}$.

We can use the [FLBRP package](http://flr-project.org/FLBRP) to calculate a large set of reference points, once we have an estimate of stock status, our *FLStock* object, and of the stock-recruitment relationship, *FLSR*. We can use them both to create the necessary object, of class *FLBRP*, after loading the package

```{r, flrbp, warnings=FALSE}
library(FLBRP)
plrp <- FLBRP(stk, sr=plsr)
summary(plrp)
```

An object of the *FLBRP* class contains a series of times series of observations on dynamics of both stock and fishery, as well as average quantities computed along a specified time period. A series of slots are designed to contain economic variables, such as `price`, variable costs (`vcost`) and fixed costs (`fcost`), which will need to be provided. When those data are available, economic reference points can also be calculated. Finally, a vector of fishing mortality values, `fbar`, has been created, by default to be of length 101 and with values going from $F=0$ to $F=4$, to be used in the calculation procedure.

We can now proceed with the calculation of reference points, through the `brp` method

```{r, brp}
plrp <- brp(plrp)
```

and the extraction of the values obtained, stored in the `refpts` slot of the object

```{r, refpts}
refpts(plrp)
```

This slot contains the estimates for a series of reference points and for a number of quantities. For example, the value stored at the intersection of `msy` and `harvest` corresponds to $F_{MSY}$, while that of `virgin` and `ssb` is the estimate of virgin stock biomass. We can, for example, extract the two standard MSY-related reference points as a named vector using

```{r, refptsextract}
pmsy <- refpts(plrp)["msy", c("harvest", "ssb"), drop=TRUE]
```

to then use them to visualize the status of the stock in relation to these points, for example

```{r, refptplot}
plot(ssb(stk) / pmsy["ssb"]) + geom_hline(aes(yintercept=1), linetype=2) +
  ylab(expression(SSB / SSB[MSY]))
```

### More information

For more information on *FLBRP*, please see the [Reference points for fisheries management with FLBRP](Reference_points_for_fisheries_management_with_FLBRP.html) tutorial.

# Forecasting under different scenarios

Our final step in this quick tour takes us to the tools in *FLR* for forecasting the future responses of the stock to alternative management decisions. This example concerns a Short Term Forecast, in which the effect of the current management, for example a maximum catch limit or TAC, on the present year is forecast, and then a series of alternative increases or decreases in catch are applied to the object representing our current understanding of the stock status and dynamics.

The [FLash package](http://flr-project.org/FLash) contains the current implementation of the forecasting mechanism in *FLR*. Please note that FLash is currently only working on 32 bit R if using Windows. A new forecasting package, called *[FLasher](http://flr-project.org/FLasher)* has been released and will only work on 64 bit R.

```{r, flash}
library(FLash)
```

The key element in the *FLash* interface is the *fwdControl* class, which allows us to specify the targets and limits that we want the stock to achieve or remain within. But first we need to make some assumptions about the future values of certain quantities, both on the biology (natural mortality, maturity) and the fishery (selectivity, proportion of discards). A common approach for this kind of forecast is to use averages over the last few years (typically three) for certain quantities. This can be done in *FLR* using the `stf` method in the *FLash* package.

```{r, stf}
library(FLash)
proj <- stf(stk)
```

Our new object extends up to 2020, and the values for certain stocks have been extended using the mean of the last three years, for example for the weight-at-age in the stock

```{r, stfstockwt}
stock.wt(proj)[, ac(2011:2020)]
```

We can now construct a `fwdControl` object that projects the stock as if catches in the current year, 2018 in this case, were those of the agreed TAC, 85,000 t, and F was then kept in 2019 and 2020 at the same levels as in 2017. The main component of a *fwdControl* object is a *data.frame* that must at least contains columns for

- *year*
- *quantity*, the quantity to be projected, for example `ssb`, `f` or `catch`
- *val*, the actual value to be obtained

```{r, fwdcontrol}
TAC <- 85000
Flevel <- fbar(stk)[,"2017"]
ctrl <- fwdControl(data.frame(year=2018:2020, quantity=c("catch", "f", "f"), val=c(TAC, Flevel, Flevel)))
```

Now the `fwd` method will carry out the projection using those three object: stock, stock-recruitment and control

```{r, fwd}
proj <- fwd(proj, control=ctrl, sr=plsr) 
```

and the results can be visualized

```{r, fwdplot}
plot(proj) + geom_vline(aes(xintercept=2017.5), linetype=2)
```

### More information

Please see the [Short Term Forecasting for advice using FLash](Short_Term_Forecasting_for_advice_using_FLash.html)  and [Short Term Forecasting for advice using FLash](Short_Term_Forecasting_for_advice_using_FLash.html) tutorials for more examples on how to use the capabilities of *FLash* to forecast population dynamics under various scenarios.

# More on FLR

Please visit [our website](http://flr-project.org) for more information on the FLR packages, tutorials and examples of use, and instructions for getting in touch if you want to report a bug, ask a question or get involved in the development of the **FLR** toolset. The project is run by a team of fisheries scientists that decided to share the R code they were writing and realised that using a common platform for Quantitative Fisheries Science in R could break the Babel tower spell that J. Schnute first identified

"The cosmic plan for confounding software languages seems to be working remarkably
well among the community of quantitative fishery scientists!"

Schnute et al. 2007


# References

L. T. Kell, I. Mosqueira, P. Grosjean, J-M. Fromentin, D. Garcia, R. Hillary, E. Jardim, S. Mardle, M. A. Pastoors, J. J. Poos, F. Scott, R. D. Scott; FLR: an open-source framework for the evaluation and development of management strategies. *ICES J Mar Sci* 2007; 64 (4): 640-646. doi: [10.1093/icesjms/fsm012](https://doi.org/10.1093/icesjms/fsm012).

# More information

* You can submit bug reports, questions or suggestions on this tutorial at <https://github.com/flr/doc/issues>.
* Alternatively, 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')`
* FLa4a: `r packageVersion('FLa4a')`
* FLBRP: `r packageVersion('FLBRP')`
* FLash: `r packageVersion('FLash')`
* **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

**Iago MOSQUEIRA**. European Commission, DG Joint Research Centre, Directorate D - Sustainable Resources, Unit D.02 Water and Marine Resources, Via E. Fermi 2749, 21027 Ispra VA, Italy. <https://ec.europa.eu/jrc/>.