Conditioning a simple single‐stock Operating Model

In this first tutorial, we set up a simple single-stock and single-fleet operating model using some sample data provided with the MixME package. You should be familiar with FLR data objects and methods before starting this tutorial, particularly the

We will cover the basic structure required for inputs to a MixME simulation and introduce some of the conditioning functions to help build the operating model. We will then assemble the remaining objects required to run a MixME simulation. For simplicity, we will implement fixed catch advice over the duration of the simulation, so that stock estimation is not required (in reality we will generate perfect observations of the stock that are ignored by the harvest control rule module... but more on that later).

Let's load the necessary libraries and load a sample dataset to condition our first operating model.

## load libraries
library(FLCore)
library(FLFishery)
library(mse)
library(stockassessment)
library(MixME)

## load example data
data("singlestock_MixME_om")

Basic Operating Model structure

In case anyone is interested, this is haddock (Melanogrammus aeglefinus) and we have stock and catch data for the years 1993 - 2019. These stock and catch data are age-structured and stored in separate FLR objects. The key stock information is stored within slots (numbers n, natural mortality m, stock mean individual weight wt and proportion mature mat) in a named FLBiol object (had), which is in turn nested in a named FLBiols object (stks). The catch information (landings numbers landings.n, landings mean individual weight landings.wt, discards numbers discards.n, discards mean individual weight discards.wt, selectivity catch.sel and catchability catch.q) are stored as a named FLCatch (had) in a named FLFishery (fleet) that is in turn nested in a named FLFisheries object (flts).

## A list of stks (FLBiols) and flts (FLFisheries)
summary(singlestock_MixME_om)

## Dimensions for each stock and fleet
summary(singlestock_MixME_om$stks)
summary(singlestock_MixME_om$flts)

This hierarchical structure is extremely convenient because we can nest multiple stocks and fleets. A few important points:

1. The name of the FLBiols object must be stks and the names of the FLFisheries object must be flts.
2. Stocks and their corresponding fleet catches must share same name. In this case both our stock and catch object are called had.
3. By default, the catchability slot in an FLCatch does not have a year dimension. This is easily modified, and MixME expects that catchability can be indexed by year.
4. A parameterised stock-recruit model must be available from the rec slot of the stock FLBiol.

## Check that stock and catch names match
names(singlestock_MixME_om$flts$fleet) %in% names(singlestock_MixME_om$stks)

## Check catchability dimensions
dim(catch.q(singlestock_MixME_om$flts$fleet$had))

## Check stock-recruit relationship
singlestock_MixME_om$stks$had@rec

The good news is that we have most of the basic structure for our operating model. This basic structure remains the same irrespective of whether we are modelling 1 or 100 stock and fleets.

There are a few steps left before we have an operating model that we can pass to MixME:

1. Estimate historic quota-share for our single fleet
2. Project stock n years into the future
3. Calculate stock numbers in initial year
4. Generate an observation error model

Estimate historic quota-share

In a multi-fleet system, we need some way of partitioning the overall target catch for each stock across the fleets that catch this stock. This proportional share of stock quota (quota-share) is not currently represented in any FLR objects, so we must generate an FLQuant containing the estimated quota-share for each catch and attach it as an attribute of the corresponding FLCatch.

In our single-stock example, quota-share will be 1 in each year (because all available quota is taken by our single fleet). In more complex mixed fishery systems, quota-share can be a tricky parameter to estimate. One convenient (but also potentially imprecise) method is to assume that the proportional share of landings for each stock across fleets corresponds to the quota-share. This is useful because we already have all the information that we need for the calculation stored in the FLFisheries object. MixME has a function that allows you to calculate the quota-share in this way for each fleet catch.

## Calculate quota-share
out <- calculateQuotashare(stks = singlestock_MixME_om$stks, 
                           flts = singlestock_MixME_om$flts, verbose = TRUE)
singlestock_MixME_om$stks <- out$stks
singlestock_MixME_om$flts <- out$flts

## Check 
attr(singlestock_MixME_om$flts$fleet$had, "quotashare")

Project stock properties into future years

MixME requires that Operating Model year dimensions run the full length of the simulation. So if we want to carry out a 20-year simulation, we need to extend the stock and fishery structures from 2019 (the current final data year) to 2039. In each of these projection years (2020 - 2039), we need to estimate the stock and fishery parameters that are not calculated dynamically during the simulation.

Future parameter values are often estimated from historical values. The easiest way to do this is to take the average of the most recent historical data. The stock parameters that we need to project are natural mortality, proportion mature, and stock mean individual weights. The fleet parameters that we need to project for each stock catch are the landings and discards mean individual weights, catch selection, proportion of catch retained (also known as the landed fraction), catchability and quota-share.

MixME has a function that will carry out this parameter projection. You might expect that some parameters covary in time and should therefore be estimated over identical time-scales. And you would be right! The function implicitly defines three groups of parameters:

• All stock parameters, landings and discards mean individual weights and landed fraction (wts.nyears)
• Catchability and catch selection (sel.nyears) [^1]
• Quota-share (qs.nyears)

For simplicity, we take the average parameter value over the most recent 3 years for all parameters.

out <- stfMixME(singlestock_MixME_om,
                method = "yearMeans", 
                nyears = 20, 
                wts.nyears = 3, 
                sel.nyears = 3, 
                qs.nyears = 3, 
                verbose = TRUE)

## Check projection
plot(out$stks$had)
plot(out$flts$fleet)

## Overwrite outputs
singlestock_MixME_om$stks <- out$stks
singlestock_MixME_om$flts <- out$flts

Calculate stock numbers in initial year

The stock numbers slot (n) contains the numbers for each age class at the beginning of the year. MixME needs to know what the stock numbers are at the beginning of the first simulation year to project the stock given some target catch.

The good news is that we have all the data we need to calculate this. We know what the numbers at the beginning of 2019 were, we know the natural mortality and fishing mortality rates in 2019, and we have a parameterised stock-recruit model to estimate recruitment. The easiest way to calculate the 2020 stock numbers is via a 1-year short-term forecast using the FLasher package. We don't actually care about the forecast target here because the 2020 numbers are calculated from existing numbers (rather than a defined target).

## Initial simulation year
iy <- 2020

## Create an arbitrary forecast target for each fishery
ctrlArgs <- lapply(1:length(singlestock_MixME_om$flts), function(x) {
  list(year = iy,
       quant = "effort",
       fishery = names(singlestock_MixME_om$flts)[x],
       value = 1)
})
## make an FCB matricx
ctrlArgs$FCB <- makeFCB(biols = singlestock_MixME_om$stks, flts = singlestock_MixME_om$flts)

## Generate effort-based FLasher::fwd forecast control
flasher_ctrl <- do.call(FLasher::fwdControl, ctrlArgs)

## carry out the projection
omfwd <- FLasher::fwd(object    = singlestock_MixME_om$stks, 
                      fishery   = singlestock_MixME_om$flts, 
                      control   = flasher_ctrl)

## assign the 
singlestock_MixME_om$stks$had@n[,ac(iy)] <- omfwd$biols$had@n[,ac(iy)]

Generate an Observation Error Model

The Observation Error Model defines how future stock and catch observations are generated. As we discuss here, these observations are critical inputs to the estimation of stock status and management advice generation. Building an Observation Error Model can get quite complex once we start considering observation error and include catch and survey indices. We will cover some of this complexity in a future tutorial, but I don't want to get too far into the weeds here.

For the moment we are constructing a very simple simulation using a constant catch rule and all we need for the moment is some structure that corresponds to the observed stock. [^2]

## convert FLBiol to FLStocks
stk_oem <- FLStocks(lapply(singlestock_MixME_om$stks, function(x) {
  xx <- as.FLStock(x, singlestock_MixME_om$flts$fleet)
  stock.n(xx)[] <- NA
  stock(xx)[]   <- NA
  harvest(xx)[] <- NA
  return(xx)
  }))

Assemble simulation inputs

Great! We have a complete Operating Model (and part of an Observation Error Model). What remains is to define the Management Procedure and simulation arguments.

Build the MixME input object

MixME has a function that, given a conditioned Operating Model, will generate a complete set of generic inputs to MixME that can then be customised to suit any simulation requirements.

I've specified two additional arguments. The 'management lag', the delay between specification and implementation of management targets, is zero, amd 'management type' is fixed catch.

singlestock_MixME_input <- makeMixME(om = singlestock_MixME_om, 
                                     catch_obs = stk_oem, 
                                     management_lag = 0, 
                                     management_type = "fixedC", 
                                     parallel = FALSE)

We need to update some of the simulation and management settings, but it is also an opportunity to look at some of important inputs that the function has created. Let's start with the management procedure methods and arguments - these are supplied as a list within each module with the overall structure:

• est (mseCtrl)
  • args (list)
    • estmethod (list)
      • stock name (character or a user-supplied function)
    • other arguments required by estmethod (list)
• phcr (mseCtrl)
  • args (list)
    • hcrpars (list)
• hcr (mseCtrl)
  • args (list)
    • hcrmethod (list)
      • stock name (character or a user-supplied function)
    • other arguments required by hcrmethod (list)
• isys (mseCtrl)
  • args (list)
    • isysmethod (list)
      • stock name (character or a user-supplied function)
    • other arguments required by isysmethod (list)
• fwd (mseCtrl)
  • args (list)
    • sr_residuals ()
    • proc_res ()

You may be thinking that forward projection module is not a part of the management procedure. You're absolutely correct! This is simply a convenient way to pass parameter noise to the operating model.

Customising management inputs

The snippet of code below shows the default arguments for the management module (ctrl_obj).

• Stock estimation has been set to perfect observation - this is fine for the basic simulation that we are constructing here.
• The harvest control rule parameterisation module is empty - again, this is fine for our basic simulation
• The harvest control rule module specifies 1 tonne of catch for each stock - we will want to update this module.
• Advice implementation system is empty - this is fine for our simulation
• The forward projection module is empty

A strength of MixME is that not every module is needed for each simulation. In this fixed catch management procedure, we don't need to parameterise a harvest control rule, catch targets don't typically need additional advice implementation (this becomes necessary when converting fishing mortality targets to catch targers, or if there is systematic deviation from the advised catch target).

## Stock estimation
singlestock_MixME_input$ctrl_obj$est@args

## Harvest control rule parameterisation
singlestock_MixME_input$ctrl_obj$phcr

## Harvest control rule (update catch target)
singlestock_MixME_input$ctrl_obj$hcr@args

## Advice implementation
singlestock_MixME_input$ctrl_obj$isys@args

## Forward projection
singlestock_MixME_input$ctrl_obj$fwd@args

The catch target is unrealistically low. Let's increase this to 1000 tonnes (assuming that individual mean weights are in kilograms).

singlestock_MixME_input$ctrl_obj$hcr@args$ctrg$had <- 1000

Customising simulation inputs

The last thing that we need to do is check that our simulation arguments are sensible. An definition of each argument is given here.

## Full set of simulation arguments
singlestock_MixME_input$args

These inputs are mostly fine. We just want to update the fishing mortality range used to calculate mean fishing mortality $\bar{F}$ (this is needed for some of the summary metrics).

## Update fbar ranges
singlestock_MixME_input$args$frange$had <- c("minfbar" = 3, "maxfbar" = 5)

Customising Observation Error Model inputs

It's also worth having a look at the arguments being supplied to the observation error model. The definitions of these arguments is given here.

singlestock_MixME_input$oem@args

It's clear that there is one argument that needs updating - if the management lag is 0, then the availability of observed catch should be instantaneous.

## Update observation arguments
singlestock_MixME_input$oem@args$catch_timing$had <- 0

Run the simulation

This dataset is now conditioned and ready to run.

res <- runMixME(om  = singlestock_MixME_input$om, 
                oem = singlestock_MixME_input$oem,
                ctrl_obj = singlestock_MixME_input$ctrl_obj,
                args     = singlestock_MixME_input$args)

We'll leave a detailed analysis the simulation outputs for a future tutorial. For now, let's just quickly check what the simulated dynamics look like.

plot_timeseries_MixME(res, quantity = "ssb")
plot_timeseries_MixME(res, quantity = "catch")
plot_timeseries_MixME(res, quantity = "uptake")
plot_timeseries_MixME(res, quantity = "fbar")

[^1]: Fleet-level variables (capacity and effort) are also projected as part of this group but we don't care about the values. Effort is overwritten and the current version of MixME largely ignores capacity. We just need to have these slots filled for some of the functions to work properly [^2]: Actually, a constant catch rule doesn't even use the generated stock observations... because it's constant catch.