This vignette shows how to condition a seasonal Operating Model, i.e. one with less than annual time steps that are less than annual, based on life-history theory.

References

Packages

Some ‘FLR’ packages are required

library(FLCore)
library(ggplotFL)

library(FLBRP)
library(FLasher)

library(mydas)

library(FLRef)

as well as R packages for plotting

library(ggplot2)
library(GGally)
library(ggpubr)

and data manipulation

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

Quick Start

The FLife packge is required for modelling life-histories, and ‘mydas’ for estimation and developing harvest control rules

The simplest way to obtain the FLR pacakges are to install them from the FLR repository

install.packages("FLife", repos = "http://flr-project.org/R")

See help(install.packages) for more details.

After installation, they can be loaded

library(FLife)
library(mydas)

Operating Model Conditioning

Annual model

Conditioned an Operating Models on life-history parameters allows more flexibility than using a data-rich stock assessment, that may never have been formally validated.

There are relationships between the life-history parameters, e.g. the von Bertalanfy growth parameters \(L_{\infty}\), \(k\) and \(t_0\), length at which \(50\%\) of individuals reach maturity (\(l_{50}\)), and of the length weight relationship \(W=aL^b\).

There is an example data set with teleost life-history parameters in FLife

An object of class "FLPar"
iters:  145 

params
               linf                   k                  t0                 l50 
45.100000(28.02114)  0.246667( 0.17297) -0.143333( 0.13590) 22.100000(11.71254) 
                  a                   b 
 0.011865( 0.00776)  3.010000( 0.15271) 
units:  NA

Figure 1 relationships between life-history parameters

If only the maximum size of a species is known then you can derive the other parameters from the life-history relationships.

par = lhPar(FLPar(c(linf = 59.1), units = "NA"))

par

An object of class "FLPar"
params
     linf         k        t0         a         b       l50       a50     ato95 
  59.1000    0.2315   -1.3667    0.0003    3.0000   31.9822    1.9989    1.0000 
     asym        bg        m1        m2        m3         s         v      sel1 
   1.0000    3.0000    0.5500   -1.6100    1.4400    0.9000 1000.0000    2.9989 
     sel2      sel3 
   1.0000 5000.0000 
units:  NA

Equilibrium dynamics

The life-history parameters can be used to create the equilibrium dynamics using the FLBRP class.

eq = lhEql(par)

plot(eq, refpts = "msy")

Figure 1 Equilibrium dynamics with \(MSY\) reference points.

Time series dynamics

For the Operating Model, the equilibrium FLBRP object needs to be converted into an FLStock object.

As an example, time series are simulated for fishing mortality, where an unexploited stock is fished until it becomes over exploited, after which a recovery plan is implemented.

fbar(eq) = refpts(eq)["msy", "harvest"] %*% FLQuant(c(rep(0.1, 60), seq(0.1, 2, length.out = 40)[-40],
    seq(2, 0.7, length.out = 11), rep(0.7, 20)))[, -1]

Create the FLStock

om = as(eq, "FLStock")

Turn it into a stochastic object so that Monte Carlo simulations can be performed

om = propagate(om, 100)

Then perform a stochastic projection

set.seed(234)
srDev = rlnoise(100, fbar(eq) %=% 0, 0.3, 0)

## project for a target fishing mortality
f = fbar(eq)[, -1]

om = ffwd(om, fbar = f, sr = eq, deviances = srDev)

## get rid of burn in period
om = window(om, start = 41)

mets = list(RMSY = function(x, y) rec(x)/refpts(y)["msy", "rec"], YMSY = function(x,
    y) catch(x)/refpts(y)["msy", "yield"], SBMSY = function(x, y) ssb(x)/sbmsy(y),
    FMSY = function(x, y) fbar(x)/fmsy(y))

fqs = FLQuants(lapply(mets, function(m) m(om, eq)))

plot(fqs, iter = 1) + ylim(c(0, NA)) + geom_hline(yintercept = 1, linetype = 2) +
    scale_x_continuous(limits = c(50, 120)) + theme_bw(16) + theme(legend.position = "none")

Figure 2 Operating Model with \(MSY\) reference points.

Seasonal Operating Model

Take an annual FLStock object and “seasonalise” it, i.e. interpolate growth, and split natural mortality by season.

om4 = seasonalise(window(iter(om, 3), end = 81))

If spawning is not year round, as is the case in temperate species, it is necessary to assign the stock-recruit relationship parameters to the spawning season

sr.par = FLPar(NA, dimnames = list(params = dimnames(params(eq))$params, season = 1:4,
    iter = 1))
sr.par[, 2] = params(eq)[, 1]
sr.par[1, 2] = sr.par[1, 2]
sr.par[2, 2] = sr.par[2, 2]
sr4 = predictModel(model = model(eq), params = sr.par)

sr4

An object of class "FLQuants": EMPTY
model:  
rec ~ a * ssb/(b + ssb)
<environment: 0x557cd434f4f0>

params:  
An object of class "FLPar"
      season
params 1     2     3     4    
     a    NA 161.2    NA    NA
     b    NA  28.6    NA    NA
units:  NA

Since the selection pattern is asymptotic, set fbar range to be oldest true age, also make years correspond to current time period.

As there is no spawning outside of \(2^{nd}\) season so set these to 0

range(om4)[c("minfbar", "maxfbar")] = 39

dimnames(om4) <- list(year = an(dimnames(om4)$year) + 1950)

mat(om4)[, , , c(1, 3:4)] = 0

Run a deterministic projection

control = as(FLQuants(f = fbar(om4)[, -1] * 0.5), "fwdControl")
om4 = fwd(om4, control = control, sr = sr4)

plot(window(om4, start = 2010, end = 2030))

Figure 3 Seasonal Operating Model

Simulate a seasonal fishery

f = fbar(om4)[, ac(2021:2030)] * 4
f[, ac(2021:2030), , 1:3] = 0

control = as(FLQuants(f = f), "fwdControl")
om4.2 = window(fwd(om4, control = control, sr = sr4), end = 2030)

plot(window(FLStocks(`All year round` = om4, `Seasonal\nFishery` = om4.2), start = 2010,
    end = 2030))

Figure 4 Seasonal Operating Model and fishery

om = annualise(om4)
om = ffwd(om, sr = eq, fbar = fbar(om)[, -1])

plot(window(FLStocks(Annual = om, Seasonal = om4), start = 2010, end = 2030))

Figure 5 Convert seasonal to an annual Operating Model

More Information

You can submit bug reports, questions or suggestions on FLife at the FLife issue page ¹, or on the FLR mailing list.
Or send a pull request to https://github.com/lauriekell/FLife/
For more information on the FLR Project for Quantitative Fisheries Science in R, visit the FLR webpage ².
The latest version of FLife can always be installed using the devtools package, by calling

library(devtools)
install_github("flr/FLife")

Software Versions

R version 4.2.1 (2022-06-23)
FLCore: 2.6.19.9020
FLPKG:
Compiled: Thu Feb 9 10:52:14 2023
Git Hash: 9511ede

Author information

Laurence KELL. laurie@seaplusplus.co.uk

Acknowledgements

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.