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.
library(wham)
library(whamMSE)
main.dir = here::here()
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
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!
wham
InputHere 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)
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)
om_with_data <- update_om_fn(om, seed = 123, random = random)
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)
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)