Generate OM and Perform Self- and Cross-test

This vignette demonstrates how to use the whamMSE (Woods Hole Assessment Model Management Strategy Evaluation) package to generate a multi-stock, multi-region operating model and perform self-tests and cross-tests. 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

The operating model is generated based on user-defined spatial and temporal resolution, historical period, projection period, biological characteristics, and fishery/survey data. The generate_basic_info() function creates a list of biological and fishery details used to build a wham input. Users can choose built-in life history templates or customize biological info from real-world assessments. See ?generate_basic_info for full details.

year_start  <- 1
year_end    <- 20

info <- generate_basic_info(
  n_stocks = 2,
  n_regions = 2,
  n_indices = 2,
  n_fleets = 2,
  n_seasons = 4,
  base.years = year_start:year_end,
  life_history = "medium",
  n_ages = 12,
  fracyr_spawn = 0.625
)

basic_info = info$basic_info
catch_info = info$catch_info
index_info = info$index_info
F_info = info$F

3. Configure Selectivity

By default, selectivity is time-invariant but can also include random effects. See the wham vignettes for more details.

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

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

4. Configure Natural Mortality

Natural mortality is age- and time-invariant by default, but random effects can be specified. See the wham vignettes for details.

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

5. Configure Numbers-at-Age (NAA)

Recruitment varies around a mean by default; alternative stock-recruitment options include “Beverton-Holt” and “Ricker”. See the wham vignettes.

sigma        <- "rec+1"
re_cor       <- "iid"
ini.opt      <- "equilibrium"
Rec_sig      <- 0.2
NAA_sig      <- 0.1

sigma_vals <- array(NAA_sig, dim = c(n_stocks, n_regions, n_ages))
sigma_vals[,,1] <- Rec_sig

NAA_re <- list(
  N1_model = rep(ini.opt, n_stocks),
  sigma = rep(sigma, n_stocks),
  cor = rep(re_cor, n_stocks),
  recruit_model = 2,
  sigma_vals = sigma_vals
)

6. Generate wham Input

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
)

7. Generate Operating Model

Use fit_wham() with do.fit = FALSE to initialize the OM without fitting data.

om <- fit_wham(input, do.fit = FALSE, do.brps = TRUE, MakeADFun.silent = TRUE)

8. Self-test

Create a function to simulate data and run the self-test.

sim_fn <- function(om, self.fit = FALSE) {
  input <- om$input
  input$data <- om$simulate(complete = TRUE)
  if (self.fit) {
    fit <- fit_wham(input, do.osa = FALSE, do.retro = FALSE, MakeADFun.silent = TRUE)
    return(fit)
  } else return(input)
}

set.seed(12345)
self_sim_fit <- sim_fn(om, self.fit = TRUE)

Check convergence:

check_convergence(self_sim_fit)

9. Create HTML Output (Optional)

plot_wham_output(self_sim_fit, out.type = "html")

10. Cross-test

Test the OM against a misspecified EM with different NAA configuration.

sigma <- "rec+1"
re_cor <- "2dar1"
ini.opt <- "equilibrium"

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

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
)

em <- fit_wham(input, do.fit = FALSE, do.brps = FALSE, MakeADFun.silent = TRUE)

Cross-test function:

sim_fn2 <- function(om, em, cross.fit = FALSE) {
  input <- em$input
  input$data <- om$simulate(complete = TRUE)
  if (cross.fit) {
    fit <- fit_wham(input, do.osa = FALSE, do.retro = FALSE, MakeADFun.silent = TRUE)
    return(fit)
  } else return(input)
}

set.seed(123)
cross_sim_fit <- sim_fn2(om, em, cross.fit = TRUE)

Check convergence:

check_convergence(cross_sim_fit)

11. 100 Replicates: Self-tests (Optional)

nsim <- 100
set.seed(8675309)
sim_input <- lapply(1:nsim, function(x) {
  input_i <- om$input
  sim <- om$simulate(complete = TRUE)
  input_i$data <- sim
  return(input_i)
})

sim_fits <- lapply(1:nsim, function(x) {
  cat(paste("model_fit:", x, "start \n"))
  out <- fit_wham(sim_input[[x]], do.osa = FALSE, MakeADFun.silent = TRUE, retro.silent = TRUE, save.sdrep = FALSE)
  cat(paste("model_fit:", x, "done \n"))
  return(out)
})

conv <- sapply(1:nsim, function(x) {
  fit <- sim_fits[[x]]
  if (length(fit) != 0) {
    if (fit$is_sdrep & !fit$na_sdrep & !fit$hessian) TRUE else FALSE
  } else FALSE
})

cat(paste("Convergence rate:", sum(conv)/nsim))

12. 100 Replicates: Cross-tests (Optional)

sim_input <- lapply(1:nsim, function(x) {
  input_i <- em$input
  sim <- om$simulate(complete = TRUE)
  input_i$data <- sim
  return(input_i)
})

cross_sim_fit <- lapply(1:nsim, function(x) {
  cat(paste("model_fit:", x, "start \n"))
  out <- fit_wham(sim_input[[x]], do.osa = FALSE, MakeADFun.silent = TRUE, retro.silent = TRUE, save.sdrep = FALSE)
  cat(paste("model_fit:", x, "done \n"))
  return(out)
})