mahiMSE User Manual
An R Package for Informing Management of U.S. South Atlantic Dolphinfish Fisheries
Tom Carruthers (tom@bluematterscience.com)
2026-06-30
1 Foreword
This research and decision support tools use operating models and MSE (Management Strategy Evaluation) simulation testing to characterize scientific uncertainty to inform robust management decision making.
2 Objective of this document
This document is intended to guide prospective users on how to:
- work with mahi operating models
- design and test management procedures (MPs)
- summarize MP performance
3 Version Notes
The package is subject to ongoing testing and development. If you find a bug or a problem please send a report to tom@bluematterscience.com so that it can be fixed!
4 Background information
Project information, links and results are available from a project webpage
5 Software
The openMSE R package was used to organize fishery data, define operating models, condition operating models to data, calculate reference points and specify and test existing and alternative management options.
The Rapid Conditioning Model was used to fit operating models to observed data.
6 In Development
A number of additional features for the mahiMSE package are in development:
| Limitation / new feature / bug | Description |
|---|---|
| TAC by fleet and area | As of 0.8.2 this is only by fleet |
| Add spatial VAST index as an option | Needs adding to the addFind function for OM building |
| Add CPUE index | Something about blue-blobs etc |
| Current VB indices are simulated | Until real fish-dependent indices can be added to the OM objects, these are simulated included the historical period |
| Current SSB indices are simulated | As above, no SSB index is actually available |
| Carry out peer-review of MSE framework | Are the current OMs, MPs and performance metrics sufficient to meet ToRs? Can the process now progress to formal MP testing and potentially adoption? |
| Alternative Indices of Abundance for MPs | Formalization of PR_tourn, MRIP, ENSO, Blue-Blob |
| Effort dynamics model | To model opportunity a model is needed to predict whether ACLs will be reached given available fishing characteristics. |
| Slick improvements | Add landings by fleet to time series plots |
| Slick improvements | Make robustness slick or add to existing |
7 Installation
Install R for Windows
Install RStudio
Install latest openMSE packages
Open RStudio and enter the following commands into the console:
install.packages('openMSE')
install.package('pak')You will need to upgrade to the very latest version of MSEtool.
MSEtool was radically overhauled to consolidate all the various features required from applications and workshops for more than 200 fisheries worldwide. MSEtool version 4.x.x is faster, smaller, more flexible and much more intuitive to use. Furthermore it includes features for seasonality that are particularly applicable for short-lived species such as mahi.
As of March 2026 v4.x.x exists on the noaa-mahi branch of the MSEtool R package repository. By the end of 2026 this will be rolled into the regular CRAN release of MSEtool (v4.0.0) and openMSE (v2.0).
To install MSEtool v4.x.x it you use the following code:
library('pak')
pkg_install('blue-matter/MSEtool@noaa-mahi') # you need the prerelease branch- Install the latest version of the mahiMSE GitHub package:
pkg_install("blue-matter/mahiMSE") # mahiMP, plotting etc- Check that the installation is successful by finding this help file (the mahiMSE vignette):
library(mahiMSE)
vignette('mahiMSE') # Not working with pak (its the document you are reading!)- Install the Reference Set and Robustness Set OM Packages
pkg_install("blue-matter/mahiRefSet") # 32 OMs
pkg_install("blue-matter/mahiRobSet") # 13 OMs
If you have any difficulties please send an inquiry including some reproducible code to: tom@bluematterscience.com
8 Quick start: Run an MSE with a Small OM
Load libraries:
library(mahiMSE)Build a historical simulation for the small (8 simulation) version of OM_29 (the reference case operating model):
Hist = Simulate(smallOM)Create two management procedures using the mahiMP function:
StatQuo = HalfEff = mahiMP # copy the mahiMP
formals(HalfEff)$rel_TAE = 1/2 # all fleets effort is halved
class(HalfEff) = class(StatQuo) = "mp" # assign class 'mp'
myMSE = Project(Hist, c('StatQuo','HalfEff')) # Project Plot projected SSB and Landings:
mahiplot(mahiResults(myMSE)) # plot annual results
#> Processing MSE Results
#> Dimensions: nsim = 8, nt = 148, np = 80, ny = 228, na = 17, nf = 8, nl = 39, nr = 5, nMSE = 1, nMP = 2
#> Extracting Spawning Stock Biomass (SSB)
#> Extracting Biomass (B)
#> Extracting Landings (Landings)
#> Calculating Catch Diff (Cdif)
#> Calculating Proportion of Landings in each age class Biomass (PropCAA)
#> Calculating Vulnerable Biomass as a retained catch rate predictor (VBr)
#> Calculating Vulnerable Biomass as a catch rate predictor (VB)
#> Calculating MSY reference points (Brel, Frel)
#> Returning mahiResults object
mahiplot_fleet(mahiResults(myMSE))
#> Processing MSE Results
#> Dimensions: nsim = 8, nt = 148, np = 80, ny = 228, na = 17, nf = 8, nl = 39, nr = 5, nMSE = 1, nMP = 2
#> Extracting Spawning Stock Biomass (SSB)
#> Extracting Biomass (B)
#> Extracting Landings (Landings)
#> Calculating Catch Diff (Cdif)
#> Calculating Proportion of Landings in each age class Biomass (PropCAA)
#> Calculating Vulnerable Biomass as a retained catch rate predictor (VBr)
#> Calculating Vulnerable Biomass as a catch rate predictor (VB)
#> Calculating MSY reference points (Brel, Frel)
#> Returning mahiResults object
mahiplot_MSY(mahiResults(myMSE))
#> Processing MSE Results
#> Dimensions: nsim = 8, nt = 148, np = 80, ny = 228, na = 17, nf = 8, nl = 39, nr = 5, nMSE = 1, nMP = 2
#> Extracting Spawning Stock Biomass (SSB)
#> Extracting Biomass (B)
#> Extracting Landings (Landings)
#> Calculating Catch Diff (Cdif)
#> Calculating Proportion of Landings in each age class Biomass (PropCAA)
#> Calculating Vulnerable Biomass as a retained catch rate predictor (VBr)
#> Calculating Vulnerable Biomass as a catch rate predictor (VB)
#> Calculating MSY reference points (Brel, Frel)
#> Returning mahiResults object
9 Getting Help
9.1 Splash Page & Trial Specifications Document (TSD)
The project splash page provides a centralized place for project status, links to supporting documentation and other help on MSE concepts. If you had to have one mahi MSE link on your desktop, this is it!
The TSD is intended to provide sufficient detail on the technical aspects of the MSE to ensure reproducibility.
9.2 Inline R Package Help
User-facing mahiMSE functions have inline R help. Along with user-facing functions of the OpenMSE package, these can be accessed with the ? operator:
?mahiMP
?mahiplot
?Simulate
?Project
All of the code for the mahiMSE R package is open source. Internal functions such as those for computing things like the impact of trip-limits on retention rate, are all available to the user.
Use the double colon operators to locate and view code:
mahiMSE::do_HCR # Internal function for calculates advice using a harvest control rule
mahiMSE::gen_ind # Internal function for generating indices of abundance for fleets
mahiMSE::do_TL # Internal function for predicting catch rate distribution and retention rate for a trip limit
9.3 Useful Package Objects
It can be challenging to keep track of the nomenclature of things like the various fleets and areas. The mahiMSE package includes a number of objects that are automatically loaded into the session to help the user.
Areas # Named vector of area names
Areas2 # Alternative area names
Fleets # Fleet names
Example_Data # An example data file recorded after one projected time step, MSEtool object of class 'data'
Example_Data2 # An example data file recorded after many projected time steps given over-exploitation, MSEtool object of class 'data'
Example_Data3 # An example data file recorded after many projected time steps given under-exploitation, MSEtool object of class 'data'
rTAC # regional, seasonal catch from last historical year for reference a 3D array [nseason, nfleet, narea]
sVAST # The spatial VAST index (areas sum to 1) by time step (quarterly), a matrix [ntimestep, narea]
oVAST # The original VAST output (data table with CVs, strata etc)
TL_info # Trip limit information. The distribution of catch rates at 2018 - 2022 stock levels. This is used
# to predict the catch rate and hence release rates within the mahiMP
ref_grid # The reference grid of operating models
cond_grid # The grid of conditioned (fitted to data) reference grid operating models
C_yfr # Catches (kg) in a 3D array [historical time period, fleet, area]
Length_classes # The length classes of the operating model (helps to specify length vectors for MP inputs)
Time_steps # All of the time steps of the historical and projection period
hTime_steps # The time steps of the historical period
pTime_steps # The time steps of the projection period
10 Mahi MSE Dimensions: Areas, Fleets, Time Steps and Length Classes
When exploring management procedures you may wish to specify advice for specific fleets and/or areas, and also control how management changes occur over time. In this section we review the dimensions of the MSE and show you how you can access these.
10.1 Areas
The operating model has five explicit areas which can be used to specify management advice:
Figure 1. The five spatial regions of the operating model.
Table 1. The discrete spatial strata.
| Code | Alternative code | Description |
|---|---|---|
| CAR+FLK | CAR+SFL | Caribbean, South Florida and the Keys |
| NCA | SAR | Off Caribbean, South Florida and the Keys |
| NCFL | SE | Wilmington down to Mid Florida |
| NED | NED | Off North coast |
| NNC+VBM | NC+NE | Wilmington to Montauk |
10.2 Fleets
Table 2. Individual fleets used to partition data for operating model conditioning and projection of fishing dynamics.
| Code | Fleet | Areas | Years | Selectivity | Description |
|---|---|---|---|---|---|
| USCom | USA commercial longline | All | All | Logistic | US pelagic longline fleet |
| RecN | Recreational North | NNC+VBM | All | Domed | Private recreational fishing |
| RecS | Recreational South | CAR+FLK, NCFL | All | Domed | Private recreational fishing |
| HireN | For hire North | NNC+VBM | All | Domed | Headboat / party boat / charter |
| HireS | For hire South | CAR+FLK, NCFL | All | Domed | Headboat / party boat / charter |
| Intl | International | CAR+FLK, NCA, NED | All | Logistic (USCom) | Assumed to have the same seasonal distribution as the USCom fleet. |
| Disc | Discard | CAR+FLK, NCA, NNC+VBM | All | Logistic (USCom) | Assumed to have same catches as UnRep fleet and the same seasonal distribution as the USCom fleet. |
| UnRep | Unreported | CAR+FLK | All | Logistic (USCom) | Assummed to have the same catches as the Disc fleet and the same seasonal distribution as the USCom fleet. |
Here is the breakdown of catches (t) in 2022 for Mahi:
nts = dim(C_yfr)[1]
last4seasons = nts-(3:0)
Cagg = round(apply(C_yfr[last4seasons,,],2:3,sum)/1E3,1) # sum over fleet and area
kable(Cagg)| CAR+FLK | NCA | NCFL | NED | NNC+VBM | |
|---|---|---|---|---|---|
| USCom | 16.5 | 0 | 69.6 | 0.0 | 18.1 |
| RecN | 0.0 | 0 | 0.0 | 0.0 | 643.5 |
| RecS | 1757.5 | 0 | 1245.2 | 0.0 | 0.0 |
| HireN | 0.0 | 0 | 0.0 | 0.0 | 202.0 |
| HireS | 286.8 | 0 | 42.1 | 0.0 | 0.0 |
| Intl | 5019.5 | 287 | 0.0 | 26.3 | 0.0 |
| Disc | 292.8 | 0 | 193.1 | 0.0 | 117.3 |
| UnRep | 1244.3 | 0 | 0.0 | 0.0 | 0.0 |
10.3 Time Steps
The mahi operating models run on a seasonal time step (1: Winter, 2: Spring, 3: Summer, 4:Autumn). You can view historical, projected and all timesteps using the data objects loaded with mahiMSE
Time_steps # all
#> [1] 1986.00 1986.25 1986.50 1986.75 1987.00 1987.25 1987.50 1987.75 1988.00 1988.25 1988.50 1988.75 1989.00 1989.25 1989.50 1989.75 1990.00 1990.25 1990.50 1990.75 1991.00 1991.25 1991.50 1991.75 1992.00 1992.25 1992.50 1992.75 1993.00 1993.25 1993.50 1993.75 1994.00 1994.25 1994.50 1994.75 1995.00 1995.25 1995.50 1995.75 1996.00 1996.25 1996.50 1996.75 1997.00 1997.25 1997.50 1997.75 1998.00 1998.25 1998.50 1998.75 1999.00 1999.25 1999.50 1999.75 2000.00 2000.25 2000.50 2000.75 2001.00 2001.25 2001.50 2001.75 2002.00 2002.25 2002.50 2002.75 2003.00 2003.25 2003.50 2003.75 2004.00 2004.25 2004.50 2004.75 2005.00 2005.25 2005.50 2005.75
#> [81] 2006.00 2006.25 2006.50 2006.75 2007.00 2007.25 2007.50 2007.75 2008.00 2008.25 2008.50 2008.75 2009.00 2009.25 2009.50 2009.75 2010.00 2010.25 2010.50 2010.75 2011.00 2011.25 2011.50 2011.75 2012.00 2012.25 2012.50 2012.75 2013.00 2013.25 2013.50 2013.75 2014.00 2014.25 2014.50 2014.75 2015.00 2015.25 2015.50 2015.75 2016.00 2016.25 2016.50 2016.75 2017.00 2017.25 2017.50 2017.75 2018.00 2018.25 2018.50 2018.75 2019.00 2019.25 2019.50 2019.75 2020.00 2020.25 2020.50 2020.75 2021.00 2021.25 2021.50 2021.75 2022.00 2022.25 2022.50 2022.75 2023.00 2023.25 2023.50 2023.75 2024.00 2024.25 2024.50 2024.75 2025.00 2025.25 2025.50 2025.75
#> [161] 2026.00 2026.25 2026.50 2026.75 2027.00 2027.25 2027.50 2027.75 2028.00 2028.25 2028.50 2028.75 2029.00 2029.25 2029.50 2029.75 2030.00 2030.25 2030.50 2030.75 2031.00 2031.25 2031.50 2031.75 2032.00 2032.25 2032.50 2032.75 2033.00 2033.25 2033.50 2033.75 2034.00 2034.25 2034.50 2034.75 2035.00 2035.25 2035.50 2035.75 2036.00 2036.25 2036.50 2036.75 2037.00 2037.25 2037.50 2037.75 2038.00 2038.25 2038.50 2038.75 2039.00 2039.25 2039.50 2039.75 2040.00 2040.25 2040.50 2040.75 2041.00 2041.25 2041.50 2041.75 2042.00 2042.25 2042.50 2042.75
hTime_steps # historical
#> [1] 1986.00 1986.25 1986.50 1986.75 1987.00 1987.25 1987.50 1987.75 1988.00 1988.25 1988.50 1988.75 1989.00 1989.25 1989.50 1989.75 1990.00 1990.25 1990.50 1990.75 1991.00 1991.25 1991.50 1991.75 1992.00 1992.25 1992.50 1992.75 1993.00 1993.25 1993.50 1993.75 1994.00 1994.25 1994.50 1994.75 1995.00 1995.25 1995.50 1995.75 1996.00 1996.25 1996.50 1996.75 1997.00 1997.25 1997.50 1997.75 1998.00 1998.25 1998.50 1998.75 1999.00 1999.25 1999.50 1999.75 2000.00 2000.25 2000.50 2000.75 2001.00 2001.25 2001.50 2001.75 2002.00 2002.25 2002.50 2002.75 2003.00 2003.25 2003.50 2003.75 2004.00 2004.25 2004.50 2004.75 2005.00 2005.25 2005.50 2005.75
#> [81] 2006.00 2006.25 2006.50 2006.75 2007.00 2007.25 2007.50 2007.75 2008.00 2008.25 2008.50 2008.75 2009.00 2009.25 2009.50 2009.75 2010.00 2010.25 2010.50 2010.75 2011.00 2011.25 2011.50 2011.75 2012.00 2012.25 2012.50 2012.75 2013.00 2013.25 2013.50 2013.75 2014.00 2014.25 2014.50 2014.75 2015.00 2015.25 2015.50 2015.75 2016.00 2016.25 2016.50 2016.75 2017.00 2017.25 2017.50 2017.75 2018.00 2018.25 2018.50 2018.75 2019.00 2019.25 2019.50 2019.75 2020.00 2020.25 2020.50 2020.75 2021.00 2021.25 2021.50 2021.75 2022.00 2022.25 2022.50 2022.75
pTime_steps # projected
#> [1] 2023.00 2023.25 2023.50 2023.75 2024.00 2024.25 2024.50 2024.75 2025.00 2025.25 2025.50 2025.75 2026.00 2026.25 2026.50 2026.75 2027.00 2027.25 2027.50 2027.75 2028.00 2028.25 2028.50 2028.75 2029.00 2029.25 2029.50 2029.75 2030.00 2030.25 2030.50 2030.75 2031.00 2031.25 2031.50 2031.75 2032.00 2032.25 2032.50 2032.75 2033.00 2033.25 2033.50 2033.75 2034.00 2034.25 2034.50 2034.75 2035.00 2035.25 2035.50 2035.75 2036.00 2036.25 2036.50 2036.75 2037.00 2037.25 2037.50 2037.75 2038.00 2038.25 2038.50 2038.75 2039.00 2039.25 2039.50 2039.75 2040.00 2040.25 2040.50 2040.75 2041.00 2041.25 2041.50 2041.75 2042.00 2042.25 2042.50 2042.75If for example, you wanted to view these four quarters of catches for both 2021 and 2022 that would be eight timesteps.
Lets sum those up across areas and tabulate them:
library(kableExtra)
dim(C_yfr) # catches in kg in a 3D array [timestep, fleet, area]
#> [1] 148 8 5
nts = dim(C_yfr)[1] # number of time steps in total 37 years, 4 seasons is 148 seasons total
Csub = C_yfr[nts-(7:0),,] # catches for last 8 time steps (2 years)
Csum = apply(Csub,1:2,sum) # sum over the first two dimensions (timestep and fleet)
kable(round(Csum/1000,0)) # tabulate to the nearest ton| USCom | RecN | RecS | HireN | HireS | Intl | Disc | UnRep | |
|---|---|---|---|---|---|---|---|---|
| 2021.0 | 2 | 0 | 149 | 0 | 13 | 506 | 24 | 92 |
| 2021.25 | 69 | 447 | 1045 | 221 | 252 | 2686 | 273 | 628 |
| 2021.5 | 23 | 124 | 620 | 87 | 124 | 1338 | 133 | 325 |
| 2021.75 | 18 | 190 | 490 | 133 | 241 | 803 | 140 | 199 |
| 2022.0 | 3 | 0 | 94 | 0 | 13 | 440 | 16 | 96 |
| 2022.25 | 75 | 332 | 1349 | 90 | 156 | 2065 | 279 | 463 |
| 2022.5 | 18 | 27 | 1242 | 68 | 130 | 1923 | 211 | 466 |
| 2022.75 | 9 | 284 | 318 | 43 | 29 | 905 | 98 | 219 |
10.4 Length Classes
Management options such as post-release mortality rate are available by length class to model the impact of alternative gears or regulations. Length class information is available in the operating model objects (more on those below) but since it is the same among OMs, by default a vector is loaded with the mahiMSE to make these easy to locate:
Length_classes
#> [1] 22.23 66.69 111.15 155.61 200.07 244.53 288.99 333.45 377.91 422.37 466.83 511.29 555.75 600.21 644.67 689.13 733.59 778.05 822.51 866.97 911.43 955.89 1000.35 1044.81 1089.27 1133.73 1178.19 1222.65 1267.11 1311.57 1356.03 1400.49 1444.95 1489.41 1533.87 1578.33 1622.79 1667.25 1711.7111 Operating Models
11.1 Operating Model Structure
MSEtool v4.x.x operating models are of class ‘om’ and there is full documentation from the R command line:
class?om
Operating models contain all information about the dynamics for the stock(s), fleets, observation (data simulation) processes and management implementation.
slotNames(smallOM)
?Stock
?Fleet
These Stock and Fleet functions are used to populate slots in the om objects and have various sub components:
smallOM@Stock$Dolphinfish@Length
#>
#> ── A `?MSEtool::length-class()` Object ──
#>
#> `Pars`:
#> `Model`:
#>
#> `Units`: not specified
#> `MeanAtAge`: "8 Sim x 17 Age x 228 Year array"
#> `CVatAge`: "1 Sim x 1 Age x 1 Year array"
#> `Dist`: "normal"
#> `TruncSD`: 2
#>
#> `Timing`: 0
#> `Random`: not specified
#> `ALK`: "8 Sim x 17 Age x 39 Class x 228 Year array"
#> `Classes`: 22.2, 66.7, 111, 156, 200, 245, 289, 333, 378, 422, 467, 511, 556, 600, 645, 689, 734, 778, …, 1670, and 171011.2 Reference Set Operating Models
The reference set is a fully orthogonal grid with two levels for each of five factors (see the TSD for more info):
| Uncertainty | Level 1 | Level 2 |
|---|---|---|
| Natural Mortality | Low (m): 0.25 per season | High (M): 0.5 per season |
| Resilience (steepness) | Low (s): 0.7 | High (S): 0.95 |
| Catch levels of Disc and Unrep | Low (c): half NOAA calculated | High (C):100% NOAA calculated |
| Recruitment Level | Low (r) 75% historical mean levels | High (R) 100% historical mean levels |
| Viscosity | Low (v): prob. stay. = 0.4 | High (V) prob. stay. = 0.6 |
These 32 operating models are organized in the general order of most (OM_1) to least (OM_32) challenging test of an MP:
| Code | M | Steep | Cat | Rec | Visc |
|---|---|---|---|---|---|
| OM_1 | 0.25 | 0.70 | 0.5 | 0.75 | 60 |
| OM_2 | 0.50 | 0.70 | 0.5 | 0.75 | 60 |
| OM_3 | 0.25 | 0.95 | 0.5 | 0.75 | 60 |
| OM_4 | 0.50 | 0.95 | 0.5 | 0.75 | 60 |
| OM_5 | 0.25 | 0.70 | 1.0 | 0.75 | 60 |
| OM_6 | 0.50 | 0.70 | 1.0 | 0.75 | 60 |
| OM_7 | 0.25 | 0.95 | 1.0 | 0.75 | 60 |
| OM_8 | 0.50 | 0.95 | 1.0 | 0.75 | 60 |
| OM_9 | 0.25 | 0.70 | 0.5 | 1.00 | 60 |
| OM_10 | 0.50 | 0.70 | 0.5 | 1.00 | 60 |
| OM_11 | 0.25 | 0.95 | 0.5 | 1.00 | 60 |
| OM_12 | 0.50 | 0.95 | 0.5 | 1.00 | 60 |
| OM_13 | 0.25 | 0.70 | 1.0 | 1.00 | 60 |
| OM_14 | 0.50 | 0.70 | 1.0 | 1.00 | 60 |
| OM_15 | 0.25 | 0.95 | 1.0 | 1.00 | 60 |
| OM_16 | 0.50 | 0.95 | 1.0 | 1.00 | 60 |
| OM_17 | 0.25 | 0.70 | 0.5 | 0.75 | 40 |
| OM_18 | 0.50 | 0.70 | 0.5 | 0.75 | 40 |
| OM_19 | 0.25 | 0.95 | 0.5 | 0.75 | 40 |
| OM_20 | 0.50 | 0.95 | 0.5 | 0.75 | 40 |
| OM_21 | 0.25 | 0.70 | 1.0 | 0.75 | 40 |
| OM_22 | 0.50 | 0.70 | 1.0 | 0.75 | 40 |
| OM_23 | 0.25 | 0.95 | 1.0 | 0.75 | 40 |
| OM_24 | 0.50 | 0.95 | 1.0 | 0.75 | 40 |
| OM_25 | 0.25 | 0.70 | 0.5 | 1.00 | 40 |
| OM_26 | 0.50 | 0.70 | 0.5 | 1.00 | 40 |
| OM_27 | 0.25 | 0.95 | 0.5 | 1.00 | 40 |
| OM_28 | 0.50 | 0.95 | 0.5 | 1.00 | 40 |
| OM_29 | 0.25 | 0.70 | 1.0 | 1.00 | 40 |
| OM_30 | 0.50 | 0.70 | 1.0 | 1.00 | 40 |
| OM_31 | 0.25 | 0.95 | 1.0 | 1.00 | 40 |
| OM_32 | 0.50 | 0.95 | 1.0 | 1.00 | 40 |
The reference set of operating models are loaded with the mahiRefSet package:
library(mahiRefSet)
OM_32@Name
#> [1] "RefOM #32: 0.5-0.95-1-1-40"11.3 Reference Case Operating Model
In many MSE processes a Reference Case operating model is identified that represents a ‘default’ set of assumptions. This operating models provides a basis for concise testing of sensitivity to alternative assumptions and quick comparison of candidate MPs. The current Reference Case operating model is OM_29 that has low M, low steepness, full historical catches, continuation of historical mean recruitment and high stock mixing.
11.4 Robustness Set Operating Models
The robustness set consists of second-order scenarios for exploring MP performance and further discriminating among MPs that perform similarly for the Reference Set. The robustness operating models are based on the Reference Case OM (OM_29).
| Code | Type | Description |
|---|---|---|
| ROM_C1 | Catches | IUU increases by 1% every year |
| ROM_C2 | Catches | Catches halved for non-US fleets in area 1 (Florida and Caribbean) |
| ROM_R1 | Recruitment | Future recruitment declines 1% per year |
| ROM_R2 | Recruitment | Future recruitment reduces by 25% after 5 years |
| ROM_R3 | Recruitment | Future recruitment reduces by 25% after 10 years |
| ROM_R4 | Recruitment | Future recruitment is 50% more variable |
| ROM_S1 | Spatial | Two percent decline in CAR and SFL, 1 percent increase in SE, 2% increase in NC and NE |
| ROM_S2 | Spatial | 50% greater variability in spatial / seasonal distribution |
| ROM_S3 | Spatial | 1% pa. increase in catchability reflecting range contraction |
| ROM_S4 | Spatial | 20% target probability of staying |
| ROM_S5 | Spatial | 80% target probability of staying |
| ROM_P1 | Productivity | 1% pa. decrease in somatic growth rate (k) |
| ROM_P2 | Productivity | 1% pa. decrease in condition factor (weight at length) |
| ROM_P3 | Productivity | 1% pa. increase in natural mortality rate (all ages) |
| ROM_F1 | Fishery | 33% shorter selectivity of all fleets (inflection point) |
| ROM_F2 | Fishery | Projected shift towards 33% younger selectivity of all fleets (inflection point) |
The Robustness Set of operating models are loaded with the mahiRobSet package:
library(mahiRobSet)
ROM_S2@Name
#> [1] "ROM_S2: 0.25-0.7-1-1-40"12 The mahiMP
The mahiMP is a function of openMSE class ‘mp’ that is compatible with the operating models of the mahiMSE package (MSEtool v4.x.x, openMSE v2.0).
The function is designed to mimic the current status quo management procedure for dolphinfish that includes:
a minimum size limit of 500mm for the RecS and HireS fleets;
a trip limit of 54 fish for the RecS, HireS, RecN and HireN fleets;
a US commercial (fleet USCom) ACL of 1,719,953 lbs whole weight and
since it may be affected by management regulations, post release mortality is also a management lever of the management procedure and by default assumes a 25% rate (one in four released fish subsequently dies). This can be changed for all length classes together for any fleet or be changed for individual length classes within fleet.
potential effort frozen at current levels (this is necessary because fleets are not reaching their catch limits, without an effort restriction the fleets would immediately ‘jump up’ in exploitation rate in the projections to attempt to fish their catch limit, which is not realistic given the fishery over the last few years)
The MP also includes a number of other possible management levers including:
Total allowable effort (i.e. for exploring possible changes in access or fleet capacity) and
maximum size limits to allow for preservation of larger spawning fish
In addition to these levers the MP can include empirical harvest control rules for TAC, minimum size and TAE. These harvest control rules can be linked to fishery-dependent indices and assigned control parameters that impact exploitation rate as stock levels decline. These empirical harvest control rules can be stacked, allowing for simultaneous management of TACs, size limits and TAE by fleet or in total.
The mahiMP comes with R help documentation that can be accessed from the command line:
?mahiMP
In this document I go into further detail on how to implement various management options and conduct closed-loop simulations to test their efficacy.
12.1 Catch Limit Control
12.1.1 TAC
The TAC argument of the mahiMP specifies the absolute annual TAC taken by each fleet (it is a vector nfleets long). This annual catch is redistributed to seasons and areas according to the proportions of the most recent year of the operating model (2022) (object rTAC). In future editions of mahiMP you will be able to specify TAC by fleet x area and fleet x area x season.
In this example we will reduce the international fleet to 1kt.
TAC_mp = mahiMP # copy mahiMP
nfleets = length(Fleets) # how many fleets are there?
formals(TAC_mp)$TAC = apply(rTAC,2,sum) # total annual catches of the various fleets in 2022
Intl_ind = match("Intl", Fleets) # which is the international fleet?
formals(TAC_mp)$TAC[Intl_ind] = 1E6 # 1m kg (1kt
class(TAC_mp) = 'mp' # assign class mp
ad = TAC_mp(Example_Data) # run for an example dataset
ad@TAC # the seasonal TAC specified
#> USCom RecN RecS HireN HireS Intl Disc UnRep
#> 74982.31 331886.22 1348700.11 90497.15 155908.78 2065400.78 278978.01 462802.6212.1.2 rel_TAC
Rather than specify the absolute TAC it is possible to specify relative TAC (relative to 2022) as either a fraction for all fleets, or a vector of fractions nfleets long. For example, if you wished to turn off the international fleet altogether:
TACrel_mp = mahiMP # copy mahiMP
rel_TAC = rep(1,length(Fleets)) # set relative effort to 100% all fleets
rel_TAC[match("Intl", Fleets)] = 0 # turn off International fleet
formals(TACrel_mp)$rel_TAC = rel_TAC # set relative TAC for fleets
class(TACrel_mp) = 'mp' # assign class mp
ad = TACrel_mp(Example_Data) # run for an example dataset
ad@TAC # the seasonal TAC specified
#> USCom RecN RecS HireN HireS Intl Disc UnRep
#> 74982.31 331886.22 1348700.11 90497.15 155908.78 0.00 278978.01 462802.6212.2 Minimum Size Limits
The current default size limits for mahiMP are:
Fleets
#> [1] "USCom" "RecN" "RecS" "HireN" "HireS" "Intl" "Disc" "Unrep"
formals(mahiMP)$SLmin
#> c(0, 0, 500, 0, 500, 0, 0, 0)Minimum size limits can be altered similarly to other mahiMP arguments. Just for purely demonstration purposes, lets add the 500mm minimum size to the northern recreational and for-hire fleets (RecN and HireN, respectively)
SLminMP = mahiMP
SLmin = rep(0,length(Fleets))
SLmin[Fleets %in% c("RecN", "RecS","HireN","HireS")] = 500
formals(SLminMP)$SLmin = SLmin
class(SLminMP) = 'mp'
ad0 = mahiMP(Example_Data) # get advice for the base mahiMP
ad = SLminMP(Example_Data) # run for an example dataset
fleet = match("RecN",Fleets) # locate the northern rec fleet
ad0@Retention[[fleet]]@MeanAtLength # length class x area
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 1 1 1 1 1
#> [2,] 1 1 1 1 1
#> [3,] 1 1 1 1 1
#> [4,] 1 1 1 1 1
#> [5,] 1 1 1 1 1
#> [6,] 1 1 1 1 1
#> [7,] 1 1 1 1 1
#> [8,] 1 1 1 1 1
#> [9,] 1 1 1 1 1
#> [10,] 1 1 1 1 1
#> [11,] 1 1 1 1 1
#> [12,] 1 1 1 1 1
#> [13,] 1 1 1 1 1
#> [14,] 1 1 1 1 1
#> [15,] 1 1 1 1 1
#> [16,] 1 1 1 1 1
#> [17,] 1 1 1 1 1
#> [18,] 1 1 1 1 1
#> [19,] 1 1 1 1 1
#> [20,] 1 1 1 1 1
#> [21,] 1 1 1 1 1
#> [22,] 1 1 1 1 1
#> [23,] 1 1 1 1 1
#> [24,] 1 1 1 1 1
#> [25,] 1 1 1 1 1
#> [26,] 1 1 1 1 1
#> [27,] 1 1 1 1 1
#> [28,] 1 1 1 1 1
#> [29,] 1 1 1 1 1
#> [30,] 1 1 1 1 1
#> [31,] 1 1 1 1 1
#> [32,] 1 1 1 1 1
#> [33,] 1 1 1 1 1
#> [34,] 1 1 1 1 1
#> [35,] 1 1 1 1 1
#> [36,] 1 1 1 1 1
#> [37,] 1 1 1 1 1
#> [38,] 1 1 1 1 1
#> [39,] 1 1 1 1 1
ad@Retention[[fleet]]@MeanAtLength # length class x area
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 0 0 0 0 0
#> [2,] 0 0 0 0 0
#> [3,] 0 0 0 0 0
#> [4,] 0 0 0 0 0
#> [5,] 0 0 0 0 0
#> [6,] 0 0 0 0 0
#> [7,] 0 0 0 0 0
#> [8,] 0 0 0 0 0
#> [9,] 0 0 0 0 0
#> [10,] 0 0 0 0 0
#> [11,] 0 0 0 0 0
#> [12,] 1 1 1 1 1
#> [13,] 1 1 1 1 1
#> [14,] 1 1 1 1 1
#> [15,] 1 1 1 1 1
#> [16,] 1 1 1 1 1
#> [17,] 1 1 1 1 1
#> [18,] 1 1 1 1 1
#> [19,] 1 1 1 1 1
#> [20,] 1 1 1 1 1
#> [21,] 1 1 1 1 1
#> [22,] 1 1 1 1 1
#> [23,] 1 1 1 1 1
#> [24,] 1 1 1 1 1
#> [25,] 1 1 1 1 1
#> [26,] 1 1 1 1 1
#> [27,] 1 1 1 1 1
#> [28,] 1 1 1 1 1
#> [29,] 1 1 1 1 1
#> [30,] 1 1 1 1 1
#> [31,] 1 1 1 1 1
#> [32,] 1 1 1 1 1
#> [33,] 1 1 1 1 1
#> [34,] 1 1 1 1 1
#> [35,] 1 1 1 1 1
#> [36,] 1 1 1 1 1
#> [37,] 1 1 1 1 1
#> [38,] 1 1 1 1 1
#> [39,] 1 1 1 1 1Lets do an MSE projection comparing these two MPs:
Hist = Simulate(smallOM)
Proj_sz = Project(Hist, c("mahiMP","SLminMP"))
mahiplot(Proj_sz)
12.3 Trip Limits
Trip level catch rate data (2018-2022) are only available for the southern recreational and for-hire fleets (RecS and HireS, respectively) in certain areas and quarters. These catch rates are used to predict how often a trip may exceed the trip limit and therefore be restricted to that quantity of landings (providing a prediction of the retention rate).
There are two approaches to this. The first is to fit a statistical distribution to the observed catch rate data and then use this to predict the distribution of catch rates as regional abundance fluctuates (and hence predict the retention of fish subject to the trip limit).
Below, you can see that the fit of a negative binomial model (red line) is sometimes not a great characterization of the data (black line) particularly in the upper tail where trip-limit impacts on retention are likely to be most important:
Rather than use a statistical distribution it is possible to use the empirical distribution directly to predict catch rates, essentially stretching the empirical distribution (increasing / decreasing the mean catch rate while keeping the coefficient of variation the same) to predict retention subject to the trip limit.
A comparison of the approaches shows key differences in the predicted fraction of trips above the 54 fish trip limit and a larger discrepancy in the predicted fraction of caugth fish discarded:
Trip limits are currently specified by fleet in the mahiMP and are straightforward to specify. If you apply it to a data object with plot=T, you can see how retention was calculated (the dark green line is the observed distribution of catch rates from 2018-2022, the red line is that empirical distribution ‘stretched’ by the ratio of vulnerable biomass now to 2018-2022 (VB ratio), the expected retention (fraction of fish kept, 1- discard rate) is included in each plot (‘Retn’):
TL10 = mahiMP
formals(TL10)$TL = c(NA, 10, 10, 10, 10, NA, NA, NA) # ten fish for all hire and rec fleets
class(TL10) = 'mp'
par(mfrow=c(2,3),mai=c(0.6,0.6,0.3,0.1))
ad1 = mahiMP(Example_Data,plot=T)
Now lets see how the trip limit of 10 fish changes the retention calculation:
par(mfrow=c(2,3),mai=c(0.6,0.6,0.3,0.1))
ad2 = TL10(Example_Data, plot=T)
12.4 Harvest Control Rules
Harvest Control Rules alter management advice in response to data regarding abundance. In a conventional data-rich stock assessment, an HCR has an independent variable reflecting stock status (e.g., estimated SSB relative to SSB(MSY)) and a dependent variable which is typically the estimated exploitation rate arising from the HCR (e.g., F relative FMSY).
Management procedures usually attempt to use empirical harvest control rules that work on the same principle but use observed data directly. For example an abundance index relative to a target level, catch or catch/index as the dependent variable.
The mahiMP allows for harvest control rules that can control fleet-specific effort or catch from indices that can be spatially aggregated (e.g., total VAST), spatial or fleet specific (vulnerable biomass, ie catch per unit effort based indices).
The default functional form of HCR in mahiMP is the ‘hockey-stick’ that has two control points. The dependent variable control points (LCP) provide the range of management advice for a particular lever. The independent variable control points (ICP) determine up to what point advice is constant at the lower LCP and above which advice is constant at the upper LCP. At intermediate levels between ICP control points, advice is a linear increase from the lower to upper LCP.
Here is an example of a hockey stock HCR (solid black line) from an example below, the inflection points of the HCR are the control points:
12.4.1 Index-Based Empirical MP with Constant Exploitation Rate
Lets start by specifying an HCR that continues to fish at the current catch/index by modifying TAC advice in proportion to an overall index of abundance. Because this is a linear increase in TAC with index it has lower ICP and LCP at the origin [0,0]. We set a max TAC to twice that of current levels (defined by mahiMP HCR_calib argument, by default last four seasons in 2022, time steps 145:148):
TAC_HCR = mahiMP
formals(TAC_HCR)$HCR_use = T # a switch that turns on/off any specified HCR
formals(TAC_HCR)$HCR_ICP = list(c(0,2)) # control points for the input index, 0 to 2x current levels (index)
formals(TAC_HCR)$HCR_LCP = list(c(0,2)) # corresponding lever control points, 0 to 2x current levels (catch)
formals(TAC_HCR)$HCR_calib = 145:148 # The reference time steps (what is current level) (the last 4 seasons, ie 2022)
formals(TAC_HCR)$HCR_fleet = NA # NA means that all fleets will have their TAC changed relative to current levels
formals(TAC_HCR)$HCR_sens = 1 # changes are proportional and not damped
formals(TAC_HCR)$HCR_lever = "TAC" # TAC is the management lever
formals(TAC_HCR)$HCR_up_max = 0.25 # maximum upward changes of 25%
formals(TAC_HCR)$HCR_down_max = 0.25 # maximum downward changes of 25%
formals(TAC_HCR)$HCR_input = 1 # The input is index #1
formals(TAC_HCR)$HCR_smooth = 0.2 # The HCR will used a filtered index that is a polynomial smoother with enp = 0.2
formals(TAC_HCR)$Index_type = "VAST" # The overall VAST index will be index #1
formals(TAC_HCR)$Index_cv = 0.1 # Future observation error is a lognormal sd of 0.1 (CV 10%)
formals(TAC_HCR)$Index_ac = 0 # No lag-1 autocorrelation in the index
formals(TAC_HCR)$Index_lag = 1 # Data available for the year before
formals(TAC_HCR)$Index_seed = 1 # Seed controlling sampling of error structures
par(mfrow=c(3,3),mai=c(0.6,0.6,0.3,0.1))
ad = TAC_HCR(Example_Data, plot=T) # Use plot to see internal calculations
12.5 Alternative Control Points
To demonstrate how control point are specified, lets modify the TAC_HCR MP and recalculate advice. Note that control points are vectors in a list - it has to be a list so that you can stack multiple HCRs, for example rules for the various fleets.
TAC_HCR2 = TAC_HCR # copy previous MP
formals(TAC_HCR2)$HCR_ICP = list(c(0.5,1)) # hockey stick control points at 50% and 100% current index levels
formals(TAC_HCR2)$HCR_LCP = list(c(0.2,1.5)) # hockey stick control points at 20% and 150% current catch per index
class(TAC_HCR2) = "mp"
par(mfrow=c(3,3),mai=c(0.6,0.6,0.3,0.1))
ad = TAC_HCR2(Example_Data, plot=T) # Use plot to see internal calculations
HCR plot of this figure shows that the HCR alone would have recommended a substantial increase in TAC (blue arrow) but this was subsequently adjusted downwards because of the TAC change constraints, namely HCR_up_max which is set to 0.25 (25%) (red arrow).
12.6 Use of Smoothers
It is desirable to apply smoothers to the independent variable (index) of an HCR to filter signal from noise. The default approach of mahiMSE is a polynomial smoother parameterized by ‘effective number of parameters’ (enp). The HCR_smooth argument sets enp as a multiple of the length of the time series, thereby keeping the ratio of polynomial parameters to data inputs the same. The value above of 0.2 would assign 20 parameters for an index time series of 100 timesteps. The higher the fraction the more parameters and the lower the smoothing.
For example:
TAC_HCR3 = TAC_HCR2 # copy previous MP
formals(TAC_HCR3)$HCR_smooth = 0.05 # effective number of parameters is 5% the length of the time series
class(TAC_HCR3) = "mp"
par(mfrow=c(3,3),mai=c(0.6,0.6,0.3,0.1))
ad = TAC_HCR3(Example_Data, plot=T) # Use plot to see internal calculations
The smoother panel of this figure shows how responsive the smoothed index is with respect to the data. The previous figure used the default value of 0.2 and tracked the index more closely over shorter time scales. You can see that the independent variable of the HCR has changed also and is no longer leading to a recommendation outside of the maximum increase in TAC (there is only a blue arrow).
Ultimately the choice of smoothing level come down to management performance: more smoothing is likely to reduce variability in advice but in doing so lower responsiveness to available biomass reducing average yields.
Like all mahiMP code the smoother can be used independently to better understand its properties:
smth = mahiMSE::Index_smooth
set.seed(3)
ndat = 100
fakedata = exp(0.2+sin(seq(0,30,length.out=ndat))*exp(rnorm(ndat,0,0.6)))
HCR_smooth = c(0.05,0.1,0.2,0.3)
smdat = sapply(HCR_smooth,function(x,fakedata)smth(fakedata,enp.mult=x),fakedata=fakedata)
plot(fakedata); matplot(smdat,type="l",lty=1,lwd=2,col=c("red","green","blue","darkgrey"),add=T)
12.7 Management Lever Sensitivity
Rather than setting abrupt limits on HCR recommendation changes, it might be desirable to dampen changes. Doing so might prevent the MP from constantly hitting max HCR change constraints (HCR_up_max, HCR_down_max) and lessen variability in recommendations for the same yield and biomass outcomes.
The functional form of the smoother is very simple. If the recommended TAC change from the HCR alone is:
\[ \delta_{y} = TAC_{y}/TAC_{y-1} \] Then the smoother adjusted change is:
\[ \hat{\delta}_{y} = exp(log(\delta_{y})\theta_{HCRsens}) \]
Here are the properties of various levels of the sensitivity argument HCR_sens:
delta = exp(seq(-1.5,0.6,length.out=1000))
HCR_sens = c(0.2,0.5, 1, 1.25, 1.5)
cols = c("red","green","blue","black","orange")
sensdat = sapply(HCR_sens,function(x,delta)exp(log(delta)*x),delta=delta)
matplot(delta, sensdat,type="l",xlab = "Change directly from HCR (delta)",
ylab="After sensitivity adjustment (delta_hat)",lty=1,lwd=2,
col=cols); grid()
legend('topleft',legend=HCR_sens,text.col=cols)
When HCR_sens has a value of 1, the output of the HCR passes the sensitivity function unchanged (the blue line is on the 1:1 ratio line). At values lower than 1 (green and red lines), the signal is dampened and recommendations change slower than recommended by the HCR. At values of HCR_sens above 1 recommendations arising from the sensitivity function are exaggerated.
Although similar signal damping can be achieved by using the change constraints or the gradient of the HCR itself, in other MSE settings such as Atlantic Bluefin Tuna, the sensitivity parameter could be adjusted to greatly reduce MP advice variability with little cost to yield and biomass outcomes and as such offers an important avenue for MP investigation.
12.8 Fleet-Specific Advice from and Fleet-Specific Data
Lets assume that you wish to alter the ACL for the RecS fleet using a vulnerable biomass index (e.g., a fishery dependent standardized CPUE index). We will change two arguments argument of the TAC_HCR function to achieve this:
TAC_HCR4 = TAC_HCR # copy previous linear increasing TAC with index (constant harvest rate) MP
formals(TAC_HCR4)$HCR_fleet = list(3) # RecS is Fleets[3]
formals(TAC_HCR4)$HCR_input = 1 # Use index 1
formals(TAC_HCR4)$Index_type = "VB3" # Vulnerable biomass index for fleet 3
class(TAC_HCR4) = "mp"
par(mfrow=c(3,3),mai=c(0.6,0.6,0.3,0.1))
ad = TAC_HCR4(Example_Data, plot=T)
12.9 Stacking HCRs
But what if you wanted to manage the commercial fleet with a commercial vulnable biomass index, the RecS fleet with a RecS vulnerable biomass index and the HireS, HireN and RecN with a VAST index rule?
Each HCR slot now has three positions. You’ll notice there are three indices but there doesn’t have to be, each of these HCRs could use the same index (if this was index 1 the argument HCR_input would be c(1,1,1).
Here is an example:
# A three HCR system by fleet or groups of fleets
TAC_HCR = mahiMP
formals(TAC_HCR)$HCR_use = c(T, T, T) # a switch that turns on/off any specified HCR
formals(TAC_HCR)$HCR_ICP = list(c(0,2), c(0,2), c(0,2)) # control points for the input index, 0 to 2x current levels (index)
formals(TAC_HCR)$HCR_LCP = list(c(0,2), c(0,2), c(0,2)) # corresponding lever control points, 0 to 2x current levels (catch)
formals(TAC_HCR)$HCR_fleet = list(1, 3, c(2,4,5)) # UScom, RecS, c(RecN, HireN, HireS)
formals(TAC_HCR)$HCR_sens = c(1,1,1) # changes are proportional and not damped
formals(TAC_HCR)$HCR_lever = c("TAC","TAC","TAC") # TAC is the management lever
formals(TAC_HCR)$HCR_up_max = c(0.25, 0.25, 0.25) # maximum upward changes of 25%
formals(TAC_HCR)$HCR_down_max = c(0.25, 0.25, 0.25) # maximum downward changes of 25%
formals(TAC_HCR)$HCR_input = c(1,2,3) # The input indices (defined below)
formals(TAC_HCR)$HCR_smooth = c(0.2,0.2,0.2) # The HCR filtered index that is a polynomial smoother with enp = 0.2
formals(TAC_HCR)$Index_type = c("VB1", "VB3","VAST") # (1) vulnB for fleet 1, (2) vulnB for fleet 3, (3) overall VAST index
formals(TAC_HCR)$Index_cv = c(0.1,0.1,0.1) # Future observation error is a lognormal sd of 0.1 (CV 10%)
formals(TAC_HCR)$Index_ac = c(0, 0, 0) # No lag-1 autocorrelation in the index
formals(TAC_HCR)$Index_lag = c(1, 1, 1) # Data available for the year before
formals(TAC_HCR)$Index_seed = c(1, 1, 1) # Seed controlling sampling of error structures
par(mfrow=c(4,4),mai=c(0.6,0.6,0.3,0.1))
ad = TAC_HCR(Example_Data, plot=T) # Use plot to see internal calculations
ad@TAC
#> USCom RecN RecS HireN HireS Intl Disc UnRep
#> 74982.31 331886.22 1348700.11 90497.15 155908.78 2065400.78 278978.01 462802.62
12.10 HCRs that Control Size Limits
Coming soon!
12.11 HCRs that Control Trip Limits
Coming soon!
13 Tuning MPs
Coming soon!
14 MSE Performance Presentation in Slick
Coming soon!
15 ECP
Coming soon!
16 Resources
16.1 splash page
For more information about how operating models were specified you can visit the project webpage
Contact (tom@bluematterscience for password)
16.2 Dolphinfish GitHub repository
All of the code for conditioning operating models can be found at a private GitHub repository (contact Cassidy Peterson for access)
16.3 Inverts R packages
The code for the R packages used in this manual are available from the public mahiMSE GitHub repository:
16.4 OpenMSE
Visit the OpenMSE website for more information on the framework including tutorial and case studies.
The Rapid Conditioning Model is also described there
Other useful openMSE apps include MERA for rapidly scoping operating models and Reference Point Calculator for investigating reference points and control rules.
16.5 Management Strategy Evaluation
The Ocean Foundation host the Harvest Strategies website which is an excellent resource describing MSE, MSE concepts and how MSE can be used in fishery management including introductory videos and interviews with experts. ‘Harvest strategy’ and Management Procedure are synonymous by the way.
17 References
Hordyk, A., Huynh, Q., Carruthers, T. 2025. OpenMSE: An open-source R package for Management strategy evaluation, available from: https://openmse.com
Huynn, Q., 2025. Rapid Conditioning model. Available from https://openmse.com/tutorial-rcm/
Punt, A.E., Butterworth, D.S., de Moor, C.L., De Oliveira, J.A.A., and Haddon, M. 2016. Management strategy evaluation: Best practices. Fish Fish. 17(2): 303–334. doi:10.1111/faf.12104.