This vignette introduces how remove one region (e.g. closed area, windfarm area, refuge habitat) in the assessment:

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 = 3,
                            n_regions = 3, # 3 regions which contains 1 closed area
                            n_indices = 3, # Assume survey is operated in closed area
                            n_fleets = 2, # closed area doesn't have fishery
                            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),
                            fleet_regions = c(1,2),
                            index_regions = c(1,2,3), # Regions 1-3 all have survey
                            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. 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)

7. Generate dataset

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

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

9. Conduct MSE with one region (e.g., closed area) and corresponding data removed.

Since one region is removed, we need to re-configure some processes:

new_n_stocks = 2
new_n_regions = 2
new_n_indices = 2

NAA_re_em <- list(N1_model=rep(ini.opt,new_n_stocks),
                  sigma=rep(sigma,new_n_stocks),
                  cor=rep(re_cor,new_n_stocks),
                  recruit_model = 2) 

# M Configuration
M_em <- list(model="constant",initial_means=array(0.2, dim = c(new_n_stocks,new_n_regions,n_ages)))

# Selectivity Configuration with one survey removed
sel_em <- list(model=rep("logistic",n_fleets+new_n_indices),
            initial_pars=c(rep(list(fleet_pars),n_fleets),rep(list(index_pars),new_n_indices)))

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,
                       reduce_region_info = list(remove_regions = c(1,1,0), # remove the 3rd region
                                                 NAA_re_em = NAA_re_em, # new NAA configuration
                                                 sel_em = sel_em, # new selectivity configuration
                                                 M_em = M_em), # new M configuration
                       seed = 123,
                       save.sdrep = TRUE,
                       save.last.em = TRUE)

mod1$em_input$data$agg_indices # double check if a index is removed from the EM
mod1$em_input$data$index_regions # double check if 2 indices in region 1 and 2

9. Conduct MSE with one region (e.g., closed area) removed but the corresponding data (e.g. survey data) re-assigned to one of the remaining regions.

Since one region is removed, we still need to re-configure some processes:

new_n_stocks = 2
new_n_regions = 2
new_n_indices = 2

NAA_re_em <- list(N1_model=rep(ini.opt,new_n_stocks),
                  sigma=rep(sigma,new_n_stocks),
                  cor=rep(re_cor,new_n_stocks),
                  recruit_model = 2) 

# M Configuration
M_em <- list(model="constant",initial_means=array(0.2, dim = c(new_n_stocks,new_n_regions,n_ages)))

# Selectivity Configuration should not change because survey from region 3 is re-assigned to region 1
sel_em <- list(model=rep("logistic",n_fleets+n_indices), 
            initial_pars=c(rep(list(fleet_pars),n_fleets),rep(list(index_pars),n_indices)))

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,
                       reduce_region_info = list(remove_regions = c(1,1,0), # remove the 3rd region
                                                 reassign = 1, # reassign survey data to region 1
                                                 NAA_re_em = NAA_re_em, # new NAA configuration
                                                 sel_em = sel_em, # new selectivity configuration
                                                 M_em = M_em), # new M configuration
                       seed = 123,
                       save.sdrep = TRUE,
                       save.last.em = TRUE)

mod1$em_input$data$agg_indices # double check if the index from region 3 is kept
mod1$em_input$data$index_regions # double check if index from region 3 is reassigned to region 1