An_introduction_to_MSE_using_FLR.Rmd

---
title: "An introduction to MSE with FLR"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
github_document:
  mathjax: TRUE
pdf_document:
  fig_width: 6 
  fig_height: 4 
tags: [FLR MSE]
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")
```
 
Management Strategy Evaluation (MSE) is a framework for evaluating the performance of Harvest Control Rules (HCRs) against prevailing uncertainties (Punt et al. 2016). This tutorial introduces the basic steps for building a single-species MSE: conditioning the operating model, setting up the observation error model, constructing a simple model-based HCR (based on the ICES MSY approach), performing the MSE simulations (including feedback), and producing performance statistics.

## 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/), [FLash](http://www.flr-project.org/FLCore/),[FLXSA](http://www.flr-project.org/FLCore/),[FLBRP](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("FLa4a","FLash","FLXSA","FLBRP","ggplotFL"), repos="http://flr-project.org/R")
```

```{r, pkgs}
# Loads all necessary packages
library(FLa4a)
library(FLash)
library(FLXSA)
library(FLBRP)
library(ggplotFL)
```

# CONDITIONING THE OPERATING MODEL

Conditioning the operating model is a key step in building an MSE analysis, as it allows operating models to be considered "plausible" in the sense that they are consistent with observed data. In this tutorial, the a4a assessment model is fitted to the data within the `ple4` `FLStock` and `ple4.index` `FLIndex` objects to produce `stk`, the operating model `FLStock` object. An mcmc method within the a4a assessment is used to obtain parameter uncertainty, reflected in the `iter` dimension of `stk`.

**Read in stock assessment data**

```{r, dat}
data(ple4)
data(ple4.index)
stk <- ple4
idx <- FLIndices(idx=ple4.index)
```

**Set up the iteration and projection window parameters**  

In this tutorial, 20 iterations are used along with a 12-year projection window. A 3-year period is used for to calculate averages needed for projections (e.g. mean weights, etc.).

```{r, parproj}
it <- 20    # iterations
y0 <- range(stk)["minyear"] # initial data year
dy <- range(stk)["maxyear"] # final data year
iy <- dy+1  # initial year of projection (also intermediate year)
ny <- 12    # number of years to project from intial year
fy <- dy+ny # final year
nsqy <- 3   # number of years to compute status quo metrics
```

**Fit stock assessment model a4a**  

Set up the operating model (including parameter uncertainty) based on fitting the a4a assessment model to data. Reference points are obtained for a "median" `stk` (`stk0`) to mimic best estimates of reference points used in the ICES MSY approach.

```{r, a4a}
# Set up the catchability submodel with a smoothing spline
# (setting up a 'list' allows for more than one index)
qmod <- list(~s(age, k=6))
# Set up the fishing mortality submodel as a tensor spline,
# which allows age and year to interact
fmod <- ~te(replace(age, age>9,9), year, k=c(6,8))
# Set up the MCMC parameters
mcsave <- 100
mcmc <- it * mcsave
# Fit the model
fit <- sca(stk, idx, fmodel = fmod, qmodel = qmod, fit = "MCMC", 
       mcmc = SCAMCMC(mcmc = mcmc, mcsave = mcsave, mcprobe = 0.4))
# Update the FLStock object
stk <- stk + fit
# Reduce to keep one iteration only for reference points
stk0 <- qapply(stk, iterMedians)
```

**Fit stock-recruit model**  

A Beverton-Holt stock-recruit model is fitted for each iteration, with residuals generated for the projection window based on the residuals from the historic period. A stock-recruit model is also fitted to the "median" stk for reference points.

```{r, sr, results="hide"}
# Fit the stock-recruit model
srbh <- fmle(as.FLSR(stk, model="bevholt"), method="L-BFGS-B", 
        lower=c(1e-6, 1e-6), upper=c(max(rec(stk)) * 3, Inf))
srbh0 <- fmle(as.FLSR(stk0, model="bevholt"), method="L-BFGS-B", 
         lower=c(1e-6, 1e-6), upper=c(max(rec(stk)) * 3, Inf))
# Generate stock-recruit residuals for the projection period
srbh.res <- rnorm(it, FLQuant(0, dimnames=list(year=iy:fy)), mean(c(apply(residuals(srbh), 6, sd))))
```
**Calculate reference points and set up the operating model for the projection window**  

Reference points based on the "median" stk, assuming (for illustrative purposes only) that Bpa=0.5Bmsy and Blim=Bpa/1.4. The stf method is applied to the operating model stk object in order to have the necessary data (mean weights, etc.) for the projection window.

```{r, refpts}
# Calculate the reference points
brp <- brp(FLBRP(stk0, srbh0))
Fmsy <- c(refpts(brp)["msy","harvest"])
msy <- c(refpts(brp)["msy","yield"])
Bmsy <- c(refpts(brp)["msy","ssb"])
Bpa <- 0.5*Bmsy
Blim <- Bpa/1.4
# Prepare the FLStock object for projections
stk <- stf(stk, fy-dy, nsqy, nsqy)
```

# SET UP OBSERVATION ERROR MODEL ELEMENTS

**Estimate the index catchabilities from the a4a fit (without simulation)**

Observation error is introduced through the index catchability-at-age  

```{r, idx}
# Set up the FLIndices object and populate it
# (note, FLIndices potentially has more than one index, hence the for loop)
idcs <- FLIndices()
for (i in 1:length(idx)){
#   Set up FLQuants and calculate mean and sd for catchability
    lst <- mcf(list(index(idx[[i]]), stock.n(stk0))) # make FLQuants same dimensions
    idx.lq <- log(lst[[1]]/lst[[2]]) # log catchability of index
    idx.qmu <- idx.qsig <- stock.n(iter(stk,1)) # create quants
    idx.qmu[] <- yearMeans(idx.lq) # allocate same mean-at-age to every year
    idx.qsig[] <- sqrt(yearVars(idx.lq)) # allocate same sd-at-age to every year
#   Build index catchability based on lognormal distribution with mean and sd calculated above
    idx.q <- rlnorm(it, idx.qmu, idx.qsig)
    idx_temp <- idx.q * stock.n(stk)
    idx_temp <- FLIndex(index=idx_temp, index.q=idx.q) # generate initial index
    range(idx_temp)[c("startf", "endf")] <- c(0, 0) # timing of index (as proportion of year)
    idcs[[i]] <- idx_temp
}
names(idcs) <- names(idx)
idx<-idcs[1]
```

# SET UP MSE LOOP

## Needed Functions  

**Observation error model**  
In this tutorial, observation error is applied to the operating model population numbers to obtain an index of abundance. This is implemented through the index catchability-at-age. Observation error is added during each year of the projection window, and is therefore dealt with more easily in a function.
```{r, oem}
o <- function(stk, idx, assessmentYear, dataYears) {
	# dataYears is a position vector, not the years themselves
	stk.tmp <- stk[, dataYears]
	# add small amount to avoid zeros
	catch.n(stk.tmp) <- catch.n(stk.tmp) + 0.1
	# Generate the indices - just data years
	idx.tmp <- lapply(idx, function(x) x[,dataYears])
  # Generate observed index
	for (i in 1:length(idx)) index(idx[[i]])[, assessmentYear] <- 
    stock.n(stk)[, assessmentYear]*index.q(idx[[i]])[, assessmentYear]
  list(stk=stk.tmp, idx=idx.tmp, idx.om=idx)
}
```

**XSA assessment model**  

The assessment model used to parameterise the HCR is XSA, through `FLXSA`. This function sets the control parameters for `FLXSA` and fits the assessment.

```{r, xsa}
xsa <- function(stk, idx){
  # Set up XSA settings
    control  <- FLXSA.control(tol = 1e-09, maxit=99, min.nse=0.3, fse=2.0,
                              rage = -1, qage = range(stk)["max"]-1, shk.n = TRUE, shk.f = TRUE,
                              shk.yrs = 5, shk.ages= 5, window = 100, tsrange = 99, tspower = 0)
  # Fit XSA
  fit <- FLXSA(stk, idx, control)
  # convergence diagnostic (quick and dirty)
  maxit <- c("maxit" = fit@control@maxit)
  # Update stk
  stk   <- transform(stk, harvest = harvest(fit), stock.n = stock.n(fit))
  return(list(stk = stk, converge = maxit))
}
```

**Control object for projections**  

The `fwd` method from `FLash` needs a control object, which is set by this function.

```{r, ctrlproj}
getCtrl <- function(values, quantity, years, it){
	dnms <- list(iter=1:it, year=years, c("min", "val", "max"))
	arr0 <- array(NA, dimnames=dnms, dim=unlist(lapply(dnms, length)))
	arr0[,,"val"] <- unlist(values)
	arr0 <- aperm(arr0, c(2,3,1))
	ctrl <- fwdControl(data.frame(year=years, quantity=quantity, val=NA))
	ctrl@trgtArray <- arr0
	ctrl
}
```

## MSE initialisation

The first year of the projection window is the intermediate year, for which, following the ICES WG timeline, already has a TAC. For this tutorial, the TAC in the final year of data is assumed to be the realised catch in `stk` for the same year, while the TAC in the intermediate year is set equal to the TAC in the final year of data. The `stk` object then needs to be projected through the intermediate year by applying the TAC in the intermediate year with `fwd`.

```{r, mseinit}
vy <- ac(iy:fy)
TAC <- FLQuant(NA, dimnames=list(TAC="all", year=c(dy,vy), iter=1:it))
TAC[,ac(dy)] <- catch(stk)[,ac(dy)]
TAC[,ac(iy)] <- TAC[,ac(dy)] #assume same TAC in the first intermediate year
ctrl <- getCtrl(c(TAC[,ac(iy)]), "catch", iy, it)
# Set up the operating model FLStock object
stk.om <- fwd(stk, control=ctrl, sr=srbh, sr.residuals = exp(srbh.res), sr.residuals.mult = TRUE)
```

## Start the MSE loop  

The MSE loop requires the following:

- observation error model to be applied to generate the index of abundance,
- the stock assessment (XSA) to be applied using this index to generate the management procedure stock object `stk.mp`,
- the resultant SSB estimate to be used in the HCR, along with the reference points obtained earlier, to derive target Fs for the TAC year (the year after the intermediate year), and
- the TAC associated with the target F to be calculated by applying `fwd` to the `stk.mp`.

The final step of the MSE loop is to apply the TAC to the operating model stock object, `stk.om`, by using `fwd`.
  
```{r, mseloop, results="hide"}
set.seed(231) # set seed to ensure comparability between different runs
for(i in vy[-length(vy)]){
  # set up simulations parameters
	ay <- an(i)
	cat(i, " > ")
  flush.console()
	vy0 <- 1:(ay-y0) # data years (positions vector) - one less than current year
	sqy <- (ay-y0-nsqy+1):(ay-y0) # status quo years (positions vector) - one less than current year

  # apply observation error
	oem <- o(stk.om, idx, i, vy0)
	stk.mp <- oem$stk
	idx.mp <- oem$idx
	idx <- oem$idx.om

  # perform assessment
  out.assess <- xsa(stk.mp, idx.mp)
  stk.mp <- out.assess$stk
  
  # apply ICES MSY-like Rule to obtain Ftrgt
  # (note this is not the ICES MSY rule, but is similar)
  flag <- ssb(stk.mp)[,ac(ay-1)]<Bpa
  Ftrgt <- ifelse(flag,ssb(stk.mp)[,ac(ay-1)]*Fmsy/Bpa,Fmsy) 

  # project the perceived stock to get the TAC for ay+1
  fsq.mp <- yearMeans(fbar(stk.mp)[,sqy]) # Use status quo years defined above
  ctrl <- getCtrl(c(fsq.mp, Ftrgt), "f", c(ay, ay+1), it)
  stk.mp <- stf(stk.mp, 2)
  gmean_rec <- c(exp(yearMeans(log(rec(stk.mp)))))
  stk.mp <- fwd(stk.mp, control=ctrl, sr=list(model="mean", params = FLPar(gmean_rec,iter=it)))
  TAC[,ac(ay+1)] <- catch(stk.mp)[,ac(ay+1)]

  # apply the TAC to the operating model stock
  ctrl <- getCtrl(c(TAC[,ac(ay+1)]), "catch", ay+1, it)
  stk.om <- fwd(stk.om, control=ctrl,sr=srbh, sr.residuals = exp(srbh.res), sr.residuals.mult = TRUE)
}
```

# PERFORMANCE STATISTICS

```{r, pstats}
#Some example performance statstics, but first isolate the projection period
stk.tmp<-window(stk.om,start=iy)
# annual probability of being below Blim
(risky<-iterSums(ssb(stk.tmp)/Blim<1)/it)
# mean probabiity of being below Blim in the first half of the projection period
mean(risky[,1:trunc(length(risky)/2)])
# ...and second half
mean(risky[,(trunc(length(risky)/2)+1):length(risky)])
# plot of SSB relative to Bmsy
boxplot(data~year,data=as.data.frame(ssb(stk.tmp)),main="SSB")
abline(h=Bmsy,col="red")
# plot of landings relative to MSY yield
boxplot(data~year,data=as.data.frame(landings(stk.tmp)),main="Landings")
abline(h=msy,col="red")
```

```{r fig1, fig.cap="Results for applying an ICES MSY-like rule, comparing the operating model to the management procedure"}
plot(FLStocks(stk.om=stk.om, stk.mp=stk.mp)) + theme(legend.position="top") + geom_vline(aes(xintercept=2009))
```

# References  
Punt, A.E., Butterworth, D.S., de Moor, C.L., De Oliveira, J.A.A. and M. Haddon (2016). Management Strategy Evaluation: Best Practices. Fish and Fisheries, 17(2): 303-334. DOI: 10.1111/faf.12104.

# 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

**Jose DE OLIVEIRA**. Centre for Environment, Fisheries and Aquaculture Science (CEFAS), Pakefield Road, Lowestoft, Suffolk NR33 0HT, United Kingdom. <http://www.cefas.co.uk/>

**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/>