Inverts User Manual
An R package for conditioning operating models and conducting MSE for informing management of hand harvested invertebrate fisheries in B.C.
Tom Carruthers (tom@bluematterscience.com)
2025-12-19
1 Foreword
Fisheries for hand-harvested invertebrates (HHI) present opportunities and challenges for science and management that are often not pertinent to fisheries for pelagic finfish.
Populations at the scale of the B.C. coastline are likely to have regional sub-population dynamics at varying spatial scales with uncertain larval dispersal and exchange of biomass. Data on HHI may come from a mix of commerical, recreational and organized scientific surveys that are often not available for earlier time periods of substantial exploitation.
While aspects of population dynamics and structure may be highly uncertain, the prescriptive nature of management measures (e.g. minimum size limits, rotational spatial closures, gear restrictions) for HHI could be highly robust and support the sustainable exploitation of productive stocks.
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 describe the decision support tools and guide prospective users through various aspects of MSE analyses providing examples on how to:
- specify and condition operating models
- calculate reference points
- modify operating models
- define and calculate management performance metrics
- design and test management procedures
- produce their own plots and diagnostics
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 Introduction
4.1 Hand Harvested Invertebrates
As part of the approach to meeting the Fish Stocks provisions (FSP) in the revised Fisheries Act, DFO’s Marine Invertebrate Section (MIS) has been adapting the Management Procedure (MP) Framework (Anderson et al. 2021), a decision support tool, to support investigation of management decisions related to conservation objectives for Pacific marine invertebrates. This project will result in a decision support tool for exploring impacts of spatial closures (Marine Spatial Planning) and other management actions, in order to better understand trade-offs for decision makers under future marine ecosystem and human use (primarily fishing) scenarios.
These libraries adapt existing tools of the MP Framework to: define reference case and robustness operating models (OMs); define quantitative performance metrics; test alternative management procedures, calculate appropriate reference points, with management options relevant for marine invertebrates; and determining additional management objectives beyond the FSP for the four case studies (Geoduck, intertidal clams, Green Sea Urchin, and Giant Red Sea Cucumber).
4.2 Background information
Project information, links and results are available from a project webpage
4.3 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.
Reference Point Calculator (Huynh et al. ) was used to explore fishery reference points
Objects and code are organized into a generic package ‘Inverts’ and a series of species-specific packages (e.g., Inverts.GD, Inverts.MC, Inverts.GSU, Inverts.SC for geoduck, manila clam, green sea urchin and sea cucumber, respectively).
Organizing the packages this way means that unnecessary objects do not need to be loaded when working on a particular species. Additionally the framework can be scaled to new species and allows for individual improvements to particular species without impacting the other packages.
| Inverts | The umbrella R package containing generic code and example objects for each case study stock |
| Inverts.GD | Data, operating models and tools specific to Geoduck |
| Inverts.MC | Data, operating models and tools specific to Manila Clam |
| Inverts.GSU | Data, operating models and tools specific to Green Sea Urchin |
| Inverts.SC | Data, operating models and tools specific to Giant Red Sea Cucumber |
Data, operating models, analyses and R packages are stored on the following GitHub repositories
| mis-assess/csrf_hh_data | Raw data, background analyses, supporting documentation for project | Private |
| mis-assess/csrf_hh_io | Project splash page | Public |
| blue-matter/Inverts | Inverts R package | Public |
| blue-matter/Inverts.GD | Inverts.GD R package | Public |
| blue-matter/Inverts.MC | Inverts.MC R package | Public |
| blue-matter/Inverts.GSU | Inverts.GSU R package | Public |
| blue-matter/Inverts.SC | Inverts.SC R package | Public |
5 Installation
Install R for Windows
Install RStudio
Install openMSE and devtools packages
Open RStudio and enter the following commands into the console:
install.packages('openMSE')
install.packages('devtools')- Install github packages including the Inverts and RPC R packages by opening RStudio and entering the following commands at the R console :
devtools::install_github("blue-matter/RPC")
devtools::install_github("blue-matter/Inverts") - Install any species specific packages
devtools::install_github("blue-matter/Inverts.GD") # Geoduck
devtools::install_github("blue-matter/Inverts.MC") # Manila clam- Check that the installation is successful by finding this help file (the Inverts vignette):
library(Inverts)
vignette('Inverts')If you have any difficulties please send an inquiry including some reproducible code to: tom@bluematterscience.com
6 Quick start: run an MSE with pre-specified objects
Load libraries:
library(Inverts)
library(Inverts.MC)Locate a stock to work with:
objs('In') # objects of class input ('In')
#> Searching for objects of class In in package: Inverts
#> Searching for objects of class In in package: Inverts.GD
#> Searching for objects of class In in package: Inverts.MC
#> [1] "In.GD.1" "In.GD.106" "In.GD.12" "In.GD.13" "In.GD.14" "In.GD.15"
#> [7] "In.GD.16" "In.GD.17" "In.GD.18" "In.GD.19" "In.GD.2" "In.GD.23"
#> [13] "In.GD.24" "In.GD.25" "In.GD.26" "In.GD.27" "In.GD.3" "In.GD.4"
#> [19] "In.GD.5" "In.GD.6" "In.GD.7" "In.GD.8" "In.MC.C" "In.MC.D"
#> [25] "In.MC.E" "In.MC.F"
about(In.MC.E) # where is this?
Condition an operating model for Manila Clam in CMA E and plot the fit:
myfit = cond.MC(In.MC.E)
plot(myfit)
Build the historical simulation:
Hist = runMSE(myfit@OM, Hist=T)Calculate historical levels relative to reference points:
RPs = Ref.Points(Hist)
Quantify current stock status:
Status(Hist)
#> $SSB_SSB0
#> [1] 0.4610 0.4896 0.6213 0.5724 0.7093 0.5927 0.5080 0.5980 0.5586 0.5436
#> [11] 0.6989 0.5919 0.5011 0.4337 0.7043 0.6705 0.5918 0.6463 0.5365 0.6874
#> [21] 0.5946 0.5952 0.6089 0.4805 0.5515 0.5904 0.4890 0.4778 0.6299 0.7176
#> [31] 0.6209 0.6538
#>
#> $SSB_SSBMSY
#> [1] 1.2762 1.2438 1.5581 1.5906 1.5521 1.6405 1.2987 1.7274 1.4393 1.4272
#> [11] 1.4731 1.8246 1.3928 1.1431 1.5741 1.6345 1.5903 1.3602 1.4067 1.4718
#> [21] 1.5815 1.4919 1.6031 1.1908 1.3907 1.4854 1.2399 1.2056 1.6559 1.6081
#> [31] 1.5366 1.6894
#>
#> $F_F01
#> [1] 0.3575 0.2505 0.1381 0.2268 0.0649 0.1634 0.2194 0.1587 0.1796 0.1931
#> [11] 0.0831 0.1666 0.2762 0.3331 0.0771 0.1046 0.1457 0.1321 0.2389 0.0800
#> [21] 0.2065 0.1275 0.1467 0.2522 0.1963 0.1395 0.2543 0.2591 0.1276 0.0715
#> [31] 0.1568 0.1234
#>
#> $F_FMSY
#> [1] 0.1596 0.1241 0.0404 0.0751 0.0265 0.0421 0.1011 0.0425 0.0770 0.0898
#> [11] 0.0360 0.0507 0.0863 0.1719 0.0251 0.0281 0.0523 0.0428 0.1102 0.0308
#> [21] 0.0487 0.0523 0.0353 0.1437 0.0932 0.0747 0.1344 0.1815 0.0318 0.0231
#> [31] 0.0400 0.0306Simulation test status quo fishing pressure / catches and reference fishing:
myMSE = Project(Hist, c('MP.MC', # status quo manila clam management
"FMSYref50", # perfect implementation of 50% FMSY fishing
"CurC")) # constant current catch levelsPlot projected SSB and exploitation relative to MSY levels:
Splot(myMSE)
Status at the end of the projection:
DFO_proj(myMSE)
Tabulation of DFO-related probabilities:
DFO_tab(myMSE)[1:3,] # ?DFO_tab to find out more
#> MP P_Cr_S P_Ct_S P_H_S POF_S STY P_Cr_L P_Ct_L P_H_L POF_L LTY AAVY
#> 2 FMSYref50 0 1 99 0 103 0 5 95 0 94 7
#> 3 CurC 0 0 100 0 31 0 0 100 0 31 0
#> 1 MP.MC 0 0 100 0 23 0 0 100 0 30 3
#> P_Reb
#> 2 91
#> 3 100
#> 1 1007 Specifying Operating Models
7.1 Objects and data
The Inverts package includes all of the generic functions for searching for objects, interrogating operating models, running MSEs and plotting results.
The package comes with reference operating models (class OM) already conditioned for each species (e.g. OM_GD_7 is the operating model for Geoduck in management area 7).
You can search the loaded Inverts packages for any object of a given class using the objs() function:
library(Inverts)
objs('OM') # lists all objects of class 'OM' in the loaded Inverts packages
#> Searching for objects of class OM in package: Inverts
#> Searching for objects of class OM in package: Inverts.GD
#> Searching for objects of class OM in package: Inverts.MC
#> [1] "OM.GD.7" "OM.GSU.12" "OM.MC.E" "OM.RSC.7F"When you load a species-specific library, you get bespoke input files (class ‘In’) and functions specific to that species:
library(Inverts.GD) # load the geoduck library
objs('In') # lists all objects of class 'In' in the loaded Inverts packages
#> Searching for objects of class In in package: Inverts
#> Searching for objects of class In in package: Inverts.GD
#> Searching for objects of class In in package: Inverts.MC
#> [1] "In.GD.1" "In.GD.106" "In.GD.12" "In.GD.13" "In.GD.14" "In.GD.15"
#> [7] "In.GD.16" "In.GD.17" "In.GD.18" "In.GD.19" "In.GD.2" "In.GD.23"
#> [13] "In.GD.24" "In.GD.25" "In.GD.26" "In.GD.27" "In.GD.3" "In.GD.4"
#> [19] "In.GD.5" "In.GD.6" "In.GD.7" "In.GD.8" "In.MC.C" "In.MC.D"
#> [25] "In.MC.E" "In.MC.F"Input object are specific to a particular region. You can plot the areas and some details of the operating models using the about() function:
about(In.GD.7) # Plot some information about the input object
This also works for operating model objects which may have more information if they have been conditioned on data (such as unfished recruitment R0 (millions) and stock depletion (recent SSB relative to unfished):
about(OM.MC.E) # Plot some information about an operating model 
If you don’t want the figure and just want the data you can alter some arguments of the about function
about(OM.GD.7, plot.fig = F, ret.dat = T) # Extract only the data
#> $Name
#> [1] "Geoduck Management Area 7"
#>
#> $Years
#> [1] 1976 2023
#>
#> $M
#> 5% 95%
#> 0.042 0.070
#>
#> $K
#> 5% 95%
#> 0.24 0.26
#>
#> $Linf
#> 5% 95%
#> 137.3 138.0
#>
#> $R0
#> 5% 95%
#> 2.69 4.59
#>
#> $Depln.
#> 5% 95%
#> 0.75 0.85Remember all of the functions are documented and can be accessed with the ? operator:
?aboutThe input objects are the smallest and most concise basis for creating operating models and MSE analyses in a reproducible way. These objects contain the parameters and data required to condition operating models using the Rapid Conditioning Model (RCM) of OpenMSE.
names(In.GD.7) # slot 1 is the name, slot 2 is the OM parameters, slot 3 contains the data
class(In.GD.7$OM) # this object is of class 'OM' for operating model, containing all parameters for the OM
class?OM # find out more about this class
class(In.GD.7$RCMdata) # this object is of class 'RCMdata', containing all fishery data used for conditioning the OM
class?RCMdata # find out more about this class
slotNames(In.GD.7$RCMdata) # list all the data types included in the object
In.GD.7$RCMdata@Chist # the historical catches7.2 Conditioning operating models
Using these input objects you can condition an operating model for this stock unit. The reason this is done at this level (we don’t provide pre-fitted OMs, other than the demonstration reference operating models in the Inverts package) is because you require control over the number of simulations to run. For example, in this script or a workshop we might wish to demonstrate functionality quickly using only a few simulations. In a setting where management advice is to be followed you might want 200+ simulations.
Each input object has up to 500 simulations already specified.
You can select a number of simulations and then condition an operating model. In this case we are using a function from the Geoduck package ‘Inverts.GD’. This conditioning function is bespoke and include features for model specification, fitting and acceptance/rejection that are specific to Geoduck. Other species have similar functions (e.g. cond.MC() for manila clam).
This code might take a few minutes to run. To speed up this calculation and all of the other parts of this tutorial we are going to specify only 24 simulations:
?cond.MC()
myfit = cond.MC(In.MC.E, sims=24)Once you have a fitted operating model the first thing you might wish to do is visualize the data, fit to data and model estimates:
plot(myfit)Alternatively, you could focus on the operating model characteristics including those parts of the operating model that simulate data and implement advice:
plot(myfit@OM)8 Calculating reference points
8.1 The Inverts Package
The Inverts package contains functions for the extraction and plotting of management reference points.
Firstly, you need to reconstruct the historical fishery in your R session using the runMSE() function with argument Hist=T. This means that the historical spool-up period of the operating model (before you were do to management projections) is calculated including a very wide range of reference points and quantities relating to reference points.
In later steps, projections will be undertaken to test management options (data quality, management procedures etc). By pre-calculating the historical stock and reference points, you don’t have to do this every time you do a projection (historical dynamics are fixed for each operating model) which makes simulation testing faster.
Hist = runMSE(myfit@OM, Hist=T)You can obtain a tabulated summary of reference points for a Hist object using the Ref.Points function:
RPs = Ref.Points(Hist)
names(RPs) # names of listed reference points
#> [1] "SSB_SSB0" "SSB_SSBMSY" "F_FMSY"
head(RPs$SSB_SSB0) # start of the biomass relative to unfished entry
#> Mean X5. X25. X50. X75. X95. P20 P40
#> 2000 0.4105 0.2829 0.3677 0.4202 0.4542 0.5285 0 0.4062
#> 2001 0.4475 0.3499 0.4052 0.4556 0.4795 0.5405 0 0.2188
#> 2002 0.4876 0.4283 0.4639 0.4946 0.5062 0.5474 0 0.0312
#> 2003 0.5123 0.4612 0.4966 0.5185 0.5300 0.5482 0 0.0000
#> 2004 0.5487 0.4890 0.5357 0.5528 0.5679 0.5913 0 0.0000
#> 2005 0.5531 0.4745 0.5319 0.5600 0.5829 0.6125 0 0.0000You can also use the internal plotting function on any of these tables if you want a single plot:
RefPlot(RPs[[3]],names(RPs)[3])
For a summary of stock status in the most recent historical year two functions are provided:
Status(Hist) # provides all the estimated values of SSB and F relative to reference points 
#> $SSB_SSB0
#> [1] 0.4610 0.4896 0.6213 0.5724 0.7093 0.5927 0.5080 0.5980 0.5586 0.5436
#> [11] 0.6989 0.5919 0.5011 0.4337 0.7043 0.6705 0.5918 0.6463 0.5365 0.6874
#> [21] 0.5946 0.5952 0.6089 0.4805 0.5515 0.5904 0.4890 0.4778 0.6299 0.7176
#> [31] 0.6209 0.6538
#>
#> $SSB_SSBMSY
#> [1] 1.2762 1.2438 1.5581 1.5906 1.5521 1.6405 1.2987 1.7274 1.4393 1.4272
#> [11] 1.4731 1.8246 1.3928 1.1431 1.5741 1.6345 1.5903 1.3602 1.4067 1.4718
#> [21] 1.5815 1.4919 1.6031 1.1908 1.3907 1.4854 1.2399 1.2056 1.6559 1.6081
#> [31] 1.5366 1.6894
#>
#> $F_F01
#> [1] 0.3575 0.2505 0.1381 0.2268 0.0649 0.1634 0.2194 0.1587 0.1796 0.1931
#> [11] 0.0831 0.1666 0.2762 0.3331 0.0771 0.1046 0.1457 0.1321 0.2389 0.0800
#> [21] 0.2065 0.1275 0.1467 0.2522 0.1963 0.1395 0.2543 0.2591 0.1276 0.0715
#> [31] 0.1568 0.1234
#>
#> $F_FMSY
#> [1] 0.1596 0.1241 0.0404 0.0751 0.0265 0.0421 0.1011 0.0425 0.0770 0.0898
#> [11] 0.0360 0.0507 0.0863 0.1719 0.0251 0.0281 0.0523 0.0428 0.1102 0.0308
#> [21] 0.0487 0.0523 0.0353 0.1437 0.0932 0.0747 0.1344 0.1815 0.0318 0.0231
#> [31] 0.0400 0.0306DFO_status(Hist) # plot estimated values of SSB and F relative to reference points on standard DFO regions 
If you would rather access all the raw data yourself, these are available in the Hist object in the following slots:
names(Hist@TSdata) # Time series data like F (Find) and spawning biomass (SSB)
#> [1] "Number" "Biomass" "VBiomass"
#> [4] "SBiomass" "Removals" "Landings"
#> [7] "Discards" "Find" "RecDev"
#> [10] "SPR" "Unfished_Equilibrium"
names(Hist@Ref$ReferencePoints) # Equilibrium reference points (parameters for most recent historical year)
#> [1] "N0" "B0" "SSB0" "SSN0" "VB0"
#> [6] "MSY" "FMSY" "SSBMSY" "BMSY" "VBMSY"
#> [11] "UMSY" "FMSY_M" "SSBMSY_SSB0" "BMSY_B0" "VBMSY_VB0"
#> [16] "RefY" "MGT" "Blow"
names(Hist@Ref$ByYear) # Equilibrium reference points calculated from parameters in each historical year
#> [1] "N0" "SN0" "B0" "SSB0" "VB0" "R0"
#> [7] "h" "MSY" "FMSY" "SSBMSY" "BMSY" "VBMSY"
#> [13] "F01_YPR" "Fmax_YPR" "F_SPR" "Fcrash" "Fmed" "SPRcrash"
names(Hist@Ref$Dynamic_Unfished) # Dynamic reference points (assuming no fishing for all years)
#> [1] "N0" "B0" "SN0" "SSB0" "VB0" "Rec"It is possible to compare biological status across various operating models. Lets fit a second one and do the historical spool up:
myfit2 = cond.MC(In.MC.C)
Hist2 = runMSE(myfit2@OM)
Now we can compare the current biomass estimate among the operating models:
DFO_status_comp(list(Hist,Hist2))
#> P(Crit.) P(Caut.) P(Healthy)
#> M_CCME 0 0 1
#> M_CCMC 0 0 18.2 The Reference Point Calculator (RPC) Package
Reference Point Calculator (Huynh et al. 2022) is an R package and App designed for the exploration of reference points and the testing of management procedures / Harvest Control Rules (HCRs). The package was developed as part of a DFO grant and allows for the calculation of a very wide range of reference points.
Many of the RPC plotting functions use objects of class ‘Hist’ that are historical reconstructions of the type calculated above (using runMSE(OM, Hist=T)).
If you want a summary of the historical stock dynamics:
hist_bio(Hist)
Also available are the yield curves (equilibrium catches at various population / exploitation levels):
# hist_YieldCurve(Hist)Here we plot the SSB history relative to MSY:
LRP_SSBMSY(Hist)
Exploitation rate relative to FMSY:
LRP_FMSY(Hist)
Spawning Potential Ratio:
LRP_SPR(Hist)
The historical recruitment:
LRP_R(Hist)
Spawning stock biomass (SSB) relative to the (mean) SSB corresponding to 50% max recruitment:
LRP_50Rmax(Hist)
Historical surplus production:
LRP_SP(Hist)
8.3 The Reference Point Calculator online App
The operating model is fully compatible with the shiny app Reference Point Calculator to provide a very comprehensive report on ALL fishery reference points ever described (including the various plots above) and also allows for the testing of management procedures.
To use the app you should save the operating model object to your computer, open the app link above and then upload your saved operating model object:
saveRDS(myfit@OM, file="C:/temp/my_fitted_OM.rds") # this is the file you upload to the RPC app
9 Testing management procedures
9.1 Running simulations
Each species-specific R package has a dedicated ‘status quo’ management procedure. These can be easily found using the list objects function objs().
objs('MP') # invertebrates MPs
#> Searching for objects of class MP in package: Inverts
#> Searching for objects of class MP in package: Inverts.GD
#> Searching for objects of class MP in package: Inverts.MC
#> [1] "MP.GD" "MP.MC"Many other MPs that are available for testing (from the openMSE libraries) can be found using the avail() function:
avail('MP') # all MPs that have been loaded included those with openMSE pacakages
#> → Searching for objects of class "MP" in package "MSEtool"
#> → Searching for objects of class "MP" in package "SAMtool"
#> → Searching for objects of class "MP" in package "DLMtool"
#> [1] "curEref" "Emp" "FMSYref" "FMSYref50" "FMSYref75"
#> [6] "NFref" "DDSS_4010" "DDSS_75MSY" "DDSS_MSY" "SCA_4010"
#> [11] "SCA_75MSY" "SCA_MSY" "SP_4010" "SP_75MSY" "SP_MSY"
#> [16] "SSS_4010" "SSS_75MSY" "SSS_MSY" "AvC" "AvC_MLL"
#> [21] "BK" "BK_CC" "BK_ML" "CC1" "CC2"
#> [26] "CC3" "CC4" "CC5" "CompSRA" "CompSRA4010"
#> [31] "CurC" "curE" "curE75" "DAAC" "DBSRA"
#> [36] "DBSRA_40" "DBSRA4010" "DCAC" "DCAC_40" "DCAC_ML"
#> [41] "DCAC4010" "DCACs" "DD" "DD4010" "DDe"
#> [46] "DDe75" "DDes" "DepF" "DTe40" "DTe50"
#> [51] "DynF" "EtargetLopt" "Fadapt" "Fdem" "Fdem_CC"
#> [56] "Fdem_ML" "Fratio" "Fratio_CC" "Fratio_ML" "Fratio4010"
#> [61] "GB_CC" "GB_slope" "GB_target" "Gcontrol" "HDAAC"
#> [66] "ICI" "ICI2" "Iratio" "Islope1" "Islope2"
#> [71] "Islope3" "Islope4" "IT10" "IT5" "Itarget1"
#> [76] "Itarget1_MPA" "Itarget2" "Itarget3" "Itarget4" "ItargetE1"
#> [81] "ItargetE2" "ItargetE3" "ItargetE4" "ITe10" "ITe5"
#> [86] "ITM" "L95target" "LBSPR" "LBSPR_MLL" "Lratio_BHI"
#> [91] "Lratio_BHI2" "Lratio_BHI3" "LstepCC1" "LstepCC2" "LstepCC3"
#> [96] "LstepCC4" "LstepCE1" "LstepCE2" "Ltarget1" "Ltarget2"
#> [101] "Ltarget3" "Ltarget4" "LtargetE1" "LtargetE4" "matlenlim"
#> [106] "matlenlim2" "MCD" "MCD4010" "minlenLopt1" "MRnoreal"
#> [111] "MRreal" "Rcontrol" "Rcontrol2" "SBT1" "SBT2"
#> [116] "slotlim" "SPmod" "SPslope" "SPSRA" "SPSRA_ML"
#> [121] "YPR" "YPR_CC" "YPR_ML"Most of these are reference MPs (e.g., curE - fishing at current effort levels, FMSYref75 - fishing at 75% of perfectly known FMSY, NFref - no fishing) or ‘data limited’ MPs (e.g. DBSRA, DCAC, AvC, BK_ML). There are however a handful of data rich MPs that are full stock assessments so that it is possible to quantify the benefit of moving towards more analysis-heavy types of management (e.g., SP_MSY, DDSS_MSY, SCA_MSY etc.).
The invertebrate-specific MPs (e.g. MP.MC) are designed to closely follow existing management measures but have sufficient flexibility to accommodate a wide range of modifications so that management alternatives can be investigated (more on that below in the ‘Specifying management options’ section)
Running MSE projections uses the Project() function, where the historical simulation object (we made Hist earlier) is projected for at least one management scenario. Note that each simulation of this projection is identical in terms of recruitment deviation and observation error - the only difference is the MP, so simulation 1 for the first MP tested can be compared directly with simulation 1 for the other MPs. By sharing the same projected conditions, the relative differences among MPs are revealed more efficiently.
Here we do a projection of our status quo management MP for manila clam and and reference MP that fishes at half of perfectly know FMSY levels:
myMSE = Project(Hist, c('MP.MC', # status quo manila clam management
"FMSYref50", # perfect implementation of 50% FMSY fishing
"CurC")) # constant current catch levelsThe object that was created (myMSE) is of class ‘MSE’ and contains a wide range of information about the operating model, reference points and the projected simulations:
class?MSE9.2 Performance metrics
Since the MSE object contains a comprehensive record of everything that was simulated in the projection, we can produce a large number of possible performance metrics. There are functions that provide outputs specific to particular management authorities. For example, there are DFO-specific plots that reference the probabilty of being in the critical, cautious and healthy zones over the short and long time periods:
DFO_tab(myMSE)[1:3,]
#> MP P_Cr_S P_Ct_S P_H_S POF_S STY P_Cr_L P_Ct_L P_H_L POF_L LTY AAVY
#> 2 FMSYref50 0 1 99 0 103 0 5 95 0 94 7
#> 3 CurC 0 0 100 0 31 0 0 100 0 31 0
#> 1 MP.MC 0 0 100 0 23 0 0 100 0 30 3
#> P_Reb
#> 2 91
#> 3 100
#> 1 100
?DFO_tab # to find out more
#> starting httpd help server ... doneOther outputs include NOAA_plot that shows the probability of not overfishing, probability of being above 50% BMSY (P50), long term yield and the probability catch changes are less than 15% among management updates:
NOAA_plot(myMSE)
#> PNOF B50 LTY VY
#> MP.MC 100 100 5.0 100
#> FMSYref50 100 100 100.0 100
#> CurC 100 100 6.2 1009.3 Plotting results
You can visualize some of the MSE outcomes using built-in openMSE plotting functions (basic but fast):
Splot(myMSE)
There are also a set of DFO figures and diagnostics:
DFO_proj(myMSE)
There are also a set of tradeoff plots that show median yield and biomass outcomes:
Cplot(myMSE)
#> Warning in MPtype(MSEobj@MPs): Some MPs are not found in environment: MP.MC
You do have access to all the data in the MSE object so can plot your own custom plots easily enough:
Year = myMSE@OM$CurrentYr[1] + 1:myMSE@proyears
matplot(Year,t(myMSE@SB_SBMSY[,1,]),type="l",lty=1, ylab="Spawning biomass relative to BMSY"); grid()
abline(h=1,col="grey",lwd=2,lty=2)
9.4 Specifying management options
For example, the MP.MC function (default management procedure for manila clam) has the following arguments:
args(MP.MC)
#> function (x, Data, reps = 1, Min.size = 35, Max.size = NaN, CEff.Mult = 1,
#> C_I.targ = 1, I.targ = 0.5, IS.targ = 0, IS.yrs = 6, IS.fac = 1,
#> TAC.calc = NaN, maxTAC = 5, minTAC = 0.1, TACdec = 0.2, TACinc = 0.1,
#> I.enp = 0.25, I_freq = c(0, 1, 0), calib_yrs = 2, HCR_CP_B = c(0,
#> 0), HCR_CP_TAC = c(0, 1), curI_2_target = 2, DR = 0,
#> Fdisc = 0.5)
#> NULLYou can find out what each of these arguments do via the in-line help:
?MC.MPThe next section describes how these generic MPs can be modified to test a wide range of alternative management options.
10 Designing new management procedures
10.1 Status quo and reference management procedures
Before testing alternative management options it is beneficial to understand the range of outcomes that might be possible. Do do this it is advisable to run projections with reference MPs. Here we will evaluate ‘no fishing’, ‘FMSY fishing’, and status quo fishing.
refMSE = Project(Hist, c('MP.MC', # status quo manila clam management
"FMSYref", # perfect implementation of FMSY fishing
"NFref")) # no fishingProjection plots let us understand how the FMSY and no fishing scenarios bracket our existing status quo management.
Splot(refMSE)
10.2 Implementing size limits
By default, the manila clam MP imposes a minimum size limit of 35mm and keeps effort at current levels. If you wanted to test an alternative minimum size limit you would alter that particular argument:
SL25 = MP.MC # copy the status quo management procedure
formals(SL25)$Min.size = 25 # set the size limit to 25mm
class(SL25) = 'MP' # assign your new function the class 'MP'Lets test it out:
SL_MSE = Project(Hist, c("MP.MC","SL25"))
Splot(SL_MSE)
Because both MPs are fishing much below FMSY the only impact of the size limit change is on the yield. Lets try this again with 25 adn 35mm size limits but with 20x the amount of fishing pressure as today:
SL25_E20 = SL25 # copy the basic operating model
formals(SL25_E20)$CEff.Mult = 20 # set the size limit to 25mm
class(SL25_E20) = 'MP' # assign your new function the class 'MP'
SL35_E20 = MP.MC
formals(SL35_E20)$CEff.Mult = 20
SL_MSE_E20 = Project(Hist, c("SL25_E20","SL35_E20"))
Splot(SL_MSE_E20)
Now we can see closer to FMSY fishing pressure, more comparable yield outcomes but somwehatever higher biomass outcomes for the higher size limit.
10.3 Implementing effort limits
Lets say you wanted to test a set of levels for a particular management measure. Here we are going to investigate various levels of current fishing pressure (this could be managed by access limits, seasonal spatial closures, license limits etc). You can see from the code above that this would be cumbersome for more than a few MPs. To avoid the mess, the Inverts package comes with a function for spawning multiple MPs from differing values of a particular argument:
Effs = c(1,2,4,8,16)
E_MPs = paste0("E",Effs)
spawnMP("MP.MC", spawnlist = list(CEff.Mult = Effs), spawnnames = E_MPs)
Based on our status quo management procedure, in just three lines of code we have successfully made a set of 5 MPs (with names E1, E2, E4, E8, E16) that fish at between 1x and 16x current effort levels.
We can now test those:
MSE_Eff = Project(Hist, E_MPs)
Splot(MSE_Eff)
These plots show the fundamental trade-off between what is taken (yield) and what is left over in the water (biomass). We can created custom trade-off plot to show this more clearly:
Exp_longterm_B_BMSY = apply(MSE_Eff@SB_SBMSY[,,30],2,mean)
Exp_longterm_Yield = apply(MSE_Eff@Catch[,,30],2,mean)
plot(Exp_longterm_B_BMSY,Exp_longterm_Yield, col="white");grid()
text(Exp_longterm_B_BMSY,Exp_longterm_Yield, MSE_Eff@MPs)
10.4 Index ratio approaches
By default, the invertebrate MPs fish at current fishing effort levels without total allowable catch (TAC) control. You may however wish to invetigate options that explicity set catch limits. The simplest way to do this are empirical management procedures that use data directly and do not require an estimation / assessment step. These are increasingly applied in the management of a very wide range of fish species and their adoption is the principle objective of MSE in most fishery settings.
The most commonly applied type of empirical management procedure aims to fish at a stable exploitation rate by setting catch limits to a fixed multiple of a relative abundance index:
Constant exploitation rates already include a degree of precaution since equilibrium sustainable exploitation rates increase exponentially as stocks decline. For example, in most cases it is necessary to fish at levels larger than two or three times FMSY to sustainably deplete a stock below half BMSY:
Nevertheless, in Canada, fisheries management often seeks to impose an additional ‘hockey stick’ harvest control rule (an algorithm that throttles exploitation rate as stock biomass declines - the green line of the plot above). Hockey stick harvest control rules require control points for the maximum exploitation rate (y axis) and the biomass (x axis) below which exploitation rate is reduced and where exploitation rate is zero.
In the figure above the harvest control rule for which Exploitation rate is a function of biomass (left panel), can be rephrased as an empirical rule based on observed catches and an observed Index of abundance (right panel). The control points are then phrased in terms of the recommended catch per index (exploitation rate) and the index levels below which catches are zero or catch per index is reduced.
The status quo MPs have the option of specifying control using such a system. The arguments for this function are explained in the help documentation:
objs('MP') # what inverts-specific MPs are available?
#> Searching for objects of class MP in package: Inverts
#> Searching for objects of class MP in package: Inverts.GD
#> Searching for objects of class MP in package: Inverts.MC
#> [1] "MP.GD" "MP.MC"
?MP.MC # help documentation for MP.GD
args(MP.MC)# a more concise list of arguments to the function
#> function (x, Data, reps = 1, Min.size = 35, Max.size = NaN, CEff.Mult = 1,
#> C_I.targ = 1, I.targ = 0.5, IS.targ = 0, IS.yrs = 6, IS.fac = 1,
#> TAC.calc = NaN, maxTAC = 5, minTAC = 0.1, TACdec = 0.2, TACinc = 0.1,
#> I.enp = 0.25, I_freq = c(0, 1, 0), calib_yrs = 2, HCR_CP_B = c(0,
#> 0), HCR_CP_TAC = c(0, 1), curI_2_target = 2, DR = 0,
#> Fdisc = 0.5)
#> NULLTo specify the index ratio approach you set the argument TAC.calc equal to ‘Ratio’. There are then four arguments that shape the hockey stick control rule:
- C_I.targ - the maximum catch / index (maximum exploitation rate)
- curI_2_target - a calibration factor that states where the most recent index observation (e.g., 2024) is in relation to a target level.
- HCR_CP_B - the position of the lower and upper index control points relative to target (e.g. 0.5 and 1)
- (optionally) HCR_CP_TAC - position of the lower and upper catch per index control points (default to 0 and 1 as seen in the figure above)
myIR = MP.MC
formals(myIR)$TAC.calc = "Ratio" # Empirical index-ratio approach (the default is 'NaN' - no TAC limits, only effort)
formals(myIR)$C_I.targ = 1 # Aiming to fish at the current catch per index level
formals(myIR)$curI_2target = 2 # The stock (as indicated by the index) is currently twice a calibration level (lets say BMSY)
formals(myIR)$HCR_CP_B = c(0.5,1)# Control points are 0.5x and 1x BMSY (calibrated according to curI_2target above)
formals(myIR)$HCR_CP_TAC = c(0,1)# Zero fishing below lower control point, current catch/index above upper control point
class(myIR) = "MP" # make sure it is of OpenMSE class 'MP'If you wanted a constant exploitation rate policy you would either specify both biomass control points at the origin:
formals(myIR)$HCR_CP_B = c(0,0) # Lower and upper biomass control points are at the origin = constant catch per index.
or make both exploitation rate control points at the same level:
formals(myIR)$HCR_CP_TAC = c(1,1)# Lower and upper exploitation rate control points are the same = constant catch per index
These constant exploitation rate policies are worth investigating particularly for short-lived species where natural processes contribute substantially to variability in stock levels, lower constant exploitation rate policies (essentially no HCR) often out-perform hockey stick control rules that have much higher catch variability for insubstantial gains in yields. Constant exploitation rate policies are a feature of management procedures adopted for a range of pelagic fish species ranging from Peruvian anchoveta to Atlantic bluefin tuna.
10.5 Index target approaches
An alternative to empirical rules aiming for constant exploitation rate are MP that aim for a particular index level. These may be more appropriate in situations where there is variability in the sustainable rate of exploitation (due to changes in, for example, somatic growth, natural mortality rate, fishery size or age selectivity).
The index target approaches make changes to the catch advice to reach the target index level. These rules only ask for a target index level ‘I.targ’:
myIT = MP.MC
formals(myIT)$TAC.calc = "Target" # Empirical index-ratio approach
formals(myIT)$I.targ = 1360 # Target index level
class(myIT) = "MP"
10.6 Index slope approaches
In cases where you wish to maintain a trajectory in stock biomass at some level, an alternative to index target MPs are index slope MPs. Similarly to the index target approaches you must specify a target slope (‘IS.targ’) and two other arguments that control the calculation of the slope (‘IS.yrs’) and the responsiveness of the MP (‘IS.fac’):
myIS = MP.MC
formals(myIS)$TAC.calc = "Slope" # Empirical index-ratio approach
formals(myIS)$IS.targ = 0 # Stable stock levels
formals(myIS)$IS.fac = 1 # TAC changes are proportional the slope change (e.g. a 10% reduction in TAC for a 10% decline in index slope below target levels)
class(myIS) = "MP"10.7 Stock assessments
Management procedures can include estimation models that aim to mimick the real stock assessment processes. These are sometimes referred to as ‘model-based MPs’.
Even in situations where stock assessments are not possible for data, logistical or management reasons, it may be desirable to evaluate their potential performance (perhaps as a yardstick for current approaches). In MSE projections, all data are simulated including types that may not be currently available. That allows for the theoretical testing of alternative more data intensive options.
The OpenMSE libraries come with a range of prespecified data-rich assessment MPs that include state-space surplus production (SP), state-space delay differential (DDSS) and statistical catch at age (SCA) models. A few derivatives are available that fish at constant FMSY (’_MSY’), constant 75% of FMSY (’_75MSY’) and use a 40-10 HCR (’_4010’, FMSY above 40% B0, zero fishing below 10% B0):
avail('MP')
?SCA_MSY
?DDSS_75MSY
?SP_4010You can however design your own data rich MPs by combining assessments (class ‘Assess’) with HCRs (class ‘HCR’):
avail('Assess')
avail('HCR')
SP_Fox_40_10 = make_MP(SP_Fox,HCR40_10) # combination of the Fox surplus production model and the 40_10 HCRYou can also make your own MPs with HCRs that have custom control points:
myMP_HCR = make_MP(SCA, HCR_segment, # HCR segment has arguments (below) for the control points
OCP_type = "SSB_SSBMSY", # independent variable (x axis) is spawning stock biomass relatie to MSY levels
Ftarget_type = "FMSY", # dependent variable (y axis) of HCR is F relative to FMSY
OCP = c(0.1,0.4), # lower and upper control points are 10 and 40 percent of SSB
relF = c(0,1))
class(myMP_HCR)10.8 Short-cuts
Running model-based MPs can be computationally intensive. An alternative is to take the control points of an HCR not from an estimation model but from the operating model directly assuming some realistic level of observation (data going into the assessment) and estimation error (error in estimates of control points of an HCR in this context). Because there is no estimation, short-cuts can provide indicative results in a small fraction of the time as the full model-based MP.
my_SC_4010 = make_MP(Shortcut, HCR40_10, # HCR segment has arguments (below) for the control points
method = "N", # stock numbers enter the HCR
N_err = c(0.1, # 10% coefficient of variation (precision)
0.2, # 20% lag-1 autocorrelation in error
1)) # Unbiased on average (1.1 would be 10% positively biased)
class(my_SC_4010)10.9 Idealized management
Short cut methods provide a way to create idealized ‘perfect information’ management (i.e. setting the errors to very low levels:
tiny = 1E5
perfect_4010 = make_MP(Shortcut, HCR40_10, # HCR segment has arguments (below) for the control points
method = "N", # stock numbers enter the HCR
N_err = c(tiny, # highly precise
tiny, # no lag
1)) # UnbiasedThere are also some specified perfect information MPs that fish at a fraction of FMSY:
- FMSYref, FMSYref75, FMSYref50
10.10 Additional TAC control parameters
When using the specified MPs such as There are a number of other arguments that apply to all methods of catch advice (TAC.calc is set to ‘Slope’,‘Target’ or ‘Rate’) that provide additional control over TAC changes. These include:
- max TAC - that highest TAC level that is recommended by the MP
- minTAC - the lowest TAC level that is recommended by the MP
- TACdec - the highest fractional TAC decline that the MP can make (e.g. 0.2 means that TAC(y) must be higher than 0.8 x TAC(y-1))
- TACinc - the highest fraction TAC increase that the MP can make (e.g. 0.3 mean that TAC(y) must be lower than 1.3 * TAC(y-1))
- I.enp - Indices are processed via a polynomial smoother to try to remove noise from signal. The degree of smoothing is controlled by this parameter. ENP stands for effective number of parameters per datum. So higher ENPs are more closely fitting (less smooth) polynomials. The default of 0.25 means that for every 4 data points the polynomial adds another parameter for smoothing.
- I_freq - the frequency that the indices will be observed in the forward projection. Let us assume there are three indices available a vector I_freq = c(1,3,0) would mean that there would be an observation every year for index 1, a new observation every three years for index 2, and no new data for index 3.
- calib_yrs - the number of years used to define ‘recent historical levels’. When arguments such as curI_2_target are used, they determine ‘current’ based on the calib_yrs number of years in the historical time period. So for example if the OMs were fitted to 2023 and calib_yrs was set to 3, the MP would use years 2021-2023 to determine ‘recent’ index and catch levels.
10.11 Other management levers
The default MPs (e.g. GC.MP, MC.MP, GSU.MP, RSU.MP) assume that current fishing effort (days per year) will be the same in the future. You can explore alternative effort scenarios with the ‘Effort’ argument.
A default discarding rate and post-release mortality rate are also assumed. The shift to a new gear might be expected to change both of these. These are easily changed:
new_gear_MP = RSU.MP
formals(new_gear_MP)$Effort = 0.9 # a 10% decline in effort is expected
formals(new_gear_MP)$DR = 0.15 # 15% of the catch (numbers) is expected to be discarded
formals(new_gear_MP)$Fdisc = 0.8 # only 20% of discards are expected to survive
class(new_gear_MP) = "MP"
11 Species-specific considerations
Some species are managed in a unique way that requires some specification of operating model and management procedure features. Those are described here.
11.1 Rotational closures
In some areas, species such as Geoduck and Cucumber are managed using rotational closures. This means that management procedures for these species are set up to try and identify the recent cycle of closures and maintain this in future years. These MPs (MP.GD an MP.RSC) can also be specified with a user defined pattern of closures using the ‘rotation’ argument:
?MP.GD
myGD = MP.GD
formals(myGD)$rotation = rep(c(1,1,0),100)[1:OM.GD.7@proyears] # twice open for every year closed for the projected years (starting, open, open, closed)A subtlety of geoduck and red sea cucumber operating models is that the observation error model is disrupted by the on/off nature of historical catch observations. You need to specify these on the Hist object before you do projections, for example:
Hist = runMSE(OM.GD.7,Hist=T)
newHist = Obs.GD(Hist, Catch_CV = 0.02, Catch_bias = 1) # catch observations are unbiased with a 2% CV
testy = Project(newHist,"MP.GD")12 Modifying and reconditioning operating models
12.1 Input object configuration
Input objects contain all of the parameters and data required to condition an operating model using the species-specific conditioning function (e.g., cond.MC, cond.GD, cond.GSU, cond.RSC).
You can list the available input objects using the objs() function:
objs('In') # all the invertebrate input objects available in this session
#> Searching for objects of class In in package: Inverts
#> Searching for objects of class In in package: Inverts.GD
#> Searching for objects of class In in package: Inverts.MC
#> [1] "In.GD.1" "In.GD.106" "In.GD.12" "In.GD.13" "In.GD.14" "In.GD.15"
#> [7] "In.GD.16" "In.GD.17" "In.GD.18" "In.GD.19" "In.GD.2" "In.GD.23"
#> [13] "In.GD.24" "In.GD.25" "In.GD.26" "In.GD.27" "In.GD.3" "In.GD.4"
#> [19] "In.GD.5" "In.GD.6" "In.GD.7" "In.GD.8" "In.MC.C" "In.MC.D"
#> [25] "In.MC.E" "In.MC.F"
about(In.MC.E) # information about the location and OM for manila clam in clam management area E
Input objects are lists that include information about the geographic location of the stock (in this case CMA), the operating model (OM) and the available data (RCMdata):
class(In.MC.E)
#> [1] "In"
names(In.MC.E)
#> [1] "CMA" "OM" "RCMdata"We can make changes to model parameters (OM) or data (RCMdata) by manipulating these input objects.
12.2 Changing parameters
Operating models are built from the parameter ranges specified in the operating model object and a ‘custom parameter’ slot in the operating model.
myOM = In.MC.E$OM
slotNames(myOM) # most of these are parameter ranges for uniform sampling of values
#> [1] "Name" "Agency" "Region" "Sponsor"
#> [5] "Latitude" "Longitude" "nsim" "proyears"
#> [9] "interval" "pstar" "maxF" "reps"
#> [13] "cpars" "seed" "Source" "Common_Name"
#> [17] "Species" "maxage" "R0" "M"
#> [21] "Msd" "h" "SRrel" "Perr"
#> [25] "AC" "Linf" "Linfsd" "K"
#> [29] "Ksd" "t0" "LenCV" "L50"
#> [33] "L50_95" "D" "a" "b"
#> [37] "Size_area_1" "Frac_area_1" "Prob_staying" "Fdisc"
#> [41] "nyears" "CurrentYr" "EffYears" "EffLower"
#> [45] "EffUpper" "Esd" "qinc" "qcv"
#> [49] "L5" "LFS" "Vmaxlen" "isRel"
#> [53] "LR5" "LFR" "Rmaxlen" "DR"
#> [57] "Spat_targ" "MPA" "Misc" "Cobs"
#> [61] "Cbiascv" "CAA_nsamp" "CAA_ESS" "CAL_nsamp"
#> [65] "CAL_ESS" "Iobs" "Btobs" "Btbiascv"
#> [69] "beta" "LenMbiascv" "Mbiascv" "Kbiascv"
#> [73] "t0biascv" "Linfbiascv" "LFCbiascv" "LFSbiascv"
#> [77] "FMSY_Mbiascv" "BMSY_B0biascv" "Irefbiascv" "Brefbiascv"
#> [81] "Crefbiascv" "Dbiascv" "Dobs" "hbiascv"
#> [85] "Recbiascv" "sigmaRbiascv" "Eobs" "Ebiascv"
#> [89] "TACFrac" "TACSD" "TAEFrac" "TAESD"
#> [93] "SizeLimFrac" "SizeLimSD"
names(myOM@cpars) # these are higher level, per simulation samples specifed by the user
#> [1] "SLarray" "CAL_binsmid" "L50" "K" "Linf"
#> [6] "Wa" "M"Operating model building first looks for the bespoke inputs of the custom parameters before settling on the parameter ranges. For example the operating models will use OM@cpars$M before they use OM@M:
validcpars() # list all the custom parameter optionslength(myOM@cpars$M) # nsim samples of the instantaneous natural mortality rate
#> [1] 500
dim(myOM@cpars$SLarray) # fishery selectivity by simulation (500), length class (45) and year (54)
#> [1] 500 45 54If you wanted to create the same OM but with an M value that is 20% lower you would do this:
myOM = In.MC.E$OM
prevM = In.MC.E$OM@cpars$M
myOM@cpars$M = prevM * 0.8If you wanted to sample your own M distribution and put it in the cpars slot you could do this:
nsim = myOM@nsim
myOM@cpars$M = rlnorm(nsim,log(0.31),0.1) # a lognormal distribution with median 0.31 and sd 0.1
hist(myOM@cpars$M )
12.3 Adding data / reconditioning operating models
The data used in conditing are stored in the slot RCMdata (rapid conditioning model data object). These data are of class RCMdata and have a defined format.
mydata = In.MC.E$RCMdata
slotNames(mydata)
#> [1] "Chist" "C_sd" "C_wt" "Ehist" "CAA"
#> [6] "CAA_ESS" "CAL" "CAL_ESS" "length_bin" "MS"
#> [11] "MS_type" "MS_cv" "Index" "I_sd" "I_wt"
#> [16] "IAA" "IAA_ESS" "IAL" "IAL_ESS" "C_eq"
#> [21] "C_eq_sd" "E_eq" "abs_I" "I_units" "I_delta"
#> [26] "age_error" "sel_block" "Misc"
class?RCMdataThe two principle sources of data in this manila clam data object are historical catches (Chist) and Catch-at-length composition data (CAL) and three indices (Density survey, ICMP and a weak index prior for testing). Chist is a vector n_years long, CAL is a matrix nyears by n_length_categories:
mydata = In.MC.E$RCMdata
myOM= In.MC.E$OM
myOM@nyears # the operating model historical time period
#> [1] 24
length(mydata@Chist) # check nyears in the OM matches Chist vector length
#> [1] 24
length(mydata@length_bin) # the number of length bins for the length composition data
#> [1] 45
dim(mydata@CAL) # check that the dimensions of the CAL data match up
#> [1] 24 45
dim(mydata@IAA) # check index age observations
#> [1] 24 16 3
dim(mydata@Index) # check that the index data dimensions match up
#> [1] 24 3We are going to add some fake data (the mean catch and indices of the historical series, an old year of CAL data) to these datasets to demonstrate how this is done.
mydata@Chist = c(mydata@Chist, mean(mydata@Chist))
mydata@Index = rbind(mydata@Index,apply(mydata@Index,2,mean,na.rm=T))
mydata@CAL = rbind(mydata@CAL, mydata@CAL[13,]) # 13th year, last year of CAL obs in the dataset is replicated
mydata@IAA = mydata@IAA[c(1:myOM@nyears,myOM@nyears),,] # replicate last year of index age observations
length(mydata@Chist)
#> [1] 25
dim(mydata@CAL)
#> [1] 25 45
dim(mydata@IAA)
#> [1] 25 16 3We have a few things left to do: we need to specify the precision of the new observations of catch and catch-at-length and specify the new number of years in the operating model.
mydata@C_sd = c(mydata@C_sd,mean(mydata@C_sd))
mydata@I_sd = rbind(mydata@I_sd, apply(mydata@I_sd, 2, mean,na.rm=T))
mydata@CAL_ESS = c(mydata@CAL_ESS, mydata@CAL_ESS[13]) # 13th year, last year of CAL age sample size
mydata@IAA_ESS = mydata@IAA_ESS[c(1:myOM@nyears,myOM@nyears),] # replicate last year of index age sample size
myOM@cpars$SLarray = myOM@cpars$SLarray[,,c(1,1:(myOM@nyears+myOM@proyears))] # add another year of size at length
myOM@nyears = myOM@nyears+1 # make sure OM has extra year
myOM@CurrentYr = myOM@CurrentYr+1
Now you have a fully specified RCMdata input object and compatible operating model, we can fit the new model and plot results:
myinput = In.MC.E
myinput$OM = myOM
myinput$RCMdata = mydata
myfit = cond.MC(myinput)
plot(myfit)
13 Resources
13.1 B.C. Hand-harvested invertebrates splash page
For more information about how operating models were specified for B.C. invertebrates you can visit the splash page
Contact (tom@bluematterscience for password)
13.2 B.C. Hand-harvested invertebrates GitHub repository
All of the code for conditioning operating models can be found at a private GitHub repository (contact Mackensie Mazur or Shannon Obradovich for access)
13.3 Inverts R packages
The code for the R packages used in this manaual are available from public GitHub repositories:
13.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.
13.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.
14 References
Hordyk, A., Huynh, Q., Carruthers, T. 2025. OpenMSE: An open-source R package for Management strategy evaluation, available from: https://openmse.com
Huynh, 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.