An example of fitting biological reference points, including stock and recruitment relationships and economic analyses.
## Load the example FLStock object from FLCore library(FLBRP) data(ple4) ## Create the corresponding FLBRP object ple4.brp<-FLBRP(ple4,nyrs=5)
nyrs specifies number of years used to calculate the mean values, e.g. wts-at-age, selection patterns etc.
the default is 3.
The FLBRP class holds several slots corresponding to the assumptions about selection patterns, biological parameters, stock recruitment relationships and economics
## Selection Patterns catch.sel(ple4.brp) discards.sel(ple4.brp) bycatch.harvest(ple4.brp) xyplot(data~age,data=catch.sel(ple4.brp),type="l") ## Mass-at-age stock.wt(ple4.brp) catch.wt(ple4.brp) discards.wt(ple4.brp) bycatch.wt(ple4.brp) xyplot(data~age,groups=qname,data=FLQuants(swt=stock.wt(ple4.brp),cwt=catch.wt(ple4.brp)),type="l") ## Biological parameters m(ple4.brp) mat(ple4.brp) xyplot(data~age|qname,data=FLQuants(sel=catch.sel(ple4.brp),dsel=discards.sel(ple4.brp), swt=stock.wt(ple4.brp), cwt =catch.wt(ple4.brp), mat= mat(ple4.brp), m = m(ple4.brp)) ,type="l",scale="free")
There are also slots corresponding to the historic observations
## historic observations harvest.obs(ple4.brp) yield.obs(ple4.brp) r.obs(ple4.brp) ssb.obs(ple4.brp) profit.obs(ple4.brp)
Once an FLBRP object has been created then equilibrium quantities at age and time series can be estimated
## estimate equilibrium quantities ## fmult gives the scaling of the selection patterns to generate F fmult=1:10*.2 ple4.brp<-brp(ple4.brp,fmult=fmult) ## Fishing mortality harvest(ple4.brp) ## check that it worked! sweep(harvest(ple4.brp),1,catch.sel(ple4.brp),"/") ## other estimated quantities-at-age stock.n(ple4.brp) catch.n(ple4.brp) landings.n(ple4.brp) discards.n(ple4.brp)
Time series
## estimates of time series yield.hat(ple4.brp) r.hat(ple4.brp) ssb.hat(ple4.brp) profit.hat(ple4.brp)
As well as calculating equilibrium quantities you can also calculate reference points based upon these. For example theyield and spawner per recruit reference points F0.1, Fmax.
## get the slots that will hold the reference points refpts(ple4.brp)
Fmsy is the same as Fmax, since the default assumed stock recruitment relationship is mean recruitment. Plotting the reference points ans expected quantities
## plot per recruit values ple4.brp<-brp(FLBRP(ple4)) plot(ple4.brp) ## plot against historical data # first set recruitment to mean ple4.brp@sr.params<-mean(ple4.brp@r.obs) # redo brp ple4.brp<-brp(ple4.brp) plot(brp(ple4.brp),obs=TRUE) ## Two sex example data(ple4sex) brp(FLBRP(ple4sex))@refpts ## changing slots when creating new.mat<-FLQuant(c(rep(0,4),rep(1,11)),dimnames=dimnames(ple4.brp@mat))
Adding new reference points, these can be calculated by adding the quantity (e.g. F, Yield, ssb)
## adding new reference points new.brp<-FLBRP(ple4,harvest=c(Fpa=0.6)) new.brp<-brp(new.brp) refpts(new.brp)
For an economic analysis, price-at-age and costs (both fixed and variable) are required. Prices are easy to get, but costs are more difficult.
## adding price ple4.brp@price[]<-c(1.81,1.81,1.81,1.93, 1.93, 2.39, 2.39, 2.39, 2.39, 3.42, 3.42, 3.42, 3.42, 3.42, 3.42) #For costs, we assume that costs = revenue at a particular level of F, i.e. the breakeven point. ## adding costs # assume that the break even point is at 4 times Fmsy fmult <-refpts(ple4.brp)["msy","harvest"]*4 ple4.brp <-brp(ple4.brp,fmult=fmult) revenue.brkeven<-apply(sweep(landings.n(ple4.brp),1,price(ple4.brp),"*"),2,sum)[1,1,drop=T] ## At breakeven costs=revenue ## assume that fixed costs are 10% of revenue at beakeven cost<-c(fixed =revenue.brkeven*.1, variable=revenue.brkeven*.9) ## recalculate reference points ple4.brp<-FLBRP(ple4) ple4.brp@cost <-cost price(ple4.brp)[] <-price ple4.brp<-brp(ple4.brp) plot(ple4.brp)
As well as yield/spawner per recruit estimates we can include stock recruitment assumptions. For example fit a Beverton and Holt realtionship and compare to the yield recruit estimates
## pass parameters to the FLBRP object ple4.brp@sr.params <-exp(mean(log(r.obs(ple4.brp)))) ple4.brp <-brp(ple4.brp) plot(ple4.brp) ## Stock recruitment relationships plot(ssb.obs(ple4.brp),r.obs(ple4.brp),xlim=c(0,max(ssb.obs(ple4.brp))),ylim=c(0,max(r.obs(ple4.brp)))) lines(lowess(ssb.obs(ple4.brp),r.obs(ple4.brp))) ple4.sr<-as.FLSR(ple4) ## Fit a Ricker stock recruitment relationship model(ple4.sr)<-ricker() ple4.sr<-transform(ple4.sr,ssb=ssb/1000,rec=rec/1000) ple4.sr<-mle(ple4.sr) plot(ple4.sr) ple4.brp.rk <-ple4.brp ple4.brp.rk@sr.model <-"ricker" ple4.brp.rk@sr.params<-c(params(ple4.sr)[1,1:2]) ple4.brp.rk <-brp(ple4.brp.rk) ## fit a Beverton and Holt stock recruitment relationship ple4.sr <-as.FLSR(ple4) ple4.sr <-transform(ple4.sr,ssb=ssb/1000,rec=rec/1000) model(ple4.sr)<-bevholt() ple4.sr <-mle(ple4.sr) ple4.brp.bh <-ple4.brp ple4.brp.bh@sr.model <-"bevholt" ple4.brp.bh@sr.params<-c(params(ple4.sr)[1,1:2]) ple4.brp.bh <-brp(ple4.brp.bh) ple4.brp <-brp(ple4.brp) ## reference points refpts(ple4.brp) refpts(ple4.brp.rk) refpts(ple4.brp.bh) ## surplus production curves plot( ple4.brp@ssb.hat/1000,ple4.brp@yield.hat/1000,type="l",ylim=c(0,150)) lines(ple4.brp.rk@ssb.hat,ple4.brp.rk@yield.hat,type="l",col="red") lines(ple4.brp.bh@ssb.hat,ple4.brp.bh@yield.hat,type="l",col="blue")
Also works for other dims e.g. for sex
## Two sex example data(ple4sex) brp(FLBRP(ple4sex))@refpts
Can also modify slots and parametesr
## changing slots when creating new.mat<-FLQuant(c(rep(0,4),rep(1,11)),dimnames=dimnames(ple4.brp@mat)) ple4.brp.newmat<-brp(FLBRP(ple4,mat=new.mat,harvest=c(Fpa=0.6))) ## MSY refpts(brp(ple4.brp))["msy",] ## Global MSY ple4.brp@catch.sel[1:9,]<-0 refpts(brp(ple4.brp))["msy",]
Or create new reference points
## adding new reference points brp(FLBRP(ple4,harvest=c(Fpa=0.6)))
In Europe there is now a move towards long-term fishery-based plans to bring all major fish stocks to rates of fishing at which MSY can be achieved. The previous framework was reactionary and based upon limit and precautionary reference points. Which aimed to ensure that spawning stock biomass (SSB) remains above a threshold value (Blim) where either recruitment is impaired or the dynamics are unknown, and that fishing mortality remains below a threshold (Flim) that would drive the stock to Blim? Precautionary reference points (Bpa and Fpa) that take into account uncertainty were used to actually trigger management action.
In some recent MSE work Fpa was then defined as a mulitple of F0.1 (i.e. 3 times F0,1) and Bpa defined as the SSB at Fpa. The following example shows how this can be implemented as part of an MSE using FLStock & FLBRP
## example of a routine to implement rule to calculate Fpa & Bpa calc.pa<-function(stk,pa.mult=3) { ## create FLBRP object res <-FLBRP(stk) ## assume geometric mean res@sr.params<-c(alpha=exp(mean(log(stk@stock.n[1,drop=T]),na.rm=T))) res@sr.model <-"mean" ## calculate reference points etc. ## set fmult=1, to save execution time res <-brp(res,fmult=1) ##PA ref pts: FPA=F0.1*mult, BPA = B(FPA) res@refpts <-rbind(res@refpts,pa=c(res@refpts["f0.1","harvest"]*pa.mult,NA,NA,NA,NA)) ##PA ref pts: Bpa=0.5BMSY, FPA= F that gives that SSB #res@refpts <-rbind(res@refpts,pa=c(NA,NA,res@refpts[1,"r"],res@refpts[1,"ssb"]/2,NA)) res <-brp(res) return(refpts(res)["pa",c("harvest","ssb")]) }