Management_procedure_biomass_dynamic.Rmd

---
title: "Biomass Dynamic Management Procedures"
author: "Laurence Kell"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
  github_document:
    mathjax: TRUE
  pdf_document:
    fig_width: 6 
    fig_height: 4 
bibliography: bibliography.bib
---


```{r, ini, echo=FALSE, results='hide', message=FALSE, warnings=FALSE, cache=FALSE}
library(knitr)
source("R/ini.R")
options(digits=3)
iFig=0
```

```{r,echo=FALSE}
sink(NULL)
warn=options()$warn
options(warn=-1)
library(ggplotFL)
library(plyr)
library(reshape)
library(mpb)  

theme_set(theme_bw())
options(digits=3)
options(warn=warn)
sink()
```

# Introduction

The `mpb` package implements biomass based methods for stock assessment and simulation testing using Management Strategy Evaluation (MSE).

The main processes influencing the dynamics of exploited populations are gains due to growth and recruitment and losses due to fishing and natural mortality. In a biomass dynamic stock assessment model recruitment, growth and natural mortality are simplified into a single production function ($P$), for example that of [@pella1969generalized].

\begin{figure}
\begin{equation}  B_{t+1}=B_{t}-C_{t}+P_{t}\end{equation}  
\caption{Biomass next year equals the biomass this year less catches and plus production.}
\end{figure}

\begin{figure}
\begin{equation}\frac{r}{p}\cdot~B(1-(\frac{B}{K})^p)\end{equation}  
\caption{Production function.}
\end{figure}

The dynamics are determined by the population growth rate ($r$, in the abscence of density dependence) and the shape of the production function ($p$). if $p=1$ then the maximum sustainable yield (MSY) is found halfway between 0 and virgin biomass ($K$); as p deccreases MSY shifts to the left.
There is seldom suffcient infomation in the catch data to estimate even these few parameters and so additional data are required, e.g. time series of relative abundance from catch per unit effort (CPUE), or surveys.

The provision of fisheries management advice requires the assessment of stock status relative to reference points, the prediction of the response of a stock to management, and checking that predictions are consistent with reality. Biomass dynamic models have been criticised as being too simplistic to capture the actual population dynamics, however, if a simple model can provide robust advice on stock status and the response of a stock to management why use anything more complicated  [@ludwig1985age]?  For example the Pella-Tomlinson model is used by the IWC to set catch limits. Neither the form of the model nor its parameters are meant to provide an accurate representation of the dynamics of the population. Rather, it has been demonstrated by simulation that when a biomass dynamic model is used as part of a management strategy with a harvest control rule (HCR) it allows the robust calculation and setting of catches limits [@butterworth1999experiences]. 

# `biodyn` Class

The main class is `biodyn`, which has methods for importing data, exporting results, fitting models, checking diagnostics, plotting, estimation of uncertainly, projection, simulating HCRs, and for the provision of advice. The robustness of the methods can be simulation tested using MSE. `biodyn` also includes slots for catch, parameters, historical stock status, reference points, diagnostics, and summary statistics (use `??biodyn` for more information)

An object can be created in various way, e.g. using the constructor 
```{r classCreate}
bd=biodyn()
```

or by coercion from another class
```{r classCoerce,eval=FALSE}
data(ple4)
bd=as(ple4,"biodyn")
```

or using an existing text file such as the input file of ASPIC and then coercing the `aspic` object into an object of the `mpb` class
```{r classCoerce2,eval=FALSE,echo=TRUE}
asp=aspic("aspic.inp")
bd =as(asp,"biodyn")
```

Objects for use in simulation can also be created 
```{r classSim,eval=FALSE}
bd=sim()
```

## Plotting

**Plots** are important for examining objects, exploring data, summarising results, checking outputs, and diagnosing problems. 
```{r plot, fig.margin=TRUE, fig.width=4, fig.height=5, fig.cap="Production function with simulated time series"}
bd=window(sim(),end=49)
plot(bd)+
  theme_bw()
```

`mpb` uses `ggplot2` as this allows the basic plots to be modified as required, for example a trajectory can be added to the plot of the production function
```{r plotProduction, fig.margin=TRUE, fig.height=2.5, fig.cap="Simulated CPUE series"}
plotProduction(bd)+
  geom_path( aes(stock,catch),
             model.frame(FLQuants(bd,"stock","catch")))+
  geom_point(aes(stock,catch),
             model.frame(FLQuants(bd,"stock","catch")))+
  theme_bw()+theme(legend.position="none")
```

# Estimation

**Fitting to data** can be done using either maximum likelihood or Monte Carlo Markov Chain (MCMC) simulations. Simulation can help to check robustness by allowing estimated values to be compared with the ones used to generate the data.

```{r fit, echo=TRUE, fig.margin=TRUE, fig.height=4, fig.cap="Simulated stock"}
bd=sim()
```

A CPUE series is needed for fitting and can be simulated using mid year biomass and adding error.

```{r fitU, fig.margin=TRUE, fig.cap="Simulated CPUE series"}
cpue=(stock(bd)[,-dims(bd)$year]+
      stock(bd)[,-1])/2
set.seed(7890)
cpue=rlnorm(1,log(cpue),.2)

ggplot(as.data.frame(cpue))+
  geom_point(aes(year,data))+
  geom_line(aes(year,data),col="salmon",
            data=as.data.frame(stock(bd)))+
  theme_bw()
```

Starting values for the parameters are required. The defaults assume that $r$ is 0.5, the production function is symetric (i.e. p=1) and the $b0$ ratio of the initial biomass to $k$ is 1. MSY should be the same order of magnitude as the catch and so carry capacity ($k$) can be calculated if a guess for $r$ is provided. 

```{r fitGuess2,size="tiny"}
params(bd)["k"]=guessK(params(bd)["r"],mean(catch(bd),na.rm=T),params(bd)["p"])
```

Parameters are also required for catchability ($q$) and the CV for the CPUE indices; if the population parameters are known then the stock can be calculated from the catch and initial values for $q$ and the CV derived.
```{r fitParams,fig.margin=TRUE,fig.width=4,fig.height=6}
setParams(bd)=cpue
params(bd)
```

Before fitting the `control` slot has to be provided with the initial guesses, upper and lower bounds (`min` and `max`), and the `phase` for each parameter.
```{r fitControl,fig.margin=TRUE,fig.width=4,fig.height=6}
setControl(bd)=params(bd)
control(bd)
```
Difficult to estimate parameters may be fixed by setting the `phase` (e.g. for $B_0$ and p) to <0, while parameters can be sequentially estimated by setting phase >0. 


## Maximum Likelihood

**Estimation** can be performed using maximum likelihood
```{r fitrun,fig.margin=TRUE,fig.width=4,fig.height=6}
control(bd)["r",1]=2  
bdHat=fit(bd,cpue)  
```

Since the true parameter values are known the fit can be checked 

```{r fitcheck, fig.margin=TRUE, fig.height=6,fig.cap="A comparison of the true and fitted time series"} 
params(bdHat)  
params(bdHat)/params(bd)

plot(as(list("True"=bd,"Hat"=bdHat),"biodyns"))+
  theme(legend.position="bottom")+
  theme_bw()
```

# Diagnostics

**Goodness of fit** diagnostics are important for replicability, by ensuring that a global solution has actually been found and that assumptions arnt violated, so when the assessment is repeated you get a similar result.

## Residuals
Patterns in residuals from the fits of the CPUE to stock abundance may indicate a violation of models assumptions. Which may result in biased estimates of parameters, reference points and stock trends. While variance estimates obtained from bootstrapping assume that residuals are Independently and Identically Distributed (i.i.d.).

The residuals are in the `diags` slot.
```{r diag,echo=TRUE}
head(bdHat@diags)
```

Checking the distribution of residuals can be done by plotting the obsevered quantiles against the predicted quantiles from the assumed distribution using Q-Q plots. These compare a sample of data (the residuals) on the vertical axis to a statistical population (e.g. from a normal distribution) on the horizontal axis. Any nonlinear patterns may imply that the data are not normally distributed i.e. $X ~ N(0,1)$, for example a systematic departure from a straight line may indicate skewness or over or under dispersion. 

```{r diagQQ, fig.width=4,fig.height=4,fig.margin=TRUE, fig.cap="Quantile-quantile plot to compare residual distribution with the normal distribution."}
rsdl=bdHat@diags
ggplot(rsdl)                                           +
  geom_point( aes(qqx,qqy))                            +
  stat_smooth(aes(qqx,qqHat),method="lm",se=T,fill="blue", alpha=0.1) +
  theme_bw()+theme(legend.position="bottom")               
```

It is assumed that an index is proportional to the stock so when plotting the observed against the fitted values the points should fall around the $y=x$ line, if they do not then the index may not be a good proxy for the stock trend.

```{r diagHat, fig.margin=TRUE, fig.height=4, figwidth=4, fig.cap="Observed CPUE verses fitted, blue line is a linear resgression fitted to points, black the y=x line."}
library(diags)

ggplot(with(rsdl, data.frame(obs=diags:::stdz(obs),hat=diags:::stdz(hat))))   +
    geom_abline(aes(slope=1,intercept=0))                     +
    geom_point( aes(obs,hat))                                 +
    stat_smooth(aes(obs,hat),method="lm", se=F)               +
    theme_bw()+theme(legend.position="bottom")                +
    xlab("Fitted") + ylab("Observed")
```

To look for systematic patterns the residuals can be plotted by year, a lowess smoother helps to identify if the proxy doesnt agree with the estimated stock trend based on the catch
```{r diagYr,fig.height=3, fig.margin=TRUE, fig.cap="Residuals by year, with lowess smoother"}
dat=transform(subset(rsdl,!is.na(residual), 
                     residual=diags::stdz(residual,na.rm=T)))

ggplot(aes(year,residual),data=dat)  +
  geom_hline(aes(yintercept=0))      +
  geom_point()                       +
  stat_smooth(method="loess",se=F)   +
  theme_bw()+theme(legend.position="bottom")                
```

It is also assumed that variance of the index does not vary with the mean, this can be checked by plotting the residuals against the fitted values. 
```{r diagVar,fig.height=3, fig.margin=TRUE, fig.cap="Plot of residuals against fitted value, to check variance relationship."}
ggplot(aes(hat, residual),
       data=subset(rsdl,!is.na(hat) & !is.na(residual)))   +
  geom_hline(aes(yintercept=0))         +
  geom_point()                          +
  stat_smooth(method="loess",se=F)      +
  theme_bw()+theme(legend.position="bottom")               
```

It is assumed that the residuals are not autocorrelated, which can be checked by plotting the residuals against each other with a lag of 1. Significant autocorrelations could be due to an increase in catchability with time, which may result in a more optimistic estimate of current stock status as a decline in the stock may be masked by an increase in catchability.
```{r diagAR, fig.width=4,fig.width=4,fig.margin=TRUE, fig.cap="Plot of autocorrelation, i.e. $residual_{t+1}$ verses $residual_{t}$."}
sum(rsdl$residual^2)

ggplot(rsdl)                                              +
  geom_point( aes(residual,residualLag))                  +
  stat_smooth(aes(residual,residualLag),method="lm",se=F) +
  geom_hline(aes(yintercept=0))     +
  xlab(expression(Residual[t]))     + 
  ylab(expression(Residual[t+1]))   +
  theme_bw()+theme(legend.position="bottom")                 
```
\newpage


## Profiles

Likelihood profiles are useful to check that you are actually at a global solution and not stuck on a small hill with your back to the mountain. They are also useful for evaluating the infomation content of the data and whether different data sets are telling you different things and you need to ask more questions to determine the truth.

The control slot can be used to produce a profile, i.e. fix a parameter or parameters for a range of values and then find the maximum likelihood by estimating the other parameters.


1D
```{r prfl, fig.margin=TRUE, fig.height=3, fig.cap="Likelihood profile for r"}
bdHat=fit(bdHat,cpue)
setControl(bdHat)=params(bdHat)
res=profile(bdHat,which='r',fixed=c('b0','p'),range=seq(0.95,1.03,.002))
ggplot(subset(res,ll<0))+
  geom_line(aes(r,ll))  +
  theme_bw()
```

```{r prflk, fig.margin=TRUE, fig.cap="Likelihood profile for k", eval=FALSE}
control(bdHat)["r","phase"]=1
bdHat=fit(bdHat,cpue)
setControl(bdHat)=params(bdHat)
res=profile(bdHat,which='k',fixed=c('b0','p'),range=seq(0.95,1.03,.002))
ggplot(subset(res,ll<0))+
  geom_line(aes(k,ll))  +
  theme_bw()
```


\newpage
# Uncertainty

A main objective of stock assessment is to estimate uncertainly in stock status. This requires estimates of distributions as well as point estimates. As an example a catch and cpue are simulated and fitted using `biodyn`.
```{r uncertainty, fig.margin=TRUE}
bd   =window(sim(),end=39)
cpue=(stock(bd)[,-dims(bd)$year]+
      stock(bd)[,-1])/2
set.seed(7890)
cpue=rlnorm(1,log(cpue),.2)
bdHat=bd

setParams( bdHat)=cpue
setControl(bdHat)=params(bdHat)
bdHat@control[3:4,"phase"]=-1
bdHat=fit(bdHat,cpue)

sims=as(list("True"=bd,"Best Fit"=bdHat),"biodyns")
```

There are various ways to estimate undercertainty in parameter estimates and quantities derived from them, i.e. use the covariance matrix provided by a maximum likelihood fit, bootstrapping, the jack knife or Bayesian methods such as Monte Carlo Markov Chain,

## Variance/Covariance Matrix

Fitting using maximum likelihood provides the covariance matrix for the parameters. Only the $r$ and $k$ are of interest, as $p$ and $b0$ were fixed and $q$ and $sigma$ are nusiance parameters, i.e. are not of immediate interest but which must be accounted for in the analysis.
```{r uncertaintyCov,fig.height=6, fig.margin=TRUE}
v=vcov(  bdHat)[c("r","k"),c("r","k"),1]
params(bdHat)[c("r","k")]
#refs=mvn(500,p,v)
```

## The Bootstrap 

The Bootstrap can be used to simulate CPUE series replicates and the model refitted.
```{r uncertaintyBoot, fig.height=4, fig.margin=TRUE, fig.cap="Bootstrapped CPUE series"}
set.seed(7890)
cpueBoot      =boot(bdHat)

sims["Bootstrap"]=fit(bdHat,cpueBoot)
```


## Jack knife


The Jack knife is a relatively quick procedure
```{r uncertaintyJackknife, fig.height=4,fig.margin=TRUE, fig.cap="Plot predicted stock trend by index"}
bdJK =fit(bdHat,FLQuant(jackknife(cpue)))

sims["Jack Knife"]=bdJK
```

```{r ,fig.height=4,fig.margin=TRUE, fig.cap=""}
plot(sims)+
  theme_bw()
```

\newpage

## Stock Status

The Precautionary Approach requires stock status to be estimated relative to reference points. The covariance matrix can be used to estimate uncertainty in derived quantities, i.e. those used for management such as $F:F_{MSY}$. 

Marginal densities for stock
```{r refBmsy,fig.margin=TRUE, fig.width=4, fig.height=3,fig.cap="Densities of Stock from different methods for estimating uncertainty.", eval=FALSE}

boot=stock(sims[["Bootstrap"]])[,39]

set.seed(7890) 
jack=randJack(500,stock(sims[[  "Best Fit"]])[,39],
                  stock(sims[["Jack Knife"]])[,39])

bnow=rbind(data.frame(Method="boot",stock=c(boot)),
           data.frame(Method="jack",stock=c(jack)))
  
ggplot(bnow)+ 
  geom_density(aes(x=stock, y=..count..), position = "stack",fill="red")+
  facet_wrap(~Method,scale="free_y",ncol=1)+
  geom_vline(aes(xintercept=c(stock(sims[["Best Fit"]])[,"39"])))+
  theme_bw()
```

Kobe Phase Plot
```{r kobe,eval=FALSE,fig.margin=TRUE,fig.width=4,fig.height=5,fig.caption="Kobe Phase Plots"}
library(kobe)

kb=rbind(data.frame(Method="Boot",kobe(sims[["Bootstrap"]], what="pts")),
         data.frame(Method="Jack",kobe(sims[["Jack Knife"]],what="pts")))
             
ggplot(kb)+ 
  geom_point(aes(stock,harvest))+
  facet_wrap(~Method,scale="free_y",ncol=1)+
  theme_bw()
```

## Projections

Once stock parameters and status has been estimated then projections need to be conducted to inform management.

```{r fdwd, fig.margin=TRUE,fig.width=4, fig.height=6,fig.cap="Projection"}
set.seed(7890)
harvest=rlnorm(100,log(harvest(bdHat))[,-dims(bdHat)$year],.1)
bdHat =fwd(bdHat,harvest=harvest)

plot(bdHat,worm=c(2,8))+    
  theme(legend.position="bottom")+
  theme_bw()
```


## Harvest Control Rules

Use simulated data to run annual, tri-annual, F bound and TAC bounded HCRs

Annual
```{r hcr1}
bd=window(sim(),end=29)
for (i in seq(28,49,1))
  bd=fwd(bd,harvest=hcr(bd,yr=i-1,hyr=i+1:2))

simHCR=as(list("Annual"=bd),"biodyns")
```

Tri-annual
```{r hcr}
bd=window(bd,end=29)
for (i in seq(28,49,3))
  bd=fwd(bd,harvest=hcr(bd,yr=i,hyr=i+1:3))
simHCR["Triennial"]=bd
```

Bound on F
```{r hcrF}
bd=window(bd,end=29)
for (i in seq(28,49,3))
  bd=fwd(bd,harvest=hcr(bd,yr=i,byr=i,hyr=i+1:3,bndF=c(0.9,1.1)))
simHCR["bound F"]=bd
```

Bound on catch
```{r hcrY}
bd=window(bd,end=29)
for (i in seq(28,49,3))
  bd=fwd(bd,catch  =hcr(bd,yr=i,byr=i-1,hyr=i+1:3,tac=TRUE,bndTac=c(0.9,1.1)))

simHCR["bound TAC"]=bd
```

```{r hcrPlot, fig.margin=TRUE,fig.width=4, fig.height=4,fig.cap="Plots of projections"}
plot(simHCR)+
  theme_bw()+
  theme(legend.position="bottom")
```

Process Error and Harvest Control Rule
```{r MC,fig.margin=TRUE,fig.width=6,fig.height=6}
set.seed(7890)
pe=rlnorm(500,FLQuant(0,dimnames=list(year=1:50)),0.5)

bd=window(sim(),end=30)
bd.=bd
bd@stock =propagate(bd@stock, 500)
bd=fwd(bd,harvest=harvest(bd)[,2:30],pe=pe)

for (i in seq(30,48,1))
  bd=fwd(bd,
         catch=hcr(bd,yr=i,hyr=i+1,tac=TRUE,bndTac=c(0.9,1.1)),
         pe   =pe)

plot(bd)+
  theme_bw()
```
\newpage


## Advice

```{r kobe2,fig.margin=TRUE,fig.width=6,fig.height=6}
library(plyr)
library(mpb)  
library(reshape)
library(kobe)

trks=kobe(simHCR[["Annual"]],what="trks")

kobePhase()+
    geom_path( aes(stock,harvest),col="blue",data=subset(trks,pctl=="50%"))
```

```{r kobe3,fig.margin=TRUE,fig.width=6,fig.height=6}
trks=ldply(simHCR,kobe,what="trks")

kobePhase()+
    geom_path( aes(stock,harvest,col=.id),data=subset(trks,pctl=="50%"))
```

## MSE

```{r mse,eval=FALSE}
mseBiodyn
```
\newpage

# References


<!-- ## MCMC -->

<!-- Monte Carlo Markov Chain -->
<!-- ```{r uncertaintyMCMC, eval=FALSE} -->
<!-- sims["MCMC"]=fit(bdHat,cpue,cmdOps=c("-mcmc 1000000, -mcsave 5000")) -->
<!-- ``` -->

<!-- Diagnostics need to be run to make sure that the MCMC has actually estimated a stationary distribution. -->

<!-- ```{r uncertaintyMCMC2, fig.height=4,fig.margin=TRUE,eval=FALSE} -->
<!-- acf(c(params(sims[["MCMC"]])["r"])) -->
<!-- ``` -->

<!-- ```{r ref,eval=FALSE,echo=FALSE} -->
<!-- bdHat@mng -->
<!-- ``` -->

<!-- ```{r ref2,eval=FALSE} -->
<!-- bdHat@mngVcov -->
<!-- ``` -->

<!-- ```{r ref3,echo=FALSE,eval=FALSE,fig.margin=TRUE,fig.width=4,fig.height=6,fig.cap=""} -->
<!-- currentState   =bdHat@mng[c("bbmsy","ffmsy"),"hat",drop=T] -->
<!-- currentStateVar=bdHat@mngVcov[c("bbmsy","ffmsy"), -->
<!--                               c("bbmsy","ffmsy"),drop=T] -->

<!-- refs=mvrnorm(100,currentState,currentStateVar) -->

<!-- ggplot(data=as.data.frame(refs))+ -->
<!--   geom_histogram(aes(x=bbmsy))+ -->
<!--   theme_bw() -->
<!-- ``` -->

<!-- likelihood components -->
<!-- ```{r prflLike, fig.margin=TRUE, fig.height=4, fig.width=4, fig.cap="Likelihood profile by data conmponent, i.e. CPUE series"} -->

<!-- bd=sim() -->

<!-- set.seed(7890) -->
<!-- Us  =FLQuants("Unbiased"     = -->
<!--                 rlnorm(1,log((stock(bd)[,-dims(bd)$year]+ -->
<!--                               stock(bd)[,-1])/2),0.2), -->
<!--               "Increase in q"= -->
<!--                 rlnorm(1,log((stock(bd)[,-dims(bd)$year]+ -->
<!--                               stock(bd)[,-1])/2),0.2)) -->

<!-- setParams( bd)=Us -->
<!-- setControl(bd)=params(bd) -->

<!-- bd@control[3:4,"phase"]=-1 -->
<!-- bd=fit(bd,index=Us) -->
<!-- bd@control[,c("min")]=bd@params*0.1 -->
<!-- bd@control[,c("val")]=bd@params -->
<!-- bd@control[,c("max")]=bd@params*10 -->

<!-- prfl=profile(bd,which='r',index=Us, -->
<!--              range=seq(0.975,1.05,.001)) -->

<!-- ggplot(prfl)+ -->
<!--   geom_path(aes(r,ll,group=index,col=index))+ -->
<!--   facet_wrap(~index,scale="free",ncol=1)          + -->
<!--   theme(legend.position="bottom")+ -->
<!--   theme_bw() -->
<!-- ``` -->

<!-- Profile Slot -->
<!-- ```{r prflADMB,echo=FALSE,eval=FALSE} -->
<!-- bd=fit(bdHat,cpue,cmdOps="-lprof") -->
<!-- prf=subset(bd@profile, param %in% c("bbmsy","ffmsy")) -->
<!-- prf=data.frame(What="Profile",t(daply(prf, .(param), with, sample(value,500,prob=p,replace=T)))) -->
<!-- names(prf)[2:3]=c("Stock","Harvest") -->
<!-- ``` -->