Management Strategy Evaluation (MSE) is used to evaluate alternative management plans using computer simulation. This requires the use of an Operating Model (OM) to simulate the actual system to be managed and to evaluate the performance of alternative candidate Management Procedures (MPs) to be applied in practice. Testing an MP should essentially be blind experiments where information about the system is limited to the data available to the stock assessment. Performance statistics based on the OM (e.g. yield, probability of stock collapse) are then used to evaluate the performance of the management against the management objectives.
Here we document a generic single species Management Strategy Evaluation where the Operating Model is the currently-used stock assessment model. The use of the assessment model as the operating model implies that the assessment models describe nature almost perfectly. However, this is an important test of a Management Procedure since if it cannot perform well when reality is as simple as that implied by an assessment model, it is unlikely to perform adequately for more realistic representations of uncertainty.
Basing an OM on the current assessment model also arguably has the lowest demands for knowledge and data and also ensures that advice based upon the MSE is consistent with current advice, potentially of importance for acceptance when developing recovery and management plans.
The OM is a mathematical-statistical model used to describe the underlying resource dynamics and to generate future data when projecting forward. The purpose of the Operating Model is not to reproduce the entire complexity of both species biology and ecology, but to define plausible hypothesis about population dynamics, then to evaluate consequences for exploitation.
The OM was based on the age-structured equation using the approach of Fromentin and Kell (2007) i.e.
N(a,t) = N(a-1,t-1)exp(-Za,t)
where N(a,t) are the number of fish of age a at time t, and Z(a,t) is the total mortality from age a–1 to age a, i.e. Z(a,t) = M(a,t) + F(a,t) where M(a,t) is the natural mortality at age a and F(a,t) is the fishing mortality at age a in year t.
The population dynamics of the OM is also based on a stock recruitment relationship to allow N in the first age to be simulated.
Operating models must be ‘conditioned’ on data, i.e. fitted to the data so that the model predictions of the data are approximately consistent with the observations. This conditioning process can lead to an undesirably narrow range of scenarios, with the result that candidate OMs are only tested against scenarios which are consistent with the stock assessment assumptions or are at least fairly likely given the observed data.
However, conditioning should not necessarily lead to a narrow range of scenarios since the MPs will be used to manage a stock in the future and should therefore be tested for problematic cases which are possible and may with hindsight prove the assumptions used within the stock assessment to be wrong, i.e. scenarios represent ‘justified concerns’ to which the MP should be robust.
The data used in this example are from the example ple4 data set
library(FLCore) data(ple4)
The stock recruitment relationship was fitted using FLSR part of FLCore
## Stock recruitment relationship # create an FLSR object ple4.sr <-as.FLSR(ple4) # scale to help in fitting ple4.sr <-transform(ple4.sr,ssb=ssb/1000,rec=rec/1000) # fit a Ricker relationship model(ple4.sr)<-ricker() ple4.sr <-mle(ple4.sr) # inspect fit and diagnostics plot(ple4.sr)
Often when performing a stock assessment using virtual population analysis (VPA) a retrospective pattern will be seen. While looking for such patterns is a commonly done in stock assessment working groups determining what causes such patterns is not so easy. They can be caused by changes in population characteristics that are assumed to be stable over time, such as selectivity, catchability, or natural mortality. But they may also be caused by conflicting signals coming in from commercial or survey indices. Such a bias can be important in affecting the performance of a harvest control rule (HCR) if it is also a function of trends in the data and if the stock assessment attempts to correct it.
As well as conducting a stock assessment FLR can be used to perform retrospective analyses
library(FLXSA) ## Performing a stock assessment data(ple4.indices) ple4<-ple4+FLXSA(ple4,ple4.indices,FLXSA.control()) ## performing a retrospective analysis retro.ple4<-retro(ple4,ple4.indices,FLXSA.control(),retro=10) xyplot(data~year,groups=qname,data=lapply(retro.ple4,fbar), type="l") xyplot(data~year,groups=qname,data=lapply(retro.ple4,"ssb"),type="l")
XSA to reduce retrospective patterns employs shrinkage, i.e. the estimate of fishing mortality in the most recent year includes a weighting equal to the recent average fishing mortlaity. However, changing fishing mortality is the whole point of a HCR and this can lead to bias in it performance.
To illustrate this effect, FLXSA was run with both low shrinkage (fse=0.75) and high shrinkage (fse=0.25), where fse is the assume standard error in F used to determine the inverse variance weighting.
## performing a retrospective analysis with high shrinkage retro.ple4<-retro(ple4,ple4.indices,FLXSA.control(fse=0.25),retro=10) xyplot(data~year,groups=qname,data=lapply(retro.ple4,fbar), type="l")
## performing a retrospective analysis with low shrinkage retro.ple4<-retro(ple4,ple4.indices,FLXSA.control(fse=0.75),retro=10) xyplot(data~year,groups=qname,data=lapply(retro.ple4,fbar), type="l")
If the bias is produced by a trend in selectivity, catchability, or natural mortality then assuming a trend of the opposite sign will also act to reduce the retrospective pattern. See the example where a bias was assumed in catchability.
indices<-ple4.indices ## function to apply a linearly increasing trend to an FLQuant retro.adjust<-function(x,obj) { if (x>0) res <-1-(sort(cumsum(rep(x, dims(obj)$year)),d=T)-x) else res <-sort(1-(sort(cumsum(rep(x, dims(obj)$year)),d=T)-x),d=T) return(obj*FLQuant(rep(res,each=dims(obj)$age),dimnames=dimnames(obj))) } ## apply a 5% increase in catcahability to allow indices indices[[1]]@index<-retro.adjust(-0.05,ple4.indices[[1]]@index) indices[[2]]@index<-retro.adjust(-0.05,ple4.indices[[2]]@index) retro.ple4<-retro(ple4,indices,FLXSA.control(),retro=10) xyplot(data~year,groups=qname,data=lapply(retro.ple4,fbar), type="l")
However, in practice distinguishing between potential causes of bias will be difficult. Therefore alternative scenarios should be run to test that the management advice is robust to what we don’t know.
The OEM essentially provides an interface between the “true world” of the operating model, and the “perceived world” of the management procedure by generating the types of data that would be available to the management procedure in a realistic manner (i.e. incorporating the errors, with appropriate statistical properties, that would be encountered in practice).
In order to estimate numbers and fishing mortality-at-age using VPA, estimates of catch-at-age and indices of abundance are required.
## set up objects nits <-100 # number of iterations now <-2001 # current year, i.e. starting year for MSE hist.yr <-1957 future.yr <-now+20 # final project year # create an object with mulitple iterations ple4 <-propagate(stf(ple4,nyrs=21),iter=nits) # create an FLBiol object, representing the biological population dynamics from the # original stock assessment biol <-as(ple4,"FLBiol") # create an FLFleet object, representing the fleet dynamics from the # original stock assessment fleet <-as(ple4,"FLFleet") fleet@capacity[] <-0.75
This is “arse about tit” since the values used in the MP should be sampled from the OM for use in the MP.
Values in the MP would include appropriate bias and noise associated with the data collection process. This is most easily done by creating a function to copy quantities that can be observed from the OM for use in the MP.
For example indices of abundance, as surveyed by a research cruise, or catch-at-age as sampled by observers or at a port of landing.
A function that assumes the index is an unbiased estimate of numbers-at-age (e.g. as given by a research survey) by copying numbers-at-age from the OM and adding a random term equivalent to measurement error.
## Research survey indices ## deviates dmns<-list(age=1:plusgroup,year=hist.yr:future.yr,iter=1:nits) deviates.index <-FLQuant(exp(rnorm( prod(unlist(lapply(dmns,length))), mean=0, sd=.30)),dimnames=dmns) OEMCPUE<-function(biol,deviates,start,end="missing",plusgroup="missing") { if (missing(end)) end<-start yrs<-start:end if (!missing(plusgroup)) index<-setPlusGroup(trim(biol,year=yrs),plusgroup)@n else index<-trim(biol,year=yrs)@n ## unbiased population estimates return(index*deviates[,ac(yrs)]) } # creating an index for historic years index<-OEMCPUE(biol,deviates.index,start=hist.yr,end=now-1,plusgroup=plusgroup) index<-FLIndex(index=index,startf=0,endf=.1) # creating an index for current year index<-window(index,end=now) index@index[,ac(now)]<-OEMCPUE(biol,deviates.index,now,plusgroup=plusgroup) # adding a bias consistent with the retrospective bias?
As well as catch-at-age data assumptions about mass, maturity and natural mortality-at-age are required, these can be modelled by an OEM again with noise and bias equivalent to the sampling regimes in operation. Although in this case no noise or bias has been added.
## Observation Error Model OEMStock<-function(biol,fleet,start,end="missing",plusgroup="missing") { if (missing(end)) end<-start yrs<-as.character(end:start) stk <-as(window(fleet@metiers[[1]]@catches[[1]], start=start, end=end),"FLStock") ## biological parameters stock.wt( stk)[,yrs] <-wt( biol)[dimnames(stk@m)$age,yrs] m( stk)[,yrs] <-m( biol)[dimnames(stk@m)$age,yrs] mat( stk)[,yrs] <-fec( biol)[dimnames(stk@m)$age,yrs] harvest.spwn(stk)[,yrs] <-spwn(biol)[dimnames(stk@m)$age,yrs] m.spwn( stk)[,yrs] <-spwn(biol)[dimnames(stk@m)$age,yrs] if (!missing(plusgroup)) stk<-setPlusGroup(stk,plusgroup) return(stk) } # data collection for stock # perfect sampling of catch and biological parameters # historic years stk <-OEMStock(biol,fleet,start=hist.yr,end=now-1,plusgroup=plusgroup) # current year stk <-window(stk,end=now) stk[,ac(now)] <-OEMStock(biol,fleet,start=now,plusgroup=plusgroup)
An MP is the combination of pre-defined data, together with an algorithm to which such data are input to provide a value for a TAC or effort control measure. The procedure will then be tested through simulation trials, to show robust performance in the presence of uncertainties.
stk<-stk[,ac(hist.yr:(now-1))]+FLXSA(stk[,ac(hist.yr:(now-1))],trim(index,year=hist.yr:(now-1)),diag.flag=FALSE)
## years for projection etc. yr.now <-as.character(stk@range["maxyear"]+1) yr.TAC <-as.character(stk@range["maxyear"]+2) yr.recruits<-as.character(stk@range["minyear"]:(stk@range["maxyear"]-1)) ## set up future years stk <-stf(stk,nyrs=2) ## ref points Fpa <-0.6 Bpa <-200000 ## F status quo Fsq <-c(fbar(stk)[,yr.now]) TAC <-50000 ## target options target<-fwdTarget(year =c(stk@range["maxyear"]-1,rep(stk@range["maxyear"],4)), quantity=c(rep( "f",2),"catch","ssb","f"), value =c(mean(Fsq),mean(pmin(Fsq,Fpa)), NA, NA, NA), min =c(NA,NA,mean(TAC*(1.0-TAC.bound)),mean(Bpa), Fmin), max =c(NA,NA,mean(TAC*(1.0+TAC.bound)), NA,mean(Fpa))) cat("\n===================","Target"," ====================\n") print(target[,1:5]) cat("================================================\n") ## targets by iter value <-array(NA,c(dim(target)[1], dims(stk)$iter)) min <-value max <-value value[1,]<-c(Fsq) value[2,]<-c(pmin(Fsq,Fpa)) value[2,]<-c(pmin(Fsq,Fpa)) min[3,]<-c(TAC*(1.0-TAC.bound)) min[4,]<-c(Bpa) min[5,]<-Fmin max[3,]<-c(TAC*(1.0-TAC.bound)) max[5,]<-c(Fpa) ## run t.<-system.time(stk<-fwd(stk,target[1:2,],value=value[1:2,], min=min[1:2,], max=max[1:2,], sr.model="mean",sr.param=exp(mean(log(stk@stock.n[1,yr.recruits]),na.rm=TRUE)))) ## get Quota TAC<-computeCatch(stk)[,yr.TAC]
Performance statistics are used to evaluate how well a specific MP achieves some or all of the general objectives for management for a particular scenario. These are based upon the OM since it is the actual performance not the perceived performance that is being used to evaluate the alternative MPs.
## Summary statistics SmryStats<-function(fleet,biol, start=-30, end=0) { fleet<-window(fleet, start=start, end=end) biol <-window(biol , start=start, end=end) cost <-function(fleet) {fleet@fcost+fleet@metiers[[1]]@vcost*fleet@effort} revenues=revenue(fleet@metiers[[1]])[[1]] costs =cost(fleet) profits =revenues-costs yield =landings(fleet,1)[[1]] catch =catch(fleet,1)[[1]] discards=discards(fleet,1)[[1]] ssb =ssb(biol) biomass =apply(biol@n*biol@wt,c(2,6),sum) mnsz =apply(biol@n*biol@wt,c(2,6),sum)/apply(biol@n,c(2,6),sum) pmat =ssb/biomass f =catch/biomass res=FLQuants(revenues=revenues, costs =costs, profits =profits, yield =yield, catch =catch, discards=discards, ssb =ssb, biomass =biomass, mnsz =mnsz, pmat =pmat, f =f) res=FLQuants(revenues=revenue(fleet@metiers[[1]])[[1]], costs =cost(fleet), profits =revenues-costs, yield =landings(fleet,1)[[1]], catch =catch(fleet,1)[[1]], discards=discards(fleet,1)[[1]], ssb =ssb(biol), biomass =apply(biol@n*biol@wt,c(2,6),sum), mnsz =apply(biol@n*biol@wt,c(2,6),sum)/apply(biol@n,c(2,6),sum), pmat =ssb/biomass, f =catch/biomass) return(res) }
SmryPlots<-function(fleet,biol,stk,fbarRange,yrs) { # fleet<-window(fleet,end=max(yrs)+2) # biol <-window(biol, end=max(yrs)+2) stk <-window(stk, end=max(yrs)+2) if (is.numeric(yrs)) yrs<-as.character(yrs) c.fleet =apply(catch.n(fleet,1,1)*catch.wt(fleet,1,1),c(2,6),sum) OM.f =quantile(apply(calcF(m(biol)[ac(fbarRange),yrs],catch.n(fleet,1,1)[ac(fbarRange),yrs],n(biol[ac(fbarRange),yrs])),c(2,6),mean),na.rm=T)[,yrs] MP.f =quantile(apply(stk@harvest[ac(fbarRange)], c(2,6),mean),na.rm=T)[,yrs] OM.rec =quantile(biol@n[1,], na.rm=T)[,yrs] MP.rec =quantile(stk@stock.n[1,], na.rm=T)[,yrs] OM.ssb =quantile(ssb(biol), na.rm=T)[,yrs] MP.ssb =quantile(ssb(stk), na.rm=T)[,yrs] OM.yield=quantile(c.fleet, na.rm=T)[,yrs] MP.yield=quantile(computeCatch(stk), na.rm=T)[,yrs] dmns <-dimnames(OM.ssb) dmns$iter<-1:10 ssb. <-FLQuant(c(OM.ssb, MP.ssb), dimnames=dmns) f. <-FLQuant(c(OM.f, MP.f), dimnames=dmns) rec. <-FLQuant(c(OM.rec, MP.rec), dimnames=dmns) yield.<-FLQuant(c(OM.yield,MP.yield),dimnames=dmns) print(xyplot(data~year|qname,groups=iter,data=FLQuants(F =f., SSB =ssb., Recruits =rec., Yield =yield.), type="l",col=c(rep("red",5),rep("blue",5)),lwd=c(1,2,3,2,.1,.1,2,3,2,1),scale="free")) }
An example of creating an OM from a WG stock assessment. The OM is created from the FLStock stock assessment object, and comprises two parts the FLBiol (corresponding to the biological dynamics) and the FLFLeet (corresponsing to the fleet dynamics).
data(ple4) ## add iterations fore Monte Carlo simulation nits<-100 ple4<-propagate(ple4,iter=nits) ## Create Operation model classes ## Fleet fleet<-as(ple4,"FLFleet") ## Biological population biol <-as(ple4,"FLBiol")
To conduct the stock assessment we need to collect data for indices of abundance and catch-at-age to perform VPA.
## Observation Error Models ## random noise corresponding to measurement error dmns<-list(age=1:plusgroup,year=biol@range["minyear"]:biol@range["maxyear"],iter=1:nits) deviates.index <-FLQuant(exp(rnorm( prod(unlist(lapply(dmns,length))), mean=0, sd=.30)),dimnames=dmns) ## CPUE index<-OEMCPUE(biol,deviates.index,start=biol@range["minyear"],end=biol@range["minyear"]) index<-FLIndex(index=index,startf=0,endf=.1) ## Catch-at-age and biological parameters stk <-OEMStock(biol,fleet,start=biol@range["minyear"],end=biol@range["minyear"]) ## Stock assessment stk <-stk+FLXSA(stk,index,diag.flag=FALSE)
## Check to see if the assessment can recreate the dynamics
##SSB xyplot(data~year,groups=qname,data=FLQuants(biol=ssb(biol),stk=quantile(window(ssb(stk),end=2001),na.rm=T)),type="l") ##F fbar.OM<-apply(calcF(biol@m,catch.n(fleet,1,1),biol@n)[ac(stk@range["minfbar"]:stk@range["maxfbar"])],2,mean) xyplot(data~year,groups=qname,data=FLQuants(biol=fbar.OM,stk=quantile(window(fbar(stk),end=2001),na.rm=T)),type="l")
library(FLCore) library(FLXSA) library(FLash) data(ple4) ## Check Retro retro.ple4<-retro(ple4,ple4.indices,FLXSA.control(),retro=10) xyplot(data~year,groups=qname,data=lapply(retro.ple4,fbar), type="l") xyplot(data~year,groups=qname,data=lapply(retro.ple4,"ssb"),type="l") ## Stock recruitment relationship ple4.sr <-as.FLSR(ple4) ple4.sr <-transform(ple4.sr,ssb=ssb/1000,rec=rec/1000) model(ple4.sr)<-ricker() ple4.sr <-mle(ple4.sr) plot(ple4.sr) sr.param <- c(params(ple4.sr))[1:2]*c(1,1/1000) sr.model <-"ricker" ## set up objects nits <-100 ple4 <-propagate(stf(ple4,nyrs=21),iter=nits) biol <-as(ple4,"FLBiol") fleet <-as(ple4,"FLFleet") fleet@capacity[] <-0.75 ## years now <-2001 future.yr <-now+19 history.yr <-dims(ple4)$minyear plusgroup <-ple4@range["plusgroup"] ## HCR options Fpa <-.60 Bpa <-200000 Fmin <-.10 TAC.bound <-0.15 GodOnlyKnows <-FALSE plusgrp <-ple4@range["plusgroup"] # Monte Carlo #sr.deviates <-FLQuant(sample(c(exp(residuals(ple4.sr))),replace=TRUE,size=nits*20*nits),dimnames=dmns) sr.deviates <-FLQuant(exp(rnorm(nits*dims(ple4)$year, mean=0.0, sd=.3)),dimnames=list(age=1,year=now:future.yr,iter=1:nits))/exp(.3^2/2.0) index.deviates <-FLQuant(exp(rnorm(dims(biol@m)$age*dims(biol@m)$year*nits, mean=0, sd=.30)), dimnames=list(age=ple4@range["min"]:plusgrp,year=history.yr:future.yr,iter=1:nits)) catch.q(fleet,1,1)<-FLQuant(1*exp(rnorm(length(c(catch.q(fleet,1)[[1]])),0,.2)),dimnames=dimnames(catch.q(fleet,1)[[1]])) ## MP ######################################################################## ## set up CPUE index index <-as(trim(biol, age=biol@range["min"]:plusgroup),"FLIndex") index@range[c("startf","endf")] <-c(0,0.1) index@index[,ac(history.yr:(now-1))]<-OEMCPUE(biol,index.deviates,start=history.yr,end=now-1,plusgroup=plusgroup) ## set up stock stk <-as(biol,"FLStock") stk <-setPlusGroup(stk,plusgrp) stk@stock.n[] <-NA stk@harvest[] <-NA units(stk@harvest)<-"f" ## set up stock yrs.assess<-ac(-30:(now-1)) stk <-OEMStock(biol,fleet,start=history.yr,end=now-1) TAC <-catch(fleet,1,1)[,ac(now)] cat(sprintf("\n===========================================\n")) for (iyr in now:future.yr) { stk@range["maxyear"] <-iyr-1 yrs.assess <-as.character((history.yr):stk@range["maxyear"]) cat(sprintf("Year =%.0f\n",iyr)) ## Observation Error Model index@index[,ac(iyr-1)]<-OEMCPUE( biol,index.deviates,iyr-1,plusgroup=plusgroup) stk[ ,ac(iyr-1)]<-OEMStock(biol,fleet, iyr-1,plusgroup=plusgroup) cat("=============","Assess!"," =============\n") if (TRUE) t.<-system.time(stk<-stk[,yrs.assess]+FLXSA(stk[,yrs.assess],trim(index,year=((history.yr):(iyr-1))),diag.flag=FALSE)) else t.<-system.time(stk<-perfectAssess(stk,biol,fleet,yrs.assess)) cat("system.time",t.,"\n") ## HCR cat("=============","HCR"," =================\n") ## years for projection etc. yr.now <-as.character(stk@range["maxyear"]+1) yr.TAC <-as.character(stk@range["maxyear"]+2) yr.recruits<-as.character((stk@range["minyear"]+1):(stk@range["maxyear"]+2)) ## set up future years stk <-stf(window(stk,end=stk@range["maxyear"]),nyrs=2) ## F status quo Fsq <-c(fbar(stk)[,yr.now]) ## target options target<-fwdTarget(year =c(stk@range["maxyear"]-1,rep(stk@range["maxyear"],4)), quantity=c(rep( "f",2),"catch","ssb","f"), value =c(mean(Fsq),mean(pmin(Fsq,Fpa)), NA, NA, NA), min =c(NA,NA,mean(TAC*(1.0-TAC.bound)),mean(Bpa), Fmin), max =c(NA,NA,mean(TAC*(1.0+TAC.bound)), NA,mean(Fpa))) cat("\n===================","Target"," ====================\n") print(target[,1:5]) cat("================================================\n") ## targets by iter value <-array(NA,c(dim(target)[1], dims(stk)$iter)) min <-value max <-value value[1,]<-c(Fsq) value[2,]<-c(pmin(Fsq,Fpa)) value[2,]<-c(pmin(Fsq,Fpa)) min[3,]<-c(TAC*(1.0-TAC.bound)) min[4,]<-c(Bpa) min[5,]<-Fmin max[3,]<-c(TAC*(1.0-TAC.bound)) max[5,]<-c(Fpa) ## run t.<-system.time(stk<-fwd(stk,target[1:2,],value=value[1:2,], min=min[1:2,], max=max[1:2,], sr.model="mean",sr.param=exp(mean(log(stk@stock.n[1,yr.recruits]),na.rm=TRUE)))) ## get Quota TAC<-computeCatch(stk)[,yr.TAC] cat("system.time",t.,"\n") ## go fish!! cat("=============","Go Fish!"," ============\n") OMcontrol <-fwdControl(year=iyr+1, fleet=1, metier=1, value=NA, min=.01, max=.75) OMtarget <-fwdTarget( year=iyr+1, value=mean(c(TAC)), quantity="catch", fleet=1, metier=1, spp=1) t.<-system.time(res<-fwd(biol, fleet, OMtarget, OMcontrol, target.value=t(matrix(TAC)),sr.model=sr.model, sr.param=sr.param, sr.residuals=sr.deviates[,as.character(iyr+1)])) landings.n(fleet,1,1)[,as.character(iyr+1)]<-res$landings.n[,as.character(iyr+1)] discards.n(fleet,1,1)[,as.character(iyr+1)]<-res$discards.n[,as.character(iyr+1)] effort(fleet)[ ,as.character(iyr+1)]<-res$effort[ ,as.character(iyr+1)] n(biol)[ ,as.character(iyr+1)]<-res$n[ ,as.character(iyr+1)] cat("system.time",t.,"\n") SmryPlots(fleet,biol,stk,stk@range[c("minfbar","maxfbar")],(history.yr):(iyr-1)) }