Sablefish MSE Example

A fully worked example of this MSE framework for Alaska Sablefish follows. This example includes the following features:

multiple operating models
multiple harvest control rules
state-dependent and state-independent recruitment functions
HCRs with stability constraints and harvest caps
HCRs with different reference point levels
Parallel MSE simulations
MSE output processing and plotting.

Defining Operating Models

A base operating model (OM) object is already available to the package to start from. As all of the OMs developed for this example vary only in terms of how recruitment is projected in the future, we will load that base OM object in, and simply modify the recruitment parameter.

nyears <- 100

data_dir <- file.path(here::here(), "data")
load(file.path(data_dir, "sable_om.rda"))

# Also load the historical recruitment timeseries
load(file.path(data_dir, "sable_assessment.rda"))
hist_recruits <- assessment$natage.female[,1]*2

OM1 will project future recruitment by resampling from the historical recruitment timeseries, using the resample_recruits function:

# Normal recruitment
om_rand_recruit <- sable_om
om_rand_recruit$name <- "Random Recruitment"
om_rand_recruit$recruitment$func <- resample_recruits
om_rand_recruit$recruitment$pars <- list(
    hist_recruits = hist_recruits,
    nyears = 10*nyears
)

OM2 projects future recruitment via a Beverton-Holt stock recruit relationship with $h=0.0.85$ and $\sigma_R = 1.20$ . As a state-dependent recruitment function, the underlying beverton_holt R function is a “function factory”, whose returned function is used internally by the MSE to compute recruitment based on current spawning biomass.

# B-H recruitment
om_bh_recruit <- sable_om
om_bh_recruit$name <- "Beverton-Holt Recruitment"
om_bh_recruit$recruitment$func <- beverton_holt
om_bh_recruit$recruitment$pars <- list(
    h = 0.85,
    R0 = 15,
    S0 = sbpr*15,
    sigR = 1.20
)

OM3 projects future recruitment as occurring in distinct regimes: first a “low” regime for 20 years, followed by a “high” regime for 5 years. The regimes each have a uniquely parameterized Beverton-Holt stock recruitment relationships that differ with respect to the $R_0$ and $S_0$ parameters.

om_bhcyclic_recruit <- sable_om
om_bhcyclic_recruit$name <- "Beverton-Holt Cyclic Recruitment"
om_bhcyclic_recruit$recruitment$func <- bevholt_regimes
om_bhcyclic_recruit$recruitment$pars <- list(
    h = 0.85,
    sbpr = sbpr,
    R0 = c(12.5, 50),
    sigR = c(1.20, 1.20),
    nyears = 10*nyears,
    regime_length = c(20, 5),
    starting_regime = 0
)

OM4 assumes that recruitment crashes for the first 20 years of the projection period, while the recruitment for the remainder of the projection period is historically resamples as in OM1.

om_immcrash_recruit <- sable_om
om_immcrash_recruit$name <- "Immediate Crash Recruitment"
om_immcrash_recruit$recruitment$func <- recruits_crash
om_immcrash_recruit$recruitment$pars <- list(
    crash_start_year = 1,
    crash_length = 20,
    crash_value = min(hist_recruits),
    hist_recruits = hist_recruits,
    nyears = 10*nyears
)

The four OM objects are placed in a master list that will later be passed to the MSE function.

om_list <- listN(om_rand_recruit, om_bh_recruit, om_bhcyclic_recruit, om_immcrash_recruit)

Defining Management Procedures

Management procedures (MPs) are specified as a complex list of values that define the specifics of the harvest control rule to be applied, the reference point targets to use, ABC-TAC and TAC-landings reductions, survey frequency, and assessment frequency. A default MP object is available through the setup_mp_options function. The eleven MPs specified here differ only in terms of their HCR function and their reference point target, so the default MP object can, largely be reused.

mp_base <- setup_mp_options()

MP1 is the NPFMC Tier 3 HCR that is currently used to manage sablefish in Alaska. It is a threshold HCR that uses a target reference point of $SPR_{40\%}$ . The HCR function (tier3) outputs in units of fishing mortality (F), as specified by “units”. No harvest cap or stability constraints are applied.

mp_f40 <- mp_base
mp_f40$name <- "F40"
mp_f40$spr_target<- c(0.40, 0.40)
mp_f40$hcr <- list(
    func = tier3,
    extra_pars = NA,
    extra_options = list(
        max_stability = NA,
        harvest_cap = NA
    ),
    units = "F"
)

MP2 and MP3 are simple variations on MP1, using different combination of reference point targets. MP2 used $SPR_{50\%}$ as the reference point target for both fishing mortality and spawning biomass, while MP3 uses $SPR_{50\%}$ as the target fishing mortality rate and $SPR_{40\%}$ as the target spawning biomass level.

mp_f50 <- mp_base
mp_f50$name <- "F50"
mp_f50$ref_points$spr_target <- c(0.5, 0.5)
mp_f50$hcr <- list(
    func = tier3,
    extra_pars = NA,
    extra_options = list(
        max_stability = NA,
        harvest_cap = NA
    ),
    units = "F"
)

mp_b40f55 <- mp_base
mp_b40f55$name <- "B40/F55"
mp_b40f55$ref_points$spr_target <- c(0.55, 0.4)
mp_b40f55$hcr <- list(
    func = tier3,
    extra_options = list(
        max_stability = NA,
        harvest_cap = NA
    ),
    units = "F"
)

MP4 and MP5 are another simple variation on MP1, where stability constraints of 5% and 10% are applied to the ABC. Otherwise, these MPs are defined exactly as MP1.

mp_5perc <- mp_base
mp_5perc$name <- "F40 +/- 5%"
mp_5perc$hcr <- list(
    func = tier3,
    extra_pars = NA,
    extra_options = list(
        max_stability = 0.05,
        harvest_cap = NA
    ),
    units = "F"
)

mp_10perc <- mp_base
mp_10perc$name <- "F40 +/- 10%"
mp_10perc$hcr <- list(
    func = tier3,
    extra_pars = NA,
    extra_options = list(
        max_stability = 0.10,
        harvest_cap = NA
    ),
    units = "F"
)

MP6 and MP7 are another simple variation on MP1, whereby a maximum permissible TAC is implemented. If the recommended TAC (via evaluation of the harvest control rule and after application of the ABC-TAC reduction) exceeds the specified harvest cap, the TAC will be reduced to the harvest cap exactly. Attainment will be applied after the harvest cap has been applied. Harvest caps are 15k, 25k mt respectively.

mp_15cap <- mp_base
mp_15cap$name <- "15k Harvest Cap"
mp_15cap$hcr <- list(
    func = tier3,
    extra_pars = NA,
    extra_options = list(
        max_stability = NA,
        harvest_cap = 15
    ),
    units = "F"
)

mp_25cap <- mp_base
mp_25cap$name <- "25k Harvest Cap"
mp_25cap$hcr <- list(
    func = tier3,
    extra_pars = NA,
    extra_options = list(
        max_stability = NA,
        harvest_cap = 25
    ),
    units = "F"
)

MP8 is a constant fishing mortality rate rule that applies $F_{50}$ at all levels of spawning biomass.

mp_f50chr <- mp_base
mp_f50chr$name <- "Constant F50"
# Set biomass reference point to some small level to simulate a constant F rule
mp_f50chr$ref_points$spr_target <- c(0.50, 0.001)
mp_f50chr$hcr <- list(
    func = chr,
    extra_pars = NA,
    extra_options = list(
        max_stability = NA,
        harvest_cap = NA
    ),
    units = "F"
)

MP9 and MP10 are approximate implementations of the harvest control rules used by other management bodies for Pacific sablefish.

MP9 is an implementation of the Pacific Fisheries Management Council’s “40-10” HCR that is used to manage sablefish along the west coast of the continental United States. This HCR traditionally operates at the TAC level, but is implemented here so as to output an ABC, as the other HCR functions do. MP9 uses a target fishing mortality rate of $SPR_{45\%}$ to approximate $F_{\text{MSY}}$ , and, unlike the other HCR functions, outputs in units of catch (“TAC”) rather than F.

MP10 is an implementation of DFO Canada’s sablefish HCR that is used for sablefish management in British Columbia, Canada. This HCR is a threshold rule with a maximum harvest rate of 5.5%, which declines towards 0 when $\text{SSB} <= 0.60B_{MSY}$ . $B_{\text{MSY}}$ is approximated as the spawning biomass resulting from fishing at a $F=F_{45\%}$ , akin to in MP9.

mp_pfmc4010 <- mp_base
mp_pfmc4010$name <- "PFMC 40-10"
mp_pfmc4010$ref_points$spr_target <- c(0.45, 0.40)
mp_pfmc4010$hcr <- list(
    func = pfmc4010,
    extra_pars = NA,
    extra_options = list(
        max_stability = NA,
        harvest_cap = NA
    ),
    units = "TAC"
)

mp_bcsable <- mp_base
mp_bcsable$name <- "British Columbia"
mp_bcsable$ref_points$spr_target <- c(0.45, 0.45)
mp_bcsable$hcr <- list(
    func = bc_sable,
    extra_pars = NA,
    extra_options = list(
        max_stability = NA,
        harvest_cap = NA
    ),
    units = "F"
)

MP11 is a simple rule that prohibits fishing at any level of biomass.

mp_f00chr <- mp_base
mp_f00chr$name <- "No Fishing"
mp_f00chr$hcr <- list(
    func = chr,
    extra_pars = NA,
    extra_options = list(
        max_stability = NA,
        harvest_cap = NA
    ),
    units = "TAC"
)
mp_f00chr$ref_points$spr_target <- c(1, 0.001)

The eleven MP objects are placed in a master list that will later be passed to the MSE function.

hcr_list <- listN(
    mp_f40, mp_f50, mp_b30f40, mp_b40f50,
    mp_5perc, mp_10perc, mp_10perc_up, mp_15perc,
    mp_15cap, mp_25cap,
    mp_f50chr,
    mp_pfmc4010, mp_bcsable,
    mp_f00chr
)

Defining the `mse_options` Object

An mse_options list object is also defined to provide additional options to the MSE function. Critically, this list object defines: the number of years to project for, the length of the OM conditioning period, the simulation year in which to begin projecting recruitment using the provided recruitment function, and whether the estimation method should be run or not.

The default values are all appropriate for this case study, and can be accessed using setup_mse_options.

mse_options_base <- setup_mse_options()
mse_options <- mse_options_base
mse_options$n_spinup_years <- 54
mse_options$recruitment_start_year <- 54
mse_options$n_proj_years <- 75

mse_options_list <- listN(mse_options)

Running the MSE

Now that all of the OMs and MPs have been appropriately defined, the full set of MSE simulation can be run. Here we will run each combination of OM and MP for 9 random simulations in parallel. The parallel computing overhead is handled internally by the MSE function and uses 2 fewer than the available number of compute cores on the machine.

nsims <- 20
seed_list <- sample(1:(1000*nsims), nsims)  # Draw 20 random seeds
model_runs <- run_mse_multiple(
    om_list, 
    hcr_list, 
    seed_list,
    nyears=100,
    mse_options_list=mse_options_list,
    diagnostics = TRUE,
    save=TRUE
)

Depending on the number of OMs, MPs, seeds, and other options, this function may take several hours to complete. A relatively fast machine, with >10 compute cores should take ~2 hours to complete all of the simulations as specified above.

Processing MSE Results

The above MSE will generate results from 44 unique models (11 HCRs across 5 OMs), each model consisting of results for 20 simulations, with each simulation having lasted 75 years. The amount of output data often makes processing MSE results tricky. Multiple helper functions are provided by the package to make processing results easier.

Before prcoessing output, create a data.frame that specifies which models in the model_runs used which combination of OM and MP object. Often this can be simply accomplished using expand.grid, as this is what is used internally by run_mse_multiple(). For readability and plotting purposes, it is often useful to assign each OM and HCR a human readable name in the same order as the objects appear in their respective lists (om_list for OMs and hcr_list for MPs).

om_names <- unlist(lapply(om_list, \(x) x$name))
hcr_names <- unlist(lapply(hcr_list, \(x) x$name))

extra_columns2 <- expand.grid(
    om = unlist(lapply(om_list, \(x) x$name)),
    hcr = unlist(lapply(hcr_list, \(x) x$name))
)

Five helper functions are predefined to facilitate working with common outputs from the MSE: spawning biomass, fishing mortality, recruitment, ABC and TAC, and landed catch. Each function is defined very similarly, taking the list of MSE output objects (model_runs) and the data.frame of OMs and HCRs that apply to each model (extra_columns). These are used internally by each function to assign the correct OM and HCR name to each model, and pull the correct data from the MSE outputs.

ssb_data <- get_ssb_biomass(model_runs, extra_columns2, sable_om$dem_params, hcr_filter=hcr_names, om_filter=om_names)

f_data <- get_fishing_mortalities(model_runs, extra_columns2, hcr_filter=hcr_names, om_filter=om_names)

abctac <- get_management_quantities(model_runs, extra_columns2, spinup_years=common_trajectory, hcr_filter=hcr_names, om_filter=om_names)

catch_data <- get_landed_catch(model_runs, extra_columns2, hcr_filter=publication_hcrs, om_filter=om_names)

These helper functions return tibbles in long format that can be easily provided to ggplot for plotting or used within other tidyverse style data processing pipelines. A plotting function is also defined alongside each data processing function above, and will return a ggplot object that can be further modified outside of the function. For each plotting function, the parameter v1 corresponds to the column name assigned to the “color” aesthetic, and the parameter v2 corresponds to the column names to facet by. Facetting by multiple variables (e.g. a facet grid) is not currently supported by default.

plot_ssb(ssb_data, v1="hcr", v2="om", v3=NA, common_trajectory=common_trajectory, show_est = FALSE)

plot_relative_ssb(ssb_data, v1="hcr", v2="om", common_trajectory = common_trajectory, base_hcr = "No Fishing")

plot_fishing_mortalities(f_data, v1="hcr", v2="om", common_trajectory = common_trajectory, show_est=FALSE)

plot_abc_tac(abctac, v1="hcr", v2="om", common_trajectory=common_trajectory)

plot_landed_catch(catch_data, v1="hcr", v2="om", common_trajectory = common_trajectory)

Computing Performance Metrics

Helper functions are also provided to compute performance metrics. Individual functions for each performance metric are available (though outputs may require additional processing), or the performance_metric_summary function can be used to quickly compute all of them. Like with the other data processing functions, a default plotting functions for performance metrics is also available.