---
title: "EBS Pollock FIMS Implementation"
execute:
echo: true
warning: false
message: false
format:
html:
embed-resources: true
lightbox: true
---
```{r setup}
#| include: false
suppressPackageStartupMessages({
library(dplyr)
library(tidyr)
library(ggplot2)
library(ggthemes)
library(gt)
library(FIMS)
})
theme_set(ggthemes::theme_few())
if (requireNamespace("FIMSdiags", quietly = TRUE)) {
library(FIMSdiags)
}
drop_all_na_cols <- function(data) {
data |>
select(where(~!all(is.na(.x))))
}
gt_report <- function(data, ...) {
data <- drop_all_na_cols(data)
table <- gt(data, ...)
integer_cols <- names(data)[tolower(names(data)) %in% c("year", "age")]
if (length(integer_cols) > 0) {
table <- table |>
fmt_number(columns = all_of(integer_cols), decimals = 0, use_seps = FALSE)
}
table
}
```
# Overview
This document summarizes a simplified FIMS implementation for EBS pollock based on inputs from the main 2024 EBS pollock stock assessment ([NPFMC SAFE 2024 EBS pollock assessment](https://files.npfmc.org/SAFE/2024/EBSpollock.pdf)). The model is a single-population catch-at-age assessment with one fishery fleet and four index series (BTS, ATS, AVO, and CPUE). Length compositions and seasonal timing are omitted.
The current FIMS configuration should be interpreted as a reduced implementation of the 2024 assessment rather than a one-for-one reproduction. Several assessment features are simplified or not yet represented here:
- In the 2024 assessment, age-1 observations from the acoustic-trawl survey and bottom trawl survey are treated as separate age-1 abundance indices, while the remaining ages are aggregated into biomass indices for those surveys. The FIMS implementation currently uses the survey biomass index series directly and handles age information through age-composition observations, rather than splitting age-1 survey signals into separate index fleets.
- The 2024 assessment links AVO selectivity to the ATS selectivity pattern. In this FIMS implementation, AVO is represented with fixed double-logistic selectivity parameters rather than sharing the estimated ATS selectivity curve.
- The 2024 assessment allows time variation in the BTS selectivity inflection point for ages 2 and older. The baseline FIMS run uses time-invariant logistic BTS selectivity.
- The 2024 assessment fishery selectivity is represented with non-parametric coefficients and regularity penalties. The baseline FIMS run uses a parametric logistic fishery selectivity curve, while the `TVselex` exploratory run uses a double-logistic fishery curve with annual variation only in the ascending inflection point.
- The 2024 assessment uses age-specific natural mortality. This FIMS implementation now fixes natural mortality by age, constant over years, at \(M_1 = 0.9\), \(M_2 = 0.45\), and \(M_{3+} = 0.3\).
- Other apparent simplifications include omitting length-composition data, omitting seasonal/subannual timing structure, using empirical stock-level weight-at-age inputs, and keeping recruitment variance treatment simpler than the production assessment configuration.
# Data Inputs
The index series include ATS (Acoustic-trawl survey), BTS (Bottom trawl survey), CPUE (catch per unit effort), and AVO (Acoustic Vessels of Opportunity; opportunistic backscatter data collection and processing).
```{r load-data}
doc_dir <- tryCatch(dirname(knitr::current_input()), error = function(e) NA_character_)
if (is.na(doc_dir) || doc_dir == ".") doc_dir <- getwd()
project_root <- if (file.exists(file.path(doc_dir, "data"))) doc_dir else file.path(doc_dir, "..")
project_root <- normalizePath(project_root, mustWork = FALSE)
input_path <- file.path(project_root, "data", "ebs_fims_data.rds")
if (!file.exists(input_path)) {
stop(
"Missing input data at: ",
normalizePath(input_path, mustWork = FALSE),
"\nRun: Rscript scripts/01_build_data.R"
)
}
payload <- readRDS(input_path)
data_ebs <- payload$data_ebs
years <- payload$years
ages <- payload$ages
```
The input summary is shown in @tbl-input-summary, and data availability by fleet and data type is shown in @tbl-data-availability.
```{r}
#| label: tbl-input-summary
#| tbl-cap: "Summary of the EBS pollock FIMS input data."
# Basic input summary
summary_tbl <- tibble::tibble(
item = c(
"Years",
"Ages",
"Fleets",
"Data types"
),
value = c(
paste0(min(years), "–", max(years), " (", length(years), ")"),
paste0(min(ages), "–", max(ages), " (", length(ages), ")"),
paste(sort(unique(data_ebs$name)), collapse = ", "),
paste(sort(unique(data_ebs$type)), collapse = ", ")
)
)
summary_tbl |>
gt_report() |>
tab_header(title = "FIMS Input Summary")
```
```{r}
#| label: tbl-data-availability
#| tbl-cap: "Available observations by fleet and data type."
# Data availability by fleet/type
availability <- data_ebs |>
filter(value != -999) |>
group_by(name, type) |>
summarize(
n_obs = n(),
years = paste0(min(timing, na.rm = TRUE), "–", max(timing, na.rm = TRUE)),
.groups = "drop"
)
availability |>
gt_report() |>
tab_header(title = "Data Availability by Fleet")
```
# Model Implementation
## Structure
- **Population**: single stock
- **Fleets**: fishery, BTS, ATS, AVO, CPUE
- **Recruitment**: Beverton-Holt with fixed-effect deviations and fixed recruitment log-SD
- **Natural mortality**: fixed by age and constant over years, with \(M_1 = 0.9\), \(M_2 = 0.45\), and \(M_{3+} = 0.3\)
- **Selectivity**: asymptotic logistic for fishery and BTS; estimated double logistic for ATS; fixed double logistic for AVO
- **Selectivity sharing**: CPUE shares fishery selectivity
- **Growth**: empirical WAA from stock-level SSB weights
- **Maturity**: logistic
## Run Definitions
Two model runs are defined in this implementation:
| Run | Description | Selectivity treatment |
|---|---|---|
| `run0` | Initial FIMS run used for the current fitted results. | Fishery and BTS use time-invariant logistic selectivity; ATS uses time-invariant double-logistic selectivity; AVO double-logistic selectivity is fixed. |
| `TVselex` | Time-varying fishery selectivity run. | Fishery selectivity is changed to a double-logistic form and only the ascending inflection point is expanded over model years; AVO double-logistic selectivity is fixed. |
The requested `TVselex` run is specified with the FIMS `DoubleLogisticSelectivity` module because the installed FIMS version used here (`r as.character(utils::packageVersion("FIMS"))`) exposes logistic and double-logistic selectivity classes, but not a double-normal selectivity class.
The `run0` parameter specifications are adjusted from the FIMS defaults before fitting as follows.
```{r}
#| label: lst-run0-parameter-specifications
#| eval: false
#| echo: true
#| code-fold: false
cfg <- create_default_configurations(data_4_model) |>
tidyr::unnest(cols = data) |>
dplyr::rows_update(
tibble::tibble(
module_name = "Selectivity",
fleet_name = c("fishery", "bts", "ats", "avo"),
module_type = c("Logistic", "Logistic", "DoubleLogistic", "DoubleLogistic")
),
by = c("module_name", "fleet_name")
) |>
tidyr::nest(.by = c(model_family, module_name, fleet_name))
pars <- create_default_parameters(cfg, data_4_model) |>
tidyr::unnest(cols = data) |>
dplyr::rows_update(
tibble::tibble(
module_name = "Recruitment",
label = c(
rep("log_devs", length((get_start_year(data_4_model) + 1):get_end_year(data_4_model))),
"log_sd"
),
time = c(
(get_start_year(data_4_model) + 1):get_end_year(data_4_model),
NA_real_
),
value = c(
rep(0, length((get_start_year(data_4_model) + 1):get_end_year(data_4_model))),
0.1
),
estimation_type = c(
rep("fixed_effects", length((get_start_year(data_4_model) + 1):get_end_year(data_4_model))),
"constant"
)
),
by = c("module_name", "label", "time")
) |>
dplyr::mutate(
value = dplyr::case_when(
module_name == "Population" & label == "log_M" & age == 1 ~ log(0.9),
module_name == "Population" & label == "log_M" & age == 2 ~ log(0.45),
module_name == "Population" & label == "log_M" & age >= 3 ~ log(0.3),
TRUE ~ value
),
estimation_type = dplyr::case_when(
module_name == "Population" & label == "log_M" ~ "constant",
TRUE ~ estimation_type
)
) |>
dplyr::rows_update(
tibble::tibble(
module_name = "Selectivity",
fleet_name = "avo",
label = c(
"inflection_point_asc",
"slope_asc",
"inflection_point_desc",
"slope_desc"
),
value = c(1.5, 2.0, 8.0, 0.1),
estimation_type = "constant"
),
by = c("module_name", "fleet_name", "label")
) |>
dplyr::mutate(
selectivity_shared_with = dplyr::case_when(
module_name == "Selectivity" & fleet_name == "cpue" ~ "fishery",
TRUE ~ NA_character_
)
)
```
The `TVselex` run starts from the same data and recruitment settings, but changes fishery selectivity to `DoubleLogistic` and makes only the ascending inflection point parameter time-varying by creating one row for that parameter in each model year. The remaining fishery double-logistic parameters are retained as time-invariant fixed effects.
```{r}
#| label: lst-tvselex-parameter-specifications
#| eval: false
#| echo: true
#| code-fold: false
cfg_tvselex <- create_default_configurations(data_4_model) |>
tidyr::unnest(cols = data) |>
dplyr::rows_update(
tibble::tibble(
module_name = "Selectivity",
fleet_name = c("fishery", "bts", "ats", "avo"),
module_type = c("DoubleLogistic", "Logistic", "DoubleLogistic", "DoubleLogistic")
),
by = c("module_name", "fleet_name")
) |>
tidyr::nest(.by = c(model_family, module_name, fleet_name))
pars_tvselex_base <- create_default_parameters(cfg_tvselex, data_4_model) |>
tidyr::unnest(cols = data) |>
dplyr::rows_update(
tibble::tibble(
module_name = "Recruitment",
label = c(
rep("log_devs", length((get_start_year(data_4_model) + 1):get_end_year(data_4_model))),
"log_sd"
),
time = c(
(get_start_year(data_4_model) + 1):get_end_year(data_4_model),
NA_real_
),
value = c(
rep(0, length((get_start_year(data_4_model) + 1):get_end_year(data_4_model))),
0.1
),
estimation_type = c(
rep("fixed_effects", length((get_start_year(data_4_model) + 1):get_end_year(data_4_model))),
"constant"
)
),
by = c("module_name", "label", "time")
) |>
dplyr::mutate(
value = dplyr::case_when(
module_name == "Population" & label == "log_M" & age == 1 ~ log(0.9),
module_name == "Population" & label == "log_M" & age == 2 ~ log(0.45),
module_name == "Population" & label == "log_M" & age >= 3 ~ log(0.3),
TRUE ~ value
),
estimation_type = dplyr::case_when(
module_name == "Population" & label == "log_M" ~ "constant",
TRUE ~ estimation_type
)
) |>
dplyr::rows_update(
tibble::tibble(
module_name = "Selectivity",
fleet_name = "avo",
label = c(
"inflection_point_asc",
"slope_asc",
"inflection_point_desc",
"slope_desc"
),
value = c(1.5, 2.0, 8.0, 0.1),
estimation_type = "constant"
),
by = c("module_name", "fleet_name", "label")
)
fishery_tvselex <- pars_tvselex_base |>
dplyr::filter(
module_name == "Selectivity",
fleet_name == "fishery",
label == "inflection_point_asc"
) |>
dplyr::select(-time) |>
tidyr::crossing(time = years) |>
dplyr::mutate(estimation_type = "fixed_effects")
pars_tvselex <- pars_tvselex_base |>
dplyr::filter(
!(
module_name == "Selectivity" &
fleet_name == "fishery" &
label == "inflection_point_asc"
)
) |>
dplyr::bind_rows(fishery_tvselex) |>
dplyr::mutate(
selectivity_shared_with = dplyr::case_when(
module_name == "Selectivity" & fleet_name == "cpue" ~ "fishery",
TRUE ~ NA_character_
)
)
input_tvselex <- initialize_fims(pars_tvselex, data_4_model)
```
# Results
## Fit to run0
The current fitted results are for `run0`. The fitted spawning biomass trajectory is shown in @fig-ssb. Fishery observed and predicted landings are compared in @fig-landings-fit, while survey and CPUE index fits are summarized in @fig-fits-by-fleet and residual patterns are shown in @fig-residuals. These index series are shown on their native arithmetic scales with separate facet scales because the magnitudes differ across surveys and CPUE, and the residuals are standardized by the input observation standard deviations on the arithmetic scale. Expected age-composition outputs for fishery, BTS, and ATS are shown in @fig-agecomp-fishery, @fig-agecomp-bts, and @fig-agecomp-ats, respectively, while selectivity is shown in @fig-selectivity-curves. The time-invariant selectivity assumptions used to initialize each gear are shown for selected years in @fig-selectivity-assumptions. A diagnostics overview, including model status and the fitted SSB and recruitment trajectories, is shown in @tbl-model-status and @fig-ssb-recruitment. The estimated stock-recruitment relationship for both fitted models is shown in @fig-stock-recruitment. Likelihood diagnostics and key parameter summaries are reported in @tbl-diagnostics, @tbl-recruitment-parameters, and @tbl-q.
```{r load-fit}
if (!exists("doc_dir")) {
doc_dir <- tryCatch(dirname(knitr::current_input()), error = function(e) NA_character_)
if (is.na(doc_dir) || doc_dir == ".") doc_dir <- getwd()
}
if (!exists("project_root")) {
project_root <- if (file.exists(file.path(doc_dir, "data"))) doc_dir else file.path(doc_dir, "..")
}
fit_path <- file.path(project_root, "outputs", "fims_fit_summary.rds")
if (!file.exists(fit_path)) {
stop(
"Missing model fit summary at: ",
normalizePath(fit_path, mustWork = FALSE),
"\nRun: Rscript scripts/02_fit_model.R"
)
}
fit <- readRDS(fit_path)
est <- fit$estimates
rep <- fit$report
# Map year index to actual years (robust to 0- or 1-based index)
map_year <- function(year_i) {
years[year_i]
}
```
```{r}
#| label: fig-ssb
#| fig-cap: "Estimated spawning biomass over time from the fitted EBS pollock FIMS model."
ssb <- est |>
filter(label == "spawning_biomass") |>
mutate(year = map_year(year_i))
ggplot(ssb, aes(x = year, y = estimated)) +
geom_line(color = "#1b9e77", linewidth = 0.8) +
scale_y_continuous(limits = c(0, NA)) +
labs(x = "Year", y = "SSB", title = "Spawning Biomass")
```
### Landings Fit
```{r}
#| label: fig-landings-fit
#| fig-cap: "Observed fishery landings and predicted fishery landings from the fitted model."
fishery_landings_expected <- rep$landings_expected[[which.max(vapply(rep$landings_expected, sum, numeric(1), na.rm = TRUE))]]
land <- tibble::tibble(
year = years,
expected = fishery_landings_expected,
observed = data_ebs |>
filter(type == "landings", name == "fishery") |>
arrange(timing) |>
pull(value)
)
ggplot(land, aes(x = year)) +
geom_line(aes(y = expected), color = "#377eb8", linewidth = 0.8) +
geom_point(aes(y = observed), color = "#377eb8", size = 1.2, alpha = 0.7) +
scale_y_continuous(limits = c(0, NA)) +
labs(x = "Year", y = "Landings", title = "Fishery Landings: Observed vs Expected")
```
### Fits by fleet
```{r}
#| label: fig-fits-by-fleet
#| fig-cap: "Observed and predicted index series for BTS, ATS, AVO, and CPUE shown on arithmetic scales with separate facet axes."
obs_index <- data_ebs |>
filter(type == "index", value != -999) |>
transmute(series = toupper(name), year = timing, observed = value, uncertainty = uncertainty)
pred_index <- purrr::map2_dfr(
rep$index_expected,
seq_along(rep$index_expected),
~tibble::tibble(series_id = .y, year = years, predicted = .x)
)
score_map <- function(pred_tbl, obs_tbl) {
obs_names <- sort(unique(obs_tbl$series))
pred_ids <- sort(unique(pred_tbl$series_id))
perms <- expand.grid(rep(list(pred_ids), length(obs_names)))
perms <- perms[apply(perms, 1, function(x) length(unique(x)) == length(x)), , drop = FALSE]
scores <- purrr::map_dfr(seq_len(nrow(perms)), function(i) {
perm <- as.integer(perms[i, ])
total_score <- 0
for (j in seq_along(obs_names)) {
obs_j <- obs_tbl |>
filter(series == obs_names[j]) |>
select(year, observed)
pred_j <- pred_tbl |>
filter(series_id == perm[j]) |>
select(year, predicted)
joined <- left_join(obs_j, pred_j, by = "year")
total_score <- total_score + sqrt(mean((joined$observed - joined$predicted)^2, na.rm = TRUE))
}
tibble::tibble(row_id = i, total_score = total_score)
})
best <- perms[scores$row_id[which.min(scores$total_score)], ]
tibble::tibble(series = obs_names, series_id = as.integer(best[1, ]))
}
index_map <- score_map(pred_index, obs_index)
fits <- pred_index |>
inner_join(index_map, by = "series_id") |>
left_join(obs_index, by = c("series", "year")) |>
mutate(series = factor(series, levels = c("BTS", "ATS", "AVO", "CPUE")))
ggplot(fits, aes(x = year)) +
geom_line(aes(y = predicted), color = "#377eb8", linewidth = 0.8) +
geom_point(aes(y = observed), color = "#377eb8", size = 1.1, alpha = 0.7, na.rm = TRUE) +
facet_wrap(~series, scales = "free_y") +
scale_y_continuous(limits = c(0, NA)) +
labs(x = "Year", y = "Index", title = "Survey and CPUE Fits on Arithmetic Scale")
```
### Residuals
```{r}
#| label: fig-residuals
#| fig-cap: "Standardized arithmetic residuals for BTS, ATS, AVO, and CPUE. Residuals are divided by the input observation standard deviation on the arithmetic scale."
log_sd_to_sd <- function(mu, log_sd) {
mu * sqrt(exp(log_sd^2) - 1)
}
resid_df <- fits |>
mutate(
obs_sd = log_sd_to_sd(observed, uncertainty),
std_residual = ifelse(is.finite(obs_sd) & obs_sd > 0, (observed - predicted) / obs_sd, NA_real_)
)
ylim_zero <- function(x) {
c(min(0, min(x, na.rm = TRUE)), max(0, max(x, na.rm = TRUE)))
}
ggplot(resid_df, aes(x = year, y = std_residual)) +
geom_hline(yintercept = 0, color = "gray50", linewidth = 0.4) +
geom_point(size = 1.1, alpha = 0.7) +
facet_wrap(~series, scales = "free_y") +
scale_y_continuous(limits = ylim_zero(resid_df$std_residual)) +
labs(x = "Year", y = "Standardized residual", title = "Survey and CPUE Standardized Residuals")
```
### Age composition fits
```{r}
#| label: fig-agecomp-fishery
#| fig-cap: "Expected fishery age-composition patterns, faceted by year."
#| fig-width: 8
#| fig-height: 11
#| out-width: "100%"
agecomp_counts <- vapply(rep$agecomp_expected, function(x) sum(x > 0, na.rm = TRUE), numeric(1))
agecomp_map <- c(
fishery = which(agecomp_counts == 900)[1],
bts = which(agecomp_counts == 630)[1],
ats = which(agecomp_counts == 285)[1]
)
normalize_agecomp <- function(data, value_col) {
value_col <- rlang::ensym(value_col)
data |>
group_by(year) |>
mutate(
.total = sum(!!value_col, na.rm = TRUE),
!!value_col := ifelse(.total > 0, !!value_col / .total, 0)
) |>
ungroup() |>
select(-.total)
}
make_agecomp_df <- function(series_id) {
tibble::tibble(
year = rep(years, each = length(ages)),
age = rep(ages, times = length(years)),
expected = rep$agecomp_expected[[series_id]]
) |>
group_by(year) |>
filter(any(expected > 0, na.rm = TRUE)) |>
ungroup() |>
normalize_agecomp(expected)
}
make_obs_agecomp_df <- function(fleet_name) {
data_ebs |>
filter(type == "age_comp", name == fleet_name, value != -999) |>
transmute(year = timing, age, observed = value) |>
normalize_agecomp(observed)
}
agecomp_fishery <- make_agecomp_df(agecomp_map[["fishery"]])
obs_agecomp_fishery <- make_obs_agecomp_df("fishery")
ggplot(agecomp_fishery, aes(x = age, y = expected)) +
geom_line(linewidth = 0.5, color = "#4daf4a") +
geom_point(
data = obs_agecomp_fishery,
aes(x = age, y = observed),
inherit.aes = FALSE,
color = "#1f78b4",
size = 0.9,
alpha = 0.8
) +
facet_wrap(~year, ncol = 5, dir = "v") +
scale_y_continuous(limits = c(0, NA)) +
labs(x = "Age", y = "Proportion", title = "Fishery Age Composition Fits")
```
```{r}
#| label: fig-agecomp-bts
#| fig-cap: "Observed and predicted BTS age-composition patterns, faceted by year."
#| fig-width: 8
#| fig-height: 11
#| out-width: "100%"
agecomp_bts <- make_agecomp_df(agecomp_map[["bts"]])
obs_agecomp_bts <- make_obs_agecomp_df("bts")
ggplot(agecomp_bts, aes(x = age, y = expected)) +
geom_line(linewidth = 0.5, color = "#4daf4a") +
geom_point(
data = obs_agecomp_bts,
aes(x = age, y = observed),
inherit.aes = FALSE,
color = "#1f78b4",
size = 0.9,
alpha = 0.8
) +
facet_wrap(~year, ncol = 4, dir = "v") +
scale_y_continuous(limits = c(0, NA)) +
labs(x = "Age", y = "Proportion", title = "BTS Age Composition Fits")
```
```{r}
#| label: fig-agecomp-ats
#| fig-cap: "Observed and predicted ATS age-composition patterns, faceted by year."
#| fig-width: 8
#| fig-height: 11
#| out-width: "100%"
agecomp_ats <- make_agecomp_df(agecomp_map[["ats"]])
obs_agecomp_ats <- make_obs_agecomp_df("ats")
ggplot(agecomp_ats, aes(x = age, y = expected)) +
geom_line(linewidth = 0.5, color = "#4daf4a") +
geom_point(
data = obs_agecomp_ats,
aes(x = age, y = observed),
inherit.aes = FALSE,
color = "#1f78b4",
size = 0.9,
alpha = 0.8
) +
facet_wrap(~year, ncol = 3, dir = "v") +
scale_y_continuous(limits = c(0, NA)) +
labs(x = "Age", y = "Proportion", title = "ATS Age Composition Fits")
```
### Selectivity curves
```{r}
#| label: fig-selectivity-curves
#| fig-cap: "Estimated selectivity-at-age curves by fleet. Curves are reconstructed from the fitted fleet-specific selectivity parameter blocks."
logistic_selectivity <- function(age, inflection_point, slope) {
1 / (1 + exp(-slope * (age - inflection_point)))
}
double_logistic <- function(age, inflection_point_asc, slope_asc, inflection_point_desc, slope_desc) {
asc <- 1 / (1 + exp(-slope_asc * (age - inflection_point_asc)))
desc <- 1 / (1 + exp(slope_desc * (age - inflection_point_desc)))
asc * desc
}
selectivity_template <- tibble::tribble(
~fleet, ~module_type,
"ats", "DoubleLogistic",
"avo", "DoubleLogistic",
"bts", "Logistic",
"fishery", "Logistic"
)
selectivity_parameters <- est |>
filter(module_name == "Selectivity") |>
distinct(parameter_id, module_type, label, estimated) |>
arrange(parameter_id) |>
mutate(block_id = cumsum(c(TRUE, diff(parameter_id) != 1L))) |>
group_by(block_id) |>
mutate(module_type = first(module_type)) |>
ungroup() |>
distinct(block_id, module_type, label, estimated) |>
group_by(block_id, module_type) |>
summarize(values = list(stats::setNames(estimated, label)), .groups = "drop") |>
left_join(
selectivity_template |>
mutate(block_id = row_number()),
by = c("block_id", "module_type")
)
selectivity_curves <- selectivity_parameters |>
rowwise() |>
mutate(
curve = list({
pars <- values
sel <- if (module_type == "Logistic") {
logistic_selectivity(ages, pars[["inflection_point"]], pars[["slope"]])
} else {
double_logistic(
ages,
pars[["inflection_point_asc"]],
pars[["slope_asc"]],
pars[["inflection_point_desc"]],
pars[["slope_desc"]]
)
}
tibble(age = ages, selectivity = sel)
})
) |>
ungroup() |>
select(fleet, curve) |>
unnest(curve)
ggplot(selectivity_curves, aes(x = age, y = selectivity)) +
geom_line(color = "#984ea3", linewidth = 0.8) +
facet_wrap(~fleet) +
scale_y_continuous(limits = c(0, NA)) +
labs(x = "Age", y = "Selectivity", title = "Selectivity at Age")
```
```{r}
#| label: fig-selectivity-assumptions
#| fig-cap: "Illustrative selectivity assumptions by gear for selected years. Fishery and BTS are initialized with asymptotic logistic selectivity, ATS uses an estimated double-logistic curve, and AVO uses a fixed double-logistic curve. In the current implementation these assumptions are time-invariant, so the curves are identical across years."
selected_years <- c(min(years), years[ceiling(length(years) / 2)], max(years))
sel_assumptions <- tidyr::crossing(
year = selected_years,
gear = c("fishery", "bts", "ats", "avo"),
age = ages
) |>
mutate(
selectivity = dplyr::case_when(
gear %in% c("fishery", "bts") ~ logistic_selectivity(
age = age,
inflection_point = 2,
slope = 1
),
TRUE ~ double_logistic(
age = age,
inflection_point_asc = dplyr::if_else(gear == "avo", 1.5, 2),
slope_asc = dplyr::if_else(gear == "avo", 2.0, 1),
inflection_point_desc = dplyr::if_else(gear == "avo", 8.0, 4),
slope_desc = dplyr::if_else(gear == "avo", 0.1, 1)
)
),
gear = factor(gear, levels = c("fishery", "bts", "ats", "avo"))
)
ggplot(sel_assumptions, aes(x = age, y = selectivity, color = factor(year))) +
geom_line(linewidth = 0.8) +
facet_wrap(~gear) +
scale_y_continuous(limits = c(0, 1)) +
labs(
x = "Age",
y = "Assumed selectivity",
color = "Year",
title = "Selected Years of Selectivity Assumptions by Gear"
)
```
## Time-varying selectivity run
```{r}
#| label: tvselex-status-data
#| include: false
tvselex_path <- file.path(project_root, "outputs", "tvselex_fit_summary.rds")
tvselex <- if (file.exists(tvselex_path)) readRDS(tvselex_path) else NULL
tvselex_status <- if (is.null(tvselex)) {
tibble::tibble(
metric = c("Run status", "Interpretation"),
value = c(
"not yet run",
"No saved TVselex output was found; run scripts/05_fit_tvselex.R to attempt this model."
)
)
} else {
tibble::tibble(
metric = c(
"Run status",
"FIMS version",
"Fixed effects",
"Random effects",
"Time-varying fishery selectivity rows",
"Max gradient",
"Total NLL",
"Terminal SB",
"Interpretation"
),
value = c(
tvselex$status,
as.character(tvselex$version),
as.character(tvselex$number_of_parameters[["fixed_effects"]]),
as.character(tvselex$number_of_parameters[["random_effects"]]),
as.character(tvselex$fishery_selectivity_rows),
ifelse(is.null(tvselex$max_gradient), "", sprintf("%.5f", tvselex$max_gradient)),
ifelse(is.null(tvselex$report$jnll), "", sprintf("%.2f", tvselex$report$jnll)),
ifelse(
is.null(tvselex$report$spawning_biomass),
"",
sprintf("%.2f", tail(tvselex$report$spawning_biomass[[1]], 1))
),
if (identical(tvselex$status, "fit completed")) {
if (is.finite(tvselex$max_gradient) && tvselex$max_gradient <= 0.01) {
"TVselex completed and appears well converged by gradient threshold."
} else {
"TVselex completed, but the maximum gradient is too high for the run to be treated as a reliable comparison model."
}
} else {
"TVselex did not produce an accepted FIMS output object; do not compare fitted trajectories or likelihood components from this attempt."
}
)
)
}
terminal_ssb <- function(report) {
if (is.null(report$spawning_biomass)) {
return("")
}
ssb_values <- if (is.list(report$spawning_biomass)) {
report$spawning_biomass[[1]]
} else {
report$spawning_biomass
}
sprintf("%.2f", tail(ssb_values, 1))
}
run0_status <- tibble::tibble(
metric = c(
"Run status",
"FIMS version",
"Fixed effects",
"Random effects",
"Time-varying fishery selectivity rows",
"Max gradient",
"Total NLL",
"Terminal SB",
"Interpretation"
),
run0 = c(
"fit completed",
as.character(fit$version),
as.character(fit$number_of_parameters[["fixed_effects"]]),
as.character(fit$number_of_parameters[["random_effects"]]),
"0",
sprintf("%.5f", fit$max_gradient),
ifelse(is.null(rep$jnll), "", sprintf("%.2f", rep$jnll)),
terminal_ssb(rep),
if (is.finite(fit$max_gradient) && fit$max_gradient <= 0.01) {
"run0 completed and appears well converged by gradient threshold."
} else {
"run0 completed, but the maximum gradient is too high for the run to be treated as a final model."
}
)
)
run_status_comparison <- run0_status |>
dplyr::full_join(
tvselex_status |>
dplyr::rename(TVselex = value),
by = "metric"
) |>
dplyr::mutate(
run0 = dplyr::coalesce(run0, ""),
TVselex = dplyr::coalesce(TVselex, "")
)
tvselex_fixed_effects <- tvselex_status$value[tvselex_status$metric == "Fixed effects"]
if (length(tvselex_fixed_effects) == 0) tvselex_fixed_effects <- "not available"
tvselex_fishery_rows <- tvselex_status$value[tvselex_status$metric == "Time-varying fishery selectivity rows"]
if (length(tvselex_fishery_rows) == 0) tvselex_fishery_rows <- "not available"
```
The `TVselex` run was attempted as a time-varying fishery selectivity alternative to `run0`. This model increases flexibility by estimating the fishery double-logistic ascending inflection point by year, while retaining time-invariant values for the other fishery selectivity parameters. In the current attempt, this expanded the model to `r tvselex_fixed_effects` fixed effects, including `r tvselex_fishery_rows` time-varying fishery selectivity rows.
The current `TVselex` formulation completed with a maximum gradient of `r tvselex_status$value[tvselex_status$metric == "Max gradient"]`, total NLL of `r tvselex_status$value[tvselex_status$metric == "Total NLL"]`, and terminal spawning biomass of `r tvselex_status$value[tvselex_status$metric == "Terminal SB"]`.
```{r}
#| label: tbl-tvselex-status
#| tbl-cap: "Status summary comparing run0 and the attempted TVselex model run."
run_status_comparison |>
gt_report() |>
tab_header(title = "Run Status Comparison")
```
```{r}
#| label: tbl-tvselex-gradient-diagnostics
#| tbl-cap: "Largest absolute gradients from the TVselex TMB objective, matched back to the fixed-effect parameter table."
if (!is.null(tvselex) && !is.null(tvselex$gradient_diagnostics)) {
tvselex$gradient_diagnostics |>
slice_max(abs_gradient, n = 15, with_ties = FALSE) |>
transmute(
parameter_index,
module = module_name,
fleet = dplyr::coalesce(fleet_name, ""),
label,
year = time,
estimate = estimated,
gradient,
abs_gradient
) |>
gt_report() |>
fmt_number(columns = any_of(c("estimate", "gradient", "abs_gradient")), decimals = 3) |>
tab_header(title = "Largest TVselex Gradients")
} else {
tibble::tibble(note = "Gradient diagnostics were not saved with the current TVselex output.") |>
gt_report() |>
tab_header(title = "Largest TVselex Gradients")
}
```
```{r}
#| label: fig-tvselex-selectivity-over-time
#| fig-cap: "Fishery selectivity-at-age over time from the TVselex run. The ascending inflection point varies annually; the remaining fishery double-logistic parameters are time-invariant."
if (!is.null(tvselex) && identical(tvselex$status, "fit completed")) {
tv_sel_est <- tvselex$estimates |>
filter(module_name == "Selectivity", module_type == "DoubleLogistic")
tv_inflection_asc <- tv_sel_est |>
filter(label == "inflection_point_asc", parameter_id >= 545) |>
arrange(parameter_id) |>
transmute(year = years[seq_len(n())], inflection_point_asc = estimated)
tv_fishery_static <- tv_sel_est |>
filter(parameter_id %in% c(546, 547, 548)) |>
select(label, estimated) |>
distinct() |>
tibble::deframe()
tvselex_selectivity <- tv_inflection_asc |>
tidyr::crossing(age = ages) |>
mutate(
selectivity = double_logistic(
age = age,
inflection_point_asc = inflection_point_asc,
slope_asc = tv_fishery_static[["slope_asc"]],
inflection_point_desc = tv_fishery_static[["inflection_point_desc"]],
slope_desc = tv_fishery_static[["slope_desc"]]
)
)
ggplot(tvselex_selectivity, aes(x = year, y = age, fill = selectivity)) +
geom_tile() +
scale_fill_viridis_c(limits = c(0, 1), name = "Selectivity") +
scale_y_continuous(breaks = ages) +
labs(
x = "Year",
y = "Age",
title = "TVselex Fishery Selectivity Over Time"
)
} else {
plot.new()
text(
0.5, 0.5,
"TVselex fit output not available.\nRun scripts/05_fit_tvselex.R to populate this figure.",
cex = 1
)
}
```
## Diagnostics
This section summarizes basic model-status diagnostics for the fitted models and provides a place to evaluate a 5-peel retrospective case when retrospective output is available.
```{r}
#| label: tbl-model-status
#| tbl-cap: "Summary diagnostics for the current fitted model."
model_status_tbl <- tibble::tibble(
metric = c(
"FIMS version",
"Max gradient",
"Fixed effects",
"Random effects",
"Retrospective output",
"Interpretation"
),
value = c(
as.character(fit$version),
sprintf("%.4f", fit$max_gradient),
as.character(fit$number_of_parameters[["fixed_effects"]]),
as.character(fit$number_of_parameters[["random_effects"]]),
paste(
"run0",
ifelse(file.exists(file.path(project_root, "outputs", "retro_5_peel_summary.rds")), "available", "not yet run"),
"; TVselex",
ifelse(file.exists(file.path(project_root, "outputs", "tvselex_retro_5_peel_summary.rds")), "available", "not yet run")
),
ifelse(
fit$max_gradient <= 0.01,
"Mode appears well converged by gradient threshold.",
"Mode should be treated as provisional; gradient is larger than a typical convergence target."
)
)
)
model_status_tbl |>
gt_report() |>
tab_header(title = "Model Status Diagnostics")
```
```{r}
#| label: fig-ssb-recruitment
#| fig-cap: "Fitted spawning biomass and expected recruitment trajectories from run0 and TVselex, with the ADMB 2024 derived output overlaid for comparison. FIMS SSB is divided by 2 before plotting."
make_diag_series <- function(estimates, model_name) {
recruitment_values <- estimates |>
filter(label == "expected_recruitment") |>
mutate(year = map_year(year_i)) |>
transmute(model = model_name, source = "FIMS", year, value = estimated, quantity = "Recruitment")
ssb_values <- estimates |>
filter(label == "spawning_biomass") |>
mutate(year = map_year(year_i)) |>
transmute(model = model_name, source = "FIMS", year, value = estimated / 2, quantity = "SSB")
bind_rows(ssb_values, recruitment_values)
}
diag_series <- make_diag_series(est, "run0")
if (!is.null(tvselex) && identical(tvselex$status, "fit completed")) {
diag_series <- bind_rows(
diag_series,
make_diag_series(tvselex$estimates, "TVselex")
)
}
derived_dir <- file.path(dirname(project_root), "results", "derived")
read_derived_diag_series <- function(file_name, model_name) {
path <- file.path(derived_dir, file_name)
if (!file.exists(path)) return(NULL)
derived <- readRDS(path)
if (is.null(derived$timeseries)) return(NULL)
derived$timeseries |>
filter(quantity %in% c("SSB", "Recruit", "Recruit_age1")) |>
transmute(
model = model_name,
source = "ADMB 2024",
year = as.integer(year),
value = value,
quantity = dplyr::case_when(
quantity %in% c("Recruit", "Recruit_age1") ~ "Recruitment",
TRUE ~ quantity
)
)
}
derived_diag_series <- bind_rows(
read_derived_diag_series("admb_2024.rds", "ADMB 2024")
)
if (nrow(derived_diag_series) > 0) {
diag_series <- bind_rows(diag_series, derived_diag_series)
}
model_colors <- c(
"run0" = "#1b9e77",
"TVselex" = "#d95f02",
"ADMB 2024" = "#1f78b4"
)
ggplot(diag_series, aes(x = year, y = value, color = model, group = model)) +
geom_line(data = filter(diag_series, source == "FIMS"), linewidth = 0.9) +
geom_line(data = filter(diag_series, source == "ADMB 2024"), linewidth = 0.9, linetype = "22") +
geom_point(data = filter(diag_series, source == "ADMB 2024"), size = 0.75, alpha = 0.65) +
facet_wrap(~quantity, scales = "free_y", ncol = 1) +
scale_color_manual(values = model_colors) +
scale_y_continuous(limits = c(0, NA)) +
labs(
x = "Year",
y = NULL,
color = "Model",
title = "Spawning Biomass and Recruitment by Model"
)
```
```{r}
#| label: fig-stock-recruitment
#| fig-cap: "Estimated stock-recruitment relationship for run0 and TVselex. Text labels show model years, with spawning biomass on the x-axis and expected recruitment on the y-axis."
make_stock_recruitment <- function(estimates, model_name) {
ssb_values <- estimates |>
filter(label == "spawning_biomass") |>
transmute(year = map_year(year_i), spawning_biomass = estimated)
recruitment_values <- estimates |>
filter(label == "expected_recruitment") |>
transmute(year = map_year(year_i), recruitment = estimated)
inner_join(ssb_values, recruitment_values, by = "year") |>
mutate(model = model_name)
}
stock_recruitment <- make_stock_recruitment(est, "run0")
if (!is.null(tvselex) && identical(tvselex$status, "fit completed")) {
stock_recruitment <- bind_rows(
stock_recruitment,
make_stock_recruitment(tvselex$estimates, "TVselex")
)
}
ggplot(
stock_recruitment,
aes(
x = spawning_biomass,
y = recruitment,
label = year,
color = model
)
) +
geom_text(size = 2.6, alpha = 0.85) +
scale_x_continuous(limits = c(0, NA)) +
scale_y_continuous(limits = c(0, NA)) +
labs(
x = "Spawning biomass",
y = "Expected recruitment",
color = "Model",
title = "Estimated Stock-Recruitment Relationship"
)
```
```{r}
#| label: retro-5-peel-data
#| include: false
retro_path <- file.path(project_root, "outputs", "retro_5_peel_summary.rds")
tvselex_retro_path <- file.path(project_root, "outputs", "tvselex_retro_5_peel_summary.rds")
retro_inputs <- list()
add_retro_run <- function(retro, run_name) {
if ("run" %in% names(retro)) {
retro |>
mutate(run = dplyr::coalesce(.data$run, run_name))
} else {
retro |>
mutate(run = run_name)
}
}
if (file.exists(retro_path)) {
retro_inputs$run0 <- add_retro_run(readRDS(retro_path), "run0")
}
if (file.exists(tvselex_retro_path)) {
retro_inputs$TVselex <- add_retro_run(readRDS(tvselex_retro_path), "TVselex")
}
if (length(retro_inputs) > 0) {
retro <- bind_rows(retro_inputs)
if ("status" %in% names(retro)) {
retro <- retro |>
filter(is.na(status) | status == "fit completed")
}
} else {
retro <- tibble::tibble()
}
plot_retro_model <- function(retro_data, model_name) {
retro_model <- retro_data |>
filter(run == model_name)
if (nrow(retro_model) > 0) {
ggplot(retro_model, aes(x = year, y = value, color = factor(peel))) +
geom_line(linewidth = 0.8) +
facet_wrap(vars(quantity), scales = "free_y", ncol = 1) +
scale_y_continuous(limits = c(0, NA)) +
labs(
x = "Year",
y = NULL,
color = "Peel",
title = paste("Five-Peel Retrospective Diagnostics:", model_name)
)
} else {
plot.new()
text(
0.5, 0.5,
paste0(
"Five-peel retrospective output not found for ",
model_name,
".\nFailed peels are omitted from this figure."
),
cex = 1
)
}
}
```
```{r}
#| label: fig-retro-run0
#| fig-cap: "Five-peel retrospective trajectories for spawning biomass and recruitment from the run0 configuration. Failed peels are omitted from the plotted trajectories."
plot_retro_model(retro, "run0")
```
```{r}
#| label: fig-retro-tvselex
#| fig-cap: "Five-peel retrospective trajectories for spawning biomass and recruitment from the TVselex configuration. Failed peels are omitted from the plotted trajectories."
plot_retro_model(retro, "TVselex")
```
```{r}
#| label: tbl-diagnostics
#| tbl-cap: "Likelihood diagnostics by model component."
diag_tbl <- est |>
filter(!is.na(likelihood)) |>
group_by(module_name, label) |>
summarize(
n = n(),
mean_ll = mean(likelihood, na.rm = TRUE),
.groups = "drop"
)
diag_tbl |>
gt_report() |>
tab_header(title = "Likelihood Diagnostics by Component")
```
# SparseNUTS Scaffold
The repository includes a first-pass scaffold for Bayesian sampling with `SparseNUTS` using the existing FIMS/TMB objective. The main entry points are `scripts/sparsenuts_framework.R`, which builds the FIMS objects and exposes helper functions for extracting the TMB objective, and `scripts/03_fit_sparse_nuts.R`, which runs a conditional mode fit followed by `SparseNUTS` sampling.
The intended MCMC workflow is to use the maximum-likelihood fit as an initialization and preconditioning step, then sample the joint posterior for estimable fixed effects with `SparseNUTS`. In this application, that would allow posterior summaries for quantities such as spawning biomass, recruitment, fishing mortality, selectivity parameters, and survey catchability, while retaining the same core FIMS/TMB model structure used for the deterministic fit.
At present this should be treated as an implementation scaffold rather than a validated Bayesian assessment workflow. The code path for building the model object and passing it to `SparseNUTS` is in place, but full end-to-end sampling has not yet been verified for this application, and posterior diagnostics such as divergences, effective sample size, split-\(\hat{R}\), and prior sensitivity have not yet been evaluated. In practice, additional work will likely be needed on parameter blocking, scaling, priors, and possibly on how random effects or transformed parameters are exposed to the sampler.
Once operational, the expected MCMC outputs would include posterior draws saved to `outputs/sparsenuts_fit.rds`, from which trace plots, marginal posterior summaries, posterior intervals for SSB and recruitment, and posterior predictive checks could be added to this report.
An example command is shown below.
```sh
Rscript scripts/03_fit_sparse_nuts.R \
--num-samples=250 \
--num-warmup=250 \
--chains=4 \
--cores=1 \
--metric=diag \
--seed=123
```
By default this writes `outputs/sparsenuts_fit.rds`.
# References
Ianelli, J. N., Fissel, B., Holsman, K., Honkalehto, T., Kotwicki, S., Monnahan, C., Siddon, E., and Stienessen, S. 2024. *Assessment of the Walleye Pollock Stock in the Eastern Bering Sea*. North Pacific Fishery Management Council SAFE report. <https://files.npfmc.org/SAFE/2024/EBSpollock.pdf>.
NOAA-FIMS. 2026. *Fisheries Integrated Modeling System (FIMS)*. Version `r fit$version`. GitHub organization and FIMS development repository. Accessed April 24, 2026. <https://github.com/NOAA-FIMS>.
Monnahan, C. *SparseNUTS: Sparse No-U-Turn MCMC Sampling for Template Model Builder*. R package version `r if (requireNamespace("SparseNUTS", quietly = TRUE)) as.character(utils::packageVersion("SparseNUTS")) else "not installed"`. Accessed April 24, 2026. <https://noaa-afsc.github.io/SparseNUTS>.
# Tables
This section contains result tables. Input-only data tables are provided in the Appendix.
```{r}
#| label: tbl-recruitment-parameters
#| tbl-cap: "Estimated recruitment parameters retained in the fitted model summary."
param_tbl <- est |>
filter(label %in% c("log_rzero", "logit_steep")) |>
transmute(run = "run0", module_name, label, estimated) |>
distinct()
if (!is.null(tvselex) && identical(tvselex$status, "fit completed")) {
param_tbl <- bind_rows(
param_tbl,
tvselex$estimates |>
filter(label %in% c("log_rzero", "logit_steep")) |>
transmute(run = "TVselex", module_name, label, estimated) |>
distinct()
)
}
param_tbl |>
tidyr::pivot_wider(names_from = run, values_from = estimated) |>
gt_report() |>
fmt_number(columns = any_of(c("run0", "TVselex")), decimals = 3) |>
tab_header(title = "Key Recruitment Parameters")
```
```{r}
#| label: tbl-q
#| tbl-cap: "Estimated survey catchability coefficients."
fleet_module_map <- tibble::tibble(
module_id = 1:5,
fleet = c("ats", "avo", "bts", "cpue", "fishery")
)
q_tbl <- est |>
filter(label == "log_q") |>
mutate(q = exp(estimated)) |>
transmute(run = "run0", module_id, q)
if (!is.null(tvselex) && identical(tvselex$status, "fit completed")) {
q_tbl <- bind_rows(
q_tbl,
tvselex$estimates |>
filter(label == "log_q") |>
mutate(q = exp(estimated)) |>
transmute(run = "TVselex", module_id, q)
)
}
q_tbl |>
left_join(fleet_module_map, by = "module_id") |>
select(module_id, fleet, run, q) |>
tidyr::pivot_wider(names_from = run, values_from = q) |>
arrange(module_id) |>
gt_report() |>
fmt_number(columns = any_of(c("run0", "TVselex")), decimals = 3) |>
tab_header(title = "Survey Catchability (q)")
```
```{r}
#| label: tbl-naa-comparison
#| tbl-cap: "Initial and terminal estimated numbers-at-age by model run."
make_naa_table <- function(estimates, run_name) {
estimates |>
filter(label == "numbers_at_age") |>
transmute(
year = map_year(year_i),
age = age_i,
run = run_name,
value = estimated
) |>
group_by(run) |>
filter(year %in% range(year, na.rm = TRUE)) |>
ungroup() |>
mutate(period = ifelse(year == min(year, na.rm = TRUE), "initial", "terminal")) |>
select(period, age, run, value)
}
naa_comparison <- make_naa_table(est, "run0")
if (!is.null(tvselex) && identical(tvselex$status, "fit completed")) {
naa_comparison <- bind_rows(
naa_comparison,
make_naa_table(tvselex$estimates, "TVselex")
)
}
naa_comparison |>
tidyr::pivot_wider(names_from = run, values_from = value) |>
arrange(factor(period, levels = c("initial", "terminal")), age) |>
gt_report() |>
fmt_number(columns = any_of(c("run0", "TVselex")), decimals = 3) |>
tab_header(title = "Initial and Terminal Numbers-at-Age")
```
```{r}
#| label: tbl-catch-at-age
#| tbl-cap: "Estimated terminal-year fishery catch-at-age by model run."
make_catch_table <- function(estimates, run_name) {
estimates |>
filter(label == "landings_numbers_at_age") |>
mutate(year = map_year(year_i)) |>
filter(year == max(year, na.rm = TRUE)) |>
transmute(
age = age_i,
run = run_name,
value = dplyr::coalesce(estimated, expected)
) |>
group_by(age, run) |>
summarize(value = sum(value, na.rm = TRUE), .groups = "drop")
}
catch_comparison <- make_catch_table(est, "run0")
if (!is.null(tvselex) && identical(tvselex$status, "fit completed")) {
catch_comparison <- bind_rows(
catch_comparison,
make_catch_table(tvselex$estimates, "TVselex")
)
}
catch_comparison |>
tidyr::pivot_wider(names_from = run, values_from = value) |>
arrange(age) |>
gt_report() |>
fmt_number(columns = any_of(c("run0", "TVselex")), decimals = 3) |>
tab_header(title = "Terminal-Year Fishery Catch-at-Age")
```
# Appendix
## Model Equations and Assumptions {#sec-model-equations}
The fitted runs use a standard catch-at-age model. Key fixed assumptions and structural simplifications are summarized in @tbl-model-assumptions.
```{r}
#| label: tbl-model-assumptions
#| tbl-cap: "Key model assumptions and fixed parameters in the current FIMS implementation."
assumptions_tbl <- tibble::tibble(
component = c(
"Natural mortality",
"Recruitment",
"Recruitment variability",
"Fishery selectivity",
"BTS selectivity",
"ATS selectivity",
"AVO selectivity",
"CPUE selectivity",
"Growth",
"Maturity",
"Length composition",
"Seasonal timing",
"Landings uncertainty"
),
current_fims_setting = c(
"Age- and year-indexed `log_M` rows are created; values are fixed over years at `M1 = 0.9`, `M2 = 0.45`, and `M3+ = 0.3`.",
"Beverton-Holt stock-recruitment with estimated `log_rzero`; steepness fixed at `h = 0.6875`.",
"Annual recruitment deviations are fixed effects for 1965-2024; recruitment `log_sd` fixed at 0.1.",
"`run0` uses time-invariant logistic selectivity; `TVselex` uses double-logistic selectivity with annual ascending inflection points.",
"Time-invariant logistic selectivity.",
"Time-invariant double-logistic selectivity.",
"Fixed double-logistic selectivity: ascending inflection 1.5, ascending slope 2.0, descending inflection 8.0, descending slope 0.1.",
"Shares fishery selectivity.",
"Empirical weight-at-age from ATS WAA, filled by age-specific means where needed.",
"Fixed maturity-at-age schedule.",
"Omitted.",
"Omitted; annual time step only.",
"Fixed log-SD of 0.05."
),
contrast_with_2024_assessment = c(
"Matches the production-assessment feature of age-specific natural mortality more closely than the previous scalar-M implementation.",
"Same general stock-recruitment family, but this implementation keeps the variance treatment simpler.",
"Simpler than the production assessment treatment of recruitment variability and uncertainty.",
"Production model uses non-parametric fishery selectivity coefficients with regularity penalties.",
"Production model allows time variation in the BTS selectivity inflection point for ages 2 and older.",
"Closer to the production model structure than the other survey selectivity simplifications.",
"Production model links AVO selectivity to ATS; this implementation fixes AVO independently.",
"Simplified sharing assumption retained in this implementation.",
"Uses empirical stock-level inputs rather than a full growth process.",
"Uses a fixed schedule rather than estimating maturity.",
"Production assessment includes length information; this implementation does not.",
"Production assessment includes more detailed timing structure.",
"Simplified placeholder uncertainty."
)
)
assumptions_tbl |>
gt_report() |>
fmt_markdown(columns = everything()) |>
tab_header(title = "Key FIMS Model Assumptions and Fixed Parameters")
```
Recruitment (Beverton-Holt):
$$
R_t = \frac{0.8 R_0 h S_{t-1}}{0.2 R_0 \phi_0 (1-h) + S_{t-1}(h-0.2)} \exp(\epsilon_t)
$$
with deviations \(\epsilon_t\) estimated as fixed effects.
Numbers-at-age dynamics:
$$
N_{a,t} = N_{a-1,t-1} \exp(-Z_{a-1,t-1})
$$
with a plus-group at the maximum age.
Initial numbers at age:
$$
N_{a,1} = R_0 \exp\left(-\sum_{j=1}^{a-1} M_j\right), \quad
N_{A,1} = \frac{R_0 \exp\left(-\sum_{j=1}^{A-1} M_j\right)}{1 - \exp(-M_A)}
$$
Total mortality:
$$
Z_{a,t} = M_a + \sum_f F_{f,t} s_{f,a}
$$
where \(M_a\) is fixed at 0.9 for age 1, 0.45 for age 2, and 0.3 for ages 3 and older, and \(s_{f,a}\) is fleet selectivity at age.
Spawning biomass:
$$
SSB_t = \sum_a N_{a,t} w_{a} m_{a}
$$
Likelihood components:
Lognormal (landings and indices):
$$
\ell_{\text{lnorm}} = \sum_t \left[-\frac{\left(\log y_t - \log \hat{y}_t\right)^2}{2\sigma_t^2} - \log y_t - \log(\sigma_t\sqrt{2\pi})\right]
$$
Multinomial (age compositions):
$$
\ell_{\text{mult}} = \sum_{t} \log\left(\frac{N_t!}{\prod_a n_{a,t}!}\right) + \sum_t \sum_a n_{a,t} \log p_{a,t}
$$
where \(y_t\) and \(\hat{y}_t\) are observed and expected landings or index, \(\sigma_t\) is the lognormal SD, \(n_{a,t}\) are age-composition counts (or effective sample sizes), and \(p_{a,t}\) are predicted age proportions.
## Data Tables
Supporting input data tables are provided in @tbl-weight-at-age, @tbl-maturity, and @tbl-survey-agecomp. Result tables are reported in the main Tables section.
```{r}
#| label: tbl-weight-at-age
#| tbl-cap: "Input weight-at-age values used in the most recent model year."
waa_tbl <- data_ebs |>
filter(type == "weight_at_age") |>
filter(timing == max(timing, na.rm = TRUE)) |>
select(age, value) |>
arrange(age)
waa_tbl |>
gt_report() |>
tab_header(title = "Weight-at-Age Input (Most Recent Year)")
```
```{r}
#| label: tbl-maturity
#| tbl-cap: "Input maturity-at-age schedule used in the model."
mat_tbl <- tibble::tibble(
age = ages,
value = c(
0, 0.008, 0.289, 0.641, 0.842,
0.901, 0.947, 0.963, 0.97, 1,
1, 1, 1, 1, 1
)
) |>
arrange(age)
mat_tbl |>
gt_report() |>
fmt_number(columns = value, decimals = 3) |>
tab_header(title = "Maturity-at-Age Input")
```
```{r}
#| label: tbl-survey-agecomp
#| tbl-cap: "Observed survey age-composition input proportions in the most recent year with data."
survey_agecomp <- data_ebs |>
filter(type == "age_comp", name %in% c("bts", "ats")) |>
filter(value != -999) |>
group_by(name) |>
filter(timing == max(timing, na.rm = TRUE)) |>
ungroup() |>
transmute(fleet = name, year = timing, age, proportion = value) |>
group_by(fleet, year) |>
mutate(proportion = proportion / sum(proportion, na.rm = TRUE)) |>
ungroup() |>
select(-year) |>
arrange(fleet, age)
survey_agecomp |>
gt_report(groupname_col = "fleet") |>
tab_header(title = "Survey Age-Composition Inputs (Most Recent Year)")
```
## Parameter Estimates
The full `run0` parameter listing is shown in @tbl-parameter-appendix. A standard-deviation column is included for parameter standard errors, but the current saved FIMS summary does not populate finite standard errors; blank entries indicate that uncertainty output was not available for that parameter in the current saved fit.
```{r}
#| label: tbl-parameter-appendix
#| tbl-cap: "Appendix table of fitted run0 parameter estimates and standard deviations, where available."
parameter_appendix <- est |>
filter(estimation_type %in% c("fixed_effects", "random_effects")) |>
distinct(
module_name,
module_type,
fleet,
label,
year_i,
age_i,
parameter_id,
estimated,
uncertainty,
estimation_type
) |>
mutate(
fleet = dplyr::coalesce(fleet, ""),
year = dplyr::if_else(!is.na(year_i), map_year(year_i), NA_real_),
age = age_i,
sd = dplyr::if_else(is.finite(uncertainty), uncertainty, NA_real_)
) |>
select(
parameter_id,
module_name,
module_type,
fleet,
label,
year,
age,
estimation_type,
estimate = estimated,
sd
) |>
arrange(module_name, fleet, label, year, age, parameter_id)
parameter_appendix |>
gt_report() |>
fmt_number(columns = any_of(c("estimate", "sd")), decimals = 3) |>
tab_header(title = "Appendix: run0 Parameter Estimates")
```
# Notes
- AVO selectivity is fixed independently with ascending inflection point 1.5, ascending slope 2.0, descending inflection point 8.0, and descending slope 0.1.
- ATS weight-at-age is used as the growth input.
- Composition likelihoods use a multinomial distribution with sample sizes taken from the 2024 EBS pollock assessment inputs.
