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).
Installation {#Installation}
The simplest way to obtain mpb is to install it from CRAN by using the following command in the R console:
install.packages("mpb", repos = "http://flr-project.org/R")
The repos options can be changed depending on personal preferences and includes options such as choosing the directories in which to install the packages see help(install.packages) for more details.
library(FLife)
library(mydas)
library(FLCore)
library(FLBRP)
library(FLAssess)
library(FLXSA)
library(ggplotFL)
library(FLasher)
library(FLBRP)
library(FLife)
library(mpb)
library(plyr)
theme_set(theme_bw())
load Operating Model conditioned on turbot life history data.
data(om)
and FLBRP with equilibrium dynamics
data(eq)
To simulation random variation in the time series, deviations around the stock recruitment relationship are modelled as a random variable.
Deviates can then used to create a stochastic time series by projecting the dynamics from year 1.
nits=dims(om)$iter
nits=100
set.seed(1234)
srDev=rlnoise(nits,fbar(om)%=%0,.3,b=0.0)
plot(srDev,iter=c(7,2,9))
Figure r iFig=iFig+1; iFig Time series of recruitment deviates
om =propagate(om,nits)
oms=FLStocks("Projection"=fwd(om,fbar=fbar(om)[,-1],deviances=srDev,sr=eq))
plot(oms[["Projection"]],iter=1:3)+
theme(legend.position="none")
Figure r iFig=iFig+1; iFig Stochastic Time series of F, SSB, recruitment and yield
om =propagate(om,nits)
oms=FLStocks("Projection"=fwd(om,fbar=fbar(om)[,-1],deviances=srDev,sr=eq))
plot(oms[["Projection"]],iter=1:3)+
theme(legend.position="none")
Figure r iFig=iFig+1; iFig Stochastic Time series of F, SSB, recruitment and yield
Observation Error Model
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 =rlnorm(nits,setPlusGroup(stock.n(eq),10)*0,.2)
\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
library(FLXSA)
mp=window(setPlusGroup(oms[["Projection"]],10),end=80)
xsaControl=FLXSA.control(tol =1e-09, maxit =150,
min.nse=0.3, fse =1.0,
rage =1, qage =6,
shk.n =TRUE, shk.f =TRUE,
shk.yrs=1, shk.ages=4,
window =10, tsrange =10,
tspower= 0,
vpa =FALSE)
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=xsaControl,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
library(mydas)
#save(oms,eq,mp,xsaControl,rf,srDev,uDev,file="/home/laurence/tmp/oms.RData")
!source('~/Desktop/flr/mydas/R/mseXSA.R')
control=FLXSA.control(tol =1e-16, maxit =150,
min.nse=0.3, fse =0.5,
rage =2, qage =8,
shk.n =TRUE, shk.f =TRUE,
shk.yrs=10, shk.ages=3,
window =10, tsrange =10,
tspower=0,
vpa =!TRUE)
oms=oms[1]
oms["Age"]=mseXSA(oms[["Projection"]],eq,
mp,xsaControl,rf=rf,
sr_deviances=srDev,u_deviances=uDev,
start=75,end=103,interval=1) #maxF=1.0)
#save(mp,oms,file="/home/laurence/tmp/tmp.RData")
plot(oms)+
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.
library(plyr)
library(dplyr)
library(reshape)
library(ggplot2)
library(FLCore)
library(ggplotFL)
library(FLasher)
library(FLBRP)
library(FLife)
library(mpb)
library(mydas)
Then the assessment is run without feedback
and compared to the OM
load("/home/laurence/tmp/oms.RData")
library(FLasher)
library(mpb)
source('~/Desktop/flr/mpb/R/hcr.R')
source('~/Desktop/flr/mpb/R/setMP.R')
mp=setMP(as(window(om,end=75),"biodyn"),
r = 0.25,
k =1000.0,
b0= 0.9,
p = -0.6)
save(om,mp,file="/home/laurence/tmp/t.RData")
plot(mp)
nits=dims(om)$iter
set.seed(1234)
srDev =rlnoise(nits,FLQuant(0,dimnames=dimnames(iter(catch(om),1))),0.3,b=0.0)
uDev =rlnoise(nits,FLQuant(0,dimnames=dimnames(iter(catch(om),1))),0.2,b=0.0)
selDev=rlnoise(nits,FLQuant(0,dimnames=dimnames(iter( m(om),1))),0.2,b=0.0)
eq=FLCore::iter(eq,seq(nits))
oms["Biomass"]=mydas:::mseMPB2(om,eq,mp,start=75,end=100,sr_deviances=srDev,
u_deviances=uDev,sel_deviances=selDev)
plot(window(oms[["Biomass"]],end=100),iter=1:3)+
theme(legend.position="none")
Figure r iFig=iFig+1; iFig Time series from the MSE of F, SSB, recruitment and yield
Empirical
control=FLPar(k1=0.5,k2=0.5,gamma=1)
oms["Emprirical"]=mydas:::mseSBTD(om,eq,control=control,
sr_deviances=srDev,u_deviances=uDev,
start=75,end=100)
plot(window(oms[["Emprirical"]],end=100),iter=1:3)+
theme(legend.position="none")
Figure r iFig=iFig+1; iFig Time series from the MSE of F, SSB, recruitment and yield
res=ldply(oms,omSmry,eq)
\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).
Installation {#Installation}
The simplest way to obtain mpb is to install it from CRAN by using the following command in the R console:
install.packages("mpb", repos = "http://flr-project.org/R")
The repos options can be changed depending on personal preferences and includes options such as choosing the directories in which to install the packages see help(install.packages) for more details.
library(FLife) library(mydas)
library(FLCore) library(FLBRP) library(FLAssess) library(FLXSA) library(ggplotFL) library(FLasher) library(FLBRP) library(FLife) library(mpb) library(plyr) theme_set(theme_bw())
load Operating Model conditioned on turbot life history data.
data(om)
and FLBRP with equilibrium dynamics
data(eq)
To simulation random variation in the time series, deviations around the stock recruitment relationship are modelled as a random variable.
Deviates can then used to create a stochastic time series by projecting the dynamics from year 1.
nits=dims(om)$iter nits=100 set.seed(1234) srDev=rlnoise(nits,fbar(om)%=%0,.3,b=0.0)
plot(srDev,iter=c(7,2,9))
Figure r iFig=iFig+1; iFig Time series of recruitment deviates
om =propagate(om,nits) oms=FLStocks("Projection"=fwd(om,fbar=fbar(om)[,-1],deviances=srDev,sr=eq))
plot(oms[["Projection"]],iter=1:3)+ theme(legend.position="none")
Figure r iFig=iFig+1; iFig Stochastic Time series of F, SSB, recruitment and yield
om =propagate(om,nits) oms=FLStocks("Projection"=fwd(om,fbar=fbar(om)[,-1],deviances=srDev,sr=eq))
plot(oms[["Projection"]],iter=1:3)+ theme(legend.position="none")
Figure r iFig=iFig+1; iFig Stochastic Time series of F, SSB, recruitment and yield
Observation Error Model
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 =rlnorm(nits,setPlusGroup(stock.n(eq),10)*0,.2)
\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
library(FLXSA) mp=window(setPlusGroup(oms[["Projection"]],10),end=80) xsaControl=FLXSA.control(tol =1e-09, maxit =150, min.nse=0.3, fse =1.0, rage =1, qage =6, shk.n =TRUE, shk.f =TRUE, shk.yrs=1, shk.ages=4, window =10, tsrange =10, tspower= 0, vpa =FALSE)
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=xsaControl,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
library(mydas) #save(oms,eq,mp,xsaControl,rf,srDev,uDev,file="/home/laurence/tmp/oms.RData")
!source('~/Desktop/flr/mydas/R/mseXSA.R') control=FLXSA.control(tol =1e-16, maxit =150, min.nse=0.3, fse =0.5, rage =2, qage =8, shk.n =TRUE, shk.f =TRUE, shk.yrs=10, shk.ages=3, window =10, tsrange =10, tspower=0, vpa =!TRUE) oms=oms[1] oms["Age"]=mseXSA(oms[["Projection"]],eq, mp,xsaControl,rf=rf, sr_deviances=srDev,u_deviances=uDev, start=75,end=103,interval=1) #maxF=1.0)
#save(mp,oms,file="/home/laurence/tmp/tmp.RData")
plot(oms)+ 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.
library(plyr) library(dplyr) library(reshape) library(ggplot2) library(FLCore) library(ggplotFL) library(FLasher) library(FLBRP) library(FLife) library(mpb) library(mydas)
Then the assessment is run without feedback
and compared to the OM
load("/home/laurence/tmp/oms.RData") library(FLasher) library(mpb) source('~/Desktop/flr/mpb/R/hcr.R') source('~/Desktop/flr/mpb/R/setMP.R') mp=setMP(as(window(om,end=75),"biodyn"), r = 0.25, k =1000.0, b0= 0.9, p = -0.6) save(om,mp,file="/home/laurence/tmp/t.RData") plot(mp)
nits=dims(om)$iter set.seed(1234) srDev =rlnoise(nits,FLQuant(0,dimnames=dimnames(iter(catch(om),1))),0.3,b=0.0) uDev =rlnoise(nits,FLQuant(0,dimnames=dimnames(iter(catch(om),1))),0.2,b=0.0) selDev=rlnoise(nits,FLQuant(0,dimnames=dimnames(iter( m(om),1))),0.2,b=0.0) eq=FLCore::iter(eq,seq(nits)) oms["Biomass"]=mydas:::mseMPB2(om,eq,mp,start=75,end=100,sr_deviances=srDev, u_deviances=uDev,sel_deviances=selDev)
plot(window(oms[["Biomass"]],end=100),iter=1:3)+ theme(legend.position="none")
Figure r iFig=iFig+1; iFig Time series from the MSE of F, SSB, recruitment and yield
Empirical
control=FLPar(k1=0.5,k2=0.5,gamma=1) oms["Emprirical"]=mydas:::mseSBTD(om,eq,control=control, sr_deviances=srDev,u_deviances=uDev, start=75,end=100)
plot(window(oms[["Emprirical"]],end=100),iter=1:3)+ theme(legend.position="none")
Figure r iFig=iFig+1; iFig Time series from the MSE of F, SSB, recruitment and yield
res=ldply(oms,omSmry,eq)
\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.