In tcarruth/MSEtool: Management Strategy Evaluation Toolkit

library(MSEtool)

knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
knitr::opts_chunk$set(dpi=85)
options(width = 650)

Introduction

MSEtool is an extension of the DLMtool package for implementing closed-loop simulation testing of management procedures, a process often referred to as management strategy evaluation (MSE). The components of the operating model for conducting MSEs for data-rich stocks is wholly contained within DLMtool, but MSEtool provides the tools for implementing additional features to facilitate model-based management procedures (MPs). Whereas DLMtool is primarily designed for testing data-limited management procedures, MSEtool provides an object-oriented framework for creating data-rich MPs which typically combine an assessment model with a harvest control rule to provide a catch recommendation. Management procedures created by MSEtool are designed to be used with the operating model of DLMtool through the runMSE function.

The user guide to the DLMtool package is accessible via:

DLMtool::userguide()

The overall structure of a closed-loop simulation for a MSE in DLMtool can be described through the classes deployed in the package. An operating model (class OM object) parameterizes the stock (biology), fleet, observation, and implementation dynamics of the MSE. To begin the MSE, the OM object generates the data to be used for a management procedure and stores them in an object of class Data. A function of class MP uses the Data object and returns an object of class Rec containing the management recommendations, e.g., some combination of spatial closures and size, effort and catch limits. From the management recommendations, the OM object is updated and the process is successively repeated as the operating model progresses over time. The output for the closed-loop simulation is returned in an object of class MSE, from which summary statistics can be calculated.

Figure 1: Design of DLMtool package for management strategy evaluations

Overview of MSEtool

To facilitate a standardized design for data-rich MPs, MSEtool creates new classes to build out the MP class. MSEtool uses functions of class Assess and HCR for an assessment model and harvest control rule, respectively. An MP class function can be made from these functions to be tested in a MSE within the DLMtool framework.

Additionally, Assess functions can be used in an assessment-only context (external to management strategy evaluations). Output of an Assess function can be used to provide standardized stock assessment reporting and diagnostics, such as Kobe plots, likelihood profiles for parameter estimates, and retrospective analyses of time series estimates of biomass, fishing mortality, recruitment, etc.

Figure 2: Design of MSEtool package for building management procedures from assessment models and harvest control rules for management strategy evaluations with the DLMtool package. Assess and HCR functions and the Assessment object are unique to MSEtool, but work within the framework of DLMtool objects. Separate functions are also available for standardized reporting of output from assessment models in a stock assessment setting.

Additional features

MSEtool provides functions for converting the results of assessment models into DLMtool OM and Data objects. The SS2OM and SS2Data functions convert SS (Stock Synthesis) assessments, while iSCAM2OM and iSCAM2Data functions are provided for iSCAM assessments.

To support multiple-area operating models, the simmov function estimates the movement matrix based on the prescribed number of areas, biomass distribution among areas in unfished conditions, and the probability of staying in each area in the next time step.

Assessment models

Assessment models in MSEtool are functions of class Assess. Currently, the following assessment functions are available:

DD_TMB (Delay-Difference model and its state-space counterpart, DD_SS)
cDD (Continuous Delay-Differential model and its state-space counterpart, cDD_SS)
SP (Surplus Production model and its state-space counterpart, SP_SS)
SCA (Statistical Catch at age model and variants SCA2 and SCA_Pope)
SSS (Simple Stock Synthesis implementation with fixed terminal depleton assumption)
VPA (Virtual Population Analysis)

These functions take a Data object as an argument and return an object of class Assessment, for example,

DD_result <- DD_TMB(Data = Red_snapper, ...)

where ... are model configuration arguments.

The Assessment object contains the raw output from TMB estimation model, optimization output from nlminb, and parameter estimates with covariance matrix. If applicable, time series estimates of total biomass, vulnerable biomass, spawning stock biomass, fishing mortality and/or exploitation rate, as well as their corresponding MSY and virgin reference points. The Assessment object also reports likelihood values and predicted values for the data used in the assessment model. The object documentation for the Assessment object can be obtained through:

class?Assessment

Descriptions of the Assess functions are provided in a separate vignettes:

browseVignettes("MSEtool")

Reporting and diagnostic functions

Assessment results and summary, along with diagnostic procedures are available for each assessment model. Generic functions summary and plot are available for Assessment objects. The summary function returns a list of estimated current stock status ($F/F_{MSY}$ , $B/B_{MSY}$ in the terminal year), input parameters to the model, derived quantities, and parameter estimates. The plot function generates a Markdown report with time series figures of data, fitted data, and estimated quantities (e.g., biomass and recruitment).

summary(DD_result)
plot(DD_result) # By default, also saves figures in a temporary directory for viewing. The directory can be changed to a user's filespace of choice

Two diagnostic functions are available: retrospective for retrospective analyses (re-running the model by sequentially removing terminal years of data) and profile for profiling the likelihood function over a grid of parameters.

retrospective(DD_result, nyr = 5) # Retrospective analysis going back 5 years from current year
profile(DD_result, R0 = seq(0.75, 1.25, 0.025), h = seq(0.95, 1, 2.5e-3)) # Joint profile over grid of R0 and steepness

Figure 3: Joint profile of R0 and steepness from Delay-Difference assessment model (DD_TMB function). Contours indicate change in likelihood values from minimum (red point).

Harvest control rules (HCRs)

Harvest control rules in MSEtool are functions of class HCR, and are designed to provide different for a given assessment model. HCR functions take an object of class Assessment and return an object of class Rec. Currently, three HCR functions are available as output controls (catch limits): HCR_MSY, HCR40_10, and HCR60_20. HCR_MSY prescribes the TAC to be the catch predicted with fishing at the estimated $F_{MSY}$, i.e., the product of $F_{MSY}$ and estimated vulnerable biomass, while HCR40_10 and HCR60_20 are ramped control rules which reduce the $F_{MSY}$ catch limit based on spawning depletion $SSB/SSB_0$.

par(mar = c(5,5,3,1))
Brel <- seq(0, 1, length.out = 200)
plot(Brel, MSEtool::HCRlin(Brel, 0.1, 0.4), xlab = expression(Estimated~~SSB/SSB[0]), ylab = "TAC adjustment factor \n(proportion of FMSY catch)", main = "40-10 harvest control rule", type = "l", col = "blue")
abline(v = c(0.1, 0.4), col = "red", lty = 2)

Figure 4: The 40-10 harvest control rule prescribes an additional reduction in the TAC based on estimated spawning depletion.

Model-based management procedures

Creating MPs from assessment models and HCRs

MSEtool comes with a set of pre-built data-rich MPs:

avail("MP", all_avail = FALSE)

These MPs use an assessment model and a harvest control rule to come up with a management recommendation (TAC).

Additionally, the make_MP function is a convenient way to stitch together additional MPs from an assessment model and HCR. This is handy if an assessment model with a different configuration than the default is needed (for example, specify dome vulnerability instead of logistic vulnerability in SCA) or a different harvest control rule is desired. The resulting function can then be passed to DLMtool::runMSE() as a management procedure.

DD_MSY <- make_MP(DD_TMB, HCR_MSY, ...)
DD_4010 <- make_MP(DD_TMB, HCR40_10)
myMSE <- DLMtool::runMSE(OM = DLMtool::testOM, MPs = c("FMSYref", "DD_MSY", "DD_4010"))
Tplot2(myMSE)

Figure 5. Trade-off plot of the MSE for the MPs that use the Delay-Difference model with two harvest control rules and the FMSY reference MP.

Model configuration arguments can be passed through ... in make_MP when the MP is made. For example, the Ricker stock-recruitment relationship can be used in the Delay-Difference model instead of the default Beverton-Holt relationship and steepness can be fixed instead of esimated by passing argument SR = "Ricker" and fix_h = TRUE, respectively:

DD_MSY_Ricker_fixsteep <- make_MP(DD_TMB, HCR_MSY, SR = "Ricker", fix_h = TRUE)

This MP will call DD_TMB(..., SR = "Ricker", fix_h = TRUE).

Diagnostic of assessment models in MSE

Three functions are designed to evaluate the performance of Assess models in MSE: prelim_AM, diagnostic_AM, and retrospective_AM (where AM stands for "assessment in MSE").

prelim_AM is designed to evaluate the configuration of the assessment model before running the closed-loop simulation. Given an operating model, this function generates data from the historical period of the MSE, applies the assessment model to those data, then returns the convergence rate of the model. A list of Assessment objects are returned and can be used for troubleshooting. Additional arguments are passed to the Assess function via ...:

prelim_AM(DLMtool::testOM, DD_TMB, ...)
#> Generating Hist object from OM object via runMSE:
#> Loading operating model
#> Optimizing for user-specified movement
#> Optimizing for user-specified depletion
#> Calculating historical stock and fishing dynamics
#> Calculating MSY reference points
#> Calculating B-low reference points
#> Calculating reference yield - best fixed F strategy
#> Returning historical simulations
#> Running DD_TMB with 48 simulations for DLMtool::testOM.
#> Assessments complete.
#> Total time to run 48 assessments: 1.2 seconds
#> 2 of 48 simulations (4.2%) did not converge.
#> See simulation number: 11 21

By default, assessment model results are not saved during the MSE. However, this feature can be turned on in two steps, as shown in the following example:

DD_MSY <- make_MP(DD_TMB, HCR_MSY, diagnostic = "full")
myMSE <- runMSE(..., MPs = "DD_MSY", PPD = TRUE)

First, the MP should be made with make_MP by setting diagnostic = "min" or diagnostic = "full". Second, the call to runMSE should be accompanied with argument PPD = TRUE.

With this combination, two diagnostic objects can be saved and reported in the Misc slot of the MSE object returned by runMSE: convergence diagnostics and the Assessment objects from each assessment call. The first can be obtained with diagnostic = "min", while both are returned with diagnostic = "full".

The diagnosic_AM function plots the convergence information, and retrospective_AM plots the assessment model's retrospective patterns (for a single simulation in the MSE).

diagnostic_AM(myMSE)
retrospective_AM(myMSE, sim = 1, MP = "DD_MSY")

Figure 6. Convergence diagnostics (return by TMB objects and nlminb function) for the DD_MSY MP during the MSE. Here, the assessment model was applied beginning in year 50 of the MSE and repeated every 4th year. In MSEtool, convergence is defined by a positive-definite Hessian matrix (top left figure). In the case of non-convergence, the previous management recommendation is used.

Figure 7. Retrospective analysis of the Delay-Difference model in the DD_MSY MP during the MSE for simulation #1. Operating model (true) values of spawning biomass (SSB), relative spawning biomass (SSB_SSBMSY), spawning biomass depletion (SSB_SSB0), fishing mortality (F), relative F (F_FMSY), and vulnerable bimoass (VB) are in dark black. Colored lines indicate model estimates sequentially over the projection years. Note that this assessment model assumes the vulnerable biomass to be the spawning biomass.

Once it is determined that the MSEtool assessment models are performing adequately (no major convergence issues), then the performance of the MPs that use these models can be evaluated within the DLMtool package.

Obtaining management recommendations

Like all other output MPs, DLMtool::TAC() can be used to obtain a catch limit:

DD_MSY <- make_MP(DD_TMB, HCR_MSY)
DD_4010 <- make_MP(DD_TMB, HCR40_10)
snapper_TAC <- TAC(Red_snapper, MPs = c("DD_MSY", "DD_4010"))
plot(snapper_TAC)

Figure 8. Stochastic TAC recommendations from the Delay-Difference model with two harvest control rules. When the 40-10 harvest control rule is applied, the prescribed TAC is close to zero because the assessment model estimated the terminal year biomass depletion to be 0.13 (close to the limit reference point).

Acknowledgements

MSEtool development is funded by the Canadian Department of Fisheries and Oceans and benefits from ongoing collaborations with the Natural Resources Defense Council and number of Canadian government scientists including Robyn Forrest, Sean Anderson and Daniel Duplisea.

The DLMtool package was developed at the University of British Columbia in collaboration with the Natural Resources Defense Council. DLMtool development has been funded by the Gordon and Betty Moore Foundation, the Packard Foundation, U.S. National Oceanic and Atmospheric Administration, Fisheries and Oceans Canada, the Walton Foundation, Resources Legacy Fund, the Natural Resources Defense Council, the United Nations Food and Agricultural Organization and the California Department of Fish and Wildlife.