This vignette introduces options for catch apportionment when the operating model is spatially explicit while a panmictic model is used as the assessment model.

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 = 2,
                            n_regions = 2,
                            n_indices = 4,
                            n_fleets = 4,
                            n_seasons = 4,
                            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

3. Specify movement type and movement rate

basic_info <- generate_NAA_where(basic_info = basic_info, move.type = 2) # "bidirectional" movement

move <- generate_move(basic_info = basic_info, move.type = 2, move.rate = c(0.2,0.1), move.re = "constant")

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

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

6. 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, 
                            move = move,
                            catch_info = catch_info, 
                            index_info = index_info, 
                            F = F_info)

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 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. Using panmictic assessment model

There are several options for catch apportionment:

Option 1 (weight_type = 1 with method = “equal”): Assign equal weights to all regions.

mod = list()

n_stocks = n_regions = 1

n_fleets = n_indices = 2

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

NAA_re_em <- list(N1_model="equilibrium",sigma="rec+1",cor="iid")

M_em <- list(model="constant",initial_means=array(0.2, dim = c(n_stocks,n_regions,n_ages)))

mod[[1]] = loop_through_fn(om = om_with_data,
                           em_info = info, 
                           random = random,
                           M_em = M_em, 
                           sel_em = sel_em,
                           NAA_re_em = NAA_re_em, 
                           move_em = NULL,
                           age_comp_em = "multinomial",
                           em.opt = list(separate.em = TRUE, separate.em.type = 1, 
                                         do.move = FALSE, est.move = FALSE),
                           aggregate_catch_info = list(n_fleets = 2, # 2 aggregated fleets
                                                       fleet_pointer = c(1,2,1,2), 
                                                       # fleet 1 (in region 1) and 3 (in region 2) will be aggregated, 
                                                       # fleet 2 (in region 1) and 4 (in region 2) will be aggregated 
                                                       use_catch_weighted_waa=TRUE, 
                                                       catch_Neff = c(100,100), # use ESS = 100 for aggregate fleet
                                                       catch_cv = c(0.1,0.1)), # use CV = 0.1 for aggregate fleet
                           aggregate_index_info = list(n_indices = 2, # 2 aggregated indices
                                                       index_pointer = c(1,2,1,2), 
                                                       # index 1 (in region 1) and 3 (in region 2) will be aggregated, 
                                                       # index 2 (in region 1) and 4 (in region 2) will be aggregated 
                                                       use_catch_weighted_waa=TRUE,
                                                       index_Neff = c(100,100), # use ESS = 100 for aggregate index
                                                       index_cv = c(0.1,0.1)), # use CV = 0.1 for aggregate index
                           catch_alloc = list(weight_type = 1, method = "equal", weight_years = 3),
                           assess_years = assess.years, 
                           assess_interval = assess.interval, 
                           base_years = base.years,
                           year.use = 20,
                           seed = 123,
                           save.sdrep = TRUE,
                           save.last.em = TRUE
                           )

Option 2 (weight_type = 2 with method = “fleet_region”): Use regional catch totals to compute weights.

mod[[2]] = loop_through_fn(om = om_with_data,
                           em_info = info, 
                           random = random,
                           M_em = M_em, 
                           sel_em = sel_em,
                           NAA_re_em = NAA_re_em, 
                           move_em = NULL,
                           age_comp_em = "multinomial",
                           em.opt = list(separate.em = TRUE, separate.em.type = 1, 
                                         do.move = FALSE, est.move = FALSE),
                           aggregate_catch_info = list(n_fleets = 2, # 2 aggregated fleets
                                                       fleet_pointer = c(1,2,1,2), 
                                                       # fleet 1 (in region 1) and 3 (in region 2) will be aggregated, 
                                                       # fleet 2 (in region 1) and 4 (in region 2) will be aggregated 
                                                       use_catch_weighted_waa=TRUE, 
                                                       catch_Neff = c(100,100), # use ESS = 100 for aggregate fleet
                                                       catch_cv = c(0.1,0.1)), # use CV = 0.1 for aggregate fleet
                           aggregate_index_info = list(n_indices = 2, # 2 aggregated indices
                                                       index_pointer = c(1,2,1,2), 
                                                       # index 1 (in region 1) and 3 (in region 2) will be aggregated, 
                                                       # index 2 (in region 1) and 4 (in region 2) will be aggregated 
                                                       use_catch_weighted_waa=TRUE,
                                                       index_Neff = c(100,100), # use ESS = 100 for aggregate index
                                                       index_cv = c(0.1,0.1)), # use CV = 0.1 for aggregate index
                           catch_alloc = list(weight_type = 2, method = "fleet_region", weight_years = 3),
                           assess_years = assess.years, 
                           assess_interval = assess.interval, 
                           base_years = base.years,
                           year.use = 20,
                           seed = 123,
                           save.sdrep = TRUE,
                           save.last.em = TRUE) 

Option 3 (weight_type = 3, with method = “index_equal”): Use survey-based regional weighting, then assigns equally among fleets in the same region

mod[[3]] = loop_through_fn(om = om_with_data,
                           em_info = info, 
                           random = random,
                           M_em = M_em, 
                           sel_em = sel_em,
                           NAA_re_em = NAA_re_em, 
                           move_em = NULL,
                           age_comp_em = "multinomial",
                           em.opt = list(separate.em = TRUE, separate.em.type = 1, 
                                         do.move = FALSE, est.move = FALSE),
                           aggregate_catch_info = list(n_fleets = 2, # 2 aggregated fleets
                                                       fleet_pointer = c(1,2,1,2), 
                                                       # fleet 1 (in region 1) and 3 (in region 2) will be aggregated, 
                                                       # fleet 2 (in region 1) and 4 (in region 2) will be aggregated 
                                                       use_catch_weighted_waa=TRUE, 
                                                       catch_Neff = c(100,100), # use ESS = 100 for aggregate fleet
                                                       catch_cv = c(0.1,0.1)), # use CV = 0.1 for aggregate fleet
                           aggregate_index_info = list(n_indices = 2, # 2 aggregated indices
                                                       index_pointer = c(1,2,1,2), 
                                                       # index 1 (in region 1) and 3 (in region 2) will be aggregated, 
                                                       # index 2 (in region 1) and 4 (in region 2) will be aggregated 
                                                       use_catch_weighted_waa=TRUE,
                                                       index_Neff = c(100,100), # use ESS = 100 for aggregate index
                                                       index_cv = c(0.1,0.1)), # use CV = 0.1 for aggregate index
                           catch_alloc = list(weight_type = 3, method = "index_equal", weight_years = 3, survey_pointer = 1), # use index 1
                           assess_years = assess.years, 
                           assess_interval = assess.interval, 
                           base_years = base.years,
                           year.use = 20,
                           seed = 123,
                           save.sdrep = TRUE,
                           save.last.em = TRUE)

Option 4 (weight_type = 4 with method = “user_defined_fleets”): Use manually specified weights for each region or each fleet

mod[[4]] = loop_through_fn(om = om_with_data,
                           em_info = info, 
                           random = random,
                           M_em = M_em, 
                           sel_em = sel_em,
                           NAA_re_em = NAA_re_em, 
                           move_em = NULL,
                           age_comp_em = "multinomial",
                           em.opt = list(separate.em = TRUE, separate.em.type = 1, 
                                         do.move = FALSE, est.move = FALSE),
                           aggregate_catch_info = list(n_fleets = 2, # 2 aggregated fleets
                                                       fleet_pointer = c(1,2,1,2), 
                                                       # fleet 1 (in region 1) and 3 (in region 2) will be aggregated, 
                                                       # fleet 2 (in region 1) and 4 (in region 2) will be aggregated 
                                                       use_catch_weighted_waa=TRUE, 
                                                       catch_Neff = c(100,100), # use ESS = 100 for aggregate fleet
                                                       catch_cv = c(0.1,0.1)), # use CV = 0.1 for aggregate fleet
                           aggregate_index_info = list(n_indices = 2, # 2 aggregated indices
                                                       index_pointer = c(1,2,1,2), 
                                                       # index 1 (in region 1) and 3 (in region 2) will be aggregated, 
                                                       # index 2 (in region 1) and 4 (in region 2) will be aggregated 
                                                       use_catch_weighted_waa=TRUE,
                                                       index_Neff = c(100,100), # use ESS = 100 for aggregate index
                                                       index_cv = c(0.1,0.1)), # use CV = 0.1 for aggregate index
                           catch_alloc = list(weight_type = 4, method = "user_defined_fleets", user_weights = c(0.4,0.4,0.1,0.1)), # user-specified
                           assess_years = assess.years, 
                           assess_interval = assess.interval, 
                           base_years = base.years,
                           year.use = 20,
                           seed = 123,
                           save.sdrep = TRUE,
                           save.last.em = TRUE) # average over the most recent 5 years

Option 5 (weight_type = 4): Weights are calculated using regional recruitment (which is equal in this case).

mod[[5]] = loop_through_fn(om = om_with_data,
                           em_info = info, 
                           random = random,
                           M_em = M_em, 
                           sel_em = sel_em,
                           NAA_re_em = NAA_re_em, 
                           move_em = NULL,
                           age_comp_em = "multinomial",
                           em.opt = list(separate.em = TRUE, separate.em.type = 1, 
                                         do.move = FALSE, est.move = FALSE),
                           aggregate_catch_info = list(n_fleets = 2, # 2 aggregated fleets
                                                       fleet_pointer = c(1,2,1,2), 
                                                       # fleet 1 (in region 1) and 3 (in region 2) will be aggregated, 
                                                       # fleet 2 (in region 1) and 4 (in region 2) will be aggregated 
                                                       use_catch_weighted_waa=TRUE, 
                                                       catch_Neff = c(100,100), # use ESS = 100 for aggregate fleet
                                                       catch_cv = c(0.1,0.1)), # use CV = 0.1 for aggregate fleet
                           aggregate_index_info = list(n_indices = 2, # 2 aggregated indices
                                                       index_pointer = c(1,2,1,2), 
                                                       # index 1 (in region 1) and 3 (in region 2) will be aggregated, 
                                                       # index 2 (in region 1) and 4 (in region 2) will be aggregated 
                                                       use_catch_weighted_waa=TRUE,
                                                       index_Neff = c(100,100), # use ESS = 100 for aggregate index
                                                       index_cv = c(0.1,0.1)), # use CV = 0.1 for aggregate index
                           catch_alloc = list(weight_type = 4, method = "user_defined_regions", user_weights = c(0.5,0.5)), # user-specified
                           assess_years = assess.years, 
                           assess_interval = assess.interval, 
                           base_years = base.years,
                           year.use = 20,
                           seed = 123,
                           save.sdrep = TRUE,
                           save.last.em = TRUE) # average over the most recent 5 years

Plot catch

fleet1catch_m1 <- mod[[1]]$om$rep$pred_catch[,1] # fleet1catch in model 1
fleet1catch_m2 <- mod[[2]]$om$rep$pred_catch[,1] # fleet1catch in model 2
fleet1catch_m3 <- mod[[3]]$om$rep$pred_catch[,1] # fleet1catch in model 3
fleet1catch_m4 <- mod[[4]]$om$rep$pred_catch[,1] # fleet1catch in model 4
fleet1catch_m5 <- mod[[5]]$om$rep$pred_catch[,1] # fleet1catch in model 5

plot(fleet1catch_m1, col = "red", type = "o", ylab = "Catch", xlab = "Year", main = "Fleet 1 Catch", ylim = c(0, 5000))
lines(fleet1catch_m2, col = "blue", type = "o")
lines(fleet1catch_m3, col = "green", type = "o")
lines(fleet1catch_m4, col = "purple", type = "o")
lines(fleet1catch_m5, col = "orange", type = "o")

legend("topright", legend=c("M 1", "M 2", "M 3", "M 4", "M 5"),
       col=c("red", "blue", "green", "purple", "orange"), lty=c(1,1), cex=0.8)
Figure 1
Figure 1