Composable management procedures and MSE

Code

if (dir.exists(file.path("..", "R")) && requireNamespace("pkgload", quietly = TRUE)) {
  pkgload::load_all("..", quiet = TRUE)
} else {
  library(decamod)
}

Purpose and scope

A management procedure (MP) combines an estimation model (EM), a harvest control rule (HCR), and the operational rules connecting them. The MSE runner evaluates complete procedures against an operating model (OM). It does not need separate SPM and GP versions of its annual loop.

This is an initial, deliberately bounded interface. It provides the existing rock-lobster HCR, SPM specifications, frozen single-year Bayesian GP EMs, and a single-region, two-season decamod OM adapter. It does not replace the existing do_mpe(), spm_model(), or rule_tacc() interfaces. It does not migrate any previous CRA2 report, rerun a full experiment, or make a scientific choice of EM.

Define components, then combine them

The constructors define specifications; they do not run fits or simulations. Numerical HCR parameters and biomass reference years are explicit because these are stock-specific scientific choices, not universal defaults.

Code

reference_years <- 2015:2024
hcr1 <- hcr_lobster(
  parameters = c(4, 0.1, 1, 1.8, 80, 0.25, 0.05, 0.05, 0, 0),
  reference_years = reference_years
)

original_spm <- em_spm("cra2_original", reference_years)
fixed_spm <- em_spm("toy_fixed", reference_years)

mp_original <- mp("original_SPM_HCR1", original_spm, hcr1)
mp_fixed <- mp("fixed_SPM_HCR1", fixed_spm, hcr1)
mp_fixed

Management procedure: fixed_SPM_HCR1
   EM: toy_fixed
   HCR: rock-lobster rule 4
   Target: AW relative_biomass
   Failure policy: hold_tacc ; extrapolation: flag 

The initial SPM configurations are intentionally distinct:

Configuration	Process error	Biomass transition	Estimate acceptance
cra2_original	Estimated; starts at 0.001; original prior retained	Original unguarded model	Finite positive estimate, retaining the original use of non-converged fits
toy_fixed	Fixed at 0.10 by default	Smooth positivity guard and penalty	Convergence code zero, finite objective/estimate, and no final guard activation
toy_recalibrated	As toy_fixed	As toy_fixed	As toy_fixed, with supplied status-preserving recalibration

The original process-error prior has not been changed during this refactor. Keeping an original specification available is not a recommendation to use it. Gradients and boundary flags are recorded but do not silently become new rejection criteria. SPM uncertainty intervals remain unavailable (NA), rather than being filled with invented values.

Recalibration coefficients are supplied from the chosen calibration artefact:

Code

recalibrated <- em_spm(
  "toy_recalibrated", reference_years,
  recalibration = c(intercept = chosen_intercept, slope = chosen_slope),
  recalibration_target = mse_target(reference_years, "end_of_year")
)

gp <- em_gp(
  model = trained_single_year_gp,
  reference_years = reference_years,
  catch_component = "legal",
  id = "single_year_CPUE_catch_GP",
  target = mse_target(reference_years, "end_of_year")
)
mp_gp <- mp("GP_HCR1", gp, hcr1)

The GP uses its stored scaler, fitted coefficients, covariance matrices, and training design. The three existing single-year GP alternatives are different saved models passed to em_gp(). Five-year models and SPM–GP hybrids are not supported by this constructor. Training occurs separately, not each simulated year. The current GP rule uses the posterior median on the original scale; using a conservative quantile would define a different MP.

Biological target and catch-component provenance

The original SPM supplies start-of-AW status in the last observation year. The saved CRA2 GP was trained against end-of-year Bfinal status. Bfinal uses the numbers remaining after the year’s dynamics, with terminal AW weights, selectivity, and legal-size indicators. Calling both simply “AW biomass” hides a real timing difference. The existing SPM recalibration was fitted against the latter target, so the example declares that response explicitly while retaining its native SPM target in native_target.

mse_target() records the definition; it does not convert estimates or retrain the model. A GP trained with target metadata cannot be relabelled at prediction. For legacy saved GPs without metadata, supply the known target explicitly to em_gp(). Otherwise timing is unspecified, and declared-target error summaries are refused. The OM records both ratio_OM_aw and ratio_OM_end, while retaining the old ratio_OM as the end-of-year alias for reproducibility.

Use summarise_mse(result, truth = "declared") for each EM’s stated estimation task, or truth = "ratio_OM_end" to score all signals as proxies for the same management target. The summary labels both the declared and scoring targets. Different timing-specific errors are not an identical-task comparison.

No numerical change to the original SPM, historical GP, or recalibration is implied by this annotation. A future common-target retraining experiment must be a separately versioned scientific decision, not a silent package correction.

Operate an MP on observations

Population-removal diagnostics

For catch-conditioned exploitation fractions, diagnose_exploitation() checks the shared decamod population report. It separates the existing 2% survival guard from genuinely overdrawn occupied cells, combines overlapping fisheries and handling deaths, and quantifies both excess removals and the guard’s adjustment in numbers and tonnes. The calculations use mid-season abundance; legacy catch predictions are converted to that timing for bookkeeping checks. Exact agreement between prescribed and predicted catch is checked separately: it does not prove that the unguarded population could supply those removals.

The om_decamod() adapter records these diagnostics automatically, along with previously/ever-guarded flags for each projection trajectory. Annual checkpoints retain the flags. summarise_mse_diagnostics(result) reports event frequencies, affected trajectories, magnitude maxima, and the number of known observations. Non-applicable or missing diagnostics remain NA, not zero. The same audit is included in validate_mse() and summarise_mse() output.

These are diagnostics only. A guard event does not remove a trajectory, reject an EM estimate, change an HCR decision, or cap catches. The adjustment is not called a realised catch shortfall: allocating a shortfall across fisheries and handling deaths would require a separate, explicit scientific policy. Historical pre-projection guard events are not included in the cumulative projection flag.

Code

audit <- diagnose_exploitation(obj$report(), data$catch_ytrf)
audit$summary
audit$cells[audit$cells$overdrawn %in% TRUE, ]
summarise_mse_diagnostics(result)

Annual observation packet

mse_observations() is the only annual data packet accepted by the built-in EMs. It contains observed catch, CPUE, CPUE uncertainty, and timing. Extra columns are stripped, and future or duplicate observations are rejected. It does not contain true biomass, OM parameters, or future recruitment-scenario labels.

Catch components must be identified explicitly. To reproduce the existing toy experiment, the SPM uses total catch across modelled fisheries, whereas the saved GP models use legal catch (commercial plus recreational in fishery 1). Calling the latter “total catch” without this qualification would obscure the different input definitions. A new training experiment can choose another component, but its prediction inputs must match its training inputs.

Code

observations <- mse_observations(
  catch = observed_seasonal_catch,
  cpue = observed_seasonal_cpue,
  year = 2025
)
decision <- run_mp(mp_gp, observations, tacc_history = c(80, 80))
decision$advice
decision$estimate

The catch table has year, season (1 = AW, 2 = SS), component, and value. The CPUE table has year, season, series, value, and sd (log-scale). The SPM needs complete seasonal catch history and the historical reference period. Single-year GP models need the current AW and SS CPUE and, where used, current seasonal catch. The SPM retains the legacy alternating AW/SS CPUE-series mapping; incompatible observation ordering is rejected explicitly.

The advice year is the next modelled year, after data through the current SS season. This is the current simulation convention. An operational calendar must also allow data processing and regulatory deadlines; this version does not claim that arbitrary lags or assessment frequencies are implemented.

run_mp() returns the advice, advice year, estimate, diagnostics, fallback reason, and updated EM state. Its TACC history includes the current TACC. Latent-year rules require at least two TACC values. The same function can be called on real observations independently of any MSE.

By default, unusable EM output holds the current TACC. GP extrapolation is flagged but does not itself reject advice, preserving the toy experiment. To evaluate an extrapolation fallback, define another MP explicitly:

Code

mp_gp_cautious <- mp("GP_HCR1_extrapolation_fallback", gp, hcr1,
                     out_of_range = "fallback")

These are operational fallback rules, not evidence that holding catch is always biologically safe. They still require closed-loop evaluation.

Evaluate the procedures

The real decamod example below is deliberately not executed when this vignette is built. Model files and trained emulators remain project artefacts, not package data, and constructing them does not launch an expensive MSE.

Code

om <- om_decamod(
  models = c(base = "runs/base.rda", declining = "runs/base.rda"),
  rec_devs = c(base = "Rdev.rda", declining = "Rdev_robustness.rda"),
  initial_tacc = 80,
  catch_rec = c(0.8726319, 7.8536874),
  region = "CRA 2",
  reference_years = reference_years
)

result <- run_mse(
  om = om,
  mps = list(mp_original, mp_fixed, mp_gp),
  scenarios = c("base", "declining"),
  iterations = 901:902,
  years = 2025:2026,
  cores = 1,
  checkpoint_dir = "results/interface_smoke"
)
result$data

Each MP–scenario–iteration combination has its own OM, catch, observation, and EM history. Paired MPs receive the same environmental random seeds, but their realised observations diverge as their decisions alter the stock. Stochastic EMs have a separate random stream. Serial and fork-parallel execution preserve numerical results, and the caller’s random-number state is restored afterwards.

With neutral defaults, the decamod adapter retains the toy experiment’s seasonal catch implementation, fixed recreational catches, and no TACC implementation error. Other modelled catch components retain their OM-projected values. Explicit robustness configurations can change these assumptions (below). Its posterior pool defaults to the existing 2,000-draw preparation, even for a tiny subset of iterations; reducing that pool can change draw identity and is not merely a speed setting.

Diagnostics, reproducibility, and extension

result$data contains advice and true OM outcomes alongside the EM estimate, objective, maximum gradient, convergence, boundary flags, positivity guards, GP range/distance diagnostics, fallbacks, and timing. Inapplicable diagnostics are NA. Model fitting state is kept separately from the immutable specification.

Annual checkpoints contain plain state and accumulated annual rows. They are written atomically and never overwritten. Resume verifies the MP specifications, OM input identities, iteration/year design, runtime versions, and package R code. Do not resume across a scientific or code change: use a new output directory. max_steps = 1 deliberately pauses after one annual checkpoint; another call with the same design resumes. An OM failure stops the run, keeping the last complete checkpoint, instead of treating an invalid simulated world as an EM failure with an ordinary catch fallback.

Custom EMs can be defined with em_function(), custom HCRs with hcr_function(), and a different OM can supply the three callbacks in mse_om(). All configuration must be explicit; callbacks must not depend on mutable globals or hidden closure data. The callback source and explicit configuration are fingerprinted, not arbitrary external globals. Store serialisable numeric/list state, not live RTMB objects or environments.

The existing CRA2 scripts remain the comparison reference. The accompanying project script verify_decamod_mp_interface.R checks all five toy EMs, one HCR, one iteration, and two years against the original runner, including a checkpoint resume. A successful numerical migration test establishes software equivalence for that design; it does not establish emulator robustness or MP suitability.

Several EMs, several HCRs

Code

ems <- list(original = original_spm, recalibrated = recalibrated, gp = gp)
hcrs <- list(rule_1 = hcr1, rule_2 = another_hcr, rule_3 = third_hcr)
procedures <- mp_grid(ems, hcrs)
result <- run_mse(om, procedures, scenarios = c("base", "declining"),
                  iterations = 901:902, years = 2025:2026,
                  checkpoint_dir = "results/new_design_smoke")
checks <- validate_mse(result)
summaries <- summarise_mse(result, truth = "ratio_OM_end")
contrasts <- compare_mse(summaries, benchmark = "original",
  draw_map = data.frame(Iteration = 901:902, chain = 2L, draw = 401:402),
  block_length = 1, n_boot = 20)

This tiny example checks software, not statistical precision. A full evaluation must supply the actual posterior chain/draw map. compare_mse() resamples whole paired trajectories in circular blocks within chains, never across chain boundaries. It keeps scenario and HCR comparisons separate. Zero paired risk differences do not establish risk equivalence. The function reports available pairs rather than filling undefined statistics (for example, proportional catch changes after zero TACC).

validate_mse() obtains every expected MP, iteration, and year from the saved manifest. Each HCR is replayed against its own realised TACC history, including historical catches for latent-year rules. Failure-policy checks retain fallback years. It does not assume there are two EMs or that any non-SPM is a GP.

Prospective robustness controls

mse_robustness() separates observation error, catch reporting, implementation, and biological recruitment changes. A single configuration can apply to all scenarios, or a named list must match the OM’s scenario names exactly.

Code

neutral <- mse_robustness()
# Illustrative SOFTWARE-TEST values, not selected scientific trial levels:
changed_index <- mse_robustness(cpue_power = 0.8, cpue_drift = 0.02)
unrecognised_noise <- mse_robustness(cpue_sd_multiplier = 1.5,
                                     cpue_reported_sd_multiplier = 1)
underreported <- mse_robustness(catch_reporting_multiplier = 0.9)
implementation <- mse_robustness(implementation_multiplier = 1.05,
                                implementation_log_sd = 0.05)
missing_ss <- mse_robustness(missing_cpue_seasons = 2,
                             missing_cpue_years = 2027)
low_recruitment <- mse_robustness(recruitment_multiplier = 0.8)

For a baseline index simulation $I$ with conditional median $\mu$ and the last historical seasonal observation $a$, projected year $k$ uses

$$I^*=a(\mu/a)^p(1+d k)m\exp\{s\log(I/\mu)\}.$$

Here $p$ is cpue_power, $d$ is cpue_drift, $m$ is cpue_multiplier, and $s$ is cpue_sd_multiplier. This tests distortion of the baseline index relationship. Drift is additional to baseline fitted creep, not an absolute replacement for the historical 1% assumption. A 2% drift setting must not be labelled “true creep = 2%”. Trial magnitudes need a scientific choice.

Implementation noise is mean-one lognormal before applying its explicit bias multiplier. The OM receives actual removals; EMs receive separately reported catch. True and reported CPUE uncertainty can also differ. Historical data and the four-year historical recruitment overlap are unchanged. No trial label, future recruitment, or true biomass enters an EM’s observation packet.

Missing CPUE exercises the stated failure policy; it does not impute data. The current SPM requires complete CPUE history, so an unrepaired missing value can continue to trigger its fallback in later years. The single-year GP can resume once its current-year observations are complete. That is an explicit operational difference, not an estimation-performance improvement to hide. Holding TACC after failure is not automatically biologically safe.

Alternative conditioned biological models remain separate model files supplied to om_decamod(). Arbitrary changes to growth, natural mortality, initial depletion, spatial structure, or the real-world regulatory calendar are not silently introduced by these controls. Data remain available through completed SS, with advice implemented in the following modelled year.

Readiness checks and installation

The development version is 1.0.0.9002. Build/install into a fresh, isolated R library before a study; do not overwrite the frozen 500-iteration installation. The opt-in tools/verify-cra2-mse.R script takes explicit CRA2 input paths and a new output directory. It checks the original SPM against fitting expressions extracted verbatim from MSE_SPM.R, runs nine MPs over a small three-scenario design with checkpoint restart, and checks prospective perturbations. It never launches the 1,000-iteration scientific study or changes frozen input files.

Raw maximum gradients remain available. Active bound sides, signed boundary gradients, and a projected-gradient maximum are additionally saved to assess constrained optima. These diagnostics do not silently change either SPM’s estimate-acceptance policy. OM exploitation-guard counts are also retained.

Single-EM plateau-height calibration

mse_plateau_design() varies parameter 5 only. run_plateau_search() executes the paired design and summarise_plateau_search() reports the response of true end-of-year B/B_R, scenario-specific results, and an explicitly weighted mixture. Only one frozen EM is permitted during calibration. Its HCRs can later be crossed with all candidate EMs without further tuning. Coarse target brackets are not adopted heights. Non-monotonic responses and targets above the zero-height result are flagged rather than hidden by an interpolator.

Code

outputs <- mse_outputs(B_R, "Source, hash, and derivation of the fixed B_R",
  size_breaks = c(seq(30, 90, 2), Inf),
  thresholds = data.frame(name = c("male_ge70", "female_ge84", "female_ge90"),
    sex = c(1, 2, 2), lower = c(70, 84, 90)), number_units = "individuals")
# Supply outputs to om_decamod(..., outputs = outputs).
design <- mse_plateau_design(gp_em, c(0, 40, 80, 120, 160, 200),
  c(4, .1, 1.2, 1.8, 80, .25, .06, .05, 0, 0), 2015:2024)
# No simulation is launched by constructing the design.

The 6% step has a 5.66% adjacent downward change, clearing the retained 5% minimum-change threshold. This is an explicit experimental rule change, not a change to rule_tacc() or to previously saved HCRs. A zero-height rule still implements the common initial-year TACC before the first advice takes effect.

Expanded outputs preserve both the fixed management reference B_R and the iteration-specific historical AW biomass reference. End-year abundance vectors are list-columns in the RDS, accompanied by boundaries and units in the manifest. Requested cut-offs must coincide with bin boundaries, including the open-ended 90 mm-plus group. No within-bin interpolation is performed. Catch-component fields distinguish commercial, recreational, and other removals. Seasonal U is the native OM fishery-1 exploitation parameter, not AW catch divided by end-year biomass. Inactive F or U parameterisations are NA rather than misleading zeros.

summarise_mse_indicators() reports explicit annual/evaluation windows with central 95% simulation ranges and finite denominators. Its window probabilities pool iteration-years; ever-crossing indicators are separate trajectory fields. These intervals are not Monte Carlo confidence intervals. Full numbers-at-size remain in the RDS rather than being flattened into the summary CSVs.