This vignette demonstrates how to use the whamMSE (Woods Hole Assessment Model Management Strategy Evaluation) package to incorporate implementation error into a management strategy evaluation (MSE). Note: all code shown here is for demonstration purposes only; the results are for illustrative use and should not be interpreted as reflecting real-world stock assessments.

1. Load “WHAM” and “whamMSE”:

library(wham)
library(whamMSE)

main.dir = here::here()

2. Generate Basic Information

year_start  <- 1  # starting year in the burn-in period
year_end    <- 20  # end year in the burn-in period
MSE_years   <- 3     # number of years in the feedback loop
# Note: no need to include MSE_years in simulation-estimation 

info <- generate_basic_info(n_stocks = 1,
                            n_regions = 1,
                            n_indices = 1,
                            n_fleets = 1,
                            n_seasons = 1,
                            base.years = year_start:year_end,
                            n_feedback_years = MSE_years,
                            life_history = "medium",
                            n_ages = 12) 

basic_info = info$basic_info # collect basic information
catch_info = info$catch_info # collect fleet catch information
index_info = info$index_info # collect survey information
F_info = info$F # collect fishing information

3. Configure Numbers-at-Age (NAA)

sigma        <- "rec+1"
re_cor       <- "iid"
ini.opt      <- "equilibrium" 
sigma_vals   <-  array(0.5, dim = c(n_stocks, n_regions, n_ages)) # NAA sigma

NAA_re <- list(N1_model=rep(ini.opt,n_stocks),
               sigma=rep(sigma,n_stocks),
               cor=rep(re_cor,n_stocks),
               recruit_model = 2,  # rec random around the mean
               sigma_vals = sigma_vals) # NAA_where must be specified in basic_info!

4. Generate wham Input

Here we use prepare_wham_input() function to generate a wham input using the basic information we set above:

input <- prepare_wham_input(basic_info = basic_info, 
                            NAA_re = NAA_re,
                            catch_info = catch_info, 
                            index_info = index_info, 
                            F = F_info)

5. Generate Operating Model

random = input$random # check what processes are random effects
input$random = NULL # so inner optimization won't change simulated RE
om <- fit_wham(input, do.fit = F, do.brps = T, MakeADFun.silent = TRUE)
# Note: do.fit must be FALSE (no modeling fitting yet)

6. Generate Dataset

om_with_data <- update_om_fn(om, seed = 123, random = random)

7. Specify Assessment Interval

assess.interval <- 3 # Note: assessment interval is 3 years, given the feedback period is 3 years, there will be only 1 assessment
base.years      <- year_start:year_end # Burn-in period
first.year      <- head(base.years,1)
terminal.year   <- tail(base.years,1)
assess.years    <- seq(terminal.year, tail(om$years,1)-assess.interval,by = assess.interval)

8. Conduct MSE with Implementation Error

In this example, we introduce implementation error to simulate deviations between catch advice and actual realized catch. Implementation error accounts for uncertainty or non-compliance when applying catch advice in management, which can influence stock dynamics in the MSE framework.

We support the following implementation error options:

  • "lognormal": Multiplies annual catch advice by a log-normal random variable drawn from exp(N(mean, sd)). This captures multiplicative error with log-normal variability.
  • "normal": Adds additive noise by multiplying annual catch advice by 1 + N(mean, sd). This is centered at 1 (no bias) by default if mean = 0.
  • "uniform": Multiplies annual catch advice by a uniform random multiplier between min and max. This simulates bounded but random variability. For example, setting min = 1.1 and max = 1.2 represents a case of consistent overreporting, where the annual realized catch is randomly drawn to be 10–20% higher than the annual catch advice.
  • "constant": Applies a fixed multiplier to the catch advice (e.g., 1.2 means 20% overshoot every year).

For this run, we used the "constant" method, with a constant_value of 1.2, meaning the realized catch was consistently 20% higher than the advised catch. This is useful to test performance under persistent overharvest scenarios.

mod = loop_through_fn(om = om_with_data,
                      em_info = info, 
                      random = random,
                      M_em = M, 
                      sel_em = sel,
                      NAA_re_em = NAA_re, 
                      age_comp_em = "multinomial",
                      em.opt = list(separate.em = FALSE, 
                                    separate.em.type = 3, 
                                    do.move = FALSE, 
                                    est.move = FALSE),
                      implementation_error = list(method = "constant",
                                                  mean = 0,
                                                  cv = NULL,
                                                  sd = NULL,
                                                  min = NULL,
                                                  max = NULL,
                                                  constant_value = 1.2),
                      assess_years = assess.years, 
                      assess_interval = assess.interval, 
                      base_years = base.years,
                      year.use = 20,
                      add.years = TRUE, 
                      seed = 123)