Environmental Consideration in MSE

This vignette demonstrates an end-to-end SPASAM.MSE workflow that includes an environmental covariate (Ecov) affecting recruitment with a lag-1 controlling effect. The Ecov time series is modeled with a latent process (here: AR1), so the “true” Ecov state is not perfectly observed.

A key MSE design choice is what varies across realizations:

Different true Ecov states across realizations
- Each replicate draws a different latent process trajectory (via Ecov_re).
- This represents variability/uncertainty in the underlying environmental forcing itself.
Same true Ecov state across realizations, but different Ecov observations
- All replicates share the same latent trajectory (Ecov_re_fixed), but observed covariate values (Ecov_obs) differ across replicates due to observation error.
- This is often what we want when the “true world” is the same, but measurements are noisy.

This vignette includes BOTH examples side-by-side.

This is a minimal demonstration: 1 stock, 1 region, 1 fleet, 1 index, 1 season. Results are illustrative only and should not be interpreted as a real-world stock assessment.

1. Load “WHAM” and “SPASAM.MSE”:

# Install (if needed):
# install.packages("remotes")
# remotes::install_github("lichengxue/SPASAM.MSE", dependencies = FALSE)
# remotes::install_github("lichengxue/wham", ref = "lab", dependencies = FALSE)

library(wham)
library(SPASAM.MSE)
library(here)

main.dir <- here::here()

2. Generate Basic Information

We define:

Base (burn-in) years: 2001–2020
Feedback years: 3 (so the model years extend to 2023)

n_stocks  <- 1
n_regions <- 1
n_fleets  <- 1
n_indices <- 1
n_seasons <- 1
n_ages    <- 12

year_start <- 2001
year_end   <- 2020
base_years <- year_start:year_end
n_feedback_years <- 3

info <- generate_basic_info(
  n_stocks = n_stocks,
  n_regions = n_regions,
  n_indices = n_indices,
  n_fleets  = n_fleets,
  n_seasons = n_seasons,
  base.years = base_years,
  n_feedback_years = n_feedback_years,
  life_history = "medium",
  n_ages = n_ages
)

basic_info <- info$basic_info
catch_info <- info$catch_info
index_info <- info$index_info
F_info     <- info$F

basic_info$years

3. Configure Numbers-at-Age (NAA)

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

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

NAA_re <- list(
  N1_model      = rep(ini.opt, n_stocks),
  recruit_pars  = 10,   # log scale mean recruitment
  sigma         = rep(sigma, n_stocks),
  cor           = rep(re_cor, n_stocks),
  recruit_model = 2,
  sigma_vals    = sigma_vals
)

4. Configure Environmental Covariate (Ecov)

4.1 Ecov observations and future assumptions

The model years are 2001–2023. We provide Ecov observations for 2000–2020 and extend the observations to 2023 by holding them at the historical mean (users can choose a different projection assumption).

n_ecov <- 1

years_ecov_obs <- 2000:2020
temp_vals <- c(
  8.172843, 7.305779, 7.860323, 6.728931, 6.062992,
  6.329001, 7.261095, 6.246041, 6.702910, 7.058778,
  8.338381, 7.509824, 8.949135, 8.299672, 7.942865,
  8.021900, 8.890488, 7.485228, 8.279137, 7.487395,
  7.774417
)

years_ecov_full <- 2000:2023
temp_mean_val   <- mean(temp_vals)
temp_future_vals <- rep(temp_mean_val, 3)  # 2021–2023
temp_vals_full <- c(temp_vals, temp_future_vals)

temp_mat <- matrix(temp_vals_full, ncol = 1)
colnames(temp_mat) <- "bt_temp"
rownames(temp_mat) <- years_ecov_full

4.2 Recruitment effect and Ecov process parameters

We use a lag-1 controlling effect on recruitment and an AR1 latent process.

Recruitment effect (lag-1 controlling)

The configuration "controlling-lag-1-linear" means the covariate affects recruitment in the next year (lag = 1) and the effect is linear on the recruitment scale used internally by WHAM.

A simple way to think about the effect is:

\[ \log(R_y) = \ldots + \beta \cdot Ecov_{y-1} \]

where \(\beta = -0.3969783\) in this example (a negative effect).

Ecov_beta_R <- array(
  -0.3969783,
  dim = c(1, 1, 1),
  dimnames = list(
    stock = "stock1",
    ecov  = "bt_temp",
    poly  = "linear"
  )
)

Ecov_process_pars <- matrix(
  c(
    7.3656398,   # mu (mean state)
    -0.5825355,  # log_sigma (process innovation SD)
    1.2154606    # transformed phi
  ),
  nrow = 3,
  ncol = 1,
  dimnames = list(c("mu", "log_sigma", "phi"), "bt_temp")
)

4.3 Build the ecov list

ecov <- list(
  label = "bt_temp",
  year  = years_ecov_full,
  mean  = temp_mat,
  logsigma = "est_1",
  use_obs = matrix(1, nrow = length(years_ecov_full), ncol = 1),
  process_model = "ar1",
  recruitment_how = matrix(
    "controlling-lag-1-linear",
    nrow = n_ecov,
    ncol = n_stocks
  )
)

5. Generate `wham` Input

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
)

# Fix Ecov recruitment effect and process parameters for this demo
input$par$Ecov_beta_R       <- Ecov_beta_R
input$par$Ecov_process_pars <- Ecov_process_pars

# Build OM without fitting
random <- input$random
input$random <- NULL

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

# Keep random-effect structure for update_om_fn()
om$random <- random

6. Example A: Different True States Across Realizations

Here we simulate twice with different seeds. We do not modify Ecov_re, so the latent state is allowed to vary across realizations (replicates).

om_A1 <- update_om_fn(om, seed = 123, random = random)
om_A2 <- update_om_fn(om, seed = 124, random = random)

# Years for plotting (use what exists; otherwise use row index)
yrs_A <- om_A1$years_Ecov
if (is.null(yrs_A) && !is.null(om_A1$input$years_Ecov)) yrs_A <- om_A1$input$years_Ecov
if (is.null(yrs_A)) yrs_A <- seq_len(nrow(om_A1$input$data$Ecov_obs))

# Observations differ across seeds
plot(yrs_A, as.numeric(om_A1$input$data$Ecov_obs[,1]),
     type = "l", lwd = 2, col = "blue",
     xlab = "Year", ylab = "Ecov_obs",
     main = "Example A: Ecov observations differ across seeds")
lines(yrs_A, as.numeric(om_A2$input$data$Ecov_obs[,1]),
      lwd = 2, col = "red")

legend("topleft",
       legend = c("Obs (seed 123)", "Obs (seed 124)"),
       col = c("blue", "red"), lwd = 2, bty = "n")

# True latent state also differs across seeds
plot(yrs_A, as.numeric(om_A1$rep$Ecov_x),
     type = "l", lwd = 2, col = "blue",
     xlab = "Year", ylab = "True Ecov state (Ecov_x)",
     main = "Example A: True Ecov state differs across seeds")
lines(yrs_A, as.numeric(om_A2$rep$Ecov_x),
      lwd = 2, col = "red")

legend("topleft",
       legend = c("True state (seed 123)", "True state (seed 124)"),
       col = c("blue", "red"), lwd = 2, bty = "n")

7. Example B: Same True State Across Realizations

In this example, we fix the latent deviations (Ecov_re) to a specified vector. Then we:

Turn OFF Ecov latent simulation (do_simulate_Ecov_re <- 0L)
Remove "Ecov_re" from the random effects list used for simulation
Build an unfitted OM object and simulate two realizations with different seeds

This guarantees:

the true state (Ecov_x) is identical across realizations,
observations (Ecov_obs) can still differ across realizations due to observation error.

# Start from the original WHAM input again (so we don't carry modifications from earlier)
inputB <- prepare_wham_input(
  basic_info = basic_info,
  NAA_re     = NAA_re,
  ecov       = ecov,
  catch_info = catch_info,
  index_info = index_info,
  F          = F_info
)

inputB$par$Ecov_beta_R       <- Ecov_beta_R
inputB$par$Ecov_process_pars <- Ecov_process_pars

# Fixed latent deviations (example values)
Ecov_re_vals <- c(
  0.006585013,  0.382720922, -0.584829153, -1.209735900,
  -0.998780628, -0.222456716, -0.990273119, -0.634984723,
  -0.232817773,  0.842631681,  0.266049570,  1.420400325,
  0.923464230,  0.612037739,  0.710521898,  1.382505929,
  0.238340928,  0.771175504,  0.182324835,  0.461804389,
  1.285226416, 1, 1, 1 # we assume future temperature is 1 degree higher than the mean
)

# Put the fixed Ecov_re into the input parameter list
# (this assumes inputB$par$Ecov_re has matching n_years_Ecov x n_ecov dims)
inputB$par$Ecov_re[,1] <- Ecov_re_vals

# Tell WHAM NOT to simulate Ecov_re (use what we set above)
inputB$data$do_simulate_Ecov_re <- 0L

# Remove "Ecov_re" from random effects so update_om_fn won't simulate it
randomB <- inputB$random
randomB <- randomB[!randomB %in% "Ecov_re"]

# Build OM without fitting
inputB$random <- NULL
omB <- fit_wham(inputB, do.fit = FALSE, do.brps = FALSE, MakeADFun.silent = TRUE)

# Simulate two realizations (same true state, potentially different observations)
om_B1 <- update_om_fn(omB, seed = 123, random = randomB)
om_B2 <- update_om_fn(omB, seed = 124, random = randomB)

# Years for plotting
yrs_B <- om_B1$years_Ecov
if (is.null(yrs_B) && !is.null(om_B1$input$years_Ecov)) yrs_B <- om_B1$input$years_Ecov
if (is.null(yrs_B)) yrs_B <- seq_len(nrow(om_B1$input$data$Ecov_obs))

# Observations differ across seeds (if observation error is simulated)
plot(yrs_B, as.numeric(om_B1$input$data$Ecov_obs[,1]),
     type = "l", lwd = 2, col = "blue",
     xlab = "Year", ylab = "Ecov_obs",
     main = "Example B: Ecov observations differ, fixed true state")
lines(yrs_B, as.numeric(om_B2$input$data$Ecov_obs[,1]),
      lwd = 2, col = "red")

legend("topleft",
       legend = c("Obs (seed 123)", "Obs (seed 124)"),
       col = c("blue", "red"), lwd = 2, bty = "n")

# True latent state should be identical across seeds
plot(yrs_B, as.numeric(om_B1$rep$Ecov_x),
     type = "l", lwd = 2, col = "blue",
     xlab = "Year", ylab = "True Ecov state (Ecov_x)",
     main = "Example B: True Ecov state is fixed across seeds")
lines(yrs_B, as.numeric(om_B2$rep$Ecov_x),
      lwd = 2, col = "red")

legend("topleft",
       legend = c("True state (seed 123)", "True state (seed 124)"),
       col = c("blue", "red"), lwd = 2, bty = "n")

8. Run MSE

Example B is often preferable for MSE when you want the same “true world” but noisy measurements.

We run a single assessment at the terminal year (2020) to provide advice for 2021–2023 (assess_interval = 3).

assess_interval <- 3
terminal_year   <- max(base_years)
assess_years <- terminal_year

# -------------------------
# EM 1: WITH Ecov in the EM
# -------------------------

mod1 <- loop_through_fn(
  om = om_B1,
  em_info = info,
  random = randomB,

  NAA_re_em   = NAA_re,
  ecov_em     = ecov,
  age_comp_em = "multinomial",

  em.opt = list(
    separate.em = FALSE,
    separate.em.type = NULL,
    do.move = FALSE,
    est.move = FALSE
  ),

  assess_years = assess_years,
  assess_interval = assess_interval,
  base_years = base_years,

  year.use = terminal_year,
  add.years = TRUE,
  seed = 123
)

mod1$converge_list # 2 = converged 

# ----------------------------
# EM 2: WITHOUT Ecov in the EM
# ----------------------------

mod2 <- loop_through_fn(
  om = om_B1,
  em_info = info,
  random = randomB,

  NAA_re_em   = NAA_re,
  # ecov_em   = ecov,   # <-- OMIT Ecov in the EM
  age_comp_em = "multinomial",

  em.opt = list(
    separate.em = FALSE,
    separate.em.type = NULL,
    do.move = FALSE,
    est.move = FALSE
  ),

  assess_years = assess_years,
  assess_interval = assess_interval,
  base_years = base_years,

  year.use = terminal_year,
  add.years = TRUE,
  seed = 123
)

mod2$converge_list # 2 = converged

9. Plot MSE Results

# 10. Plot MSE results: EM with Ecov vs. without Ecov

# Years corresponding to indices 20:23
yrs <- basic_info$years[20:23]

# ------------------------
# Spawning Stock Biomass
# ------------------------
plot(
  yrs,
  mod1$om$rep$SSB[20:23],
  type = "l",
  lwd  = 2,
  col  = "red",
  xlab = "Year",
  ylab = "SSB",
  main = "MSE comparison: SSB",
  xaxt = "n"          # turn off default x-axis
)

lines(
  yrs,
  mod2$om$rep$SSB[20:23],
  lwd = 2,
  col = "blue",
  lty = 2
)

axis(1, at = yrs, labels = yrs)  # integer-only x-axis

legend(
  "topleft",
  legend = c("EM with Ecov", "EM without Ecov"),
  col = c("red", "blue"),
  lty = c(1, 2),
  lwd = 2,
  bty = "n"
)

# ------------------------
# Predicted Catch
# ------------------------
plot(
  yrs,
  mod1$om$rep$pred_catch[20:23],
  type = "l",
  lwd  = 2,
  col  = "red",
  xlab = "Year",
  ylab = "Predicted Catch",
  main = "MSE comparison: Catch",
  xaxt = "n"          # turn off default x-axis
)

lines(
  yrs,
  mod2$om$rep$pred_catch[20:23],
  lwd = 2,
  col = "blue",
  lty = 2
)

axis(1, at = yrs, labels = yrs)  # integer-only x-axis

legend(
  "topleft",
  legend = c("EM with Ecov", "EM without Ecov"),
  col = c("red", "blue"),
  lty = c(1, 2),
  lwd = 2,
  bty = "n"
)