---
title: "OMMP16 model fit review"
bibliography: references.bib
format:
html:
embed-resources: true
toc: true
toc-depth: 3
toc-location: left
number-sections: true
lightbox: true
theme: default
mainfont: system-ui
code-fold: true
code-tools: true
favicon: ../images/favicon.svg
include-in-header:
text: |
<link rel="icon" href="../images/favicon.svg" type="image/svg+xml">
include-before-body: ../includes/embedded-lightbox.html
execute:
warning: false
message: false
editor_options:
chunk_output_type: console
---
# Overview
::: {.callout-note}
## Historical checkpoint
This page records the model-fit review undertaken around OMMP16. It is not the
final ESC31 assessment report: later scientific decisions, accepted MCMCs, the
nine-cell grid, and the 2,000-draw projections are documented in the
[ESC31 workflow](../ESC31/1_intro.html). In particular, the selected ESC31
base combines the 20% baseline CPUE CV with the raw year-specific GAM22
estimation CV; it does not use the scaled-CV sensitivity described below.
:::
This report presents the OMMP16 operating model fit prepared from the 2026 input
files. It includes the maximum-likelihood model setup, estimated parameters,
selectivity hyper-parameters, model-fit plots, residual diagnostics, and
population-dynamics plots.
The CCSBT data page describes the Commission data holdings used to
support assessment inputs, including catch, catch-and-effort, catch-at-size, and
related monitoring data [@CCSBTData2026]. The latest CCSBT ESC data-input
papers found at the time of drafting were the 2025 ESC30 updates, including
CSIRO close-kin mark-recapture and gene-tagging program updates
[@Farley2025CKMR; @Preece2025GT], the GAM22 CPUE input index
[@ItohTakahasi2025GAM22], and fisheries indicators [@Patterson2025Indicators].
```{r}
#| label: setup
#| include: false
knitr::knit_meta_add(list(rmarkdown::html_dependency_jquery()))
library(tidyverse)
library(sbt)
library(RTMB)
library(patchwork)
library(kableExtra)
library(ggridges)
library(scales)
theme_set(theme_bw())
fig_width <- 9
fig_height <- 6
`%||%` <- function(x, y) {
if (is.null(x) || length(x) == 0L) y else x
}
format_numeric <- function(x, digits = 4) {
ifelse(is.finite(x), signif(x, digits), x)
}
format_decimal <- function(x, digits = 3) {
out <- rep(NA_character_, length(x))
finite <- is.finite(x)
x[finite & abs(x) < 0.5 * 10^(-digits)] <- 0
out[finite] <- formatC(x[finite], format = "f", digits = digits, big.mark = ",")
out[!finite & !is.na(x)] <- as.character(x[!finite & !is.na(x)])
out
}
format_table <- function(x, digits = 4) {
x |>
mutate(across(where(is.numeric), ~ format_numeric(.x, digits = digits)))
}
format_decimal_table <- function(x, digits = 3) {
x |>
mutate(across(where(is.numeric), ~ format_decimal(.x, digits = digits)))
}
replace_missing <- function(x, missing = "–") {
x <- as.data.frame(x)
x[] <- lapply(x, function(column) {
column <- as.character(column)
column[is.na(column) | column == "NA"] <- missing
column
})
as_tibble(x)
}
format_residual_source <- function(x) {
x <- sub("^.*[:][:]", "", x)
sub("\\(\\)$", "", x)
}
parameter_labels <- function(x) {
recode(
as.character(x),
B0 = "B₀",
M0 = "M₀",
M4 = "M₄",
M10 = "M₁₀",
M30 = "M₃₀",
.default = as.character(x)
)
}
numeric_table_align <- function(x) {
ifelse(vapply(x, is.numeric, logical(1)), "r", "l")
}
format_residual_table <- function(x, digits = 3) {
out <- x
numeric_columns <- names(out)[vapply(out, is.numeric, logical(1))]
for (column in setdiff(numeric_columns, "N")) {
out[[column]] <- format_decimal(out[[column]], digits = digits)
}
if ("N" %in% names(out)) {
out$N <- ifelse(is.na(out$N), NA_character_, as.character(out$N))
}
out
}
zero_y <- function(p) {
p +
scale_y_continuous(limits = c(0, NA), expand = expansion(mult = c(0, 0.05)))
}
pad_limits <- function(x, pad = 0.25) {
x <- x[is.finite(x)]
if (!length(x)) return(c(NA_real_, NA_real_))
range(x) + c(-pad, pad)
}
year_limits <- function(x) {
pad_limits(x, pad = 0.25)
}
age_limits <- function(data) {
c(data$min_age - 0.25, data$max_age + 0.25)
}
bin_limits <- function(x, step = NULL) {
x <- sort(unique(x[is.finite(x)]))
if (!length(x)) return(c(NA_real_, NA_real_))
if (is.null(step)) {
dx <- diff(x)
step <- if (length(dx)) min(dx[dx > 0], na.rm = TRUE) else 1
}
range(x) + c(-0.5, 0.5) * step
}
selectivity_fig_height <- function(years) {
n_year <- length(unique(years[is.finite(years)]))
min(32, max(10, 6 + 0.85 * n_year))
}
fit_observed_color <- "#D55E00"
fit_expected_color <- "#0072B2"
get_report_aerial_cov <- function(model_data) {
if (!is.null(model_data$aerial_cov)) return(model_data$aerial_cov)
aerial_env <- new.env(parent = emptyenv())
get("data", envir = asNamespace("utils"))("aerial_cov", package = "sbt", envir = aerial_env)
aerial_env$aerial_cov
}
plot_index_ci <- function(df, y_label) {
ggplot(df, aes(x = .data$year, y = .data$obs)) +
geom_linerange(
aes(ymin = .data$lower, ymax = .data$upper, color = "Observed"),
linewidth = 0.35,
alpha = 0.65
) +
geom_point(aes(color = "Observed"), size = 1.8) +
geom_line(aes(y = .data$pred, color = "Predicted"), linewidth = 0.6) +
labs(x = "Year", y = y_label, color = NULL) +
scale_color_manual(values = c("Observed" = fit_observed_color, "Predicted" = fit_expected_color)) +
scale_x_continuous(limits = year_limits(df$year), breaks = pretty_breaks()) +
scale_y_continuous(limits = c(0, NA), expand = expansion(mult = c(0, 0.05)))
}
plot_aerial_fit_ci <- function(data, object) {
par_list <- object$env$parList(object$env$last.par.best)
aerial_tau <- exp(par_list$par_log_aerial_tau)
aerial_cov <- get_report_aerial_cov(data)
rep <- object$report(object$env$last.par.best)
df <- data$aerial_survey |>
transmute(
year = .data$Year,
obs = .data$Index,
pred = rep$aerial_pred,
sigma = sqrt(diag(aerial_cov)[seq_len(n())] + aerial_tau^2)
) |>
mutate(
lower = ifelse(.data$obs > 0, exp(log(.data$obs) - 1.96 * .data$sigma), NA_real_),
upper = ifelse(.data$obs > 0, exp(log(.data$obs) + 1.96 * .data$sigma), NA_real_)
)
plot_index_ci(df, "Aerial survey")
}
plot_catch_input_output <- function(data, object) {
fsh <- c("LL1", "LL2", "LL3", "LL4", "Indonesia", "Australia")
yrs1 <- data$first_yr:data$last_yr
yrs2 <- data$first_yr_catch:data$last_yr
table_year <- function(x, years) {
labels <- as.character(x)
label_years <- suppressWarnings(as.integer(labels))
if (all(is.finite(label_years)) && all(label_years %in% years)) label_years else years[as.integer(x)]
}
df_obs <- as.data.frame.table(data$catch_obs_ysf, responseName = "obs")
names(df_obs)[seq_len(3)] <- c("Var1", "Var2", "Var3")
df_obs <- df_obs |>
filter(.data$obs > 0) |>
transmute(
Year = table_year(.data$Var1, yrs2),
Season = paste("Season:", as.integer(.data$Var2)),
Fishery = fsh[as.integer(.data$Var3)],
obs = .data$obs
)
df_pred <- as.data.frame.table(
object$report(object$env$last.par.best)$catch_pred_ysf,
responseName = "pred"
)
names(df_pred)[seq_len(3)] <- c("Var1", "Var2", "Var3")
df_pred <- df_pred |>
transmute(
Year = table_year(.data$Var1, yrs1),
Season = paste("Season:", as.integer(.data$Var2)),
Fishery = fsh[as.integer(.data$Var3)],
pred = .data$pred
) |>
right_join(df_obs, by = join_by("Year", "Season", "Fishery")) |>
mutate(Fishery = factor(.data$Fishery, levels = fsh))
ggplot(df_pred, aes(x = .data$Year)) +
geom_point(aes(y = .data$obs / 1000), color = fit_observed_color, size = 1.8) +
geom_line(aes(y = .data$pred / 1000, group = interaction(.data$Fishery, .data$Season)),
color = fit_expected_color, linewidth = 0.6
) +
facet_wrap(Fishery ~ Season, scales = "free_y") +
labs(x = "Year", y = "Catch (thousands of tonnes)") +
scale_x_continuous(limits = year_limits(c(data$first_yr, data$last_yr)), breaks = pretty_breaks()) +
scale_y_continuous(limits = c(0, NA), expand = expansion(mult = c(0, 0.05))) +
theme(legend.position = "none")
}
population_plot_color <- "#00BFC4"
plot_initial_numbers_report <- function(data, object) {
rep <- object$report(object$env$last.par.best)
tibble(
age = data$min_age:data$max_age,
value = as.numeric(rep$number_ysa[1, 1, ]) / 1e6
) |>
ggplot(aes(x = .data$age, y = .data$value)) +
geom_line(color = population_plot_color, linewidth = 0.6) +
geom_point(color = population_plot_color, size = 1.8) +
labs(x = "Age", y = "Initial numbers (millions)") +
scale_x_continuous(limits = age_limits(data), breaks = pretty_breaks()) +
scale_y_continuous(limits = c(0, NA), expand = expansion(mult = c(0, 0.05)))
}
plot_natural_mortality_report <- function(data, object) {
rep <- object$report(object$env$last.par.best)
mortality_points <- tibble(
age = data$min_age:data$max_age,
value = as.numeric(rep$M_a),
Parameter = case_when(
.data$age == 0 ~ "M₀",
.data$age == 4 ~ "M₄",
.data$age == 10 ~ "M₁₀",
.data$age == 30 ~ "M₃₀",
TRUE ~ NA_character_
)
)
ggplot(mortality_points, aes(x = .data$age, y = .data$value)) +
geom_line(color = population_plot_color, linewidth = 0.6) +
geom_point(color = population_plot_color, size = 1.8) +
geom_point(
data = filter(mortality_points, !is.na(.data$Parameter)),
aes(color = .data$Parameter),
size = 3
) +
labs(x = "Age", y = "Natural mortality", color = "Parameter") +
scale_color_manual(
values = c("M₀" = "#D55E00", "M₄" = "#009E73", "M₁₀" = "#CC79A7", "M₃₀" = "#E69F00"),
breaks = c("M₀", "M₄", "M₁₀", "M₃₀")
) +
scale_x_continuous(limits = age_limits(data), breaks = pretty_breaks()) +
scale_y_continuous(limits = c(0, NA), expand = expansion(mult = c(0, 0.05)))
}
plot_rec_devs_report <- function(data, object) {
years <- data$first_yr:data$last_yr
rec_devs <- object$env$last.par.best[names(object$par) %in% "par_rdev_y"]
df <- tibble(
year = years,
value = as.numeric(rec_devs),
prior = ifelse(seq_along(years) > length(years) - 3, "AR1 prior", "Independent normal prior")
)
ggplot(df, aes(x = .data$year, y = .data$value)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black", linewidth = 0.4) +
geom_line(color = population_plot_color, linewidth = 0.45) +
geom_point(aes(color = .data$prior), size = 1.7) +
labs(x = "Year", y = "Recruitment deviate", color = "Recruitment-deviation prior") +
scale_color_manual(values = c("Independent normal prior" = population_plot_color, "AR1 prior" = fit_observed_color)) +
scale_x_continuous(limits = year_limits(df$year), breaks = pretty_breaks()) +
theme(legend.position = "none")
}
plot_recruitment_report <- function(data, object) {
rep <- object$report(object$env$last.par.best)
df <- tibble(
year = data$first_yr:(data$last_yr + 1),
value = as.numeric(rep$number_ysa[, 1, 1]) / 1e6,
prior = ifelse(seq_along(year) > length(year) - 3, "AR1 prior", "Independent normal prior")
)
df |>
ggplot(aes(x = .data$year, y = .data$value)) +
geom_hline(yintercept = rep$R0 / 1e6, color = "black", linetype = "dashed", linewidth = 0.4) +
geom_line(color = population_plot_color, linewidth = 0.6) +
geom_point(aes(color = .data$prior), size = 1.3) +
labs(x = "Year", y = "Recruitment (millions)") +
scale_color_manual(values = c("Independent normal prior" = population_plot_color, "AR1 prior" = fit_observed_color)) +
scale_x_continuous(limits = year_limits(df$year), breaks = pretty_breaks()) +
scale_y_continuous(limits = c(0, NA), expand = expansion(mult = c(0, 0.05))) +
theme(legend.position = "none")
}
plot_total_genetics_report <- function(data, object, hsp_plot_list) {
if (inherits(hsp_plot_list, "error")) {
stop(conditionMessage(hsp_plot_list), call. = FALSE)
}
if (is.null(hsp_plot_list$total_matches) || is.null(hsp_plot_list$total_matches$data)) {
stop("Total POP/HSP match panel not found.", call. = FALSE)
}
df <- as_tibble(hsp_plot_list$total_matches$data) |>
transmute(
dtype = as.character(.data$dtype),
obs = .data$obs,
med = .data$med,
lq = .data$lq,
uq = .data$uq
)
rep <- object$report(object$env$last.par.best)
if (!is.null(rep$gt_prob) && !is.null(data$gt_obs) && nrow(data$gt_obs)) {
gt_prob <- pmin(pmax(as.numeric(rep$gt_prob), 0), 1)
gt_med <- sum(data$gt_obs$Nsam * gt_prob, na.rm = TRUE)
gt_sd <- sqrt(sum(data$gt_obs$Nsam * gt_prob * (1 - gt_prob), na.rm = TRUE))
df <- bind_rows(
df,
tibble(
dtype = "GT",
obs = sum(data$gt_obs$Nmatch, na.rm = TRUE),
med = gt_med,
lq = max(0, gt_med - 1.96 * gt_sd),
uq = gt_med + 1.96 * gt_sd
)
)
}
df <- df |> mutate(dtype = factor(.data$dtype, levels = c("GT", "POPs", "HSPs")))
ggplot(df, aes(x = .data$dtype)) +
geom_linerange(aes(ymin = .data$lq, ymax = .data$uq), color = fit_expected_color, linewidth = 0.45) +
geom_point(aes(y = .data$med), color = fit_expected_color, size = 2.2) +
geom_point(aes(y = .data$obs), color = fit_observed_color, size = 2.2) +
labs(x = "Data set", y = "Matches") +
scale_y_continuous(limits = c(0, NA), expand = expansion(mult = c(0, 0.05)))
}
harmonize_ck_colors <- function(p) {
if (!inherits(p, "ggplot")) return(p)
for (i in seq_along(p$layers)) {
layer_color <- p$layers[[i]]$aes_params$colour %||%
p$layers[[i]]$aes_params$color %||%
p$layers[[i]]$aes_params$col %||%
NA_character_
if (identical(layer_color, "dark salmon")) {
p$layers[[i]]$aes_params$colour <- fit_observed_color
p$layers[[i]]$aes_params$color <- fit_observed_color
p$layers[[i]]$aes_params$col <- fit_observed_color
}
if (identical(layer_color, "dark orchid 4")) {
p$layers[[i]]$aes_params$colour <- fit_expected_color
p$layers[[i]]$aes_params$color <- fit_expected_color
p$layers[[i]]$aes_params$col <- fit_expected_color
}
}
p
}
remove_expected_points <- function(p) {
if (!inherits(p, "ggplot")) return(p)
keep <- vapply(p$layers, function(layer) {
geom_class <- class(layer$geom)[1]
layer_color <- layer$aes_params$colour %||%
layer$aes_params$color %||%
layer$aes_params$col %||%
NA_character_
!(geom_class == "GeomPoint" && identical(layer_color, fit_expected_color))
}, logical(1))
p$layers <- p$layers[keep]
p
}
harmonize_composition_fit_colors <- function(p) {
if (!inherits(p, "ggplot")) return(p)
for (i in seq_along(p$layers)) {
geom_class <- class(p$layers[[i]]$geom)[1]
if (geom_class == "GeomPoint") {
p$layers[[i]]$aes_params$colour <- fit_observed_color
p$layers[[i]]$aes_params$color <- fit_observed_color
p$layers[[i]]$aes_params$col <- fit_observed_color
}
if (geom_class %in% c("GeomLine", "GeomPath", "GeomStep")) {
p$layers[[i]]$aes_params$colour <- fit_expected_color
p$layers[[i]]$aes_params$color <- fit_expected_color
p$layers[[i]]$aes_params$col <- fit_expected_color
}
if (geom_class %in% c("GeomRibbon", "GeomArea")) {
p$layers[[i]]$aes_params$fill <- fit_expected_color
}
}
p
}
plot_removal_selectivity <- function(data, object, fishery_name = "LL4") {
fsh <- c("LL1", "LL2", "LL3", "LL4", "Indonesian", "Australian", "CPUE")
ages <- data$min_age:data$max_age
yrs <- data$first_yr:data$last_yr
rem <- c(data$removal_switch_f, 0)
fyr <- c(data$first_yr_catch_f, 1969)
df <- as.data.frame.table(
object$report(object$env$last.par.best)$hrate_fya,
responseName = "value"
) |>
mutate(
fishery = fsh[as.integer(.data$Var1)],
removal = rem[as.integer(.data$Var1)],
first_yr_catch = fyr[as.integer(.data$Var1)],
year = yrs[as.integer(.data$Var2)],
age = ages[as.integer(.data$Var3)]
) |>
filter(
.data$fishery == fishery_name,
.data$removal == 1,
.data$year >= .data$first_yr_catch
) |>
group_by(.data$fishery, .data$year) |>
mutate(
total = sum(.data$value, na.rm = TRUE),
value = ifelse(.data$total > 0, .data$value / .data$total, NA_real_)
) |>
ungroup() |>
filter(is.finite(.data$value), .data$value > 0)
if (!nrow(df)) {
return(blank_plot(
paste("Selectivity:", fishery_name),
paste(fishery_name, "has no removal-at-age values to display.")
))
}
ggplot(df, aes(x = .data$age, y = .data$year, height = .data$value, group = .data$year)) +
geom_density_ridges(stat = "identity", alpha = 0.55, rel_min_height = 0) +
facet_wrap(~fishery) +
labs(x = "Age", y = "Year") +
scale_x_continuous(limits = c(0, NA), expand = expansion(mult = c(0, 0.05))) +
scale_y_reverse(breaks = pretty_breaks())
}
blank_plot <- function(title, message) {
ggplot() +
annotate("text", x = 0, y = 0, label = str_wrap(message, width = 80), hjust = 0) +
xlim(0, 1) +
ylim(-0.2, 0.2) +
labs(title = title) +
theme_void()
}
make_plot <- function(label, code) {
output <- character()
plot <- tryCatch(
{
output <- capture.output(p <- code())
if (inherits(p, "ggplot")) invisible(ggplot_build(p))
p
},
error = function(e) {
structure(
blank_plot(label, conditionMessage(e)),
residual_error = conditionMessage(e)
)
}
)
list(label = label, plot = plot, output = output, error = attr(plot, "residual_error"))
}
read_cpue_time_varying_cv <- function(cpue_base, cpue_cv, cv_column = "cv_scaled") {
if (!all(c("year", cv_column) %in% names(cpue_cv))) {
stop("CPUE CV file must contain columns `year` and `", cv_column, "`.", call. = FALSE)
}
cpue_cv <- cpue_cv |>
transmute(Year = .data$year, CV = .data[[cv_column]])
missing_cv <- anti_join(select(cpue_base, Year), cpue_cv, by = "Year")
extra_cv <- anti_join(cpue_cv, select(cpue_base, Year), by = "Year")
if (nrow(missing_cv) || nrow(extra_cv)) {
stop(
"CPUE CV years do not align with cpue.csv. Missing years: ",
paste(missing_cv$Year, collapse = ", "),
". Extra years: ",
paste(extra_cv$Year, collapse = ", "),
".",
call. = FALSE
)
}
cpue_sens <- cpue_base |>
left_join(cpue_cv, by = "Year") |>
mutate(CV = ifelse(.data$CV > 0 & is.finite(.data$CV), .data$CV, NA_real_))
if (any(is.na(cpue_sens$CV))) {
stop("CPUE CV values must be positive and finite.", call. = FALSE)
}
cpue_sens
}
summarize_fit <- function(label, data_i, object_i, opt_i = NULL) {
par_list <- object_i$env$parList(object_i$env$last.par.best)
gradient <- tryCatch(object_i$gr(object_i$env$last.par.best), error = function(e) NA_real_)
max_gradient <- if (all(is.na(gradient))) NA_real_ else max(abs(gradient), na.rm = TRUE)
tibble(
Model = label,
`Convergence code` = opt_i$convergence %||% NA_integer_,
NLL = opt_i$objective %||% object_i$fn(object_i$env$last.par.best),
`Max gradient` = max_gradient,
B0 = exp(par_list$par_log_B0),
M0 = exp(par_list$par_log_m0),
M4 = exp(par_list$par_log_m4),
M10 = exp(par_list$par_log_m10),
M30 = exp(par_list$par_log_m30),
h = exp(par_list$par_log_h),
psi = exp(par_list$par_log_psi)
)
}
has_sensitivity_fit <- function(x) {
is.list(x) && !is.null(x$data) && !is.null(x$obj) && !is.null(x$opt)
}
restore_sensitivity_fit <- function(x) {
if (!has_sensitivity_fit(x)) return(x)
object_ok <- tryCatch(is.finite(x$obj$fn(x$opt$par)), error = function(e) FALSE)
if (!object_ok) {
x$obj <- MakeADFun(func = cmb(sbt_model, x$data), parameters = x$parameters, map = x$map)
}
x$obj$env$last.par.best <- x$opt$par
x
}
plot_pops_residuals_report <- function(data, obj) {
osa_res <- oneStepPredict(
obj = obj,
observation.name = "pop_nP",
method = "cdf",
trace = FALSE,
discrete = TRUE
)
df <- data.frame(data$pop_obs) |>
mutate(
ReleaseYear = .data$Cohort + data$first_yr - 1L,
CaptureYear = .data$CaptureYear + data$first_yr - 1L,
CaptureType = factor(.data$CaptureSwitch, levels = c(0, 1), labels = c("Direct age", "Length only")),
resid = as.numeric(osa_res$residual),
Sign = ifelse(.data$resid >= 0, "Positive", "Negative")
)
finite_resid <- df$resid[is.finite(df$resid)]
s <- sdnr(finite_resid)
subtitle <- sprintf(
"SDNR %.2f (95%% CI %.2f-%.2f); MAR %.2f",
s$SDNR, s$LCI, s$HCI, mar(finite_resid)
)
plot_df <- df |>
filter(is.finite(.data$resid)) |>
group_by(.data$CaptureType, .data$ReleaseYear, .data$CaptureYear) |>
summarise(
mean_resid = mean(.data$resid),
n = n(),
.groups = "drop"
)
p <- ggplot(plot_df, aes(x = .data$CaptureYear, y = .data$ReleaseYear, fill = .data$mean_resid)) +
geom_tile(color = "white", linewidth = 0.25) +
facet_wrap(~CaptureType, ncol = 1) +
scale_fill_gradient2(
low = "#2166AC",
mid = "white",
high = "#B2182B",
midpoint = 0,
limits = c(-3, 3),
oob = squish,
name = "Mean OSA residual"
) +
scale_x_continuous(breaks = pretty_breaks()) +
scale_y_continuous(breaks = pretty_breaks()) +
labs(
x = "Adult capture year",
y = "Juvenile release cohort",
subtitle = subtitle,
caption = "One-step-ahead (OSA) residuals for binomial POP observations from oneStepPredict on pop_nP, averaged within release-cohort and adult-capture-year cells."
)
attr(p, "residual_data") <- df
attr(p, "plot_data") <- df |> filter(is.finite(.data$resid))
attr(p, "residual_source") <- "oneStepPredict"
p
}
residual_data <- function(plot) {
dat <- attr(plot, "plot_data")
if (is.null(dat)) dat <- attr(plot, "residual_data")
if (is.null(dat) && inherits(plot, "ggplot")) dat <- plot$data
if (is.null(dat) || inherits(dat, "waiver")) return(NULL)
as.data.frame(dat)
}
composition_stat_row <- function(plot) {
summary <- attr(plot, "composition_summary")
if (is.null(summary) || !is.data.frame(summary) || !nrow(summary)) {
return(tibble(
`Francis SDNR` = NA_real_,
`Francis N multiplier` = NA_real_,
`McAllister-Ianelli harmonic N` = NA_real_,
`McAllister-Ianelli median N` = NA_real_,
`Mean input N` = NA_real_,
`McAllister-Ianelli/Input N` = NA_real_
))
}
tibble(
`Francis SDNR` = summary$francis_sdnr[1],
`Francis N multiplier` = summary$francis_n_multiplier[1],
`McAllister-Ianelli harmonic N` = summary$harmonic_eff_n[1],
`McAllister-Ianelli median N` = summary$median_eff_n[1],
`Mean input N` = summary$mean_input_n[1],
`McAllister-Ianelli/Input N` = summary$harmonic_eff_n_over_mean_input_n[1]
)
}
residual_stat_row <- function(label, plot, error = NULL) {
comp <- composition_stat_row(plot)
if (!is.null(error)) {
return(bind_cols(tibble(
`Data type` = label,
N = NA_integer_,
SDNR = NA_real_,
`SDNR lower 95%` = NA_real_,
`SDNR upper 95%` = NA_real_,
MAR = NA_real_,
Source = NA_character_
), comp))
}
dat <- residual_data(plot)
if (is.null(dat)) {
return(bind_cols(tibble(
`Data type` = label,
N = NA_integer_,
SDNR = NA_real_,
`SDNR lower 95%` = NA_real_,
`SDNR upper 95%` = NA_real_,
MAR = NA_real_,
Source = NA_character_
), comp))
}
residual_column <- intersect(c("resid", "residual"), names(dat))[1]
if (is.na(residual_column)) {
return(bind_cols(tibble(
`Data type` = label,
N = NA_integer_,
SDNR = NA_real_,
`SDNR lower 95%` = NA_real_,
`SDNR upper 95%` = NA_real_,
MAR = NA_real_,
Source = attr(plot, "residual_source") %||% NA_character_
), comp))
}
x <- dat[[residual_column]]
finite_x <- x[is.finite(x)]
if (!length(finite_x)) {
return(bind_cols(tibble(
`Data type` = label,
N = 0L,
SDNR = NA_real_,
`SDNR lower 95%` = NA_real_,
`SDNR upper 95%` = NA_real_,
MAR = NA_real_,
Source = attr(plot, "residual_source") %||% NA_character_
), comp))
}
s <- sdnr(finite_x)
bind_cols(tibble(
`Data type` = label,
N = length(finite_x),
SDNR = s$SDNR,
`SDNR lower 95%` = s$LCI,
`SDNR upper 95%` = s$HCI,
MAR = mar(finite_x),
Source = attr(plot, "residual_source") %||% "oneStepPredict"
), comp)
}
```
# Input data
```{r}
#| label: read-inputs
data_loc <- file.path("..", "ESC31", "csv_2026")
length_mean <- read_csv(file.path(data_loc, "mean_length.csv"), show_col_types = FALSE)
length_sd <- read_csv(file.path(data_loc, "sd_length.csv"), show_col_types = FALSE)
catch <- read_csv(file.path(data_loc, "catch.csv"), show_col_types = FALSE)
catch_UA <- read_csv(file.path(data_loc, "catch_UA.csv"), show_col_types = FALSE)
scenarios_surface <- read_csv(file.path(data_loc, "scenarios_surface.csv"), show_col_types = FALSE)
scenarios_LL1 <- read_csv(file.path(data_loc, "scenarios_LL1.csv"), show_col_types = FALSE)
POPs <- read_csv(file.path(data_loc, "POPs.csv"), show_col_types = FALSE)
HSPs <- read_csv(file.path(data_loc, "HSPs.csv"), show_col_types = FALSE)
GTs <- read_csv(file.path(data_loc, "GTs.csv"), show_col_types = FALSE)
troll <- read_csv(file.path(data_loc, "trolling_index.csv"), show_col_types = FALSE)
cpue <- read_csv(file.path(data_loc, "cpue.csv"), show_col_types = FALSE)
cpue_cv <- read_csv(file.path(data_loc, "CV_out_10.csv"), show_col_types = FALSE)
age_freq <- read_csv(file.path(data_loc, "age_freq.csv"), show_col_types = FALSE)
length_freq <- read_csv(file.path(data_loc, "lf_assessment.csv"), show_col_types = FALSE)
```
```{r}
#| label: create-data
data_in <- list(
last_yr = 2025,
age_increase_M = 25,
length_m50 = 150,
length_m95 = 180,
catch_surf_case = 1,
catch_LL1_case = 1,
length_mean = length_mean,
length_sd = length_sd,
catch = catch,
catch_UA = catch_UA,
scenarios_surf = scenarios_surface,
scenarios_LL1 = scenarios_LL1,
POPs = POPs,
HSPs = HSPs,
GTs = GTs,
troll = troll,
cpue = cpue,
age_freq = age_freq,
length_freq = length_freq,
removal_switch_f = c(0, 0, 0, 1, 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, 1957, 1961, 1965, 1969, 1973, 1977, 1981, 1985, 1989, 1993, 1997, 2001, 2006, 2007, 2008, 2011, 2014, 2017, 2020, 2023),
sel_LL2_yrs = c(1969, 2001, 2005, 2008, 2011, 2014, 2017, 2020, 2023),
sel_LL3_yrs = c(1954, 1961, 1965, 1969, 1970, 1971, 2005, 2006, 2007),
sel_LL4_yrs = c(1953),
sel_Ind_yrs = c(1976, 1997, 1999, 2002, 2004, 2006, 2008, 2010, 2012:2021),
sel_Aus_yrs = c(1952, 1969, 1973, 1977, 1981, 1985, 1989, 1993, 1997:2025),
sel_CPUE_yrs = c(1969, 1973, 1977, 1981, 1985, 1989, 1993, 1997, 2001, 2006, 2007, 2008, 2011, 2014, 2017, 2020, 2023),
af_switch = 1,
lf_switch = 1,
lf_minbin = c(1, 1, 1, 11, 6),
cpue_switch = 1,
cpue_a1 = 5,
cpue_a2 = 17,
aerial_switch = 4,
aerial_tau = 0.3,
troll_switch = 0,
pop_switch = 1,
hsp_switch = 1,
hsp_false_negative = 0.6840729,
gt_switch = 1,
tag_switch = 1,
tag_var_factor = 1.82
)
data <- get_data(data_in = data_in)
```
```{r}
#| label: tbl-input-summary
input_summary <- tibble(
`Input` = c("Assessment years", "Fisheries", "CPUE observations", "CPUE length-frequency rows", "Length-frequency rows", "Age-frequency rows", "Gene-tagging rows", "HSP rows", "POP rows", "Tag cohorts"),
`Value` = c(
paste(data$first_yr, data$last_yr, sep = "-"),
data$n_fishery,
length(data$cpue_obs),
ifelse(is.null(data$cpue_lfs), 0, nrow(data$cpue_lfs)),
nrow(data$lf_obs),
nrow(data$af_obs),
nrow(data$gt_obs),
nrow(data$hsp_obs),
nrow(data$pop_obs),
data$n_K
)
)
kable(
replace_missing(input_summary),
caption = "Summary of model-ready OMMP16 input data."
)
```
The plots below show the input data before fitting. For time-varying biological
inputs and composition data, the figures focus on the most recent available
year and two earlier recent years to keep the review compact.
```{r}
#| label: prepare-input-plots
#| include: false
fishery_lookup <- c(
"1" = "LL1",
"2" = "LL2",
"3" = "LL3",
"4" = "LL4",
"5" = "Indonesia",
"6" = "Australia",
"7" = "CPUE"
)
recent_years <- function(year, n = 3) {
sort(tail(sort(unique(year[!is.na(year)])), n))
}
growth_plot_years <- sort(unique(c(
data$first_yr,
round(mean(c(data$first_yr, data$last_yr))),
data$last_yr
)))
catch_plot_data <- catch |>
pivot_longer(-.data$Year, names_to = "Fishery", values_to = "Catch")
index_plot_data <- bind_rows(
cpue |>
transmute(Year = .data$Year, Series = "CPUE", Value = .data$CPUE),
data$aerial_survey |>
transmute(Year = .data$Year, Series = "Aerial survey", Value = .data$Index),
troll |>
transmute(Year = .data$Year, Series = "Troll", Value = .data$Median),
GTs |>
transmute(Year = .data$RecYear, Series = "Gene tagging", Value = .data$Nmatch)
)
length_bin_columns <- setdiff(names(length_freq), c("Fishery", "Year", "N"))
length_snapshot <- length_freq |>
mutate(
Fishery = recode(as.character(.data$Fishery), !!!fishery_lookup),
Year = as.integer(.data$Year)
) |>
filter(.data$Year %in% recent_years(.data$Year)) |>
pivot_longer(all_of(length_bin_columns), names_to = "Length", values_to = "Proportion") |>
mutate(
Length = as.numeric(.data$Length),
Year = factor(.data$Year)
) |>
group_by(.data$Fishery, .data$Year, .data$Length) |>
summarise(Proportion = sum(.data$Proportion, na.rm = TRUE), .groups = "drop") |>
group_by(.data$Fishery, .data$Year) |>
filter(sum(.data$Proportion, na.rm = TRUE) > 0) |>
ungroup()
age_bin_columns <- setdiff(names(age_freq), c("Fishery", "Year", "N", "MinAge", "MaxAge"))
age_snapshot <- age_freq |>
mutate(
Fishery = recode(as.character(.data$Fishery), !!!fishery_lookup),
Year = as.integer(.data$Year)
) |>
group_by(.data$Fishery) |>
filter(.data$Year %in% recent_years(.data$Year)) |>
ungroup() |>
pivot_longer(all_of(age_bin_columns), names_to = "Age", values_to = "Proportion") |>
mutate(
Age = as.numeric(.data$Age),
Year = factor(.data$Year)
) |>
filter(.data$Age >= .data$MinAge, .data$Age <= .data$MaxAge) |>
group_by(.data$Fishery, .data$Year, .data$Age) |>
summarise(Proportion = sum(.data$Proportion, na.rm = TRUE), .groups = "drop") |>
group_by(.data$Fishery, .data$Year) |>
filter(sum(.data$Proportion, na.rm = TRUE) > 0) |>
ungroup()
paly_years <- as.integer(dimnames(data$paly)[[3]])
paly_nonzero_years <- paly_years[apply(data$paly, 3, sum) > 0]
paly_plot_years <- sort(unique(c(
min(paly_nonzero_years),
max(paly_nonzero_years)
)))
paly_plot_data <- as_tibble(
as.data.frame.table(
data$paly[, , as.character(paly_plot_years), drop = FALSE],
responseName = "Probability"
)
) |>
transmute(
LengthBin = as.integer(.data$Var1),
Age = as.integer(as.character(.data$Var2)),
Year = factor(as.integer(as.character(.data$Var3)), levels = paly_plot_years),
Probability = as.numeric(.data$Probability)
)
paly_length_lookup <- paly_plot_data |>
mutate(
YearInteger = as.integer(as.character(.data$Year)),
YearIndex = .data$YearInteger - data$first_yr + 1L,
AgeIndex = .data$Age - data$min_age + 1L,
LengthAtAge = data$length_mu_ysa[cbind(.data$YearIndex, 1L, .data$AgeIndex)]
) |>
group_by(.data$LengthBin) |>
summarise(
Length = ifelse(sum(.data$Probability, na.rm = TRUE) > 0,
weighted.mean(.data$LengthAtAge, w = .data$Probability, na.rm = TRUE),
NA_real_
),
.groups = "drop"
)
paly_plot_data <- paly_plot_data |>
left_join(paly_length_lookup, by = join_by(LengthBin)) |>
filter(is.finite(.data$Length), .data$Age >= 5) |>
mutate(LengthIndex = .data$LengthBin)
paly_age_breaks <- sort(unique(paly_plot_data$Age))
paly_length_axis <- paly_plot_data |>
distinct(.data$LengthIndex, .data$Length) |>
arrange(.data$LengthIndex)
paly_length_axis_index <- unique(round(seq(1, nrow(paly_length_axis), length.out = min(6, nrow(paly_length_axis)))))
paly_y_breaks <- paly_length_axis$LengthIndex[paly_length_axis_index]
paly_probability_limits <- c(0, max(paly_plot_data$Probability, na.rm = TRUE))
```
```{r}
#| label: fig-input-length-at-age
#| fig-cap: "Mean length at age by season for the first, middle, and final input years."
#| fig-width: 9
#| fig-height: 5
plot_length_at_age(data = data, years = growth_plot_years)
```
```{r}
#| label: fig-input-weight-at-age
#| fig-cap: "Weight at age by fishery for the first, middle, and final input years."
#| fig-width: 9
#| fig-height: 5
plot_weight_at_age(data = data, years = growth_plot_years)
```
```{r}
#| label: fig-input-paly
#| fig-cap: "Conditional probability of age given length and year (`paly`) for the first and final non-zero PALY years. Length is shown in centimeters using PALY probabilities and model length-at-age inputs."
#| fig-width: 10
#| fig-height: 5.8
ggplot(paly_plot_data, aes(x = .data$Age, y = .data$LengthIndex, fill = .data$Probability)) +
geom_tile(width = 1, height = 1, color = "grey82", linewidth = 0.18, alpha = 0.86) +
facet_wrap(~Year, nrow = 1) +
labs(x = "Age", y = "Length (cm)", fill = "Probability") +
scale_x_continuous(limits = bin_limits(paly_plot_data$Age, step = 1), breaks = paly_age_breaks, expand = expansion(mult = 0)) +
scale_y_continuous(
limits = bin_limits(paly_plot_data$LengthIndex, step = 1),
breaks = paly_y_breaks,
labels = function(x) round(paly_length_axis$Length[match(x, paly_length_axis$LengthIndex)]),
expand = expansion(mult = 0)
) +
scale_fill_viridis_c(option = "magma", direction = -1, limits = paly_probability_limits, oob = squish) +
coord_cartesian(clip = "off") +
theme(panel.grid = element_blank())
```
```{r}
#| label: fig-input-catch
#| fig-cap: "Observed catch inputs by fishery."
#| fig-width: 9
#| fig-height: 7
ggplot(catch_plot_data, aes(x = .data$Year, y = .data$Catch)) +
geom_line(color = fit_expected_color, linewidth = 0.45) +
geom_point(color = fit_expected_color, size = 0.9) +
facet_wrap(~Fishery, scales = "free_y", ncol = 2) +
labs(x = "Year", y = "Catch") +
scale_x_continuous(limits = year_limits(catch_plot_data$Year), breaks = pretty_breaks()) +
scale_y_continuous(limits = c(0, NA), expand = expansion(mult = c(0, 0.05)))
```
```{r}
#| label: fig-input-indices
#| fig-cap: "Input index and monitoring series."
#| fig-width: 9
#| fig-height: 6
ggplot(index_plot_data, aes(x = .data$Year, y = .data$Value)) +
geom_line(color = fit_expected_color, linewidth = 0.45) +
geom_point(color = fit_expected_color, size = 1.2) +
facet_wrap(~Series, scales = "free_y", ncol = 2) +
labs(x = "Year", y = "Input value") +
scale_x_continuous(limits = year_limits(index_plot_data$Year), breaks = pretty_breaks()) +
scale_y_continuous(limits = c(0, NA), expand = expansion(mult = c(0, 0.05)))
```
```{r}
#| label: fig-input-length-frequency
#| fig-cap: "Observed length-frequency proportions for the three most recent years with length-frequency inputs."
#| fig-width: 11
#| fig-height: 9
ggplot(length_snapshot, aes(x = .data$Length, y = .data$Proportion, color = .data$Year)) +
geom_line(linewidth = 0.5) +
geom_point(size = 0.8) +
facet_wrap(~Fishery, scales = "free_y", ncol = 2) +
labs(x = "Length bin", y = "Proportion", color = "Year") +
scale_x_continuous(breaks = pretty_breaks(n = 5)) +
scale_y_continuous(limits = c(0, NA), expand = expansion(mult = c(0, 0.05)))
```
```{r}
#| label: fig-input-age-frequency
#| fig-cap: "Observed age-frequency proportions for the three most recent years with age-frequency inputs in each fishery."
#| fig-width: 9
#| fig-height: 5.5
ggplot(age_snapshot, aes(x = .data$Age, y = .data$Proportion, color = .data$Year)) +
geom_line(linewidth = 0.6) +
geom_point(size = 1.4) +
facet_wrap(~Fishery, scales = "free_y", ncol = 1) +
labs(x = "Age", y = "Proportion", color = "Year") +
scale_x_continuous(breaks = pretty_breaks(n = 6)) +
scale_y_continuous(limits = c(0, NA), expand = expansion(mult = c(0, 0.05)))
```
# Model setup
```{r}
#| label: setup-parameters
parameters <- get_parameters(data = data)
parameters$par_log_h <- log(0.72)
parameters$par_log_psi <- log(1.75)
b0_start_multiplier <- 1.05
parameters$par_log_B0 <- parameters$par_log_B0 + log(b0_start_multiplier)
map <- get_map(parameters = parameters)
data$priors <- get_priors(parameters = parameters)
obj <- MakeADFun(func = cmb(sbt_model, data), parameters = parameters, map = map)
bounds <- get_bounds(obj = obj, parameters = parameters)
```
The initial B₀ value is increased by 5% relative to the default
`get_parameters()` value before the objective is created. This small upward
nudge gives the optimizer a less boundary-like starting biomass without changing
the model structure.
```{r}
#| label: tbl-starting-values
starting_values <- tibble(
Parameter = parameter_labels(c("Steepness", "Psi", "B0", "M0", "M4", "M10", "M30")),
`Starting value` = c(
exp(parameters$par_log_h),
exp(parameters$par_log_psi),
exp(parameters$par_log_B0),
exp(parameters$par_log_m0),
exp(parameters$par_log_m4),
exp(parameters$par_log_m10),
exp(parameters$par_log_m30)
)
)
kable(
starting_values |>
mutate(`Starting value` = format_decimal(.data$`Starting value`, digits = 3)) |>
replace_missing(),
align = numeric_table_align(starting_values),
caption = "Key starting values used to initialize the OMMP16 fit."
)
```
```{r}
#| label: tbl-estimated-parameter-names
estimated_parameter_names <- tibble(`Estimated parameter` = names(obj$par)) |>
count(.data$`Estimated parameter`, name = "Number of active parameters") |>
arrange(.data$`Estimated parameter`)
estimated_parameter_names <- bind_rows(
estimated_parameter_names,
tibble(
`Estimated parameter` = "Total",
`Number of active parameters` = length(obj$par)
)
)
kable(
replace_missing(estimated_parameter_names),
align = c("l", "r"),
caption = "Active parameter blocks in the maximum-likelihood optimization, including the number of estimated parameters in each block."
)
```
# Maximum-likelihood fit
```{r}
#| label: run-mle
#| results: hide
control <- list(eval.max = 10000, iter.max = 10000)
initial_nll <- obj$fn(obj$par)
opt <- nlminb(
start = obj$par,
objective = obj$fn,
gradient = obj$gr,
hessian = obj$he,
lower = bounds$lower,
upper = bounds$upper,
control = control
)
opt <- nlminb(
start = opt$par,
objective = obj$fn,
gradient = obj$gr,
hessian = obj$he,
lower = bounds$lower,
upper = bounds$upper,
control = control
)
obj$par <- opt$par
obj$env$last.par.best <- opt$par
final_nll <- obj$fn(opt$par)
max_gradient <- max(abs(obj$gr(opt$par)))
```
```{r}
#| label: estimability-check
#| include: false
estimability_result <- tryCatch(
{
hessian <- obj$he(opt$par)
check_estimability(obj = obj, h = hessian)
},
error = function(e) e
)
estimability_summary <- if (inherits(estimability_result, "error")) {
tibble(`Estimability result` = "Failed", Detail = conditionMessage(estimability_result))
} else if (length(estimability_result$WhichBad) == 0L) {
tibble(`Estimability result` = "All parameters are estimable", Detail = paste("Minimum Hessian eigenvalue:", signif(min(estimability_result$Eigen$values), 4)))
} else {
bad_params <- estimability_result$BadParams |>
filter(.data$Param_check != "OK") |>
distinct(.data$Param) |>
pull(.data$Param)
tibble(`Estimability result` = "Some parameters are not estimable", Detail = paste(bad_params, collapse = ", "))
}
```
```{r}
#| label: tbl-optimisation-summary
optimisation_summary <- tibble(
Section = "Optimization",
Statistic = c("Initial negative log-likelihood", "Final negative log-likelihood", "Convergence code", "Convergence message", "Maximum absolute gradient"),
Value = c(initial_nll, final_nll, opt$convergence, opt$message, max_gradient)
)
estimability_table <- estimability_summary |>
pivot_longer(everything(), names_to = "Statistic", values_to = "Value") |>
mutate(Section = "Estimability", .before = 1)
optimisation_estimability_summary <- bind_rows(
optimisation_summary,
estimability_table
)
kable(
format_table(optimisation_estimability_summary) |>
replace_missing(),
caption = "Maximum-likelihood optimization and estimability summary. Estimability was checked using the AD Hessian at the maximum-likelihood estimate."
)
```
In @tbl-optimisation-summary, the initial negative log-likelihood is the
objective value at the starting parameter vector and the final negative
log-likelihood is the value after the two `nlminb()` passes. A convergence code
of 0 from [`nlminb()`](https://stat.ethz.ch/R-manual/R-devel/library/stats/html/nlminb.html)
means the optimizer reported successful convergence; the convergence message is
the corresponding optimizer stopping message. The maximum absolute gradient is a
local first-order check on the optimized parameter vector, with smaller values
indicating a flatter objective at the solution. The estimability rows report the
`check_estimability()` result and the associated detail from the AD Hessian at
the maximum-likelihood estimate.
::: {.callout-warning}
## Historical fit is not an accepted ESC31 result
This reconstruction stopped with convergence code
`r opt$convergence` and maximum absolute gradient
`r format(max_gradient, scientific = TRUE, digits = 3)`. It therefore fails the
current ESC31 numerical acceptance gates and is retained only as an OMMP16
diagnostic checkpoint. The accepted base fit and posterior are reported in
[the ESC31 base-model page](../ESC31/2_base.html).
:::
# Parameter tables
```{r}
#| label: tbl-estimated-parameters
estimated_parameters <- make_parameter_table(
obj = obj,
data = data,
exclude_recruitment_devs = TRUE,
include_selectivity_devs = FALSE,
include_details = TRUE
) |>
filter(.data$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(.data$Parameter))
kable(
format_decimal_table(estimated_parameters, digits = 3) |>
replace_missing(),
align = numeric_table_align(estimated_parameters),
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))
```
```{r}
#| label: tbl-selectivity-hyperparameters
par_list <- obj$env$parList(obj$env$last.par.best)
fishery_names <- c("LL1", "LL2", "LL3", "LL4", "Indonesia", "Australia", "CPUE")
selectivity_hyperparameters <- tibble(
Fishery = fishery_names,
`ρ year` = sel_rho_from_par(par_list$par_sel_rho_y),
`ρ age` = sel_rho_from_par(par_list$par_sel_rho_a),
`σ` = exp(par_list$par_log_sel_sigma)
)
kable(
format_table(selectivity_hyperparameters) |>
replace_missing(),
caption = "Selectivity hyper-parameters by fishery on natural scale. Status columns are omitted because these hyper-parameters are fixed by the active parameter map for this fit."
)
```
# Model fit plots
```{r}
#| label: make-fit-plots
#| include: false
hsp_plot_list <- tryCatch(
plot_hsps(data = data, object = obj, return_list = TRUE),
error = function(e) e
)
hsp_panel <- function(name, label) {
make_plot(label, function() {
if (inherits(hsp_plot_list, "error")) {
stop(conditionMessage(hsp_plot_list), call. = FALSE)
}
if (is.null(hsp_plot_list[[name]])) {
stop("HSP panel not found: ", name, call. = FALSE)
}
p <- hsp_plot_list[[name]]
p <- harmonize_ck_colors(p)
if (name %in% c("pop_by_adult_capture_age", "pop_by_adult_capture_year")) {
p <- remove_expected_points(p)
}
if (name == "pop_by_juvenile_cohort") {
x <- data$pop_obs[, "Cohort"]
p <- p + scale_x_continuous(
limits = pad_limits(x),
breaks = pretty_breaks(),
labels = function(z) z + data$first_yr - 1
)
}
if (name == "pop_by_adult_capture_year") {
x <- data$pop_obs[, "CaptureYear"]
p <- p + scale_x_continuous(
limits = pad_limits(x),
breaks = pretty_breaks(),
labels = function(z) z + data$first_yr - 1
)
}
if (name == "pop_by_adult_capture_age") {
p <- p + scale_x_continuous(limits = age_limits(data), breaks = pretty_breaks())
}
if (name == "hsp_by_cohort_pair") {
p <- p + scale_x_continuous(limits = pad_limits(data$hsp_obs[, "cmax"]), breaks = pretty_breaks())
}
if (name == "hsp_by_initial_cohort") {
p <- p + scale_x_continuous(limits = pad_limits(data$hsp_obs[, "cmin"]), breaks = pretty_breaks())
}
if (inherits(p, "ggplot")) {
p <- p + labs(y = "Number of matches")
}
if (inherits(p, "ggplot")) p <- zero_y(p)
p
})
}
fit_plots <- list(
natural_mortality = make_plot("Natural mortality", function() plot_natural_mortality_report(data = data, object = obj)),
recruitment = make_plot("Recruitment", function() plot_recruitment_report(data = data, object = obj)),
recruitment_deviates = make_plot("Recruitment deviates", function() plot_rec_devs_report(data = data, object = obj)),
initial_numbers = make_plot("Initial numbers", function() plot_initial_numbers_report(data = data, object = obj)),
catch = make_plot("Catch input and output", function() plot_catch_input_output(data = data, object = obj)),
selectivity_ll1 = make_plot("Selectivity: LL1", function() plot_selectivity(data = data, object = obj, fisheries = "LL1", years = data$sel_LL1_yrs)),
selectivity_ll2 = make_plot("Selectivity: LL2", function() plot_selectivity(data = data, object = obj, fisheries = "LL2", years = data$sel_LL2_yrs)),
selectivity_ll3 = make_plot("Selectivity: LL3", function() plot_selectivity(data = data, object = obj, fisheries = "LL3", years = data$sel_LL3_yrs)),
selectivity_indonesian = make_plot("Selectivity: Indonesian", function() plot_selectivity(data = data, object = obj, fisheries = "Indonesian", years = data$sel_Ind_yrs)),
selectivity_australian = make_plot("Selectivity: Australian", function() plot_selectivity(data = data, object = obj, fisheries = "Australian", years = data$sel_Aus_yrs)),
selectivity_cpue = make_plot("Selectivity: CPUE", function() plot_selectivity(data = data, object = obj, fisheries = "CPUE", years = data$sel_CPUE_yrs)),
cpue_fit = make_plot("CPUE fit", function() plot_cpue(data = data, object = obj)),
aerial_fit = make_plot("Aerial survey fit", function() plot_aerial_fit_ci(data = data, object = obj)),
troll_fit = make_plot("Troll index fit", function() plot_troll(data = data, object = obj)),
gt_fit = make_plot("Gene-tagging fit", function() plot_gt(data = data, object = obj)),
hsp_pop_cohort = hsp_panel("pop_by_juvenile_cohort", "POPs by juvenile cohort"),
hsp_pop_age = hsp_panel("pop_by_adult_capture_age", "POPs by adult capture age"),
hsp_pop_year = hsp_panel("pop_by_adult_capture_year", "POPs by adult capture year"),
hsp_pair = hsp_panel("hsp_by_cohort_pair", "HSPs by cohort pair"),
hsp_initial = hsp_panel("hsp_by_initial_cohort", "HSPs by initial cohort"),
hsp_total = make_plot("Total POP, HSP, and GT matches", function() plot_total_genetics_report(data = data, object = obj, hsp_plot_list = hsp_plot_list)),
lf_ll1 = make_plot("LL1 length composition fit", function() harmonize_composition_fit_colors(plot_lf(data = data, object = obj, fishery = "LL1"))),
lf_ll2 = make_plot("LL2 length composition fit", function() harmonize_composition_fit_colors(plot_lf(data = data, object = obj, fishery = "LL2"))),
lf_ll3 = make_plot("LL3 length composition fit", function() harmonize_composition_fit_colors(plot_lf(data = data, object = obj, fishery = "LL3"))),
lf_ll4 = make_plot("LL4 length composition fit", function() harmonize_composition_fit_colors(plot_lf(data = data, object = obj, fishery = "LL4"))),
lf_cpue = make_plot("CPUE length composition fit", function() harmonize_composition_fit_colors(plot_lf(data = data, object = obj, fishery = "CPUE"))),
af_indonesia = make_plot("Indonesian age composition fit", function() harmonize_composition_fit_colors(plot_af(data = data, object = obj, fishery = "Indonesian"))),
af_australia = make_plot("Australian age composition fit", function() harmonize_composition_fit_colors(plot_af(data = data, object = obj, fishery = "Australian")))
)
```
## Fishery and index fits
```{r}
#| label: fig-catch-fit
#| fig-cap: "Input and output catch by fishery and season. Orange points are input catch; blue lines are model output catch. Catch is an input, not a fitted observation."
#| fig-width: 9
#| fig-height: 6
fit_plots$catch$plot
```
```{r}
#| label: fig-cpue-fit
#| fig-cap: "Observed and predicted CPUE index values with 95% intervals from the input CPUE standard deviations plus the fitted model observation error. Orange points are observed values and blue lines are expected values."
#| fig-width: 9
#| fig-height: 5
fit_plots$cpue_fit$plot
```
```{r}
#| label: fig-aerial-fit
#| fig-cap: "Observed and predicted aerial survey index values with 95% intervals from the aerial covariance matrix plus the fitted model observation error. Orange points are observed values and blue lines are expected values."
#| fig-width: 9
#| fig-height: 5
fit_plots$aerial_fit$plot
```
```{r}
#| label: fig-troll-fit
#| fig-cap: "Observed and predicted troll index values with 95% intervals from the input troll standard deviations plus the fixed troll observation error. Orange points are observed values and blue lines are expected values. The troll index is switched off in this fit, so this plot is diagnostic only and does not contribute to the objective function."
#| fig-width: 9
#| fig-height: 5
fit_plots$troll_fit$plot
```
```{r}
#| label: fig-gt-fit
#| fig-cap: "Observed and expected gene-tagging (GT) matches by release year with 95% binomial intervals. Orange points are observed matches, the solid blue line is expected matches, and dashed blue lines are the lower and upper 95% binomial bounds."
#| fig-width: 9
#| fig-height: 5
fit_plots$gt_fit$plot
```
```{r}
#| label: fig-pop-juvenile-cohort
#| fig-cap: "Observed and predicted parent-offspring-pair matches (POP) by juvenile cohort. Orange points are observed matches and blue lines and intervals are expected matches."
#| fig-width: 9
#| fig-height: 5
fit_plots$hsp_pop_cohort$plot
```
```{r}
#| label: fig-pop-adult-age
#| fig-cap: "Observed and predicted parent-offspring-pair matches (POP) by adult capture age. Orange points are observed matches and blue lines and intervals are expected matches."
#| fig-width: 9
#| fig-height: 5
fit_plots$hsp_pop_age$plot
```
```{r}
#| label: fig-pop-adult-year
#| fig-cap: "Observed and predicted parent-offspring-pair matches (POP) by adult capture year. Orange points are observed matches and blue lines and intervals are expected matches."
#| fig-width: 9
#| fig-height: 5
fit_plots$hsp_pop_year$plot
```
```{r}
#| label: fig-hsp-cohort-pair
#| fig-cap: "Observed and predicted half-sibling-pair matches (HSP) by cohort pair. Orange points are observed matches and blue lines and intervals are expected matches."
#| fig-width: 11
#| fig-height: 7
fit_plots$hsp_pair$plot
```
```{r}
#| label: fig-hsp-initial-cohort
#| fig-cap: "Observed and predicted half-sibling-pair matches (HSP) by initial cohort. Orange points are observed matches and blue lines and intervals are expected matches."
#| fig-width: 9
#| fig-height: 5
fit_plots$hsp_initial$plot
```
```{r}
#| label: fig-close-kin-total
#| fig-cap: "Observed and predicted total parent-offspring-pair (POP), half-sibling-pair (HSP), and gene-tagging (GT) matches. Orange points are observed matches and blue points and intervals are expected matches."
#| fig-width: 7
#| fig-height: 5
fit_plots$hsp_total$plot
```
## Composition fits
```{r}
#| label: fig-lf-ll1
#| fig-cap: "Length-composition fit for LL1. Orange points are observed proportions and blue dashed lines are expected proportions."
#| fig-width: 10
#| fig-height: 10
fit_plots$lf_ll1$plot
```
```{r}
#| label: fig-lf-ll2
#| fig-cap: "Length-composition fit for LL2. Orange points are observed proportions and blue dashed lines are expected proportions."
#| fig-width: 10
#| fig-height: 6
fit_plots$lf_ll2$plot
```
```{r}
#| label: fig-lf-ll3
#| fig-cap: "Length-composition fit for LL3. Orange points are observed proportions and blue dashed lines are expected proportions."
#| fig-width: 10
#| fig-height: 9.5
fit_plots$lf_ll3$plot
```
```{r}
#| label: fig-lf-ll4
#| fig-cap: "Length-composition fit for LL4. Orange points are observed proportions and blue dashed lines are expected proportions."
#| fig-width: 10
#| fig-height: 6
fit_plots$lf_ll4$plot
```
```{r}
#| label: fig-lf-cpue
#| fig-cap: "Length-composition fit for the CPUE length-frequency data. Orange points are observed proportions and blue dashed lines are expected proportions."
#| fig-width: 10
#| fig-height: 10
fit_plots$lf_cpue$plot
```
```{r}
#| label: fig-af-indonesia
#| fig-cap: "Age-composition fit for the Indonesian fishery. Orange points are observed proportions and blue dashed lines are expected proportions."
#| fig-width: 10
#| fig-height: 6
fit_plots$af_indonesia$plot
```
```{r}
#| label: fig-af-australia
#| fig-cap: "Age-composition fit for the Australian fishery. Orange points are observed proportions and blue dashed lines are expected proportions."
#| fig-width: 10
#| fig-height: 8
fit_plots$af_australia$plot
```
## Selectivity fits
Selectivity panels are restricted to the configured `sel_*_yrs` values for each
fishery. The default plot spans every model year, which creates large gaps
because selectivity is only configured for selected years.
```{r}
#| label: fig-selectivity-ll1
#| fig-cap: "Estimated selectivity-at-age for LL1. Dashed vertical lines mark the minimum and maximum ages over which selectivity is estimated; crosses mark years with composition observations shown in @fig-lf-ll1."
#| fig-width: 9
#| fig-height: !expr selectivity_fig_height(data$sel_LL1_yrs)
fit_plots$selectivity_ll1$plot
```
```{r}
#| label: fig-selectivity-ll2
#| fig-cap: "Estimated selectivity-at-age for LL2. Dashed vertical lines mark the minimum and maximum ages over which selectivity is estimated; crosses mark years with composition observations shown in @fig-lf-ll2."
#| fig-width: 9
#| fig-height: !expr selectivity_fig_height(data$sel_LL2_yrs)
fit_plots$selectivity_ll2$plot
```
```{r}
#| label: fig-selectivity-ll3
#| fig-cap: "Estimated selectivity-at-age for LL3. Dashed vertical lines mark the minimum and maximum ages over which selectivity is estimated; crosses mark years with composition observations shown in @fig-lf-ll3."
#| fig-width: 9
#| fig-height: !expr selectivity_fig_height(data$sel_LL3_yrs)
fit_plots$selectivity_ll3$plot
```
```{r}
#| label: fig-selectivity-indonesian
#| fig-cap: "Estimated selectivity-at-age for the Indonesian fishery. Dashed vertical lines mark the minimum and maximum ages over which selectivity is estimated; crosses mark years with composition observations shown in @fig-af-indonesia."
#| fig-width: 9
#| fig-height: !expr selectivity_fig_height(data$sel_Ind_yrs)
fit_plots$selectivity_indonesian$plot
```
```{r}
#| label: fig-selectivity-australian
#| fig-cap: "Estimated selectivity-at-age for the Australian fishery. Dashed vertical lines mark the minimum and maximum ages over which selectivity is estimated; crosses mark years with composition observations shown in @fig-af-australia."
#| fig-width: 9
#| fig-height: !expr selectivity_fig_height(data$sel_Aus_yrs)
fit_plots$selectivity_australian$plot
```
```{r}
#| label: fig-selectivity-cpue
#| fig-cap: "Estimated selectivity-at-age for the CPUE fleet. Dashed vertical lines mark the minimum and maximum ages over which selectivity is estimated; crosses mark years with composition observations shown in @fig-lf-cpue. The high value at the first displayed age is an edge effect from the lower bound of the CPUE selectivity window rather than evidence for an additional younger-age mode."
#| fig-width: 9
#| fig-height: !expr selectivity_fig_height(data$sel_CPUE_yrs)
fit_plots$selectivity_cpue$plot
```
# Residual diagnostics
One-step-ahead (OSA) residuals provide randomized quantile-style diagnostics
for non-Gaussian observations and are useful for assessing stock-assessment fit
across observation types [@DunnSmyth1996; @StewartMonnahan2025]. 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 [@Francis2011; @McAllisterIanelli1997].
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 @tbl-residual-statistics.
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.
```{r}
#| label: make-residual-plots
#| include: false
residual_plots <- list(
cpue = make_plot("CPUE index", function() plot_cpue_residuals(data, obj)),
aerial = make_plot("Aerial survey", function() plot_aerial_residuals(data, obj)),
troll = make_plot("Troll survey", function() plot_troll_residuals(data, obj)),
gt = make_plot("Gene tagging", function() plot_gt_residuals(data, obj)),
hsp = make_plot("Half-sibling pairs", function() plot_hsps_residuals(data, obj)),
pop = make_plot("Close-kin POPs", function() plot_pops_residuals_report(data, obj)),
tags = make_plot("Conventional tags", function() {
p <- plot_tags_residuals(data, obj)
attr(p, "residual_source") <- "compResidual"
p
}),
lf_ll1 = make_plot("Length composition: LL1", function() plot_lf_residuals(data, obj, fishery = "LL1")),
lf_ll2 = make_plot("Length composition: LL2", function() plot_lf_residuals(data, obj, fishery = "LL2")),
lf_ll3 = make_plot("Length composition: LL3", function() plot_lf_residuals(data, obj, fishery = "LL3")),
lf_cpue = make_plot("Length composition: CPUE", function() plot_lf_residuals(data, obj, fishery = "CPUE")),
af_indonesia = make_plot("Age composition: Indonesian", function() plot_af_residuals(data, obj, fishery = "Indonesian")),
af_australia = make_plot("Age composition: Australian", function() plot_af_residuals(data, obj, fishery = "Australian"))
)
residual_statistics <- map_dfr(
residual_plots,
~ residual_stat_row(.x$label, .x$plot, .x$error)
)
```
```{r}
#| label: tbl-residual-statistics
kable(
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(.data$Source)) |>
format_residual_table(digits = 3) |>
replace_missing(),
align = c("l", rep("r", 5), "l", rep("r", 6)),
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, " " = 2, "Francis" = 2, "McAllister-Ianelli" = 4))
```
```{r}
#| label: fig-cpue-residuals
#| fig-cap: "One-step-ahead (OSA) residuals for the CPUE index, derived with oneStepPredict on the lognormal CPUE observation."
#| fig-width: 9
#| fig-height: 5
residual_plots$cpue$plot
```
```{r}
#| label: fig-aerial-residuals
#| fig-cap: "One-step-ahead (OSA) residuals for the aerial survey index, derived with oneStepPredict on the lognormal aerial-survey observation."
#| fig-width: 9
#| fig-height: 5
residual_plots$aerial$plot
```
```{r}
#| label: fig-troll-residuals
#| fig-cap: "Residual diagnostics for the troll index, derived as likelihood-scale standardized lognormal residuals. The troll index is switched off in this fit, so this panel records that status rather than a fitted residual pattern."
#| fig-width: 9
#| fig-height: 5
residual_plots$troll$plot
```
```{r}
#| label: fig-gt-residuals
#| fig-cap: "One-step-ahead (OSA) residuals for gene-tagging recaptures, derived with oneStepPredict on the binomial GT recapture observation."
#| fig-width: 9
#| fig-height: 5
residual_plots$gt$plot
```
```{r}
#| label: fig-hsp-residuals
#| fig-cap: "One-step-ahead (OSA) residuals for half-sibling-pair observations, derived with oneStepPredict on the binomial HSP observation."
#| fig-width: 9
#| fig-height: 7
residual_plots$hsp$plot
```
```{r}
#| label: fig-pop-residuals
#| fig-cap: "One-step-ahead (OSA) residual map for close-kin POP observations, derived with oneStepPredict on the binomial POP observation and averaged within release-cohort and adult-capture-year cells."
#| fig-width: 11
#| fig-height: 7
residual_plots$pop$plot
```
```{r}
#| label: fig-tag-residuals
#| fig-cap: "Residual diagnostics for conventional tag recaptures, derived with compResidual for the Dirichlet-multinomial recapture categories."
#| fig-width: 12
#| fig-height: 8
residual_plots$tags$plot
```
```{r}
#| label: fig-lf-ll1-residuals
#| fig-cap: "Length-composition residual diagnostics for LL1, derived with oneStepPredict where finite and the compResidual fallback otherwise."
#| fig-width: 11
#| fig-height: 6
residual_plots$lf_ll1$plot
```
```{r}
#| label: fig-lf-ll2-residuals
#| fig-cap: "Length-composition residual diagnostics for LL2, derived with oneStepPredict where finite and the compResidual fallback otherwise."
#| fig-width: 11
#| fig-height: 6
residual_plots$lf_ll2$plot
```
```{r}
#| label: fig-lf-ll3-residuals
#| fig-cap: "Length-composition residual diagnostics for LL3, derived with oneStepPredict where finite and the compResidual fallback otherwise."
#| fig-width: 11
#| fig-height: 6
residual_plots$lf_ll3$plot
```
```{r}
#| label: fig-lf-cpue-residuals
#| fig-cap: "Length-composition residual diagnostics for CPUE length-frequency data, derived with oneStepPredict where finite and the compResidual fallback otherwise."
#| fig-width: 11
#| fig-height: 6
residual_plots$lf_cpue$plot
```
```{r}
#| label: fig-af-indonesia-residuals
#| fig-cap: "Age-composition residual diagnostics for the Indonesian fishery, derived with oneStepPredict where finite and the compResidual fallback otherwise."
#| fig-width: 11
#| fig-height: 6
residual_plots$af_indonesia$plot
```
```{r}
#| label: fig-af-australia-residuals
#| fig-cap: "Age-composition residual diagnostics for the Australian fishery, derived with oneStepPredict where finite and the compResidual fallback otherwise."
#| fig-width: 11
#| fig-height: 6
residual_plots$af_australia$plot
```
## Population dynamics
```{r}
#| label: fig-initial-numbers
#| fig-cap: "Initial numbers at age from the maximum-likelihood fit."
#| fig-width: 8
#| fig-height: 5
fit_plots$initial_numbers$plot
```
```{r}
#| label: fig-natural-mortality
#| fig-cap: "Estimated natural mortality at age from the maximum-likelihood fit."
#| fig-width: 8
#| fig-height: 5
fit_plots$natural_mortality$plot
```
```{r}
#| label: fig-recruitment-deviates
#| fig-cap: "Recruitment deviates from the maximum-likelihood fit. The black dashed horizontal line marks zero; orange points identify the final three recruitment deviations penalized with the AR1 prior, and blue points identify deviations penalized with independent normal priors."
#| fig-width: 9
#| fig-height: 5
fit_plots$recruitment_deviates$plot
```
```{r}
#| label: fig-recruitment
#| fig-cap: "Recruitment trajectory from the maximum-likelihood fit. The black dashed horizontal line is unfished recruitment; orange points identify the final three recruitment values associated with the AR1 recruitment-deviation prior, and blue points identify earlier recruitment values."
#| fig-width: 9
#| fig-height: 5
fit_plots$recruitment$plot
```
# CPUE sensitivity fits
The historical sensitivity fits were configured in `script_OMMP16.R` and saved
to `sensitivity/ommp16_sensitivities.rds` when run. That result file is not
retained in the current checkout. The figures below therefore display labelled
unavailable panels so the historical report remains reproducible; they are not
unfinished ESC31 outputs and should not be regenerated for current sign-off.
The accepted current sensitivity set is reported in
[the ESC31 sensitivity page](../ESC31/3_sens.html).
```{r}
#| label: read-sensitivity-results
#| include: false
sensitivity_results_path <- file.path("sensitivity", "ommp16_sensitivities.rds")
sensitivity_results <- if (file.exists(sensitivity_results_path)) {
readRDS(sensitivity_results_path)
} else {
NULL
}
sens_q2008 <- restore_sensitivity_fit(sensitivity_results$q2008 %||% NULL)
sens_cpue_tv_cv <- restore_sensitivity_fit(sensitivity_results$cpue_tv_cv %||% NULL)
cpue_tv_cv_input <- sensitivity_results$cpue_tv_cv_input %||%
read_cpue_time_varying_cv(cpue_base = cpue, cpue_cv = cpue_cv)
cpue_tv_cv_mean <- mean(cpue_tv_cv_input$CV, na.rm = TRUE)
sensitivity_summary <- bind_rows(
summarize_fit("Base OMMP16", data, obj, opt),
if (has_sensitivity_fit(sens_q2008)) {
summarize_fit("Sensitivity 1: CPUE q split from 2008", sens_q2008$data, sens_q2008$obj, sens_q2008$opt)
},
if (has_sensitivity_fit(sens_cpue_tv_cv)) {
summarize_fit("Sensitivity 2: time-varying CPUE CV", sens_cpue_tv_cv$data, sens_cpue_tv_cv$obj, sens_cpue_tv_cv$opt)
}
)
```
```{r}
#| label: tbl-sensitivity-status
sensitivity_status <- tibble(
`Output file` = sensitivity_results_path,
Status = ifelse(
file.exists(sensitivity_results_path),
"Historical output file found",
"Historical output not retained; see ESC31/3_sens.html"
)
)
kable(
replace_missing(sensitivity_status),
caption = "Sensitivity-output status for the report."
)
```
```{r}
#| label: tbl-sensitivity-summary
kable(
sensitivity_summary |>
mutate(across(c(B0, M0, M4, M10, M30, h, psi), ~ format_decimal(.x, digits = 3))) |>
format_decimal_table(digits = 3) |>
replace_missing(),
align = c("l", rep("r", ncol(sensitivity_summary) - 1)),
caption = "Maximum-likelihood comparison for the base model and available CPUE sensitivity fits. Convergence code 0 indicates normal nlminb convergence."
)
```
## Sensitivity 1: CPUE q split from 2008
Sensitivity 1 adds a CPUE catchability split in 2008. This checks whether a
post-2008 CPUE scale shift changes the core population trajectory or key
estimated parameters.
```{r}
#| label: fig-sensitivity-q2008-cpue
#| fig-cap: "CPUE index fit for Sensitivity 1, with a catchability split from 2008. Orange points are observed values and blue lines are expected values."
#| fig-width: 9
#| fig-height: 5
sensitivity_q2008_cpue_plot <- make_plot("Sensitivity 1 CPUE fit", function() {
if (!has_sensitivity_fit(sens_q2008)) {
stop("Run the Sensitivity 1 section in script_OMMP16.R to create ", sensitivity_results_path, call. = FALSE)
}
plot_cpue(data = sens_q2008$data, object = sens_q2008$obj)
})
sensitivity_q2008_cpue_plot$plot
```
## Sensitivity 2: time-varying CPUE CV
Sensitivity 2 uses the year-specific CPUE CVs in `CV_out_10.csv`, taking the
`cv_scaled` column as the model CV while leaving the CPUE index values from
`cpue.csv` unchanged. The mean scaled CV across years is
`r format_decimal(cpue_tv_cv_mean, digits = 3)`, matching the intended 0.20
target. The scaling preserves the relative annual pattern in `cv_raw`, but
forces the average uncertainty level to 20% using
`cv_scaled = cv_raw * 0.20 / mean(cv_raw)`.
```{r}
#| label: fig-sensitivity-cpue-cv-input
#| fig-cap: "Time-varying CPUE CV used in Sensitivity 2. Values are read from the cv_scaled column of CV_out_10.csv and aligned to cpue.csv by year."
#| fig-width: 9
#| fig-height: 4.5
cpue_tv_cv_input |>
ggplot(aes(x = .data$Year, y = .data$CV)) +
geom_line(color = fit_expected_color, linewidth = 0.6) +
geom_point(color = fit_expected_color, size = 1.8) +
labs(x = "Year", y = "CPUE CV") +
scale_x_continuous(limits = year_limits(cpue_tv_cv_input$Year), breaks = pretty_breaks()) +
scale_y_continuous(limits = c(0, NA), expand = expansion(mult = c(0, 0.05)))
```
```{r}
#| label: fig-sensitivity-cpue-tv-cv-fit
#| fig-cap: "CPUE index fit for Sensitivity 2, with time-varying CPUE CVs from CV_out_10.csv. Orange points are observed values and blue lines are expected values."
#| fig-width: 9
#| fig-height: 5
sensitivity_cpue_tv_cv_plot <- make_plot("Sensitivity 2 CPUE fit", function() {
if (!has_sensitivity_fit(sens_cpue_tv_cv)) {
stop("Run the Sensitivity 2 section in script_OMMP16.R to create ", sensitivity_results_path, call. = FALSE)
}
plot_cpue(data = sens_cpue_tv_cv$data, object = sens_cpue_tv_cv$obj)
})
sensitivity_cpue_tv_cv_plot$plot
```
## Sensitivity comparison
```{r}
#| label: fig-sensitivity-biomass-comparison
#| fig-cap: "Relative spawning-biomass comparison for the base OMMP16 fit and the two CPUE sensitivity fits."
#| fig-width: 9
#| fig-height: 5.5
sensitivity_biomass_plot <- make_plot("Sensitivity biomass comparison", function() {
if (!has_sensitivity_fit(sens_q2008) || !has_sensitivity_fit(sens_cpue_tv_cv)) {
stop("Run both CPUE sensitivity sections in script_OMMP16.R to create ", sensitivity_results_path, call. = FALSE)
}
plot_biomass_spawning(
data_list = list(data, sens_q2008$data, sens_cpue_tv_cv$data),
object_list = list(obj, sens_q2008$obj, sens_cpue_tv_cv$obj),
labels = c("Base OMMP16", "CPUE q split from 2008", "Time-varying CPUE CV"),
relative = TRUE
)
})
sensitivity_biomass_plot$plot
```
# Grid MCMC results
This section records the prototype grid configured in `script_OMMP16.R`. Its
historical output directory is not retained in this checkout, so the table
shows the expected legacy file layout rather than current results. It must not
be confused with the later accepted ESC31 nine-cell grid, diagnostics, balanced
2,000-draw posterior, and direct-\(M\) grid reported in
[the ESC31 grid page](../ESC31/4_grid.html).
```{r}
#| label: tbl-grid-mcmc-status
grid_values_report <- expand.grid(h = c(0.6, 0.7, 0.8), psi = c(1.5, 1.75, 2))
grid_mcmc_dirs <- list.dirs("grid_mcmc", recursive = FALSE, full.names = TRUE)
grid_mcmc_dir <- if (length(grid_mcmc_dirs)) {
grid_mcmc_dirs[which.max(file.info(grid_mcmc_dirs)$mtime)]
} else {
NA_character_
}
grid_mcmc_status <- grid_values_report |>
mutate(
Cell = row_number(),
`h start` = .data$h,
`psi start` = .data$psi,
`Expected output` = if (is.na(grid_mcmc_dir)) {
file.path("grid_mcmc", "RUN_ID", paste0("grid", .data$Cell, ".rda"))
} else {
file.path(grid_mcmc_dir, paste0("grid", .data$Cell, ".rda"))
},
Status = ifelse(
!is.na(grid_mcmc_dir) & file.exists(.data$`Expected output`),
"Historical output file found",
"Historical output not retained"
)
) |>
select(Cell, `h start`, `psi start`, `Expected output`, Status)
kable(
grid_mcmc_status |>
format_decimal_table(digits = 3) |>
replace_missing(),
align = c("r", "r", "r", "l", "l"),
caption = "Historical OMMP16 prototype grid cells and expected legacy output files. These are not the accepted ESC31 grid artifacts."
)
```
# M10 profile and selectivity sensitivity
The follow-up fixed-`M10` profile was run to identify which data and
selectivity components resist moving natural mortality at age 10 below 0.1. The
profile fixes `M10`, re-optimizes the remaining active maximum-likelihood
parameters, and then decomposes the objective into likelihood and prior
components. The profile minimum was near `M10 = 0.125`; values below 0.1 were
meaningfully worse even after removing the direct `M10` prior contribution.
The fixed-profile results are summarized in @tbl-m10-profile-summary, with
objective-component deltas in @tbl-m10-component-deltas and fishery-level
composition/selectivity deltas in @tbl-m10-fishery-deltas.
```{r}
#| label: read-m10-profile-results
#| include: false
m10_profile_summary <- read_csv("m10_profile_combined_summary.csv", show_col_types = FALSE)
m10_component_deltas <- read_csv("m10_profile_component_deltas_vs_min.csv", show_col_types = FALSE)
m10_fishery_deltas <- read_csv("m10_profile_fishery_deltas_0.100_vs_0.125.csv", show_col_types = FALSE)
m10_tuning_summary <- read_csv("m10_hyper_tuning_combined_summary.csv", show_col_types = FALSE)
```
```{r}
#| label: tbl-m10-profile-summary
m10_profile_table <- m10_profile_summary |>
filter(.data$status == "ok") |>
transmute(
`Fixed M10` = .data$M10_fixed,
`NLL` = .data$final_nll,
`Delta NLL` = .data$delta_nll,
`Delta NLL without M10 prior` = .data$delta_nll_without_m10_prior,
`Convergence code` = .data$convergence,
`Max gradient` = .data$max_gradient
)
kable(
m10_profile_table |>
format_decimal_table(digits = 3) |>
replace_missing(),
align = c("r", "r", "r", "r", "r", "r"),
caption = "Fixed-M10 profile summary. Delta NLL is measured relative to the best fixed-M10 profile point. Convergence code 0 indicates normal nlminb convergence; code 1 indicates the evaluation limit was reached."
)
```
The profile results in @tbl-m10-profile-summary show that `M10` is not being
held high by the direct `M10` prior. At `M10 = 0.100`, the total delta NLL was
about 1.82 and the delta NLL after removing the `M10` prior was about 1.89. At
`M10 = 0.085`, those deltas were about 4.35 and 4.38, respectively. The
resistance to lower `M10` therefore comes from the likelihood and selectivity
penalty rather than the `M10` prior.
```{r}
#| label: tbl-m10-component-deltas
m10_component_delta_table <- m10_component_deltas |>
filter(abs(.data$M10_fixed - 0.100) < 1e-8) |>
transmute(
Component = .data$component,
`Delta versus M10 = 0.125` = .data$delta_vs_min
) |>
arrange(desc(.data$`Delta versus M10 = 0.125`))
kable(
m10_component_delta_table |>
format_decimal_table(digits = 3) |>
replace_missing(),
align = c("l", "r"),
caption = "Objective-component deltas at fixed M10 = 0.100 relative to the profile minimum near M10 = 0.125. Positive values penalize lower M10; negative values favor lower M10."
)
```
The component decomposition in @tbl-m10-component-deltas points to the
selectivity penalty and composition fits as the main resistance to
`M10 = 0.100`. Tagging, POP, recruitment, CPUE index, and CPUE length
frequencies move slightly in the opposite direction and therefore do not explain
the high `M10`.
```{r}
#| label: tbl-m10-fishery-deltas
m10_fishery_delta_table <- m10_fishery_deltas |>
transmute(
Component = .data$component,
Fishery = .data$fishery,
`Delta versus M10 = 0.125` = .data$delta_vs_0125
) |>
arrange(desc(.data$`Delta versus M10 = 0.125`))
kable(
m10_fishery_delta_table |>
format_decimal_table(digits = 3) |>
replace_missing(),
align = c("l", "l", "r"),
caption = "Fishery-level selectivity and composition deltas at M10 = 0.100 relative to M10 = 0.125. Positive values penalize lower M10."
)
```
The fishery-level decomposition in @tbl-m10-fishery-deltas does not support
Indonesia selectivity as the main driver. The largest penalty is LL1
selectivity, followed by Indonesia age composition and LL1 length composition.
Indonesia selectivity itself slightly favors lower `M10`, and CPUE length
frequencies also slightly favor lower `M10`. CPUE selectivity is a secondary
positive penalty, but it is much smaller than LL1 selectivity.
```{r}
#| label: tbl-m10-hyper-tuning
m10_tuning_table <- m10_tuning_summary |>
arrange(.data$final_nll) |>
transmute(
Variant = .data$label,
`M10` = .data$M10_parameter,
`NLL` = .data$final_nll,
`Delta from baseline` = .data$delta_from_best_baseline,
`Max gradient` = .data$max_gradient,
`Below 0.1` = .data$below_0_1,
`LL1 rho year` = .data$LL1_rho_y,
`LL1 rho age` = .data$LL1_rho_a,
`LL1 sigma` = .data$LL1_sigma,
`CPUE sigma` = .data$CPUE_sigma,
`Indonesia sigma` = .data$Indonesia_sigma
)
kable(
m10_tuning_table |>
format_decimal_table(digits = 3) |>
replace_missing(),
align = c("l", rep("r", 10)),
caption = "Fixed selectivity-hyperparameter tuning runs used to test whether M10 can be reduced below 0.1. Delta from baseline is measured relative to the best baseline continuation in this tuning set."
)
```
The tuning results in @tbl-m10-hyper-tuning show that the best low-`M10`
fixed-hyperparameter run used LL1 `rho_y = rho_a = 0.98` and LL1 `sigma =
0.75`. That run produced `M10 = 0.089` after a longer continuation, so low
`M10` is numerically attainable. The cost is large: the objective was about 185
NLL units above the baseline continuation and the final maximum gradient was
still about 0.52. Adding CPUE or Indonesia to the hyperparameter tuning did not
improve the tradeoff. The specific longer LL1 run used for the selectivity
comparison is listed in @tbl-ll1-tuned-summary.
```{r}
#| label: make-ll1-tuned-object
#| include: false
make_object_from_parameter_list <- function(parameter_list) {
# The saved tuning object predates the explicit fixed overdispersion
# parameters. They were constants in that run, so restore the equivalent
# current parameter values without estimating them or changing the fit.
fixed_overdispersion <- c("pop_od", "hsp_od", "gt_od")
for (parameter_name in fixed_overdispersion) {
if (is.null(parameter_list[[parameter_name]])) {
parameter_list[[parameter_name]] <- parameters[[parameter_name]]
}
}
data_i <- data
data_i$priors <- get_priors(parameters = parameter_list)
object_i <- MakeADFun(
func = cmb(sbt_model, data_i),
parameters = parameter_list,
map = get_map(parameter_list),
silent = TRUE
)
object_i$env$tracemgc <- FALSE
object_i$env$inner.control$trace <- FALSE
object_i$par <- object_i$par
object_i$env$last.par.best <- object_i$par
object_i
}
ll1_tuned_results <- readRDS("m10_hyper_tuning_results_best_ll1_long.rds")
ll1_tuned_parameters <- ll1_tuned_results$details[[1]]$parameters
ll1_tuned_obj <- make_object_from_parameter_list(ll1_tuned_parameters)
ll1_tuned_summary <- ll1_tuned_results$summary |>
transmute(
Variant = .data$label,
`M10` = .data$M10_parameter,
`NLL` = .data$final_nll,
`Max gradient` = .data$max_gradient,
`LL1 rho year` = .data$LL1_rho_y,
`LL1 rho age` = .data$LL1_rho_a,
`LL1 sigma` = .data$LL1_sigma
)
```
```{r}
#| label: tbl-ll1-tuned-summary
kable(
ll1_tuned_summary |>
format_decimal_table(digits = 3) |>
replace_missing(),
align = c("l", rep("r", 6)),
caption = "Best low-M10 LL1 fixed-hyperparameter tuning case used for the selectivity comparison."
)
```
```{r}
#| label: fig-ll1-selectivity-tuning
#| fig-cap: "LL1 selectivity-at-age in the maximum-likelihood fit and the best low-M10 fixed-hyperparameter test. The low-M10 test fixes LL1 rho_y and rho_a at 0.98 and LL1 sigma at 0.75, producing M10 = 0.089 but with a large objective penalty."
#| fig-width: 15
#| fig-height: 8.5
baseline_ll1_selectivity <- plot_selectivity(
data = data,
object = obj,
fisheries = "LL1",
years = data$sel_LL1_yrs
) +
ggtitle("Maximum-likelihood fit")
tuned_ll1_selectivity <- plot_selectivity(
data = data,
object = ll1_tuned_obj,
fisheries = "LL1",
years = data$sel_LL1_yrs
) +
ggtitle("LL1 rho = 0.98, sigma = 0.75")
baseline_ll1_selectivity | tuned_ll1_selectivity
```
Visually, the low-`M10` LL1 test allows LL1 selectivity to move much more
freely across change years and ages. The lower `M10` is therefore not achieved
by a small, localized adjustment to Indonesia selectivity; it is achieved by
relaxing the LL1 selectivity penalty enough for the LL1 selectivity surface to
absorb tension that is otherwise expressed as higher age-10 natural mortality.
