Processing MSE Output Data

MSEs output large amounts of data – numbers-at-age across years, landed catch by year, fishing mortality by year, recommended harvest levels by year, etc. – and this data output quickly grows as the number of operating models (OMs), management procedures (MPs), and simulations grows. SablefishMSE provides several wrapper functions to facilitate easier processing of large MSE outputs and computation of common derived quantities (spawning biomass, total biomass, total catch, etc.)

MSE Output Data

All three of the MSE functions (run_mse, run_mse_parallel, and run_mse_multiple) output the same types of data in a similar format.

run_mse - data is output as a named list of multidimensional arrays of maximum dimension [n_proj_years+n_spinup_years, nages, nsexes, nregions, nfleets].
run_mse_parallel - data is output as a named list of multidimensional arrays the same as run_mse but with an added dimension referring to the simulation
run_mse_multiple - data is output at a list of completed MSE simulations from run_mse_parallel, with each list element referring to one combination of OM, MP, and MSE option.

Description of Output Data

The following data is output with every call to run_mse and is consequently available with calls to run_mse_parallel and run_mse_multiple:

land_caa - landed catch-at-age from the OM ([nyears, nages, nsexes, nregions, nfleets]). Note that this is in units of weight rather than numbers.
disc_caa - discarded catch-at-age from the OM ([nyears, nages, nsexes, nregions, nfleets]). Note that this is in units of weight rather than numbers.
caa - total catch-at-age (landed + discarded) from the OM ([nyears, nages, nsexes, nregions, nfleets]). Note that this is in units of weight rather than numbers.
faa - fishing mortality -at-age from the OM ([nyears, nages, nsexes, nregions, nfleets])
faa_est - fishing mortality -at-age estimated from the EM ([nyears, nages, nsexes, nregions, nfleets])
naa - numbers-at-age from the OM ([nyears, nages, nsexes, nregions])
naa_est - numbers-at-age estimated from the EM ([nyears, nages, nsexes, nregions])
hcr_f - the fishing mortality recommended by the HCR ([n_proj_years, 1])
abc - the Acceptable Biological Catch resulting from applying the HCR ([n_proj_years, 1])
tac - the Total Allowable Catch resulting from applying the HCR and subsequent stability constraint or harvest caps ([n_proj_years, 1])
exp_land - the expected landings after accounting for partial attainment of the TAC ([n_proj_years, 1])
global_rec_devs - recruitment deviates from the OM ([nyears, 1])

Additionally, if diagnostics=TRUE for run_mse_parallel or run_mse_multiple, the following outputs are also available:

survey_obs - observation data generated by the OM
model_outs - completed TMB model objects

The functions discussed below were not developed with the intent to be used on the survey_obs or model_outs output data types.

Reading MSE Output from a File

If run_mse_multiple is used, users can optionally specify save=TRUE to save MSE outputs to a file. These outputs can be read back into R using the provided get_saved_model_runs function. Specifying the name of the OMS and MPs to read back from file as a vector will filter the files read back in to only those that match the specified OMS and MPs. If no filters are provided, all saved model runs will be read back in:

om_filter <- c("OM1", "OM2")
hcr_filter <- c("MP1", "MP2")

model_runs <- get_saved_model_runs(om_filter=om_filter, hcr_filter=hcr_filter)

This function returns a list of two objects: model_runs, which is a list of the MSE output data for each model run, and extra_columns, which is a data.frame containing the OM and MP names corresponding to each model run.

Depending on the number of model runs being read back in, this function often takes several minutes to complete. For very large numbers of models, it is not possible to read all model runs back into R at once. In such scenarios, existing MSE processing functions (see below) can be used to process data file-by-file, and then combine and summarize across files at the end.

Existing MSE Processing Functions

Several processing functions are available through the SablefishMSE package allowing for easy computation of common fisheries quantities.

get_ssb_biomass - calculates spawning biomass (SSB) and total biomass for every year across all model runs and simulations.
get_fishing_mortalities - calculates fleet-specific and total fishing mortality for every year across all model runs and simulations. Total fishing mortality is assumed to be the maximum age-specific fishing mortality aggregated across fleets.
get_recruits - calculates recruitment for every year across all model runs and simulations. This function does not handle recruitment deviates.
get_landed_catch - calculates fleet-specific and total landed catch for every year across all model runs and simulations.
get_dynamic_economic_value - calculates annual economic value of the fixed gear fleet using a dynamic price model that assumes prices decline with total catch volume (see Goethel et al. 2025 for details on the price model).
get_average age - calculates average age of the population for every year across all model runs and simulations.
get_abi - calculates population ABI (Griffiths et al. 2023) for every year across all model runs and simulations.
get_management_quantities - aggregates ABC, TAC, and expected landings data for every year across all model runs and simulations. Historical ABC, TAC, and realized landings are hard coded into the function.

All functions generally accept the same parameters, and return a long format tibble. The first two parameters, model_runs and extra_columns, are outputs of get_saved_model_runs and are required for all functions. If model_runs=NULL, the function will attempt to iteratively read in and apply the underlying the process function to, each saved model run file as appropriate. When model_runs=NULL, the extra_columns parameter must stll be provided as a data.frame containing the OM and MP names corresponding to each model run file; this is most easily accomplished as extra_columns = expand.grid(om=om_filter, hcr=hcr_filter).

The om_filter and hcr_filter parameters can be used to specify which OMs and MPs to read in and process.

Additional parameters may be required for some functions.

For example:

om_filter <- c("OM1", "OM2")
hcr_filter <- c("MP1", "MP2")

mse_runs <- get_saved_model_runs(om_filter=om_filter, hcr_filter=hcr_filter)
model_runs <- mse_runs$model_runs
extra_columns <- mse_runs$extra_columns

# Calculate SSB and total biomass across all model runs and simulations
ssb_biomass <- get_ssb_biomass(model_runs=model_runs, extra_columns=extra_columns, om1$dem_params, om_filter=om_filter, hcr_filter=hcr_filter)

# Alternatively, if `get_saved_model_runs` is not used
rm(model_runs, extra_columns)
mse_runs <- NULL

# Itertaively read in and process each model run file as appropriate
ssb_biomass <- get_ssb_biomass(model_runs=NULL, extra_columns=NULL, om1$dem_params, om_filter=om_filter, hcr_filter=hcr_filter)

Performance Metrics

Many functions are available for calculating common and sablefish-specific performance metrics from MSE output data. All functions calculate their performance metric across the projection period only, and across all model runs.

average_catch - average annual catch
total_catch - total catch (sum of catch over all years)
prop_years_catch - proportion of projection years in which catch exceeded some threshold value
average_ssb - average annual spawning biomass
prop_low_biomass - proportion of projection years in which SSB was $<0.35\text{SSB}_{1960}$
time_below_bref - average number of projection years in which SSB was $<B_{ref} (where $B_{ref}$ is the MP specific biomass reference point).
average_age - average age of an individual in the population
average_abi - average ABI of the population (Griffiths et al. 2023)
average_annual_catch_variation - annual year-to-year variation in catch, calculated as the average of the absolute value of the year-to-year change in catch divided by the average catch across all years
average_proportion_catch_large - average annual proportion of the landed catch >9yo
average_proportion_biomass_old - average annual proportion of the biomass >21yo
annual_value - average annual value of the catch assuming a fixed relative price per age class
dynamic_annual_value - average annual value of the catch assuming prices fluctuate with landed catch volume

Two resiliency metrics can also be calculated, but should only be calculated under crash recruitment OM scenarios, as they internally rely on knowing over what years a recruitment crash occurs.

biomass_crash_time - average time required for the SSB to decline to $<0.175\text{SSB}_{1960}$. This metric is only calculated over the first 30 years of the projection period.
biomass_recovery_time - average time required for the SSB to recover to $>0.35\text{SSB}_{1960}$. This metric is calculated starting year 31 of the projection period.

For additional information on any of these functions, see their function documentation.

`performance_metric_summary` Function

A summary function is available to compute multiple performance metrics and concatenate them into a single summary tibble object: performance_metric_summary.

performance_metric_summary(
    model_runs=model_runs, 
    extra_columns=extra_columns, 
    dem_params=om$dem_params, 
    ref_naa=ref_naa,
    hcr_filter=c("MP1", "MP2"),
    om_filter=c("OM1", "OM2"),
    interval_widths=c(0.50, 0.95), 
    # List of metrics to calculate
    metric_list = c("avg_catch", "avg_variation", "avg_ssb", "avg_age", "prop_years_lowssb")
)