Skip to contents

Simulation and projections

Source: vignettes/projection.qmd

Scope

This vignette separates two related tasks:

  1. simulating observations from a fitted historical likelihood; and
  2. projecting population dynamics and simulated monitoring data through the CCSBT candidate target procedure (CTP).

Start projections from one accepted fit

Production projections require an accepted sbt_fit containing a compatible, validated posterior. The bundled package optimum is a deliberately restricted API fixture, and the archived MCMC and projection files shipped with earlier package versions are not scientific projection inputs.

The executable example below uses the fixture only to demonstrate historical simulation. The complete projection recipe is not evaluated during the website build because it requires an accepted posterior and is computationally expensive.

Historical simulation

Load the package fixture through the current staged fit lifecycle. The RTMB objective is transient and is rebuilt from the stored data and configuration.

library(sbt)
library(tidyverse)

theme_set(theme_bw())

article_extdata <- function(file) {
  candidates <- c(
    file.path("inst", "extdata", file),
    file.path("..", "inst", "extdata", file),
    system.file("extdata", file, package = "sbt")
  )
  candidates <- candidates[nzchar(candidates) & file.exists(candidates)]
  if (!length(candidates)) {
    stop("Could not find extdata file: ", file, call. = FALSE)
  }
  candidates[[1]]
}

fixture <- new.env(parent = emptyenv())
load(article_extdata("opt.rda"), envir = fixture)

fit <- sbt_fit(
  fixture$data,
  metadata = list(run = "bundled conditional API fixture")
)
fit <- sbt_add_parameters(fit, parameters = fixture$parameters)
fit <- sbt_add_map(fit, map = fixture$map)
fit <- sbt_add_bounds(fit, bounds = fixture$bounds)
fit <- sbt_build_object(fit)
obj <- sbt_obj(fit)
invisible(obj$fn(fixture$opt$par))
fit <- sbt_add_optimisation(
  fit,
  opt = fixture$opt,
  estimability = check_estimability(obj = obj),
  diagnostics = list(optimizer = list(method = "nlminb", n_passes = 2L)),
  optimizer = "nlminb",
  check = FALSE
)
obj <- sbt_obj(fit)
data <- fit$data

Data passed through RTMB’s OBS() mechanism can be simulated from their historical likelihood (Kristensen et al. 2016; Kristensen 2026). The CPUE example in Figure 1 draws from the fitted observation model. It demonstrates the API, not the adequacy of the conditional fixture.

set.seed(102)
plot(
  data$cpue_year,
  data$cpue_obs,
  type = "p",
  pch = 16,
  col = "firebrick",
  xlab = "Year",
  ylab = "CPUE"
)
for (iteration in seq_len(10L)) {
  lines(
    data$cpue_year,
    exp(obj$simulate()$cpue_log_obs),
    col = grDevices::adjustcolor("grey30", alpha.f = 0.35)
  )
}
points(data$cpue_year, data$cpue_obs, pch = 16, col = "firebrick")
Figure 1: Observed CPUE and ten simulations from the conditional API fixture.

The historical likelihood currently provides simulation for:

  • cpue_log_obs
  • troll_log_obs
  • aerial_log_obs
  • gt_nrec
  • hsp_nK
  • pop_nP

Age- and length-composition components do not all provide simulation methods, so a whole-model checkConsistency() call is not yet available. The component-specific OSA residual functions remain available:


Diagnostic Summary:
  SDNR: 1.83349  (95% CI 1.54128-2.26348; Target: ~1.0)
  MAR:  1.35812  (Target: ~0.67)
Figure 2: One-step-ahead CPUE residual API output for the conditional fixture.

Current projection workflow

The projection workflow keeps the accepted fit intact, adapts its normalized posterior only at the external SparseNUTS boundary, and uses a fresh objective for mutable projection work.

fit <- sbt_fit_read("accepted-base-model.sbt.rds", strict = TRUE)
posterior <- as_tmbfit(fit)
obj <- sbt_obj(fit, fresh = TRUE)

posterior_draws <- SparseNUTS::extract_samples(posterior)
set.seed(44)
iters <- sample.int(nrow(posterior_draws), size = 100L)

first_yr <- 2022L
last_yr <- 2035L
projection_years <- first_yr:last_yr

as_tmbfit() is needed because the low-level projection helpers call SparseNUTS::extract_samples(). The package-owned fit$mcmc remains the portable source of truth.

Recruitment and selectivity

Generate recruitment and selectivity inputs in exactly the same draw order as iters. The example below pins the first projected recruitment deviate to its fitted 2022 value, bootstraps future innovations, retains fitted selectivity through 2025, and then repeats the mean selectivity from 2016–2025.

set.seed(102)
projected_recruitment <- project_rec_devs(
  data = fit$data,
  obj = obj,
  mcmc = posterior,
  first_yr = first_yr,
  last_yr = last_yr,
  samp_years = fit$data$first_yr:(first_yr - 1L),
  iters = iters,
  pin_year = 2022L,
  option = "auto",
  bootstrap = TRUE
)

projected_selectivity <- project_selectivity(
  data = fit$data,
  obj = obj,
  mcmc = posterior,
  first_yr = first_yr,
  last_yr = last_yr,
  samp_years = 2016L:2025L,
  iters = iters,
  option = "mean_last",
  n_years = 10L,
  retain_fitted_through = 2025L
)

project_rec_devs() changes the caller’s random-number state, so set the seed explicitly. project_selectivity(option = "mean_last") is deterministic; option = "lognormal" also consumes random numbers.

CTP calendar

project_ctp_schedule() makes the TAC-update and data-availability assumptions visible before a costly run. Table 1 shows the current ESC31 example: historical dynamics through 2025, four fixed future TAC years, and a CTP update every three years thereafter.

ctp_schedule <- project_ctp_schedule(
  first_yr = 2022L,
  last_yr = 2035L,
  data_last_yr = 2025L,
  fixed_catch_n_years = 4L,
  tac_schedule = 3L,
  tac_calculation_lag = 2L,
  catch_data_lag = 3L,
  cpue_data_lag = 3L,
  gt_data_lag = 4L,
  ckmr_data_lag = 8L
)
knitr::kable(ctp_schedule)
Table 1: Example CTP TAC and data-availability schedule.
year tac_change tac_calculation_year catch_data_year cpue_data_year gt_data_year ckmr_data_year ctp_update
2026 hardwired 2022 NA NA NA NA FALSE
2027 hardwired 2025 2024 2024 2023 2019 FALSE
2028 hardwired NA 2025 2025 2024 2020 FALSE
2029 hardwired NA 2026 2026 2025 2021 FALSE
2030 Yes 2028 2027 2027 2026 2022 TRUE
2031 No NA 2028 2028 2027 2023 FALSE
2032 No NA 2029 2029 2028 2024 FALSE
2033 Yes 2031 2030 2030 2029 2025 TRUE
2034 No NA 2031 2031 2030 2026 FALSE
2035 No NA 2032 2032 2031 2027 FALSE

Run and cache the CTP projection

The high-level runner validates every selected historical biological state, simulates future monitoring observations, applies the CTP at scheduled update years, and checks projected dynamics. It records nominal TAC separately from biological removals, so any fleet-specific removal multipliers must be supplied explicitly.

projection <- run_projections(
  data = fit$data,
  object = obj,
  mcmc = posterior,
  first_yr = first_yr,
  last_yr = last_yr,
  n_iter = length(iters),
  iters = iters,
  seed = 102L,
  cores = 4L,
  fixed_catch_n_years = 4L,
  fixed_projection_tac = c(22670.91, 23647, 23647, 23647),
  ctp_tac_schedule = 3L,
  ctp_tac_calculation_lag = 2L,
  ctp_catch_data_lag = 3L,
  ctp_cpue_data_lag = 3L,
  ctp_gt_data_lag = 4L,
  ctp_ckmr_data_lag = 8L,
  projection_removal_multiplier_f = c(
    LL1 = 1.11,
    LL2 = 1,
    LL3 = 1,
    LL4 = 1,
    Indonesia = 1,
    Australia = 1.20
  ),
  rdev_y = projected_recruitment$proj_rdev_y[
    , as.character(projection_years), drop = FALSE
  ],
  sel_fya = projected_selectivity,
  cache_file = "projection-cache.rds",
  overwrite = FALSE,
  cache_only = FALSE,
  verbose = TRUE
)

If exact fleet allocation proportions are part of the scientific contract, supply them with projection_tac_split_yf or projection_tac_split_f; otherwise the runner uses its documented recent-allocation rule.

Compatible caches are checked against the projection inputs and internal scientific contract before reuse. New caches are written atomically. A cache hit is convenient, but acceptance still requires the saved diagnostics and reported outputs to be reviewed.

Paired scenarios

For scenario comparisons, reuse the same posterior rows, recruitment, selectivity, projection seed, TAC assumptions, and allocation. Change only the scenario input. This common-random-number design makes differences attributable to the scenario rather than Monte Carlo noise.

Gene-tag release years and recapture or harvest-sampling years are separate controls:

gt_skip_projection <- run_projections(
  data = fit$data,
  object = obj,
  mcmc = posterior,
  first_yr = first_yr,
  last_yr = last_yr,
  n_iter = length(iters),
  iters = iters,
  seed = 102L,
  cores = 4L,
  fixed_catch_n_years = 4L,
  fixed_projection_tac = c(22670.91, 23647, 23647, 23647),
  ctp_tac_schedule = 3L,
  ctp_tac_calculation_lag = 2L,
  ctp_catch_data_lag = 3L,
  ctp_cpue_data_lag = 3L,
  ctp_gt_data_lag = 4L,
  ctp_ckmr_data_lag = 8L,
  gt_skip_years = c(2026L, 2028L, 2030L, 2032L, 2034L),
  gt_skip_recapture_years = c(2027L, 2029L, 2031L, 2033L, 2035L),
  projection_removal_multiplier_f = c(
    LL1 = 1.11,
    LL2 = 1,
    LL3 = 1,
    LL4 = 1,
    Indonesia = 1,
    Australia = 1.20
  ),
  rdev_y = projected_recruitment$proj_rdev_y[
    , as.character(projection_years), drop = FALSE
  ],
  sel_fya = projected_selectivity,
  cache_file = "projection-gt-skip-cache.rds"
)

The exact skip vectors above are illustrative. Record the chosen release and recapture schedule alongside every result so that “alternate-year gene tagging” is unambiguous.

Reviewing projection results

At minimum, review:

  • the MLE and posterior acceptance records in the source fit;
  • the selected-draw plan and historical biological-state diagnostics;
  • CTP status and input snapshots for every draw and update year;
  • finite, feasible population trajectories;
  • nominal TACs separately from biological removals;
  • paired TAC and stock-status differences by year; and
  • Monte Carlo uncertainty, especially for short screening runs.

The returned ctp_total_tac_iy matrix contains draw-by-year nominal TACs, while ctp_total_removal_iy contains the corresponding removals after multipliers. For a paired comparison, summarize the within-draw differences rather than comparing two unrelated marginal summaries:

tac_difference <- gt_skip_projection$ctp_total_tac_iy -
  projection$ctp_total_tac_iy

data.frame(
  Year = projection$projection_years,
  Median_difference = apply(tac_difference, 2L, median),
  Lower_95 = apply(tac_difference, 2L, quantile, probs = 0.025),
  Upper_95 = apply(tac_difference, 2L, quantile, probs = 0.975)
)

A 100-draw run is useful for screening and checking implementation behaviour. It is not automatically adequate for final tail probabilities or management advice; those require a draw-count sensitivity check and the full acceptance workflow.

References

Kristensen, Kasper. 2026. RTMB: ’R’ Bindings for ’TMB’. https://doi.org/10.32614/CRAN.package.RTMB.
Kristensen, Kasper, Anders Nielsen, Casper W. Berg, Hans Skaug, and Bradley M. Bell. 2016. “TMB: Automatic Differentiation and Laplace Approximation.” Journal of Statistical Software 70 (5): 1–21. https://doi.org/10.18637/jss.v070.i05.