In flr/mydas: Tools for conducting Management Strategy Evaluation for data limited stocks
library(knitr)
## Global options
opts_chunk$set(cache =TRUE,
cache.path='cache/mse/',
echo =TRUE,
eval =TRUE,
prompt =FALSE,
comment =NA,
message =FALSE,
warning =FALSE,
tidy =FALSE,
fig.height=6,
fig.width =8,
fig.path ='tex/mse/')
iFig=0
Introduction
mpb is an R package for conducting Management Strategy Evaluation (MSE) and simulating a variety of management procedures (MPs). An MP is the combination of pre-defined data, together with an algorithm to which the data are input to provide a value for a TAC or effort control measure. In this vignette the FLife package is used to condition an Operating Model (OM) using life history parameters and relationships. Both packages are part of FLR (@kell2007flr).
library(plyr)
library(FLCore)
library(FLBRP)
library(FLAssess)
library(FLXSA)
library(ggplotFL)
library(FLasher)
library(FLBRP)
library(FLife)
library(mpb)
theme_set(theme_bw())
#source('~/Desktop/flr/mse/R/msy.R')
FLife
The FLife package is used to create a stock. The first steps are to load the example teleost dataset and select the parameters for albacore.
data(teleost)
teleost
alb=lhPar(teleost[,"Thunnus alalunga"])
alb
The lhPar method is then used to derive the parameters for natural mortality-at=age, based on @gislason2008does, and default parameters and relationships for selection pattern and stock recruitment.
The default parameters can be changed, e.g. by changing a parameter. sl is the standard deviation for the lefthand limb of the double normal selection pattern, here we change it from 2 to 1 to make it steeper.
alb["sel2"]=1
Equilibrium dynamics
The parameters are then used by lhEql to simulate the equilibrium dynamics by combining the spawner/yield per recruit relationships with a stock recruiment relationship.
eql=lhEql(alb)
sel<-function(x)
catch.sel(x)%/%fapex(catch.sel(x))
ggplot(FLQuants(eql,"m","catch.sel"=sel,"mat","catch.wt"))+
geom_line(aes(age,data))+
facet_wrap(~qname,scale="free")+
scale_x_continuous(limits=c(0,20))+
guides(colour=guide_legend(title="Species",title.position="top"))
Figure r iFig=iFig+1; iFig Vectors of m, selection pattern, maturity and weight-at-age.
Estimate equilibrium dynamics and reference points
plot(eql)
Figure r iFig=iFig+1; iFig Expected, equilibrium, dynamics and reference points.
Time series
To go from equilibrium to time series dynamics the FLBRP object created by lhEql can be coerced to an FLStock object.
First change the F time series so that it represents a time series where the stock was origionally lightly exploited, F increased until the stock was overfished and then fishing pressure was reduced to ensure spawning stock biomass was greater than $B_{MSY}$.
fbar(eql)=refpts(eql)["msy","harvest"]*FLQuant(c(rep(.1,19),
seq(.1,2,length.out=40),
seq(2,.7,length.out=11)[-1],
rep(.7,61)))[,1:105]
om=as(eql,"FLStock")
Stochastic dynamics
To simulation random variation in the time series, deviations around the stock recruitment relationship was modelled as a random variable.
nits=100
set.seed(1234)
srDev=FLife:::rlnoise(nits,iter(fbar(eql),1)*0,.3,b=0.0)
plot(srDev)+
geom_point(aes(year,data,col=iter),data=as.data.frame(iter(srDev,c(7,1))))
Figure r iFig=iFig+1; iFig Time series of recruitment deviates
While to generate data for use in the MP, random measurement error was added to the simulated catch per unit effort (CPUE).
set.seed(3321)
uDev =rlnoise(nits,setPlusGroup(stock.n(eql),20)*0,.2,b=0)
These deviates were then used to create a stochastic time series by projecting the dynamics from year 1.
om =propagate(om,nits)
oms=FLStocks("Projection"=fwd(om,fbar=fbar(om)[,-1],
sr=eql,deviances=srDev))
plot(oms[["Projection"]])+
geom_line(aes(year,data,col=iter),
data=as.data.frame(FLQuants(iter(oms[["Projection"]],c(7,1)),"Rec"=rec,"F"=fbar,"SSB"=ssb,"Catch"=catch),drop=TRUE))+
theme(legend.position="none")
Figure r iFig=iFig+1; iFig Stochastic Time series of F, SSB, recruitment and yield
\newpage
Management Procedures
Feedback control
Management of a fish stocks is done using feedback control. The stock is assessed using historical data which is used estimate current stock status and then to project the stock forward under alternative management regulations for a variety of hypotheses and system dynamics. This procedure is then repeated in subsequent year to monitor and adjust the impact of management. MSE does this my simulating a MP. These can either be model based or empirical, i.e. based on a stock assessment or data alone.
In the mpb package there are a variety of MP, e.g. age, biomass and empirical based.
Harvest Control Rule
library(kobe)
hcr= data.frame(stock =c(0.0 ,0.1 , 0.6,2.0),
harvest=c(0.01,0.01, 0.7,0.7))
kobePhase()+
geom_line(aes(stock,harvest),data=hcr,col="orange",size=2)
Figure r iFig=iFig+1; iFig Hockey stick harvest control rule.
Age Based
In this example the MP is based on an Virtual Population Analysis (VPA).
First the control settings are checked by running FLXSA on data simulated by the OM without error and feedback. Ideally there should be no bias in the estimates from the stock assessment
mp=window(setPlusGroup(oms[["Projection"]],20),end=80)
##Assessment
control=FLXSA.control(tol =1e-16, maxit =150,
min.nse=0.3, fse =0.5,
rage =2, qage =10,
shk.n =TRUE, shk.f =TRUE,
shk.yrs=10, shk.ages=10,
window =10, tsrange =10,
tspower=0,
vpa =!TRUE)
idx=FLIndex(index=stock.n(mp)%*%uDev[,dimnames(stock.n(mp))$year])
range(idx)[c("plusgroup","startf","endf")]=c(NA,0.1,.2)
xsa=FLXSA(mp,idx,
control=FLXSA.control(),diag.flag=FALSE)
range(xsa)[c("min","max","plusgroup")]=range(mp)[c("min","max","plusgroup")]
mp=mp+xsa
sr=fmle(as.FLSR(mp,model="bevholt"),control=list(silent=TRUE))
rf=FLBRP(mp,sr)
plot(FLStocks("Stock\nAssessment"=mp,
"Operating\nModel" =window(oms[["Projection"]],end=80)))
Before running the MSE, i.e. using XSA as part of a feedback control procedure, the current reference points need to be estimated.
Then the MSE can be run using the mseXSA function
#source('~/Desktop/flr/FLBRP/R/fwd-setup.R')
oms[["Age"]]=mseXSA(oms[["Projection"]],eql, #OM
mp,control,rf=rf, #MP
srDev=srDev,uDev=uDev, #Random deviates for OM
start=75,end=103,maxF=1.0) #year range
plot(oms[["Age"]])+
geom_line(aes(year,data,col=iter),
data=as.data.frame(FLQuants(iter(oms[["Age"]],c(7,12,19)),"Rec"=rec,"F"=fbar,"SSB"=ssb,"Catch"=catch),drop=TRUE))+
theme(legend.position="none")
Figure r iFig=iFig+1; iFig Time series from the MSE of F, SSB, recruitment and yield
Biomass Based
In mpb there is a biomass dynamic stock assessment, designed to be used as an MP.
First the control object has to be set, i.e. setting best guess, bounds and any priors for parameters.
mp =as(window(oms[["Projection"]],start=20,end=75),"biodyn")
mp@indices=FLQuants("1"=(stock(oms[["Projection"]][,20:74])+
stock(oms[["Projection"]][,21:75]))/2.0)
params( mp)["r"]=.25
mp=fwd( mp,catch=catch(mp))
setParams( mp)=mp@indices[[1]]
setControl(mp)=params(mp)
control( mp)["r",2:4]=c(.05,0.25,1.0)
control( mp)["q1",]=c(-1,.1,1,10)
Then the assessment is run without feedback
mp=fit(mp)
and compared to the OM
setControl(mp)=params(mp)
dat=plot(FLQuants(fit(mp),"Biomass"=stock,"F"=function(x) catch(x)/stock(x)[,dimnames(catch(x))$year]))$data
plot(FLQuants(window(oms[["Projection"]],start=20,end=75),"Biomass"=stock,"F"=function(x) catch(x)/stock(x)))+
geom_line( aes(year,`50%`),data=dat,fill="blue",col="blue")+
geom_ribbon(aes(year,ymax=`75%`,ymin=`25%`),data=dat,alpha=.25,fill="blue",col="blue")
Figure r iFig=iFig+1; iFig Comparision of estimates and simulated time series of harvest rate and stock biomass.
#source('~/Desktop/flr/mpb/R/hcr.R')
setControl(mp)=params(mp)
oms[["Biomass"]]=
mseMPB(window(oms[["Projection"]],start=20,end=103),eql,mp,srDev=srDev,uDev=uDev,start=75,end=103)
plot(window(oms[["Biomass"]],end=100))+
geom_line(aes(year,data,col=iter),
data=as.data.frame(FLQuants(window(iter(oms[["Biomass"]],c(7,12,19)),end=100),"Rec"=rec,"F"=fbar,"SSB"=ssb,"Catch"=catch),drop=TRUE))+
theme(legend.position="none")
Figure r iFig=iFig+1; iFig Time series from the MSE of F, SSB, recruitment and yield
oms[["Biomass2"]]=
mseMPB(window(oms[["Projection"]],start=20,end=103),eql,mp,srDev=srDev,uDev=uDev,ftar=0.5,start=75,end=103)
```{r biodyn-mse-plot-2, echo=FALSE, eval=FALSE
norfolk }
plot(window(oms[["Biomass2"]],end=100))+
geom_line(aes(year,data,col=iter),
data=as.data.frame(FLQuants(window(iter(oms[["Biomass"]],c(7,12,19)),end=100),"Rec"=rec,"F"=fbar,"SSB"=ssb,"Catch"=catch),drop=TRUE))+
theme(legend.position="none")
**Figure `r iFig=iFig+1; iFig`** Time series from the MSE of F, SSB, recruitment and yield
## Empirical
```r
#source('~/Desktop/flr/mydas/R/mse.R')
#source('~/Desktop/flr/mydas/R/sbt-mse.R')
oms[["Emprirical"]]=mydas:::mseEMP(oms[["Projection"]],eql,srDev=srDev,uDev=uDev,start=75,end=103)
plot(window(oms[["Emprirical"]],end=100))+
geom_line(aes(year,data,col=iter),
data=as.data.frame(FLQuants(iter(window(oms[["Emprirical"]],end=100),c(7,1)),"Rec"=rec,"F"=fbar,"SSB"=ssb,"Catch"=catch),drop=TRUE))
Figure r iFig=iFig+1; iFig Time series from the MSE of F, SSB, recruitment and yield
\newpage
Software Versions
r version$version.string- FLCore:
r packageVersion('FLCore')
- FLPKG:
r # packageVersion('FLPKG')
- Compiled:
r date()
- Git Hash:
r system("git log --pretty=format:'%h' -n 1", intern=TRUE)
Author information
Laurence Kell. laurie@seaplusplus.es
Acknowledgements
This vignette and many of the methods documented in it were developed under the MyDas project funded by the Irish exchequer and EMFF 2014-2020. The overall aim of MyDas is to develop and test a range of assessment models and methods to establish Maximum Sustainable Yield (MSY) reference points (or proxy MSY reference points) across the spectrum of data-limited stocks.
References {#References}
flr/mydas documentation built on Jan. 19, 2024, 10:33 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(knitr) ## Global options opts_chunk$set(cache =TRUE, cache.path='cache/mse/', echo =TRUE, eval =TRUE, prompt =FALSE, comment =NA, message =FALSE, warning =FALSE, tidy =FALSE, fig.height=6, fig.width =8, fig.path ='tex/mse/') iFig=0
Introduction
mpb is an R package for conducting Management Strategy Evaluation (MSE) and simulating a variety of management procedures (MPs). An MP is the combination of pre-defined data, together with an algorithm to which the data are input to provide a value for a TAC or effort control measure. In this vignette the FLife package is used to condition an Operating Model (OM) using life history parameters and relationships. Both packages are part of FLR (@kell2007flr).
library(plyr) library(FLCore) library(FLBRP) library(FLAssess) library(FLXSA) library(ggplotFL) library(FLasher) library(FLBRP) library(FLife) library(mpb) theme_set(theme_bw()) #source('~/Desktop/flr/mse/R/msy.R')
FLife
The FLife package is used to create a stock. The first steps are to load the example teleost dataset and select the parameters for albacore.
data(teleost)
teleost
alb=lhPar(teleost[,"Thunnus alalunga"]) alb
The lhPar method is then used to derive the parameters for natural mortality-at=age, based on @gislason2008does, and default parameters and relationships for selection pattern and stock recruitment.
The default parameters can be changed, e.g. by changing a parameter. sl is the standard deviation for the lefthand limb of the double normal selection pattern, here we change it from 2 to 1 to make it steeper.
alb["sel2"]=1
Equilibrium dynamics
The parameters are then used by lhEql to simulate the equilibrium dynamics by combining the spawner/yield per recruit relationships with a stock recruiment relationship.
eql=lhEql(alb)
sel<-function(x) catch.sel(x)%/%fapex(catch.sel(x)) ggplot(FLQuants(eql,"m","catch.sel"=sel,"mat","catch.wt"))+ geom_line(aes(age,data))+ facet_wrap(~qname,scale="free")+ scale_x_continuous(limits=c(0,20))+ guides(colour=guide_legend(title="Species",title.position="top"))
Figure r iFig=iFig+1; iFig Vectors of m, selection pattern, maturity and weight-at-age.
Estimate equilibrium dynamics and reference points
plot(eql)
Figure r iFig=iFig+1; iFig Expected, equilibrium, dynamics and reference points.
Time series
To go from equilibrium to time series dynamics the FLBRP object created by lhEql can be coerced to an FLStock object.
First change the F time series so that it represents a time series where the stock was origionally lightly exploited, F increased until the stock was overfished and then fishing pressure was reduced to ensure spawning stock biomass was greater than $B_{MSY}$.
fbar(eql)=refpts(eql)["msy","harvest"]*FLQuant(c(rep(.1,19), seq(.1,2,length.out=40), seq(2,.7,length.out=11)[-1], rep(.7,61)))[,1:105] om=as(eql,"FLStock")
Stochastic dynamics
To simulation random variation in the time series, deviations around the stock recruitment relationship was modelled as a random variable.
nits=100 set.seed(1234) srDev=FLife:::rlnoise(nits,iter(fbar(eql),1)*0,.3,b=0.0)
plot(srDev)+ geom_point(aes(year,data,col=iter),data=as.data.frame(iter(srDev,c(7,1))))
Figure r iFig=iFig+1; iFig Time series of recruitment deviates
While to generate data for use in the MP, random measurement error was added to the simulated catch per unit effort (CPUE).
set.seed(3321) uDev =rlnoise(nits,setPlusGroup(stock.n(eql),20)*0,.2,b=0)
These deviates were then used to create a stochastic time series by projecting the dynamics from year 1.
om =propagate(om,nits) oms=FLStocks("Projection"=fwd(om,fbar=fbar(om)[,-1], sr=eql,deviances=srDev))
plot(oms[["Projection"]])+ geom_line(aes(year,data,col=iter), data=as.data.frame(FLQuants(iter(oms[["Projection"]],c(7,1)),"Rec"=rec,"F"=fbar,"SSB"=ssb,"Catch"=catch),drop=TRUE))+ theme(legend.position="none")
Figure r iFig=iFig+1; iFig Stochastic Time series of F, SSB, recruitment and yield
\newpage
Management Procedures
Feedback control
Management of a fish stocks is done using feedback control. The stock is assessed using historical data which is used estimate current stock status and then to project the stock forward under alternative management regulations for a variety of hypotheses and system dynamics. This procedure is then repeated in subsequent year to monitor and adjust the impact of management. MSE does this my simulating a MP. These can either be model based or empirical, i.e. based on a stock assessment or data alone.
In the mpb package there are a variety of MP, e.g. age, biomass and empirical based.
Harvest Control Rule
library(kobe) hcr= data.frame(stock =c(0.0 ,0.1 , 0.6,2.0), harvest=c(0.01,0.01, 0.7,0.7)) kobePhase()+ geom_line(aes(stock,harvest),data=hcr,col="orange",size=2)
Figure r iFig=iFig+1; iFig Hockey stick harvest control rule.
Age Based
In this example the MP is based on an Virtual Population Analysis (VPA).
First the control settings are checked by running FLXSA on data simulated by the OM without error and feedback. Ideally there should be no bias in the estimates from the stock assessment
mp=window(setPlusGroup(oms[["Projection"]],20),end=80) ##Assessment control=FLXSA.control(tol =1e-16, maxit =150, min.nse=0.3, fse =0.5, rage =2, qage =10, shk.n =TRUE, shk.f =TRUE, shk.yrs=10, shk.ages=10, window =10, tsrange =10, tspower=0, vpa =!TRUE) idx=FLIndex(index=stock.n(mp)%*%uDev[,dimnames(stock.n(mp))$year]) range(idx)[c("plusgroup","startf","endf")]=c(NA,0.1,.2) xsa=FLXSA(mp,idx, control=FLXSA.control(),diag.flag=FALSE) range(xsa)[c("min","max","plusgroup")]=range(mp)[c("min","max","plusgroup")] mp=mp+xsa sr=fmle(as.FLSR(mp,model="bevholt"),control=list(silent=TRUE)) rf=FLBRP(mp,sr)
plot(FLStocks("Stock\nAssessment"=mp, "Operating\nModel" =window(oms[["Projection"]],end=80)))
Before running the MSE, i.e. using XSA as part of a feedback control procedure, the current reference points need to be estimated.
Then the MSE can be run using the mseXSA function
#source('~/Desktop/flr/FLBRP/R/fwd-setup.R') oms[["Age"]]=mseXSA(oms[["Projection"]],eql, #OM mp,control,rf=rf, #MP srDev=srDev,uDev=uDev, #Random deviates for OM start=75,end=103,maxF=1.0) #year range
plot(oms[["Age"]])+ geom_line(aes(year,data,col=iter), data=as.data.frame(FLQuants(iter(oms[["Age"]],c(7,12,19)),"Rec"=rec,"F"=fbar,"SSB"=ssb,"Catch"=catch),drop=TRUE))+ theme(legend.position="none")
Figure r iFig=iFig+1; iFig Time series from the MSE of F, SSB, recruitment and yield
Biomass Based
In mpb there is a biomass dynamic stock assessment, designed to be used as an MP.
First the control object has to be set, i.e. setting best guess, bounds and any priors for parameters.
mp =as(window(oms[["Projection"]],start=20,end=75),"biodyn") mp@indices=FLQuants("1"=(stock(oms[["Projection"]][,20:74])+ stock(oms[["Projection"]][,21:75]))/2.0) params( mp)["r"]=.25 mp=fwd( mp,catch=catch(mp)) setParams( mp)=mp@indices[[1]] setControl(mp)=params(mp) control( mp)["r",2:4]=c(.05,0.25,1.0) control( mp)["q1",]=c(-1,.1,1,10)
Then the assessment is run without feedback
mp=fit(mp)
and compared to the OM
setControl(mp)=params(mp) dat=plot(FLQuants(fit(mp),"Biomass"=stock,"F"=function(x) catch(x)/stock(x)[,dimnames(catch(x))$year]))$data plot(FLQuants(window(oms[["Projection"]],start=20,end=75),"Biomass"=stock,"F"=function(x) catch(x)/stock(x)))+ geom_line( aes(year,`50%`),data=dat,fill="blue",col="blue")+ geom_ribbon(aes(year,ymax=`75%`,ymin=`25%`),data=dat,alpha=.25,fill="blue",col="blue")
Figure r iFig=iFig+1; iFig Comparision of estimates and simulated time series of harvest rate and stock biomass.
#source('~/Desktop/flr/mpb/R/hcr.R') setControl(mp)=params(mp) oms[["Biomass"]]= mseMPB(window(oms[["Projection"]],start=20,end=103),eql,mp,srDev=srDev,uDev=uDev,start=75,end=103)
plot(window(oms[["Biomass"]],end=100))+ geom_line(aes(year,data,col=iter), data=as.data.frame(FLQuants(window(iter(oms[["Biomass"]],c(7,12,19)),end=100),"Rec"=rec,"F"=fbar,"SSB"=ssb,"Catch"=catch),drop=TRUE))+ theme(legend.position="none")
Figure r iFig=iFig+1; iFig Time series from the MSE of F, SSB, recruitment and yield
oms[["Biomass2"]]= mseMPB(window(oms[["Projection"]],start=20,end=103),eql,mp,srDev=srDev,uDev=uDev,ftar=0.5,start=75,end=103)
```{r biodyn-mse-plot-2, echo=FALSE, eval=FALSE norfolk } plot(window(oms[["Biomass2"]],end=100))+ geom_line(aes(year,data,col=iter), data=as.data.frame(FLQuants(window(iter(oms[["Biomass"]],c(7,12,19)),end=100),"Rec"=rec,"F"=fbar,"SSB"=ssb,"Catch"=catch),drop=TRUE))+ theme(legend.position="none")
**Figure `r iFig=iFig+1; iFig`** Time series from the MSE of F, SSB, recruitment and yield
## Empirical
```r
#source('~/Desktop/flr/mydas/R/mse.R')
#source('~/Desktop/flr/mydas/R/sbt-mse.R')
oms[["Emprirical"]]=mydas:::mseEMP(oms[["Projection"]],eql,srDev=srDev,uDev=uDev,start=75,end=103)
plot(window(oms[["Emprirical"]],end=100))+ geom_line(aes(year,data,col=iter), data=as.data.frame(FLQuants(iter(window(oms[["Emprirical"]],end=100),c(7,1)),"Rec"=rec,"F"=fbar,"SSB"=ssb,"Catch"=catch),drop=TRUE))
Figure r iFig=iFig+1; iFig Time series from the MSE of F, SSB, recruitment and yield
\newpage
Software Versions
r version$version.string- FLCore:
r packageVersion('FLCore') - FLPKG:
r # packageVersion('FLPKG') - Compiled:
r date() - Git Hash:
r system("git log --pretty=format:'%h' -n 1", intern=TRUE)
Author information
Laurence Kell. laurie@seaplusplus.es
Acknowledgements
This vignette and many of the methods documented in it were developed under the MyDas project funded by the Irish exchequer and EMFF 2014-2020. The overall aim of MyDas is to develop and test a range of assessment models and methods to establish Maximum Sustainable Yield (MSY) reference points (or proxy MSY reference points) across the spectrum of data-limited stocks.
References {#References}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.