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.

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 information

3. 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. Increase Observation Error and Reduce Sample Size for Survey Indices in the Feedback Period

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 file

7. 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

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. Run 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), # 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 analysis

11. A Worse Case Scenario (No Data for Survey 1 in the Future)

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

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 analysis
par(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")