In laurieKell/mydas-pkg: Tools for conducting Management Strategy Evaluation for data limited stocks
library(knitr)
#source("R/ini.R")
library(knitr)
## Global options
opts_chunk$set(echo =!TRUE,
eval =TRUE,
cache =TRUE,
cache.path="cache/oem/",
prompt =FALSE,
comment =NA,
message =FALSE,
tidy =FALSE,
warnings=FALSE,
fig.height=4.5,
fig.width =6,
fig.path ="tex/oem-")
iFig=0
Introduction {#Introduction}
In Management Strategy Evaluation (MSE) an Operating Model (OM) is used to simulate resource dynamics to evaluate the performance of a Management Procedure (MP). Where the MP is the combination of pre-defined data, together with an algorithm to which such data are input to provide a value for a management control measure.
The link between the OM and the MP is the Observation Error Model (OEM), which generates fishery-dependent or independent resource monitoring data. The OEM models the uncertainties due to sampling and limited data.
Installation
To run the code in this vignette a number of packages need to be installed, both from CRAN and
FLR, where tutorials are available.
Load Libraries
CRAN
All the examples make extensive use of the packages of Hadley Wickham
library(plyr)
library(reshape)
library(ggplot2)
FLR
The FLR packages can be installed from www.flr-project.org
install.packages(c("FLCore","FLFishery","FLasher","FLBRP","mpb","FLife"),
repos="http://flr-project.org/R")
library(FLCore)
library(FLasher)
library(FLBRP)
library(FLife)
Devtools
The MyDas package has various methods for conditioning OMs, generating pseudo data and simulation testing data-limited stock assessment methods. As mydas is still being developed, it is currently best installed from its GitHub repository. It can be installed using the devtools package, (RTools)[https://cran.r-project.org/bin/windows/Rtools/] also needs to be installed, see (guide)[https://www.coursera.org/lecture/data-scientists-tools/installing-rtools-ecVeq].
install.packages("devtools",dependencies=TRUE)
The mydas pakage can then be installed.
library(devtools)
devtools::install_github("lauriekell/mydas-pkg")
library(mydas)
Example
Operating Model
Create an Operating Model by specifying the life histoy parameters and filling in any missing values
lh=FLPar(c(linf= 59.1, k=0.28, t0=-0.4, a=0.01111,b=3.15,a50=4.0,l50=43.25),units="NA")
lh=lhPar(lh)
Then create an FLBRP object that represents the equilibrium, i.e. expected dynamics and reference points.
eq=lhEql(lh)
Simulate a exploitation history and current depletion level.
The stock was originally underexploited, then effort increased until the stock was over fished at 2$F_{MSY}$. At which point a recovery plan was implemented, where F was reduced to $F_{MSY}$ in 1 generation time, after which the stock was exploited at $F_{MSY}$
gTime=c(round(gt(eq)))
fbar(eq)=refpts(eq)["msy","harvest"]%*%FLQuant(c(rep(.1,19),
seq(.1,2,length.out=30),
seq(2,1.0,length.out=gTime)[-1],
rep(1.0,61)))[,1:105]
plot(fbar(eq))
Figure r iFig=iFig+1; iFig Exploitation History.
om=as(eq,"FLStock")
om=fwd(om,f=fbar(om)[,-1], sr=eq)
plot(om)
Catch per unit effort
An unbiased index of total biomass
plot(stock(om))
Figure r iFig=iFig+1; iFig Stock biomass
plot(FLQuants(om,Survey=stock,
CPUE=function(x) catch(x)/fbar(x)))
Figure r iFig=iFig+1; iFig Comparison between unbiased survey and catch per unit effort
library(plyr)
dat=as.data.frame(FLQuants(om,Survey=stock,
CPUE=function(x) catch(x)/fbar(x)))
dat=ddply(dat,.(qname), transform, index=data/mean(data))
ggplot(dat)+geom_line(aes(year,index,col=qname))+
theme(legend.position="bottom")
Figure r iFig=iFig+1; iFig Comparison between unbiased survey and catch per unit effort
Effort creep
Model an increase in catchabilty, e.g. due to an increase in catchability.
trend<-function(object,bias=0.02)
object%*%FLQuant(cumprod(1+rep(bias,dim(object)[2])),dimnames=dimnames(object))
plot(FLQuants(stock(om)%=%1,trend(stock(om)%=%1)))
Figure r iFig=iFig+1; iFig 2% increase per year in catchability.
plot(FLQuants(catch(om)%/%fbar(om),trend(catch(om)%/%fbar(om))))
Figure r iFig=iFig+1; iFig Effect of a 2% increase per year in catchability on CPUE.
Hyperstabilty
Hyperstability is when catch rate stays stable while the actual fish population declines.
hyperstability<-function(object,omega=1,ref=apply(object,c(1,3:6),mean))
ref%*%((object%/%ref)^omega)
om2=window(om,start=25,end=50)
plot(catch(om2)%/%fbar(om2))+
geom_line(aes(year,data),
data=as.data.frame(hyperstability(catch(om2)%/%fbar(om2),omega=0.5)),col="red")
Figure r iFig=iFig+1; iFig Effect of hyperstabilty (red), the stock appears to decline less.
It is possible to develop indices tailored to different case studies
plot(FLQuants(om,"Stock" =stock,
"Adults" =ssb,
"Recruits" =rec,
"Juveniles"=function(x) stock(x)-ssb(x),
"CPUE" =function(x) catch(x)%/%fbar(x)))
Figure r iFig=iFig+1; iFig Different indices of abundance.
Uncertainty
cv=rlnorm(100,log(stock(om)),0.3)
plot(cv, iter=1)
Figure r iFig=iFig+1; iFig Measurement error, the red line is a single realisation.
cv=rlnoise(100,log(stock(om)),0.3,0.8)
plot(cv, iter=1)
Figure r iFig=iFig+1; iFig Measurement error with auto-correlation, the red line is a single realisation.
set.seed(1234)
u =FLQuants("Unbiased" =rlnorm(100,log(stock(om)),.3),
"AR" =rlnoise(100,log(stock(om)),.3,b=.7),
"Hyperstability" =rlnorm(100,log(hyperstability(stock(om),0.5)),.3),
"Trend" =rlnorm(100,log(trend(stock(om),0.01)),.3),
"Hetroscedascity"=rlnorm(100,log(stock(om)),.3)*trend(stock(om)%=%1,0.01),
"Juvenile" =rlnorm(100,log(stock(om)-ssb(om)),.3),
"Mature" =rlnorm(100,log(ssb(om)),.3),
"Numbers" =rlnorm(100,log(apply(stock.n(om),2:6,sum)),.3))
u=FLQuants(llply(u,function(x) x/mean(x)))
u=ldply(u,as.data.frame)
u.=ddply(u,.(year,.id), with, quantile(data,na.rm=T))
ggplot()+
geom_line(aes(year,data,col=factor(iter)),
data=subset(u,iter%in%c(2,11)))+
geom_ribbon(aes(year,ymin=`25%`,ymax=`75%`),data=u.,col="grey",alpha=.5)+
facet_wrap(~.id,ncol=2)+
theme_bw()+theme(legend.position="none")
Figure r iFig=iFig+1; iFig Example of OEM for indices of abundance.
Mean length
Create an index of mean length
mnLn=apply(wt2len(stock.wt(om),lh)%*%stock.n(om),c(2,6),sum)%/%apply(stock.n(om),c(2.,6),sum)
and compare with F
plot(FLQuants(mnLn,fbar(om)))
Figure r iFig=iFig+1; iFig Comparison of mean length and fishing mortality.
First create an Age Length Key that will be used to slice up numbers-at-age
ak=alk(lh)
Then use the ALK to sample lengths
lfd=lenSample(catch.n(om)[,20:65],ak,nsample=500)
ggplot(melt(lfd[,seq(1,45,10)]))+
geom_histogram(aes(length,weight=value),binwidth=1)+
facet_grid(year~iter,scale="free")+
xlab("Length (cm)")+ylab("Frequency")+
coord_cartesian(xlim=c(0,mean(lh["linf"])))
Figure r iFig=1; iFig Observation error model for turbot length samples.
- You can submit bug reports, questions or suggestions on
FLPKG at the FLPKG issue page ^[https://github.com/flr/FLPKG/issues], or on the FLR mailing list.
- Or send a pull request to https://github.com/flr/FLPKG/
- For more information on the FLR Project for Quantitative Fisheries Science in R, visit the FLR webpage ^[http://flr-project.org].
- The latest version of
FLPKG can always be installed using the devtools package, by calling
library(devtools)
install_github('flr/FLPKG')
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}
laurieKell/mydas-pkg documentation built on Nov. 8, 2019, 10:58 p.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) #source("R/ini.R")
library(knitr) ## Global options opts_chunk$set(echo =!TRUE, eval =TRUE, cache =TRUE, cache.path="cache/oem/", prompt =FALSE, comment =NA, message =FALSE, tidy =FALSE, warnings=FALSE, fig.height=4.5, fig.width =6, fig.path ="tex/oem-") iFig=0
Introduction {#Introduction}
In Management Strategy Evaluation (MSE) an Operating Model (OM) is used to simulate resource dynamics to evaluate the performance of a Management Procedure (MP). Where the MP is the combination of pre-defined data, together with an algorithm to which such data are input to provide a value for a management control measure.
The link between the OM and the MP is the Observation Error Model (OEM), which generates fishery-dependent or independent resource monitoring data. The OEM models the uncertainties due to sampling and limited data.
Installation
To run the code in this vignette a number of packages need to be installed, both from CRAN and FLR, where tutorials are available.
Load Libraries
CRAN
All the examples make extensive use of the packages of Hadley Wickham
library(plyr) library(reshape) library(ggplot2)
FLR
The FLR packages can be installed from www.flr-project.org
install.packages(c("FLCore","FLFishery","FLasher","FLBRP","mpb","FLife"), repos="http://flr-project.org/R")
library(FLCore) library(FLasher) library(FLBRP) library(FLife)
Devtools
The MyDas package has various methods for conditioning OMs, generating pseudo data and simulation testing data-limited stock assessment methods. As mydas is still being developed, it is currently best installed from its GitHub repository. It can be installed using the devtools package, (RTools)[https://cran.r-project.org/bin/windows/Rtools/] also needs to be installed, see (guide)[https://www.coursera.org/lecture/data-scientists-tools/installing-rtools-ecVeq].
install.packages("devtools",dependencies=TRUE)
The mydas pakage can then be installed.
library(devtools) devtools::install_github("lauriekell/mydas-pkg")
library(mydas)
Example
Operating Model
Create an Operating Model by specifying the life histoy parameters and filling in any missing values
lh=FLPar(c(linf= 59.1, k=0.28, t0=-0.4, a=0.01111,b=3.15,a50=4.0,l50=43.25),units="NA") lh=lhPar(lh)
Then create an FLBRP object that represents the equilibrium, i.e. expected dynamics and reference points.
eq=lhEql(lh)
Simulate a exploitation history and current depletion level.
The stock was originally underexploited, then effort increased until the stock was over fished at 2$F_{MSY}$. At which point a recovery plan was implemented, where F was reduced to $F_{MSY}$ in 1 generation time, after which the stock was exploited at $F_{MSY}$
gTime=c(round(gt(eq))) fbar(eq)=refpts(eq)["msy","harvest"]%*%FLQuant(c(rep(.1,19), seq(.1,2,length.out=30), seq(2,1.0,length.out=gTime)[-1], rep(1.0,61)))[,1:105] plot(fbar(eq))
Figure r iFig=iFig+1; iFig Exploitation History.
om=as(eq,"FLStock") om=fwd(om,f=fbar(om)[,-1], sr=eq) plot(om)
Catch per unit effort
An unbiased index of total biomass
plot(stock(om))
Figure r iFig=iFig+1; iFig Stock biomass
plot(FLQuants(om,Survey=stock, CPUE=function(x) catch(x)/fbar(x)))
Figure r iFig=iFig+1; iFig Comparison between unbiased survey and catch per unit effort
library(plyr) dat=as.data.frame(FLQuants(om,Survey=stock, CPUE=function(x) catch(x)/fbar(x))) dat=ddply(dat,.(qname), transform, index=data/mean(data)) ggplot(dat)+geom_line(aes(year,index,col=qname))+ theme(legend.position="bottom")
Figure r iFig=iFig+1; iFig Comparison between unbiased survey and catch per unit effort
Effort creep
Model an increase in catchabilty, e.g. due to an increase in catchability.
trend<-function(object,bias=0.02) object%*%FLQuant(cumprod(1+rep(bias,dim(object)[2])),dimnames=dimnames(object))
plot(FLQuants(stock(om)%=%1,trend(stock(om)%=%1)))
Figure r iFig=iFig+1; iFig 2% increase per year in catchability.
plot(FLQuants(catch(om)%/%fbar(om),trend(catch(om)%/%fbar(om))))
Figure r iFig=iFig+1; iFig Effect of a 2% increase per year in catchability on CPUE.
Hyperstabilty
Hyperstability is when catch rate stays stable while the actual fish population declines.
hyperstability<-function(object,omega=1,ref=apply(object,c(1,3:6),mean)) ref%*%((object%/%ref)^omega)
om2=window(om,start=25,end=50) plot(catch(om2)%/%fbar(om2))+ geom_line(aes(year,data), data=as.data.frame(hyperstability(catch(om2)%/%fbar(om2),omega=0.5)),col="red")
Figure r iFig=iFig+1; iFig Effect of hyperstabilty (red), the stock appears to decline less.
It is possible to develop indices tailored to different case studies
plot(FLQuants(om,"Stock" =stock, "Adults" =ssb, "Recruits" =rec, "Juveniles"=function(x) stock(x)-ssb(x), "CPUE" =function(x) catch(x)%/%fbar(x)))
Figure r iFig=iFig+1; iFig Different indices of abundance.
Uncertainty
cv=rlnorm(100,log(stock(om)),0.3) plot(cv, iter=1)
Figure r iFig=iFig+1; iFig Measurement error, the red line is a single realisation.
cv=rlnoise(100,log(stock(om)),0.3,0.8) plot(cv, iter=1)
Figure r iFig=iFig+1; iFig Measurement error with auto-correlation, the red line is a single realisation.
set.seed(1234) u =FLQuants("Unbiased" =rlnorm(100,log(stock(om)),.3), "AR" =rlnoise(100,log(stock(om)),.3,b=.7), "Hyperstability" =rlnorm(100,log(hyperstability(stock(om),0.5)),.3), "Trend" =rlnorm(100,log(trend(stock(om),0.01)),.3), "Hetroscedascity"=rlnorm(100,log(stock(om)),.3)*trend(stock(om)%=%1,0.01), "Juvenile" =rlnorm(100,log(stock(om)-ssb(om)),.3), "Mature" =rlnorm(100,log(ssb(om)),.3), "Numbers" =rlnorm(100,log(apply(stock.n(om),2:6,sum)),.3)) u=FLQuants(llply(u,function(x) x/mean(x))) u=ldply(u,as.data.frame) u.=ddply(u,.(year,.id), with, quantile(data,na.rm=T)) ggplot()+ geom_line(aes(year,data,col=factor(iter)), data=subset(u,iter%in%c(2,11)))+ geom_ribbon(aes(year,ymin=`25%`,ymax=`75%`),data=u.,col="grey",alpha=.5)+ facet_wrap(~.id,ncol=2)+ theme_bw()+theme(legend.position="none")
Figure r iFig=iFig+1; iFig Example of OEM for indices of abundance.
Mean length
Create an index of mean length
mnLn=apply(wt2len(stock.wt(om),lh)%*%stock.n(om),c(2,6),sum)%/%apply(stock.n(om),c(2.,6),sum)
and compare with F
plot(FLQuants(mnLn,fbar(om)))
Figure r iFig=iFig+1; iFig Comparison of mean length and fishing mortality.
First create an Age Length Key that will be used to slice up numbers-at-age
ak=alk(lh)
Then use the ALK to sample lengths
lfd=lenSample(catch.n(om)[,20:65],ak,nsample=500)
ggplot(melt(lfd[,seq(1,45,10)]))+ geom_histogram(aes(length,weight=value),binwidth=1)+ facet_grid(year~iter,scale="free")+ xlab("Length (cm)")+ylab("Frequency")+ coord_cartesian(xlim=c(0,mean(lh["linf"])))
Figure r iFig=1; iFig Observation error model for turbot length samples.
- You can submit bug reports, questions or suggestions on
FLPKGat theFLPKGissue page ^[https://github.com/flr/FLPKG/issues], or on the FLR mailing list. - Or send a pull request to https://github.com/flr/FLPKG/
- For more information on the FLR Project for Quantitative Fisheries Science in R, visit the FLR webpage ^[http://flr-project.org].
- The latest version of
FLPKGcan always be installed using thedevtoolspackage, by calling
library(devtools) install_github('flr/FLPKG')
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.