Decision esc31_2026_retain_frozen_1996_fishery7_length_composition_v1 explicitly retains the checksum-bound 1996 Fishery 7 length composition already used by the accepted base, sensitivities, grids, and projection inputs. The pre-sign-off audit found 23 duplicated Japanese size-data keys at the complete year/month/adjusted-area/length level. Correct aggregation changes only the 1996 composition by more than \(10^{-6}\), with a maximum absolute bin-proportion change of 0.000529. This small, localized discrepancy does not justify rebuilding the ESC31 assessment chain. No accepted assessment input or fitted artifact is changed by this decision.
At the next planned data refresh, raw frequencies must first be summed over the complete key and then normalized once within each spatiotemporal cell. That corrected construction must enter only as part of the next deliberate data revision, with new checksums, provenance, and downstream validation. Retaining the frozen ESC31 vector does not endorse duplicate records; it preserves the scientific identity of the assessment already reviewed and accepted. The duplicate-safe input-build work is tracked in GitHub issue 6.
LL4 is now a standard selectivity-based fishery with a length-composition likelihood and exact total-catch conditioning. LL4 has one selectivity block beginning in 1953, estimates ages 8–21, and uses the same selectivity hyperparameters as LL3 (rho_year = 0.5, rho_age = 0.5, sigma = 0.75).
A preventive harvest wall begins at 0.85, with strength 10, ceiling 0.90, and scale 0.01. Accepted fits and draws must still satisfy the biological-state, catch-accounting, and zero-continuation gates. The model estimates natural mortality at \(M_{10}\) and \(M_{30}\) and uses the fixed length-scaled curve \(M_a=M_{10}(L_{10}/L_a)\) through age 25 before increasing linearly to \(M_{30}\) at age 30. The explicit exponent parameter \(m_c\) is fixed at -1. Every retained posterior draw must have finite, positive mortality at every age; no mortality exception is allowed.
The fixed base selectivity settings are deliberate ESC31 choices: LL2 age/year/sigma is 0.50/0.70/0.50, Indonesia is 0.98/0.98/0.12, and CPUE is 0.90/0.95/0.23. The 1976 Indonesian change node is retained. The base CPUE uncertainty deliberately combines fixed log-SD 0.2 with the raw annual GAM22 CVs. These documented departures from OMMP16 Table 1 and its constant-20% reference case are accepted without another base refit.
For 2025, posterior median total reproductive output relative to its MSY reference point is TRO/TRO_MSY = 0.987 (95% credible interval 0.779–1.237), with a 54.6% posterior probability of being below one. Fishing mortality is clearly below its reference point: F/F_MSY = 0.553 (0.466–0.649), with no retained draw above one.
The four-chain MCMC is accepted: maximum R-hat is 1.0087, minimum bulk and tail ESS are 1,157 and 1,178, with 0 divergences and 0 maximum-treedepth hits. The fixed LL2, Indonesian, and CPUE selectivity settings, retained 1976 Indonesian node, and combined CPUE uncertainty are also accepted as deliberate base choices; no refit is required for those documented deviations. The main reservations are structural rather than numerical. The largest sensitivity shifts are lower recent relative TRO under the CPUE catchability change from 2008 and higher relative TRO under direct four-anchor mortality. Estimating the length-mortality exponent and relaxing Indonesian selectivity remain close to the base. All sensitivities are accepted, with the explicit two-divergence NoPOPHSP exception, and all nine selected grid cells pass their production gates. The balanced grid posterior, direct-get_M 108-fit MLE grid, and resample have now also been built and independently validated. The four 2,000-draw projection arms were completed and accepted under the preceding projection contract. The 3 August format-9 update changes only the staging of not-yet-decided future catch and records realized catch separately. The completed fail-closed paths had no catch shortfall, zero harvest penalties, and maximum raw harvest below 0.57, so their exact fingerprint-bound format-8 outputs remain accepted for reporting. Future projection runs use format 9. The sensitivity page evaluates both an estimated mortality exponent and the direct four-anchor mortality alternative.
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(kableExtra)
library(scales)
library(sbt)
library(readxl)
library(SparseNUTS)
theme_set(theme_bw())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)
model_file <- file.path(run_dir, "esc31_base.rds")
# One canonical sbt_fit owns the complete base run. A new MLE deliberately
# replaces any stale downstream state; a completed MCMC augments this file.
base_mcmc_file <- model_file
esc31_env_flag <- function(name) {
tolower(trimws(Sys.getenv(name, "false"))) %in%
c("1", "true", "yes", "on")
}
# Set TRUE only when replacing the saved MLE with a fresh optimization.
# This flag is independent of Quarto's execution cache.
run_new_mle <- FALSE || esc31_env_flag("ESC31_RUN_BASE_MLE")
# Set TRUE only when replacing the saved M-parameter profiles.
run_new_m_profiles <-
FALSE || esc31_env_flag("ESC31_RUN_BASE_M_PROFILES")
m_profile_file <- file.path(supporting_run_dir, "esc31_base_m_profiles.rds")
render_m_profile_results <- run_new_m_profiles ||
esc31_env_flag("ESC31_RENDER_BASE_M_PROFILES") ||
file.exists(m_profile_file)
# Set TRUE only when replacing the saved posterior with a fresh MCMC run.
run_new_mcmc <- FALSE || esc31_env_flag("ESC31_RUN_BASE_MCMC")
render_mcmc_results <- run_new_mcmc ||
esc31_env_flag("ESC31_RENDER_BASE_MCMC") ||
file.exists(base_mcmc_file)
base_loo_file <- file.path(
supporting_run_dir,
"esc31_base_loo.rds"
)
# Set TRUE only when replacing the saved base PSIS-LOO calculation.
run_new_base_loo <- FALSE || esc31_env_flag("ESC31_RUN_BASE_LOO")
render_base_loo_results <- run_new_base_loo ||
esc31_env_flag("ESC31_RENDER_BASE_LOO") ||
file.exists(base_loo_file)
base_loo_composition_pattern_text <-
"No corrected base LOO result was available."
base_loo_n_pattern_text <- ""
model_code_signature <- "sbt_model_base_2026_07_21"The previous-assessment comparison is generated from an explicit V1 MLE rather than the bundled ADMB output. The data, switches, parameterization, parameter map, and optimizer bounds match the sbt_vs_admb package vignette. The only specification change is that steepness is fixed at the current value h = 0.70 instead of the legacy grid-cell value h = 0.72.
previous_v1_model_file <- file.path(
supporting_run_dir,
"esc31_previous_assessment_v1_h070.sbt.rds"
)
run_new_previous_v1_mle <- esc31_env_flag("ESC31_RUN_PREVIOUS_V1_MLE")
previous_v1_release <- sbtdata::sbt_load_data("2023")
previous_v1_csv <- previous_v1_release$data_csv1
previous_v1_labrep <- previous_v1_release$data_labrep1
previous_v1_par <- previous_v1_release$data_par1
previous_v1_data_in <- list(
last_yr = 2022,
age_increase_M = 25,
length_m50 = 150,
length_m95 = 180,
catch_surf_case = 1,
catch_LL1_case = 1,
length_mean = previous_v1_release$length_mean,
length_sd = previous_v1_release$length_sd,
catch = previous_v1_release$catch,
catch_UA = previous_v1_release$catch_UA,
scenarios_surf = previous_v1_csv[["scenarios_surface"]],
scenarios_LL1 = previous_v1_csv[["scenarios_LL1"]],
POPs_v1 = previous_v1_release$POPs_v1,
HSPs = previous_v1_release$HSPs,
GTs = previous_v1_release$GTs,
aerial_survey = previous_v1_release$aerial_survey,
aerial_cov = previous_v1_release$aerial_cov,
troll = previous_v1_release$troll,
cpue = previous_v1_release$cpue,
age_freq = previous_v1_release$age_freq,
length_freq = previous_v1_release$length_freq,
tag_reporting = previous_v1_release$tag_reporting,
tag_releases = previous_v1_release$tag_releases,
tag_recaptures = previous_v1_release$tag_recaptures,
sel_min_age_f = c(2, 2, 2, 8, 6, 0),
sel_max_age_f = c(17, 9, 17, 22, 25, 7),
sel_end_f = c(1, 0, 1, 1, 1, 0),
sel_change_sd_fy = t(as.matrix(
previous_v1_csv[["sel_change_sd"]][, -1]
)),
sel_smooth_sd_f = previous_v1_labrep$sel.smooth.sd,
pop_switch = 1,
hsp_switch = 1,
hsp_false_negative = 0.7467647,
gt_switch = 1,
cpue_switch = 1,
cpue_a1 = 5,
cpue_a2 = 17,
aerial_switch = 4,
aerial_tau = previous_v1_labrep$tau.aerial,
troll_switch = 0,
lf_minbin = c(1, 1, 1, 11),
tag_switch = 1,
tag_var_factor = 1.82
)
previous_v1_data <- sbt::get_data_v1(data_in = previous_v1_data_in)
# The V1 comparison has no explicit parameter priors. An empty list preserves
# that historical target while allowing the staged fit to distinguish it from
# unresolved current-model defaults.
previous_v1_data$priors <- list()
previous_v1_parameters <- list(
par_log_B0 = previous_v1_par$ln_B0,
par_log_psi = log(previous_v1_par$psi),
par_log_m0 = log(previous_v1_par$m0),
par_log_m4 = log(previous_v1_par$m4),
par_log_m10 = log(previous_v1_par$m10),
par_log_m30 = log(previous_v1_par$m30),
par_log_h = log(0.70),
par_log_sigma_r = log(previous_v1_labrep$sigma.r),
par_log_cpue_q = previous_v1_par$lnq,
par_log_cpue_sigma = log(previous_v1_par$sigma_cpue),
par_log_cpue_omega = log(previous_v1_par$cpue_omega),
par_log_aerial_tau = log(previous_v1_par$tau_aerial),
par_log_aerial_sel = previous_v1_par$ln_sel_aerial,
par_log_troll_tau = log(previous_v1_par$tau_troll),
par_log_hsp_q = previous_v1_par$lnqhsp,
par_logit_hstar_i = qlogis(exp(
previous_v1_par$par_log_hstar_i
)),
par_log_tag_H_factor = log(previous_v1_par$tag_H_factor),
par_rdev_y = previous_v1_par$Reps,
par_sels_init_i = previous_v1_par$par_sels_init_i,
par_sels_change_i = previous_v1_par$par_sels_change_i
)
previous_v1_map <- list(
par_log_psi = factor(NA),
par_log_m0 = factor(NA),
par_log_m10 = factor(NA),
par_log_h = factor(NA),
par_log_sigma_r = factor(NA),
par_log_cpue_sigma = factor(NA),
par_log_cpue_omega = factor(NA),
par_log_troll_tau = factor(NA),
par_log_aerial_tau = factor(NA),
par_log_aerial_sel = factor(rep(NA, 2)),
par_log_hsp_q = factor(NA),
par_log_tag_H_factor = factor(NA)
)
previous_v1_get_bounds <- function(obj, parameters) {
lower <- rep(-Inf, length(obj$par))
upper <- rep(Inf, length(obj$par))
lower[grep("par_log_m0", names(obj$par))] <- log(0.2)
upper[grep("par_log_m0", names(obj$par))] <- log(0.55)
lower[grep("par_log_m4", names(obj$par))] <-
parameters$par_log_m10
upper[grep("par_log_m4", names(obj$par))] <- log(
0.333 * exp(parameters$par_log_m10) +
0.667 * exp(parameters$par_log_m0)
)
lower[grep("par_log_m10", names(obj$par))] <- log(0.029)
upper[grep("par_log_m10", names(obj$par))] <- log(0.21)
lower[grep("par_log_m30", names(obj$par))] <- log(0.2)
upper[grep("par_log_m30", names(obj$par))] <- log(0.7)
upper[grep("par_log_troll_tau", names(obj$par))] <- log(0.4)
lower[grep("par_rdev_y", names(obj$par))] <- -5
upper[grep("par_rdev_y", names(obj$par))] <- 5
names(lower) <- names(obj$par)
names(upper) <- names(obj$par)
list(lower = lower, upper = upper)
}if (run_new_previous_v1_mle || !file.exists(previous_v1_model_file)) {
previous_v1_control <- list(eval.max = 10000, iter.max = 10000)
previous_v1_fit <- sbt::sbt_fit(
data = previous_v1_data,
control = previous_v1_control,
metadata = list(
comparison = "base_page_previous_assessment_v1",
source = "sbt_vs_admb vignette specification",
fixed_h = 0.70
),
model = "sbt_model_v1"
)
previous_v1_fit <- sbt_add_parameters(
previous_v1_fit,
previous_v1_parameters
)
previous_v1_fit <- sbt_add_map(previous_v1_fit, previous_v1_map)
previous_v1_fit <- sbt_build_object(previous_v1_fit)
previous_v1_fit <- sbt_add_bounds(
previous_v1_fit,
previous_v1_get_bounds(
sbt_obj(previous_v1_fit),
previous_v1_fit$parameters
)
)
previous_v1_fit <- sbt_optimise(
previous_v1_fit,
n_passes = 3L,
check = TRUE,
check_args = list(gradient_tolerance = 0.01, cores = 1L)
)
sbt_fit_validation(
previous_v1_fit,
scope = "mle",
require_pass = TRUE
)
sbt_fit_save(previous_v1_fit, previous_v1_model_file, overwrite = TRUE)
} else {
previous_v1_fit <- sbt_fit_read(
previous_v1_model_file,
strict = TRUE,
rebuild = FALSE,
verify_validation = FALSE
)
}
previous_v1_report <- sbt::sbt_fit_report(previous_v1_fit)
if (!isTRUE(all.equal(
as.numeric(previous_v1_report$par_h)[1],
0.70,
tolerance = 1e-10
))) {
stop("The cached previous-assessment V1 fit does not use h = 0.70.")
}kable(
data.frame(
Objective = previous_v1_fit$fit$opt$objective,
`Maximum gradient` =
previous_v1_fit$fit$diagnostics$max_gradient,
Estimability = previous_v1_fit$fit$estimability$status,
`Fixed psi` = as.numeric(previous_v1_report$par_psi)[1],
`Fixed h` = as.numeric(previous_v1_report$par_h)[1],
check.names = FALSE
),
digits = 8
) |>
kable_styling(full_width = FALSE)| Objective | Maximum gradient | Estimability | Fixed psi | Fixed h |
|---|---|---|---|---|
| 6050.071 | 1.4e-07 | estimable | 1.75 | 0.7 |
The V1 comparison uses the same parameter treatment as the selected previous-assessment ADMB grid cell. Structural grid parameters and externally conditioned observation scales are fixed, while population scale, the two within-cell mortality anchors, CPUE catchability, annual recruitment deviations, and selectivity effects are estimated. The base-page comparison differs from the package vignette only by fixing h = 0.70 rather than the legacy cell’s h = 0.72.
previous_v1_hstar_count <- sum(
names(previous_v1_fit$fit$opt$par) == "par_logit_hstar_i"
)
if (!is.null(previous_v1_fit$map$par_logit_hstar_i) ||
previous_v1_hstar_count !=
length(previous_v1_fit$parameters$par_logit_hstar_i)) {
stop("The previous-assessment hstar parameter treatment is inconsistent.")
}
kable(
data.frame(
Treatment = c("Fixed", "Estimated"),
Parameters = c(
paste(
"psi, M0, M10, h, sigma_r, CPUE sigma and omega,",
"troll tau, aerial tau and selectivity, HSP q,",
"and the tag harvest-rate factor"
),
paste(
"B0, M4, M30, CPUE q,", previous_v1_hstar_count,
"hstar effects, annual recruitment",
"deviations, and initial and time-varying selectivity effects"
)
),
check.names = FALSE
)
) |>
kable_styling(full_width = FALSE)| Treatment | Parameters |
|---|---|
| Fixed | psi, M0, M10, h, sigma_r, CPUE sigma and omega, troll tau, aerial tau and selectivity, HSP q, and the tag harvest-rate factor |
| Estimated | B0, M4, M30, CPUE q, 17 hstar effects, annual recruitment deviations, and initial and time-varying selectivity effects |
The hstar treatment was checked against both the executable V1 map and the archived ADMB declaration. The hstar vector is absent from the fixed-parameter map, all 17 par_logit_hstar_i coefficients are active in the fitted parameter vector, and lnhstar_s1_ya is declared in phase 1 in the ADMB template. The hstar effects were therefore estimated; the separate tag harvest-rate factor was fixed.
The complete data setup used by the production workflow is shown below. It loads the checksum-verified 2026 release from sbtdata, defines the model switches, and records the scientific contract. The original direct-file reads remain commented for audit while the submitted files stay in csv_2026. The same straight-line code remains in esc31_inputs.R for the sensitivity, grid, projection, and command-line workflows.
source(file.path(esc_dir, "esc31_workflow.R"))
source(file.path(esc_dir, "esc31_base_posterior.R"))assessment_data <- sbtdata::sbt_load_data("2026")
length_mean <- assessment_data$length_mean
length_sd <- assessment_data$length_sd
catch <- assessment_data$catch
catch_UA <- assessment_data$catch_UA
scenarios_surface <- assessment_data$scenarios_surface
scenarios_LL1 <- assessment_data$scenarios_LL1
POPs <- assessment_data$POPs
HSPs <- assessment_data$HSPs
GTs <- assessment_data$GTs
troll <- assessment_data$troll
cpue <- assessment_data$cpue
age_freq <- assessment_data$age_freq
length_freq <- assessment_data$length_freq
cpue_gam22_cv <- assessment_data$cpue_gam22_cv
# Previous direct-file reads are retained temporarily for audit and comparison.
# data_loc <- esc31_data_dir(esc_dir)
# length_mean <- read_csv(file.path(data_loc, "mean_length.csv"))
# length_sd <- read_csv(file.path(data_loc, "sd_length.csv"))
# catch <- read_csv(file.path(data_loc, "catch.csv"))
# catch_UA <- read_csv(file.path(data_loc, "catch_UA.csv"))
# scenarios_surface <- read_csv(file.path(data_loc, "scenarios_surface.csv"))
# scenarios_LL1 <- read_csv(file.path(data_loc, "scenarios_LL1.csv"))
# POPs <- read_csv(file.path(data_loc, "POPs.csv"))
# HSPs <- read_csv(file.path(data_loc, "HSPs.csv"))
# GTs <- read_csv(file.path(data_loc, "GTs.csv"))
# troll <- read_csv(file.path(data_loc, "trolling_index.csv"))
# cpue <- read_csv(file.path(data_loc, "cpue.csv"))
# age_freq <- read_csv(file.path(data_loc, "age_freq.csv"))
# length_freq <- read_csv(file.path(data_loc, "lf_assessment.csv"))
# cpue_gam22_cv <- read_excel(file.path(data_loc, "CV_GAM22_20260713.xlsx"), sheet = "Index", range = "E3:H60")
# names(cpue_gam22_cv) <- c("Year", "index_mean", "cv_raw", "cv_scaled")
# cpue_gam22_cv$Year <- as.integer(cpue_gam22_cv$Year)
# cpue_gam22_cv$CV <- cpue_gam22_cv$cv_raw
# cpue$CV <- cpue_gam22_cv$CV[match(cpue$Year, cpue_gam22_cv$Year)]data_in <- list(
last_yr = 2025,
M_switch = 2L, age_increase_M = 25,
length_m50 = 150, length_m95 = 180,
catch_surf_case = 1L,
catch_LL1_case = 1L,
length_mean = length_mean, length_sd = length_sd,
catch = catch,
catch_UA = catch_UA,
scenarios_surf = scenarios_surface,
scenarios_LL1 = scenarios_LL1,
paly = assessment_data$paly,
POPs = POPs,
HSPs = HSPs,
GTs = GTs,
troll = troll,
cpue = cpue,
age_freq = age_freq,
length_freq = length_freq,
aerial_survey = assessment_data$aerial_survey,
aerial_cov = assessment_data$aerial_cov,
tag_reporting = assessment_data$tag_reporting,
tag_releases = assessment_data$tag_releases,
tag_recaptures = assessment_data$tag_recaptures,
harvest_wall_strength = 10,
harvest_wall_onset = 0.85,
harvest_wall_ceiling = 0.9,
harvest_wall_scale = 0.01,
removal_switch_f = c(0, 0, 0, 0, 0, 0),
sel_min_age_f = c(2, 2, 2, 8, 6, 0, 4),
sel_max_age_f = c(17, 9, 17, 21, 25, 7, 17),
sel_end_f = c(1, 0, 1, 1, 1, 0, 1),
sel_LL1_yrs = c(1952, seq(1957, 2001, 4), 2006:2008, seq(2011, 2023, 3)),
sel_LL2_yrs = c(1969, 2001, 2005, 2008, seq(2011, 2023, 3)),
sel_LL3_yrs = c(1954, 1961, 1965, 1969:1971, 2005:2007),
sel_LL4_yrs = 1953,
sel_Ind_yrs = c(1976, 1997, 1999, seq(2002, 2010, 2), 2012:2021),
sel_Aus_yrs = c(1952, seq(1969, 1993, 4), 1997:2025),
sel_CPUE_yrs = c(seq(1969, 2001, 4), 2006:2008, seq(2011, 2023, 3)),
af_switch = 1L,
lf_switch = 1L, lf_minbin = c(1, 1, 1, 11, 6),
cpue_switch = 1L, cpue_a1 = 5, cpue_a2 = 17, cpue_sel_fishery = 7L,
cpue_lf_switch = 1L, cpue_lf_sel_fishery = 7L,
aerial_switch = 4L,
troll_switch = 0L,
pop_switch = 1L,
hsp_switch = 1L, hsp_false_negative = 0.6840729,
gt_switch = 1L,
tag_switch = 1L, tag_var_factor = 2.4
)The base conditioning includes the NCNM/UAM additions stored in the 2026 sbtdata release (sourced from catch_UA.csv). The NoUAM sensitivity removes those additions before rebuilding the data. This is separate from the retired ADMB blanket catch-underreporting switch, which is not part of the current sbt or ESC31 input interface.
base_cpue_cv_spec <- list(
implementation = "gam22_raw_plus_default_cpue_sigma_v1",
file = "CV_GAM22_20260713.xlsx",
sheet = "Index",
block = "B=1000",
cv_column = "cv_raw",
combination = "sigma2_total = default_cpue_log_sigma^2 + log1p(cv_raw^2)",
decision_id = "esc31_2026_cpue_sigma_020_plus_gam22_raw_v1"
)
workbook_md5 <- unname(tools::md5sum(
file.path(esc31_data_dir(esc_dir), "CV_GAM22_20260713.xlsx")
))
ll34_conditioning_spec <- list(
status = "approved",
approved_configuration =
"ll3_ll4_separate_standard_selectivity_single_ll4_block_ll3_hyperparameters",
ll3_treatment =
"standard_selectivity_based_removal_with_soft_lf_likelihood",
ll4_treatment =
"standard_selectivity_based_removal_with_soft_lf_likelihood",
ll3_ll4_lumping = FALSE,
catch_accounting =
"retain_full_observed_ll4_catch_by_standard_conditioning",
ll4_selectivity = list(
n_time_blocks = 1L,
block_start_years = 1953L,
estimated_ages = 8:21,
n_length_compositions = 38L,
terminal_extension = "age_21_value_extended_through_terminal_age",
parameterization = "14_mean_normalized_log_age_effects",
parameter_map = "all_14_ll4_age_effects_estimated",
prior = list(
implementation = "fixed_separable_ar1",
inheritance = "copy_ll3_fixed_hyperparameters",
source_fishery = "LL3",
target_fishery = "LL4",
rho_year = 0.5,
rho_age = 0.5,
sigma = 0.75,
review_status = "accepted_after_mle_shape_review"
)
),
ll4_length_composition_likelihood = "multinomial_soft_likelihood",
combined_age_specific_harvest_limit = 0.9,
accepted_continuation_penalty = 0,
positive_abundance_required = TRUE,
legacy_admb_configuration = "ll3_ll4_selectivity_soft_lf_likelihood",
rejected_configuration = "lumped_ll3_ll4_exact_sliced_removal",
superseded_base_configuration =
"ll3_standard_ll4_separate_full_exact_sliced_removal",
direct_sliced_ll4_role = "formal_sensitivity_pending_implementation",
workshop_reference = "OM WS 2023 Tokyo report, agenda item 5, paragraphs 19 and 27(a)",
decision_id =
"esc31_2026_ll3_ll4_separate_standard_single_block_ll3_hyper_v2"
)
harvest_wall_spec <- list(
status = "approved",
implementation = "normalized_squared_softplus_harvest_wall_v1",
strength = 10,
onset = 0.85,
ceiling = 0.9,
scale = 0.01,
aggregation = "sum_over_raw_combined_year_season_age_harvest",
role = "preventive_objective_regularization",
scope = "global_model_fit_and_posterior",
continuation_accounting = "separate_unchanged_existing_posfun",
catch_accounting = "unchanged_exact_conditioned_catch",
decision_id = "esc31_2026_preventive_harvest_wall_a10_v1"
)
natural_mortality_spec <- list(
status = "approved",
M_switch = 2L,
implementation = "m10_scaled_length_allometry_v1",
curve = paste(
"M[a] = M10 * (L[a] / L10)^mc through age_increase_M;",
"linear on the natural scale to M30 at max_age"
),
length_reference = "first_year_first_season_mean_length_at_age",
age_increase_M = 25L,
parameterization = c(
M10 = "exp(par_log_m10)",
M30 = "exp(par_log_m30)",
mc = "par_mc",
M0 = "derived_from_M10_length_and_fixed_mc",
M4 = "derived_from_M10_length_and_fixed_mc"
),
estimated_parameters = c("par_log_m10", "par_log_m30"),
fixed_parameters = c(mc = -1),
parameter_map = c(par_mc = "fixed"),
optimizer_bounds = c(M10 = "[0.029, 0.21]", M30 = "[0.20, 0.70]"),
exponent_convention = "fixed_mc_minus1_in_direct_length_ratio",
prior_scale = "direct_log_m10_m30",
change_of_variables_adjustment = "none",
posterior_mortality_requirement =
"finite_positive_at_every_age_for_every_draw_no_exceptions",
supersedes = c(
"esc31_2026_direct_log_m4_v1",
"esc31_2026_two_age2_mortality_draws_v1"
),
decision_id = "esc31_2026_m10_scaled_length_m_mc_minus1_v1"
)
observation_error_spec <- list(
status = "approved",
implementation = "fixed_aerial_tag_sdnr_calibration_v1",
aerial = list(
parameter = "par_log_aerial_tau",
additional_log_scale_sd = 0.59,
treatment = "fixed_by_get_map"
),
conventional_tags = list(
data_field = "tag_var_factor",
variance_factor = 2.4,
treatment = "fixed_model_data"
),
residual_diagnostic = list(
target = "SDNR approximately one",
conventional_tag_seed = 123L
),
role = "downweight_aerial_and_conventional_tag_likelihoods",
decision_id = "esc31_2026_aerial_tag_sdnr_calibration_v1"
)
review_method_decisions <- list(
troll = list(
status = "accepted",
decision_id = "esc31_2026_troll_base_aerial_tau_059_v1",
interpretation = "trolling_index_only",
aerial_tau = 0.59,
future_larger_tau_requires_new_sensitivity = TRUE
),
conventional_tag_osa = list(
status = "accepted",
decision_id =
"esc31_2026_conventional_tag_compResidual_osa_fallback_v1",
method = "compResidual::resDirM",
seed = 123L,
role = "diagnostic_only_no_fit_change"
)
)
base_model_contract <- list(
data_in = data_in,
cpue_cv_input = cpue_gam22_cv,
audit_metadata = list(workbook_md5 = workbook_md5),
scientific_config = c(
base_cpue_cv_spec,
list(
file_md5 = workbook_md5,
aligned_cv_md5 = esc31_object_md5(cpue_gam22_cv),
years = cpue_gam22_cv$Year,
ll34_conditioning = ll34_conditioning_spec,
harvest_wall = harvest_wall_spec,
natural_mortality = natural_mortality_spec,
observation_error = observation_error_spec
)
)
)base_fit <- sbt_fit(
data_in = data_in,
metadata = list(
assessment = "ESC31",
model_code_signature = model_code_signature,
raw_input_md5 = base_run_identity$raw_input_md5,
workflow_git_sha = base_run_identity$workflow_git_sha,
sbt_git_sha = base_run_identity$sbt_git_sha,
run_identity = base_run_identity
)
)These are all external and package data objects supplied through sbt_fit(data_in = data_in) to get_data(). Dimensions are shown before the package converts them into model arrays. Table 3 is therefore the source-and-dimensions inventory for the assessment inputs, rather than a model-results table.
input_tables <- list(
"mean_length.csv" = length_mean,
"sd_length.csv" = length_sd,
"catch.csv" = catch,
"catch_UA.csv" = catch_UA,
"scenarios_surface.csv" = scenarios_surface,
"scenarios_LL1.csv" = scenarios_LL1,
"POPs.csv" = POPs,
"HSPs.csv" = HSPs,
"GTs.csv" = GTs,
"trolling_index.csv" = troll,
"age_freq.csv" = age_freq,
"lf_assessment.csv" = length_freq,
"cpue.csv + CV_GAM22_20260713.xlsx" = cpue,
"paly" = assessment_data$paly,
"aerial_survey" = assessment_data$aerial_survey,
"aerial_cov" = assessment_data$aerial_cov,
"tag_reporting" = assessment_data$tag_reporting,
"tag_releases" = assessment_data$tag_releases,
"tag_recaptures" = assessment_data$tag_recaptures
)
input_summary <- tibble(
Input = names(input_tables),
Source = rep("sbtdata 2026 release", length(input_tables)),
Dimensions = map_chr(input_tables, ~ paste(dim(.x), collapse = " x "))
) |>
kable() |>
kable_styling(full_width = FALSE)The 2026 sbtdata release retains the submitted GAM22 index values from cpue.csv and the raw year-specific GAM estimation uncertainty from the validated area-weighted B=1000 workbook. That uncertainty is added to the historical 20% CPUE uncertainty in log-variance space. The simulation-based index_scaled and cv_scaled values are not used by the model.
cpue_default_log_sigma <- exp(
get_parameters(data = base_fit$data)$par_log_cpue_sigma
)
cpue_uncertainty <- cpue |>
mutate(
raw_gam_log_sd = sqrt(log1p(CV^2)),
combined_log_sd = sqrt(raw_gam_log_sd^2 + cpue_default_log_sigma^2),
combined_cv = sqrt(expm1(combined_log_sd^2))
)
cpue_uncertainty |>
filter(Year >= max(Year) - 9L) |>
transmute(
Year,
`Submitted GAM22 index` = CPUE,
`Raw GAM22 CV` = CV,
`Combined CPUE CV used by model` = combined_cv,
`Combined log-scale SD used by model` = combined_log_sd
) |>
kable(digits = 3, caption = "Recent GAM22 index and observation uncertainty") |>
kable_styling(full_width = FALSE)| Year | Submitted GAM22 index | Raw GAM22 CV | Combined CPUE CV used by model | Combined log-scale SD used by model |
|---|---|---|---|---|
| 2016 | 1.247 | 0.097 | 0.225 | 0.222 |
| 2017 | 1.256 | 0.094 | 0.224 | 0.221 |
| 2018 | 1.393 | 0.089 | 0.222 | 0.219 |
| 2019 | 1.633 | 0.094 | 0.224 | 0.221 |
| 2020 | 1.272 | 0.093 | 0.223 | 0.220 |
| 2021 | 1.315 | 0.099 | 0.226 | 0.223 |
| 2022 | 2.151 | 0.104 | 0.228 | 0.225 |
| 2023 | 1.723 | 0.119 | 0.236 | 0.233 |
| 2024 | 2.009 | 0.153 | 0.255 | 0.251 |
| 2025 | 2.761 | 0.219 | 0.301 | 0.294 |
The raw GAM22 CVs range from 1.28% to 21.87%. After combining their log variances with the package-default time-invariant CPUE log-scale sigma of 0.20, the ordinary-scale CVs used by the model range from 20.24% to 30.10%, with mean 21.34%.
The following figures show a representative spread of the biological, fishery, index, monitoring, and composition inputs before fitting. They are input summaries rather than goodness-of-fit plots: model expectations and posterior intervals are introduced later under Model Fit Plots. The composition panels use the three most recent available years within each fishery so that every active series is represented.
The inputs continue the CCSBT data-exchange and catch-review series (Commission for the Conservation of Southern Bluefin Tuna 2026; CCSBT Secretariat 2023a, 2023b). The NCNM additions are documented by Edwards and Hoyle (2023); the Indonesian length/age and Australian otolith inputs by Indonesia (2023) and Australia (2023); and the trolling index by Itoh (2023). The conventional-tag treatment follows the established recovery and reporting/shedding framework (Brownie et al. 1985; Hampton and Fournier 1997; Eveson et al. 2015), while the close-kin inputs build on the SBT CKMR design and its current sampling update (Davies et al. 2012; Farley et al. 2023).
plot_length_at_age(
data = base_fit$data,
years = input_growth_plot_years
)
plot_weight_at_age(
data = base_fit$data,
years = input_growth_plot_years
)
plot_paly(data = base_fit$data)
ggplot(input_catch_plot_data, aes(x = Year, y = Catch)) +
geom_line(color = input_data_color, linewidth = 0.45) +
geom_point(color = input_data_color, size = 0.9) +
facet_wrap(vars(Fishery), scales = "free_y", ncol = 2) +
scale_x_continuous(
limits = input_year_limits(input_catch_plot_data$Year),
breaks = pretty_breaks()
) +
scale_y_continuous(
labels = comma,
limits = c(0, NA),
expand = expansion(mult = c(0, 0.05))
) +
labs(x = "Year", y = "Catch (tonnes)")
ggplot(input_index_plot_data, aes(x = Year, y = Value)) +
geom_line(color = input_data_color, linewidth = 0.45) +
geom_point(color = input_data_color, size = 1.2) +
facet_wrap(vars(Series), scales = "free_y", ncol = 2) +
scale_x_continuous(
limits = input_year_limits(input_index_plot_data$Year),
breaks = pretty_breaks()
) +
scale_y_continuous(
labels = comma,
limits = c(0, NA),
expand = expansion(mult = c(0, 0.05))
) +
labs(x = "Year", y = "Input value")
plot_tags_recaptures_pooled(base_fit)
ggplot(
input_length_snapshot,
aes(x = Length, y = Proportion, color = Year)
) +
geom_line(linewidth = 0.5) +
geom_point(size = 0.8) +
facet_wrap(vars(Fishery), scales = "free_y", ncol = 2) +
scale_x_continuous(breaks = pretty_breaks(n = 5)) +
scale_y_continuous(
limits = c(0, NA),
expand = expansion(mult = c(0, 0.05))
) +
labs(x = "Length bin (cm)", y = "Proportion", color = "Year")
ggplot(
input_age_snapshot,
aes(x = Age, y = Proportion, color = Year)
) +
geom_line(linewidth = 0.6) +
geom_point(size = 1.3) +
facet_wrap(vars(Fishery), scales = "free_y", ncol = 1) +
scale_x_continuous(breaks = pretty_breaks(n = 7)) +
scale_y_continuous(
limits = c(0, NA),
expand = expansion(mult = c(0, 0.05))
) +
labs(x = "Age", y = "Proportion", color = "Year")
The base model is set up as an ESC31 implementation of the previous CCSBT stock-assessment structure. The main reference case is the 2023 SBT assessment (Hillary et al. 2023), with supporting inputs from the published GAM22 CPUE update through 2024 (Itoh and Takahashi 2025), gene-tagging update (Preece and Bradford 2023), CKMR sampling and kin-finding update (Farley et al. 2023), and fisheries-indicator summaries (Patterson and Woodhams 2023).
The current annual GAM22 index and B=1000 estimation CVs through 2025 come from the supplied CV_GAM22_20260713.xlsx workbook in the raw-input manifest. The 2025 endpoint is therefore sourced to that dated ESC31 workbook rather than being implied to occur in the published through-2024 paper.
The base CPUE observation error retains the package-default time-invariant log-scale sigma and adds the raw year-specific GAM22 estimation CV inside the CPUE likelihood. The workbook CV is converted to log variance, giving total log-scale variance default_cpue_log_sigma^2 + log1p(cv_raw^2). The resulting lower CPUE OSA residual SDNR is accepted.
Natural mortality estimates \(M_{10}\) and \(M_{30}\) on the log scale. From age 0 through age 25, \(M_a=M_{10}(L_a/L_{10})^{m_c}\) using the fixed first-year, first-season mean length-at-age curve. Mortality then increases linearly on the natural scale to \(M_{30}\) at age 30. The exponent \(m_c\) is an explicit model parameter fixed at -1; \(M_0\) and \(M_4\) are derived values. Every retained posterior draw must have finite, positive mortality at every age; there is no mortality exception.
Fixing \(m_c=-1\) makes mortality inversely proportional to length through age 25. This length-inverse form follows the empirical body-size scaling described by Lorenzen (1996), the explicit inverse-length model evaluated by Lorenzen (2000), and the generalized length-inverse mortality synthesis of Lorenzen (2022). In this assessment, \(m_c\) is a fixed structural assumption, whereas \(M_{10}\) and \(M_{30}\) are estimated.
base_fit <- sbt_add_parameters(base_fit)
base_fit <- sbt_add_parameters(
base_fit,
overrides = list(
par_log_psi = log(base_psi),
par_log_h = log(base_h),
par_log_B0 =
base_fit$parameters$par_log_B0 + log(base_b0_start_multiplier)
)
)base_fit <- sbt_add_map(base_fit)The explicit parameter priors below are the priors returned by get_priors() for the ESC31 base model and printed with priors_to_math(). Recruitment-deviation and selectivity priors are separate structured prior terms in the model objective.
base_fit <- sbt_add_priors(base_fit)
cat(priors_to_math(base_fit$data$priors), sep = "\n\n")par_log_psi: \(\log(\psi) \sim \mathrm{Normal}\left(\log(1.75),\,0.122^2\right)\).
par_log_m10: \(\log(M_{10}) \sim \mathrm{Normal}\left(\log(0.1),\,0.6^2\right)\).
par_log_m30: \(\log(M_{30}) \sim \mathrm{Normal}\left(\log(0.4574),\,1.5^2\right)\).
par_mc: \(m_c \sim \mathrm{Normal}\left(-1,\,0.3^2\right)\).
par_log_h: \(\log(h) \sim \mathrm{Normal}\left(\log(0.55),\,0.5455^2\right)\).
par_log_cpue_omega: \(\log(\omega_{\mathrm{CPUE}}) \sim \mathrm{Normal}\left(\log(0.875),\,0.1143^2\right)\).
par_log_sigma_r: \(\log(\sigma_R) \sim \mathrm{Normal}\left(\log(0.6),\,1^2\right)\).
par_cpue_creep: \(c_{\mathrm{CPUE}} \sim \mathrm{Normal}\left(0.005,\,0.01^2\right)\).
par_sel_rho_y: \(\mathrm{logit}\left((\rho_{\mathrm{sel},y}+1)/2\right) \sim \mathrm{Normal}\left(3.664,\,1^2\right)\).
par_sel_rho_a: \(\mathrm{logit}\left((\rho_{\mathrm{sel},a}+1)/2\right) \sim \mathrm{Normal}\left(3.664,\,1^2\right)\).
par_log_sel_sigma: \(\log(\sigma_{\mathrm{sel}}) \sim \mathrm{Normal}\left(\log(0.5),\,1^2\right)\).
The canonical base MLE uses the staged fitted-model lifecycle. Objective construction and bounds are resolved automatically, and the two sequential optimizer passes start from RTMB’s retained best state. A failed diagnostic gate leaves the expensive result on disk for inspection but cannot pass the acceptance check below.
base_fit <- sbt_optimise(
base_fit,
n_passes = 2L,
control = list(eval.max = 10000, iter.max = 1000),
check = TRUE,
check_args = list(
gradient_tolerance = 0.01
)
)
sbt_fit_save(base_fit, model_file, overwrite = TRUE)
sbt_fit_validation(base_fit, scope = "mle", require_pass = TRUE)The complete portable base run is saved to runs/esc31_base.rds. The same sbt_fit contains the MLE, MCMC draws, diagnostics, provenance, and embedded posterior plot summaries; there is no second base-posterior fit.
base_harvest_wall_identity <- esc31_fit_harvest_wall_identity(base_fit)
state_summary <- base_fit$fit$diagnostics$biological_state_mle$summary[1L, ]
state_draw <- base_fit$fit$diagnostics$biological_state_mle$per_draw[1L, ]
ll4_hyperparameters <- c(
rho_year = sel_rho_from_par(base_fit$parameters$par_sel_rho_y)[[4L]],
rho_age = sel_rho_from_par(base_fit$parameters$par_sel_rho_a)[[4L]],
sigma = exp(base_fit$parameters$par_log_sel_sigma[[4L]])
)
mle_acceptance_table <- tibble(
Check = c(
"Final negative log-likelihood",
"Optimizer convergence code",
"Maximum absolute gradient",
"Estimability",
"Biological-state gate",
"Maximum raw seasonal harvest at age",
"Minimum populated numbers-at-age",
"Maximum all-fleet scaled catch error",
"Preventive harvest-wall settings",
"Harvest-wall objective contribution",
"Continuation objective contribution",
"LL4 standard-fleet contract",
"LL4 selectivity structure",
"LL4 rho-year / rho-age / sigma",
"LL4 length-composition evidence",
"Maximum scaled LL4 catch error",
"LL4 peak selectivity"
),
Requirement = c(
"Finite", "0", "<= 0.01", "Estimable", "Pass", "<= 0.9",
">= -1e-8", "<= 1e-8", "10 / 0.85 / 0.90 / 0.01", "Report",
"<= 1e-8", "Pass", "1 block; 1953; ages 8--21; 14 parameters",
"0.5 / 0.5 / 0.75", "Review", "<= 1e-8", "Review"
),
Result = c(
format_decimal(base_fit$fit$opt$objective, 3),
as.character(base_fit$fit$opt$convergence),
format_decimal(base_fit$fit$diagnostics$max_gradient, 6),
base_fit$fit$estimability$message,
ifelse(isTRUE(state_summary$passes), "Pass", "Fail"),
format_decimal(state_summary$max_raw_harvest_rate, 6),
format_decimal(state_summary$min_number, 3),
paste0(
formatC(state_draw$max_catch_relative_error, format = "e", digits = 3)
),
paste0(
format_decimal(base_harvest_wall_identity$executable$strength, 2), " / ",
format_decimal(base_harvest_wall_identity$executable$onset, 2), " / ",
format_decimal(base_harvest_wall_identity$executable$ceiling, 2), " / ",
format_decimal(base_harvest_wall_identity$executable$scale, 2)
),
format_decimal(
base_fit$fit$diagnostics$harvest_wall_mle$objective_contribution,
6
),
format_decimal(
base_fit$fit$diagnostics$harvest_wall_mle$
continuation_objective_contribution,
10
),
ifelse(
isTRUE(base_fit$fit$diagnostics$ll4_base_mle$passes),
"Pass",
"Fail"
),
paste0(
base_fit$fit$diagnostics$ll4_base_mle$n_time_blocks, " block from ",
base_fit$fit$diagnostics$ll4_base_mle$block_start_year, "; ages ",
base_fit$fit$diagnostics$ll4_base_mle$minimum_estimated_age, "--",
base_fit$fit$diagnostics$ll4_base_mle$maximum_estimated_age, "; ",
base_fit$fit$diagnostics$ll4_base_mle$n_active_selectivity_parameters,
" parameters"
),
paste(format_decimal(ll4_hyperparameters, 2), collapse = " / "),
paste0(
base_fit$fit$diagnostics$ll4_base_mle$n_length_compositions,
" compositions; effective N ",
format_decimal(
base_fit$fit$diagnostics$ll4_base_mle$total_effective_sample_size,
3
),
"; NLL ",
format_decimal(
base_fit$fit$diagnostics$ll4_base_mle$length_composition_nll,
3
)
),
formatC(
base_fit$fit$diagnostics$ll4_base_mle$maximum_relative_catch_error,
format = "e",
digits = 3
),
paste0(
"Age ",
base_fit$fit$diagnostics$ll4_base_mle$maximum_selectivity_age,
"; value ",
format_decimal(
base_fit$fit$diagnostics$ll4_base_mle$maximum_selectivity,
3
)
)
)
)
mle_acceptance_table |> kable(align = c("l", "l", "l"))| Check | Requirement | Result |
|---|---|---|
| Final negative log-likelihood | Finite | 7,616.638 |
| Optimizer convergence code | 0 | 0 |
| Maximum absolute gradient | <= 0.01 | 0.000000 |
| Estimability | Estimable | All 1,559 active fixed-effect parameters are estimable. |
| Biological-state gate | Pass | Pass |
| Maximum raw seasonal harvest at age | <= 0.9 | 0.617533 |
| Minimum populated numbers-at-age | >= -1e-8 | 1,696.516 |
| Maximum all-fleet scaled catch error | <= 1e-8 | 2.533e-16 |
| Preventive harvest-wall settings | 10 / 0.85 / 0.90 / 0.01 | 10.00 / 0.85 / 0.90 / 0.01 |
| Harvest-wall objective contribution | Report | 0.000000 |
| Continuation objective contribution | <= 1e-8 | 0.0000000000 |
| LL4 standard-fleet contract | Pass | Pass |
| LL4 selectivity structure | 1 block; 1953; ages 8–21; 14 parameters | 1 block from 1953; ages 8–21; 14 parameters |
| LL4 rho-year / rho-age / sigma | 0.5 / 0.5 / 0.75 | 0.50 / 0.50 / 0.75 |
| LL4 length-composition evidence | Review | 38 compositions; effective N 334.995; NLL 275.878 |
| Maximum scaled LL4 catch error | <= 1e-8 | 2.058e-16 |
| LL4 peak selectivity | Review | Age 11; value 5.172 |
estimated_parameters <- make_parameter_table(
fit = base_fit,
exclude_recruitment_devs = TRUE,
include_selectivity_devs = FALSE,
include_details = TRUE
)
estimated_parameters_report <- estimated_parameters |>
filter(Estimated) |>
select(any_of(c("Parameter", "Model scale", "Value", "Lower", "Upper", "Prior", "Prior par1", "Prior par2"))) |>
rename(
Distribution = Prior,
par1 = `Prior par1`,
par2 = `Prior par2`
) |>
mutate(Parameter = parameter_labels(Parameter))
estimated_parameters_report |>
format_decimal_table(digits = 3) |>
replace_missing() |>
kable(
align = numeric_table_align(estimated_parameters_report),
caption = "Estimated non-recruitment, non-annual-selectivity model parameters. Log-scale parameters are shown on natural scale in the Value column."
) |>
add_header_above(c(" " = 5, "Prior" = 3))| Parameter | Model scale | Value | Lower | Upper | Distribution | par1 | par2 |
|---|---|---|---|---|---|---|---|
| B0 | 15.675 | 6,422,189.198 | 100,000.000 | 1,000,000,000.000 | - | - | - |
| CPUE q | -0.043 | 0.957 | 0.000 | 10,000.000 | - | - | - |
| M10 | -2.232 | 0.107 | 0.029 | 0.210 | normal | 0.100 | 0.600 |
| M30 | -0.782 | 0.458 | 0.200 | 0.700 | normal | 0.457 | 1.500 |
base_mcmc_sampler_table <- tibble(
Setting = c(
"Chains", "Warmup iterations per chain", "Retained iterations per chain",
"Total retained draws", "Metric", "Chain initialization",
"Target acceptance probability", "Maximum tree depth", "Seed",
"Optimization inside sampler"
),
Value = c(
4L,
150L,
750L,
3000L,
"Dense SparseNUTS inverse-Hessian metric",
"Exact source-MLE last.par.best",
0.999,
13L,
73001L,
"Skipped"
)
)
base_mcmc_sampler_table |> kable(align = c("l", "l"))| Setting | Value |
|---|---|
| Chains | 4 |
| Warmup iterations per chain | 150 |
| Retained iterations per chain | 750 |
| Total retained draws | 3000 |
| Metric | Dense SparseNUTS inverse-Hessian metric |
| Chain initialization | Exact source-MLE last.par.best |
| Target acceptance probability | 0.999 |
| Maximum tree depth | 13 |
| Seed | 73001 |
| Optimization inside sampler | Skipped |
The posterior uses four chains with 150 SparseNUTS warmup and 750 retained iterations per chain. Every chain starts at the exact MLE mode and the retained 3,000 draws must pass every production diagnostic gate.
if (run_new_mcmc) {
config <- base_page_mcmc_contract$sampler
base_fit <- sbt_mcmc(
base_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,
seed = config$seed,
skip_optimization = config$skip_optimization,
control = list(
adapt_delta = config$adapt_delta,
max_treedepth = config$max_treedepth
),
refresh = config$refresh
)
base_fit <- esc31_finalize_base_mcmc(
fit = base_fit,
contract = base_page_mcmc_contract,
file = base_mcmc_file
)
} else {
base_fit <- esc31_load_base_mcmc(
base_mcmc_file, base_fit, esc_dir
)$fit
}if (nrow(posterior_parameter_summary)) {
posterior_parameter_summary |>
format_decimal_table(digits = 4) |>
replace_missing() |>
kable(
align = c("l", rep("r", 4)),
caption = paste(
"Maximum-likelihood estimates and posterior summaries for six",
"reported scalar quantities. M0 and M4 are derived from M10, fixed",
"mc = -1, and length at age. Intervals are equal-tailed 95%",
"credible intervals from all",
format(length(base_retained_draws[, , 1L]), big.mark = ","),
"retained draws."
)
)
} else {
knitr::asis_output(
"*No saved base posterior was available for this render.*"
)
}| Parameter | MLE | Median | Lower 95% | Upper 95% |
|---|---|---|---|---|
| B0 | 6,422,189.1984 | 6,309,202.8164 | 5,109,888.7095 | 7,755,606.3267 |
| M0 | 0.3657 | 0.3674 | 0.3447 | 0.3913 |
| M4 | 0.1617 | 0.1624 | 0.1524 | 0.1730 |
| M10 | 0.1074 | 0.1079 | 0.1012 | 0.1149 |
| M30 | 0.4577 | 0.4638 | 0.3803 | 0.5626 |
| CPUE q | 0.9575 | 0.9612 | 0.9056 | 1.0192 |
The comparison below shows all four estimated scalar parameters other than recruitment deviations and selectivity parameters. M10 and M30 have explicit scalar priors, so both prior and posterior densities are shown. B0 and CPUE catchability have no explicit scalar prior in the base specification, so their panels show the posterior only. Fixed parameters are not shown. Recruitment deviations and selectivity surfaces use separate structured prior terms and are reported elsewhere. Posterior density lines are drawn only from the minimum to the maximum retained posterior draw. Prior densities are initially evaluated across their central 99.8% range, but each panel displays the prior only to one retained-posterior range beyond the posterior limits. This display-only clipping keeps a diffuse prior from determining the x-axis while leaving the prior used in fitting unchanged.
if (base_mcmc_accepted) {
base_prior_posterior_plot <- plot_prior_posterior(
base_fit,
pars = c(
"par_log_B0",
"par_log_m10",
"par_log_m30",
"par_log_cpue_q"
),
include_without_prior = TRUE,
labels = c(
par_log_B0 = "B0",
par_log_m10 = "M10",
par_log_m30 = "M30",
par_log_cpue_q = "CPUE catchability"
),
scale = "natural",
ncol = 2L
)
base_prior_posterior_limits <- base_prior_posterior_plot$data |>
filter(as.character(.data$Distribution) == "Posterior") |>
group_by(.data$Parameter) |>
summarise(
posterior_min = min(.data$Value),
posterior_max = max(.data$Value),
posterior_span = .data$posterior_max - .data$posterior_min,
display_min = .data$posterior_min - .data$posterior_span,
display_max = .data$posterior_max + .data$posterior_span,
.groups = "drop"
)
base_prior_posterior_plot$data <-
base_prior_posterior_plot$data |>
left_join(base_prior_posterior_limits, by = "Parameter") |>
filter(
as.character(.data$Distribution) != "Prior" |
(.data$Value >= .data$display_min &
.data$Value <= .data$display_max)
) |>
select(
-.data$posterior_min,
-.data$posterior_max,
-.data$posterior_span,
-.data$display_min,
-.data$display_max
)
base_prior_posterior_plot +
scale_y_continuous(
limits = c(0, NA),
expand = expansion(mult = c(0, 0.05))
)
}
Pareto-smoothed importance-sampling leave-one-out cross-validation (PSIS-LOO) approximates the predictive density for each omitted observation from the existing posterior draws (Vehtari et al. 2017). The Pareto shape diagnostic assesses the stability of that importance-sampling approximation (Vehtari et al. 2024). LOO-IC is -2 * elpd_loo, so smaller values indicate better expected out-of-sample predictive performance when models are evaluated against the same observations, effective weights, likelihood definitions, and pointwise likelihood partition. Because only the base model is evaluated here, its LOO-IC is a reference value rather than a model-ranking result.
Length-frequency rows with zero input effective N remain unchanged in the model data and prediction/output structures. They contribute no composition likelihood and are excluded only from the pointwise LOO matrix. For multinomial composition data, pointwise likelihoods use the same gamma-function expression as the fitted objective so that fractional effective bin counts are evaluated consistently. Any genuinely non-finite Pareto shape diagnostic remains in the numerical result but is omitted from plot_loo() rather than being assigned an artificial finite coordinate. The SE column in the LOO summary describes variation across pointwise predictive contributions; it is distinct from the Monte Carlo SE of the importance-sampling estimate discussed below.
base_active_parameter_count <- length(base_fit$mcmc$par_names)
if (base_active_parameter_count != length(base_fit$mcmc$mle$est)) {
stop(
"The posterior and MLE parameter counts do not agree.",
call. = FALSE
)
}
if (base_loo_available) {
base_loo_summary <- base_loo$estimates |>
as.data.frame() |>
rownames_to_column("Quantity") |>
as_tibble() |>
mutate(
Quantity = recode(
.data$Quantity,
elpd_loo = "Expected log predictive density (elpd_loo)",
p_loo = "Effective number of parameters (p_loo)",
looic = "LOO information criterion (LOO-IC)"
)
) |>
rename(`Pointwise SE` = SE)
base_loo_summary |>
format_decimal_table(digits = 2) |>
add_row(
Quantity = "Active estimated model parameters",
Estimate = format(
base_active_parameter_count,
big.mark = ",",
scientific = FALSE,
trim = TRUE
),
`Pointwise SE` = "\u2014",
.before = 2L
) |>
kable(align = c("l", "r", "r"))
} else {
knitr::asis_output(
paste0(
"*No saved base LOO calculation was available. Set ",
"`ESC31_RUN_BASE_LOO=true` while rendering the accepted posterior ",
"to create it.*"
)
)
}| Quantity | Estimate | Pointwise SE |
|---|---|---|
| Expected log predictive density (elpd_loo) | -7,732.99 | 433.95 |
| Active estimated model parameters | 1,559 | — |
| Effective number of parameters (p_loo) | 296.30 | 22.36 |
| LOO information criterion (LOO-IC) | 15,465.97 | 867.91 |
The raw count is 1,559 active estimated parameters, including annual recruitment deviations and selectivity effects. The effective count, \(p_\mathrm{LOO}\), measures predictive complexity after accounting for model structure and priors, so it is not expected to equal the raw parameter count.
The base input contains 213 LF rows. Its 51 zero-N rows remain in model I/O, while the 162 positive-N LF likelihood terms enter PSIS-LOO. No positive-N LF term has a non-finite Pareto shape estimate.
An earlier provisional cache reported four infinite LF diagnostics (LL3 in 2005–2006 and LL4 in 1986 and 1990). In each case the effective count in every length bin was below 0.5. Numeric likelihood reporting had dispatched those fractional counts to the discrete base-R multinomial implementation, which rounds bin counts to integers; every bin therefore became zero and produced an identically zero reported contribution. The fitted RTMB objective already used the continuous gamma-function extension, so this was a LOO post-processing error and does not require a base-model refit or new MCMC. The version-3 LOO cache uses the same expression in both paths and rejects the provisional cache.
base_loo_composition_pattern_text <-
"No corrected base LOO result was available."
base_loo_n_pattern_text <- ""
if (base_loo_available) {
base_loo_pareto <- tibble(
Point = rownames(base_loo$pointwise),
k = as.numeric(base_loo_k)
) |>
mutate(
Component = sub(
"_[0-9]+$",
"",
sub("^lp_", "", .data$Point)
),
Component = recode(
.data$Component,
af = "Age frequencies",
lf = "Length frequencies",
cpue_lf = "CPUE length frequencies",
cpue = "CPUE",
aerial = "Aerial survey",
troll = "Troll survey",
tags = "Conventional tags",
pop = "POP",
hsp = "HSP",
gt = "Gene tagging",
.default = .data$Component
)
) |>
group_by(.data$Component) |>
summarise(
Points = n(),
`k > reliability threshold` = sum(
!is.finite(.data$k) |
.data$k > base_loo_threshold
),
`k >= 1` = sum(
!is.finite(.data$k) |
.data$k >= 1
),
`Non-finite k` = sum(!is.finite(.data$k)),
.groups = "drop"
) |>
arrange(desc(.data$`k > reliability threshold`), .data$Component)
component_diagnostic <- function(component) {
base_loo_pareto |>
filter(.data$Component == component) |>
slice_head(n = 1L)
}
af_diagnostic <- component_diagnostic("Age frequencies")
lf_diagnostic <- component_diagnostic("Length frequencies")
cpue_lf_diagnostic <- component_diagnostic(
"CPUE length frequencies"
)
gt_diagnostic <- component_diagnostic("Gene tagging")
base_loo_composition_pattern_text <- paste0(
"The corrected pattern is concentrated in age frequencies: **",
af_diagnostic$`k > reliability threshold`,
" of ",
af_diagnostic$Points,
" AF points** exceed the 0.7 threshold (**",
af_diagnostic$`k >= 1`,
" have k >= 1**), compared with **",
lf_diagnostic$`k > reliability threshold`,
" of ",
lf_diagnostic$Points,
" catch-LF points** (**",
lf_diagnostic$`k >= 1`,
" have k >= 1**), **",
cpue_lf_diagnostic$`k > reliability threshold`,
" of ",
cpue_lf_diagnostic$Points,
" CPUE-LF points**, and **",
gt_diagnostic$`k > reliability threshold`,
" of ",
gt_diagnostic$Points,
" gene-tagging points**. No other component exceeds the threshold, ",
"and every corrected Pareto k estimate is finite."
)
composition_k_n_correlation <- function(component, input_n) {
component_points <- grepl(
paste0("^lp_", component, "_[0-9]+$"),
rownames(base_loo$pointwise)
)
component_index <- as.integer(
sub(
"^.*_([0-9]+)$",
"\\1",
rownames(base_loo$pointwise)[component_points]
)
)
stats::cor(
input_n[component_index],
base_loo_k[component_points],
method = "spearman"
)
}
af_k_n_correlation <- composition_k_n_correlation(
"af", base_fit$data$af_n
)
lf_k_n_correlation <- composition_k_n_correlation(
"lf", base_fit$data$lf_n
)
cpue_lf_k_n_correlation <- composition_k_n_correlation(
"cpue_lf", base_fit$data$cpue_n
)
base_loo_n_pattern_text <- paste0(
"Within composition types, the Spearman association between input ",
"effective N and Pareto k is **",
format_decimal(af_k_n_correlation, 2),
" for AF**, **",
format_decimal(lf_k_n_correlation, 2),
" for catch LF**, and **",
format_decimal(cpue_lf_k_n_correlation, 2),
" for CPUE LF**. Larger input N is therefore associated with greater ",
"influence for the two length-composition series, but it does not ",
"explain the broad AF pattern."
)
base_loo_pareto |>
kable(align = c("l", "r", "r", "r", "r"))
}| Component | Points | k > reliability threshold | k >= 1 | Non-finite k |
|---|---|---|---|---|
| Age frequencies | 88 | 49 | 14 | 0 |
| Length frequencies | 162 | 15 | 2 | 0 |
| CPUE length frequencies | 57 | 6 | 0 | 0 |
| Gene tagging | 8 | 1 | 0 | 0 |
| Aerial survey | 1 | 0 | 0 | 0 |
| CPUE | 57 | 0 | 0 | 0 |
| Conventional tags | 69 | 0 | 0 | 0 |
| HSP | 116 | 0 | 0 | 0 |
| POP | 9005 | 0 | 0 | 0 |
if (base_loo_available) {
loo_fishery_labels <- c(
"LL1", "LL2", "LL3", "LL4", "Indonesian", "Australian", "CPUE"
)
loo_input_year <- function(year_index) {
base_fit$data$first_yr + as.integer(year_index) - 1L
}
collapse_loo_years <- function(years) {
years <- sort(unique(years[!is.na(years)]))
if (length(years)) paste(years, collapse = ", ") else "—"
}
base_loo_high_k_years <- tibble(
Point = rownames(base_loo$pointwise),
k = as.numeric(base_loo_k)
) |>
mutate(
Point_index = as.integer(sub("^.*_([0-9]+)$", "\\1", .data$Point)),
Component = sub(
"_[0-9]+$",
"",
sub("^lp_", "", .data$Point)
),
Series = case_when(
.data$Component == "af" ~ paste(
loo_fishery_labels[
as.integer(base_fit$data$af_fishery[.data$Point_index])
],
"age frequencies"
),
.data$Component == "lf" ~ paste(
loo_fishery_labels[
as.integer(base_fit$data$lf_fishery[.data$Point_index])
],
"length frequencies"
),
.data$Component == "cpue_lf" ~ "CPUE length frequencies",
.data$Component == "gt" ~ "GT",
TRUE ~ NA_character_
),
Year = case_when(
.data$Component == "af" ~ loo_input_year(
base_fit$data$af_year[.data$Point_index]
),
.data$Component == "lf" ~ loo_input_year(
base_fit$data$lf_year[.data$Point_index]
),
.data$Component == "cpue_lf" ~ loo_input_year(
base_fit$data$cpue_lf_years[.data$Point_index]
),
.data$Component == "gt" ~ as.integer(
base_fit$data$gt_obs$RecYear[.data$Point_index]
),
TRUE ~ NA_integer_
)
) |>
filter(is.finite(.data$k), .data$k > base_loo_threshold)
if (anyNA(base_loo_high_k_years$Series) ||
anyNA(base_loo_high_k_years$Year)) {
stop("A high-k point could not be matched to its data series and year.")
}
base_loo_high_k_years |>
group_by(.data$Series) |>
summarise(
`High-k points` = n(),
`Years with k > 0.7` = collapse_loo_years(.data$Year),
`Years with k >= 1.0` = collapse_loo_years(
.data$Year[.data$k >= 1]
),
`Maximum k` = max(.data$k),
`Year of maximum k` = .data$Year[which.max(.data$k)][1L],
.groups = "drop"
) |>
arrange(desc(.data$`Maximum k`)) |>
mutate(`Maximum k` = format_decimal(.data$`Maximum k`, 3L)) |>
kable(align = c("l", "r", "l", "l", "r", "r")) |>
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover", "condensed"),
full_width = TRUE
)
}| Series | High-k points | Years with k > 0.7 | Years with k >= 1.0 | Maximum k | Year of maximum k |
|---|---|---|---|---|---|
| Australian age frequencies | 40 | 1965, 1966, 1967, 1968, 1969, 1970, 1972, 1973, 1976, 1978, 1985, 1987, 1989, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2020, 2021, 2022, 2023, 2024, 2025 | 1966, 1973, 1978, 1989, 2010, 2011, 2012, 2013, 2018, 2022 | 1.234 | 1978 |
| LL3 length frequencies | 4 | 1959, 1960, 1961, 1962 | 1960 | 1.131 | 1960 |
| Indonesian age frequencies | 9 | 1999, 2003, 2008, 2013, 2014, 2015, 2017, 2018, 2020 | 2013, 2014, 2017, 2018 | 1.129 | 2018 |
| LL1 length frequencies | 11 | 1952, 1997, 1998, 2006, 2007, 2008, 2010, 2011, 2013, 2024, 2025 | 2006 | 1.064 | 2006 |
| GT | 1 | 2018 | — | 0.993 | 2018 |
| CPUE length frequencies | 6 | 1985, 2006, 2007, 2008, 2018, 2022 | — | 0.929 | 2018 |
The corrected pattern is concentrated in age frequencies: 49 of 88 AF points exceed the 0.7 threshold (14 have k >= 1), compared with 15 of 162 catch-LF points (2 have k >= 1), 6 of 57 CPUE-LF points, and 1 of 8 gene-tagging points. No other component exceeds the threshold, and every corrected Pareto k estimate is finite.
Within composition types, the Spearman association between input effective N and Pareto k is 0.06 for AF, 0.69 for catch LF, and 0.54 for CPUE LF. Larger input N is therefore associated with greater influence for the two length-composition series, but it does not explain the broad AF pattern.
The Monte Carlo SE of elpd_loo is NA. This is expected from the loo implementation when any Pareto \(k\) exceeds its sample-size-dependent reliability threshold. For these 3,000 retained draws the threshold is 0.700, and 71 pointwise estimates exceed it. The reported LOO-IC can still be inspected, but its importance-sampling approximation and SE require caution; exact leave-one-out refits, moment matching, or K-fold validation are appropriate follow-up options for influential observations (Vehtari et al. 2024).
Pareto \(k\) diagnoses whether importance sampling can reliably approximate refitting after leaving out one point; it is not, by itself, an estimate that a data set has been given too much statistical weight (Vehtari et al. 2024). Here one AF or LF point is an entire annual composition, not an individual fish or length/age bin. Leaving out a whole composition can therefore move the posterior appreciably, particularly where composition data inform recruitment or selectivity.
The base model uses multinomial composition likelihoods with fixed input effective sample sizes. Multinomial sampling assumes independent fish, whereas schooling, clustered sampling, ageing error, and other process or observation errors can produce correlated, overdispersed compositions (Francis 2011; Hulson et al. 2011). This is a reason to review composition weighting and likelihood choice, but not to lower effective N automatically until a high-\(k\) row has been checked for data error, residual or posterior-predictive misfit, and conflict with abundance data. Changing effective N also changes the likelihood scale, so raw LOO-IC values should not be used to select among different weighting schemes as though they were scores for the same observation model. Appropriate sensitivity work includes exact leave-one-composition-out or K-fold refits for problematic rows, externally justified effective-N alternatives, and composition likelihoods that allow overdispersion or correlation. Francis’s logistic-normal formulation is one such alternative (Francis 2014); the package now supports it for age frequencies, but the accepted base remains unchanged.
if (base_loo_available) {
plot_loo(
x = base_loo,
data = base_fit$data,
exclude = "pop"
)
}
if (base_loo_available) {
plot_loo(
x = base_loo,
data = base_fit$data,
include = "pop"
)
}
The following overview links the recent observations represented in the accepted base-model inputs to the exact PSIS-LOO diagnostics from the accepted posterior. The year is the recapture year for GT, the offspring cohort for POP, and the newer cohort in each HSP pair; conventional-tag release cohorts predate the displayed period. A grey marker denotes a series that ended before the corresponding year. Troll-survey observations are available but are not fitted because the base model has troll_switch = 0; consequently, there is no troll likelihood contribution or Pareto-\(k\) diagnostic. Marker area is scaled within each row and is not comparable among data types. Fitted rows are ordered by mean Pareto \(k\), from lowest at the bottom to highest at the top; the unfitted troll row is kept separately at the bottom. Pareto \(k\) diagnoses the reliability of the leave-one-out approximation; it does not measure the direction or magnitude of a data set’s effect on stock status.
plot_catch_input_output(base_fit)
plot_cpue(base_fit, posterior = base_plot_posterior)
plot_aerial(base_fit)
plot_gt(base_fit)
plot_tags_fits_pooled(base_fit)
close_kin_plots <- plot_hsps(base_fit, return_list = TRUE)
close_kin_plots$pop_by_juvenile_cohort
close_kin_plots$pop_by_adult_capture_age
close_kin_plots$pop_by_adult_capture_year
close_kin_plots$hsp_by_cohort_pair
close_kin_plots$hsp_by_initial_cohort
close_kin_plots$total_matches
composition_selectivity_blocks(
plot_lf(base_fit, posterior = base_plot_posterior, fishery = "LL1"),
base_fit,
"LL1"
)
composition_selectivity_blocks(
plot_lf(base_fit, posterior = base_plot_posterior, fishery = "LL2"),
base_fit,
"LL2"
)
composition_selectivity_blocks(
plot_lf(base_fit, posterior = base_plot_posterior, fishery = "LL3"),
base_fit,
"LL3"
)
composition_selectivity_blocks(
plot_lf(base_fit, posterior = base_plot_posterior, fishery = "LL4"),
base_fit,
"LL4"
)
composition_selectivity_blocks(
plot_lf(base_fit, posterior = base_plot_posterior, fishery = "CPUE"),
base_fit,
"CPUE"
)
composition_selectivity_blocks(
plot_af(base_fit, posterior = base_plot_posterior, fishery = "Indonesian"),
base_fit,
"Indonesian"
)
composition_selectivity_blocks(
plot_af(base_fit, posterior = base_plot_posterior, fishery = "Australian"),
base_fit,
"Australian"
)
plot_selectivity(base_fit, fisheries = "LL1")
plot_selectivity(base_fit, fisheries = "LL2")
plot_selectivity(base_fit, fisheries = "LL3")
plot_selectivity(base_fit, fisheries = "LL4")
plot_selectivity(base_fit, fisheries = "Indonesian")
plot_selectivity(base_fit, fisheries = "Australian")
plot_selectivity(base_fit, fisheries = "CPUE")
One-step-ahead (OSA) residuals provide randomized quantile-style diagnostics for non-Gaussian observations (Dunn and Smyth 1996) and are useful for assessing stock-assessment fit across observation types. For composition data, the Francis and McAllister-Ianelli diagnostics summarize whether the input sample sizes are broadly consistent with the dispersion in the composition residuals.
The scalar likelihood components use oneStepPredict() because each observation can be evaluated directly as a univariate OSA residual. The composition and tag recapture likelihoods are constrained multivariate observations, so compResidual is used when the residuals need to respect the simplex and Dirichlet-multinomial structure or when oneStepPredict() returns non-finite composition residuals. The residual source used for each data type is also reported in Table 12. Decision esc31_2026_conventional_tag_compResidual_osa_fallback_v1 accepts this likelihood-matched conventional-tag fallback as the production diagnostic; it does not change the fitted likelihood or any model estimate.
For normalized residuals, SDNR values near 1 indicate that the spread of the residuals is broadly consistent with the assumed observation error. High SDNR values indicate residuals that are more variable than expected, which can point to lack of fit or observation errors that are too small. Low SDNR values indicate residuals that are less variable than expected, which can point to observation errors that are too large or an overly down-weighted data set. MAR is the median absolute residual; values near 0.67 are expected for standard normal residuals, with larger values indicating larger typical residuals and smaller values indicating smaller typical residuals.
The CPUE and composition likelihoods are deliberately down-weighted in this base model. Their SDNRs are therefore expected to be below one and are not a request to increase their weights until the SDNRs equal one. Those diagnostics remain useful for finding temporal or age/length structure in the residuals; the intended weighting changes their interpretation, not the need to inspect their pattern.
p_cpue_residuals <- plot_cpue_residuals(base_fit)
Diagnostic Summary:
SDNR: 0.79539 (95% CI 0.67151-0.97575; Target: ~1.0)
MAR: 0.60767 (Target: ~0.67)p_cpue_residuals
p_aerial_residuals <- plot_aerial_residuals(base_fit)
Diagnostic Summary:
SDNR: 0.99859 (95% CI 0.75942-1.45851; Target: ~1.0)
MAR: 0.45206 (Target: ~0.67)
Residuals summarised: 20 aerial survey observationsp_aerial_residuals
p_gt_residuals <- plot_gt_residuals(base_fit)
Diagnostic Summary:
SDNR: 0.84899 (95% CI 0.56133-1.72792; Target: ~1.0)
MAR: 0.60246 (Target: ~0.67)p_gt_residuals
p_hsp_residuals <- plot_hsps_residuals(base_fit)
Diagnostic Summary:
SDNR: 1.05670 (95% CI 0.93600-1.21342; Target: ~1.0)
MAR: 0.67470 (Target: ~0.67)p_hsp_residuals
p_pop_residuals <- plot_pops_residuals(base_fit)
Diagnostic Summary:
SDNR: 0.99429 (95% CI 0.97998-1.00903; Target: ~1.0)
MAR: 0.65799 (Target: ~0.67)p_pop_residuals
tag_osa_decision <- review_method_decisions$conventional_tag_osa
stopifnot(
identical(tag_osa_decision$status, "accepted"),
identical(
tag_osa_decision$decision_id,
"esc31_2026_conventional_tag_compResidual_osa_fallback_v1"
),
identical(tag_osa_decision$method, "compResidual::resDirM")
)
p_tag_residuals <- plot_tags_residuals(
base_fit,
seed = tag_osa_decision$seed
)
Diagnostic Summary:
SDNR: 1.00933 (95% CI 0.92508-1.11059; Target: ~1.0)
MAR: 0.57735 (Target: ~0.67)
Residuals summarised: 232 recapture categoriesp_tag_residuals
p_lf_ll1_residuals <- plot_lf_residuals(base_fit, fishery = "LL1")
Diagnostic Summary:
SDNR: 0.54418 (95% CI 0.52686-0.56269; Target: ~1.0)
MAR: 0.34754 (Target: ~0.67)
Francis mean length SDNR: 1.08408 (N multiplier: 0.851)
McAllister-Ianelli N: harmonic 219.3, median 415.7, mean input 69.0 (harmonic/input: 3.179)
Residuals summarised: 1776 length-composition cells across 74 compositionsp_lf_ll1_residuals
p_lf_ll2_residuals <- plot_lf_residuals(base_fit, fishery = "LL2")
Diagnostic Summary:
SDNR: 0.63413 (95% CI 0.60138-0.67067; Target: ~1.0)
MAR: 0.41887 (Target: ~0.67)
Francis mean length SDNR: 0.78429 (N multiplier: 1.626)
McAllister-Ianelli N: harmonic 82.0, median 167.0, mean input 8.8 (harmonic/input: 9.278)
Residuals summarised: 648 length-composition cells across 27 compositions
Non-finite residual cells dropped before plotting: 48
Composition diagnostics based on 29 compositionsp_lf_ll2_residuals
p_lf_ll3_residuals <- plot_lf_residuals(base_fit, fishery = "LL3")
Diagnostic Summary:
SDNR: 0.70838 (95% CI 0.66519-0.75760; Target: ~1.0)
MAR: 0.47688 (Target: ~0.67)
Francis mean length SDNR: 0.92095 (N multiplier: 1.179)
McAllister-Ianelli N: harmonic 78.2, median 140.5, mean input 34.3 (harmonic/input: 2.283)
Residuals summarised: 456 length-composition cells across 19 compositions
Non-finite residual cells dropped before plotting: 48
Composition diagnostics based on 21 compositionsp_lf_ll3_residuals
p_lf_ll4_residuals <- plot_lf_residuals(base_fit, fishery = "LL4")
Diagnostic Summary:
SDNR: 0.66794 (95% CI 0.63007-0.71069; Target: ~1.0)
MAR: 0.43576 (Target: ~0.67)
Francis mean length SDNR: 0.82127 (N multiplier: 1.483)
McAllister-Ianelli N: harmonic 58.0, median 63.6, mean input 8.8 (harmonic/input: 6.578)
Residuals summarised: 532 length-composition cells across 38 compositionsp_lf_ll4_residuals
p_lf_cpue_residuals <- plot_lf_residuals(base_fit, fishery = "CPUE")
Diagnostic Summary:
SDNR: 0.61895 (95% CI 0.59394-0.64618; Target: ~1.0)
MAR: 0.36456 (Target: ~0.67)
Francis mean length SDNR: 0.75987 (N multiplier: 1.732)
McAllister-Ianelli N: harmonic 255.3, median 295.1, mean input 80.1 (harmonic/input: 3.187)
Residuals summarised: 1083 length-composition cells across 57 compositionsp_lf_cpue_residuals
p_af_indonesia_residuals <- plot_af_residuals(base_fit, fishery = "Indonesian")
Diagnostic Summary:
SDNR: 0.61883 (95% CI 0.58630-0.65522; Target: ~1.0)
MAR: 0.40651 (Target: ~0.67)
Francis mean age SDNR: 1.04836 (N multiplier: 0.910)
McAllister-Ianelli N: harmonic 223.7, median 224.4, mean input 86.6 (harmonic/input: 2.583)
Residuals summarised: 624 age-composition cells across 26 compositionsp_af_indonesia_residuals
p_af_australia_residuals <- plot_af_residuals(base_fit, fishery = "Australian")
Diagnostic Summary:
SDNR: 0.77326 (95% CI 0.72501-0.82844; Target: ~1.0)
MAR: 0.41629 (Target: ~0.67)
Francis mean age SDNR: 0.67260 (N multiplier: 2.210)
McAllister-Ianelli N: harmonic 65.4, median 174.4, mean input 27.6 (harmonic/input: 2.370)
Residuals summarised: 434 age-composition cells across 62 compositionsp_af_australia_residuals
residual_plots <- list(
`CPUE index` = p_cpue_residuals,
`Aerial survey` = p_aerial_residuals,
`Gene tagging` = p_gt_residuals,
`Half-sibling pairs` = p_hsp_residuals,
`Close-kin POPs` = p_pop_residuals,
`Conventional tags` = p_tag_residuals,
`Length composition: LL1` = p_lf_ll1_residuals,
`Length composition: LL2` = p_lf_ll2_residuals,
`Length composition: LL3` = p_lf_ll3_residuals,
`Length composition: LL4` = p_lf_ll4_residuals,
`Length composition: CPUE` = p_lf_cpue_residuals,
`Age composition: Indonesian` = p_af_indonesia_residuals,
`Age composition: Australian` = p_af_australia_residuals
)
residual_statistics <- imap_dfr(
residual_plots,
~ residual_stat_row(.y, .x)
)
residual_statistics |>
rename(
Estimate = SDNR,
`Lower 95%` = `SDNR lower 95%`,
`Upper 95%` = `SDNR upper 95%`,
`SDNR` = `Francis SDNR`,
`N multiplier` = `Francis N multiplier`,
`Harmonic N` = `McAllister-Ianelli harmonic N`,
`Median N` = `McAllister-Ianelli median N`,
`Effective/input N` = `McAllister-Ianelli/Input N`
) |>
mutate(Source = format_residual_source(Source)) |>
select(-Source, Source) |>
format_residual_table(digits = 3) |>
replace_missing() |>
kable(
align = c("l", rep("r", 11), "l"),
caption = "Residual diagnostic statistics by data type. SDNR is the standard deviation of normalized residuals; MAR is the median absolute residual. Francis and McAllister-Ianelli diagnostics are reported for age- and length-composition data where available."
) |>
add_header_above(c(" " = 2, "SDNR" = 3, " " = 1, "Francis" = 2, "McAllister-Ianelli" = 4, " " = 1))| Data type | N | Estimate | Lower 95% | Upper 95% | MAR | SDNR | N multiplier | Harmonic N | Median N | Mean input N | Effective/input N | Source |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| CPUE index | 57 | 0.795 | 0.672 | 0.976 | 0.608 | - | - | - | - | - | - | oneStepPredict |
| Aerial survey | 20 | 0.999 | 0.759 | 1.459 | 0.452 | - | - | - | - | - | - | oneStepPredict |
| Gene tagging | 8 | 0.849 | 0.561 | 1.728 | 0.602 | - | - | - | - | - | - | oneStepPredict |
| Half-sibling pairs | 116 | 1.057 | 0.936 | 1.213 | 0.675 | - | - | - | - | - | - | oneStepPredict |
| Close-kin POPs | 9005 | 0.994 | 0.980 | 1.009 | 0.658 | - | - | - | - | - | - | oneStepPredict |
| Conventional tags | 232 | 1.009 | 0.925 | 1.111 | 0.577 | - | - | - | - | - | - | compResidual |
| Length composition: LL1 | 1776 | 0.544 | 0.527 | 0.563 | 0.348 | 1.084 | 0.851 | 219.274 | 415.689 | 68.982 | 3.179 | compResidual |
| Length composition: LL2 | 648 | 0.634 | 0.601 | 0.671 | 0.419 | 0.784 | 1.626 | 82.019 | 167.001 | 8.840 | 9.278 | compResidual |
| Length composition: LL3 | 456 | 0.708 | 0.665 | 0.758 | 0.477 | 0.921 | 1.179 | 78.245 | 140.489 | 34.271 | 2.283 | compResidual |
| Length composition: LL4 | 532 | 0.668 | 0.630 | 0.711 | 0.436 | 0.821 | 1.483 | 57.993 | 63.571 | 8.816 | 6.578 | compResidual |
| Length composition: CPUE | 1083 | 0.619 | 0.594 | 0.646 | 0.365 | 0.760 | 1.732 | 255.289 | 295.078 | 80.103 | 3.187 | compResidual |
| Age composition: Indonesian | 624 | 0.619 | 0.586 | 0.655 | 0.407 | 1.048 | 0.910 | 223.690 | 224.445 | 86.600 | 2.583 | compResidual |
| Age composition: Australian | 434 | 0.773 | 0.725 | 0.828 | 0.416 | 0.673 | 2.210 | 65.420 | 174.449 | 27.606 | 2.370 | compResidual |
sdnr_plot_data <- residual_statistics |>
filter(is.finite(SDNR)) |>
mutate(
`Data type` = factor(
`Data type`,
levels = rev(residual_statistics$`Data type`)
)
)
ggplot(sdnr_plot_data, aes(y = `Data type`)) +
geom_vline(
xintercept = 1,
color = "#c65d36",
linewidth = 0.7,
linetype = "dashed"
) +
geom_segment(
aes(
x = `SDNR lower 95%`,
xend = `SDNR upper 95%`,
yend = `Data type`
),
color = "#234e70",
linewidth = 0.8,
na.rm = TRUE
) +
geom_point(aes(x = SDNR), color = "#234e70", size = 2.7) +
scale_x_continuous(limits = c(0, NA), expand = expansion(mult = c(0, 0.05))) +
labs(x = "SDNR (95% confidence interval)", y = NULL) +
theme(panel.grid.major.y = element_blank())
The aerial-survey tau and conventional-tag variance factor were calibrated together and are fixed at 0.59 and 2.40, respectively, for production fitting and sampling. At the calibrated MLE, their SDNRs are approximately 0.999 and 1.009. The randomized conventional-tag diagnostic uses seed 123, without changing the caller’s random-number state, so the reported result is reproducible. Both likelihood controls remain fixed in the posterior.
The gene-tagging SDNR is based on only eight observations, and its 95% interval includes one. The binomial gene-tagging likelihood is therefore retained unchanged. The available beta-binomial alternative adds overdispersion and would lower the SDNR further, so no GT variance or effective-sample-size adjustment is applied.
plot_initial_numbers(base_fit, color = fit_expected_color) +
geom_line(
data = filter(
initial_numbers_comparison,
Assessment == "Current assessment"
),
aes(Age, Value, linetype = Assessment),
colour = assessment_line_colours[["Current assessment"]],
linewidth = 0.8,
inherit.aes = FALSE
) +
geom_line(
data = filter(
initial_numbers_comparison,
Assessment == "Previous assessment"
),
aes(Age, Value, linetype = Assessment),
colour = assessment_line_colours[["Previous assessment"]],
linewidth = 0.8,
inherit.aes = FALSE
) +
scale_linetype_manual(
values = assessment_linetypes,
breaks = names(assessment_linetypes)
) +
labs(linetype = NULL) +
guides(linetype = assessment_line_guide) +
theme(legend.position = "bottom")
plot_natural_mortality(
base_fit,
posterior = base_plot_posterior,
color = fit_expected_color
) +
geom_line(
data = filter(mortality_comparison, Assessment == "Current assessment"),
aes(Age, Value, linetype = Assessment),
colour = assessment_line_colours[["Current assessment"]],
linewidth = 0.8,
inherit.aes = FALSE
) +
geom_line(
data = filter(mortality_comparison, Assessment == "Previous assessment"),
aes(Age, Value, linetype = Assessment),
colour = assessment_line_colours[["Previous assessment"]],
linewidth = 0.8,
inherit.aes = FALSE
) +
scale_linetype_manual(
values = assessment_linetypes,
breaks = names(assessment_linetypes)
) +
labs(linetype = NULL) +
guides(linetype = assessment_line_guide) +
theme(legend.position = "bottom")
mortality_report <- base_dynamics_report
mortality_minimum_index <- which.min(mortality_report$M_a)
mortality_summary <- tibble(
Quantity = c("M0", "M4", "M10", "M30", "Minimum M"),
Age = c(
"0", "4", "10", "30",
as.character(base_fit$data$min_age + mortality_minimum_index - 1L)
),
Value = c(
mortality_report$par_m0,
mortality_report$par_m4,
mortality_report$par_m10,
mortality_report$par_m30,
mortality_report$M_a[[mortality_minimum_index]]
)
)
mortality_summary |>
format_decimal_table(digits = 6) |>
kable(align = c("l", "r", "r"))| Quantity | Age | Value |
|---|---|---|
| M0 | 0 | 0.365742 |
| M4 | 4 | 0.161674 |
| M10 | 10 | 0.107361 |
| M30 | 30 | 0.457715 |
| Minimum M | 25 | 0.089302 |
The following profiles re-optimize all remaining parameters while fixing each estimated log-mortality parameter in turn. The vertical dashed line is the MLE and the horizontal dashed line is a two-unit increase in the objective. Profile results are stored separately from Quarto’s execution cache and are invalidated whenever the accepted base-fit scientific identity, optimizer result, profile settings, or profile-cache scientific contract changes. Package and function versions remain audit metadata rather than cross-machine invalidators. Set ESC31_RUN_BASE_M_PROFILES=true to recalculate them, or ESC31_RENDER_BASE_M_PROFILES=true to display a compatible saved result.
m_profile_parameters <- c(
"par_log_m10", "par_log_m30"
)
m_profile_labels <- c(
par_log_m10 = expression(M[10]),
par_log_m30 = expression(M[30])
)
m_profile_settings <- list(
h = 1e-4,
ytol = 2,
ystep = 0.1,
slice = FALSE,
adaptive = TRUE
)
m_profile_identity <- esc31_object_md5(list(
contract = "esc31_natural_mortality_profile_cache_v2",
base_fit_scientific_signature = esc31_fit_scientific_signature(base_fit),
settings = m_profile_settings,
parameters = m_profile_parameters
))
if (render_m_profile_results) {
if (run_new_m_profiles) {
profile_obj <- sbt_obj(base_fit, fresh = TRUE)
profile_bounds <- base_fit$bounds
m_profiles <- setNames(
lapply(m_profile_parameters, function(parameter) {
parameter_index <- match(parameter, names(profile_obj$par))
if (is.na(parameter_index)) {
stop("The base MLE does not estimate `", parameter, "`.",
call. = FALSE)
}
sbtprofile(
obj = profile_obj,
name = parameter,
h = m_profile_settings$h,
ytol = m_profile_settings$ytol,
ystep = m_profile_settings$ystep,
parm.range = c(
profile_bounds$lower[[parameter_index]],
profile_bounds$upper[[parameter_index]]
),
slice = m_profile_settings$slice,
adaptive = m_profile_settings$adaptive,
trace = FALSE
)
}),
m_profile_parameters
)
saveRDS(
list(
cache_schema_version = 2L,
scientific_identity = m_profile_identity,
profiles = m_profiles
),
m_profile_file
)
} else {
if (!file.exists(m_profile_file)) {
stop("The saved natural-mortality profiles are missing.", call. = FALSE)
}
m_profile_cache <- readRDS(m_profile_file)
m_profile_cache_has_current_identity <-
is.list(m_profile_cache) &&
identical(m_profile_cache$cache_schema_version, 2L) &&
identical(m_profile_cache$scientific_identity, m_profile_identity)
m_profile_cache_is_legacy <-
is.list(m_profile_cache) &&
is.null(m_profile_cache$cache_schema_version) &&
is.character(m_profile_cache$identity) &&
length(m_profile_cache$identity) == 1L
m_profile_payload_is_compatible <-
is.list(m_profile_cache) &&
(m_profile_cache_has_current_identity || m_profile_cache_is_legacy) &&
identical(names(m_profile_cache$profiles), m_profile_parameters) &&
all(vapply(m_profile_parameters, function(parameter) {
profile <- m_profile_cache$profiles[[parameter]]
mle_x <- base_fit$fit$opt$par[
match(parameter, names(base_fit$fit$opt$par))
]
is.data.frame(profile) &&
all(c(parameter, "value") %in% names(profile)) &&
all(is.finite(profile[[parameter]])) &&
all(is.finite(profile$value)) &&
min(profile[[parameter]]) < mle_x &&
max(profile[[parameter]]) > mle_x &&
isTRUE(all.equal(
profile[[parameter]][which.min(profile$value)],
unname(mle_x),
tolerance = 1e-10
)) &&
isTRUE(all.equal(
min(profile$value),
unname(base_fit$fit$opt$objective),
tolerance = 1e-8
))
}, logical(1)))
if (!m_profile_payload_is_compatible) {
stop("The saved natural-mortality profiles are stale.", call. = FALSE)
}
profile_obj <- sbt_obj(base_fit, fresh = TRUE)
m_profiles <- m_profile_cache$profiles
}
}if (render_m_profile_results) {
plot_profile(
obj = profile_obj,
x = m_profiles[["par_log_m10"]],
xlab = m_profile_labels[["par_log_m10"]],
rescale = TRUE
) +
theme(legend.position = "bottom")
}
if (render_m_profile_results) {
plot_profile(
obj = profile_obj,
x = m_profiles[["par_log_m30"]],
xlab = m_profile_labels[["par_log_m30"]],
rescale = TRUE
) +
theme(legend.position = "bottom")
}
if (render_m_profile_results) {
m_profile_summary <- imap_dfr(
m_profiles,
function(profile, parameter) {
profile_x <- as.numeric(profile[[parameter]])
profile_delta <- as.numeric(profile$value) - min(profile$value)
mle_x <- base_fit$fit$opt$par[
match(parameter, names(base_fit$fit$opt$par))
]
inside <- is.finite(profile_delta) & profile_delta <= 2
tibble(
Parameter = parameter_labels(sub("par_log_m", "M", parameter)),
MLE = exp(mle_x),
`Lower profile value` = exp(min(profile_x[inside])),
`Upper profile value` = exp(max(profile_x[inside])),
`Lower side closed` = any(profile_delta[profile_x < mle_x] >= 2),
`Upper side closed` = any(profile_delta[profile_x > mle_x] >= 2)
)
}
)
m_profile_summary |>
mutate(across(where(is.numeric), ~ format_decimal(.x, 4))) |>
kable(align = c("l", "r", "r", "r", "c", "c"))
} else {
m_profile_summary <- tibble()
}| Parameter | MLE | Lower profile value | Upper profile value | Lower side closed | Upper side closed |
|---|---|---|---|---|---|
| M10 | 0.1074 | 0.1004 | 0.1144 | TRUE | TRUE |
| M30 | 0.4577 | 0.3730 | 0.5617 | TRUE | TRUE |
plot_rec_devs(base_fit, posterior = base_plot_posterior) +
geom_line(
data = filter(
recruitment_deviate_comparison,
Assessment == "Current assessment"
),
aes(Year, Value, linetype = Assessment),
colour = assessment_line_colours[["Current assessment"]],
linewidth = 0.8,
inherit.aes = FALSE
) +
geom_line(
data = filter(
recruitment_deviate_comparison,
Assessment == "Previous assessment"
),
aes(Year, Value, linetype = Assessment),
colour = assessment_line_colours[["Previous assessment"]],
linewidth = 0.8,
inherit.aes = FALSE
) +
scale_linetype_manual(
values = assessment_linetypes,
breaks = names(assessment_linetypes)
) +
labs(linetype = NULL) +
guides(linetype = assessment_line_guide) +
theme(legend.position = "bottom")
plot_recruitment(base_fit, posterior = base_plot_posterior) +
geom_line(
data = filter(recruitment_comparison, Assessment == "Current assessment"),
aes(Year, Value, linetype = Assessment),
colour = assessment_line_colours[["Current assessment"]],
linewidth = 0.8,
inherit.aes = FALSE
) +
geom_line(
data = filter(recruitment_comparison, Assessment == "Previous assessment"),
aes(Year, Value, linetype = Assessment),
colour = assessment_line_colours[["Previous assessment"]],
linewidth = 0.8,
inherit.aes = FALSE
) +
scale_linetype_manual(
values = assessment_linetypes,
breaks = names(assessment_linetypes)
) +
labs(linetype = NULL) +
guides(linetype = assessment_line_guide) +
theme(legend.position = "bottom")
base_biomass_report <- base_dynamics_report
absolute_spawning_biomass_comparison <- bind_rows(
data.frame(
Year = seq.int(
base_fit$data$first_yr,
length.out = length(base_biomass_report$spawning_biomass_y)
),
Spawning_biomass = as.numeric(base_biomass_report$spawning_biomass_y),
Assessment = "Current assessment"
),
data.frame(
Year = seq.int(
previous_v1_fit$data$first_yr,
length.out = length(previous_v1_report$spawning_biomass_y)
),
Spawning_biomass = as.numeric(previous_v1_report$spawning_biomass_y),
Assessment = "Previous assessment"
)
)
ggplot(
absolute_spawning_biomass_comparison,
aes(Year, Spawning_biomass, colour = Assessment, linetype = Assessment)
) +
geom_line(linewidth = 0.9) +
scale_colour_manual(values = c(
"Current assessment" = "#0072B2",
"Previous assessment" = "black"
)) +
scale_linetype_manual(values = c(
"Current assessment" = "solid",
"Previous assessment" = "dashed"
)) +
scale_y_continuous(
limits = c(0, NA),
labels = label_comma(),
expand = expansion(mult = c(0, 0.04))
) +
labs(
x = NULL,
y = "TRO",
colour = NULL,
linetype = NULL
) +
theme(legend.position = "bottom")
spawning_biomass_comparison <- bind_rows(
data.frame(
Year = seq.int(
base_fit$data$first_yr,
length.out = length(base_biomass_report$spawning_biomass_y)
),
Relative_SB = base_biomass_report$spawning_biomass_y /
base_biomass_report$B0,
Assessment = "Current assessment"
),
data.frame(
Year = seq.int(
previous_v1_fit$data$first_yr,
length.out = length(previous_v1_report$spawning_biomass_y)
),
Relative_SB = previous_v1_report$spawning_biomass_y /
previous_v1_report$par_B0,
Assessment = "Previous assessment"
)
)
ggplot(
spawning_biomass_comparison,
aes(Year, Relative_SB, colour = Assessment, linetype = Assessment)
) +
geom_line(linewidth = 0.9) +
scale_colour_manual(values = c(
"Current assessment" = "#0072B2",
"Previous assessment" = "black"
)) +
scale_linetype_manual(values = c(
"Current assessment" = "solid",
"Previous assessment" = "dashed"
)) +
scale_y_continuous(
limits = c(0, NA),
expand = expansion(mult = c(0, 0.04))
) +
labs(x = NULL, y = "Relative TRO", colour = NULL, linetype = NULL) +
theme(legend.position = "bottom")
The legacy equilibrium MSY calculation is applied separately to every fitted catch year. Each annual calculation holds that year’s fishery selectivities, mean weights, and model-predicted catch allocation fixed while solving for the six fully selected fishing mortalities that maximize equilibrium yield. These are equilibrium reference points, not a dynamic forecast or a catch recommendation.
The annual MLE calculation is retained as a reproducible point estimate. When the MCMC chunks are enabled, the posterior calculation uses all 3,000 retained draws across all fitted catch years. Both calculations are shown below.
base_msy <- run_msy(base_fit, posterior = FALSE)base_posterior_msy <- NULL
if (base_mcmc_accepted) {
base_posterior_msy_file <- file.path(
run_dir, "msy", "base_posterior_msy.rds"
)
base_posterior_msy_years <- seq.int(
base_fit$data$first_yr_catch,
base_fit$data$last_yr
)
base_posterior_msy_identity <- esc31_base_msy_identity(
base_fit,
base_posterior_msy_years
)
base_posterior_msy_cache <- if (file.exists(base_posterior_msy_file)) {
tryCatch(
readRDS(base_posterior_msy_file),
error = function(error) NULL
)
} else {
NULL
}
if (!is.null(base_posterior_msy_cache) &&
!esc31_base_msy_cache_is_current(
base_posterior_msy_cache,
base_posterior_msy_identity,
base_fit
)) {
stop(
"The current posterior MSY cache is missing or stale. Run ",
"`Rscript scripts/run-esc31-base-msy.R .` before rendering ",
"posterior results.",
call. = FALSE
)
}
if (!is.null(base_posterior_msy_cache)) {
posterior_msy_objectives <- base_posterior_msy_cache$result$summary
required_msy_objectives <- c("phase1_objective", "objective")
if (!is.data.frame(posterior_msy_objectives) ||
length(setdiff(
required_msy_objectives,
names(posterior_msy_objectives)
)) ||
any(!is.finite(as.matrix(
posterior_msy_objectives[required_msy_objectives]
))) ||
any(posterior_msy_objectives$phase1_objective >= 0) ||
any(posterior_msy_objectives$objective >= 0)) {
stop(
"The posterior MSY cache contains a false near-zero-yield optimizer ",
"solution and cannot be presented.",
call. = FALSE
)
}
base_posterior_msy <- base_posterior_msy_cache$result
}
}All 74 fitted catch years from 1952 to 2025 passed both optimization phases and the finite-result checks at the MLE. The posterior calculation retained all 222,000 draw-year cells; every cell passed the biological-state, optimization, and finite-result checks. TRO_MSY is on the model’s normalized reproductive-output scale and is unitless; MSY and total age-2+ biomass at MSY are in tonnes.
The following table is the current counterpart of Table 4 in the 2023 assessment. Every quantity is calculated over the same accepted balanced 2,000-draw nine-cell MCMC grid. Relative age-10-plus biomass follows the legacy definition: season-1 numbers aged 10 and older, weighted using the LL1 weight-at-age schedule. The 2026 state has no 2026 fitted catch, so its MSY quantities use the latest fitted catch year, 2025; the reference year is shown explicitly.
reference_set_status_file <- file.path(
esc_dir,
"report_data", "table4_grid_stock_status.csv"
)
if (!file.exists(reference_set_status_file)) {
stop(
"The web stock-status tables are missing. Run ",
"`Rscript scripts/build-esc31-web-status-tables.R .`.",
call. = FALSE
)
}
reference_set_status <- readr::read_csv(
reference_set_status_file,
show_col_types = FALSE
)
reference_set_status |>
transmute(
`State year` = .data$State_year,
`MSY reference year` = .data$MSY_reference_year,
`Relative TRO` = format_interval(
.data$Relative_TRO_median,
.data$Relative_TRO_lower,
.data$Relative_TRO_upper,
3L
),
`P(relative TRO < 0.2)` = format_decimal(
.data$Probability_relative_TRO_below_0_2,
3L
),
`Relative B10+` = format_interval(
.data$Relative_B10_plus_median,
.data$Relative_B10_plus_lower,
.data$Relative_B10_plus_upper,
3L
),
`F/F_MSY` = format_interval(
.data$F_F_MSY_median,
.data$F_F_MSY_lower,
.data$F_F_MSY_upper,
3L
),
`TRO/TRO_MSY` = format_interval(
.data$TRO_TRO_MSY_median,
.data$TRO_TRO_MSY_lower,
.data$TRO_TRO_MSY_upper,
3L
),
`MSY (t)` = format_interval(
.data$MSY_t_median,
.data$MSY_t_lower,
.data$MSY_t_upper,
0L
)
) |>
kable(align = rep("r", 8L))| State year | MSY reference year | Relative TRO | P(relative TRO < 0.2) | Relative B10+ | F/F_MSY | TRO/TRO_MSY | MSY (t) |
|---|---|---|---|---|---|---|---|
| 2023 | 2023 | 0.255 (0.189-0.338) | 0.068 | 0.231 (0.171-0.310) | 0.567 (0.397-0.805) | 0.929 (0.590-1.495) | 32,984 (26,813-40,131) |
| 2025 | 2025 | 0.271 (0.201-0.359) | 0.022 | 0.236 (0.174-0.319) | 0.553 (0.383-0.811) | 0.985 (0.629-1.587) | 32,529 (26,400-39,381) |
| 2026 | 2025 | 0.281 (0.209-0.373) | 0.011 | 0.231 (0.171-0.313) | 0.553 (0.383-0.811) | 0.985 (0.629-1.587) | 32,529 (26,400-39,381) |
The complete numeric values are available as a machine-readable CSV, with source identities recorded in the shared stock-status provenance file.
if (base_posterior_msy_available) {
base_posterior_msy_year_summary |>
slice_tail(n = 5L) |>
transmute(
Year = as.character(year),
`MSY (t)` = format_interval(MSY_median, MSY_lower, MSY_upper, 0),
`TRO_MSY/TRO_0` = format_interval(
Bmsy_B0_median, Bmsy_B0_lower, Bmsy_B0_upper, 3
),
`Age-2+ biomass at MSY (t)` = format_interval(
TBmsy_median, TBmsy_lower, TBmsy_upper, 0
),
`F_MSY (ages 2-15)` = format_interval(
Fmsy_a215_median, Fmsy_a215_lower, Fmsy_a215_upper, 4
),
`Fitted TRO/TRO_MSY` = format_interval(
B_Bmsy_median, B_Bmsy_lower, B_Bmsy_upper, 3
),
`Fitted F/F_MSY` = format_interval(
F_Fmsy_median, F_Fmsy_lower, F_Fmsy_upper, 3
)
) |>
kable(
align = c("r", rep("r", 6)),
caption = "Posterior medians and equal-tailed 95% credible intervals for equilibrium MSY reference points in the five most recent fitted catch years."
)
} else {
base_msy_summary |>
slice_tail(n = 5L) |>
transmute(
Year = as.character(year),
`MSY (t)` = format_decimal(MSY, digits = 0),
`TRO_MSY/TRO_0` = format_decimal(Bmsy_B0, digits = 3),
`Age-2+ biomass at MSY (t)` = format_decimal(TBmsy, digits = 0),
`F_MSY (ages 2-15)` = format_decimal(Fmsy_a215, digits = 4),
`Fitted TRO/TRO_MSY` = format_decimal(B_Bmsy, digits = 3),
`Fitted F/F_MSY` = format_decimal(F_Fmsy, digits = 3)
) |>
kable(
align = c("r", rep("r", 6)),
caption = "MLE equilibrium MSY reference points for the five most recent fitted catch years."
)
}| Year | MSY (t) | TRO_MSY/TRO_0 | Age-2+ biomass at MSY (t) | F_MSY (ages 2-15) | Fitted TRO/TRO_MSY | Fitted F/F_MSY |
|---|---|---|---|---|---|---|
| 2021 | 32,985 (27,042-40,330) | 0.274 (0.272-0.276) | 465,626 (381,197-565,787) | 0.0923 (0.0878-0.0963) | 0.876 (0.688-1.107) | 0.565 (0.495-0.644) |
| 2022 | 31,743 (25,956-38,475) | 0.273 (0.271-0.276) | 459,205 (375,924-558,258) | 0.0920 (0.0878-0.0961) | 0.906 (0.711-1.145) | 0.575 (0.502-0.657) |
| 2023 | 32,901 (26,979-40,003) | 0.275 (0.272-0.278) | 466,840 (383,453-569,447) | 0.0915 (0.0864-0.0961) | 0.927 (0.731-1.175) | 0.566 (0.491-0.650) |
| 2024 | 32,864 (26,977-40,071) | 0.275 (0.272-0.278) | 467,496 (383,598-570,270) | 0.0909 (0.0859-0.0953) | 0.955 (0.751-1.204) | 0.570 (0.491-0.656) |
| 2025 | 32,451 (26,663-39,862) | 0.275 (0.272-0.278) | 466,470 (382,885-570,011) | 0.0920 (0.0868-0.0965) | 0.987 (0.779-1.237) | 0.553 (0.466-0.649) |
if (base_posterior_msy_available) {
ggplot(base_posterior_msy_year_summary, aes(year)) +
geom_ribbon(
aes(ymin = MSY_lower, ymax = MSY_upper),
fill = fit_expected_color,
alpha = 0.2
) +
geom_line(
aes(y = MSY_median, color = "Posterior MSY"),
linewidth = 0.75
) +
geom_line(
data = base_msy_summary,
aes(x = year, y = MSY),
color = "grey30",
linetype = "dotted",
linewidth = 0.6,
inherit.aes = FALSE
) +
geom_line(
aes(y = observed_catch, color = "Observed catch"),
linewidth = 0.7
) +
scale_color_manual(
values = c(
`Observed catch` = fit_observed_color,
`Posterior MSY` = fit_expected_color
)
) +
scale_y_continuous(
labels = label_comma(), limits = c(0, NA),
expand = expansion(mult = c(0, 0.05))
) +
labs(x = NULL, y = "Tonnes", color = NULL) +
theme(legend.position = "top")
} else {
base_msy_summary |>
select(
Year = year,
`Observed catch` = observed_catch,
`Equilibrium MSY` = MSY
) |>
pivot_longer(-Year, names_to = "Series", values_to = "Tonnes") |>
ggplot(aes(Year, Tonnes, color = Series)) +
geom_line(linewidth = 0.7) +
scale_color_manual(
values = c(
`Observed catch` = fit_observed_color,
`Equilibrium MSY` = fit_expected_color
)
) +
scale_y_continuous(
labels = label_comma(), limits = c(0, NA),
expand = expansion(mult = c(0, 0.05))
) +
labs(x = NULL, y = "Tonnes", color = NULL) +
theme(legend.position = "top")
}
msy_status_breaks <- c("TRO/TRO_MSY", "F/F_MSY")
msy_status_labels <- expression(TRO/TRO[MSY], F/F[MSY])
if (base_posterior_msy_available) {
posterior_status <- bind_rows(
base_posterior_msy_year_summary |>
transmute(
Year = year, Ratio = "TRO/TRO_MSY",
Median = B_Bmsy_median, Lower = B_Bmsy_lower, Upper = B_Bmsy_upper
),
base_posterior_msy_year_summary |>
transmute(
Year = year, Ratio = "F/F_MSY",
Median = F_Fmsy_median, Lower = F_Fmsy_lower, Upper = F_Fmsy_upper
)
)
mle_status <- base_msy_summary |>
select(
Year = year,
`TRO/TRO_MSY` = B_Bmsy,
`F/F_MSY` = F_Fmsy
) |>
pivot_longer(-Year, names_to = "Ratio", values_to = "Value")
ggplot(posterior_status, aes(Year, Median, color = Ratio, fill = Ratio)) +
geom_hline(yintercept = 1, linetype = "dashed", color = "grey35") +
geom_ribbon(
aes(ymin = Lower, ymax = Upper),
alpha = 0.16,
color = NA
) +
geom_line(linewidth = 0.75) +
geom_line(
data = mle_status,
aes(Year, Value, color = Ratio),
linetype = "dotted",
linewidth = 0.6,
inherit.aes = FALSE
) +
scale_color_manual(
values = c(
`TRO/TRO_MSY` = fit_expected_color,
`F/F_MSY` = fit_observed_color
),
breaks = msy_status_breaks,
labels = msy_status_labels
) +
scale_fill_manual(
values = c(
`TRO/TRO_MSY` = fit_expected_color,
`F/F_MSY` = fit_observed_color
),
breaks = msy_status_breaks,
labels = msy_status_labels
) +
scale_y_continuous(
limits = c(0, NA), expand = expansion(mult = c(0, 0.05))
) +
labs(x = NULL, y = "Ratio", color = NULL, fill = NULL) +
theme(legend.position = "top")
} else {
base_msy_summary |>
select(
Year = year,
`TRO/TRO_MSY` = B_Bmsy,
`F/F_MSY` = F_Fmsy
) |>
pivot_longer(-Year, names_to = "Ratio", values_to = "Value") |>
ggplot(aes(Year, Value, color = Ratio)) +
geom_hline(yintercept = 1, linetype = "dashed", color = "grey35") +
geom_line(linewidth = 0.7) +
scale_color_manual(
values = c(
`TRO/TRO_MSY` = fit_expected_color,
`F/F_MSY` = fit_observed_color
),
breaks = msy_status_breaks,
labels = msy_status_labels
) +
scale_y_continuous(
limits = c(0, NA), expand = expansion(mult = c(0, 0.05))
) +
labs(x = NULL, y = "Ratio", color = NULL) +
theme(legend.position = "top")
}
The appendix reports the saved posterior’s diagnostics. Failed gates are marked explicitly; a failed posterior remains available for inspection but is not used downstream. The sampler-parameter plot includes warmup and marks its boundary. The pairs plot uses retained draws and shows the five slowest parameters with trace plots on the diagonal. All 750 retained iterations in each of the four chains are included in the parameter diagnostics.
if (nrow(base_mcmc_gate_table)) {
base_mcmc_gate_table |>
mutate(
Result = cell_spec(
Result,
bold = TRUE,
color = ifelse(Result == "Pass", "#006100", "#9C0006"),
background = ifelse(Result == "Pass", "#C6EFCE", "#FFC7CE")
)
) |>
kable(align = c("l", "r", "r", "l", "c"), escape = FALSE)
} else {
knitr::asis_output("*No saved base MCMC was available for this render.*")
}| Gate | Requirement | Observed | Limiting variable | Result |
|---|---|---|---|---|
| MLE convergence | 0 | 0 | - | Pass |
| Finite MLE objective | Finite | 7,616.638 | - | Pass |
| Maximum MLE gradient | <= 0.01 | 0.000000 | - | Pass |
| Estimability | Estimable | Estimable | - | Pass |
| Maximum R-hat | < 1.01 | 1.0087 | par_log_sel_1[78] | Pass |
| Minimum bulk ESS | >= 400 | 1,157 | lp__ | Pass |
| Minimum tail ESS | >= 400 | 1,178 | par_rdev_y[24] | Pass |
| Divergences | 0 | 0 | - | Pass |
| Maximum-treedepth hits | 0 | 0 | - | Pass |
| Biological draws expected | 3000 | 3000 | - | Pass |
| Biological draws evaluated | 3000 | 3000 | - | Pass |
| Non-finite biological draws | 0 | 0 | - | Pass |
| Invalid biological draws | 0 | 0 | - | Pass |
| Maximum raw seasonal harvest | <= 0.9 | 0.878696 | chain 3, iteration 732, year 1960, season 1, age 10 | Pass |
| Minimum populated numbers-at-age | > 0 | 490.556 | chain 2, iteration 503, year 2011, season 2, age 29 | Pass |
| Minimum natural mortality | > 0 in every draw | 0.079444 | chain 4, iteration 319, age 25 | Pass |
| Maximum continuation penalty | <= 1e-8 | 0.000000 | - | Pass |
if (nrow(base_monitor_table)) {
base_monitor_table |>
mutate(
`R-hat` = format_decimal(`R-hat`, 4),
`Bulk ESS` = format_decimal(`Bulk ESS`, 0),
`Tail ESS` = format_decimal(`Tail ESS`, 0)
) |>
kable(align = c("l", "r", "r", "r"))
} else {
knitr::asis_output("*No parameter diagnostics are available.*")
}| Parameter | R-hat | Bulk ESS | Tail ESS |
|---|---|---|---|
| par_log_B0 | 1.0020 | 5,025 | 2,272 |
| par_log_m10 | 1.0010 | 5,074 | 2,540 |
| par_log_m30 | 0.9997 | 3,822 | 2,106 |
| par_log_cpue_q | 1.0005 | 5,365 | 1,849 |
| par_log_sel_1[78] | 1.0087 | 5,306 | 2,141 |
| par_rdev_y[24] | 1.0026 | 1,352 | 1,178 |
| lp__ | 1.0014 | 1,157 | 1,608 |
if (base_mcmc_available) {
mcmc <- as_tmbfit(base_fit)
SparseNUTS::plot_sampler_params(mcmc, plot = FALSE) +
geom_vline(
xintercept = base_fit$mcmc$warmup + 0.5,
linetype = "dashed",
color = "grey25"
)
}
if (base_mcmc_available) {
pairs_rtmb(fit = mcmc, order = "slow", pars = 1:5)
}
The length-based natural-mortality model is parsimonious and biologically interpretable. It estimates only \(M_{10}\) and \(M_{30}\), fixes \(m_c=-1\), and derives the younger-age mortality curve from observed length at age. This retains the broad shape of the previous assessment while avoiding four freely estimated mortality anchors. At the MLE, \(M_{10}\) is 0.1074 and \(M_{30}\) is 0.4577; the all-age minimum is 0.0893. The inverse-length part of the curve has a clear empirical basis (Lorenzen 1996, 2000, 2022). The two-unit profile ranges are 0.1004-0.1144 for M10 and 0.3730-0.5617 for M30.
The numerical behaviour is strong. The MLE has objective 7,616.638, maximum absolute gradient 2.6e-10, normal optimizer convergence, and an estimable Hessian for all 1,559 active effects. The biological constraints are comfortably inactive: maximum raw seasonal harvest is 0.618, the minimum populated number-at-age is 1,696.5, the continuation penalty is zero to reported precision, and the preventive wall contributes only 2.57e-21 objective units. Thus, the solution is not being held at an artificial harvest boundary.
The current and refitted previous-assessment total reproductive output trajectories tell a recognizably similar long-term story. At their latest overlapping year, 2023, relative TRO is 0.242 in the current fit and 0.281 in the comparison fit. This continuity is useful: the revised model changes the level and uncertainty without producing an unexplained reversal of the stock history.
The deliberate CPUE and composition down-weighting is also a defensible feature of this base, not a residual-calibration failure. The CPUE SDNR is 0.795 and the available composition SDNRs range from 0.544 to 0.773. Those values reflect the intended conservative influence of these data. Aerial-survey and conventional-tag uncertainty are separately calibrated and fixed, while the remaining scalar diagnostics provide useful independent checks.
The four-chain posterior passes the specified gates: maximum R-hat is 1.0087, minimum bulk ESS is 1,157, minimum tail ESS is 1,178, and there are no divergences or maximum-treedepth hits. These diagnostics meet the production criteria recommended for rank-normalized MCMC assessment (Vehtari et al. 2021) and the additional stock-assessment checks described by Monnahan (2024). All retained draws must also pass the full biological-state contract, so sampler diagnostics alone are not being treated as sufficient. For 2025, posterior median TRO/TRO_MSY is 0.987 (95% interval 0.779-1.237) and median F/F_MSY is 0.553 (0.466-0.649).
Fixing \(m_c=-1\) transfers uncertainty from estimation into model structure. The resulting curve is plausible and smooth, but the base deliberately does not estimate the exponent. The natural-mortality profiles show the conditional information for \(M_{10}\) and \(M_{30}\), not uncertainty about the functional form itself. The sensitivity page now separates two questions: one model estimates the length exponent while retaining the base curve, and the direct four-anchor model tests a different mortality parameterization. The former, not the four-anchor comparison, resolves the exponent-estimation reservation.
The purposeful down-weighting also has a cost: CPUE and composition data exert less leverage on recent abundance, recruitment, and selectivity. Their low SDNRs should not be tuned back to one, but temporal runs in CPUE or survey residuals and coherent age/length patterns in compositions still indicate possible structural mismatch. Similarly, choosing aerial and tag error scales to bring their SDNRs near one is transparent and practical, but those same SDNRs are no longer independent validation of the chosen weights.
The model remains high-dimensional and contains strong structural choices, including the one-block LL4 selectivity and fixed selectivity hyperparameters. All MCMC chains start at the same MLE to avoid biologically infeasible warmup states. Independent random streams, rank plots, effective sample sizes, and full-state checks mitigate that choice, but identical starts give a weaker test for undiscovered posterior modes than genuinely dispersed feasible starts.
The completed sensitivities qualify, but do not overturn, the base result. For the 2026 endpoint, posterior median relative TRO is 0.282 in the 3,000-draw base mid-cell posterior. The earlier reference-set table reports 0.281 from the separately identified balanced 2,000-draw nine-cell posterior. For comparison, the same endpoint median is 0.239 when CPUE catchability changes from 2008 and 0.339 under the direct four-anchor mortality alternative; these are the largest shifts among the accepted cases. Estimating the length-mortality exponent (0.287) and relaxing Indonesian selectivity (0.279) remain close to the base, while excluding both close-kin likelihoods gives 0.266 and is accepted only under the documented two-divergence NoPOPHSP exception. The sensitivity intervals overlap, but the spread confirms that CPUE catchability and mortality parameterization remain material structural uncertainties.
Finally, the annual MSY calculations are equilibrium reference points based on each year’s selectivity, weights, and fitted catch allocation. They are not a dynamic projection, TAC recommendation, or substitute for the projection analysis. The sensitivity set, all nine MCMC grid cells, the balanced 2,000-draw grid posterior, and the direct-get_M 108-fit MLE grid and resample are accepted. Status close to a reference-point boundary should be interpreted from the full credible interval and sensitivities, not from the MLE or posterior median alone.
On balance, I consider this a strong and viable ESC31 base model. Its mortality curve is simpler and more biologically coherent than the superseded direct-log anchor model, its MLE is numerically clean and biologically interior, and its long-term trajectory remains explainable relative to the previous assessment. The important caveats are explicit rather than hidden: fixed mortality shape, intentional CPUE and composition down-weighting, calibrated survey/tag weights, and the demonstrated sensitivity to the CPUE-2008 and direct-mortality alternatives. The accepted sensitivity set and nine-cell grid support using it as the conditioning base. The grid products have now been built and checked. The four 2,000-draw projection arms were completed and reviewed under the preceding contract. Their exact fingerprint-bound format-8 outputs remain accepted for reporting; format 9 applies to future projection runs.