Introduction

This vignette, conditions an Operating Model based on life history relationships. Parameters are extracted from FishBase using the FishLife and SPMpriors packages, and the Operating Model built using the FLife.

The code is generic and will work for any species. However, conditioning an Operating Model requires expert local knowledge and the input of stakeholders before it should be used to provide advice.

The life history parameters used are for Lake Tanganyika sprat (Stolothrissa tanganicae). Sprat spawns several times a year, annual mortality is high and so it has a short life span. As juveniles grow, they move progressively closer to shore before moving offshore and start to appear in pelagic catches.

In this vignette for demonstration the Operating Model is age-structured and is
annual. Another vignette develops a seasonal Operating Model, with multiple spawning episodes per year. Each spawning episode can generate a “subcohort” which can have different growth, maturity, and natural mortality.

The Operating Model is used to generate pseudo data, for simulation testing, by the Observation Error Model. Data sets simulated are catch ($C_t$), an index of relative abundance ($I_t$), and length frequency data.

These data are used to develop length based-indicators (LBIs) and an index for relative harvest rate ($H_t = C_t/I_t$).

To identify which indicators can track changes in stock status, the simulated inducators are compared to the “true” values from the Operating Model, i.e. by comparing

Relative index of abundance $I_t$ verses vulnerable biomass $V_t$
length based-indicators such as $L_{mean}$ verses Fishing Mortality $F_t$
relative harvest rate $H_t$ verses $F_t$

The best indicators can then be used in an empirical Harvest Control Rule (HCR).

Packages

Package required

FLR

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

library(FLCandy)

FishLife

library(rfishbase)
library(SPMpriors)
library(FishLife)

Plotting and Data wrangling

library(plyr)
library(dplyr)
library(reshape)

library(ggplot2)
theme_set(theme_bw())

library(ggpubr)
library(ggcorrplot)
library(GGally)

Help on individual functions is available, e.g.

?flmvn_traits

Operating Model Conditioning

Life History Parameters

Life history parameters are used to condition the Operating Model, these are the Von Bertalanffy (1957) growth parameters ($k$ and $L_{\infty}$), age at maturity ($a\_{50}$, natural mortality ($M$), the steepness of the stock recruitment relationship ($h$), and recruitment variability (sigR).

First extract life history parameters from Fishbase using the FishLife package

S.tanganicae=flmvn_traits(Genus="Stolothrissa",Species="tanganicae",Plot=FALSE)

     [,1]           [,2]         
[1,] "K"            "M"          
[2,] "Winfinity"    "Loo"        
[3,] "tmax"         "tm"         
[4,] "Lm"           "Temperature"
[5,] "ln_margsd"    "rho"        
[6,] "logitbound_h" "ln_r"

Create an FLPar object

par=FLPar(c(cast(S.tanganicae[[1]][,c("trait","mu.sp")],~trait)[,c("K","Loo","tm","M","h","sigR")]))

Rename parameters, so they are compatible with the FLR package naming convention

dimnames(par)[[1]]=c("k","linf","a50","m","s","sigmaR")

The correlations between parameters and their variability are modelled by the covariance matrix

cov=S.tanganicae[[2]][[1]][[1]][c("K","Loo","tm","M","h"),c("K","Loo","tm","M","h")]
dimnames(cov)=list(c("k","linf","a50","m","s"),
                   c("k","linf","a50","m","s"))

Some parameters, however, are missing, i.e. $t_0$, and can be set manually [bug in lhPar, as iy doesnt like the parameter $m$]

par=rbind(lhPar(rbind(par[c("k","linf","a50","s")],FLPar("t0"=-0.1))),par["m"])

Individuals mature at age 1 as this is an annual model

par["a50"]=1

Natural mortality rates ($M$) vary with body size and age, therefore, take the single value from FishBase and use the ‘generalized length-inverse mortality’ paradigm (GLIM, Lorenzen 2022) to model M-at-age, i.e.

$log(M)= m_1 + m_2 \times log(L/L_{50}_{\infty}) + m_3 \times log(k)$

par["m2"]=-1.30
par["m3"]= 1.08

m1 setb the average levl of M so scale m1 so that $M$ at $l_{50}$ is equal to the ‘FishBase’ value

$ m_1 = log(M) - m_2log(L_{50}/L_{}) - m_3 log(k)$

mGlim<-function(object,par) {
  #M= exp(a + b ln(L/L∞) + c*lnK)
  exp(par["m1"]%+%(par["m2"]%*%log(object%/%par["linf"]))%+%(par["m3"]%*%log(par["k"])))}

par["m1"]=log(par["m"])%-%(par["m2"]%*%log(par["l50"]%/%par["linf"]))%-%(par["m3"]%*%log(par["k"]))

An object of class "FLPar"
params
     linf         k        t0         a         b       l50       a50     ato95 
   9.9650    2.1581   -0.1000    0.0003    3.0000    7.8310    1.0000    1.0000 
     asym        bg        m1        m2        m3         s         v      sel1 
   1.0000    3.0000    0.0608   -1.3000    1.0800    0.9141 1000.0000    1.6141 
     sel2      sel3         m 
   1.0000 5000.0000    3.3364 
units:  NA

Plot the vectors-at-age for growth, maturity and matural mortality, …

Figure 1 Vectors of quantities-at-age.

Selectivity and maturity-at-age

Figure 2 Selectivity and maturity-at-age.

Monte Carlo Simulation

Use the covariance matrix to simulate uncertainty

cov

            k      linf      a50        m         s
k     0.02900 -0.003638 -0.02244  0.01385 -0.001889
linf -0.00364  0.015468 -0.00754  0.00102 -0.000969
a50  -0.02244 -0.007536  0.14871 -0.03597 -0.002483
m     0.01385  0.001024 -0.03597  0.07366  0.010097
s    -0.00189 -0.000969 -0.00248  0.01010  0.008367

Simulate uncertainty for $k$ using the variance

nits=101
CV  =0.20

par=rlmvn(nits,par[dimnames(cov)[[1]]],diag(cov)^0.5*c(0.1,0,0,0,0),cov2cor(cov))

Figure 3 Monte Carlo simulation of $k$.

Use lhPar to generate missing parameters values

par=rbind(lhPar(rbind(par[c("k","linf","a50","s")],propagate(FLPar("t0"=-0.1),nits))),par["m"])

Set up Lorenzen (2022) $M$ model.

par["m1"]= 0.28
par["m2"]=-1.30
par["m3"]= 1.08
par["m1"]=log(par["m"])%-%(par["m2"]%*%log(par["l50"]%/%par["linf"]))%-%(par["m3"]%*%log(par["k"]))

Equilibrium Dynamics

Create an https://flr-project.org/doc/Reference_points_for_fisheries_management_with_FLBRP.html object t0 models the equilibrium dynamics

mGlimFn=function(object,par) {
  #lnM= a + b ln(L/L∞) + c*lnK
  len=wt2len(stock.wt(object),par)
  
  exp(par["m1"]%+%(par["m2"]%*%log(len%/%par["linf"]))%+%(par["m3"]%*%log(par["k"])))}

Create an FLBRP object with natural mortality based on Lorenzen.

eq=lhEql(par,m=mGlimFn)

Figure 4 Production functions.

Time Series

Simulate projections using FLasher

Simulate a stock which is exploited sustainably, i.e. by taking $M$ as a proxy for $F_{MSY}$

fbar(eq)[]=par["m"]

Convert the FLBRP object to an FLStock object

om=as(brp(eq),"FLStock")

Use FLasher to make projections

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

Simulate trends in $F$

trnd=FLQuant(c(t(maply(c(seq(2.5,1.0,length.out=51),
                         seq(1,  0.5,length.out=51)[-1]), function(x)
  seq(1,x,length.out=21)))),dimnames=list(year=seq(21),iter=seq(101)))

Plot using ggplotFL

plot(trnd,iter=c(20,70))

Figure 5 Simulated trends in $F$, condidence intervals, with two example trajectories (red & blue).

Simulate a stock

omTrend=fwd(om,f=fbar(om)[,ac(51:71)]%*%trnd,sr=eq)

Plot

plot(window(omTrend,start=50,end=70), iter=c(20,70))

[1] "maxfbar has been changed to accomodate new plusgroup"

Simulate trends and a sudden change in $F$

kick=FLQuant(rep(c(trnd[,21]),each=21),dimnames=list(year=seq(21),iter=seq(101)))

oms=FLStocks("Constant"=om)
oms[["Trend"]]   =fwd(om,f=fbar(om)[,ac(51:71)]%*%trnd,sr=eq)
oms[["Kick"]]    =fwd(om,f=fbar(om)[,ac(51:71)]%*%kick,sr=eq)
oms=FLStocks(llply(oms, window, start=31, end=71))

[1] "maxfbar has been changed to accomodate new plusgroup"
[1] "maxfbar has been changed to accomodate new plusgroup"
[1] "maxfbar has been changed to accomodate new plusgroup"

Figure 6 Operating Model, for the changes in fishing mortality.

Indicators

Figure 7 Indices generated by the Observation Error Model, for trends in fishing mortality.

Classification Skill

Truncate the time series to the region of interest, i.e. where $F$ changes.

Create an index of abundance with measurement error, with a CV of 30%

Combine into an FLQuants object

and plot

plot(flqs, iter=c(17,23))

Figure 8 Time series of SSB and relative index of abundance.

source("/home/laurie/Desktop/active/startup.R")
sv(oms)

Compare

Figure 9 plot of index of abundance verses SSB.

More Information

For more information on the FLR Project for Quantitative Fisheries Science in R, visit the FLR webpage ¹.

Author information

Laurence Kell. laurie@seaplusplus.co.uk

Acknowledgements

Software Versions

R version 4.2.1 (2022-06-23)

knitr: 1.43
FLCore: 2.6.19.9105
FLBRP: 2.5.9.9022
FLasher: 0.7.1.9221
FLife: 3.4.0
ggplotFL: 2.7.0.9132
FLCandy: 0.1.0
rfishbase: 3.1.9.99
SPMpriors: 1.0.4
FishLife: 3.0.0
plyr: 1.8.9
dplyr: 1.1.4
reshape: 0.8.9
ggplot2: 3.4.4
ggpubr: 0.6.0
ggcorrplot: 0.1.3
GGally: 2.1.2

Compiled: Mon Jan 22 13:58:22 2024

References

Lorenzen, Kai. 2022. “Size-and Age-Dependent Natural Mortality in Fish Populations: Biology, Models, Implications, and a Generalized Length-Inverse Mortality Paradigm.” Fisheries Research 255: 106454.

Von Bertalanffy, L. 1957. “Quantitative Laws in Metabolism and Growth.” Quarterly Review of Biology, 217–31.

http://flr-project.org ↩︎

01: Operating Model Conditioning

FAO: Small Pelagics Fisheries

L Kell

22 January, 2024

Introduction

Packages

FLR

FishLife

Plotting and Data wrangling

Operating Model Conditioning

Life History Parameters

Monte Carlo Simulation

Equilibrium Dynamics

Time Series

Indicators

Classification Skill

More Information

Author information

Acknowledgements

Software Versions

References

References