FLife

The FLife package is used to create a stock.

turbot=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")

The lhPar method is then used to derive the parameters for natural mortality-at=age, based on @gislason2008does, and default parameters and relationships for selection pattern and stock recruitment.

turbot=lhPar(turbot)

turbot

An object of class "FLPar"
params
     linf         k        t0         a         b     ato95       a50 
  59.1000    0.2800   -0.4000    0.0111    3.1500    1.0000    4.0000 
     asym        bg        a1        sl        sr         s         v 
   1.0000    3.1500    5.0000    1.0000 5000.0000    0.9000 1000.0000 
      l50 
  43.2500 
units:  cm 

The default parameters can be changed, e.g. by changing a parameter. sl is the standard deviation for the lefthand limb of the double normal selection pattern, here we change it from 2 to 1 to make it steeper.

Equilibrium dynamics

The parameters are then used by lhEql to simulate the equilibrium dynamics by combining the spawner/yield per recruit relationships with a stock recruiment relationship.

eq=lhEql(turbot)

range(eq)[c("minfbar","maxfbar")]=ceiling(mean(turbot["a1"]))

Figure 1 Vectors of m, selection pattern, maturity and weight-at-age.

Estimate equilibrium dynamics and reference points

Figure 2 Expected, equilibrium, dynamics and reference points.

Time series

To go from equilibrium to time series dynamics the FLBRP object created by lhEql can be coerced to an FLStock object.

First change the F time series so that it represents a time series where the stock was origionally lightly exploited, F increased until the stock was overfished and then fishing pressure was reduced to ensure spawning stock biomass was greater than \(B_{MSY}\).

fbar(eq)=refpts(eq)["msy","harvest"]%*%FLQuant(c(rep(.1,19),
                                              seq(.1,2,length.out = 30),
                                              seq(2,1.0,length.out = 10),
                                              rep(1.0,61)))[,1:105]
om=as(eq,"FLStock")

om=fwd(om,fbar=fbar(om)[,-1],sr=eq)

Stochastic dynamics

To simulation random variation in the time series, deviations around the stock recruitment relationship was modelled as a random variable.

While to generate data for use in the MP, random measurement error was added to the simulated catch per unit effort (CPUE).

nits=10
set.seed(3321)
uDev =rlnorm(nits,setPlusGroup(stock.n(eq),20)*0,.2)

These deviates were then used to create a stochastic time series by projecting the dynamics from year 1.

set.seed(1234)
srDev=rlnoise(nits,fbar(om)%=%0,.3,b=0.0)

Figure 3 Time series of recruitment deviates

om =propagate(om,nits)
oms=FLStocks("Projection"=fwd(om,fbar=fbar(om)[,-1],residuals=srDev,sr=eq))

Figure 4 Stochastic Time series of F, SSB, recruitment and yield

Back to Top

Feedback control

Management of a fish stocks is done using feedback control. The stock is assessed using historical data which is used estimate current stock status and then to project the stock forward under alternative management regulations for a variety of hypotheses and system dynamics. This procedure is then repeated in subsequent year to monitor and adjust the impact of management. MSE does this my simulating a MP. These can either be model based or empirical, i.e. based on a stock assessment or data alone.

In the mpb package there are a variety of MP, e.g. age, biomass and empirical based.

library(kobe)

hcr= data.frame(stock  =c(0.0 ,0.1 , 0.6,2.0), 
                harvest=c(0.01,0.01, 0.7,0.7))
kobePhase()+
  geom_line(aes(stock,harvest),data=hcr,col="orange",size=2)

Age Based

In this example the MP is based on an Virtual Population Analysis (VPA).

First the control settings are checked by running FLXSA on data simulated by the OM without error and feedback. Ideally there should be no bias in the estimates from the stock assessment

library(FLXSA)
mp=window(setPlusGroup(oms[["Projection"]],20),end=80)

xsaControl=FLXSA.control(tol    =1e-09, maxit   =150, 
                         min.nse=0.3,   fse     =1.0, 
                         rage   =1,     qage    =6, 
                         shk.n  =TRUE,  shk.f   =TRUE, 
                         shk.yrs=1,     shk.ages=4, 
                         window =10,    tsrange =10, 
                         tspower= 0,
                         vpa    =FALSE)

idx=FLIndex(index=stock.n(mp)%*%uDev[,dimnames(stock.n(mp))$year])
range(idx)[c("plusgroup","startf","endf")]=c(NA,0.1,.2)

xsa=FLXSA(mp,idx,
          control=xsaControl,diag.flag=FALSE)
range(xsa)[c("min","max","plusgroup")]=range(mp)[c("min","max","plusgroup")]
mp=mp+xsa

sr=fmle(as.FLSR(mp,model="bevholt"),control=list(silent=TRUE))
rf=FLBRP(mp,sr)

plot(FLStocks("Stock\nAssessment"=mp,
              "Operating\nModel" =window(oms[["Projection"]],end=80)))

Before running the MSE, i.e. using XSA as part of a feedback control procedure, the current reference points need to be estimated.

Then the MSE can be run using the mseXSA function

oms["Age"]=mseXSA(oms[["Projection"]],eq, #OM
                    mp,control,rf=rf,        #MP
                    srDev=srDev,uDev=uDev,   #Random deviates for OM 
                    start=75,end=103,maxF=1.0)        #year range

Figure 6 Time series from the MSE of F, SSB, recruitment and yield

Biomass Based

In mpb there is a biomass dynamic stock assessment, designed to be used as an MP.

First the control object has to be set, i.e. setting best guess, bounds and any priors for parameters.

Then the assessment is run without feedback

and compared to the OM

source('~/Desktop/flr/mpb/R/setMP.R')

## MP
mp=setMP(as(window(om,end=75),"biodyn"),
         r =   0.25,
         k =1000.0,
         b0=   0.9,
         p =  -0.6)

Figure 7 Comparision of estimates and simulated time series of harvest rate and stock biomass.

save(om,eq,mp,file="/home/laurence/Desktop/tmp/om.RData")

library(FLasher)
library(mpb)

source('~/Desktop/sea++/mydas/pkg/R/hcr.R')
source('~/Desktop/sea++/mydas/pkg/R/mseMPB.R')

nits=dims(om)$iter
set.seed(1234)
srDev =rlnoise(nits,FLQuant(0,dimnames=dimnames(iter(catch(om),1))),0.3,b=0.0)
uDev  =rlnoise(nits,FLQuant(0,dimnames=dimnames(iter(catch(om),1))),0.2,b=0.0)
selDev=rlnoise(nits,FLQuant(0,dimnames=dimnames(iter(    m(om),1))),0.2,b=0.0)

eq=FLCore::iter(eq,seq(nits))

oms["Biomass"]=mseMPB2(om,eq,mp,start=75,end=100,ftar=0.5,srDev=srDev,uDev=uDev,selDev=selDev)

==75, 

78, 

81, 

84, 

87, 

90, 

93, 

96, 

==

Figure 8 Time series from the MSE of F, SSB, recruitment and yield

Empirical

save(om,eq,control,srDev,uDev,file="/home/laurence/Desktop/tmp/test.RData")

source('~/Desktop/sea++/mydas/pkg/R/hcrSBTD.R')
source('~/Desktop/sea++/mydas/pkg/R/mseSBTD.R')

control=FLPar(k1=0.5,k2=0.5,gamma=1)
oms["Emprirical"]=mseSBTD(om,eq,control=control,
                            srDev=srDev,uDev=uDev,
                            start=75,end=100)

==75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, ==

Figure 9 Time series from the MSE of F, SSB, recruitment and yield

res=ldply(oms,omSmry,eq)

Back to Top

Software Versions

• R version 3.5.1 (2018-07-02)
• FLCore: 2.6.10.9004
• FLPKG:
• Compiled: Mon Dec 10 10:12:15 2018
• Git Hash: c308b3b

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.