This vignette demonstrates how to use the whamMSE (Woods
Hole Assessment Model Management Strategy Evaluation) package to reduce
data availability in the assessment (e.g., deal with data limited
situations in the future) in 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 <- 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 informationn_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)))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,
selectivity = sel,
M = M,
NAA_re = NAA_re,
catch_info = catch_info,
index_info = index_info,
F = F_info)agg_index_sigma = input$data$agg_index_sigma
agg_index_sigma[21:50,] = 1 # Increase CV for both survey indices in the feedback period
index_Neff = input$data$index_Neff
index_Neff[21:50,] = 30 # Decrease ESS for both survey indices in the feedback period
input <- update_input_index_info(input, agg_index_sigma, index_Neff) # Update input filerandom = 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)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), # Must have this!
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 analysisinput <- 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)
agg_index_sigma = input$data$agg_index_sigma
agg_index_sigma[21:50,] = 1 # Increase CV for both survey indices in the feedback period
index_Neff = input$data$index_Neff
index_Neff[21:50,] = 30 # Decrease ESS for both survey indices in the feedback period
remove_agg = TRUE # remove a aggregate index for some years
remove_agg_pointer = 1 # target on index 1
remove_agg_years = 21:50 # all feedback years
remove_paa = TRUE # Also remove age comp for that index
remove_paa_pointer = 1 # target on index 1
remove_paa_years = 21:50 # all feedback years
input <- update_input_index_info(input, agg_index_sigma, index_Neff,
remove_agg, remove_agg_pointer, remove_agg_years,
remove_paa, remove_paa_pointer, remove_paa_years) # Update input file
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)
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 = "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), # Must have this!
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 analysispar(mfrow = c(1,2))
SSB_s1_m1 <- mod1$om$rep$SSB[,1] # Stock 1 SSB in model 1
plot(SSB_s1_m1, col = "red", type = "o", ylab = "SSB", xlab = "Year")
SSB_s1_m2 <- mod2$om$rep$SSB[,1] # Stock 1 SSB in model 2
plot(SSB_s1_m2, col = "blue", type = "o", ylab = "SSB", xlab = "Year")