Exploring simulation outputs

This is the third tutorial in a series on conditioning and running simple MixME models. In the first tutorial, we set up a simple single-stock operating model and assembled all the structures needed to run a MixME simulation. In the second tutorial, we set up a simple mixed fishery operating model and ran a simulation with a fixed catch management procedure. This operating model comprising two stocks (cod and haddock) that are each exploited by the same two demersal otter trawl fleets (OTB_A and OTB_B).

In this tutorial, we re-run this simulation and explore the results of this simulation. Specifically, we will first evaluate the quality of the simulation (checking if there are any errors in the simulation). Then, we will explore fleet and stock dynamics, and see if our management objectives have been met. Finally, we will alter the catch target and explore how this impacts the dynamics.

As previously, we first load the necessary libraries, but we now load pre-assembled inputs to MixME. These inputs are identical to the objects we assembled in the previous tutorial; specifically, a mixed fishery operating model, an observation error model for the two stocks, a fixed catch management procedure specifying catches of 1000 tonnes for each stock, and a 20-year simulation.

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

## load example data
data("mixedfishery_MixME_input")

Run the simulation

The loaded inputs are already configured and ready to run, so we can simply plug them into runMixME.

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

Understanding MixME outputs

An overview of the structure of MixME simulation outputs is given here. To quickly recap, om contains the stock and fishery structure that is now updated with the landings, discards and stock numbers and fleet effort in each simulation year. ctrl_obj and args are the management procedure and the simulation arguments taken from the input object. tracking records summary metrics for the modelled stock and fleet dynamics, as well as metrics describing the observed and perceived state of the system by the management procedure, and simulation performance and diagnostics statistics.

Collectively, the outputs from a MixME simulation contains all the necessary information to reconstruct much of the perceived state that drives management advice generation. This means that diagnosing strange simulation dynamics is possible without re-running and disassembling potentially time-consuming and computationally-costly simulations.

Check simulation quality

The first thing we should check is the quality of the simulation. Were there any years when the management procedure failed to generate advice? Were there any years where the forward projection module failed to optimise fleet effort? Are there any over-quota catches that, in a simulation without implementation error, could suggest that effort optimisation was poor?

This is where the tracking object proves to be very useful. We can see that there were no management procedure failures.¹ iterfail is a matrix of years by iterations (or replicates) where advice successes and failures are recorded as '0' and '1', respectively. When an advice failure occurs, a '1' is recorded in each subsequent year of that iteration, regardless of whether or not advice failure persists.

## Check for advice failure
apply(res$tracking$iterfail, 1, mean)

Additionally, there were no effort optimisation failures. optim records the optimised objective function value and the convergence status in each years and iteration. Convergence is '0' in all years, denoting successful convergence. If we encountered convergence problems, we might want to check message for the convergence message returned by the optimisation algorithm. message also contains the 'rescale' metric - which is consistently NA here. This is an optional procedure (disabled by default) where effort is re-scaled down if quota over-shoot exceeds some value to bring the overall over-shoot to zero.

## Check for effort optimisation failure
res$tracking$optim

## Check effort optimisation message
res$tracking$message

Consistent with this, over-quota catches are smaller than 1 kg. overquota records the over-quota catch for each stock by each fleet in every year and iteration. On balance, the simulation has performed well!

## Check over-quota catches
round(res$tracking$overquota,3)["cod",,,1]
round(res$tracking$overquota,3)["had",,,1]

## What is the maximum overshoot?
max(res$tracking$overquota, na.rm = TRUE)

Explore stock and fleet dynamics

Stock and fleet dynamics can be rapidly evaluated by plotting time-series summaries with plot_timeseries_MixME, specifying the specific property of interest to generate figures for stock spawning biomass, overall fishing mortality, total catch and advice uptake, and fleet effort. The figures display data that are mostly aggregated at the total-stock-level because the goal is to provide a quick high-level view of the overall dynamics of the system. These figures use ggplot syntax and can therefore be further customised by users.

It is immediately evident that the low catch target of 1000 tonnes for cod and haddock allows for rapid rebuilding of the two stocks to unprecedented high levels. This coincides with a massive decrease in the level of fishing mortality and catches compared to historical levels.

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

At the same time, we can see that 100 % of cod target catch is taken, whereas the percentage of haddock target catch that is taken drops from 60 % at the beginning of the simulation to a little over 10 % by 2039. This strongly suggests that cod is the effort-limiting 'choke' stock - this is the first quota that is fully consumed by the fleet - for both fleets in the system - and a quick look at the tracking object confirms this. In more complex the mixed fishery systems, the figure of % uptake at the stock-level might be insufficient to identify the 'choke' stock because fleets might have different 'chock' stocks and no stock achieves 100 % catch target uptake.

The characterisation of 'choke' stocks for each fleet is one of the fundamental objectives of mixed fisheries simulations, and MixME makes this information very accessible.

plot_timeseries_MixME(res, quantity = "uptake")

## Check choke stock
res$tracking$choke

The reason haddock % uptake drops over time is that the increasing availability of cod leads to a higher cod catch per unit effort for both fleets. This means that a progressively lower effort is needed to achieve the cod catch target - with negative implications for haddock utilisation. We can visualise this effort time-series, but we need to begin plotting from 2020 because historical effort levels are so much higher!

plot_timeseries_MixME(res, quantity = "effort", minyr = 2020)

It is often useful to calculate a 'risk' metric for each stock - a measure of risk that stock safety or long-term viability is compromised by implementing a given management procedure. MixME defines risk as the proportion of iterations (read replicates) in each year where stock spawning biomass falls below the biological limit reference point ($B_{lim}$). We haven't yet defined $B_{lim}$ for our two stocks but we can do this now.

The figure for risk is not exceptionally informative because we only have one replicate, so the only two possible values are '0' or '1'. However, we can see that $B_{lim}$ is not breached at any point in the simulation.

## Add reference points
hcrpars <- list(cod = c(Blim = 107000), 
                had = c(Blim = 9227))
res$ctrl_obj$phcr <- mseCtrl(args   = list(hcrpars = hcrpars))

## Plot risk time-series
plot_timeseries_MixME(res, quantity = "risk")

It might also useful to breakdown aspects of harvesting by fleet. We can do this using built-in summary functions. We can also extract the properties of interest from the Operating Model or tracking object, but this is easier!

summary_catch_MixME(res, minyr = 2020, byfleet = TRUE)
summary_uptake_MixME(res, minyr = 2020, byfleet = TRUE)

Changing management parameters

Constraining cod and haddock catches to 1000 tonnes each allowed both stocks to massively increase over the course of the simulation but resulted in increasing underutilisation problems for haddock. What if we increase the catch target for cod? Will this alleviate the haddock underutilisation problem?

Updating the management procedure is as simple as changing the catch target in the simulation inputs and rerunning MixME.

input <- mixedfishery_MixME_input
input$ctrl_obj$hcr@args$ctrg$cod <- 10000
input$ctrl_obj$hcr@args$ctrg$had <- 1000

res2 <- runMixME(om  = input$om, 
                oem = input$oem,
                ctrl_obj = input$ctrl_obj,
                args     = input$args)

The simulation runs without issues (great!), and we can see that the choking patterns have changed (interesting!). Now, haddock chokes both fleets initially, but cod becomes the dominant choke stock in the latter part of the simulation.

## Check for advice failure
apply(res2$tracking$iterfail, 1, mean)

## Check for effort optimisation failure
res2$tracking$optim[2,,]

## Check for effort optimisation message
res2$tracking$message

## Check maximum overshoot
max(res2$tracking$overquota, na.rm = TRUE)

## Check choke stock
res2$tracking$choke

## plot stock properties
plot_timeseries_MixME(res2, quantity = "ssb")
plot_timeseries_MixME(res2, quantity = "fbar")
plot_timeseries_MixME(res2, quantity = "catch")
plot_timeseries_MixME(res2, quantity = "uptake")

Perhaps we should be fishing the system even harder!

input$ctrl_obj$hcr@args$ctrg$cod <- 5000000
input$ctrl_obj$hcr@args$ctrg$had <- 1000000

res3 <- runMixME(om  = input$om, 
                 oem = input$oem,
                 ctrl_obj = input$ctrl_obj,
                 args     = input$args)

We've clearly gone too far! We start to see effort optimisation problems as we fish the stocks to extinction.

## Check for advice failure
apply(res3$tracking$iterfail, 1, mean)

## Check for effort optimisation failure
res3$tracking$optim[2,,]

## Check for effort optimisation message
res3$tracking$message

## Check maximum overshoot
max(res3$tracking$overquota, na.rm = TRUE)

## Check choke stock
res3$tracking$choke

## plot stock properties
plot_timeseries_MixME(res3, quantity = "ssb")
plot_timeseries_MixME(res3, quantity = "fbar")
plot_timeseries_MixME(res3, quantity = "catch")
plot_timeseries_MixME(res3, quantity = "uptake")
plot_timeseries_MixME(res3, quantity = "effort", minyr = 2020, maxyr = 2038)

Overall, we can see that it is easy to run, modify and interpret MixME simulations. However, we have been working with a very simple fixed catch management procedure and are not considering sources of uncertainty in the system. In the next set of tutorials, we will look at setting up more complex management procedures.

Notes

¹: This is not surprising. The management procedure is returning a fixed catch target for each stock. There is no stock-estimation or short-term forecast to go awry.