library(knitr)
library(knitr) ## Global options opts_chunk$set(cache =TRUE, echo =FALSE, eval =TRUE, prompt =FALSE, comment =NA, message =FALSE, warning =FALSE, tidy =TRUE, fig.height=6, fig.width =8, fig.path ="tex/simtest/srar-", cache.path="cache/simtest/sra/")
options(digits=3) iFig=0
Simuation testing of catch based methods.
ICES defines category 1 stocks as those with quantitative assessments, i.e.where full analytical assessments and forecasts can be conducted. This includes methods based on production models. While category 4 stocks are those for which only reliable catch data are available.
An Operating Model is used to simulate data to first test a biomass dynamic stock assessment model using catch and catch per unit effort data, and then catch only are used with a Stock Reduction Analysis (SRA).
To follow this tutorial you should have installed the following packages:
# Load packages library(ggplot2) library(plyr) library(reshape) library(FLCore) library(ggplotFL) library(FLBRP) library(FLasher) library(FLife) library(mydas)
Turbot
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), 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)
plot(FLQuants(om, "f" = function(x) fbar(x)%/%refpts(eq)["msy","harvest"], "ssb" = function(x) ssb(x)%/%refpts( eq)["msy","ssb"], "catch"=function(x) landings(x)%/%refpts(eq)["msy","yield"], "rec" = function(x) rec(x)%/%refpts( eq)["msy","rec"])) + geom_hline(aes(yintercept=1),col="red",linetype=2)+ theme_bw()
Figure r iFig=iFig+1; iFig
Time series relative to MSY benchmarks.
Set up biomass dynamic model, by setting starting parameters etc.
Get priors
library(popbio) prior=popdyn(lh) prior
library(mpb) ## Bug, need to add to NAMESPACE setMP=mpb:::setMP p=mpb:::p om=window(om,end=55) mp=setMP(as(om,"biodyn"), r =median(prior["r"],na.rm=T), k =median(prior["v"],na.rm=T), b0=0.8, p =median(p(prior["bmsy"]/prior["v"]),na.rm=TRUE))
mp=fit(mp,stock(mp)[,20:54])
plot(mp)
Figure r iFig=iFig+1; iFig
Biomass dynamic model fit to the Operating Model
plot(as(list("MP"=mp,"OM"=as(om,"biodyn")),"biodyns"))
Figure r iFig=iFig+1; iFig
Biomass dynamic model fit to the Operating Model
sra=mp dplIndex =window(stock(mp,0.5)%/%params(sra)["k"],start=20,end=54) dplIndex[,ac(c(22:52))]=NA ## change starting values and bounds for catchability control(sra)["q1",2:4]=c(100,1000,10000) sra=fit(sra,dplIndex)
plot(as(list("SRA"=sra,"MP"=mp,"OM"=as(om,"biodyn")),"biodyns"))
Figure r iFig=iFig+1; iFig
Biomass dynamic model fit to the Operating Model
growth<-FLife::vonB
dplIndex[,"54"]=0.2 sra2=fit(sra,dplIndex)
plot(as(list("SRA-2"=sra2,"SRA"=sra,"MP"=mp,"OM"=as(om,"biodyn")),"biodyns"))
Figure r iFig=iFig+1; iFig
Biomass dynamic model fit to the Operating Model
r version$version.string
r packageVersion('FLCore')
r packageVersion('FLasher')
r date()
This document is licensed under the Creative Commons Attribution-ShareAlike 4.0 International license.
Laurence KELL. laurie@seaplusplus.co.uk
This vignette and 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.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.