This vignette demonstrates how to use the whamMSE (Woods Hole Assessment Model Management Strategy Evaluation) package to incorporate an environmental covariate 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   <- 3     # 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 = 1,
                            n_fleets = 1,
                            n_seasons = 1,
                            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. 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!

4. Environmental Covariate on Recruitment

In wham, environmental covariates can be incorporated into various processes (e.g., natural mortality, recruitment, movement, and survey catchability), and can be applied in the Operating Model (OM), Estimation Model (EM), as well as during projection. See the wham vignettes for more details.

Here, we demonstrate how to incorporate an environmental covariate into recruitment in the OM, which controls fish population dynamics. The same environmental covariate is also included in the EM to ensure consistency during the feedback period. Note that environmental covariates can also be applied in projection, but we do not use that option in this example.

# Create a pseudo ecov time series
years = basic_info$years
env.dat <- data.frame(Year = years, mean = runif(length(years), -1, 1), sigma = runif(length(years), 0, 0.2))

ecov <- list(
    label = "Ecov_x",
    mean = as.matrix(env.dat$mean),
    logsigma = as.matrix(log(env.dat$sigma)),
    year = env.dat$Year,
    use_obs = matrix(1, ncol=1, nrow=dim(env.dat)[1]), # use all obs (=1)
    process_model = "rw", # "rw" or "ar1"
    recruitment_how = matrix("controlling-lag-1-linear",1,1)) #matrix for number of stocks (1) and number of Ecovs (1)

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, 
                            NAA_re = NAA_re,
                            ecov = ecov,
                            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

assess.interval <- 3
base.years      <- year_start:year_end
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

During the MSE feedback period, the environmental covariate time series will be automatically truncated to ensure that each assessment model only uses covariate data available up to the most recent assessment year. This prevents the EM from accessing future covariate values beyond the current assessment year, even if the full time series was specified in the OM.

mod = loop_through_fn(om = om_with_data,
                      em_info = info, 
                      random = random,
                      NAA_re_em = NAA_re, 
                      ecov_em = ecov,
                      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,
                      add.years = TRUE)