Single Stock Multi-Fleets MSE
Chengxue Li
Last update: 2025-06-24
This vignette introduces how to run MSE for a single stock with multiple fleets and surveys:
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 <- 30 # 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 = 2,
n_fleets = 2,
n_seasons = 1,
catch_info = list(catch_cv = 0.1, catch_Neff = 100, use_agg_catch = 1, use_catch_paa = 1),
index_info = list(index_cv = 0.1, index_Neff = 100, fracyr_indices = 0.5, q = 0.2,
use_indices = 1, use_index_paa = 1, units_indices = 2, units_index_paa = 2),
base.years = year_start:year_end,
n_feedback_years = MSE_years,
life_history = "medium",
n_ages = 12,
fracyr_spawn = 0.5)
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 information3. Configure selecitvity and natural mortality
n_stocks <- as.integer(basic_info['n_stocks'])
n_regions <- as.integer(basic_info['n_regions'])
n_fleets <- as.integer(basic_info['n_fleets'])
n_indices <- as.integer(basic_info['n_indices'])
n_ages <- as.integer(basic_info['n_ages'])
# Selectivity Configuration
fleet_pars <- c(5,1)
index_pars <- c(2,1)
sel <- list(model=rep("logistic",n_fleets+n_indices),
initial_pars=c(rep(list(fleet_pars),n_fleets),rep(list(index_pars),n_indices)))
# M Configuration
M <- list(model="constant",initial_means=array(0.2, dim = c(n_stocks,n_regions,n_ages)))4. 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!5. 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,
selectivity = sel,
M = M,
NAA_re = NAA_re,
catch_info = catch_info,
index_info = index_info,
F = F_info)6. Change observation error for survey indices in the feedback period
agg_index_sigma = input$data$agg_index_sigma
agg_index_sigma[21:50,] = 0.5 # Increase CV for both survey indices in the feedback period
index_Neff = input$data$index_Neff
index_Neff[21:50,] = 50 # Decrease ESS for both survey indices in the feedback period
input <- update_input_index_info(input, agg_index_sigma, index_Neff) # Update input file7. 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)8. Generate dataset
om_with_data <- update_om_fn(om, seed = 123, random = random)9. Specify assessment interval and assessment year in the feedback loop
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)10. Use the EM same as OM to conduct MSE
mod1 = 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),
assess_years = assess.years,
assess_interval = assess.interval,
base_years = base.years,
year.use = 20,
update_index_info = list(agg_index_sigma = agg_index_sigma,index_Neff = index_Neff),
# Assuming we know the true observation error of the data
# agg_index_sigma and index_Neff can be different than the true values in the OM
add.years = TRUE,
# assessment will use 20 years of data from historical period + new years in the feedback period
seed = 123,
save.sdrep = TRUE,
save.last.em = TRUE,
do.retro = FALSE, # Perform retrospective analysis
do.osa = FALSE) # Perform OSA residual analysis11. Use the EM with semi selfweighting likelihoods for age composition data
mod2 = 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 = "logistic-normal-miss0", # semi self-weighting
em.opt = list(separate.em = FALSE,
separate.em.type = 3,
do.move = FALSE,
est.move = FALSE),
assess_years = assess.years,
assess_interval = assess.interval,
base_years = base.years,
year.use = 20,
update_index_info = list(agg_index_sigma = agg_index_sigma,index_Neff = index_Neff),
# Assuming we know the true observation error of the data
# agg_index_sigma and index_Neff can be different than the true values in the OM
add.years = TRUE,
# assessment will use 20 years of data from historical period + new years in the feedback period
seed = 123,
save.sdrep = TRUE,
save.last.em = TRUE,
do.retro = FALSE, # Perform retrospective analysis
do.osa = FALSE) # Perform OSA residual analysis12. Performance
par(mfrow = c(1,3))
plot(mod1$om$rep$SSB, type = "o", col = "red", main = "SSB")
lines(mod2$om$rep$SSB, type = "o", col = "blue")
legend("topright", legend=c("Multinomial", "Logistic-normal"),
col=c("red", "blue"), lty=c(1,1), cex=0.8)
plot(mod1$om$rep$Fbar[,3], type = "o", col = "red", main = "F", ylim = c(0,1))
lines(mod2$om$rep$Fbar[,3], type = "o", col = "blue")
plot(rowSums(mod1$om$rep$pred_catch), type = "o", col = "red", main = "Catch")
lines(rowSums(mod2$om$rep$pred_catch), type = "o", col = "blue")
Figure 1