MLE and posterior acceptance are reported separately. Every sensitivity MLE must pass the numerical and biological gates before it is shown. A posterior is included only when its embedded MCMC matches the current scientific and sampler identities and passes the stated acceptance decision. In the accepted production artifacts summarized here, all posteriors meet the standard diagnostics except NoPOPHSP, for which the explicitly reviewed 12,000-draw run is accepted with two divergent transitions. Its otherwise passing diagnostics and complete biological-state validation remain visible below. A final isolated recovery run removed the divergences but failed the unchanged R-hat and ESS gates, so it was rejected and did not replace the canonical fit. A deliberately partial render may omit one or more artifacts; the tables and figures below then report those posteriors as pending rather than accepted. The conventional-tag-exclusion sensitivity has an accepted 7,200-draw posterior after its longer production run passed every standard gate.
These sensitivities inherit the current base treatment: LL3 and LL4 remain separate standard selectivity-based fisheries with soft length-composition likelihoods. LL4 has one time-invariant selectivity block beginning in 1953, with 14 estimated age effects for ages 8–21 and terminal extension of the age-21 effect. Its fixed selectivity hyperparameters inherit LL3’s values: rho_year = 0.5, rho_age = 0.5, and sigma = 0.75. LL3 and LL4 are not lumped, and LL4 retains its complete observed catch through the standard exact catch-conditioning calculation. The length-based base MLE and posterior have been reviewed; this page keeps each data change, MLE, MCMC definition, and available diagnostic in one place.
Accepted MLEs and posterior draws require combined raw seasonal harvest at age no greater than 0.9, positive abundance, exact catch accounting, and zero continuation penalty within the numerical acceptance tolerance. All fits also use the approved global preventive harvest wall: a normalized squared-softplus contribution with strength 10, onset 0.85, ceiling 0.9, and scale 0.01. The wall does not change catches or population dynamics, and the strict state gates and feasibility continuation remain unchanged.
options(readr.show_col_types = FALSE)
invisible(knitr::knit_meta_add(list(rmarkdown::html_dependency_jquery())))
library(knitr)
library(tidyverse)
library(RTMB)
library(TMBhelper)
library(scales)
library(sbt)
theme_set(theme_bw())
`%||%` <- function(x, y) {
if (is.null(x) || length(x) == 0L) y else x
}if (basename(getwd()) == "ESC31") {
esc_dir <- normalizePath(".")
} else if (dir.exists("ESC31")) {
esc_dir <- normalizePath("ESC31")
} else {
stop("Could not find the ESC31 directory.", call. = FALSE)
}
run_dir <- file.path(esc_dir, "runs")
dir.create(run_dir, recursive = TRUE, showWarnings = FALSE)
supporting_run_dir <- file.path(esc_dir, "supporting")
dir.create(supporting_run_dir, recursive = TRUE, showWarnings = FALSE)
sensitivity_run_dir <- file.path(run_dir, "sens")
dir.create(sensitivity_run_dir, recursive = TRUE, showWarnings = FALSE)
previous_v1_model_file <- file.path(
supporting_run_dir,
"esc31_previous_assessment_v1_h070.sbt.rds"
)
if (!file.exists(previous_v1_model_file)) {
stop(
"The refitted V1 comparison is unavailable; render 2_base.qmd first.",
call. = FALSE
)
}
previous_v1_fit <- sbt_fit_read(
previous_v1_model_file,
strict = TRUE,
rebuild = FALSE,
verify_validation = FALSE
)
previous_v1_report <- sbt_fit_report(previous_v1_fit)
model_file <- file.path(run_dir, "esc31_base.rds")
sensitivity_files <- c(
no_uam = file.path(sensitivity_run_dir, "esc31_sens_no_uam.sbt.rds"),
drop_5yrs = file.path(sensitivity_run_dir, "esc31_sens_drop_5yrs.sbt.rds"),
cpue_omega_075 = file.path(sensitivity_run_dir, "esc31_sens_cpue_omega_075.sbt.rds"),
q2008 = file.path(sensitivity_run_dir, "esc31_sens_q2008.sbt.rds"),
ll1_terminal_3yr = file.path(sensitivity_run_dir, "esc31_sens_ll1_terminal_3yr.sbt.rds"),
no_pop_hsp = file.path(sensitivity_run_dir, "esc31_sens_no_pop_hsp.sbt.rds"),
no_hsp = file.path(sensitivity_run_dir, "esc31_sens_no_hsp.sbt.rds"),
troll = file.path(sensitivity_run_dir, "esc31_sens_troll.sbt.rds"),
cpue_ll1_sel = file.path(sensitivity_run_dir, "esc31_sens_cpue_ll1_sel.sbt.rds"),
constant_cpue_cv = file.path(sensitivity_run_dir, "esc31_sens_constant_cpue_cv.sbt.rds"),
# This established final registry key now contains the direct four-anchor
# comparison because length-based M is the selected base.
length_m = file.path(sensitivity_run_dir, "esc31_sens_length_m.sbt.rds"),
# Append new targets so every established registry position and seed remains
# unchanged.
indo_sel = file.path(sensitivity_run_dir, "esc31_sens_indo_sel.sbt.rds"),
estimate_m_slope = file.path(
sensitivity_run_dir,
"esc31_sens_estimate_m_slope.sbt.rds"
),
# Append the conventional-tag exclusion so all existing registry positions
# and sampler seeds remain unchanged.
no_tags = file.path(
sensitivity_run_dir,
"esc31_sens_no_tags.sbt.rds"
)
)
sensitivity_mle_files <- sensitivity_files
sensitivity_mcmc_registry_keys <- c("base", names(sensitivity_files))
rejected_sensitivity_mle_keys <- character()
valid_sensitivity_mcmc_keys <- setdiff(
sensitivity_mcmc_registry_keys,
rejected_sensitivity_mle_keys
)
sensitivity_mcmc_registry_files <- c(
base = model_file,
sensitivity_files
)
sensitivity_mcmc_files <- sensitivity_mcmc_registry_files[
valid_sensitivity_mcmc_keys
]
source(file.path(esc_dir, "esc31_workflow.R"))
source(file.path(esc_dir, "esc31_inputs.R"), local = TRUE)
source(file.path(esc_dir, "esc31_base_posterior.R"), local = TRUE)
esc31_require_ll34_conditioning_approval(
base_model_contract,
stage = "sensitivity MLE and MCMC"
)
sensitivity_harvest_wall_contract <- harvest_wall_contract(
strength = 10,
onset = 0.85,
ceiling = 0.9,
scale = 0.01
)
sensitivity_harvest_wall_recomputation_tolerance <- 1e-8
sensitivity_harvest_wall_decision <-
base_model_contract$scientific_config$harvest_wall
stopifnot(
esc31_harvest_wall_is_approved(base_model_contract),
identical(
sensitivity_harvest_wall_decision$decision_id,
"esc31_2026_preventive_harvest_wall_a10_v1"
),
identical(
sensitivity_harvest_wall_contract$implementation,
sensitivity_harvest_wall_decision$implementation
),
identical(
unname(unlist(
sensitivity_harvest_wall_contract[
c("strength", "onset", "ceiling", "scale")
],
use.names = FALSE
)),
c(10, 0.85, 0.9, 0.01)
)
)
if (!exists("refit_model", inherits = FALSE)) refit_model <- FALSE
# Direct-render controls. The data changes, MLEs, MCMCs, and diagnostics all
# remain defined in this QMD. The command-line workers only set these same
# controls when one model is run independently.
rerun_sensitivities <- FALSE
run_sensitivity_fits <- TRUE
run_sensitivity_mcmc <- FALSE
sensitivity_mcmc_keys <- "all"
worker_fit_flag <- trimws(Sys.getenv("ESC31_RUN_SENSITIVITY_FITS", ""))
if (nzchar(worker_fit_flag)) {
run_sensitivity_fits <- tolower(worker_fit_flag) %in%
c("1", "true", "yes", "on")
}
worker_mcmc_flag <- trimws(Sys.getenv("ESC31_RUN_SENSITIVITY_MCMC", ""))
if (nzchar(worker_mcmc_flag)) {
run_sensitivity_mcmc <- tolower(worker_mcmc_flag) %in%
c("1", "true", "yes", "on")
}
worker_mcmc_keys <- trimws(Sys.getenv("ESC31_SENSITIVITY_MCMC_KEYS", ""))
if (nzchar(worker_mcmc_keys)) {
sensitivity_mcmc_keys <- trimws(strsplit(
worker_mcmc_keys,
",",
fixed = TRUE
)[[1]])
sensitivity_mcmc_keys <- sensitivity_mcmc_keys[
nzchar(sensitivity_mcmc_keys)
]
}
if (run_sensitivity_mcmc && !length(sensitivity_mcmc_keys)) {
stop(
"Set ESC31_SENSITIVITY_MCMC_KEYS to one or more model keys or `all`.",
call. = FALSE
)
}
rejected_requested_mcmc_keys <- intersect(
sensitivity_mcmc_keys,
rejected_sensitivity_mle_keys
)
if (length(rejected_requested_mcmc_keys)) {
stop(
"MCMC is prohibited for non-executable sensitivity key(s): ",
paste(rejected_requested_mcmc_keys, collapse = ", "),
".",
call. = FALSE
)
}
unknown_sensitivity_mcmc_keys <- setdiff(
sensitivity_mcmc_keys,
c("all", valid_sensitivity_mcmc_keys)
)
if (length(unknown_sensitivity_mcmc_keys)) {
stop(
"Unknown sensitivity MCMC key(s): ",
paste(unknown_sensitivity_mcmc_keys, collapse = ", "),
call. = FALSE
)
}
if (run_sensitivity_mcmc) {
esc31_require_ll34_conditioning_approval(
base_model_contract,
stage = "sensitivity MCMC"
)
run_sensitivity_fits <- FALSE
}expected_base_identity <- esc31_base_run_identity(
esc_dir,
base_model_contract = base_model_contract
)
base_fit_candidate <- if (file.exists(model_file)) {
tryCatch(
sbt_fit_read(
model_file,
strict = TRUE,
rebuild = FALSE,
verify_validation = FALSE
),
error = function(e) NULL
)
} else {
NULL
}
if (refit_model || !esc31_base_fit_is_current(base_fit_candidate, expected_base_identity)) {
if (run_sensitivity_mcmc) {
stop(
"The base MLE cache must be finalized before launching sensitivity MCMC workers.",
call. = FALSE
)
}
mod_qmd <- file.path(esc_dir, "2_base.qmd")
quarto::quarto_render(mod_qmd, quiet = TRUE)
}
base_fit <- sbt_fit_read(
model_file,
strict = TRUE,
rebuild = FALSE,
verify_validation = FALSE
)
data_in <- base_model_contract$data_informat_decimal <- function(x, digits = 3) {
out <- rep(NA_character_, length(x))
finite <- is.finite(x)
x[finite & abs(x) < 0.5 * 10^(-digits)] <- 0
out[finite] <- formatC(x[finite], format = "f", digits = digits, big.mark = ",")
out[!finite & !is.na(x)] <- as.character(x[!finite & !is.na(x)])
out
}
format_scientific <- function(x, digits = 2) {
out <- rep(NA_character_, length(x))
finite <- is.finite(x)
out[finite] <- formatC(x[finite], format = "e", digits = digits)
out[!finite & !is.na(x)] <- as.character(x[!finite & !is.na(x)])
out
}
replace_missing <- function(x, missing = "-") {
x <- as.data.frame(x)
x[] <- lapply(x, function(column) {
column <- as.character(column)
column[is.na(column) | column == "NA"] <- missing
column
})
as_tibble(x)
}
plot_cpue_osa_quiet <- function(fit) {
capture.output(p <- plot_cpue_residuals(fit))
p
}
summarize_fit <- function(label, fit) {
par_list <- fit$parameters
report <- sbt_fit_report(fit)
opt_i <- fit$fit$opt
nll <- opt_i$objective
n_parameters <- length(opt_i$par)
tibble(
Model = label,
`Convergence code` = opt_i$convergence,
`Penalized objective` = nll,
`AIC-style score` = 2 * nll + 2 * n_parameters,
`Max gradient` = fit$fit$diagnostics$max_gradient %||% NA_real_,
B0 = exp(par_list$par_log_B0),
M0 = as.numeric(report$par_m0),
M4 = as.numeric(report$par_m4),
M10 = as.numeric(report$par_m10),
M30 = as.numeric(report$par_m30),
h = exp(par_list$par_log_h),
psi = exp(par_list$par_log_psi)
)
}
has_sensitivity_fit <- function(x) {
inherits(x, "sbt_fit")
}
sensitivity_mle_thresholds <- list(
convergence = 0L,
max_gradient = 0.01,
estimability = "estimable",
biological_state = biological_state_contract()
)
sensitivity_optimizer_contract <- list(
schema_version = 2L,
implementation = "esc31_staged_bounded_nlminb_newton_v2",
primary = list(
method = "bounded_nlminb_exact_hessian",
max_passes = 1L,
control = base_fit$control,
transition = list(
policy = "single_exact_pass_then_scaled_v1",
certified = "accept",
finite_nonworsening_uncertified = "retain_best_then_scaled",
invalid_or_worsening = "rollback_then_scaled",
branch_on_message = FALSE
)
),
scaled = list(
method = "hessian_diagonal_scaled_bounded_gradient_nlminb",
max_passes = 6L,
control = list(
eval.max = 15000L,
iter.max = 5000L,
rel.tol = 1e-10,
x.tol = 1e-8
),
hessian_diagonal_floor = 1e-8,
scale_min = 1e-4,
scale_max = 10
),
newton = list(
method = "bounded_spd_newton_armijo",
enabled = TRUE,
minimum_reciprocal_condition = 1e-12,
initial_alpha = 1,
fraction_to_boundary = 0.995,
backtrack_factor = 0.5,
maximum_backtracks = 12L,
armijo_c1 = 1e-4,
require_gradient_improvement = TRUE
),
certification = list(
method = "bounded_gradient_nlminb",
max_passes = 1L,
control = list(eval.max = 2000L, iter.max = 1000L)
),
acceptance = c(
sensitivity_mle_thresholds,
list(
require_finite_parameters = TRUE,
require_finite_objective = TRUE,
require_finite_gradient = TRUE,
require_in_bounds = TRUE,
bound_tolerance = 1e-10
)
)
)
sensitivity_optimizer_contract_signature <- esc31_object_md5(
sensitivity_optimizer_contract
)
sensitivity_optimizer_provenance_passes <- function(fit) {
if (!has_sensitivity_fit(fit)) return(FALSE)
fit_diagnostics <- fit$fit$diagnostics
max_gradient <- fit_diagnostics$max_gradient %||% NA_real_
optimizer_history <- fit_diagnostics$optimizer_history
optimizer_summary <- fit_diagnostics$optimizer_summary
required_history_fields <- c(
"pass", "stage", "attempt", "objective", "max_abs_gradient",
"gradient_l2", "convergence", "within_bounds", "finite_gradient",
"accepted_as_best", "previous_objective",
"relative_objective_decrease", "gradient_ratio",
"relative_parameter_step", "transition", "transition_reason", "message"
)
primary_history <- if (
is.data.frame(optimizer_history) && "stage" %in% names(optimizer_history)
) {
optimizer_history[
optimizer_history$stage == "primary bounded nlminb",
,
drop = FALSE
]
} else {
data.frame()
}
allowed_primary_transitions <- unname(unlist(
sensitivity_optimizer_contract$primary$transition[
c(
"certified", "finite_nonworsening_uncertified",
"invalid_or_worsening"
)
],
use.names = FALSE
))
identical(
fit$optimization$method,
sensitivity_optimizer_contract$implementation
) &&
identical(
fit$provenance$metadata$optimizer_contract_signature,
sensitivity_optimizer_contract_signature
) &&
identical(
fit_diagnostics$optimizer_contract_signature,
sensitivity_optimizer_contract_signature
) &&
identical(
fit_diagnostics$optimizer_contract,
sensitivity_optimizer_contract
) &&
is.data.frame(optimizer_history) &&
nrow(optimizer_history) >= 1L &&
all(required_history_fields %in% names(optimizer_history)) &&
!anyNA(optimizer_history$transition) &&
all(nzchar(optimizer_history$transition)) &&
!any(optimizer_history$transition == "pending") &&
nrow(primary_history) == 1L &&
is.list(optimizer_summary) &&
identical(optimizer_summary$schema_version, 2L) &&
isTRUE(optimizer_summary$certified) &&
identical(
optimizer_summary$contract_signature,
sensitivity_optimizer_contract_signature
) &&
identical(
as.integer(optimizer_summary$total_nlminb_calls),
as.integer(fit$fit$diagnostics$optimizer$n_passes)
) &&
identical(
as.integer(optimizer_summary$final_convergence),
as.integer(fit$fit$opt$convergence)
) &&
isTRUE(all.equal(
as.numeric(optimizer_summary$final_objective),
as.numeric(fit$fit$opt$objective),
tolerance = 1e-10
)) &&
isTRUE(all.equal(
as.numeric(optimizer_summary$final_max_abs_gradient),
as.numeric(max_gradient),
tolerance = 1e-10
)) &&
identical(as.integer(optimizer_summary$primary_passes), 1L) &&
is.character(optimizer_summary$primary_transition) &&
length(optimizer_summary$primary_transition) == 1L &&
optimizer_summary$primary_transition %in% allowed_primary_transitions &&
is.character(optimizer_summary$primary_transition_reason) &&
length(optimizer_summary$primary_transition_reason) == 1L &&
nzchar(optimizer_summary$primary_transition_reason) &&
identical(
optimizer_summary$primary_transition,
primary_history$transition[[1L]]
) &&
identical(
optimizer_summary$primary_transition_reason,
primary_history$transition_reason[[1L]]
)
}
sensitivity_harvest_wall_report_diagnostics <- function(data_fit, report) {
fail <- function(reason) {
list(
schema_version = 1L,
decision_id = sensitivity_harvest_wall_decision$decision_id,
contract = sensitivity_harvest_wall_contract,
recomputation_tolerance =
sensitivity_harvest_wall_recomputation_tolerance,
passes = FALSE,
reason = reason
)
}
if (!is.list(data_fit) || !is.list(report)) {
return(fail("The fit data or model report is unavailable."))
}
required_data <- c(
"harvest_wall_strength", "harvest_wall_onset",
"harvest_wall_ceiling", "harvest_wall_scale"
)
required_report <- c(
"hrate_raw_ysa", "harvest_wall_penalty_ysa", "lp_harvest_wall",
"harvest_wall_strength", "harvest_wall_onset",
"harvest_wall_ceiling", "harvest_wall_scale"
)
missing_data <- setdiff(required_data, names(data_fit))
missing_report <- setdiff(required_report, names(report))
if (length(missing_data) || length(missing_report)) {
return(fail(paste0(
"Missing wall data/report fields: ",
paste(c(missing_data, missing_report), collapse = ", "), "."
)))
}
raw_harvest <- report$hrate_raw_ysa
reported_penalty <- report$harvest_wall_penalty_ysa
reported_total <- report$lp_harvest_wall
reported_scalars <- lapply(
report[c(
"harvest_wall_strength", "harvest_wall_onset",
"harvest_wall_ceiling", "harvest_wall_scale"
)],
as.numeric
)
if (!is.numeric(raw_harvest) || !length(raw_harvest) ||
any(!is.finite(raw_harvest)) ||
!is.numeric(reported_penalty) ||
length(reported_penalty) != length(raw_harvest) ||
any(!is.finite(reported_penalty)) ||
!identical(dim(reported_penalty), dim(raw_harvest)) ||
!is.numeric(reported_total) || length(reported_total) != 1L ||
!is.finite(reported_total) ||
any(lengths(reported_scalars) != 1L) ||
any(!is.finite(unlist(reported_scalars, use.names = FALSE)))) {
return(fail("The reported wall arrays or scalars are non-finite or malformed."))
}
reported_contract <- harvest_wall_contract(
strength = reported_scalars$harvest_wall_strength,
onset = reported_scalars$harvest_wall_onset,
ceiling = reported_scalars$harvest_wall_ceiling,
scale = reported_scalars$harvest_wall_scale
)
data_contract <- harvest_wall_contract(
strength = data_fit$harvest_wall_strength,
onset = data_fit$harvest_wall_onset,
ceiling = data_fit$harvest_wall_ceiling,
scale = data_fit$harvest_wall_scale
)
recomputed_penalty <- get_harvest_wall_penalty(
raw_harvest,
contract = sensitivity_harvest_wall_contract
)
maximum_penalty_error <- max(abs(
as.numeric(reported_penalty) - as.numeric(recomputed_penalty)
))
total_penalty_error <- abs(
as.numeric(reported_total) - sum(as.numeric(recomputed_penalty))
)
maximum_raw_harvest <- max(as.numeric(raw_harvest))
raw_harvest_within_ceiling <- maximum_raw_harvest <=
sensitivity_mle_thresholds$biological_state$hrate_limit +
sensitivity_mle_thresholds$biological_state$hrate_tolerance
data_contract_matches <- identical(
data_contract,
sensitivity_harvest_wall_contract
)
reported_contract_matches <- identical(
reported_contract,
sensitivity_harvest_wall_contract
)
nonnegative_penalty <- all(
as.numeric(reported_penalty) >=
-sensitivity_harvest_wall_recomputation_tolerance
)
passes <- data_contract_matches && reported_contract_matches &&
raw_harvest_within_ceiling && nonnegative_penalty &&
maximum_penalty_error <=
sensitivity_harvest_wall_recomputation_tolerance &&
total_penalty_error <= sensitivity_harvest_wall_recomputation_tolerance
list(
schema_version = 1L,
decision_id = sensitivity_harvest_wall_decision$decision_id,
contract = sensitivity_harvest_wall_contract,
recomputation_tolerance = sensitivity_harvest_wall_recomputation_tolerance,
raw_harvest_cells = length(raw_harvest),
cells_above_onset = sum(
as.numeric(raw_harvest) > sensitivity_harvest_wall_contract$onset
),
maximum_raw_harvest = maximum_raw_harvest,
total_wall_objective = as.numeric(reported_total),
maximum_cell_penalty = max(as.numeric(reported_penalty)),
maximum_penalty_recomputation_error = maximum_penalty_error,
total_penalty_recomputation_error = total_penalty_error,
data_contract_matches = data_contract_matches,
reported_contract_matches = reported_contract_matches,
raw_harvest_within_ceiling = raw_harvest_within_ceiling,
nonnegative_penalty = nonnegative_penalty,
payload_checksum = esc31_object_md5(list(
raw_harvest = raw_harvest,
reported_penalty = reported_penalty
)),
passes = passes,
reason = if (passes) "pass" else "One or more harvest-wall gates failed."
)
}
sensitivity_harvest_wall_fit_diagnostics <- function(fit) {
if (!has_sensitivity_fit(fit)) {
return(sensitivity_harvest_wall_report_diagnostics(NULL, NULL))
}
tryCatch(
sensitivity_harvest_wall_report_diagnostics(
fit$data,
sbt_fit_report(fit)
),
error = function(e) {
failed <- sensitivity_harvest_wall_report_diagnostics(NULL, NULL)
failed$reason <- paste(
"Harvest-wall diagnostic reconstruction failed:",
conditionMessage(e)
)
failed
}
)
}
sensitivity_mle_passes <- function(fit) {
if (!has_sensitivity_fit(fit)) return(FALSE)
validation <- tryCatch(
sbt_fit_validation(
fit,
scope = "mle",
verify = !identical(sensitivity_fit_validation_mode, "skip"),
require_pass = TRUE
),
error = function(e) NULL
)
max_gradient <- fit$fit$diagnostics$max_gradient %||% NA_real_
biological_state <- fit$fit$diagnostics$biological_state_mle
wall <- fit$fit$diagnostics$harvest_wall_mle
stored_wall_passes <- is.list(wall) &&
is.list(wall$identity) &&
identical(
wall$identity$decision,
sensitivity_harvest_wall_decision
) &&
identical(
wall$identity$executable,
sensitivity_harvest_wall_contract
) &&
is.finite(wall$maximum_raw_harvest) &&
wall$maximum_raw_harvest <=
sensitivity_mle_thresholds$biological_state$hrate_limit +
sensitivity_mle_thresholds$biological_state$hrate_tolerance &&
is.finite(wall$maximum_recomputation_error) &&
wall$maximum_recomputation_error <=
sensitivity_harvest_wall_recomputation_tolerance &&
is.finite(wall$maximum_continuation_penalty) &&
abs(wall$maximum_continuation_penalty) <=
sensitivity_mle_thresholds$biological_state$penalty_tolerance
!is.null(validation) &&
identical(as.integer(fit$fit$opt$convergence),
sensitivity_mle_thresholds$convergence) &&
length(max_gradient) == 1L &&
is.finite(max_gradient) &&
max_gradient <= sensitivity_mle_thresholds$max_gradient &&
identical(
fit$fit$estimability$status,
sensitivity_mle_thresholds$estimability
) &&
esc31_biological_state_diagnostics_pass(biological_state, 1L) &&
stored_wall_passes
}
sensitivity_mcmc_default_config <- list(
num_samples = 750L,
num_warmup = 150L,
chains = 4L,
cores = 4L,
metric = "dense",
init = "last.par.best",
adapt_delta = 0.999,
max_treedepth = 13L,
refresh = 200L,
skip_optimization = TRUE
)
# Key-specific implementation controls must remain provenance-bound. The
# NoPOPHSP fixed-effect Hessian is positive definite under RTMB's exact AD
# derivatives, but RTMB::sdreport()'s finite-difference Hessian is not. Passing
# the exact-Hessian covariance preserves the approved dense metric and avoids a
# false preconditioner failure without changing the sampler settings.
sensitivity_mcmc_overrides <- list(
no_pop_hsp = list(
dense_metric_precision_source =
"RTMB exact AD Hessian at source MLE",
# The 1,200-draw rerun still had maximum R-hat 1.025406 and minimum bulk
# ESS 243.3. Preserve its seed and transition controls, but extend the
# retained chain substantially so the effective sample size has a
# reasonable prospect of clearing 400 with margin.
num_samples = 3000L
),
# The first annual-terminal-LL1 run narrowly missed only the strict R-hat
# gate (1.010519 versus < 1.01). Retain its seed and transition controls and
# add 20% more retained draws.
ll1_terminal_3yr = list(
num_samples = 900L
),
# The first No-HSP run also had otherwise healthy diagnostics and missed
# only the strict R-hat gate (1.010443 versus < 1.01). Retain its seed and
# transition controls and add 20% more retained draws.
no_hsp = list(
num_samples = 900L
),
# The first 750-draw Constant-CPUE-CV run had otherwise healthy diagnostics
# but missed the strict R-hat gate narrowly (1.01069 versus < 1.01). Retain
# the seed and transition controls and add 20% more retained draws.
constant_cpue_cv = list(
num_samples = 900L
),
# The first Troll run had maximum R-hat 1.014088 with otherwise healthy
# ESS and no sampler pathologies. Preserve its seed and use longer chains.
troll = list(
num_samples = 1200L
),
# The first estimated-slope run missed only the R-hat gate narrowly at
# 1.010244. Preserve its seed and add 20% more retained draws.
estimate_m_slope = list(
num_samples = 900L
),
# The first No-tags run had no sampler pathologies, but maximum R-hat was
# 1.017 and minimum bulk ESS was 219.7 for M10. Preserve its seed and
# transition controls, and extend the retained chains to provide adequate
# effective-sample-size margin.
no_tags = list(
num_samples = 1800L
)
)
sensitivity_mcmc_thresholds <- list(
max_rhat = 1.01,
min_bulk_ess = 400,
min_tail_ess = 400,
divergences = 0L,
max_treedepth_hits = 0L,
biological_state = biological_state_contract()
)
no_pop_hsp_mcmc_acceptance <- list(
schema_version = 1L,
status = "accepted_with_explicit_divergence_exception",
decision_id =
"esc31_2026_no_pop_hsp_accept_two_divergences_12000_draws_v1",
decision_date = "2026-07-23",
model_key = "no_pop_hsp",
retained_draws_per_chain = 3000L,
chains = 4L,
retained_draws_total = 12000L,
accepted_divergences = 2L,
standard_thresholds_otherwise_unchanged = TRUE,
biological_state_exception = FALSE,
rationale = paste(
"Explicit assessment decision to accept the two divergent transitions",
"after the longer NoPOPHSP run passed R-hat, bulk ESS, tail ESS,",
"maximum-treedepth, MLE, and complete retained-draw biological-state",
"checks."
)
)
# State rescans are deterministic across worker counts, so this performance
# control is deliberately excluded from every scientific identity.
sensitivity_state_diagnostic_cores <- local({
configured <- trimws(Sys.getenv("ESC31_STATE_DIAGNOSTIC_CORES", ""))
if (nzchar(configured)) {
if (!grepl("^[1-9][0-9]*$", configured)) {
stop(
"`ESC31_STATE_DIAGNOSTIC_CORES` must be one positive integer.",
call. = FALSE
)
}
cores <- suppressWarnings(as.integer(configured))
if (is.na(cores)) {
stop(
"`ESC31_STATE_DIAGNOSTIC_CORES` is outside the supported integer range.",
call. = FALSE
)
}
cores
} else {
detected <- suppressWarnings(parallel::detectCores(logical = TRUE))
if (length(detected) != 1L || is.na(detected) || detected < 1L) {
detected <- sensitivity_mcmc_default_config$cores
}
min(16L, as.integer(detected))
}
})
# Expensive validation is performed once when a fit is completed and stored in
# the portable `sbt_fit`. Normal report and cache-check paths reuse that
# payload-bound record. Set `ESC31_FIT_VALIDATION_MODE=full` for a deliberate
# end-to-end audit or `skip` to trust the stored record without recalculating
# its compact checksums.
sensitivity_fit_validation_mode <- local({
configured <- tolower(trimws(Sys.getenv(
"ESC31_FIT_VALIDATION_MODE", "auto"
)))
if (!configured %in% c("auto", "full", "skip")) {
stop(
"`ESC31_FIT_VALIDATION_MODE` must be `auto`, `full`, or `skip`.",
call. = FALSE
)
}
configured
})
sensitivity_mcmc_requested <- function(key) {
key %in% valid_sensitivity_mcmc_keys &&
run_sensitivity_mcmc &&
("all" %in% sensitivity_mcmc_keys || key %in% sensitivity_mcmc_keys)
}
sensitivity_mcmc_config <- function(key) {
if (length(key) != 1L || is.na(key) ||
!key %in% valid_sensitivity_mcmc_keys) {
stop("The requested sensitivity MCMC key is not executable.", call. = FALSE)
}
config <- utils::modifyList(
sensitivity_mcmc_default_config,
sensitivity_mcmc_overrides[[key]] %||% list()
)
config$seed <- 73000L + match(key, sensitivity_mcmc_registry_keys)
config
}
sensitivity_mcmc_identity <- function(key, source_fit, fit_identity_signature,
sampler_config) {
if (length(key) != 1L || is.na(key) ||
!key %in% valid_sensitivity_mcmc_keys) {
stop("Cannot create an MCMC identity for a non-executable sensitivity.",
call. = FALSE)
}
esc31_workflow_identity(
esc_dir,
workflow_files = character(),
sbt_entry_points = esc31_sbt_entry_points("sensitivities"),
config = list(
stage = "sensitivity_mcmc",
mcmc_workflow_version =
"esc31_sensitivity_mcmc_v13_length_base_check_mcmc",
model_key = key,
base_run_signature =
base_fit$provenance$metadata$run_identity$signature,
base_fit_scientific_signature = esc31_fit_scientific_signature(base_fit),
fit_identity_signature = fit_identity_signature,
fit_scientific_signature = esc31_fit_scientific_signature(source_fit),
harvest_wall_contract = sensitivity_harvest_wall_contract,
harvest_wall_decision = sensitivity_harvest_wall_decision,
harvest_wall_recomputation_tolerance =
sensitivity_harvest_wall_recomputation_tolerance,
mle_acceptance_thresholds = sensitivity_mle_thresholds,
sampler = sampler_config,
thresholds = sensitivity_mcmc_thresholds
)
)
}
sensitivity_mcmc_standard_passes <- function(diagnostics) {
is.data.frame(diagnostics) && nrow(diagnostics) == 1L &&
isTRUE(diagnostics$passes)
}
sensitivity_mcmc_divergence_exception_passes <- function(diagnostics) {
state <- sensitivity_mcmc_thresholds$biological_state
is.data.frame(diagnostics) && nrow(diagnostics) == 1L &&
identical(
as.character(diagnostics$Model),
no_pop_hsp_mcmc_acceptance$model_key
) &&
identical(
as.integer(diagnostics$chains),
no_pop_hsp_mcmc_acceptance$chains
) &&
identical(
as.integer(diagnostics$samples_per_chain),
no_pop_hsp_mcmc_acceptance$retained_draws_per_chain
) &&
identical(
as.integer(diagnostics$state_draws_expected),
no_pop_hsp_mcmc_acceptance$retained_draws_total
) &&
identical(
as.integer(diagnostics$state_draws_evaluated),
no_pop_hsp_mcmc_acceptance$retained_draws_total
) &&
identical(
as.integer(diagnostics$divergences),
no_pop_hsp_mcmc_acceptance$accepted_divergences
) &&
identical(as.integer(diagnostics$convergence), 0L) &&
is.finite(diagnostics$objective) &&
isTRUE(diagnostics$estimable) &&
is.finite(diagnostics$max_rhat) &&
diagnostics$max_rhat < sensitivity_mcmc_thresholds$max_rhat &&
is.finite(diagnostics$min_bulk_ess) &&
diagnostics$min_bulk_ess >= sensitivity_mcmc_thresholds$min_bulk_ess &&
is.finite(diagnostics$min_tail_ess) &&
diagnostics$min_tail_ess >= sensitivity_mcmc_thresholds$min_tail_ess &&
identical(as.integer(diagnostics$max_treedepth_hits), 0L) &&
isTRUE(diagnostics$state_passes) &&
identical(as.integer(diagnostics$state_invalid_draws), 0L) &&
identical(as.integer(diagnostics$state_non_finite_draws), 0L) &&
identical(as.integer(diagnostics$state_invalid_cells), 0L) &&
is.finite(diagnostics$state_max_raw_harvest) &&
diagnostics$state_max_raw_harvest <=
state$hrate_limit + state$hrate_tolerance &&
is.finite(diagnostics$state_min_number) &&
diagnostics$state_min_number > state$number_tolerance &&
is.finite(diagnostics$state_max_harvest_penalty) &&
abs(diagnostics$state_max_harvest_penalty) <= state$penalty_tolerance
}
sensitivity_mcmc_passes <- function(diagnostics) {
sensitivity_mcmc_standard_passes(diagnostics) ||
sensitivity_mcmc_divergence_exception_passes(diagnostics)
}
sensitivity_mcmc_acceptance_label <- function(diagnostics) {
if (sensitivity_mcmc_standard_passes(diagnostics)) return("Pass")
if (sensitivity_mcmc_divergence_exception_passes(diagnostics)) {
return("Accepted (2 divergences)")
}
"Fail"
}
sensitivity_mcmc_identity_matches <- function(saved, expected, fit,
source_fit) {
if (!is.list(saved) || !is.list(saved$config) ||
!is.list(expected) || !is.list(expected$config)) {
return(FALSE)
}
normalized <- saved
binding_field <- "fit_identity_signature"
expected_binding <- expected$config[[binding_field]]
if (!is.null(expected_binding)) {
source_identity <-
source_fit$provenance$metadata$workflow_run_identity %||%
source_fit$provenance$metadata$run_identity
source_mcmc_binding <-
source_fit$provenance$metadata$mcmc_source_fit_identity_signature
saved_binding <- saved$config[[binding_field]]
if (!is.list(source_identity) ||
!is.character(source_identity$signature) ||
length(source_identity$signature) != 1L ||
!is.character(saved_binding) || length(saved_binding) != 1L ||
is.na(saved_binding) ||
!saved_binding %in% c(
source_identity$signature,
source_mcmc_binding,
source_fit$provenance$metadata$
mcmc_source_fit_workflow_signature,
expected_binding
)) {
return(FALSE)
}
normalized$config[[binding_field]] <- expected_binding
}
esc31_mcmc_identity_equivalent(
normalized,
expected,
fit = fit,
base_fit = base_fit
)
}
read_sensitivity_mcmc_fit <- function(file, source_fit, expected_identity) {
if (!file.exists(file)) return(NULL)
expected_key <- expected_identity$config$model_key
if (length(expected_key) != 1L || is.na(expected_key) ||
!expected_key %in% valid_sensitivity_mcmc_keys) return(NULL)
fit <- tryCatch(
sbt_fit_read(
file,
strict = TRUE,
rebuild = FALSE,
verify_validation = FALSE
),
error = function(e) NULL
)
if (is.null(fit) || is.null(fit$mcmc)) return(NULL)
diagnostics <- fit$fit$diagnostics$mcmc
saved_state <- fit$fit$diagnostics$biological_state_posterior
saved_identity <- fit$provenance$metadata$mcmc_run_identity
cache_checks <- c(
scientific_signature = identical(
esc31_fit_scientific_signature(fit),
esc31_fit_scientific_signature(source_fit)
),
mcmc_identity = sensitivity_mcmc_identity_matches(
saved_identity,
expected_identity,
fit = fit,
source_fit = source_fit
),
diagnostic_shape =
is.data.frame(diagnostics) && nrow(diagnostics) == 1L,
diagnostic_run_signature =
is.data.frame(diagnostics) && nrow(diagnostics) == 1L &&
identical(
as.character(diagnostics$run_signature),
saved_identity$signature
),
sampler_settings = identical(
esc31_mcmc_sampler_configuration(fit),
expected_identity$config$sampler
),
acceptance_thresholds = identical(
esc31_mcmc_acceptance_thresholds(fit),
expected_identity$config$thresholds
)
)
if (!all(cache_checks)) {
if (identical(
tolower(trimws(Sys.getenv("ESC31_DEBUG_MCMC_CACHE", ""))),
"true"
)) {
message(
"Rejected sensitivity MCMC cache `", expected_key, "` before ",
"recalculation: ",
paste(names(cache_checks)[!cache_checks], collapse = ", "),
"."
)
if (!cache_checks[["mcmc_identity"]] &&
is.list(saved_identity) && is.list(saved_identity$config)) {
debug_identity <- saved_identity
for (field in c(
"base_fit_scientific_signature", "fit_scientific_signature"
)) {
if (!is.null(expected_identity$config[[field]])) {
debug_identity$config[[field]] <-
expected_identity$config[[field]]
}
}
message(
"MCMC scientific-identity difference for `", expected_key, "`: ",
paste(
as.character(all.equal(
esc31_workflow_scientific_identity(debug_identity),
esc31_workflow_scientific_identity(expected_identity),
check.attributes = FALSE
)),
collapse = "; "
),
"."
)
}
}
return(NULL)
}
recalculated_fit <- tryCatch(
check_mcmc(
fit,
cores = sensitivity_state_diagnostic_cores,
stop_on_failure = FALSE,
mode = sensitivity_fit_validation_mode
),
error = function(e) NULL
)
validation_record <- if (is.null(recalculated_fit)) {
NULL
} else {
tryCatch(
sbt_fit_validation(
recalculated_fit,
scope = "mcmc",
verify = !identical(sensitivity_fit_validation_mode, "skip")
),
error = function(e) NULL
)
}
if (is.null(recalculated_fit) || is.null(validation_record) ||
!isTRUE(all.equal(
unclass(diagnostics),
unclass(recalculated_fit$fit$diagnostics$mcmc),
check.attributes = FALSE,
tolerance = 1e-12
)) ||
!isTRUE(all.equal(
saved_state,
recalculated_fit$fit$diagnostics$biological_state_posterior,
check.attributes = FALSE,
tolerance = 1e-10
))) {
if (identical(
tolower(trimws(Sys.getenv("ESC31_DEBUG_MCMC_CACHE", ""))),
"true"
)) {
message(
"Rejected sensitivity MCMC cache `", expected_key,
"` after embedded fit validation."
)
}
return(NULL)
}
recalculated_fit
}
sensitivity_dense_metric_qinv <- function(key, source_fit, sampling_object,
config) {
precision_source <- config$dense_metric_precision_source %||% NULL
if (is.null(precision_source)) return(NULL)
expected_source <- "RTMB exact AD Hessian at source MLE"
if (!identical(key, "no_pop_hsp") ||
!identical(config$metric, "dense") ||
!identical(precision_source, expected_source)) {
stop(
"The exact-Hessian dense metric is approved only for `no_pop_hsp`.",
call. = FALSE
)
}
sampling_par <- source_fit$fit$opt$par
if (length(sampling_par) != length(sampling_object$par) ||
!identical(names(sampling_par), names(sampling_object$par))) {
stop(
"The exact-Hessian metric parameters do not match the source MLE.",
call. = FALSE
)
}
sampling_object$par <- sampling_par
sampling_object$env$last.par.best <- sampling_par
objective_at_mle <- sampling_object$fn(sampling_par)
if (!isTRUE(all.equal(
as.numeric(objective_at_mle),
as.numeric(source_fit$fit$opt$objective),
tolerance = 1e-10
))) {
stop(
"The exact-Hessian metric object does not reproduce the source MLE.",
call. = FALSE
)
}
dense_hessian <- sampling_object$he(sampling_par)
if (!is.matrix(dense_hessian) ||
!identical(dim(dense_hessian), rep(length(sampling_par), 2L)) ||
any(!is.finite(dense_hessian))) {
stop("The exact dense-metric Hessian is invalid.", call. = FALSE)
}
dense_hessian <- 0.5 * (dense_hessian + t(dense_hessian))
dense_chol <- tryCatch(
chol(dense_hessian),
error = function(e) {
stop(
"The exact dense-metric Hessian is not positive definite: ",
conditionMessage(e),
call. = FALSE
)
}
)
dense_qinv <- chol2inv(dense_chol)
dense_qinv <- 0.5 * (dense_qinv + t(dense_qinv))
if (any(!is.finite(dense_qinv)) || any(diag(dense_qinv) <= 0)) {
stop("The exact dense-metric covariance is invalid.", call. = FALSE)
}
tryCatch(
chol(dense_qinv),
error = function(e) {
stop(
"The exact dense-metric covariance is not positive definite: ",
conditionMessage(e),
call. = FALSE
)
}
)
dense_qinv
}
ensure_sensitivity_mcmc_fit <- function(key, source_fit, file,
fit_identity_signature) {
if (length(key) != 1L || is.na(key) ||
!key %in% valid_sensitivity_mcmc_keys) {
stop("MCMC cannot be launched for a non-executable sensitivity key.",
call. = FALSE)
}
if (run_sensitivity_mcmc && !sensitivity_mcmc_requested(key)) return(NULL)
if (identical(key, "base")) {
return(esc31_load_base_mcmc(
file,
source_fit,
esc_dir,
validation_mode = sensitivity_fit_validation_mode
)$fit)
}
config <- sensitivity_mcmc_config(key)
expected_identity <- sensitivity_mcmc_identity(
key,
source_fit,
fit_identity_signature,
config
)
cached <- read_sensitivity_mcmc_fit(file, source_fit, expected_identity)
if (!is.null(cached) &&
(!sensitivity_mcmc_requested(key) ||
sensitivity_mcmc_passes(cached$fit$diagnostics$mcmc))) {
return(cached)
}
if (!sensitivity_mcmc_requested(key)) return(NULL)
if (!sensitivity_mle_passes(source_fit)) {
stop(
"The source MLE does not pass the convergence, maximum-gradient, and ",
"estimability gates for `", key, "`.",
call. = FALSE
)
}
if (!identical(key, "base") &&
!sensitivity_optimizer_provenance_passes(source_fit)) {
stop(
"The source sensitivity MLE does not match the current optimizer ",
"contract for `", key, "`.",
call. = FALSE
)
}
sampling_object <- sbt_obj(source_fit, fresh = TRUE)
dense_metric_qinv <- sensitivity_dense_metric_qinv(
key,
source_fit,
sampling_object,
config
)
mcmc_fit <- sbt_mcmc(
source_fit,
init = config$init,
check = FALSE,
num_samples = config$num_samples,
num_warmup = config$num_warmup,
chains = config$chains,
cores = config$cores,
metric = config$metric,
Qinv = dense_metric_qinv,
seed = config$seed,
skip_optimization = config$skip_optimization,
control = list(
adapt_delta = config$adapt_delta,
max_treedepth = config$max_treedepth
),
refresh = config$refresh
)
mcmc_fit$provenance$metadata <- utils::modifyList(
mcmc_fit$provenance$metadata,
list(
mcmc_run_identity = expected_identity,
mcmc_source_fit_identity_signature = fit_identity_signature,
mcmc_source_fit_workflow_signature = fit_identity_signature
)
)
# Preserve the assessment-specific biological-state contract together with
# the five scalar sampler gates that `check_mcmc()` evaluates. Without this
# binding, a newly sampled fit would pass its first check but appear stale
# when the complete threshold identity is verified on the next load.
mcmc_fit$mcmc$settings$acceptance_thresholds <-
sensitivity_mcmc_thresholds
source_mle_validation <- sbt_fit_validation(
source_fit,
scope = "mle",
require_pass = TRUE
)
if (!identical(
mcmc_fit$validation$records$mle,
source_mle_validation
) ||
!identical(esc31_mcmc_sampler_configuration(mcmc_fit), config) ||
!identical(
esc31_mcmc_acceptance_thresholds(mcmc_fit),
sensitivity_mcmc_thresholds
)) {
stop(
"The staged MCMC result does not preserve its source-MLE validation ",
"and production sampler contract for `", key, "`.",
call. = FALSE
)
}
sbt_fit_validate(mcmc_fit)
# Preserve an expensive sampler result before running posterior diagnostics.
sbt_fit_save(mcmc_fit, file, overwrite = TRUE)
mcmc_fit <- check_mcmc(
mcmc_fit,
cores = sensitivity_state_diagnostic_cores,
stop_on_failure = FALSE,
mode = "auto"
)
diagnostics <- mcmc_fit$fit$diagnostics$mcmc
sbt_fit_save(mcmc_fit, file, overwrite = TRUE)
if (!sensitivity_mcmc_passes(diagnostics)) {
stop(
"Sensitivity MCMC diagnostics do not satisfy the acceptance decision for `",
key, "`; the saved fit is retained for inspection and will be resampled ",
"on the next requested run.",
call. = FALSE
)
}
mcmc_fit
}
sensitivity_workflow_identity <- function(sensitivity_id, scientific_config) {
if (!is.character(sensitivity_id) || length(sensitivity_id) != 1L ||
is.na(sensitivity_id) || !nzchar(sensitivity_id)) {
stop("`sensitivity_id` must be one non-empty string.", call. = FALSE)
}
if (!is.list(scientific_config) || is.null(names(scientific_config)) ||
any(!nzchar(names(scientific_config)))) {
stop("`scientific_config` must be a named list.", call. = FALSE)
}
esc31_workflow_identity(
esc_dir,
workflow_files = character(),
sbt_entry_points = esc31_sbt_entry_points("sensitivities"),
config = list(
stage = "sensitivity",
fit_workflow_version =
"esc31_mle_sensitivity_v12_length_base_standard_ll4_ll3_hyper",
sensitivity_id = sensitivity_id,
base_fit_scientific_signature = esc31_fit_scientific_signature(base_fit),
base_run_signature =
base_fit$provenance$metadata$run_identity$signature,
harvest_wall_contract = sensitivity_harvest_wall_contract,
harvest_wall_decision = sensitivity_harvest_wall_decision,
harvest_wall_recomputation_tolerance =
sensitivity_harvest_wall_recomputation_tolerance,
optimizer = sensitivity_optimizer_contract,
mle_acceptance_thresholds = sensitivity_mle_thresholds,
scientific_config = scientific_config,
fit_seed = NA_integer_
)
)
}
sensitivity_metadata <- function(sensitivity_id, scientific_config, data_sens,
parameter_overrides = list(),
estimated_parameters = character(),
reference_obj = sbt_obj(base_fit)) {
setup <- prepare_sensitivity_model(
data_sens = data_sens,
parameter_overrides = parameter_overrides,
estimated_parameters = estimated_parameters,
reference_obj = reference_obj
)
scientific_config <- utils::modifyList(
scientific_config,
setup$scientific_config
)
metadata <- utils::modifyList(
list(
sensitivity_id = sensitivity_id,
base_fit_scientific_signature = esc31_fit_scientific_signature(base_fit),
base_run_signature =
base_fit$provenance$metadata$run_identity$signature,
optimizer_contract_signature = sensitivity_optimizer_contract_signature,
workflow_run_identity = sensitivity_workflow_identity(
sensitivity_id,
scientific_config
)
),
scientific_config
)
attr(metadata, "expected_model_data") <- setup$data
metadata
}
sensitivity_workflow_identity_equivalent <- function(saved, expected) {
if (!is.list(saved) || !is.list(expected) ||
!is.list(saved$config) || !is.list(expected$config)) {
return(FALSE)
}
base_signature_field <- "base_fit_scientific_signature"
if (!is.null(expected$config[[base_signature_field]])) {
if (!esc31_saved_fit_signature_matches(
saved$config[[base_signature_field]], base_fit)) {
return(FALSE)
}
saved$config[[base_signature_field]] <-
expected$config[[base_signature_field]]
}
scientific_identity <- function(identity) {
identity <- esc31_workflow_scientific_identity(identity)
if (is.null(identity)) return(NULL)
config <- identity$config$scientific_config
if (is.list(config)) {
config$transformed_data_md5 <- NULL
config$model_data_md5 <- NULL
# The seed records how optimisation was launched; it is not part of the
# converged MLE or the payload-bound posterior identity.
config$seeded_parameters_md5 <- NULL
identity$config$scientific_config <- config
}
identity
}
saved_scientific <- scientific_identity(saved)
expected_scientific <- scientific_identity(expected)
!is.null(saved_scientific) && !is.null(expected_scientific) &&
identical(saved_scientific, expected_scientific)
}
sensitivity_model_data_equivalent <- function(saved, expected) {
isTRUE(all.equal(
saved,
expected,
# Reconstructing the unchanged inputs after their move to sbtdata can
# differ at floating-point round-off only (approximately 2e-16 here).
tolerance = 1e-14,
check.attributes = TRUE
))
}
sensitivity_metadata_matches <- function(fit, expected) {
if (!has_sensitivity_fit(fit)) return(FALSE)
metadata <- fit$provenance$metadata
expected_model_data <- attr(expected, "expected_model_data", exact = TRUE)
!is.null(expected_model_data) &&
sensitivity_model_data_equivalent(fit$data, expected_model_data) &&
all(vapply(names(expected), function(name) {
if (identical(name, "workflow_run_identity")) {
sensitivity_workflow_identity_equivalent(
metadata[[name]],
expected[[name]]
)
} else if (identical(name, "base_fit_scientific_signature")) {
esc31_saved_fit_signature_matches(metadata[[name]], base_fit)
} else if (name %in% c(
"transformed_data_md5", "model_data_md5", "seeded_parameters_md5"
)) {
TRUE
} else {
identical(metadata[[name]], expected[[name]])
}
}, logical(1)))
}
read_saved_sensitivity_fit <- function(file, expected_metadata, required_diagnostics = character()) {
if (!file.exists(file)) return(NULL)
fit <- tryCatch(
sbt_fit_read(
file,
strict = TRUE,
rebuild = FALSE,
verify_validation = FALSE
),
error = function(e) NULL
)
initial_checks <- c(
readable = !is.null(fit),
metadata = !is.null(fit) &&
sensitivity_metadata_matches(fit, expected_metadata),
optimizer_provenance = !is.null(fit) &&
sensitivity_optimizer_provenance_passes(fit),
mle = !is.null(fit) && sensitivity_mle_passes(fit)
)
if (!all(initial_checks)) {
if (identical(
tolower(trimws(Sys.getenv("ESC31_DEBUG_MCMC_CACHE", ""))),
"true"
)) {
message(
"Rejected saved sensitivity fit `", basename(file), "`: ",
paste(names(initial_checks)[!initial_checks], collapse = ", "),
"."
)
}
return(NULL)
}
optimizer_diagnostics <- fit$fit$diagnostics
if (!identical(
optimizer_diagnostics$optimizer_contract_signature,
sensitivity_optimizer_contract_signature
) ||
!identical(
optimizer_diagnostics$optimizer_contract,
sensitivity_optimizer_contract
) ||
!identical(
esc31_object_md5(optimizer_diagnostics$optimizer_contract),
sensitivity_optimizer_contract_signature
)) return(NULL)
if (length(setdiff(required_diagnostics, names(fit$fit$diagnostics)))) return(NULL)
compatible_error <- NULL
compatible <- tryCatch(
{
sbt_obj(fit)
TRUE
},
error = function(e) {
compatible_error <<- conditionMessage(e)
FALSE
}
)
if (!compatible && identical(
tolower(trimws(Sys.getenv("ESC31_DEBUG_MCMC_CACHE", ""))),
"true"
)) {
message(
"Rejected saved sensitivity fit `", basename(file),
"` during objective reconstruction: ", compatible_error, "."
)
}
if (compatible) fit else NULL
}
seed_sensitivity_parameters <- function(parameters_sens, reference_obj) {
reference_parameters <- reference_obj$env$parList(reference_obj$env$last.par.best)
for (name in intersect(names(parameters_sens), names(reference_parameters))) {
current <- parameters_sens[[name]]
fitted <- reference_parameters[[name]]
if (is.matrix(current) && is.matrix(fitted)) {
current_rows <- rownames(current)
fitted_rows <- rownames(fitted)
current_cols <- colnames(current)
fitted_cols <- colnames(fitted)
if (!is.null(current_rows) && !is.null(fitted_rows) &&
!is.null(current_cols) && !is.null(fitted_cols)) {
shared_rows <- intersect(current_rows, fitted_rows)
shared_cols <- intersect(current_cols, fitted_cols)
current[shared_rows, shared_cols] <- fitted[shared_rows, shared_cols, drop = FALSE]
} else {
shared_rows <- seq_len(min(nrow(current), nrow(fitted)))
shared_cols <- seq_len(min(ncol(current), ncol(fitted)))
current[shared_rows, shared_cols] <- fitted[shared_rows, shared_cols, drop = FALSE]
}
parameters_sens[[name]] <- current
} else if (length(fitted) == 1L && length(current) > 1L) {
parameters_sens[[name]][] <- as.numeric(fitted)
} else {
n <- min(length(current), length(fitted))
parameters_sens[[name]][seq_len(n)] <- as.numeric(fitted)[seq_len(n)]
}
}
parameters_sens
}
sensitivity_map_contract <- list(
generator = "sbt::get_map",
ll4_selectivity = "one_active_1_x_14_block_starting_1953",
ll4_hyperparameters = "fixed_values_inherited_from_ll3",
mortality = "map follows M_switch; length exponent fixed at minus one"
)
prepare_sensitivity_model <- function(data_sens, parameter_overrides = list(),
estimated_parameters = character(),
reference_obj = sbt_obj(base_fit)) {
if (!is.list(parameter_overrides) ||
(length(parameter_overrides) &&
(is.null(names(parameter_overrides)) ||
any(!nzchar(names(parameter_overrides)))))) {
stop("`parameter_overrides` must be a named list.", call. = FALSE)
}
if (!is.character(estimated_parameters) || anyNA(estimated_parameters) ||
any(!nzchar(estimated_parameters)) || anyDuplicated(estimated_parameters)) {
stop("`estimated_parameters` must contain unique parameter names.",
call. = FALSE)
}
transformed_data_md5 <- esc31_object_md5(data_sens)
parameters_sens <- get_parameters(data = data_sens)
parameters_sens <- seed_sensitivity_parameters(parameters_sens, reference_obj)
for (name in names(parameter_overrides)) {
if (!name %in% names(parameters_sens)) {
stop("Unknown sensitivity parameter override: ", name, call. = FALSE)
}
parameters_sens[[name]] <- parameter_overrides[[name]]
}
map_sens <- get_map(parameters = parameters_sens)
for (name in estimated_parameters) {
if (!name %in% names(map_sens) || !all(is.na(map_sens[[name]]))) {
stop(
"An explicitly estimated sensitivity parameter must be fixed by the ",
"default map: ", name, ".",
call. = FALSE
)
}
map_sens[[name]] <- NULL
}
data_sens$priors <- get_priors(parameters = parameters_sens)
scientific_config <- list(
transformed_data_md5 = transformed_data_md5,
model_data_md5 = esc31_object_md5(data_sens),
seeded_parameters_md5 = esc31_object_md5(parameters_sens),
parameter_overrides = parameter_overrides,
map_contract = sensitivity_map_contract,
parameter_map_md5 = esc31_object_md5(map_sens)
)
if (length(estimated_parameters)) {
scientific_config$estimated_parameters <- estimated_parameters
}
list(
data = data_sens,
parameters = parameters_sens,
map = map_sens,
scientific_config = scientific_config
)
}
sensitivity_point_diagnostics <- function(obj, par, bounds,
bound_tolerance = 0) {
expected_names <- names(obj$par)
finite_parameters <- is.numeric(par) &&
length(par) == length(expected_names) &&
identical(names(par), expected_names) &&
all(is.finite(par))
within_bounds <- finite_parameters &&
all(par >= bounds$lower - bound_tolerance) &&
all(par <= bounds$upper + bound_tolerance)
objective <- if (finite_parameters) {
tryCatch(as.numeric(obj$fn(par)), error = function(e) Inf)
} else {
Inf
}
finite_objective <- length(objective) == 1L && is.finite(objective)
if (!finite_objective) objective <- Inf
gradient <- if (finite_parameters && finite_objective) {
tryCatch(as.numeric(obj$gr(par)), error = function(e) numeric())
} else {
numeric()
}
finite_gradient <- length(gradient) == length(par) &&
all(is.finite(gradient))
max_abs_gradient <- if (finite_gradient) max(abs(gradient)) else Inf
gradient_l2 <- if (finite_gradient) sqrt(sum(gradient^2)) else Inf
list(
valid = finite_parameters && within_bounds && finite_objective &&
finite_gradient,
finite_parameters = finite_parameters,
within_bounds = within_bounds,
finite_objective = finite_objective,
finite_gradient = finite_gradient,
objective = objective,
gradient = gradient,
max_abs_gradient = max_abs_gradient,
gradient_l2 = gradient_l2
)
}
sensitivity_optimizer_certified <- function(opt, point_diagnostics,
gradient_tolerance) {
is.list(opt) &&
length(opt$convergence) == 1L &&
is.finite(opt$convergence) &&
identical(as.integer(opt$convergence), 0L) &&
isTRUE(point_diagnostics$valid) &&
is.finite(point_diagnostics$max_abs_gradient) &&
point_diagnostics$max_abs_gradient <= gradient_tolerance
}
sensitivity_optimizer_progress <- function(previous_par, previous_point,
candidate_par, candidate_point) {
previous_objective <- previous_point$objective %||% NA_real_
candidate_objective <- candidate_point$objective %||% NA_real_
relative_objective_decrease <- if (
length(previous_objective) == 1L && is.finite(previous_objective) &&
length(candidate_objective) == 1L && is.finite(candidate_objective)
) {
(previous_objective - candidate_objective) /
max(1, abs(previous_objective))
} else {
NA_real_
}
previous_gradient <- previous_point$max_abs_gradient %||% NA_real_
candidate_gradient <- candidate_point$max_abs_gradient %||% NA_real_
gradient_ratio <- if (
length(previous_gradient) == 1L && is.finite(previous_gradient) &&
previous_gradient >= 0 &&
length(candidate_gradient) == 1L && is.finite(candidate_gradient) &&
candidate_gradient >= 0
) {
candidate_gradient / max(previous_gradient, .Machine$double.eps)
} else {
NA_real_
}
relative_parameter_step <- if (
is.numeric(previous_par) && is.numeric(candidate_par) &&
length(previous_par) == length(candidate_par) &&
identical(names(previous_par), names(candidate_par)) &&
all(is.finite(previous_par)) && all(is.finite(candidate_par))
) {
max(abs(candidate_par - previous_par) / pmax(1, abs(previous_par)))
} else {
NA_real_
}
list(
previous_objective = as.numeric(previous_objective),
relative_objective_decrease = relative_objective_decrease,
gradient_ratio = gradient_ratio,
relative_parameter_step = relative_parameter_step
)
}
fit_sensitivity <- function(
data_sens, sensitivity_id, label,
reference_obj = sbt_obj(base_fit),
max_passes = sensitivity_optimizer_contract$primary$max_passes,
parameter_overrides = list(),
estimated_parameters = character(),
gradient_tolerance = sensitivity_optimizer_contract$acceptance$max_gradient,
diagnostics = list(), metadata = list()) {
attr(metadata, "expected_model_data") <- NULL
required_metadata <- c(
"sensitivity_id", "base_fit_scientific_signature", "base_run_signature",
"optimizer_contract_signature", "workflow_run_identity"
)
if (!all(required_metadata %in% names(metadata)) ||
!identical(metadata$sensitivity_id, sensitivity_id) ||
!identical(
metadata$optimizer_contract_signature,
sensitivity_optimizer_contract_signature
) ||
!identical(
metadata$workflow_run_identity$config$optimizer,
sensitivity_optimizer_contract
) ||
!identical(
metadata$workflow_run_identity$config$harvest_wall_contract,
sensitivity_harvest_wall_contract
) ||
!identical(
metadata$workflow_run_identity$config$harvest_wall_decision,
sensitivity_harvest_wall_decision
) ||
!identical(
metadata$workflow_run_identity$config$harvest_wall_recomputation_tolerance,
sensitivity_harvest_wall_recomputation_tolerance
)) {
stop("`metadata` must be the matching per-sensitivity metadata.", call. = FALSE)
}
if (length(max_passes) != 1L || !is.numeric(max_passes) ||
!is.finite(max_passes) || max_passes != as.integer(max_passes) ||
!identical(
as.integer(max_passes),
sensitivity_optimizer_contract$primary$max_passes
)) {
stop(
"`max_passes` must match the hashed sensitivity optimizer contract.",
call. = FALSE
)
}
if (length(gradient_tolerance) != 1L ||
!is.numeric(gradient_tolerance) || !is.finite(gradient_tolerance) ||
gradient_tolerance <= 0 ||
!identical(
as.numeric(gradient_tolerance),
as.numeric(sensitivity_optimizer_contract$acceptance$max_gradient)
)) {
stop(
"`gradient_tolerance` must match the hashed sensitivity optimizer contract.",
call. = FALSE
)
}
if (!is.list(diagnostics) || !is.list(metadata)) {
stop("`diagnostics` and `metadata` must be lists.", call. = FALSE)
}
reserved_diagnostics <- c(
"initial_nll", "final_nll", "max_gradient", "gradient_l2",
"optimizer_contract", "optimizer_contract_signature",
"optimizer_history", "optimizer_summary", "biological_state_mle",
"harvest_wall_mle", "esc31_harvest_wall_mle"
)
supplied_reserved <- intersect(names(diagnostics), reserved_diagnostics)
if (length(supplied_reserved)) {
stop(
"`diagnostics` may not replace optimizer diagnostics: ",
paste(supplied_reserved, collapse = ", "), ".",
call. = FALSE
)
}
setup <- prepare_sensitivity_model(
data_sens = data_sens,
parameter_overrides = parameter_overrides,
estimated_parameters = estimated_parameters,
reference_obj = reference_obj
)
setup_matches_metadata <- all(vapply(
names(setup$scientific_config),
function(name) identical(metadata[[name]], setup$scientific_config[[name]]),
logical(1)
))
if (!setup_matches_metadata) {
stop(
"The current transformed data, parameters, or map do not match metadata for `",
sensitivity_id, "`.",
call. = FALSE
)
}
data_sens <- setup$data
parameters_sens <- setup$parameters
map_sens <- setup$map
fit <- sbt_fit(
data_sens,
control = sensitivity_optimizer_contract$primary$control,
metadata = utils::modifyList(
list(
assessment = "ESC31",
label = label
),
metadata
),
optimizer = sensitivity_optimizer_contract$implementation
)
fit <- sbt_add_parameters(fit, parameters_sens)
fit <- sbt_add_map(fit, map_sens)
fit <- sbt_add_priors(fit)
fit <- sbt_build_object(fit)
obj_sens <- sbt_obj(fit)
bounds_sens <- fit$bounds
optimizer_start <- obj_sens$par
bound_tolerance <- sensitivity_optimizer_contract$acceptance$bound_tolerance
initial_diagnostics <- sensitivity_point_diagnostics(
obj_sens,
optimizer_start,
bounds_sens,
bound_tolerance = bound_tolerance
)
if (!isTRUE(initial_diagnostics$valid)) {
stop(
label,
" has non-finite or out-of-bounds starting parameters, objective, or gradient.",
call. = FALSE
)
}
initial_nll <- initial_diagnostics$objective
history <- data.frame(
pass = integer(),
stage = character(),
attempt = integer(),
objective = numeric(),
max_abs_gradient = numeric(),
gradient_l2 = numeric(),
convergence = integer(),
iterations = integer(),
function_evaluations = integer(),
gradient_evaluations = integer(),
elapsed_seconds = numeric(),
finite_parameters = logical(),
within_bounds = logical(),
finite_objective = logical(),
finite_gradient = logical(),
accepted_as_best = logical(),
previous_objective = numeric(),
relative_objective_decrease = numeric(),
gradient_ratio = numeric(),
relative_parameter_step = numeric(),
transition = character(),
transition_reason = character(),
scale_min = numeric(),
scale_max = numeric(),
nonpositive_hessian_diagonal = integer(),
hessian_reciprocal_condition = numeric(),
alpha = numeric(),
armijo_pass = logical(),
gradient_improvement = logical(),
message = character(),
stringsAsFactors = FALSE
)
append_optimizer_history <- function(
stage, attempt, point, result = NULL, accepted_as_best = FALSE,
elapsed_seconds = NA_real_, scale = NULL,
nonpositive_hessian_diagonal = NA_integer_,
hessian_reciprocal_condition = NA_real_, alpha = NA_real_,
armijo_pass = NA, gradient_improvement = NA,
progress = list(
previous_objective = NA_real_,
relative_objective_decrease = NA_real_,
gradient_ratio = NA_real_,
relative_parameter_step = NA_real_
),
transition = "pending", transition_reason = "not classified",
message = NULL) {
evaluations <- result$evaluations %||%
c("function" = NA_integer_, "gradient" = NA_integer_)
convergence <- result$convergence %||% NA_integer_
iterations <- result$iterations %||% NA_integer_
function_evaluations <- if (
!is.null(names(evaluations)) && "function" %in% names(evaluations)
) evaluations[["function"]] else NA_integer_
gradient_evaluations <- if (
!is.null(names(evaluations)) && "gradient" %in% names(evaluations)
) evaluations[["gradient"]] else NA_integer_
result_message <- paste(
as.character(message %||% result$message %||% stage),
collapse = "; "
)
history <<- rbind(history, data.frame(
pass = nrow(history),
stage = as.character(stage),
attempt = as.integer(attempt),
objective = as.numeric(point$objective),
max_abs_gradient = as.numeric(point$max_abs_gradient),
gradient_l2 = as.numeric(point$gradient_l2),
convergence = as.integer(convergence),
iterations = as.integer(iterations),
function_evaluations = as.integer(function_evaluations),
gradient_evaluations = as.integer(gradient_evaluations),
elapsed_seconds = as.numeric(elapsed_seconds),
finite_parameters = isTRUE(point$finite_parameters),
within_bounds = isTRUE(point$within_bounds),
finite_objective = isTRUE(point$finite_objective),
finite_gradient = isTRUE(point$finite_gradient),
accepted_as_best = isTRUE(accepted_as_best),
previous_objective = as.numeric(progress$previous_objective),
relative_objective_decrease = as.numeric(
progress$relative_objective_decrease
),
gradient_ratio = as.numeric(progress$gradient_ratio),
relative_parameter_step = as.numeric(
progress$relative_parameter_step
),
transition = as.character(transition),
transition_reason = as.character(transition_reason),
scale_min = if (is.null(scale)) NA_real_ else min(scale),
scale_max = if (is.null(scale)) NA_real_ else max(scale),
nonpositive_hessian_diagonal = as.integer(
nonpositive_hessian_diagonal
),
hessian_reciprocal_condition = as.numeric(
hessian_reciprocal_condition
),
alpha = as.numeric(alpha),
armijo_pass = as.logical(armijo_pass),
gradient_improvement = as.logical(gradient_improvement),
message = as.character(result_message),
stringsAsFactors = FALSE
))
invisible(NULL)
}
set_last_optimizer_transition <- function(transition, transition_reason) {
if (!nrow(history)) {
stop("Cannot classify an empty optimizer history.", call. = FALSE)
}
updated_history <- history
updated_history$transition[[nrow(updated_history)]] <- transition
updated_history$transition_reason[[nrow(updated_history)]] <-
transition_reason
history <<- updated_history
invisible(NULL)
}
append_optimizer_history(
stage = "start",
attempt = 0L,
point = initial_diagnostics,
accepted_as_best = TRUE,
transition = "seed",
transition_reason = "seeded sensitivity parameters",
message = "seeded sensitivity parameters"
)
best <- list(par = optimizer_start, diagnostics = initial_diagnostics)
opt_sens <- NULL
final_stage <- NA_character_
nlminb_calls <- 0L
primary_passes <- 0L
scaled_passes <- 0L
certification_passes <- 0L
newton_metadata <- list(
attempted = FALSE,
accepted = FALSE,
status = "not required"
)
run_nlminb_stage <- function(stage, attempt, stage_control,
use_hessian = FALSE, scale = NULL,
nonpositive_hessian_diagonal = NA_integer_) {
previous_best <- best
start <- previous_best$par
result <- NULL
timing <- system.time({
result <- tryCatch(
{
if (isTRUE(use_hessian)) {
nlminb(
start = start,
objective = obj_sens$fn,
gradient = obj_sens$gr,
hessian = obj_sens$he,
lower = bounds_sens$lower,
upper = bounds_sens$upper,
control = stage_control
)
} else if (!is.null(scale)) {
nlminb(
start = start,
objective = obj_sens$fn,
gradient = obj_sens$gr,
scale = scale,
lower = bounds_sens$lower,
upper = bounds_sens$upper,
control = stage_control
)
} else {
nlminb(
start = start,
objective = obj_sens$fn,
gradient = obj_sens$gr,
lower = bounds_sens$lower,
upper = bounds_sens$upper,
control = stage_control
)
}
},
error = function(e) e
)
})
nlminb_calls <<- nlminb_calls + 1L
if (inherits(result, "error")) {
failed <- list(
finite_parameters = FALSE,
within_bounds = FALSE,
finite_objective = FALSE,
finite_gradient = FALSE,
objective = Inf,
max_abs_gradient = Inf,
gradient_l2 = Inf
)
progress <- sensitivity_optimizer_progress(
previous_best$par,
previous_best$diagnostics,
NULL,
failed
)
append_optimizer_history(
stage = stage,
attempt = attempt,
point = failed,
elapsed_seconds = timing[["elapsed"]],
scale = scale,
nonpositive_hessian_diagonal = nonpositive_hessian_diagonal,
progress = progress,
transition = "retain_previous_best",
transition_reason = "optimizer_error",
message = conditionMessage(result)
)
return(list(result = NULL, diagnostics = failed, accepted = FALSE))
}
point <- sensitivity_point_diagnostics(
obj_sens,
result$par,
bounds_sens,
bound_tolerance = bound_tolerance
)
if (isTRUE(point$finite_objective)) result$objective <- point$objective
progress <- sensitivity_optimizer_progress(
previous_best$par,
previous_best$diagnostics,
result$par,
point
)
accepted <- isTRUE(point$valid) &&
point$objective <= previous_best$diagnostics$objective
if (accepted) {
best <<- list(par = result$par, diagnostics = point)
}
default_transition_reason <- if (accepted) {
"finite_nonworsening_candidate"
} else if (!isTRUE(point$valid)) {
"invalid_or_nonfinite_candidate"
} else {
"objective_worsened"
}
append_optimizer_history(
stage = stage,
attempt = attempt,
point = point,
result = result,
accepted_as_best = accepted,
elapsed_seconds = timing[["elapsed"]],
scale = scale,
nonpositive_hessian_diagonal = nonpositive_hessian_diagonal,
progress = progress,
transition = if (accepted) {
"stage_candidate_retained"
} else {
"retain_previous_best"
},
transition_reason = default_transition_reason
)
list(result = result, diagnostics = point, accepted = accepted)
}
for (attempt in seq_len(max_passes)) {
primary_passes <- primary_passes + 1L
previous_best <- best
primary_fit <- NULL
primary_timing <- system.time({
primary_fit <- tryCatch(
sbt_optimise(
fit,
n_passes = 1L,
control = sensitivity_optimizer_contract$primary$control,
check = FALSE
),
error = function(error) error
)
})
nlminb_calls <- nlminb_calls + 1L
if (inherits(primary_fit, "error")) {
failed <- list(
finite_parameters = FALSE,
within_bounds = FALSE,
finite_objective = FALSE,
finite_gradient = FALSE,
objective = Inf,
max_abs_gradient = Inf,
gradient_l2 = Inf
)
append_optimizer_history(
stage = "primary bounded nlminb",
attempt = attempt,
point = failed,
elapsed_seconds = primary_timing[["elapsed"]],
transition = "retain_previous_best",
transition_reason = "optimizer_error",
message = conditionMessage(primary_fit)
)
stage_result <- list(
result = NULL,
diagnostics = failed,
accepted = FALSE
)
} else {
fit <- primary_fit
obj_sens <- sbt_obj(fit)
bounds_sens <- fit$bounds
primary_opt <- fit$fit$opt
primary_point <- sensitivity_point_diagnostics(
obj_sens,
primary_opt$par,
bounds_sens,
bound_tolerance = bound_tolerance
)
if (isTRUE(primary_point$finite_objective)) {
primary_opt$objective <- primary_point$objective
}
primary_progress <- sensitivity_optimizer_progress(
previous_best$par,
previous_best$diagnostics,
primary_opt$par,
primary_point
)
primary_accepted <- isTRUE(primary_point$valid) &&
primary_point$objective <= previous_best$diagnostics$objective
if (primary_accepted) {
best <- list(par = primary_opt$par, diagnostics = primary_point)
}
append_optimizer_history(
stage = "primary bounded nlminb",
attempt = attempt,
point = primary_point,
result = primary_opt,
accepted_as_best = primary_accepted,
elapsed_seconds = primary_timing[["elapsed"]],
progress = primary_progress,
transition = if (primary_accepted) {
"stage_candidate_retained"
} else {
"retain_previous_best"
},
transition_reason = if (primary_accepted) {
"finite_nonworsening_candidate"
} else if (!isTRUE(primary_point$valid)) {
"invalid_or_nonfinite_candidate"
} else {
"objective_worsened"
},
message = paste(
"sbt_optimise primary pass:",
primary_opt$message %||% "completed"
)
)
stage_result <- list(
result = primary_opt,
diagnostics = primary_point,
accepted = primary_accepted
)
}
primary_certified <- isTRUE(stage_result$accepted) &&
sensitivity_optimizer_certified(
stage_result$result,
stage_result$diagnostics,
gradient_tolerance
)
primary_transition <- sensitivity_optimizer_contract$primary$transition
if (primary_certified) {
set_last_optimizer_transition(
primary_transition$certified,
"strict_numeric_certification_passed"
)
opt_sens <- stage_result$result
final_stage <- "primary bounded nlminb"
break
} else if (isTRUE(stage_result$accepted)) {
normal_convergence <- is.list(stage_result$result) &&
length(stage_result$result$convergence) == 1L &&
is.numeric(stage_result$result$convergence) &&
is.finite(stage_result$result$convergence) &&
identical(as.integer(stage_result$result$convergence), 0L)
set_last_optimizer_transition(
primary_transition$finite_nonworsening_uncertified,
if (normal_convergence) {
"maximum_gradient_gate_failed"
} else {
"nonzero_numeric_convergence"
}
)
} else {
set_last_optimizer_transition(
primary_transition$invalid_or_worsening,
history$transition_reason[[nrow(history)]]
)
}
}
if (is.null(opt_sens)) {
scaled_contract <- sensitivity_optimizer_contract$scaled
for (attempt in seq_len(scaled_contract$max_passes)) {
scaling_hessian <- tryCatch(
obj_sens$he(best$par),
error = function(e) e
)
valid_hessian <- is.matrix(scaling_hessian) &&
identical(dim(scaling_hessian), c(length(best$par), length(best$par))) &&
all(is.finite(scaling_hessian))
if (!valid_hessian) {
hessian_message <- if (inherits(scaling_hessian, "error")) {
conditionMessage(scaling_hessian)
} else {
"The scaling Hessian was non-finite or had the wrong dimensions."
}
append_optimizer_history(
stage = "scaled Hessian unavailable",
attempt = attempt,
point = best$diagnostics,
transition = "advance_to_newton",
transition_reason = "invalid_or_nonfinite_scaling_hessian",
message = hessian_message
)
break
}
hessian_diagonal <- diag(scaling_hessian)
parameter_scale <- 1 / sqrt(pmax(
abs(hessian_diagonal),
scaled_contract$hessian_diagonal_floor
))
parameter_scale <- pmin(
pmax(parameter_scale, scaled_contract$scale_min),
scaled_contract$scale_max
)
scaled_passes <- scaled_passes + 1L
stage_result <- run_nlminb_stage(
stage = "scaled bounded gradient nlminb",
attempt = attempt,
stage_control = scaled_contract$control,
scale = parameter_scale,
nonpositive_hessian_diagonal = sum(hessian_diagonal <= 0)
)
scaled_certified <- isTRUE(stage_result$accepted) &&
sensitivity_optimizer_certified(
stage_result$result,
stage_result$diagnostics,
gradient_tolerance
)
if (scaled_certified) {
set_last_optimizer_transition(
"accept",
"strict_numeric_certification_passed"
)
opt_sens <- stage_result$result
final_stage <- "scaled bounded gradient nlminb"
break
} else if (isTRUE(stage_result$accepted)) {
set_last_optimizer_transition(
"continue_scaled",
"finite_nonworsening_uncertified"
)
} else {
set_last_optimizer_transition(
"retry_scaled_from_previous_best",
history$transition_reason[[nrow(history)]]
)
}
}
}
if (is.null(opt_sens) && isTRUE(sensitivity_optimizer_contract$newton$enabled)) {
newton_contract <- sensitivity_optimizer_contract$newton
newton_metadata <- list(
attempted = TRUE,
accepted = FALSE,
status = "Hessian validation pending"
)
pre_newton <- best$diagnostics
newton_hessian <- tryCatch(
obj_sens$he(best$par),
error = function(e) e
)
valid_newton_hessian <- is.matrix(newton_hessian) &&
identical(dim(newton_hessian), c(length(best$par), length(best$par))) &&
all(is.finite(newton_hessian))
if (valid_newton_hessian) {
newton_hessian <- (newton_hessian + t(newton_hessian)) / 2
hessian_cholesky <- tryCatch(
chol(newton_hessian),
error = function(e) e
)
hessian_rcond <- tryCatch(
as.numeric(rcond(newton_hessian)),
error = function(e) NA_real_
)
hessian_acceptable <- !inherits(hessian_cholesky, "error") &&
length(hessian_rcond) == 1L && is.finite(hessian_rcond) &&
hessian_rcond >= newton_contract$minimum_reciprocal_condition
} else {
hessian_cholesky <- NULL
hessian_rcond <- NA_real_
hessian_acceptable <- FALSE
}
if (!hessian_acceptable) {
newton_metadata$status <- paste0(
"skipped: Hessian was not finite, SPD, and sufficiently conditioned; rcond=",
format(hessian_rcond, digits = 6)
)
append_optimizer_history(
stage = "Newton correction skipped",
attempt = 0L,
point = pre_newton,
hessian_reciprocal_condition = hessian_rcond,
transition = "stop_newton",
transition_reason = "hessian_not_spd_or_sufficiently_conditioned",
message = newton_metadata$status
)
} else {
newton_step <- tryCatch(
as.numeric(backsolve(
hessian_cholesky,
forwardsolve(
t(hessian_cholesky),
matrix(pre_newton$gradient, ncol = 1L)
)
)),
error = function(e) numeric()
)
direction <- -newton_step
descent_slope <- if (length(direction) == length(best$par) &&
all(is.finite(direction))) {
sum(pre_newton$gradient * direction)
} else {
NA_real_
}
valid_direction <- length(direction) == length(best$par) &&
all(is.finite(direction)) && any(direction != 0) &&
is.finite(descent_slope) && descent_slope < 0
if (!valid_direction) {
newton_metadata$status <- "skipped: Newton direction was not finite and descending"
append_optimizer_history(
stage = "Newton correction skipped",
attempt = 0L,
point = pre_newton,
hessian_reciprocal_condition = hessian_rcond,
transition = "stop_newton",
transition_reason = "newton_direction_not_finite_and_descending",
message = newton_metadata$status
)
} else {
upper_limited <- direction > 0 & is.finite(bounds_sens$upper)
lower_limited <- direction < 0 & is.finite(bounds_sens$lower)
feasible_limits <- c(
(bounds_sens$upper[upper_limited] - best$par[upper_limited]) /
direction[upper_limited],
(bounds_sens$lower[lower_limited] - best$par[lower_limited]) /
direction[lower_limited]
)
feasible_limits <- feasible_limits[
is.finite(feasible_limits) & feasible_limits >= 0
]
maximum_feasible_alpha <- if (length(feasible_limits)) {
min(feasible_limits)
} else {
Inf
}
initial_alpha <- newton_contract$initial_alpha
if (is.finite(maximum_feasible_alpha)) {
initial_alpha <- min(
initial_alpha,
newton_contract$fraction_to_boundary * maximum_feasible_alpha
)
}
if (!is.finite(initial_alpha) || initial_alpha <= 0) {
newton_metadata$status <- "skipped: Newton direction had no positive bounded step"
append_optimizer_history(
stage = "Newton correction skipped",
attempt = 0L,
point = pre_newton,
hessian_reciprocal_condition = hessian_rcond,
transition = "stop_newton",
transition_reason = "no_positive_bounded_newton_step",
message = newton_metadata$status
)
} else {
for (backtrack in 0:newton_contract$maximum_backtracks) {
alpha <- initial_alpha * newton_contract$backtrack_factor^backtrack
candidate <- best$par + alpha * direction
candidate_point <- sensitivity_point_diagnostics(
obj_sens,
candidate,
bounds_sens,
bound_tolerance = 0
)
armijo_pass <- isTRUE(candidate_point$valid) &&
candidate_point$objective <=
pre_newton$objective +
newton_contract$armijo_c1 * alpha * descent_slope
gradient_improvement <- isTRUE(candidate_point$finite_gradient) &&
candidate_point$max_abs_gradient < pre_newton$max_abs_gradient
gradient_requirement_pass <-
!isTRUE(newton_contract$require_gradient_improvement) ||
gradient_improvement
accepted_newton <- armijo_pass && gradient_requirement_pass
newton_progress <- sensitivity_optimizer_progress(
best$par,
pre_newton,
candidate,
candidate_point
)
append_optimizer_history(
stage = "bounded Newton line search",
attempt = backtrack + 1L,
point = candidate_point,
accepted_as_best = accepted_newton,
hessian_reciprocal_condition = hessian_rcond,
alpha = alpha,
armijo_pass = armijo_pass,
gradient_improvement = gradient_improvement,
progress = newton_progress,
transition = if (accepted_newton) {
"retain_newton_candidate_then_certify"
} else {
"continue_newton_backtracking"
},
transition_reason = if (accepted_newton) {
"armijo_and_gradient_improvement_passed"
} else if (!armijo_pass) {
"armijo_check_failed"
} else {
"gradient_improvement_check_failed"
},
message = if (accepted_newton) {
"accepted SPD Newton correction"
} else {
"rejected Newton candidate"
}
)
if (accepted_newton) {
best <- list(par = candidate, diagnostics = candidate_point)
newton_metadata <- list(
attempted = TRUE,
accepted = TRUE,
status = "accepted; bounded nlminb certification required",
hessian_reciprocal_condition = hessian_rcond,
alpha = alpha,
max_abs_step = max(abs(alpha * direction)),
pre_objective = pre_newton$objective,
post_objective = candidate_point$objective,
pre_max_abs_gradient = pre_newton$max_abs_gradient,
post_max_abs_gradient = candidate_point$max_abs_gradient,
armijo_c1 = newton_contract$armijo_c1,
descent_slope = descent_slope
)
break
}
}
if (!isTRUE(newton_metadata$accepted)) {
newton_metadata$status <-
"no bounded Newton candidate passed Armijo and gradient checks"
set_last_optimizer_transition(
"stop_newton",
"no_newton_candidate_passed_safeguards"
)
}
}
}
}
if (isTRUE(newton_metadata$accepted)) {
certification_contract <- sensitivity_optimizer_contract$certification
for (attempt in seq_len(certification_contract$max_passes)) {
certification_passes <- certification_passes + 1L
stage_result <- run_nlminb_stage(
stage = "bounded gradient nlminb certification",
attempt = attempt,
stage_control = certification_contract$control
)
certification_passed <- isTRUE(stage_result$accepted) &&
sensitivity_optimizer_certified(
stage_result$result,
stage_result$diagnostics,
gradient_tolerance
)
if (certification_passed) {
set_last_optimizer_transition(
"accept",
"strict_numeric_certification_passed"
)
opt_sens <- stage_result$result
final_stage <- "bounded gradient nlminb certification"
break
} else if (isTRUE(stage_result$accepted)) {
set_last_optimizer_transition(
"certification_failed",
"finite_nonworsening_uncertified"
)
} else {
set_last_optimizer_transition(
"certification_failed_with_rollback",
history$transition_reason[[nrow(history)]]
)
}
}
}
}
if (is.null(opt_sens)) {
stop(
label,
" did not meet the staged optimizer contract after ",
primary_passes, " primary pass(es), ", scaled_passes,
" scaled pass(es), and ", certification_passes,
" certification pass(es); best maximum gradient = ",
signif(best$diagnostics$max_abs_gradient, 6),
", best objective = ", signif(best$diagnostics$objective, 10), ".",
call. = FALSE
)
}
final_diagnostics <- sensitivity_point_diagnostics(
obj_sens,
opt_sens$par,
bounds_sens,
bound_tolerance = bound_tolerance
)
if (!sensitivity_optimizer_certified(
opt_sens,
final_diagnostics,
gradient_tolerance
)) {
stop(
label,
" failed strict final optimizer certification after the staged fit.",
call. = FALSE
)
}
if (any(opt_sens$par < bounds_sens$lower) ||
any(opt_sens$par > bounds_sens$upper)) {
stop(label, " has final parameters outside the optimizer bounds.",
call. = FALSE)
}
opt_sens$objective <- final_diagnostics$objective
obj_sens$par <- opt_sens$par
obj_sens$env$last.par.best <- opt_sens$par
final_nll <- final_diagnostics$objective
max_gradient <- final_diagnostics$max_abs_gradient
estimability <- tryCatch(
{
final_hessian <- obj_sens$he(opt_sens$par)
if (!is.matrix(final_hessian) ||
!identical(
dim(final_hessian),
c(length(opt_sens$par), length(opt_sens$par))
) || any(!is.finite(final_hessian))) {
stop("The final Hessian is non-finite or has invalid dimensions.")
}
check_estimability(obj = obj_sens, h = final_hessian)
},
error = function(e) e
)
if (inherits(estimability, "error")) {
stop(label, " estimability check failed: ", conditionMessage(estimability), call. = FALSE)
}
if (length(estimability$WhichBad) > 0L) {
stop(label, " estimability check found non-estimable parameters.", call. = FALSE)
}
biological_state_mle <- diagnose_sbt_states(
obj_sens,
data = data_sens,
mle_par = opt_sens$par,
cores = 1L
)
if (!esc31_biological_state_diagnostics_pass(biological_state_mle, 1L)) {
stop(label, " failed the raw harvest or biological-state gate.",
call. = FALSE)
}
harvest_wall_mle <- sensitivity_harvest_wall_report_diagnostics(
data_sens,
obj_sens$report(opt_sens$par)
)
if (!isTRUE(harvest_wall_mle$passes)) {
stop(
label,
" failed the approved preventive harvest-wall gate: ",
harvest_wall_mle$reason,
call. = FALSE
)
}
symmetric_final_hessian <- (final_hessian + t(final_hessian)) / 2
final_hessian_cholesky <- !inherits(
try(chol(symmetric_final_hessian), silent = TRUE),
"try-error"
)
final_hessian_rcond <- tryCatch(
as.numeric(rcond(symmetric_final_hessian)),
error = function(e) NA_real_
)
at_lower_bound <- is.finite(bounds_sens$lower) &
abs(opt_sens$par - bounds_sens$lower) <= bound_tolerance
at_upper_bound <- is.finite(bounds_sens$upper) &
abs(opt_sens$par - bounds_sens$upper) <= bound_tolerance
primary_history <- history[
history$stage == "primary bounded nlminb",
,
drop = FALSE
]
if (nrow(primary_history) != 1L) {
stop(
"The single-pass primary optimizer contract produced an invalid history.",
call. = FALSE
)
}
optimizer_summary <- list(
schema_version = 2L,
contract_signature = sensitivity_optimizer_contract_signature,
certified = TRUE,
final_stage = final_stage,
primary_passes = primary_passes,
primary_transition = primary_history$transition[[1L]],
primary_transition_reason = primary_history$transition_reason[[1L]],
scaled_passes = scaled_passes,
newton = newton_metadata,
certification_passes = certification_passes,
total_nlminb_calls = nlminb_calls,
final_convergence = as.integer(opt_sens$convergence),
final_message = opt_sens$message %||% final_stage,
initial_objective = initial_nll,
final_objective = final_nll,
final_max_abs_gradient = max_gradient,
final_gradient_l2 = final_diagnostics$gradient_l2,
final_hessian_cholesky = final_hessian_cholesky,
final_hessian_reciprocal_condition = final_hessian_rcond,
at_lower_bound = sum(at_lower_bound),
at_upper_bound = sum(at_upper_bound),
bound_violations = 0L
)
fit <- sbt_add_optimisation(
fit,
opt = opt_sens,
estimability = estimability,
diagnostics = utils::modifyList(
list(
optimizer = list(
method = sensitivity_optimizer_contract$implementation,
n_passes = nlminb_calls
),
initial_nll = initial_nll,
final_nll = final_nll,
max_gradient = max_gradient,
gradient_l2 = final_diagnostics$gradient_l2,
optimizer_contract = sensitivity_optimizer_contract,
optimizer_contract_signature = sensitivity_optimizer_contract_signature,
optimizer_history = history,
optimizer_summary = optimizer_summary,
biological_state_mle = biological_state_mle,
esc31_harvest_wall_mle = harvest_wall_mle
),
diagnostics
),
optimizer = sensitivity_optimizer_contract$implementation,
check = FALSE
)
check_mle(
fit,
gradient_tolerance = gradient_tolerance,
cores = 1L,
mode = "full"
)
}The first 12 rows below reproduce the sensitivity table in the OMMP16 meeting report exactly, preserving its row order, column headings, codes, wording, grid scope, and priorities. The four blue rows are assessment comparisons added after OMMP16 and were not part of the original table.
ommp16_sensitivity_specification <- tribble(
~`Test name`, ~Code, ~`Conditioning and projection notes`, ~Grid, ~Priority,
"No UAM", "noUAM", "Remove NCNM catches from conditioning and projections.", "Mid cell", "H",
"CPUE_Drop5", "Drop_5yrs", "Eliminate the last 5 years of CPUE Series", "Mid cell", "H",
"Omega75", "cpueom75", "Power function for biomass-CPUE relationship with power = 0.75", "Mid cell", "H",
"Upq2008", "Cpueupq", "Estimate CPUE 2008 change in q for LL1 fleet", "Mid cell", "H",
"LL1_sel", "LL1_sel", "Allow the terminal 3-years to be flexibly estimated to evaluate impact on year-class uncertainty and magnitude", "Mid cell", "M",
"Indo_sel", "Indo_sel", "Restore bi-modality in Indonesian selectivity, more flexibility", "Mid cell", "H",
"NoPOP&HSP", "NoPOPHSP", "Exclude both close-kin data (Parent-Offspring and Half-Sibling Pairs)", "Mid cell", "H",
"No HSP", "NoHSP", "Exclude half-sibling-pair close-kin data", "Mid cell", "H",
"GTI", "Troll", "Includes the grid-type trolling index as additional recruitment index. Increase CV of aerial survey to preclude aerial survey dominating the fit, given apparent conflicts in the data", "Mid cell", "H",
"Projections", "Alloc", "Alternative allocations", "Full grid", "H",
"OldGrid", "Og", "E.g., 108 (or other) dimension MLE with updated data and values for M0 and M10 (and maybe psi)", "Full grid", "H",
"CPUE_fleet", "U_Fleet", "Evaluate impact relative to assuming LL1 selectivity for CPUE predictions (as in previous approach)", "Mid cell", "H",
"CPUE_CV20", "constant_cpue_cv", "Restore the former time-invariant CPUE uncertainty with zero input log-SD and fixed additional log-scale sigma = 0.2", "Mid cell", "Additional",
"Length_m", "length_m", "Estimate M0, M4, M10, and M30 using the direct four-anchor natural-mortality parameterization", "Mid cell", "Additional",
"Estimate_m_slope", "estimate_m_slope", "Estimate the length-based natural-mortality exponent while retaining the selected base parameterization", "Mid cell", "Additional",
"No conventional tags", "no_tags", "Exclude the conventional-tag likelihood while retaining all other data and model settings", "Mid cell", "Additional"
)
ommp16_original_row_count <- 12L
additional_sensitivity_rows <- seq.int(
ommp16_original_row_count + 1L,
nrow(ommp16_sensitivity_specification)
)
ommp16_sensitivity_specification |>
kable(align = c("l", "l", "l", "c", "c")) |>
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover", "condensed"),
full_width = TRUE
) |>
kableExtra::column_spec(1L, bold = TRUE) |>
kableExtra::column_spec(5L, bold = TRUE) |>
kableExtra::row_spec(
additional_sensitivity_rows,
background = "#DDEBF7"
)| Test name | Code | Conditioning and projection notes | Grid | Priority |
|---|---|---|---|---|
| No UAM | noUAM | Remove NCNM catches from conditioning and projections. | Mid cell | H |
| CPUE_Drop5 | Drop_5yrs | Eliminate the last 5 years of CPUE Series | Mid cell | H |
| Omega75 | cpueom75 | Power function for biomass-CPUE relationship with power = 0.75 | Mid cell | H |
| Upq2008 | Cpueupq | Estimate CPUE 2008 change in q for LL1 fleet | Mid cell | H |
| LL1_sel | LL1_sel | Allow the terminal 3-years to be flexibly estimated to evaluate impact on year-class uncertainty and magnitude | Mid cell | M |
| Indo_sel | Indo_sel | Restore bi-modality in Indonesian selectivity, more flexibility | Mid cell | H |
| NoPOP&HSP | NoPOPHSP | Exclude both close-kin data (Parent-Offspring and Half-Sibling Pairs) | Mid cell | H |
| No HSP | NoHSP | Exclude half-sibling-pair close-kin data | Mid cell | H |
| GTI | Troll | Includes the grid-type trolling index as additional recruitment index. Increase CV of aerial survey to preclude aerial survey dominating the fit, given apparent conflicts in the data | Mid cell | H |
| Projections | Alloc | Alternative allocations | Full grid | H |
| OldGrid | Og | E.g., 108 (or other) dimension MLE with updated data and values for M0 and M10 (and maybe psi) | Full grid | H |
| CPUE_fleet | U_Fleet | Evaluate impact relative to assuming LL1 selectivity for CPUE predictions (as in previous approach) | Mid cell | H |
| CPUE_CV20 | constant_cpue_cv | Restore the former time-invariant CPUE uncertainty with zero input log-SD and fixed additional log-scale sigma = 0.2 | Mid cell | Additional |
| Length_m | length_m | Estimate M0, M4, M10, and M30 using the direct four-anchor natural-mortality parameterization | Mid cell | Additional |
| Estimate_m_slope | estimate_m_slope | Estimate the length-based natural-mortality exponent while retaining the selected base parameterization | Mid cell | Additional |
| No conventional tags | no_tags | Exclude the conventional-tag likelihood while retaining all other data and model settings | Mid cell | Additional |
OMMP16 noted that some sensitivities might be run over the full MCMC grid where time allowed. In the current workflow, Projections and OldGrid are handled on the grid and projection pages rather than as mid cell fits. The Troll implementation adds the trolling index but retains aerial tau at the base value of 0.59; it does not implement the requested increase because the base tau was already calibrated to improve aerial-survey SDNR. The recorded decision accepts this as the trolling-index-only sensitivity; any larger-tau case is future work requiring an explicit numerical specification. OMMP16 did not give numeric hyperparameters for Indo_sel; its working definition is therefore stated explicitly in that section below.
Each tab below gives the sensitivity definition and its immediate MLE TRO comparison with the base model. Run-specific diagnostics that examine other quantities are collected below the tabset.
Every call uses the same staged package lifecycle. The sensitivity-specific data and parameter map are constructed first, after which the executable path is:
fit <- sbt_fit(sensitivity_data)
fit <- sbt_add_parameters(fit, sensitivity_parameters)
fit <- sbt_add_map(fit, sensitivity_map)
fit <- sbt_add_priors(fit)
fit <- sbt_build_object(fit)
fit <- sbt_optimise(fit, n_passes = 1L, check = FALSE)The ordinary first pass is therefore the package sbt_optimise() workflow. Only a fit that remains uncertified after that pass enters the recorded scaled-gradient and safeguarded-Newton recovery policy; its final result is returned to the same portable lifecycle with sbt_add_optimisation().
This sensitivity removes all NCNM additions in catch_UA.csv before the catch array is constructed.
uam_columns <- setdiff(names(data_in$catch_UA), "Year")
uam_removed_total <- sum(unlist(data_in$catch_UA[uam_columns]), na.rm = TRUE)
no_uam_input <- data_in
no_uam_input$catch_UA <- no_uam_input$catch_UA |>
mutate(across(all_of(uam_columns), ~ .x * 0))
no_uam_data <- get_data(data_in = no_uam_input)
no_uam_details <- list(catch_UA_zeroed = TRUE, removed_total = uam_removed_total)
no_uam_expected <- sensitivity_metadata(
"noUAM",
no_uam_details,
data_sens = no_uam_data
)
no_uam_fit <- if (rerun_sensitivities) NULL else {
read_saved_sensitivity_fit(sensitivity_files[["no_uam"]], no_uam_expected)
}
if (is.null(no_uam_fit)) {
if (!run_sensitivity_fits) {
stop("No UAM fit is unavailable and run_sensitivity_fits is FALSE.", call. = FALSE)
}
no_uam_fit <- fit_sensitivity(
data_sens = no_uam_data,
sensitivity_id = "noUAM",
label = "No UAM",
metadata = no_uam_expected
)
sbt_fit_save(no_uam_fit, sensitivity_files[["no_uam"]], overwrite = TRUE)
}plot_sensitivity_tro(
fits = list(base_fit, no_uam_fit),
labels = c("Base", "No UAM")
)
This sensitivity removes the 2021-2025 abundance-index observations. The CPUE length compositions remain in the composition likelihood.
cpue_drop_years <- tail(sort(unique(as.integer(data_in$cpue$Year))), 5L)
drop_5yrs_input <- data_in
drop_5yrs_input$cpue <- drop_5yrs_input$cpue |>
filter(!.data$Year %in% cpue_drop_years)
drop_5yrs_data <- get_data(data_in = drop_5yrs_input)
drop_5yrs_details <- list(cpue_drop_years = cpue_drop_years)
drop_5yrs_expected <- sensitivity_metadata(
"Drop_5yrs",
drop_5yrs_details,
data_sens = drop_5yrs_data
)
drop_5yrs_fit <- if (rerun_sensitivities) NULL else {
read_saved_sensitivity_fit(sensitivity_files[["drop_5yrs"]], drop_5yrs_expected)
}
if (is.null(drop_5yrs_fit)) {
if (!run_sensitivity_fits) {
stop("Drop_5yrs fit is unavailable and run_sensitivity_fits is FALSE.", call. = FALSE)
}
drop_5yrs_fit <- fit_sensitivity(
data_sens = drop_5yrs_data,
sensitivity_id = "Drop_5yrs",
label = "Drop CPUE 2021-2025",
metadata = drop_5yrs_expected
)
sbt_fit_save(drop_5yrs_fit, sensitivity_files[["drop_5yrs"]], overwrite = TRUE)
}plot_sensitivity_tro(
fits = list(base_fit, drop_5yrs_fit),
labels = c("Base", "Drop CPUE 2021-2025")
)
This sensitivity fixes the biomass-CPUE power parameter at 0.75.
cpue_omega <- 0.75
cpue_omega_075_data <- get_data(data_in = data_in)
cpue_omega_075_details <- list(cpue_omega = cpue_omega)
cpue_omega_075_overrides <- list(par_log_cpue_omega = log(cpue_omega))
cpue_omega_075_expected <- sensitivity_metadata(
"cpueom75",
cpue_omega_075_details,
data_sens = cpue_omega_075_data,
parameter_overrides = cpue_omega_075_overrides
)
cpue_omega_075_fit <- if (rerun_sensitivities) NULL else {
read_saved_sensitivity_fit(
sensitivity_files[["cpue_omega_075"]],
cpue_omega_075_expected
)
}
if (is.null(cpue_omega_075_fit)) {
if (!run_sensitivity_fits) {
stop("cpueom75 fit is unavailable and run_sensitivity_fits is FALSE.", call. = FALSE)
}
cpue_omega_075_fit <- fit_sensitivity(
data_sens = cpue_omega_075_data,
sensitivity_id = "cpueom75",
label = "CPUE power 0.75",
parameter_overrides = cpue_omega_075_overrides,
metadata = cpue_omega_075_expected
)
sbt_fit_save(
cpue_omega_075_fit,
sensitivity_files[["cpue_omega_075"]],
overwrite = TRUE
)
}plot_sensitivity_tro(
fits = list(base_fit, cpue_omega_075_fit),
labels = c("Base", "CPUE power 0.75")
)
This sensitivity estimates separate catchability blocks for 1969-2007 and 2008 onward.
cpue_q_years <- 2008L
q2008_input <- data_in
q2008_input$cpue_q_yrs <- cpue_q_years
q2008_data <- get_data(data_in = q2008_input)
q2008_details <- list(cpue_q_years = cpue_q_years)
q2008_expected <- sensitivity_metadata(
"Cpueupq",
q2008_details,
data_sens = q2008_data
)
q2008_fit <- if (rerun_sensitivities) NULL else {
read_saved_sensitivity_fit(sensitivity_files[["q2008"]], q2008_expected)
}
if (is.null(q2008_fit)) {
if (!run_sensitivity_fits) {
stop("Cpueupq fit is unavailable and run_sensitivity_fits is FALSE.", call. = FALSE)
}
q2008_fit <- fit_sensitivity(
data_sens = q2008_data,
sensitivity_id = "Cpueupq",
label = "CPUE q split from 2008",
metadata = q2008_expected
)
sbt_fit_save(q2008_fit, sensitivity_files[["q2008"]], overwrite = TRUE)
}plot_sensitivity_tro(
fits = list(base_fit, q2008_fit),
labels = c("Base", "CPUE q split from 2008")
)
The base model already has a 2023 LL1 node. This sensitivity adds separate nodes for 2024 and 2025 so all three terminal years are estimated annually.
ll1_terminal_years <- seq.int(data_in$last_yr - 2L, data_in$last_yr)
ll1_terminal_input <- data_in
ll1_terminal_input$sel_LL1_yrs <- sort(unique(c(
ll1_terminal_input$sel_LL1_yrs,
ll1_terminal_years
)))
ll1_terminal_data <- get_data(data_in = ll1_terminal_input)
ll1_terminal_details <- list(ll1_terminal_years = ll1_terminal_years)
ll1_terminal_expected <- sensitivity_metadata(
"LL1_sel",
ll1_terminal_details,
data_sens = ll1_terminal_data
)
ll1_terminal_fit <- if (rerun_sensitivities) NULL else {
read_saved_sensitivity_fit(
sensitivity_files[["ll1_terminal_3yr"]],
ll1_terminal_expected
)
}
if (is.null(ll1_terminal_fit)) {
if (!run_sensitivity_fits) {
stop("LL1_sel fit is unavailable and run_sensitivity_fits is FALSE.", call. = FALSE)
}
ll1_terminal_fit <- fit_sensitivity(
data_sens = ll1_terminal_data,
sensitivity_id = "LL1_sel",
label = "Annual LL1 selectivity 2023-2025",
metadata = ll1_terminal_expected
)
sbt_fit_save(
ll1_terminal_fit,
sensitivity_files[["ll1_terminal_3yr"]],
overwrite = TRUE
)
}plot_sensitivity_tro(
fits = list(base_fit, ll1_terminal_fit),
labels = c("Base", "Annual LL1 selectivity 2023-2025")
)
OMMP16 defines this high-priority sensitivity as restoring bimodality in Indonesian selectivity by allowing more flexibility. The report does not give an exact numeric prior. For this review MLE, the Indonesian separable-AR1 hyperparameters are changed from the base values (rho_year = 0.98, rho_age = 0.98, sigma = 0.12) to the earlier exploratory relaxation (rho_year = 0.5, rho_age = 0.5, sigma = 0.5). All other model data, parameters, maps, and acceptance gates are unchanged. The final Indonesian selectivity block still begins in 2021 and therefore applies to 2021–2025, as requested in the OMMP16 report.
Following review of the MLE selectivity shape, this 0.5/0.5/0.5 configuration is the accepted working definition for the Indo_sel sensitivity and is now a production MCMC target. It is appended after the established sensitivity registry so none of the existing model seeds change. It remains additional to, rather than part of, the base-only grid prerequisite.
indo_sel_hyperparameters <- c(
rho_year = 0.5,
rho_age = 0.5,
sigma = 0.5
)
indo_sel_data <- get_data(data_in = data_in)
indo_sel_overrides <- list(
par_sel_rho_y = base_fit$parameters$par_sel_rho_y,
par_sel_rho_a = base_fit$parameters$par_sel_rho_a,
par_log_sel_sigma = base_fit$parameters$par_log_sel_sigma
)
indo_sel_overrides$par_sel_rho_y[[5L]] <- sel_rho_to_par(
indo_sel_hyperparameters[["rho_year"]]
)
indo_sel_overrides$par_sel_rho_a[[5L]] <- sel_rho_to_par(
indo_sel_hyperparameters[["rho_age"]]
)
indo_sel_overrides$par_log_sel_sigma[[5L]] <- log(
indo_sel_hyperparameters[["sigma"]]
)
indo_sel_details <- list(
role = "production sensitivity",
decision_source = "OMMP16 sensitivity table",
target = "restore bimodality with more flexible Indonesian selectivity",
selected_hyperparameters = as.list(indo_sel_hyperparameters),
selection_basis = "accepted review of earlier ESC31 exploratory relaxation",
final_selectivity_block = 2021:2025,
production_mcmc_registry = TRUE,
grid_prerequisite = FALSE
)
indo_sel_expected <- sensitivity_metadata(
"Indo_sel",
indo_sel_details,
data_sens = indo_sel_data,
parameter_overrides = indo_sel_overrides
)
indo_sel_fit <- if (rerun_sensitivities) NULL else {
read_saved_sensitivity_fit(
sensitivity_mle_files[["indo_sel"]],
indo_sel_expected
)
}
if (is.null(indo_sel_fit)) {
if (!run_sensitivity_fits) {
stop(
"The Indo_sel MLE is unavailable and ",
"run_sensitivity_fits is FALSE.",
call. = FALSE
)
}
indo_sel_fit <- fit_sensitivity(
data_sens = indo_sel_data,
sensitivity_id = "Indo_sel",
label = "More flexible Indonesian selectivity",
parameter_overrides = indo_sel_overrides,
metadata = indo_sel_expected
)
sbt_fit_save(
indo_sel_fit,
sensitivity_mle_files[["indo_sel"]],
overwrite = TRUE
)
}
stopifnot(
sensitivity_mle_passes(indo_sel_fit),
identical(
max(as.integer(colnames(indo_sel_fit$data$sel_change_year_fy))[
indo_sel_fit$data$sel_change_year_fy[5L, ] != 0L
]),
2021L
)
)indo_sel_state <- indo_sel_fit$fit$diagnostics$biological_state_mle$summary[1L, ]
indo_sel_state_draw <-
indo_sel_fit$fit$diagnostics$biological_state_mle$per_draw[1L, ]
indo_sel_wall <- sensitivity_harvest_wall_fit_diagnostics(indo_sel_fit)
indo_sel_diagnostics <- tibble(
`rho year` = indo_sel_hyperparameters[["rho_year"]],
`rho age` = indo_sel_hyperparameters[["rho_age"]],
Sigma = indo_sel_hyperparameters[["sigma"]],
Objective = as.numeric(indo_sel_fit$fit$opt$objective),
Convergence = as.integer(indo_sel_fit$fit$opt$convergence),
`Max gradient` = as.numeric(indo_sel_fit$fit$diagnostics$max_gradient),
Estimable = identical(indo_sel_fit$fit$estimability$status, "estimable"),
`Max raw harvest` = as.numeric(indo_sel_state$max_raw_harvest_rate),
`Min abundance` = as.numeric(indo_sel_state$min_number),
`Max catch error` = as.numeric(
indo_sel_state_draw$max_catch_relative_error
),
`Max continuation` = as.numeric(indo_sel_state$max_harvest_penalty),
`Wall objective` = as.numeric(indo_sel_wall$total_wall_objective),
`MLE gate` = if (sensitivity_mle_passes(indo_sel_fit)) "Pass" else "Fail"
)
indo_sel_diagnostics |>
mutate(
Objective = format_decimal(.data$Objective, 6L),
`Max gradient` = format_decimal(.data$`Max gradient`, 9L),
`Max raw harvest` = format_decimal(.data$`Max raw harvest`, 6L),
`Min abundance` = format_decimal(.data$`Min abundance`, 3L),
`Max catch error` = format_scientific(.data$`Max catch error`, 3L),
`Max continuation` = format_scientific(.data$`Max continuation`, 3L),
`Wall objective` = format_scientific(.data$`Wall objective`, 3L)
) |>
replace_missing() |>
kable(align = c(rep("r", 6L), "l", rep("r", 5L), "l"))| rho year | rho age | Sigma | Objective | Convergence | Max gradient | Estimable | Max raw harvest | Min abundance | Max catch error | Max continuation | Wall objective | MLE gate |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.5 | 0.5 | 0.5 | 8,129.683961 | 0 | 0.000000000 | TRUE | 0.579650 | 2,897.195 | 3.019e-16 | 0.000e+00 | 1.700e-24 | Pass |
plot_sensitivity_tro(
fits = list(base_fit, indo_sel_fit),
labels = c("Base", "More flexible Indo_sel")
)
Earlier no-wall NoPOPHSP trials were boundary-dominated. That conclusion is superseded by the approved global preventive wall, which applies to the base model and every sensitivity rather than being introduced only for NoPOPHSP. An isolated exact-switch A10 pilot passed all MLE and biological-state gates, but its artifact is pilot-only: the production fit below must be rebuilt or loaded under the current sensitivity identity and pass the same strict gates as every other model.
A final isolated recovery test retained 3,000 draws in each of four chains after 1,000 warmup iterations, used adapt_delta = 0.9995, seed 73107, and the same exact-AD-Hessian dense metric. It eliminated divergent transitions and all 12,000 retained states passed, but maximum R-hat was 1.033873, minimum bulk ESS was 82.440, and minimum tail ESS was 186.027, all for par_log_m10. That candidate therefore failed the unchanged standard gates and was not promoted. The canonical result below remains the previously reviewed 12,000-draw run accepted under its explicit two-divergence decision.
no_pop_hsp_input <- data_in
no_pop_hsp_input$pop_switch <- 0L
no_pop_hsp_input$hsp_switch <- 0L
no_pop_hsp_data <- get_data(data_in = no_pop_hsp_input)
no_pop_hsp_details <- list(pop_switch = 0L, hsp_switch = 0L)
no_pop_hsp_expected <- sensitivity_metadata(
"NoPOPHSP",
no_pop_hsp_details,
data_sens = no_pop_hsp_data
)
no_pop_hsp_fit <- if (rerun_sensitivities) NULL else {
read_saved_sensitivity_fit(
sensitivity_files[["no_pop_hsp"]],
no_pop_hsp_expected
)
}
if (is.null(no_pop_hsp_fit)) {
if (!run_sensitivity_fits) {
stop(
"NoPOPHSP fit is unavailable and run_sensitivity_fits is FALSE.",
call. = FALSE
)
}
no_pop_hsp_fit <- fit_sensitivity(
data_sens = no_pop_hsp_data,
sensitivity_id = "NoPOPHSP",
label = "No POP or HSP",
metadata = no_pop_hsp_expected
)
sbt_fit_save(
no_pop_hsp_fit,
sensitivity_files[["no_pop_hsp"]],
overwrite = TRUE
)
}
stopifnot(
identical(as.integer(no_pop_hsp_fit$data$pop_switch), 0L),
identical(as.integer(no_pop_hsp_fit$data$hsp_switch), 0L),
sensitivity_mle_passes(no_pop_hsp_fit)
)
no_pop_hsp_historical_evidence <- tibble::tribble(
~case, ~status, ~objective, ~max_raw_harvest, ~max_gradient,
~wall_objective, ~interpretation,
"Best unguarded homotopy solution", "Superseded pre-wall", 6962.885469443,
0.899999877, 0.902653748, 0,
"False convergence, high gradient, and non-estimability at the raw-harvest boundary.",
"Exact-switch global A10 certification", "Pilot-only pass", 6962.88550863527,
0.7903891490, 0.00000000023587, 0.0000026475,
"Independent isolated pilot; it demonstrates viability but is never reused as a production result."
)no_pop_hsp_historical_evidence |>
mutate(
objective = format_decimal(.data$objective, 6L),
max_raw_harvest = format_decimal(.data$max_raw_harvest, 9L),
max_gradient = format_decimal(.data$max_gradient, 9L),
wall_objective = format_decimal(.data$wall_objective, 9L)
) |>
replace_missing() |>
kable(
col.names = c(
"Case", "Status", "Objective", "Max raw harvest", "Max gradient",
"Wall objective", "Interpretation"
),
align = c("l", "l", rep("r", 4L), "l")
)| Case | Status | Objective | Max raw harvest | Max gradient | Wall objective | Interpretation |
|---|---|---|---|---|---|---|
| Best unguarded homotopy solution | Superseded pre-wall | 6,962.885469 | 0.899999877 | 0.902653748 | 0.000000000 | False convergence, high gradient, and non-estimability at the raw-harvest boundary. |
| Exact-switch global A10 certification | Pilot-only pass | 6,962.885509 | 0.790389149 | 0.000000000 | 0.000002648 | Independent isolated pilot; it demonstrates viability but is never reused as a production result. |
plot_sensitivity_tro(
fits = list(base_fit, no_pop_hsp_fit),
labels = c("Base", "No POP or HSP")
)
This sensitivity removes the HSP likelihood while retaining POPs and all other base-model inputs and settings.
no_hsp_input <- data_in
no_hsp_input$hsp_switch <- 0L
no_hsp_data <- get_data(data_in = no_hsp_input)
no_hsp_details <- list(pop_switch = 1L, hsp_switch = 0L)
no_hsp_expected <- sensitivity_metadata(
"NoHSP",
no_hsp_details,
data_sens = no_hsp_data
)
no_hsp_fit <- if (rerun_sensitivities) NULL else {
read_saved_sensitivity_fit(sensitivity_files[["no_hsp"]], no_hsp_expected)
}
if (is.null(no_hsp_fit)) {
if (!run_sensitivity_fits) {
stop("NoHSP fit is unavailable and run_sensitivity_fits is FALSE.", call. = FALSE)
}
no_hsp_fit <- fit_sensitivity(
data_sens = no_hsp_data,
sensitivity_id = "NoHSP",
label = "No HSP",
metadata = no_hsp_expected
)
sbt_fit_save(no_hsp_fit, sensitivity_files[["no_hsp"]], overwrite = TRUE)
}plot_sensitivity_tro(
fits = list(base_fit, no_hsp_fit),
labels = c("Base", "No HSP")
)
The conventional-tag tables are retained because they define the release and recapture array dimensions used by the model, but tag_switch is set to zero. This removes the complete conventional-tag likelihood while leaving all other data, parameters, maps, and acceptance gates unchanged.
no_tags_input <- data_in
no_tags_input$tag_switch <- 0L
no_tags_data <- get_data(data_in = no_tags_input)
no_tags_details <- list(
tag_switch = 0L,
excluded_likelihood = "conventional tags",
implementation = "sbt_tag_switch_zero_v1",
production_mcmc_registry = TRUE,
mcmc_status = "deferred after accepted MLE"
)
no_tags_expected <- sensitivity_metadata(
"No_tags",
no_tags_details,
data_sens = no_tags_data
)
no_tags_fit <- if (rerun_sensitivities) NULL else {
read_saved_sensitivity_fit(
sensitivity_files[["no_tags"]],
no_tags_expected
)
}
if (is.null(no_tags_fit)) {
if (!run_sensitivity_fits) {
stop(
"The no-conventional-tags MLE is unavailable and ",
"run_sensitivity_fits is FALSE.",
call. = FALSE
)
}
no_tags_fit <- fit_sensitivity(
data_sens = no_tags_data,
sensitivity_id = "No_tags",
label = "No conventional tags",
metadata = no_tags_expected
)
sbt_fit_save(
no_tags_fit,
sensitivity_files[["no_tags"]],
overwrite = TRUE
)
}
stopifnot(
identical(as.integer(no_tags_fit$data$tag_switch), 0L),
sensitivity_mle_passes(no_tags_fit)
)plot_sensitivity_tro(
fits = list(base_fit, no_tags_fit),
labels = c("Base", "No conventional tags")
)
OMMP16 requests a compound sensitivity that activates the grid-type trolling index and increases the aerial-survey CV, but the report does not give a numeric value. The executable case below fixes aerial tau at 0.59, exactly the base level, and therefore isolates activation of the trolling index rather than implementing an unspecified additional increase. The base tau was already calibrated to 0.59 to bring aerial-survey SDNR close to one. Decision esc31_2026_troll_base_aerial_tau_059_v1 accepts this as the defined trolling-index-only sensitivity. Any larger-tau test would require a new numerical specification and a separately identified sensitivity.
troll_decision <- review_method_decisions$troll
stopifnot(
identical(troll_decision$status, "accepted"),
identical(
troll_decision$decision_id,
"esc31_2026_troll_base_aerial_tau_059_v1"
),
identical(troll_decision$interpretation, "trolling_index_only"),
isTRUE(troll_decision$future_larger_tau_requires_new_sensitivity)
)
troll_aerial_tau <- troll_decision$aerial_tau
troll_input <- data_in
troll_input$troll_switch <- 1L
troll_input$aerial_tau <- troll_aerial_tau
troll_data <- get_data(data_in = troll_input)
troll_details <- list(troll_switch = 1L, aerial_tau = troll_aerial_tau)
troll_overrides <- list(par_log_aerial_tau = log(troll_aerial_tau))
troll_expected <- sensitivity_metadata(
"Troll",
troll_details,
data_sens = troll_data,
parameter_overrides = troll_overrides
)
troll_fit <- if (rerun_sensitivities) NULL else {
read_saved_sensitivity_fit(sensitivity_files[["troll"]], troll_expected)
}
if (is.null(troll_fit)) {
if (!run_sensitivity_fits) {
stop("Troll fit is unavailable and run_sensitivity_fits is FALSE.", call. = FALSE)
}
troll_fit <- fit_sensitivity(
data_sens = troll_data,
sensitivity_id = "Troll",
label = "GTI with base aerial tau 0.59",
parameter_overrides = troll_overrides,
metadata = troll_expected
)
sbt_fit_save(troll_fit, sensitivity_files[["troll"]], overwrite = TRUE)
}plot_sensitivity_tro(
fits = list(base_fit, troll_fit),
labels = c("Base", "GTI with base aerial tau 0.59")
)
This sensitivity changes only the selectivity used for the CPUE abundance index. CPUE length compositions continue to use the separately estimated fishery-7 selectivity.
u_fleet_cpue_sel_fishery <- 1L
u_fleet_cpue_lf_sel_fishery <- 7L
u_fleet_input <- data_in
u_fleet_input$cpue_sel_fishery <- u_fleet_cpue_sel_fishery
u_fleet_data <- get_data(data_in = u_fleet_input)
u_fleet_details <- list(
cpue_sel_fishery = u_fleet_cpue_sel_fishery,
cpue_lf_sel_fishery = u_fleet_cpue_lf_sel_fishery
)
u_fleet_expected <- sensitivity_metadata(
"U_Fleet",
u_fleet_details,
data_sens = u_fleet_data
)
u_fleet_fit <- if (rerun_sensitivities) NULL else {
read_saved_sensitivity_fit(sensitivity_files[["cpue_ll1_sel"]], u_fleet_expected)
}
if (is.null(u_fleet_fit)) {
if (!run_sensitivity_fits) {
stop("U_Fleet fit is unavailable and run_sensitivity_fits is FALSE.", call. = FALSE)
}
u_fleet_fit <- fit_sensitivity(
data_sens = u_fleet_data,
sensitivity_id = "U_Fleet",
label = "CPUE index with LL1 selectivity",
metadata = u_fleet_expected
)
sbt_fit_save(u_fleet_fit, sensitivity_files[["cpue_ll1_sel"]], overwrite = TRUE)
}plot_sensitivity_tro(
fits = list(base_fit, u_fleet_fit),
labels = c("Base", "CPUE index with LL1 selectivity")
)
The revised base combines a 20% baseline CV with the raw year-specific GAM22 estimation CV by adding their log variances. This sensitivity restores the former time-invariant CPUE uncertainty exactly: the raw CPUE table has zero input log-SD and the fixed additional log-scale sigma is 0.2. This is an ordinary-scale CV of approximately 0.202; setting an ordinary CV column to 0.2 would silently change the historical likelihood by about 0.98%.
constant_cpue_sigma <- 0.2
constant_cpue_cv_input <- data_in
constant_cpue_cv_input$cpue$CV <- NULL
constant_cpue_cv_data <- get_data(data_in = constant_cpue_cv_input)
constant_cpue_cv_details <- list(
cpue_input_uncertainty = "raw CPUE with zero input log-SD",
fixed_cpue_sigma = constant_cpue_sigma,
implementation = "legacy_constant_cpue_sigma_v1"
)
constant_cpue_cv_overrides <- list(
par_log_cpue_sigma = log(constant_cpue_sigma)
)
constant_cpue_cv_expected <- sensitivity_metadata(
"Constant_CV_0.2",
constant_cpue_cv_details,
data_sens = constant_cpue_cv_data,
parameter_overrides = constant_cpue_cv_overrides
)
constant_cpue_cv_fit <- if (rerun_sensitivities) NULL else {
read_saved_sensitivity_fit(
sensitivity_files[["constant_cpue_cv"]],
constant_cpue_cv_expected
)
}
if (is.null(constant_cpue_cv_fit)) {
if (!run_sensitivity_fits) {
stop("Constant CPUE CV fit is unavailable and run_sensitivity_fits is FALSE.", call. = FALSE)
}
constant_cpue_cv_fit <- fit_sensitivity(
data_sens = constant_cpue_cv_data,
sensitivity_id = "Constant_CV_0.2",
label = "Constant CPUE CV 0.2",
parameter_overrides = constant_cpue_cv_overrides,
metadata = constant_cpue_cv_expected
)
sbt_fit_save(
constant_cpue_cv_fit,
sensitivity_files[["constant_cpue_cv"]],
overwrite = TRUE
)
}plot_sensitivity_tro(
fits = list(base_fit, constant_cpue_cv_fit),
labels = c("Base", "Constant CPUE CV 0.2")
)
The selected base estimates \(M_{10}\) and \(M_{30}\) under the length-based curve with \(m_c=-1\). This comparison restores the package’s direct log-anchor parameterization and estimates \(M_0\), \(M_4\), \(M_{10}\), and \(M_{30}\). The established length_m registry key is retained so its final registry position, file name, and seed remain stable; the scientific identity records the actual M_switch = 1 data and parameter map.
direct_m_input <- data_in
direct_m_input$M_switch <- 1L
direct_m_data <- get_data(data_in = direct_m_input)
base_mortality_report <- sbt_fit_report(base_fit)
direct_m_overrides <- list(
par_log_m0 = log(as.numeric(base_mortality_report$par_m0)),
par_log_m4 = log(as.numeric(base_mortality_report$par_m4))
)
direct_m_details <- list(
registry_key = "length_m",
comparison = "direct four-anchor mortality against length-based base",
M_switch = 1L,
estimated_anchors = c(0L, 4L, 10L, 30L),
base_length_exponent = -1,
implementation = "sbt_direct_log_mortality_v1"
)
length_m_expected <- sensitivity_metadata(
"Direct_M_anchors",
direct_m_details,
data_sens = direct_m_data,
parameter_overrides = direct_m_overrides
)
length_m_fit <- if (rerun_sensitivities) NULL else {
read_saved_sensitivity_fit(
sensitivity_files[["length_m"]],
length_m_expected
)
}
if (is.null(length_m_fit)) {
if (!run_sensitivity_fits) {
stop(
"Direct-anchor mortality fit is unavailable and ",
"run_sensitivity_fits is FALSE.",
call. = FALSE
)
}
length_m_fit <- fit_sensitivity(
data_sens = direct_m_data,
sensitivity_id = "Direct_M_anchors",
label = "Direct four-anchor M",
parameter_overrides = direct_m_overrides,
metadata = length_m_expected
)
sbt_fit_save(
length_m_fit,
sensitivity_files[["length_m"]],
overwrite = TRUE
)
}
stopifnot(
identical(as.integer(length_m_fit$data$M_switch), 1L),
sensitivity_mle_passes(length_m_fit)
)plot_sensitivity_tro(
fits = list(base_fit, length_m_fit),
labels = c("Length-based base", "Direct four-anchor M")
)
The selected base deliberately fixes the length-based natural-mortality exponent \(m_c\) at -1 and estimates \(M_{10}\) and \(M_{30}\). This separate sensitivity retains the same length-based curve and all other base settings but estimates \(m_c\). The package prior is \(m_c \sim N(-1, 0.3^2)\) and its existing bounds are -2 to -0.25. The base specification is unchanged.
estimate_m_slope_data <- get_data(data_in = data_in)
estimate_m_slope_details <- list(
comparison = "estimate length-based natural-mortality exponent",
M_switch = 2L,
estimated_parameter = "par_mc",
base_exponent = -1,
prior = list(distribution = "normal", mean = -1, sd = 0.3),
bounds = c(lower = -2, upper = -0.25),
implementation = "sbt_length_m_estimated_exponent_v1"
)
estimate_m_slope_expected <- sensitivity_metadata(
"Estimate_M_slope",
estimate_m_slope_details,
data_sens = estimate_m_slope_data,
estimated_parameters = "par_mc"
)
estimate_m_slope_fit <- if (rerun_sensitivities) NULL else {
read_saved_sensitivity_fit(
sensitivity_files[["estimate_m_slope"]],
estimate_m_slope_expected
)
}
if (is.null(estimate_m_slope_fit)) {
if (!run_sensitivity_fits) {
stop(
"The estimated-M-slope fit is unavailable and ",
"run_sensitivity_fits is FALSE.",
call. = FALSE
)
}
estimate_m_slope_fit <- fit_sensitivity(
data_sens = estimate_m_slope_data,
sensitivity_id = "Estimate_M_slope",
label = "Estimated length-based M slope",
estimated_parameters = "par_mc",
metadata = estimate_m_slope_expected
)
sbt_fit_save(
estimate_m_slope_fit,
sensitivity_files[["estimate_m_slope"]],
overwrite = TRUE
)
}
estimate_m_slope_mc <- as.numeric(
sbt_fit_report(estimate_m_slope_fit)$par_mc
)
stopifnot(
identical(as.integer(estimate_m_slope_fit$data$M_switch), 2L),
is.finite(estimate_m_slope_mc),
estimate_m_slope_mc >= -2,
estimate_m_slope_mc <= -0.25,
sensitivity_mle_passes(estimate_m_slope_fit)
)The MLE estimates \(m_c\) as -0.969.
plot_sensitivity_tro(
fits = list(base_fit, estimate_m_slope_fit),
labels = c("Base: slope fixed at -1", "Estimated M-at-age slope")
)
The following plots examine sensitivity-specific mechanisms rather than repeat the base-versus-sensitivity TRO comparison shown in the tabs above.
p_sensitivity_q2008_posteriors
The Indonesian-selectivity comparison shows representative blocks from the base model, the more-flexible sensitivity, and the refitted V1 model. The terminal panel is the block applied throughout 2021–2025.
p_sensitivity_indo_sel
p_sensitivity_no_tags_m_at_age
p_sensitivity_no_tags_m_parameters
p_sensitivity_length_m_at_age
p_sensitivity_estimate_m_slope_at_age
accepted_sensitivity_keys <- setdiff(
names(sensitivity_files),
rejected_sensitivity_mle_keys
)
sensitivity_fits <- list(
no_uam = no_uam_fit,
drop_5yrs = drop_5yrs_fit,
cpue_omega_075 = cpue_omega_075_fit,
q2008 = q2008_fit,
ll1_terminal_3yr = ll1_terminal_fit,
no_pop_hsp = no_pop_hsp_fit,
no_hsp = no_hsp_fit,
troll = troll_fit,
cpue_ll1_sel = u_fleet_fit,
constant_cpue_cv = constant_cpue_cv_fit,
length_m = length_m_fit,
indo_sel = indo_sel_fit,
estimate_m_slope = estimate_m_slope_fit,
no_tags = no_tags_fit
)
sensitivity_labels <- c(
no_uam = "No UAM",
drop_5yrs = "Drop CPUE 2021-2025",
cpue_omega_075 = "CPUE power 0.75",
q2008 = "CPUE q split from 2008",
ll1_terminal_3yr = "Annual LL1 selectivity 2023-2025",
no_pop_hsp = "No POP or HSP",
no_hsp = "No HSP",
troll = "GTI with base aerial tau 0.59",
cpue_ll1_sel = "CPUE index with LL1 selectivity",
constant_cpue_cv = "Constant CPUE CV 0.2",
length_m = "Direct four-anchor M",
indo_sel = "More flexible Indonesian selectivity",
estimate_m_slope = "Estimated M-at-age slope",
no_tags = "No conventional tags"
)
stopifnot(
identical(names(sensitivity_fits), accepted_sensitivity_keys),
sensitivity_mle_passes(no_pop_hsp_fit)
)
cpue_sensitivity_keys <- c(
"drop_5yrs", "cpue_omega_075", "q2008", "cpue_ll1_sel",
"constant_cpue_cv"
)
cpue_sensitivity_summary <- bind_rows(
summarize_fit("Base", base_fit),
imap_dfr(
sensitivity_fits[cpue_sensitivity_keys],
~ summarize_fit(sensitivity_labels[[.y]], .x)
)
)
mle_review_sources <- c(list(base = base_fit), sensitivity_fits)
mle_review_labels <- c(base = "Base", sensitivity_labels)
mle_review_summary <- imap_dfr(
mle_review_sources,
function(fit, key) {
state <- fit$fit$diagnostics$biological_state_mle$summary[1L, ]
state_draw <- fit$fit$diagnostics$biological_state_mle$per_draw[1L, ]
wall <- sensitivity_harvest_wall_fit_diagnostics(fit)
tibble(
Model = unname(mle_review_labels[[key]]),
Objective = as.numeric(fit$fit$opt$objective),
`Active parameters` = length(fit$fit$opt$par),
Convergence = as.integer(fit$fit$opt$convergence),
`Max gradient` =
as.numeric(fit$fit$diagnostics$max_gradient %||% NA_real_),
Estimable = identical(fit$fit$estimability$status, "estimable"),
`Max raw harvest` = as.numeric(state$max_raw_harvest_rate),
`Min abundance` = as.numeric(state$min_number),
`Max catch error` = as.numeric(state_draw$max_catch_relative_error),
`Max continuation` = as.numeric(state$max_harvest_penalty),
`Wall objective` = as.numeric(wall$total_wall_objective),
`MLE gate` = if (sensitivity_mle_passes(fit)) "Pass" else "Fail"
)
}
)
stopifnot(all(mle_review_summary$`MLE gate` == "Pass"))mle_fit_sources <- if (exists("mcmc_worker_sources", inherits = FALSE)) {
mcmc_worker_sources
} else {
c(list(base = base_fit), sensitivity_fits)
}
mle_fit_identity_signatures <- if (
exists("mcmc_worker_identity_signatures", inherits = FALSE)
) {
mcmc_worker_identity_signatures
} else {
c(
base = base_fit$provenance$metadata$run_identity$signature,
no_uam = no_uam_expected$workflow_run_identity$signature,
drop_5yrs = drop_5yrs_expected$workflow_run_identity$signature,
cpue_omega_075 = cpue_omega_075_expected$workflow_run_identity$signature,
q2008 = q2008_expected$workflow_run_identity$signature,
ll1_terminal_3yr = ll1_terminal_expected$workflow_run_identity$signature,
no_pop_hsp = no_pop_hsp_expected$workflow_run_identity$signature,
no_hsp = no_hsp_expected$workflow_run_identity$signature,
troll = troll_expected$workflow_run_identity$signature,
cpue_ll1_sel = u_fleet_expected$workflow_run_identity$signature,
constant_cpue_cv =
constant_cpue_cv_expected$workflow_run_identity$signature,
length_m = length_m_expected$workflow_run_identity$signature,
indo_sel = indo_sel_expected$workflow_run_identity$signature,
estimate_m_slope =
estimate_m_slope_expected$workflow_run_identity$signature,
no_tags = no_tags_expected$workflow_run_identity$signature
)
}
active_sensitivity_mcmc_files <- sensitivity_mcmc_files[
names(mle_fit_sources)
]
stopifnot(
identical(
sensitivity_mcmc_registry_keys,
c("base", names(sensitivity_files))
),
identical(
valid_sensitivity_mcmc_keys,
sensitivity_mcmc_registry_keys
),
length(mle_fit_sources) >= 1L,
all(names(mle_fit_sources) %in% sensitivity_mcmc_registry_keys),
identical(names(mle_fit_sources), names(active_sensitivity_mcmc_files)),
identical(names(mle_fit_sources), names(mle_fit_identity_signatures)),
all(vapply(mle_fit_sources, inherits, logical(1), what = "sbt_fit"))
)
fits <- Map(
function(key, source_fit, file, fit_identity_signature) {
ensure_sensitivity_mcmc_fit(
key = key,
source_fit = source_fit,
file = file,
fit_identity_signature = fit_identity_signature
)
},
key = names(mle_fit_sources),
source_fit = unname(mle_fit_sources),
file = unname(active_sensitivity_mcmc_files),
fit_identity_signature = unname(mle_fit_identity_signatures)
)
sensitivity_mcmc_fits <- setNames(fits, names(mle_fit_sources))sensitivity_mcmc_fits <- Map(
function(key, fit, file) {
if (is.null(fit)) return(NULL)
acceptance_record <- if (
identical(key, no_pop_hsp_mcmc_acceptance$model_key) &&
sensitivity_mcmc_divergence_exception_passes(
fit$fit$diagnostics$mcmc
)
) {
utils::modifyList(
no_pop_hsp_mcmc_acceptance,
list(
mcmc_run_signature =
fit$provenance$metadata$mcmc_run_identity$signature,
posterior_payload_checksum =
as.character(
fit$fit$diagnostics$mcmc$posterior_payload_checksum
),
state_checksum =
as.character(fit$fit$diagnostics$mcmc$state_checksum)
)
)
} else {
NULL
}
fit_changed <- FALSE
if (!is.null(acceptance_record)) {
if (!identical(
fit$fit$diagnostics$mcmc_acceptance,
acceptance_record
)) {
fit$fit$diagnostics$mcmc_acceptance <- acceptance_record
fit_changed <- TRUE
}
if (!identical(
fit$provenance$metadata$mcmc_acceptance,
acceptance_record
)) {
fit$provenance$metadata$mcmc_acceptance <- acceptance_record
fit_changed <- TRUE
}
}
summaries_before <- fit$posterior_summaries
plot_summary_available <- isTRUE(tryCatch(
{
get_posterior_summary(fit, name = "plots")
TRUE
},
error = function(error) FALSE
))
if (!plot_summary_available) {
fit <- add_posterior_summaries(fit)
}
if (fit_changed ||
!identical(summaries_before, fit$posterior_summaries)) {
sbt_fit_save(fit, file, overwrite = TRUE)
}
fit
},
key = names(sensitivity_mcmc_fits),
fit = sensitivity_mcmc_fits,
file = unname(sensitivity_mcmc_files)
)
sensitivity_mcmc_fits <- setNames(
sensitivity_mcmc_fits,
names(mle_fit_sources)
)source(file.path(esc_dir, "esc31_grid_contract.R"), local = TRUE)
grid_prerequisite_file <- esc31_grid_prerequisite_acceptance_file(esc_dir)
grid_prerequisite_keys <-
esc31_grid_specification$launch_prerequisites$required_posterior_keys
grid_prerequisite_fits <- sensitivity_mcmc_fits[grid_prerequisite_keys]
grid_prerequisite_files <- sensitivity_mcmc_files[grid_prerequisite_keys]
grid_prerequisites_available <-
identical(names(grid_prerequisite_fits), grid_prerequisite_keys) &&
identical(names(grid_prerequisite_files), grid_prerequisite_keys) &&
all(vapply(grid_prerequisite_fits, function(fit) {
inherits(fit, "sbt_fit") && !is.null(fit$mcmc) &&
sensitivity_mcmc_passes(
fit$fit$diagnostics$mcmc
)
}, logical(1)))
refresh_grid_prerequisite_acceptance <- isTRUE(get0(
"refresh_grid_prerequisite_acceptance",
inherits = FALSE,
ifnotfound = FALSE
))
if ((run_sensitivity_mcmc || refresh_grid_prerequisite_acceptance) &&
grid_prerequisites_available) {
grid_prerequisite_acceptance <-
esc31_build_grid_prerequisite_acceptance(
esc_dir = esc_dir,
fits = grid_prerequisite_fits,
files = grid_prerequisite_files,
base_run_signature =
base_fit$provenance$metadata$run_identity$signature
)
esc31_atomic_save_grid_prerequisite_acceptance(
grid_prerequisite_acceptance,
grid_prerequisite_file
)
} else if (run_sensitivity_mcmc && file.exists(grid_prerequisite_file)) {
unlink(grid_prerequisite_file)
}comparable_aic <- cpue_sensitivity_summary |>
filter(.data$Model != sensitivity_labels[["drop_5yrs"]]) |>
pull(.data$`AIC-style score`)
aic_best <- min(comparable_aic, na.rm = TRUE)
aic_worst <- max(comparable_aic, na.rm = TRUE)
cpue_sensitivity_summary |>
mutate(
`AIC-style score` = kableExtra::cell_spec(
format_decimal(.data$`AIC-style score`, digits = 3),
background = case_when(
.data$Model != sensitivity_labels[["drop_5yrs"]] &
.data$`AIC-style score` == aic_best ~ "#C6EFCE",
.data$Model != sensitivity_labels[["drop_5yrs"]] &
.data$`AIC-style score` == aic_worst ~ "#FFC7CE",
TRUE ~ "transparent"
),
color = case_when(
.data$Model != sensitivity_labels[["drop_5yrs"]] &
.data$`AIC-style score` == aic_best ~ "#006100",
.data$Model != sensitivity_labels[["drop_5yrs"]] &
.data$`AIC-style score` == aic_worst ~ "#9C0006",
TRUE ~ "#000000"
),
bold = .data$Model != sensitivity_labels[["drop_5yrs"]] &
.data$`AIC-style score` %in% c(aic_best, aic_worst)
),
`Convergence code` = if_else(
is.na(`Convergence code`),
NA_character_,
as.character(as.integer(`Convergence code`))
),
across(where(is.numeric), ~ format_decimal(.x, digits = 3))
) |>
replace_missing() |>
kable(
align = c("l", rep("r", ncol(cpue_sensitivity_summary) - 1)),
escape = FALSE
)| Model | Convergence code | Penalized objective | AIC-style score | Max gradient | B0 | M0 | M4 | M10 | M30 | h | psi |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Base | 0 | 7,616.638 | 18,351.276 | 0.000 | 6,422,189.198 | 0.366 | 0.162 | 0.107 | 0.458 | 0.700 | 1.750 |
| Drop CPUE 2021-2025 | 0 | 7,616.967 | 18,351.934 | 0.000 | 6,319,130.368 | 0.368 | 0.163 | 0.108 | 0.460 | 0.700 | 1.750 |
| CPUE power 0.75 | 0 | 7,628.842 | 18,375.685 | 0.000 | 6,596,395.572 | 0.366 | 0.162 | 0.108 | 0.487 | 0.700 | 1.750 |
| CPUE q split from 2008 | 0 | 7,611.083 | 18,342.165 | 0.000 | 6,486,131.526 | 0.361 | 0.160 | 0.106 | 0.473 | 0.700 | 1.750 |
| CPUE index with LL1 selectivity | 0 | 7,624.576 | 18,367.151 | 0.000 | 6,415,834.105 | 0.367 | 0.162 | 0.108 | 0.456 | 0.700 | 1.750 |
| Constant CPUE CV 0.2 | 0 | 7,615.753 | 18,349.507 | 0.001 | 6,433,605.006 | 0.366 | 0.162 | 0.107 | 0.457 | 0.700 | 1.750 |
mle_review_summary |>
mutate(
Objective = format_decimal(.data$Objective, 6L),
`Max gradient` = format_decimal(.data$`Max gradient`, 8L),
Estimable = if_else(.data$Estimable, "Yes", "No"),
`Max raw harvest` = format_decimal(.data$`Max raw harvest`, 6L),
`Min abundance` = format_decimal(.data$`Min abundance`, 3L),
`Max catch error` = format_scientific(.data$`Max catch error`, 2L),
`Max continuation` = format_scientific(.data$`Max continuation`, 2L),
`Wall objective` = format_scientific(.data$`Wall objective`, 2L),
`MLE gate` = kableExtra::cell_spec(
.data$`MLE gate`,
background = if_else(.data$`MLE gate` == "Pass", "#C6EFCE", "#FFC7CE"),
color = if_else(.data$`MLE gate` == "Pass", "#006100", "#9C0006"),
bold = TRUE
)
) |>
replace_missing() |>
kable(
escape = FALSE,
align = c("l", rep("r", 4L), "c", rep("r", 5L), "c")
)| Model | Objective | Active parameters | Convergence | Max gradient | Estimable | Max raw harvest | Min abundance | Max catch error | Max continuation | Wall objective | MLE gate |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Base | 7,616.637822 | 1559 | 0 | 0.00000000 | Yes | 0.617533 | 1,696.516 | 2.53e-16 | 0.00e+00 | 2.57e-21 | Pass |
| No UAM | 7,616.449823 | 1559 | 0 | 0.00000000 | Yes | 0.618829 | 1,699.370 | 2.75e-16 | 0.00e+00 | 3.32e-21 | Pass |
| Drop CPUE 2021-2025 | 7,616.967176 | 1559 | 0 | 0.00000165 | Yes | 0.618054 | 1,737.077 | 3.02e-16 | 0.00e+00 | 2.85e-21 | Pass |
| CPUE power 0.75 | 7,628.842371 | 1559 | 0 | 0.00003724 | Yes | 0.578886 | 1,509.379 | 3.11e-16 | 0.00e+00 | 1.80e-24 | Pass |
| CPUE q split from 2008 | 7,611.082570 | 1560 | 0 | 0.00001233 | Yes | 0.590348 | 2,385.866 | 2.58e-16 | 0.00e+00 | 1.12e-23 | Pass |
| Annual LL1 selectivity 2023-2025 | 7,611.966180 | 1591 | 0 | 0.00000001 | Yes | 0.619410 | 1,702.603 | 2.24e-16 | 0.00e+00 | 3.73e-21 | Pass |
| No POP or HSP | 6,931.080345 | 1559 | 0 | 0.00000330 | Yes | 0.692190 | 505.499 | 3.20e-16 | 0.00e+00 | 7.83e-15 | Pass |
| No HSP | 7,445.149591 | 1559 | 0 | 0.00000286 | Yes | 0.622431 | 1,598.656 | 2.31e-16 | 0.00e+00 | 6.83e-21 | Pass |
| GTI with base aerial tau 0.59 | 7,638.201775 | 1559 | 0 | 0.00078151 | Yes | 0.627841 | 1,667.129 | 3.49e-16 | 0.00e+00 | 2.02e-20 | Pass |
| CPUE index with LL1 selectivity | 7,624.575617 | 1559 | 0 | 0.00000132 | Yes | 0.627057 | 1,341.578 | 2.92e-16 | 0.00e+00 | 1.72e-20 | Pass |
| Constant CPUE CV 0.2 | 7,615.753311 | 1559 | 0 | 0.00065463 | Yes | 0.618915 | 1,653.463 | 2.24e-16 | 0.00e+00 | 3.38e-21 | Pass |
| Direct four-anchor M | 7,608.440612 | 1561 | 0 | 0.00000132 | Yes | 0.592350 | 719.966 | 2.70e-16 | 0.00e+00 | 2.45e-23 | Pass |
| More flexible Indonesian selectivity | 8,129.683961 | 1559 | 0 | 0.00000000 | Yes | 0.579650 | 2,897.195 | 3.02e-16 | 0.00e+00 | 1.70e-24 | Pass |
| Estimated M-at-age slope | 7,616.622002 | 1560 | 0 | 0.00011990 | Yes | 0.618526 | 1,636.017 | 2.71e-16 | 0.00e+00 | 3.13e-21 | Pass |
| No conventional tags | 6,998.932988 | 1559 | 0 | 0.00002639 | Yes | 0.607662 | 2,386.215 | 7.93e-16 | 0.00e+00 | 3.56e-22 | Pass |
The noUAM fit zeros the NCNM additions in catch_UA.csv before constructing the model catch array. Drop_5yrs removes CPUE index observations for 2021, 2022, 2023, 2024, 2025; the separate CPUE length-frequency rows remain in the composition likelihood because OMMP16 specified the CPUE series rather than those rows. cpueom75 fixes the biomass-CPUE exponent at 0.75, and Cpueupq estimates separate catchability blocks beginning in 2008.
LL1_sel adds annual LL1 selectivity nodes for 2023-2025. Indo_sel applies the accepted 0.5/0.5/0.5 selectivity-hyperparameter working definition and is now an additional MCMC target. The NoPOPHSP transformation switches off both close-kin likelihoods exactly and is now an ordinary MLE and MCMC target under the approved global A10 wall. NoHSP switches off only HSPs. no_tags sets tag_switch = 0, excluding the complete conventional-tag likelihood while retaining the tag tables required to define model dimensions; its MLE and accepted MCMC are reported here. The modified Troll case activates the grid-type trolling index and fixes aerial-survey tau at 0.59, exactly the base value. It therefore does not implement the increase requested by OMMP16. Because the base tau was already adjusted to improve aerial SDNR, this sensitivity isolates only the trolling-index effect and is accepted under decision esc31_2026_troll_base_aerial_tau_059_v1. A future additional aerial-CV change would be a new sensitivity. U_Fleet uses LL1 selectivity for the CPUE abundance-index prediction while retaining fishery-7 selectivity for the CPUE length-composition likelihood.
The revised base reads the raw GAM22 B=1000 estimation CV from the Index sheet in CV_GAM22_20260713.xlsx, adds it to a 20% baseline in log-variance space, and fixes the additional CPUE log-SD at 0.2. The constant-CV sensitivity removes the GAM22 component and restores the historical raw-CPUE plus fixed-log-SD-0.2 implementation exactly.
The final length_m registry slot now switches from the length-based base to the direct four-anchor mortality parameterization. It estimates M0, M4, M10, and M30 while preserving that slot’s path and seed. The appended estimate_m_slope sensitivity retains the base length-based parameterization but estimates \(m_c\); it does not change the base decision to fix \(m_c=-1\).
cpue_diagnostic_keys <- cpue_sensitivity_keys
cpue_diagnostic_fits <- c(
list(Base = base_fit),
sensitivity_fits[cpue_diagnostic_keys]
)
names(cpue_diagnostic_fits) <- c(
"Base",
unname(sensitivity_labels[cpue_diagnostic_keys])
)
cpue_fit_plots <- imap(cpue_diagnostic_fits, function(fit, label) {
plot_cpue(fit) + ggtitle(label)
})
patchwork::wrap_plots(cpue_fit_plots, ncol = 2, guides = "collect") &
theme(legend.position = "bottom")
sbt::plot_biomass_spawning(
c(list(base_fit), unname(sensitivity_fits)),
labels = c(
"Base",
unname(sensitivity_labels[names(sensitivity_fits)])
),
posterior = FALSE
)
comparison_reports <- lapply(mle_fit_sources, sbt_fit_report)
comparison_labels <- c(base = "Base", sensitivity_labels)
spawning_biomass_comparison <- imap_dfr(
mle_fit_sources,
function(fit, key) {
report <- comparison_reports[[key]]
tibble(
Year = seq.int(
fit$data$first_yr,
length.out = length(report$spawning_biomass_y)
),
TRO = as.numeric(report$spawning_biomass_y),
Model = unname(comparison_labels[[key]])
)
}
)
recruitment_comparison <- imap_dfr(
mle_fit_sources,
function(fit, key) {
report <- comparison_reports[[key]]
tibble(
Year = seq.int(
fit$data$first_yr,
length.out = length(report$recruitment_y)
),
Recruitment = as.numeric(report$recruitment_y) / 1e6,
Model = unname(comparison_labels[[key]])
)
}
)
mortality_comparison <- imap_dfr(
mle_fit_sources,
function(fit, key) {
tibble(
Age = fit$data$age_a,
`Natural mortality` = as.numeric(comparison_reports[[key]]$M_a),
Model = unname(comparison_labels[[key]])
)
}
)
previous_spawning_biomass <- tibble(
Year = seq.int(
previous_v1_fit$data$first_yr,
length.out = length(previous_v1_report$spawning_biomass_y)
),
TRO = as.numeric(previous_v1_report$spawning_biomass_y)
)
previous_recruitment <- tibble(
Year = seq.int(
previous_v1_fit$data$first_yr,
length.out = length(previous_v1_report$recruitment_y)
),
Recruitment = as.numeric(previous_v1_report$recruitment_y) / 1e6
)
previous_mortality <- tibble(
Age = previous_v1_fit$data$min_age:previous_v1_fit$data$max_age,
`Natural mortality` = as.numeric(previous_v1_report$M_a)
)ggplot(
spawning_biomass_comparison,
aes(.data$Year, .data$TRO, colour = .data$Model)
) +
geom_line(linewidth = 0.65, alpha = 0.9) +
geom_line(
data = previous_spawning_biomass,
colour = "black",
linetype = "dashed",
linewidth = 0.9
) +
scale_y_continuous(
labels = label_comma(),
limits = c(0, NA),
expand = expansion(mult = c(0, 0.04))
) +
labs(
x = NULL,
y = "TRO",
colour = NULL
) +
theme(legend.position = "bottom")
ggplot(
recruitment_comparison,
aes(.data$Year, .data$Recruitment, colour = .data$Model)
) +
geom_line(linewidth = 0.65, alpha = 0.9) +
geom_line(
data = previous_recruitment,
colour = "black",
linetype = "dashed",
linewidth = 0.9
) +
scale_y_continuous(
limits = c(0, NA),
expand = expansion(mult = c(0, 0.04))
) +
labs(x = NULL, y = "Recruitment (millions)", colour = NULL) +
theme(legend.position = "bottom")
ggplot(
mortality_comparison,
aes(.data$Age, .data$`Natural mortality`, colour = .data$Model)
) +
geom_line(linewidth = 0.7, alpha = 0.9) +
geom_line(
data = previous_mortality,
colour = "black",
linetype = "dashed",
linewidth = 0.9
) +
scale_y_continuous(
limits = c(0, NA),
expand = expansion(mult = c(0, 0.04))
) +
labs(x = "Age", y = "Natural mortality", colour = NULL) +
theme(legend.position = "bottom")
mortality_comparison |>
filter(.data$Model %in% c("Base", "Direct four-anchor M")) |>
ggplot(aes(
.data$Age,
.data$`Natural mortality`,
colour = .data$Model
)) +
geom_line(linewidth = 0.9) +
geom_line(
data = previous_mortality,
colour = "black",
linetype = "dashed",
linewidth = 0.9
) +
scale_y_continuous(
limits = c(0, NA),
expand = expansion(mult = c(0, 0.04))
) +
labs(x = "Age", y = "Natural mortality", colour = NULL) +
theme(legend.position = "bottom")
sensitivity_harvest_wall_summary <- imap_dfr(
mle_fit_sources,
function(fit, key) {
wall <- sensitivity_harvest_wall_fit_diagnostics(fit)
tibble(
Model = if (identical(key, "base")) "Base" else sensitivity_labels[[key]],
`Max raw harvest` = wall$maximum_raw_harvest %||% NA_real_,
`Cells above 0.85` = wall$cells_above_onset %||% NA_integer_,
`Wall objective` = wall$total_wall_objective %||% NA_real_,
`Max recomputation error` = max(
wall$maximum_penalty_recomputation_error %||% NA_real_,
wall$total_penalty_recomputation_error %||% NA_real_,
na.rm = TRUE
),
`Wall gate` = if (isTRUE(wall$passes)) "Pass" else "Fail"
)
}
)
stopifnot(all(sensitivity_harvest_wall_summary$`Wall gate` == "Pass"))
sensitivity_harvest_wall_summary |>
mutate(
`Max raw harvest` = format_decimal(.data$`Max raw harvest`, 6L),
`Wall objective` = format_decimal(.data$`Wall objective`, 8L),
`Max recomputation error` =
format_decimal(.data$`Max recomputation error`, 10L)
) |>
replace_missing() |>
kable(align = c("l", rep("r", 4L), "l"))| Model | Max raw harvest | Cells above 0.85 | Wall objective | Max recomputation error | Wall gate |
|---|---|---|---|---|---|
| Base | 0.617533 | 0 | 0.00000000 | 0.0000000000 | Pass |
| No UAM | 0.618829 | 0 | 0.00000000 | 0.0000000000 | Pass |
| Drop CPUE 2021-2025 | 0.618054 | 0 | 0.00000000 | 0.0000000000 | Pass |
| CPUE power 0.75 | 0.578886 | 0 | 0.00000000 | 0.0000000000 | Pass |
| CPUE q split from 2008 | 0.590348 | 0 | 0.00000000 | 0.0000000000 | Pass |
| Annual LL1 selectivity 2023-2025 | 0.619410 | 0 | 0.00000000 | 0.0000000000 | Pass |
| No POP or HSP | 0.692190 | 0 | 0.00000000 | 0.0000000000 | Pass |
| No HSP | 0.622431 | 0 | 0.00000000 | 0.0000000000 | Pass |
| GTI with base aerial tau 0.59 | 0.627841 | 0 | 0.00000000 | 0.0000000000 | Pass |
| CPUE index with LL1 selectivity | 0.627057 | 0 | 0.00000000 | 0.0000000000 | Pass |
| Constant CPUE CV 0.2 | 0.618915 | 0 | 0.00000000 | 0.0000000000 | Pass |
| Direct four-anchor M | 0.592350 | 0 | 0.00000000 | 0.0000000000 | Pass |
| More flexible Indonesian selectivity | 0.579650 | 0 | 0.00000000 | 0.0000000000 | Pass |
| Estimated M-at-age slope | 0.618526 | 0 | 0.00000000 | 0.0000000000 | Pass |
| No conventional tags | 0.607662 | 0 | 0.00000000 | 0.0000000000 | Pass |
The base model and all fourteen registered sensitivities use the same four-chain, 150-warmup, dense-metric sampling structure and transition controls. Retained chain lengths are increased for models whose first production run narrowly missed a diagnostic, as shown below. The conventional-tag-exclusion posterior uses 1,800 retained draws per chain after its initial shorter run missed the R-hat and bulk-ESS gates. Its final run passes every standard acceptance gate. NoPOPHSP remains in its original registry slot so later seeds remain unchanged. The standard acceptance rule requires maximum rank-normalized R-hat below 1.01, minimum bulk and tail ESS of at least 400, an estimable MLE, and no divergences, maximum-treedepth hits, or biologically invalid retained draws. The longer NoPOPHSP run has 12,000 retained draws and passes every one of those gates except for two divergent transitions. Decision esc31_2026_no_pop_hsp_accept_two_divergences_12000_draws_v1 explicitly accepts those two divergences for this sensitivity only. It does not change the standard rule for another sensitivity or any grid cell. Every retained draw is checked once at fit completion for raw combined seasonal harvest no greater than 0.9, positive population states, exact catch accounting, and zero continuation penalty within numerical tolerance. The resulting validation record is stored in the same sbt_fit and is bound to its model, posterior, sampler diagnostics, scientific contract, thresholds, and diagnostic checksums.
The final seed-73107 recovery candidate is not included in the table or posterior comparison because it failed maximum R-hat and both ESS gates, despite zero divergences and a complete passing state audit.
The direct-anchor, Indonesian-selectivity, estimated-slope, and conventional-tag-exclusion posteriors are reported with the other sensitivities. Sensitivity results do not gate the MCMC or MLE grids; no grid is launched from this page.
An ordinary report render never launches a sampler. It loads only compatible posteriors already embedded in the canonical sbt_fit files and reuses their payload-bound validation records. ESC31_FIT_VALIDATION_MODE=full deliberately repeats the complete audit, while skip trusts the stored record without recalculating its compact checksums. Missing or incompatible posteriors remain pending. ESC31_RUN_SENSITIVITY_MCMC=true is reserved for an explicit MCMC worker run.
mcmc_model_labels <- c(base = "Base", sensitivity_labels)
stopifnot(
identical(names(mcmc_model_labels), sensitivity_mcmc_registry_keys)
)
sensitivity_mcmc_contract <- imap_dfr(
mcmc_model_labels[sensitivity_mcmc_registry_keys],
function(model_label, key) {
config <- sensitivity_mcmc_config(key)
tibble(
Model = model_label,
Chains = config$chains,
Warmup = config$num_warmup,
`Retained per chain` = config$num_samples,
`Total retained` = config$chains * config$num_samples,
Metric = stringr::str_to_title(config$metric),
`Adapt delta` = config$adapt_delta,
`Max treedepth` = config$max_treedepth
)
}
)
extended_mcmc_rows <- which(
sensitivity_mcmc_contract$`Retained per chain` !=
sensitivity_mcmc_default_config$num_samples
)
sensitivity_mcmc_contract |>
mutate(
across(
c(Chains, Warmup, `Retained per chain`, `Total retained`, `Max treedepth`),
~ format_decimal(.x, digits = 0L)
),
`Adapt delta` = format_decimal(.data$`Adapt delta`, digits = 3L)
) |>
kable(align = c("l", rep("r", 4L), "c", "r", "r")) |>
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover", "condensed"),
full_width = TRUE
) |>
kableExtra::row_spec(extended_mcmc_rows, background = "#DDEBF7")| Model | Chains | Warmup | Retained per chain | Total retained | Metric | Adapt delta | Max treedepth |
|---|---|---|---|---|---|---|---|
| Base | 4 | 150 | 750 | 3,000 | Dense | 0.999 | 13 |
| No UAM | 4 | 150 | 750 | 3,000 | Dense | 0.999 | 13 |
| Drop CPUE 2021-2025 | 4 | 150 | 750 | 3,000 | Dense | 0.999 | 13 |
| CPUE power 0.75 | 4 | 150 | 750 | 3,000 | Dense | 0.999 | 13 |
| CPUE q split from 2008 | 4 | 150 | 750 | 3,000 | Dense | 0.999 | 13 |
| Annual LL1 selectivity 2023-2025 | 4 | 150 | 900 | 3,600 | Dense | 0.999 | 13 |
| No POP or HSP | 4 | 150 | 3,000 | 12,000 | Dense | 0.999 | 13 |
| No HSP | 4 | 150 | 900 | 3,600 | Dense | 0.999 | 13 |
| GTI with base aerial tau 0.59 | 4 | 150 | 1,200 | 4,800 | Dense | 0.999 | 13 |
| CPUE index with LL1 selectivity | 4 | 150 | 750 | 3,000 | Dense | 0.999 | 13 |
| Constant CPUE CV 0.2 | 4 | 150 | 900 | 3,600 | Dense | 0.999 | 13 |
| Direct four-anchor M | 4 | 150 | 750 | 3,000 | Dense | 0.999 | 13 |
| More flexible Indonesian selectivity | 4 | 150 | 750 | 3,000 | Dense | 0.999 | 13 |
| Estimated M-at-age slope | 4 | 150 | 900 | 3,600 | Dense | 0.999 | 13 |
| No conventional tags | 4 | 150 | 1,800 | 7,200 | Dense | 0.999 | 13 |
empty_sensitivity_mcmc_row <- function(model, status) {
tibble(
Model = model,
Status = status,
`Max R-hat` = NA_real_,
`Min bulk ESS` = NA_real_,
`Min tail ESS` = NA_real_,
Divergences = NA_integer_,
`Max treedepth hits` = NA_integer_,
`MLE convergence` = NA_integer_,
Estimable = NA,
`State draws expected` = NA_integer_,
`State draws checked` = NA_integer_,
`Invalid state draws` = NA_integer_,
`Non-finite state draws` = NA_integer_,
`Max raw harvest` = NA_real_,
`Min abundance` = NA_real_,
`Max state penalty` = NA_real_
)
}
sensitivity_mcmc_summary <- imap_dfr(
mcmc_model_labels[sensitivity_mcmc_registry_keys],
function(model_label, key) {
if (key %in% rejected_sensitivity_mle_keys) {
return(empty_sensitivity_mcmc_row(model_label, "MLE rejected"))
}
fit <- sensitivity_mcmc_fits[[key]]
if (is.null(fit)) {
status <- if (run_sensitivity_mcmc) {
"Not run"
} else {
"Pending current MCMC"
}
return(empty_sensitivity_mcmc_row(model_label, status))
}
diagnostics <- fit$fit$diagnostics$mcmc
tibble(
Model = model_label,
Status = sensitivity_mcmc_acceptance_label(diagnostics),
`Max R-hat` = diagnostics$max_rhat,
`Min bulk ESS` = diagnostics$min_bulk_ess,
`Min tail ESS` = diagnostics$min_tail_ess,
Divergences = as.integer(diagnostics$divergences),
`Max treedepth hits` = as.integer(diagnostics$max_treedepth_hits),
`MLE convergence` = as.integer(diagnostics$convergence),
Estimable = isTRUE(diagnostics$estimable),
`State draws expected` = as.integer(diagnostics$state_draws_expected),
`State draws checked` = as.integer(diagnostics$state_draws_evaluated),
`Invalid state draws` = as.integer(diagnostics$state_invalid_draws),
`Non-finite state draws` =
as.integer(diagnostics$state_non_finite_draws),
`Max raw harvest` = diagnostics$state_max_raw_harvest,
`Min abundance` = diagnostics$state_min_number,
`Max state penalty` = diagnostics$state_max_harvest_penalty
)
}
)format_mcmc_cell <- function(
value, failed, digits = 0L, accepted_exception = FALSE) {
label <- ifelse(
is.na(value),
"-",
formatC(value, format = "f", digits = digits, big.mark = ",")
)
kableExtra::cell_spec(
label,
background = ifelse(
!is.na(accepted_exception) & accepted_exception,
"#FFF2CC",
ifelse(!is.na(failed) & failed, "#FFC7CE", "transparent")
),
color = ifelse(
!is.na(accepted_exception) & accepted_exception,
"#7F6000",
ifelse(!is.na(failed) & failed, "#9C0006", "#000000")
),
bold = (!is.na(failed) & failed) |
(!is.na(accepted_exception) & accepted_exception)
)
}
sensitivity_mcmc_summary |>
mutate(
accepted_divergence_exception =
.data$Status == "Accepted (2 divergences)",
Status = kableExtra::cell_spec(
.data$Status,
background = case_when(
.data$Status == "Pass" ~ "#C6EFCE",
.data$Status == "Accepted (2 divergences)" ~ "#FFF2CC",
.data$Status == "Fail" ~ "#FFC7CE",
.data$Status == "MLE rejected" ~ "#E7E6E6",
TRUE ~ "transparent"
),
color = case_when(
.data$Status == "Pass" ~ "#006100",
.data$Status == "Accepted (2 divergences)" ~ "#7F6000",
.data$Status == "Fail" ~ "#9C0006",
.data$Status == "MLE rejected" ~ "#7F6000",
TRUE ~ "#000000"
),
bold = .data$Status %in%
c("Pass", "Accepted (2 divergences)", "Fail", "MLE rejected")
),
`Max R-hat` = format_mcmc_cell(
`Max R-hat`,
`Max R-hat` >= sensitivity_mcmc_thresholds$max_rhat,
digits = 3L
),
`Min bulk ESS` = format_mcmc_cell(
`Min bulk ESS`,
`Min bulk ESS` < sensitivity_mcmc_thresholds$min_bulk_ess
),
`Min tail ESS` = format_mcmc_cell(
`Min tail ESS`,
`Min tail ESS` < sensitivity_mcmc_thresholds$min_tail_ess
),
Divergences = format_mcmc_cell(
Divergences,
Divergences > sensitivity_mcmc_thresholds$divergences &
!.data$accepted_divergence_exception,
accepted_exception = .data$accepted_divergence_exception
),
`Max treedepth hits` = format_mcmc_cell(
`Max treedepth hits`,
`Max treedepth hits` > sensitivity_mcmc_thresholds$max_treedepth_hits
),
`MLE convergence` = format_mcmc_cell(
`MLE convergence`,
`MLE convergence` != 0L
),
Estimable = kableExtra::cell_spec(
ifelse(is.na(Estimable), "-", ifelse(Estimable, "Yes", "No")),
background = ifelse(!is.na(Estimable) & !Estimable, "#FFC7CE", "transparent"),
color = ifelse(!is.na(Estimable) & !Estimable, "#9C0006", "#000000"),
bold = !is.na(Estimable) & !Estimable
),
`State draws checked` = format_mcmc_cell(
`State draws checked`,
is.na(`State draws expected`) |
`State draws checked` != `State draws expected`
),
`State draws expected` = format_mcmc_cell(
`State draws expected`,
`State draws expected` <= 0L
),
`Invalid state draws` = format_mcmc_cell(
`Invalid state draws`,
`Invalid state draws` > 0L
),
`Non-finite state draws` = format_mcmc_cell(
`Non-finite state draws`,
`Non-finite state draws` > 0L
),
`Max raw harvest` = format_mcmc_cell(
`Max raw harvest`,
`Max raw harvest` >
sensitivity_mcmc_thresholds$biological_state$hrate_limit +
sensitivity_mcmc_thresholds$biological_state$hrate_tolerance,
digits = 4L
),
`Min abundance` = format_mcmc_cell(
`Min abundance`,
`Min abundance` <=
sensitivity_mcmc_thresholds$biological_state$number_tolerance,
digits = 3L
),
`Max state penalty` = format_mcmc_cell(
`Max state penalty`,
abs(`Max state penalty`) >
sensitivity_mcmc_thresholds$biological_state$penalty_tolerance,
digits = 3L
)
) |>
select(-accepted_divergence_exception) |>
kable(
escape = FALSE,
align = c("l", "l", rep("r", 6), "c", rep("r", 7))
)| Model | Status | Max R-hat | Min bulk ESS | Min tail ESS | Divergences | Max treedepth hits | MLE convergence | Estimable | State draws expected | State draws checked | Invalid state draws | Non-finite state draws | Max raw harvest | Min abundance | Max state penalty |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Base | Pass | 1.009 | 1,157 | 1,178 | 0 | 0 | 0 | Yes | 3,000 | 3,000 | 0 | 0 | 0.8787 | 490.556 | 0.000 |
| No UAM | Pass | 1.009 | 1,058 | 1,407 | 0 | 0 | 0 | Yes | 3,000 | 3,000 | 0 | 0 | 0.8730 | 378.619 | 0.000 |
| Drop CPUE 2021-2025 | Pass | 1.009 | 925 | 857 | 0 | 0 | 0 | Yes | 3,000 | 3,000 | 0 | 0 | 0.8713 | 520.517 | 0.000 |
| CPUE power 0.75 | Pass | 1.008 | 1,271 | 1,397 | 0 | 0 | 0 | Yes | 3,000 | 3,000 | 0 | 0 | 0.8837 | 434.136 | 0.000 |
| CPUE q split from 2008 | Pass | 1.010 | 1,148 | 1,394 | 0 | 0 | 0 | Yes | 3,000 | 3,000 | 0 | 0 | 0.8622 | 434.374 | 0.000 |
| Annual LL1 selectivity 2023-2025 | Pass | 1.006 | 1,173 | 977 | 0 | 0 | 0 | Yes | 3,600 | 3,600 | 0 | 0 | 0.8691 | 496.067 | 0.000 |
| No POP or HSP | Accepted (2 divergences) | 1.006 | 735 | 2,094 | 2 | 0 | 0 | Yes | 12,000 | 12,000 | 0 | 0 | 0.8870 | 95.790 | 0.000 |
| No HSP | Pass | 1.010 | 1,076 | 1,556 | 0 | 0 | 0 | Yes | 3,600 | 3,600 | 0 | 0 | 0.8815 | 479.011 | 0.000 |
| GTI with base aerial tau 0.59 | Pass | 1.007 | 1,613 | 1,592 | 0 | 0 | 0 | Yes | 4,800 | 4,800 | 0 | 0 | 0.8833 | 543.472 | 0.000 |
| CPUE index with LL1 selectivity | Pass | 1.008 | 966 | 761 | 0 | 0 | 0 | Yes | 3,000 | 3,000 | 0 | 0 | 0.8698 | 354.700 | 0.000 |
| Constant CPUE CV 0.2 | Pass | 1.007 | 1,391 | 1,414 | 0 | 0 | 0 | Yes | 3,600 | 3,600 | 0 | 0 | 0.8739 | 549.147 | 0.000 |
| Direct four-anchor M | Pass | 1.010 | 945 | 1,387 | 0 | 0 | 0 | Yes | 3,000 | 3,000 | 0 | 0 | 0.8728 | 149.360 | 0.000 |
| More flexible Indonesian selectivity | Pass | 1.008 | 1,127 | 910 | 0 | 0 | 0 | Yes | 3,000 | 3,000 | 0 | 0 | 0.8775 | 558.173 | 0.000 |
| Estimated M-at-age slope | Pass | 1.009 | 1,168 | 1,162 | 0 | 0 | 0 | Yes | 3,600 | 3,600 | 0 | 0 | 0.8731 | 309.740 | 0.000 |
| No conventional tags | Pass | 1.007 | 688 | 421 | 0 | 0 | 0 | Yes | 7,200 | 7,200 | 0 | 0 | 0.8733 | 456.635 | 0.000 |
completed_mcmc_diagnostics <- sensitivity_mcmc_summary |>
filter(.data$Status %in%
c("Pass", "Accepted (2 divergences)", "Fail"))
if (nrow(completed_mcmc_diagnostics)) {
mcmc_diagnostic_plot_data <- completed_mcmc_diagnostics |>
select(
.data$Model, .data$Status, `Max R-hat`, `Min bulk ESS`,
`Min tail ESS`, .data$Divergences, `Max treedepth hits`
) |>
pivot_longer(
cols = -c(.data$Model, .data$Status),
names_to = "Metric",
values_to = "Value"
) |>
mutate(
Metric = factor(
.data$Metric,
levels = c(
"Max R-hat", "Min bulk ESS", "Min tail ESS", "Divergences",
"Max treedepth hits"
)
)
)
mcmc_diagnostic_thresholds <- tibble(
Metric = factor(
c(
"Max R-hat", "Min bulk ESS", "Min tail ESS", "Divergences",
"Max treedepth hits"
),
levels = levels(mcmc_diagnostic_plot_data$Metric)
),
Threshold = c(1.01, 400, 400, 0, 0)
)
ggplot(
mcmc_diagnostic_plot_data,
aes(x = .data$Value, y = forcats::fct_rev(factor(.data$Model)), color = .data$Status)
) +
geom_vline(
data = mcmc_diagnostic_thresholds,
aes(xintercept = .data$Threshold),
linetype = "dashed",
color = "grey35",
inherit.aes = FALSE
) +
geom_point(size = 2.2) +
facet_wrap(vars(.data$Metric), scales = "free_x", ncol = 2) +
scale_color_manual(values = c(
Pass = "#0072B2",
`Accepted (2 divergences)` = "#E69F00",
Fail = "#D55E00"
)) +
labs(x = NULL, y = NULL, color = NULL) +
theme(legend.position = "bottom")
} else {
knitr::asis_output(
"*No current posterior is available for diagnostic comparison.*"
)
}
This is the current counterpart of the stock-status sensitivity table in the 2023 assessment. It reports every accepted posterior, rather than only the selected examples discussed in the assessment-paper text. The relative TRO and relative age-10-plus biomass columns use the 2026 model state. Equilibrium MSY quantities use 2025, the latest fitted catch year. Each row uses all retained draws from that accepted posterior.
sensitivity_status_file <- file.path(
esc_dir,
"report_data", "table5_sensitivity_stock_status.csv"
)
if (!file.exists(sensitivity_status_file)) {
stop(
"The web stock-status tables are missing. Run ",
"`Rscript scripts/build-esc31-web-status-tables.R .`.",
call. = FALSE
)
}
sensitivity_stock_status <- readr::read_csv(
sensitivity_status_file,
show_col_types = FALSE
)
stopifnot(
identical(
sensitivity_stock_status$Model_key,
names(sensitivity_mcmc_fits)[vapply(
sensitivity_mcmc_fits,
function(fit) inherits(fit, "sbt_fit") && !is.null(fit$mcmc),
logical(1)
)]
)
)
format_status_interval <- function(median, lower, upper, digits) {
paste0(
format_decimal(median, digits),
" (", format_decimal(lower, digits),
"-", format_decimal(upper, digits), ")"
)
}
sensitivity_stock_status |>
transmute(
Model = .data$Model,
`Relative TRO (2026)` = format_status_interval(
.data$Relative_TRO_median,
.data$Relative_TRO_lower,
.data$Relative_TRO_upper,
3L
),
`Relative B10+ (2026)` = format_status_interval(
.data$Relative_B10_plus_median,
.data$Relative_B10_plus_lower,
.data$Relative_B10_plus_upper,
3L
),
`F/F_MSY (2025)` = format_status_interval(
.data$F_F_MSY_median,
.data$F_F_MSY_lower,
.data$F_F_MSY_upper,
3L
),
`TRO/TRO_MSY (2025)` = format_status_interval(
.data$TRO_TRO_MSY_median,
.data$TRO_TRO_MSY_lower,
.data$TRO_TRO_MSY_upper,
3L
),
`MSY (t; 2025)` = format_status_interval(
.data$MSY_t_median,
.data$MSY_t_lower,
.data$MSY_t_upper,
0L
)
) |>
kable(align = c("l", rep("r", 5L)))| Model | Relative TRO (2026) | Relative B10+ (2026) | F/F_MSY (2025) | TRO/TRO_MSY (2025) | MSY (t; 2025) |
|---|---|---|---|---|---|
| Base | 0.282 (0.222-0.353) | 0.232 (0.180-0.294) | 0.553 (0.466-0.649) | 0.987 (0.779-1.237) | 32,451 (26,663-39,862) |
| No UAM | 0.290 (0.233-0.362) | 0.240 (0.189-0.301) | 0.498 (0.424-0.588) | 1.014 (0.811-1.265) | 31,825 (26,315-38,504) |
| Drop CPUE 2021-2025 | 0.270 (0.214-0.338) | 0.223 (0.176-0.281) | 0.592 (0.491-0.710) | 0.952 (0.758-1.188) | 32,336 (26,708-39,587) |
| CPUE power 0.75 | 0.292 (0.233-0.366) | 0.243 (0.191-0.307) | 0.521 (0.441-0.615) | 1.023 (0.811-1.276) | 33,787 (27,334-41,498) |
| CPUE q split from 2008 | 0.239 (0.182-0.305) | 0.192 (0.143-0.250) | 0.611 (0.511-0.733) | 0.831 (0.633-1.062) | 32,423 (26,587-39,954) |
| Annual LL1 selectivity 2023-2025 | 0.283 (0.225-0.356) | 0.230 (0.181-0.293) | 0.570 (0.479-0.673) | 0.981 (0.778-1.233) | 31,835 (26,114-38,887) |
| No POP or HSP | 0.266 (0.201-0.349) | 0.216 (0.159-0.292) | 0.586 (0.476-0.716) | 0.931 (0.701-1.229) | 32,606 (26,684-39,950) |
| No HSP | 0.287 (0.226-0.363) | 0.238 (0.184-0.304) | 0.548 (0.459-0.654) | 1.007 (0.793-1.279) | 32,655 (26,608-39,735) |
| GTI with base aerial tau 0.59 | 0.287 (0.229-0.357) | 0.246 (0.194-0.308) | 0.629 (0.540-0.735) | 1.024 (0.817-1.271) | 31,686 (26,251-38,301) |
| CPUE index with LL1 selectivity | 0.294 (0.233-0.363) | 0.242 (0.190-0.304) | 0.536 (0.450-0.634) | 1.027 (0.815-1.271) | 32,597 (26,589-39,760) |
| Constant CPUE CV 0.2 | 0.283 (0.224-0.352) | 0.232 (0.182-0.292) | 0.545 (0.459-0.643) | 0.989 (0.781-1.236) | 32,557 (26,742-40,189) |
| Direct four-anchor M | 0.339 (0.266-0.432) | 0.258 (0.204-0.331) | 0.471 (0.385-0.565) | 1.153 (0.911-1.461) | 33,435 (27,070-41,387) |
| More flexible Indonesian selectivity | 0.279 (0.224-0.345) | 0.230 (0.183-0.288) | 0.539 (0.458-0.639) | 0.975 (0.783-1.201) | 32,246 (26,704-39,756) |
| Estimated M-at-age slope | 0.287 (0.228-0.362) | 0.234 (0.184-0.295) | 0.546 (0.460-0.652) | 1.004 (0.795-1.268) | 32,751 (26,838-40,262) |
| No conventional tags | 0.273 (0.214-0.345) | 0.226 (0.175-0.288) | 0.543 (0.459-0.645) | 0.950 (0.742-1.207) | 32,503 (26,655-39,911) |
The complete numeric values, retained-draw counts, acceptance labels, and source identities are available in the machine-readable Table 5 CSV and the shared stock-status provenance file.
passing_mcmc_fits <- sensitivity_mcmc_fits[vapply(
sensitivity_mcmc_fits,
function(fit) {
!is.null(fit) && sensitivity_mcmc_passes(
fit$fit$diagnostics$mcmc
)
},
logical(1)
)]
if (length(passing_mcmc_fits)) {
sbt::plot_biomass_spawning(
unname(passing_mcmc_fits),
labels = unname(mcmc_model_labels[names(passing_mcmc_fits)])
)
} else {
knitr::asis_output(
"*No current posterior satisfies an accepted MCMC decision.*"
)
}
cpue_osa_plots <- imap(cpue_diagnostic_fits, function(fit, label) {
plot_cpue_osa_quiet(fit) + ggtitle(label)
})
patchwork::wrap_plots(cpue_osa_plots, ncol = 2)
The accepted working Indo_sel definition responds to the OMMP16 direction to restore bimodality through greater flexibility. OMMP16 gives no numeric hyperparameters, so the selected 0.5/0.5/0.5 values and the resulting posterior should remain visible as an assessment assumption. This appended target does not alter the base-only grid-prerequisite decision.
The expanded stock-status table closes the management-quantity interpretation for both appended targets. Relative to the base medians, Indo_sel changes relative TRO by -0.002, relative B10+ by -0.002, F/F_MSY by -0.014, TRO/TRO_MSY by -0.011, and MSY by -205 t. The estimated-slope case changes the same medians by +0.005, +0.002, -0.007, +0.017, and +300 t, respectively. Their 95% credible intervals overlap the base intervals broadly, and these shifts are small relative to the within-posterior uncertainty shown in the table. Neither appended case identifies a management-relevant change that would justify refitting the accepted base; they remain explicit accepted sensitivity interpretations rather than new base assumptions.