In flr/mydas: 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

Installation

Examples

More information

References

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.

Back to Top

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","mydas","FLife"), 
             repos="http://flr-project.org/R")

library(FLCore)
library(FLasher)
library(FLBRP)
library(FLife)

library(plyr)
library(reshape)

library(mydas) 

Back to Top

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)
eq=lhEql(lh)

gTime=c(round(gt(eq)))
fbar(eq)=refpts(eq)["msy","harvest"]%*%FLQuant(c(rep(.1,19),
                                              seq(.1,2,length.out=30)[-30],
                                              seq(2,1.0,length.out=gTime)[-1],
                                              rep(1.0,61)))[,1:105]

om=as(eq,"FLStock")
om=fwd(om,f=fbar(om)[,-1], sr=eq)

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=invAlk(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(len,weight=value),binwidth=1)+
  facet_grid(year~iter,scale="free")+
  xlab("Length (cm)")+ylab("Frequency")+
  coord_cartesian(xlim=c(20,mean(lh["linf"])*1.25))

Figure r iFig=1; iFig Observation error model for turbot length samples.

Back to Top

More information

• 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. [e-mail]

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}

Back to Top

flr/mydas documentation built on Jan. 19, 2024, 10:33 a.m.