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, sable_om, is already available in 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 this base OM object and simply modify the recruitment parameter(s).

library(SablefishMSE)
data(sable_om)

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

# Normal recruitment
data(om_rand_recruit)
om_rand_recruit$recruitment

## $func
## function(hist_recruits, nyears, seed){
##     set.seed(seed)
##     return(sample(hist_recruits, size=nyears, replace=TRUE))
## }
## 
## $pars
## $pars$hist_recruits
##     1960     1961     1962     1963     1964     1965     1966     1967 
## 28.02120 29.71580 31.69360 33.56320 34.66580 35.07580 33.41820 30.67780 
##     1968     1969     1970     1971     1972     1973     1974     1975 
## 27.83260 25.40080 23.22660 21.31440 19.81614 18.79288 18.57922  9.06294 
##     1976     1977     1978     1979     1980     1981     1982     1983 
## 10.31614 10.85754 13.65474 83.98160 41.82600 19.33290 74.28580 36.24560 
##     1984     1985     1986     1987     1988     1989     1990     1991 
## 14.32186 15.99548 23.86060 10.12922  7.20310  8.38934 11.69706 22.77260 
##     1992     1993     1994     1995     1996     1997     1998     1999 
##  7.86546 24.32680  6.83952  7.27954 12.16022 20.26580 10.05634 36.35860 
##     2000     2001     2002     2003     2004     2005     2006     2007 
## 15.98118 16.15600 42.62800 13.30080  9.89198 11.78916  7.59538 10.23582 
##     2008     2009     2010     2011     2012     2013     2014     2015 
##  9.69452 15.73100 21.31020 10.30400 11.22174  4.69700  7.48688 14.40154 
##     2016     2017     2018     2019     2020     2021     2022     2023 
## 49.54060 21.69120 95.37040 86.48080 40.60020 75.09140 42.68520 26.71460 
## 
## $pars$nyears
## [1] 1100

OM2 (om_bh_recruit) projects future recruitment via a Beverton-Holt stock recruit relationship with \(h=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.

data(om_bh_recruit)
om_bh_recruit$recruitment

## $func
## function(h, R0, S0, sigR, seed){
##     # note that the set.seed() call needs to happen
##     # outside of the returned function, or else there
##     # will be no random variability in recruitment draws
##     set.seed(seed)
##     function(ssb, y){
##         bh <- (4*R0*h*ssb)/((1-h)*R0*(S0/R0) + (5*h - 1)*ssb)
##         return(
##             bh #* exp(rnorm(1, mean=-sigR*sigR/2, sd=sigR)) # lognormal
##         )
##     }
## }
## 
## $pars
## $pars$h
## [1] 0.85
## 
## $pars$R0
## [1] 15
## 
## $pars$S0
## [1] 328.6342
## 
## $pars$sigR
## [1] 1.2

OM3 (om_bhcyclic_recruit) 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. This type of recruitment function is available with the provided bevholt_regimes function.

data(om_bhcyclic_recruit)
om_bhcyclic_recruit$recruitment

## $func
## function(h, sbpr, R0s, sigRs, nyears, regime_length, starting_regime, seed){
##     set.seed(seed)
##     function(ssb, y){
##         regs <- rep(NA, nyears)
##         curr_regime <- starting_regime
##         y1 <- 1
##         while(y1 < nyears){
##             reg_len <- regime_length[curr_regime+1]
##             regs[y1:(y1+reg_len-1)] <- curr_regime+1
##             curr_regime <- !curr_regime
##             y1 <- y1+reg_len
##         }
##         R0s <- R0s[regs]
##         sigRs <- sigRs[regs]
##         # print(regs)
##         # print(R0s)
##         # print(sigRs)
## 
##         R0 <- R0s[y]
##         sigR <- sigRs[y]
##         S0 <- sbpr*R0
##         bh <- (4*R0*h*ssb)/((1-h)*R0*(S0/R0) + (5*h - 1)*ssb)
##         return(
##             bh #* exp(rnorm(1, mean=-sigR*sigR/2, sd=sigR)) # lognormal
##         )
##     }
## }
## 
## $pars
## $pars$h
## [1] 0.85
## 
## $pars$sbpr
## [1] 21.90895
## 
## $pars$R0
## [1] 12.5 50.0
## 
## $pars$sigR
## [1] 1.2 1.2
## 
## $pars$nyears
## [1] 1100
## 
## $pars$regime_length
## [1] 20  5
## 
## $pars$starting_regime
## [1] 0

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

om_list <- afscOM::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()
mp_base

## $hcr
## $hcr$func
## NULL
## 
## $hcr$extra_pars
## [1] NA
## 
## $hcr$extra_options
## $hcr$extra_options$max_stability
## [1] NA
## 
## $hcr$extra_options$harvest_cap
## [1] NA
## 
## 
## $hcr$units
## NULL
## 
## 
## $ref_points
## $ref_points$spr_target
## [1] 0.4
## 
## $ref_points$rp_start_age
## [1] 1
## 
## $ref_points$rp_hyperallometry
## [1] 1
## 
## 
## $management
## $management$abc_tac_reduction
## [1] 1
## 
## $management$tac_land_reduction
## [1] 1
## 
## 
## $survey_frequency
## [1] 1
## 
## $assessment_frequency
## [1] 1

MP1 (mp_f40) is an implementation of 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.

data(mp_f40)
mp_f40

## $hcr
## $hcr$func
## function(ref_pts, naa, dem_params, avgrec, cutoff_age=1, hyperallometry=1){
##     nages <- afscOM::get_model_dimensions(dem_params$sel)$nages
##     a <- cutoff_age-1
##     ssb <- apply(naa[,a:nages,1,]*(dem_params$waa[,a:nages,1,,drop=FALSE]^hyperallometry)*dem_params$mat[,a:nages,1,,drop=FALSE], 1, sum)
##     return(
##         npfmc_tier3_F(ssb, ref_pts$Bref, ref_pts$Fref)
##     )
## }
## 
## $hcr$extra_pars
## [1] NA
## 
## $hcr$extra_options
## $hcr$extra_options$max_stability
## [1] NA
## 
## $hcr$extra_options$harvest_cap
## [1] NA
## 
## 
## $hcr$units
## [1] "F"
## 
## 
## $ref_points
## $ref_points$spr_target
## [1] 0.4
## 
## $ref_points$rp_start_age
## [1] 1
## 
## $ref_points$rp_hyperallometry
## [1] 1
## 
## 
## $management
## $management$abc_tac_reduction
## [1] 1
## 
## $management$tac_land_reduction
## [1] 1
## 
## 
## $survey_frequency
## [1] 1
## 
## $assessment_frequency
## [1] 1
## 
## $name
## [1] "F40"

MP2 (mp_f50) is a simple variation on MP1, where the target reference points are based on \(SPR_{50\%}\) instead of \(SPR_{40\%}\).

data(mp_f50)
mp_f50

## $hcr
## $hcr$func
## function(ref_pts, naa, dem_params, avgrec, cutoff_age=1, hyperallometry=1){
##     nages <- afscOM::get_model_dimensions(dem_params$sel)$nages
##     a <- cutoff_age-1
##     ssb <- apply(naa[,a:nages,1,]*(dem_params$waa[,a:nages,1,,drop=FALSE]^hyperallometry)*dem_params$mat[,a:nages,1,,drop=FALSE], 1, sum)
##     return(
##         npfmc_tier3_F(ssb, ref_pts$Bref, ref_pts$Fref)
##     )
## }
## 
## $hcr$extra_pars
## [1] NA
## 
## $hcr$extra_options
## $hcr$extra_options$max_stability
## [1] NA
## 
## $hcr$extra_options$harvest_cap
## [1] NA
## 
## 
## $hcr$units
## [1] "F"
## 
## 
## $ref_points
## $ref_points$spr_target
## [1] 0.5 0.5
## 
## $ref_points$rp_start_age
## [1] 1
## 
## $ref_points$rp_hyperallometry
## [1] 1
## 
## 
## $management
## $management$abc_tac_reduction
## [1] 1
## 
## $management$tac_land_reduction
## [1] 1
## 
## 
## $survey_frequency
## [1] 1
## 
## $assessment_frequency
## [1] 1
## 
## $name
## [1] "F50"

MP3 (mp_5perc) is another simple variation on MP1, where a stability constraint of 5% are applied to the ABC. Otherwise, these MPs are defined exactly as MP1. The stability constraints limit the year-to-year change in ABC to the specified percentage, such that the ABC cannot increase or decrease by more than the specified percentage from one year to the next. Stability constraints are implemented in catch units, and are thus applied after the tier3 HCR function has output a recommended F.

MP4 (mp_25cap) is another simple variation on MP1, whereby a maximum permissible TAC is implemented. If the recommended TAC (via evaluation of the HCR function) 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.

data(mp_5perc)
data(mp_25cap)

mp_5perc

## $hcr
## $hcr$func
## function(ref_pts, naa, dem_params, avgrec, cutoff_age=1, hyperallometry=1){
##     nages <- afscOM::get_model_dimensions(dem_params$sel)$nages
##     a <- cutoff_age-1
##     ssb <- apply(naa[,a:nages,1,]*(dem_params$waa[,a:nages,1,,drop=FALSE]^hyperallometry)*dem_params$mat[,a:nages,1,,drop=FALSE], 1, sum)
##     return(
##         npfmc_tier3_F(ssb, ref_pts$Bref, ref_pts$Fref)
##     )
## }
## 
## $hcr$extra_pars
## [1] NA
## 
## $hcr$extra_options
## $hcr$extra_options$max_stability
## [1] 0.05
## 
## $hcr$extra_options$harvest_cap
## [1] NA
## 
## 
## $hcr$units
## [1] "F"
## 
## 
## $ref_points
## $ref_points$spr_target
## [1] 0.4
## 
## $ref_points$rp_start_age
## [1] 1
## 
## $ref_points$rp_hyperallometry
## [1] 1
## 
## 
## $management
## $management$abc_tac_reduction
## [1] 1
## 
## $management$tac_land_reduction
## [1] 1
## 
## 
## $survey_frequency
## [1] 1
## 
## $assessment_frequency
## [1] 1
## 
## $name
## [1] "F40 +/- 5%"

mp_25cap

## $hcr
## $hcr$func
## function(ref_pts, naa, dem_params, avgrec, cutoff_age=1, hyperallometry=1){
##     nages <- afscOM::get_model_dimensions(dem_params$sel)$nages
##     a <- cutoff_age-1
##     ssb <- apply(naa[,a:nages,1,]*(dem_params$waa[,a:nages,1,,drop=FALSE]^hyperallometry)*dem_params$mat[,a:nages,1,,drop=FALSE], 1, sum)
##     return(
##         npfmc_tier3_F(ssb, ref_pts$Bref, ref_pts$Fref)
##     )
## }
## 
## $hcr$extra_pars
## [1] NA
## 
## $hcr$extra_options
## $hcr$extra_options$max_stability
## [1] NA
## 
## $hcr$extra_options$harvest_cap
## [1] 25
## 
## 
## $hcr$units
## [1] "F"
## 
## 
## $ref_points
## $ref_points$spr_target
## [1] 0.4
## 
## $ref_points$rp_start_age
## [1] 1
## 
## $ref_points$rp_hyperallometry
## [1] 1
## 
## 
## $management
## $management$abc_tac_reduction
## [1] 1
## 
## $management$tac_land_reduction
## [1] 1
## 
## 
## $survey_frequency
## [1] 1
## 
## $assessment_frequency
## [1] 1
## 
## $name
## [1] "25k Harvest Cap"

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

hcr_list <- afscOM::listN(
    mp_f40, mp_f50, 
    mp_5perc, mp_25cap
)

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 <- 25

mse_options

## $n_proj_years
## [1] 25
## 
## $n_spinup_years
## [1] 54
## 
## $recruitment_start_year
## [1] 54
## 
## $run_estimation
## [1] TRUE

mse_options_list <- afscOM::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 20 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 40 unique models (10 HCRs across 4 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).

# Quick Way to get the names of the OMs and HCRs in the same order as they appear in the lists
om_names <- unlist(lapply(om_list, \(x) x$name))
hcr_names <- unlist(lapply(hcr_list, \(x) x$name))

extra_columns <- expand.grid(
    om = om_names,
    hcr = hcr_names
)

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_columns, sable_om$dem_params, hcr_filter=hcr_names, om_filter=om_names)

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

abctac <- get_management_quantities(model_runs, extra_columns, spinup_years=mse_options$n_spinup_years, hcr_filter=hcr_names, om_filter=om_names)

catch_data <- get_landed_catch(model_runs, extra_columns, hcr_filter=hcr_names, 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.

common_trajectory <- mse_options$n_spinup_years # the number of years of historical data before the recruitment and HCR functions begin applying

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

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.