Check endyr before you start

The OM and the EM must share a terminal year and a projection horizon; the MSE fills catch by row position across both models, so mismatched years are refused rather than run.

Then check endyr against the catch series. The MSE takes over at endyr + 1, and every year after that is a projection year whose catch the control rule sets, so catch recorded past endyr is blanked to NA at setup and run_mse() warns and names the years.

That warning almost always means endyr has fallen behind. A workbook can hold catch, index and composition data for years the model does not fit; the likelihood scores only Year <= endyr, so the data are added and endyr is left where it was. Conditioning the assessment on those years and projecting over them are different questions, and the warning is where you choose. Move endyr up and refit both models to condition on them; leave it and the MSE simulates those years instead.

mse1 <- run_mse(
  om = om, em = em,
  nsim = 10,
  assessment_period = 2,
  # Fishery = fleet 1 (annual), BTS = fleet 2 (biannual)
  sampling_period = c(1, 2),
  seed = 666 # default; per-sim seed = seed + sim
)

Reproducibility

run_mse() accepts two seed arguments, both with default seed = 666:

• seed controls the random draws inside each simulated trajectory (recruitment, observation error, composition resampling). Each simulation worker calls set.seed(seed + sim), so simulation 1 uses seed 667, simulation 2 uses 668, and so on. Re-running the same run_mse() call with the same seed reproduces the trajectory bit-for-bit.
• regenerate_seed (defaults to seed) controls the optional regenerate_past = TRUE step that re-fits the EM to OM-simulated historical data. Set this separately if you want to vary past-data realizations while holding the projection draws fixed.

You do not need to call set.seed() before run_mse(), the seed arguments fully control reproducibility. (You may still want set.seed() for plotting code that draws random colours or jitter.)

Reproducibility holds within a package version, not across one. In 5.9.0 the observation draws moved from R into the TMB model, which displaces the random-number stream, so a seeded run_mse() or self_test() will not reproduce a run from 5.8.x. 5.18.0 is the same kind of boundary: it seeds each assessment’s draw on its own year and stops shortening the operating model’s projection horizon. Both move the stream for every seeded run_mse(simulate_data = TRUE).

5.13.0 is a boundary for two kinds of run rather than all of them: one using a scalar cap, which that release rewrote, and one whose control rule reads reference points from an estimated stock-recruit curve, which it corrected. A run with neither reproduces across 5.13.0. simulate_data = FALSE is unaffected by all three.

Two further points when comparing runs. The stream position is data-dependent: a covariance survey fleet redraws until every row is positive, and a multinomial composition draw stops early once the sample is exhausted, so changing an unrelated covariance or sample size moves everything drawn after it. And any model with an estimated suitability now varies its stomach contents, which it previously did not.

For management-facing runs, record packageVersion("Rceattle") and the seed alongside the saved .rds so a later reviewer can rerun the exact same trajectories from the same package tag.

Different reference points

Here the EM uses a more aggressive rule (F30% target, F10% limit, and Plimit = 0, i.e. no SSB cutoff) to contrast with the Tier 3 run above.

em2 <- Rceattle::fit_mod(
  data_list = nrdata,
  estimateMode = 0, # Estimate
  initMode = 1,      # Assume unfished equilibrium
  M1Fun = build_M1(
    M1_model = "sex_age_invariant",
    linkages = list(
      M1 = linkage_spec(
        formula = ~ 1,
        priors = list(
          "(Intercept)" = lognormal(log(0.06), 0.05)
        )
      )
    )
  ),
  HCR = build_hcr(HCR = "NPFMC", # Tier3 HCR
                  Ftarget = 0.30,# F30%
                  Flimit = 0.10, # F10%
                  Plimit = 0,    # No limit
                  Alpha = 0.1),
  fit_control = fit_control(verbose = 1, phase = TRUE)
)

mse2 <- run_mse(
  om = om, em = em2, 
  nsim = 10, 
  assessment_period = 2,
  sampling_period = c(1, 2) # Fishery = fleet 1, BTS = fleet 2
)

Triennial bottom-trawl survey

sampling_period = c(1, 3) samples the fishery every year but the bottom-trawl survey (BTS) only every third year, testing the cost of a sparser survey.

mse3 <- run_mse(
  om = om, em = em, 
  nsim = 10, 
  assessment_period = 2,
  sampling_period = c(1, 3) # Fishery = fleet 1, BTS = fleet 2
)

4-year assessment cycle and reduced survey coverage

Here the EM is re-fitted only every fourth year (assessment_period = 4) and fut_sample = 0.5 halves future sampling effort on both data streams: composition sample sizes are halved and the index log-SD is doubled (Log_sd * 1 / fut_sample). A lower-cost monitoring scenario.

mse4 <- run_mse(
  om = om, em = em, 
  nsim = 10, 
  fut_sample = 0.5,
  assessment_period = 4,
  sampling_period = c(1, 2) # Fishery = fleet 1, BTS = fleet 2
)

A missed assessment, and a buffer while advice is stale

The scenarios above all set a period: an assessment every second or fourth year, for as long as the projection runs. A skipped assessment is a different thing (one gap in an otherwise regular cycle) and assessment_period also takes the schedule itself, as a vector of the years an assessment is completed.

biennial <- seq(om$data_list$endyr + 2, om$data_list$projyr, by = 2)
miss_yr  <- biennial[ceiling(length(biennial) / 3)]

mse6 <- run_mse(
  om = om, em = em,
  nsim = 10,
  assessment_period = setdiff(biennial, miss_yr), # every year but miss_yr
  sampling_period = c(1, 2)
)

The distinction matters for the answer, not just the bookkeeping. A triennial period is a permanent policy and its cost compounds; one missed assessment is a single episode of stale advice, and how much it costs depends on the state the stock was in when the gap opened.

To pair the gap with a precautionary reduction, give catch_mult a data.frame of Year, Species and mult. It applies only in the pairs listed, so the reduction can be held for the years advice is stale and lifted once the next assessment lands. Species is the species number, and any pair the table omits is multiplied by 1.

Be careful which years those are. The assessment in year Y sets catch for Y+1 onward, so a missed assessment leaves the two years after it on the previous assessment’s advice, not the missed year itself, which the previous assessment sets in the baseline too.

schedule <- setdiff(biennial, miss_yr)

# The years the missed assessment would have set: from the year after it,
# through the next assessment that was completed.
gap_yrs <- (miss_yr + 1):min(schedule[schedule > miss_yr])

# 10% reduction, held for exactly those years.
buffer <- expand.grid(Year = gap_yrs, Species = seq_len(om$data_list$nspp))
buffer$mult <- 0.90

mse7 <- run_mse(
  om = om, em = em,
  nsim = 10,
  assessment_period = schedule,
  sampling_period = c(1, 2),
  catch_mult = buffer
)

Run these against a baseline that differs only in the schedule (same om, same em, same seed) and compare them as paired within-replicate differences.

Two schedules on the same seed share their recruitment deviations, drawn once per replicate before the first assessment. Each assessment’s observation draws are seeded on that assessment’s own year rather than on wherever the previous assessments left the random stream, so an assessment in year Y starts from the same place in every schedule that assesses Y, and a divergence at one assessment does not displace every assessment after it.

Common random numbers are not complete, though, and it is worth knowing where they stop. One sim_mod() call draws every year in the assessment interval under that assessment’s seed, so a year sitting inside a longer interval is drawn under a different seed than the same year in a schedule that assessed it directly. Two schedules that differ in which years they assess realize different observation error in the years between, even where the stock is in the same state. The paired difference is therefore mostly, but not entirely, the schedule, size the replicate count with that in mind rather than assuming all the noise has cancelled.

Once two schedules genuinely diverge (different catch taken, different years surveyed) their draws diverge too. That is a real difference between the runs rather than noise, and no seeding scheme removes it.

Make the last assessment year reach the projection horizon. Catch is filled only up to the last assessment, so a schedule that stops short leaves the trailing years at NA, and mse_summary() summarises over the whole projection, understating average catch, catch IAV and P(Closed) with nothing in the table to say so. A biennial cycle over an odd number of projection years lands here every time: seq(2027, 2050, by = 2) ends at 2049 against a 2050 horizon. run_mse() warns when it happens.

One caveat to take into any write-up: catch_mult multiplies catch, not ABC. A real ABC reduction changes removals only to the extent the fishery attains ABC, and for a stock whose catch sits well below ABC it changes them hardly at all. Either scale the multiplier by recent attainment, 1 - (1 - mult) * attainment, or report the unscaled result as an upper bound on the effect of the reduction.

Climate-linked MSE

To test a strategy under environmental change, link a process in the OM to a covariate. Here recruitment (via R0) tracks a steadily rising environmental index, so the OM’s productivity shifts over time while the EM keeps assuming stationary recruitment.

# Generate random env data
projyrs <- nrdata$styr:nrdata$projyr
projnyrs <- length(projyrs)

# Env data
# A fixed-effect covariate must cover every model year, so build the table
# over styr:projyr rather than merging into whatever years are already there
# (merge() inner-joins, which drops the rest and leaves EnvIndex NA).
nrdata$env_data <- data.frame(Year = projyrs,
                              EnvIndex = seq(0, 1, length.out = projnyrs))

# Fit OM
om2 <- Rceattle::fit_mod(
  data_list = nrdata,
  estimateMode = 0,  # Estimate
  initMode = 1,      # Assume unfished equilibrium
  M1Fun = build_M1(
    M1_model = "sex_age_invariant",
    linkages = list(
      M1 = linkage_spec(
        formula = ~ 1,
        priors = list(
          "(Intercept)" = lognormal(log(0.06), 0.05)
        )
      )
    )
  ),
  recFun = build_srr(
    srr_fun = "mean",    # Mean R with covariates
    proj_mean_rec = FALSE,
    linkages = list(
      R0 = linkage_spec(formula = ~ EnvIndex)
    )),
  fit_control = fit_control(phase = TRUE)
)

plot_biomass(om2, incl_proj = TRUE)

mse5 <- run_mse(
  om = om2, em = em, 
  nsim = 10, 
  fut_sample = 0.5,
  assessment_period = 4,
  sampling_period = c(1, 2) # Fishery = fleet 1, BTS = fleet 2
)