An evaluation of ABC and TAC projections for the Bering Sea and Aleutian Island groundfish

library(here)
library(tidyverse)
library(scales)
library(GGally)

theme_set(theme_minimal(base_size = 12))

data_root <- if (file.exists(file.path("..", "data", "BSAI_OFL_ABC_TAC.csv"))) {
  file.path("..", "data")
} else {
  "data"
}

df <- read_csv(file.path(data_root, "BSAI_OFL_ABC_TAC.csv"), show_col_types = FALSE)

#goaraw <- read_csv(file.path(data_root, "GOA_OFL_ABC_TAC.csv"), show_col_types = FALSE)
goaraw <- read_csv(file.path(data_root, "GOA_OFL_ABC_TAC_specs.csv"), show_col_types = FALSE)

mainspp <- c(
  "Pollock",
  "Yellowfin sole",
  "Pacific cod",
  "Atka mackerel",
  "Northern rock sole",
  "Flathead sole",
  "Pacific ocean perch"
)

dfOY <- df %>%
  filter(OY == 1) %>%
  group_by(Year = ProjYear, lag, Species, Order) %>%
  summarise(
    ABC = sum(ABC, na.rm = TRUE),
    OFL = sum(OFL, na.rm = TRUE),
    TAC = sum(TAC, na.rm = TRUE),
    .groups = "drop"
  )

species_order <- dfOY %>% distinct(Species, Order)

main_share <- dfOY %>%
  filter(lag == 1) %>%
  summarise(share = sum(ABC[Species %in% mainspp], na.rm = TRUE) / sum(ABC, na.rm = TRUE)) %>%
  pull(share)

cv_tbl <- dfOY %>%
  filter(lag == 1, Species %in% mainspp) %>%
  group_by(Species) %>%
  summarise(
    cv_ABC = sd(ABC, na.rm = TRUE) / mean(ABC, na.rm = TRUE),
    cv_TAC = sd(TAC, na.rm = TRUE) / mean(TAC, na.rm = TRUE),
    .groups = "drop"
  )

cv_max <- cv_tbl %>% arrange(desc(cv_ABC)) %>% slice(1)
cv_min_tac <- cv_tbl %>% arrange(cv_TAC) %>% slice(1)

goaOY <- goaraw %>%
  filter(OY == 1, Area == "Total") %>%
  group_by(Year = ProjYear, lag, Species) %>%
  summarise(
    ABC = sum(ABC, na.rm = TRUE),
    OFL = sum(OFL, na.rm = TRUE),
    TAC = sum(TAC, na.rm = TRUE),
    .groups = "drop"
  )

goa_main_spp <- goaOY %>%
  filter(lag == 1) %>%
  group_by(Species) %>%
  summarise(mean_ABC = mean(ABC, na.rm = TRUE), .groups = "drop") %>%
  arrange(desc(mean_ABC)) %>%
  slice_head(n = 7) %>%
  pull(Species)

goa_cv_tbl <- goaOY %>%
  filter(lag == 1, Species %in% goa_main_spp) %>%
  group_by(Species) %>%
  summarise(
    cv_ABC = if_else(mean(ABC, na.rm = TRUE) > 0, sd(ABC, na.rm = TRUE) / mean(ABC, na.rm = TRUE), NA_real_),
    cv_TAC = if_else(mean(TAC, na.rm = TRUE) > 0, sd(TAC, na.rm = TRUE) / mean(TAC, na.rm = TRUE), NA_real_),
    .groups = "drop"
  )

goa_cv_max <- goa_cv_tbl %>% arrange(desc(cv_ABC)) %>% slice(1)
goa_cv_min_tac <- goa_cv_tbl %>% arrange(cv_TAC) %>% slice(1)

dat <- dfOY %>%
  filter(Species %in% mainspp) %>%
  mutate(
    Stock = fct_reorder(Species, -Order),
    Species = fct_reorder(Species, Order),
    lag = if_else(lag == 2, "two", "one")
  ) %>%
  select(Stock, Species, Year, TAC, ABC, lag, Order) %>%
  pivot_wider(names_from = lag, values_from = c(TAC, ABC)) %>%
  mutate(
    ABC_deltat = (ABC_one - ABC_two) / 1000,
    TAC_deltat = (TAC_one - TAC_two) / 1000,
    ABC_delta = if_else(ABC_two > 0, (ABC_one - ABC_two) / ABC_two, NA_real_),
    TAC_delta = if_else(ABC_two > 0, (TAC_one - TAC_two) / ABC_two, NA_real_),
    TAC_sign = if_else(TAC_delta < 0, "-", "+"),
    ABC_sign = if_else(ABC_delta < 0, "-", "+")
  )

mean_deltas <- dat %>%
  summarise(
    mean_abc = mean(ABC_delta, na.rm = TRUE),
    mean_tac = mean(TAC_delta, na.rm = TRUE)
  )

mean_by_species <- dat %>%
  group_by(Species) %>%
  summarise(
    mean_abc = mean(ABC_delta, na.rm = TRUE),
    mean_tac = mean(TAC_delta, na.rm = TRUE),
    .groups = "drop"
  )

max_abc <- mean_by_species %>% arrange(desc(mean_abc)) %>% slice(1)
min_abc <- mean_by_species %>% arrange(mean_abc) %>% slice(1)
min_tac <- mean_by_species %>% arrange(mean_tac) %>% slice(1)

rollover <- dfOY %>%
  filter(Year > 2000, Species %in% mainspp) %>%
  mutate(application = if_else(lag == 1, "Used", "Projected")) %>%
  select(Year, Species, ABC, TAC, application) %>%
  pivot_wider(names_from = application, values_from = c(ABC, TAC)) %>%
  group_by(Species) %>%
  arrange(Year) %>%
  mutate(
    Projected_ABC = lag(ABC_Used),
    Projected_TAC = lag(TAC_Used)
  ) %>%
  ungroup() %>%
  mutate(
    ABC_delta = if_else(Projected_ABC > 0, (ABC_Used - Projected_ABC) / Projected_ABC, NA_real_),
    TAC_delta = if_else(Projected_TAC > 0, (TAC_Used - Projected_TAC) / Projected_TAC, NA_real_)
  ) %>%
  left_join(species_order, by = "Species")

rollover_means <- rollover %>%
  summarise(
    mean_abc = mean(ABC_delta, na.rm = TRUE),
    mean_tac = mean(TAC_delta, na.rm = TRUE)
  )

model_err <- dat %>%
  filter(Year > 2000) %>%
  mutate(
    ABC_err_model = if_else(ABC_one > 0, abs(ABC_one - ABC_two) / ABC_one, NA_real_),
    TAC_err_model = if_else(TAC_one > 0, abs(TAC_one - TAC_two) / TAC_one, NA_real_)
  ) %>%
  select(Species, Year, ABC_err_model, TAC_err_model)

rollover_err <- rollover %>%
  mutate(
    ABC_err_roll = if_else(ABC_Used > 0, abs(ABC_Used - Projected_ABC) / ABC_Used, NA_real_),
    TAC_err_roll = if_else(TAC_Used > 0, abs(TAC_Used - Projected_TAC) / TAC_Used, NA_real_)
  ) %>%
  select(Species, Year, ABC_err_roll, TAC_err_roll)

comp_err <- model_err %>%
  inner_join(rollover_err, by = c("Species", "Year"))

comp_summary <- comp_err %>%
  summarise(
    model_better_abc = mean(ABC_err_model < ABC_err_roll, na.rm = TRUE),
    model_better_tac = mean(TAC_err_model < TAC_err_roll, na.rm = TRUE),
    mean_abc_model = mean(ABC_err_model, na.rm = TRUE),
    mean_abc_roll = mean(ABC_err_roll, na.rm = TRUE),
    mean_tac_model = mean(TAC_err_model, na.rm = TRUE),
    mean_tac_roll = mean(TAC_err_roll, na.rm = TRUE)
  )

comp_by_species <- comp_err %>%
  group_by(Species) %>%
  summarise(
    mean_abc_model = mean(ABC_err_model, na.rm = TRUE),
    mean_abc_roll = mean(ABC_err_roll, na.rm = TRUE),
    mean_tac_model = mean(TAC_err_model, na.rm = TRUE),
    mean_tac_roll = mean(TAC_err_roll, na.rm = TRUE),
    .groups = "drop"
  )

species_abc_adv <- sum(comp_by_species$mean_abc_model < comp_by_species$mean_abc_roll, na.rm = TRUE)
species_tac_adv <- sum(comp_by_species$mean_tac_model < comp_by_species$mean_tac_roll, na.rm = TRUE)

Introduction

Annual ABC and TAC decisions in the Bering Sea and Aleutian Islands (BSAI) groundfish fisheries set the management ceilings that balance conservation, stability, and economic planning. Multi-year advice provides predictable expectations, but it is built on forecasts that can change as new data arrive, which can create abrupt revisions in advice.

This work evaluates whether two-year model projections provide a meaningful improvement over a simple rollover of the previous year’s values. Quantifying when projections track final outcomes, and when they do not, helps interpret the reliability of advice, supports communication with stakeholders, and identifies where interim updates may add value.

This document compares final (lag 1) ABC and TAC values to their two-year projections (lag 2) using OY = 1 records. The seven main species listed above account for about 88% of total ABC, so most of the signal comes from them.

Methods

Data were compiled for the 1986–2025 period for the Bering Sea and Aleutian Islands (BSAI) region, consistent with the management framework described in the BSAI fishery management plan and the BSAI groundfish amendments action summary. The dataset includes ABCs, TACs, and OFLs set by management for the next year and the year after (lag 1 and lag 2). We first evaluated overall variability across the main stocks for both ABC and TAC, then examined patterns in the Pollock stock in more detail.

As an exploratory extension, we attempted to replace the ATTACH catch function with a DSEM-based approach that maps ABC inputs to expected TACs for key species while modeling cross-species dependence and temporal dynamics. This provides a direct comparison to the static SUR-style framework in ATTACH and helps assess whether a dynamic model can better replicate management outcomes.

Performance measures

For each species \(s\) and projection year \(t\), let \(X_{t,1}\) be the final (lag 1) value and \(X_{t,2}\) the two-year projection (lag 2), where \(X\) is either ABC or TAC. We computed percent differences as:

\[ \Delta X_{t} = \frac{X_{t,1} - X_{t,2}}{X_{t,2}}. \]

To keep ABC and TAC deltas on a common scale, TAC percent changes were scaled by the two-year ABC:

\[ \Delta TAC_{t} = \frac{TAC_{t,1} - TAC_{t,2}}{ABC_{t,2}}. \]

Absolute changes in thousand tons are:

\[ \Delta X_{t}^{(kt)} = \frac{X_{t,1} - X_{t,2}}{1000}. \]

Interannual variability was summarized using the coefficient of variation:

\[ CV(X) = \frac{sd(X)}{mean(X)}. \]

For the projection-versus-rollover comparison, we used absolute percent error against the final value:

\[ APE_{model}(X) = \frac{\lvert X_{t,1} - X_{t,2} \rvert}{X_{t,1}}, \quad APE_{roll}(X) = \frac{\lvert X_{t,1} - X_{t-1,1} \rvert}{X_{t,1}}. \]

Dynamic structural equation model setup

This section sets up dynamic structural equation models (DSEMs) using the dsem package. DSEMs blend structural equation modeling with state-space time-series dynamics to estimate both contemporaneous effects and temporal persistence, allowing mechanistic interpretation of multivariate ecological systems (Thorson et al. 2024). Below, the DSEM work is shown separately for BSAI and GOA. In both regions, the dependent variables are TACs for Pollock and Pacific cod, with contemporaneous ABC inputs drawn from a region-specific set of major species.

Seemingly unrelated regression (SUR) provides a useful static reference point for understanding dynamic structural equation models (DSEM). In SUR, multiple regression equations are estimated jointly, allowing for contemporaneous correlation among equation-specific error terms while excluding any dynamic feedback or causal coupling across outcomes (Zellner 1962). DSEM can be viewed as a strict generalization of this framework: correlated disturbances are retained, but are embedded within a dynamic state-space formulation that allows for autoregressive behavior, latent processes, and explicit structural pathways across variables and time (Asparouhov, Hamaker, and Muthén 2018). From this perspective, SUR corresponds to a single-time-slice, purely contemporaneous special case of DSEM in which temporal dependence and structural effects are suppressed. This connection clarifies how joint modeling gains efficiency through shared stochastic components, while DSEM extends that logic by separating persistent process variation from transient noise and propagating information forward through time.

In the BSAI context, the ATTACH model implemented in the catchfunction package uses a two-step regression workflow (ABC → TAC, then TAC → catch) with SUR linking TAC-stage errors across species, and an ensemble of alternative error structures in the catch stage to reflect correlated shocks and the ecosystem cap (Faig and Haynie 2020). This provides a static, cross-equation baseline that motivates the move to DSEM when temporal dynamics and lagged effects are of interest.

For clarity and understanding, we constructed similar log-log TAC~ABC regressions for the main stocks and pooled “Others”. Here, TACs for each stock were estimated statistically using a log-linear model. For \(j = 1, 2, \ldots, n\) stocks, the general model for stock \(i\) took the form:

\[ \ln(\mathrm{TAC}_{i,t}) = \alpha_i + \beta_i \ln(\mathrm{ABC}_{i,t}) + \sum_{j \ne i}^{n} \beta_{ij} \ln(\mathrm{ABC}_{j,t}) + \sum_{k=1}^{m} \beta_k I_{k,t} + \varepsilon_{i,t}. \]

where \(\alpha_i\) is the stock-specific intercept for species \(i\), \(\beta_i\) is the elasticity of the TAC of species \(i\) with respect to its own ABC, and \(\beta_{ij}\) is the elasticity for the ABC of species \(j \neq i\) in the TAC equation for species \(i\). The effect (\(\beta_k\)) of \(k = 1, 2, \ldots, m\) events or policy changes (e.g., changes in management, area closures, or implementation of catch share programs) on TAC was also estimated, where \(I_{k,t}\) is an indicator variable for event \(k\) in year \(t\), and \(\varepsilon_{i,t}\) denotes the residual error for the prediction in year \(t\).

Results

Patterns in catch advice and limits

BSAI totals by year

Here we show simple annual BSAI totals for TAC, ABC, and OFL (lag 1, OY = 1), summed across species.

bsaiyr <- dfOY %>%
  filter(lag == 1) %>%
  group_by(Year) %>%
  summarise(
    TAC = sum(TAC, na.rm = TRUE),
    ABC = sum(ABC, na.rm = TRUE),
    OFL = sum(OFL, na.rm = TRUE),
    n_TAC = sum(!is.na(TAC)),
    n_ABC = sum(!is.na(ABC)),
    n_OFL = sum(!is.na(OFL)),
    .groups = "drop"
  )

plot_bsai <- bsaiyr %>%
  select(Year, TAC, ABC, OFL) %>%
  pivot_longer(-Year, names_to = "Measure", values_to = "value") %>%
  mutate(
    value = na_if(value, 0),
    Measure = factor(Measure, levels = c("TAC", "ABC", "OFL"))
  )

ggplot(plot_bsai, aes(x = Year, y = value / 1000, color = Measure)) +
  geom_line(linewidth = 1) +
  scale_y_continuous(labels = comma) +
  labs(x = NULL, y = "thousand t", color = NULL) +
  theme(legend.position = "top")

Figure 1: BSAI totals by year (sum across species): TAC, ABC, and OFL.

BSAI Pollock

Pollock provides a clear example of how ABC, OFL, and TAC move together over time (see Figure 2). OFL generally sits above ABC, and TAC sits below ABC.

pollock_ts <- dfOY %>%
  filter(lag == 1, Species == "Pollock") %>%
  pivot_longer(cols = c(ABC, OFL, TAC), names_to = "Measure", values_to = "value") %>%
  mutate(value = na_if(value, 0))

ggplot(pollock_ts, aes(x = Year, y = value / 1000, color = Measure)) +
  geom_line(linewidth = 1) +
  scale_y_continuous(labels = comma) +
  scale_color_brewer(palette = "Dark2") +
  labs(x = NULL, y = "thousand t", color = NULL)

Figure 2: BSAI Pollock ABC, OFL, and TAC over time (lag 1, OY = 1).

GOA totals by year

The GOA harvest specifications file was munged into the same long format as the BSAI data. Here we show annual totals from GOA Area == "Total" rows for TAC, ABC, and OFL. Because the GOA spreadsheet contains some n/a fields (especially for OFL/ABC in early years), we report totals as the sum of available numeric entries.

goayr <- goaraw %>%
  filter(lag == 1, OY == 1, Area == "Total") %>%
  group_by(Year = ProjYear) %>%
  summarise(
    TAC = sum(TAC, na.rm = TRUE),
    ABC = sum(ABC, na.rm = TRUE),
    OFL = sum(OFL, na.rm = TRUE),
    n_TAC = sum(!is.na(TAC)),
    n_ABC = sum(!is.na(ABC)),
    n_OFL = sum(!is.na(OFL)),
    .groups = "drop"
  )

plot_goa <- goayr %>%
  select(Year, TAC, ABC, OFL) %>%
  pivot_longer(-Year, names_to = "Measure", values_to = "value") %>%
  mutate(
    value = na_if(value, 0),
    Measure = factor(Measure, levels = c("TAC", "ABC", "OFL"))
  )

ggplot(plot_goa, aes(x = Year, y = value / 1000, color = Measure)) +
  geom_line(linewidth = 1) +
  scale_y_continuous(labels = comma) +
  labs(x = NULL, y = "thousand t", color = NULL) +
  theme(legend.position = "top")

Figure 3: GOA totals by year (sum across all rows): TAC, ABC, and OFL.

GOA Pollock

GOA Pollock follows the same comparison format as Figure 4 (Figure 4), showing annual ABC, OFL, and TAC trajectories.

goa_pollock_ts <- goaraw %>%
  filter(lag == 1, OY == 1, Area == "Total", Species == "Pollock") %>%
  group_by(Year = ProjYear) %>%
  summarise(
    ABC = sum(ABC, na.rm = TRUE),
    OFL = sum(OFL, na.rm = TRUE),
    TAC = sum(TAC, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  pivot_longer(cols = c(ABC, OFL, TAC), names_to = "Measure", values_to = "value") %>%
  mutate(value = na_if(value, 0))

ggplot(goa_pollock_ts, aes(x = Year, y = value / 1000, color = Measure)) +
  geom_line(linewidth = 1) +
  scale_y_continuous(labels = comma) +
  scale_color_brewer(palette = "Dark2") +
  labs(x = NULL, y = "thousand t", color = NULL)

Figure 4: GOA Pollock ABC, OFL, and TAC over time (lag 1, OY = 1, Area = Total).

Variability across species

BSAI

Figure 5 shows the interannual CVs for ABC and TAC. Pacific ocean perch has the highest ABC variability (CV about 0.95), while Pollock has the lowest TAC variability (CV about 0.12), indicating tighter management consistency for the largest stock.

cv_long <- cv_tbl %>%
  pivot_longer(cols = c(cv_ABC, cv_TAC), names_to = "Measure", values_to = "CV") %>%
  mutate(
    Measure = recode(Measure, cv_ABC = "ABC", cv_TAC = "TAC"),
    Species = fct_reorder(Species, CV)
  )

ggplot(cv_long, aes(x = Species, y = CV, fill = Measure)) +
  geom_col(position = position_dodge(width = 0.8)) +
  scale_y_continuous(labels = percent_format(accuracy = 1)) +
  coord_flip() +
  labs(x = NULL, y = "Coefficient of variation", fill = NULL) +
  theme(legend.position = "top")

Figure 5: Interannual variability (CV) of ABC and TAC for the main BSAI species, 1986-2026.

GOA

For GOA, we calculated the same CV metrics for the seven largest species by mean ABC (lag 1, OY = 1, Area == "Total"). Figure 6 shows that Other Flounder has the highest ABC variability (CV about 1.18), while Sablefish has the lowest TAC variability (CV about 0.33).

goa_cv_long <- goa_cv_tbl %>%
  pivot_longer(cols = c(cv_ABC, cv_TAC), names_to = "Measure", values_to = "CV") %>%
  mutate(
    Measure = recode(Measure, cv_ABC = "ABC", cv_TAC = "TAC"),
    Species = fct_reorder(Species, CV)
  )

ggplot(goa_cv_long, aes(x = Species, y = CV, fill = Measure)) +
  geom_col(position = position_dodge(width = 0.8)) +
  scale_y_continuous(labels = percent_format(accuracy = 1)) +
  coord_flip() +
  labs(x = NULL, y = "Coefficient of variation", fill = NULL) +
  theme(legend.position = "top")

Figure 6: Interannual variability (CV) of ABC and TAC for selected GOA species, 1986-2024 (top seven by mean ABC; lag 1, OY = 1, Area = Total).

Two-year vs final differences (percent)

Figure 7 summarizes the percent difference between final and two-year values, scaled by the two-year ABC to keep ABC and TAC deltas on the same scale. Across species, the average ABC adjustment is about 4.8%, while TAC adjustments are smaller at 2.4%. Atka mackerel has the largest positive mean ABC revision (13.0%), while Northern rock sole shows the only negative mean ABC adjustment (-1.7%); Pollock has the most negative mean TAC change (-3.7%).

delta_long_pct <- dat %>%
  filter(Year > 1990) %>%
  select(Year, Species, Order, ABC_delta, TAC_delta) %>%
  pivot_longer(cols = c(ABC_delta, TAC_delta), names_to = "Measure", values_to = "Delta") %>%
  mutate(
    Measure = recode(Measure, ABC_delta = "ABC", TAC_delta = "TAC"),
    Species = fct_reorder(Species, Order)
  )

ggplot(delta_long_pct, aes(x = Year, y = Delta, fill = Delta)) +
  geom_col() +
  facet_grid(Species ~ Measure, scales = "free_y") +
  scale_y_continuous(labels = percent_format(accuracy = 1)) +
  scale_fill_gradient2(
    low = muted("red"),
    mid = "white",
    high = muted("blue"),
    midpoint = 0,
    na.value = "grey90"
  ) +
  labs(x = NULL, y = "Percent change", fill = "Change") +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    strip.text.y = element_text(angle = 0)
  )

Figure 7: Percent changes between final (lag 1) and two-year (lag 2) values, scaled by the two-year ABC.

Two-year vs final differences (kt)

The heatmaps in Figure 8 highlight the absolute magnitude of adjustments in thousands of tons. Large swings cluster in the biggest stocks (notably Pollock), while smaller species show less absolute movement even when percent changes are sizable.

complete_years <- dat %>%
  group_by(Year) %>%
  summarise(miss = sum(is.na(ABC_deltat) | is.na(TAC_deltat)), .groups = "drop") %>%
  filter(miss == 0) %>%
  pull(Year)

delta_long_kt <- dat %>%
  filter(Year > 1990, Year %in% complete_years) %>%
  select(Year, Species, Order, ABC_deltat, TAC_deltat) %>%
  pivot_longer(cols = c(ABC_deltat, TAC_deltat), names_to = "Measure", values_to = "Delta") %>%
  mutate(
    Measure = recode(Measure, ABC_deltat = "ABC", TAC_deltat = "TAC"),
    Species = fct_reorder(Species, Order)
  )

ggplot(delta_long_kt, aes(x = Year, y = Species, fill = Delta)) +
  geom_tile() +
  facet_grid(Measure ~ .) +
  scale_fill_gradient2(
    low = muted("red"),
    mid = "white",
    high = muted("blue"),
    midpoint = 0,
    na.value = "grey90",
    labels = comma
  ) +
  labs(x = NULL, y = NULL, fill = "Change (kt)") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Figure 8: Absolute changes between final and two-year values (thousand tons).

The results shown in Figure 9 is provided to show contrast among other stocks with pollock removed.

delta_long_kt_nopollock <- delta_long_kt %>%
  filter(Species != "Pollock")

ggplot(delta_long_kt_nopollock, aes(x = Year, y = Species, fill = Delta)) +
  geom_tile() +
  facet_grid(Measure ~ .) +
  scale_fill_gradient2(
    low = muted("red"),
    mid = "white",
    high = muted("blue"),
    midpoint = 0,
    na.value = "grey90",
    labels = comma
  ) +
  labs(x = NULL, y = NULL, fill = "Change (kt)") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Figure 9: Absolute changes between final and two-year values (thousand tons), excluding Pollock.

Rollover scenario

To gauge a simple alternative, Figure 10 compares final values to a naive rollover that uses the previous year’s final value as the “projection.” The average rollover changes are small but positive (2.7% for ABC and 2.0% for TAC), suggesting a mild upward drift when viewed against a strict carry-forward baseline.

rollover_long <- rollover %>%
  select(Year, Species, Order, ABC_delta, TAC_delta) %>%
  pivot_longer(cols = c(ABC_delta, TAC_delta), names_to = "Measure", values_to = "Delta") %>%
  mutate(
    Measure = recode(Measure, ABC_delta = "ABC", TAC_delta = "TAC"),
    Species = fct_reorder(Species, Order)
  )

ggplot(rollover_long, aes(x = Year, y = Delta, fill = Delta)) +
  geom_col() +
  facet_grid(Species ~ Measure, scales = "free_y") +
  scale_y_continuous(labels = percent_format(accuracy = 1)) +
  scale_fill_gradient2(
    low = muted("red"),
    mid = "white",
    high = muted("blue"),
    midpoint = 0,
    na.value = "grey90"
  ) +
  labs(x = NULL, y = "Rollover change", fill = "Change") +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    strip.text.y = element_text(angle = 0)
  ) + ggthemes::theme_few()

Figure 10: Rollover comparison: final values vs prior-year final values used as a naive projection.

Projection vs rollover accuracy

Using the overlap period (2001 onward), we compared how close the two-year model projections (lag 2) and a naive rollover (previous year’s lag 1 value) came to the final values used (lag 1). On average, the model-based projections show slightly lower absolute percent error for ABC (13.1% vs 13.8%), and a similarly small improvement for TAC (11.5% vs 12.0%). However, the model is closer than the rollover in only 52% of ABC cases and 35% of TAC cases, and the mean advantage appears in 3 of 7 stocks for ABC and 4 of 7 stocks for TAC. This suggests there is enough information to compare the two approaches, but the performance edge of two-year projections over a rollover baseline is modest and uneven—more evident for ABC than for TAC. Table 1 summarizes the mean absolute percent errors by species for each approach.

tbl_proj_roll <- comp_by_species %>%
  mutate(
    abc_diff = mean_abc_model - mean_abc_roll,
    tac_diff = mean_tac_model - mean_tac_roll
  ) %>%
  select(
    Species,
    mean_abc_model,
    mean_abc_roll,
    abc_diff,
    mean_tac_model,
    mean_tac_roll,
    tac_diff
  ) %>%
  mutate(across(-Species, ~ percent(.x, accuracy = 0.1)))

knitr::kable(
  tbl_proj_roll,
  col.names = c(
    "Species",
    "ABC model",
    "ABC rollover",
    "ABC model - rollover",
    "TAC model",
    "TAC rollover",
    "TAC model - rollover"
  ),
  align = "lrrrrrr"
)

Table 1: Mean absolute percent error for two-year projections vs rollover baselines (2001+).

Species	ABC model	ABC rollover	ABC model - rollover	TAC model	TAC rollover	TAC model - rollover
Atka mackerel	16.1%	19.8%	-3.7%	15.5%	15.0%	0.6%
Flathead sole	5.5%	5.5%	0.0%	17.3%	17.7%	-0.4%
Northern rock sole	17.7%	18.4%	-0.7%	10.7%	10.7%	0.0%
Pacific cod	10.4%	9.6%	0.9%	9.1%	8.6%	0.5%
Pacific ocean perch	9.7%	13.5%	-3.8%	7.6%	12.5%	-4.9%
Pollock	18.0%	17.0%	1.0%	8.4%	7.1%	1.3%
Yellowfin sole	14.1%	12.7%	1.4%	11.8%	12.4%	-0.6%

Projected vs used time series

Figure 11 contrasts the two-year projections with the final values used in the next year. The projected series are generally smoother, with the largest departures in higher-variability stocks (Pacific ocean perch and Flathead sole), which aligns with the elevated CVs in Figure 5.

ts_dat <- dfOY %>%
  filter(Year > 2004, Species %in% mainspp) %>%
  mutate(Series = if_else(lag == 1, "Used (lag 1)", "Projected (lag 2)")) %>%
  select(Year, Species, ABC, Series) %>%
  left_join(species_order, by = "Species") %>%
  mutate(Species = fct_reorder(Species, Order))

ggplot(ts_dat, aes(x = Year, y = ABC / 1000, color = Series)) +
  geom_line(linewidth = 0.9) +
  facet_wrap(~Species, scales = "free_y") +
  scale_y_continuous(labels = comma) +
  scale_color_brewer(palette = "Dark2") +
  labs(x = NULL, y = "ABC (thousand t)", color = NULL)

Figure 11: Projected (lag 2) vs used (lag 1) ABC time series by species.

TAC cross-species relationships

BSAI

Figure 12 shows pairwise relationships among BSAI TAC values for the main species, with all remaining species pooled as “Others,” using points with a smooth fit in the lower triangle.

tac_pairs_spp <- setdiff(mainspp, c("Pacific ocean perch", "Flathead sole"))

tac_main <- dfOY %>%
  filter(lag == 1, Species %in% tac_pairs_spp) %>%
  mutate(TAC = TAC / 1000) %>%
  select(Year, Species, TAC) %>%
  pivot_wider(names_from = Species, values_from = TAC)

tac_other <- dfOY %>%
  filter(lag == 1, !Species %in% tac_pairs_spp) %>%
  group_by(Year) %>%
  summarise(Others = sum(TAC, na.rm = TRUE) / 1000, .groups = "drop")

tac_wide <- tac_main %>%
  left_join(tac_other, by = "Year") %>%
  select(Year, all_of(tac_pairs_spp), Others)

spp_labels <- c(tac_pairs_spp, "Others")
spp_cols <- gsub(" ", "_", spp_labels)
tac_wide <- tac_wide %>% rename_with(~ spp_cols, all_of(spp_labels))

# Draw lower-triangle pairs panels with points plus a LOESS trend.
lower_points_smooth <- function(data, mapping, ...) {
  ggplot(data = data, mapping = mapping) +
    geom_point(alpha = 0.5, size = 0.7) +
    geom_smooth(method = "loess", se = FALSE, linewidth = 0.6, color = "steelblue")
}

GGally::ggpairs(
  tac_wide,
  columns = spp_cols,
  columnLabels = spp_labels,
  lower = list(continuous = lower_points_smooth),
  diag = list(continuous = "densityDiag"),
  upper = list(continuous = "cor")
) +
  ggthemes::theme_few(base_size = 9)

Figure 12: BSAI pairs plot of TAC across the main species plus Others (lag 1, OY = 1).

GOA

Figure 13 shows the analogous GOA TAC relationships using the top five species by total TAC (lag 1, OY = 1, Area == "Total"), with all remaining GOA species pooled as “Others.”

goa_tac_pairs_spp <- goaraw %>%
  filter(lag == 1, OY == 1, Area == "Total") %>%
  group_by(Species) %>%
  summarise(total_TAC = sum(TAC, na.rm = TRUE), .groups = "drop") %>%
  filter(is.finite(total_TAC), total_TAC > 0) %>%
  arrange(desc(total_TAC)) %>%
  slice_head(n = 5) %>%
  pull(Species)

goa_tac_main <- goaraw %>%
  filter(lag == 1, OY == 1, Area == "Total", Species %in% goa_tac_pairs_spp) %>%
  group_by(Year = ProjYear, Species) %>%
  summarise(TAC = sum(TAC, na.rm = TRUE) / 1000, .groups = "drop") %>%
  pivot_wider(names_from = Species, values_from = TAC)

goa_tac_other <- goaraw %>%
  filter(lag == 1, OY == 1, Area == "Total", !Species %in% goa_tac_pairs_spp) %>%
  group_by(Year = ProjYear) %>%
  summarise(Others = sum(TAC, na.rm = TRUE) / 1000, .groups = "drop")

goa_tac_wide <- goa_tac_main %>%
  left_join(goa_tac_other, by = "Year") %>%
  select(Year, all_of(goa_tac_pairs_spp), Others)

goa_spp_labels <- c(goa_tac_pairs_spp, "Others")
goa_spp_cols <- gsub("[^A-Za-z0-9]+", "_", goa_spp_labels)
goa_tac_wide <- goa_tac_wide %>% rename_with(~ goa_spp_cols, all_of(goa_spp_labels))

GGally::ggpairs(
  goa_tac_wide,
  columns = goa_spp_cols,
  columnLabels = goa_spp_labels,
  lower = list(continuous = lower_points_smooth),
  diag = list(continuous = "densityDiag"),
  upper = list(continuous = "cor")
) +
  ggthemes::theme_few(base_size = 9)

Figure 13: GOA pairs plot of TAC across top five species (by total TAC) plus Others (lag 1, OY = 1, Area = Total).

TAC-ABC regressions

As noted, we fit multiple regressions with log(ABC) for the focal stock, separate log(ABC) terms for each other stock group, and a centered year term to proxy policy indicators.

# Convert arbitrary stock labels into safe column-name fragments.
safe_names <- function(x) {
  x %>%
    gsub("[^A-Za-z0-9]+", "_", ., perl = TRUE) %>%
    gsub("^_|_$", "", ., perl = TRUE)
}

# Fit stock-specific log(TAC) regressions using own-stock ABC, other-stock ABC,
# and a centered year effect; return coefficient summaries and fitted values.
fit_tac_abc_regressions <- function(data, main_groups) {
  group_names <- c(main_groups, "Other")

  reg_base <- data %>%
    mutate(Group = if_else(Species %in% main_groups, Species, "Other")) %>%
    group_by(Year, Group) %>%
    summarise(
      ABC = sum(ABC, na.rm = TRUE),
      TAC = sum(TAC, na.rm = TRUE),
      .groups = "drop"
    ) %>%
    pivot_wider(names_from = Group, values_from = c(ABC, TAC), values_fill = 0)

  names(reg_base) <- safe_names(names(reg_base))
  group_safe <- safe_names(group_names)

  reg_fits <- purrr::map2(group_names, group_safe, function(group, group_safe_i) {
    abc_col <- paste0("ABC_", group_safe_i)
    tac_col <- paste0("TAC_", group_safe_i)
    other_groups <- setdiff(group_names, group)
    other_cols <- paste0("ABC_", safe_names(other_groups))

    df <- reg_base %>%
      transmute(
        Year,
        ABC_i = .data[[abc_col]],
        TAC_i = .data[[tac_col]],
        across(all_of(other_cols), ~ .x, .names = "{.col}"),
        Year_c = Year - mean(Year)
      ) %>%
      filter(TAC_i > 0) %>%
      filter(if_all(all_of(c("ABC_i", other_cols)), ~ .x > 0)) %>%
      mutate(
        log_TAC_i = log(TAC_i),
        log_ABC_i = log(ABC_i),
        across(all_of(other_cols), ~ log(.x), .names = "log_{.col}")
      )

    rhs <- c("log_ABC_i", paste0("log_", other_cols), "Year_c")
    model <- lm(as.formula(paste("log_TAC_i ~", paste(rhs, collapse = " + "))), data = df)
    list(species = group, data = df, model = model)
  })

  reg_results <- purrr::map_dfr(reg_fits, function(x) {
    coefs <- coef(x$model)
    out <- tibble(
      Species = x$species,
      n = nrow(x$data),
      Intercept = unname(coefs["(Intercept)"]),
      log_ABC = unname(coefs["log_ABC_i"]),
      Year_c = unname(coefs["Year_c"]),
      R2 = summary(x$model)$r.squared
    )

    for (g in group_names) {
      col_name <- paste0("log_ABC_", safe_names(g))
      out[[col_name]] <- if (g == x$species) NA_real_ else unname(coefs[col_name])
    }

    out
  }) %>%
    mutate(order = match(Species, group_names)) %>%
    arrange(order) %>%
    select(-order)

  reg_pred <- purrr::map_dfr(reg_fits, function(x) {
    pred_log <- predict(x$model, newdata = x$data)
    tibble(
      Year = x$data$Year,
      Species = x$species,
      Observed = x$data$TAC_i,
      Predicted = exp(pred_log)
    )
  })

  list(
    results = reg_results,
    pred = reg_pred,
    group_names = group_names
  )
}

# Format the regression coefficient summary as a gt table.
render_regression_table <- function(reg_results, group_names) {
  abc_cols <- paste0("log_ABC_", safe_names(group_names))

  reg_results %>%
    gt::gt(rowname_col = "Species") %>%
  gt::cols_label(
    n = "n",
    Intercept = "Intercept",
    log_ABC = "log(ABC)",
    Year_c = "Year (centered)",
    R2 = "R2",
    !!!setNames(paste0("log(ABC ", group_names, ")"), abc_cols)
  ) %>%
  gt::fmt_number(columns = c(Intercept, log_ABC, all_of(abc_cols), Year_c, R2), decimals = 3) %>%
  gt::fmt_number(columns = n, decimals = 0)
}

# Plot observed TAC against fitted regression predictions by stock.
plot_regression_predictions <- function(reg_pred) {
  ggplot(reg_pred, aes(x = Year)) +
    geom_line(aes(y = Predicted, color = "Predicted"), linewidth = 1) +
    geom_point(aes(y = Observed, color = "Observed"), size = 1.2, alpha = 0.8) +
    facet_grid(Species ~ ., scales = "free_y") +
    scale_color_manual(values = c(Predicted = "steelblue", Observed = "grey40")) +
    scale_y_continuous(labels = comma) +
    labs(x = NULL, y = "TAC", color = NULL)
}

bsai_reg_main5 <- setdiff(mainspp, c("Pacific ocean perch", "Flathead sole"))
bsai_reg <- fit_tac_abc_regressions(
  data = dfOY %>% filter(lag == 1),
  main_groups = bsai_reg_main5
)

goa_reg_main5 <- goaraw %>%
  filter(lag == 1, OY == 1, Area == "Total") %>%
  group_by(Species) %>%
  summarise(total_TAC = sum(TAC, na.rm = TRUE), .groups = "drop") %>%
  filter(is.finite(total_TAC), total_TAC > 0) %>%
  arrange(desc(total_TAC)) %>%
  slice_head(n = 5) %>%
  pull(Species)

goa_reg <- fit_tac_abc_regressions(
  data = goaOY %>% filter(lag == 1),
  main_groups = goa_reg_main5
)

BSAI

For BSAI, we retain the existing five-stock specification: Pollock, Pacific cod, Yellowfin sole, Atka mackerel, Northern rock sole, plus a pooled “Other” category.

render_regression_table(bsai_reg$results, bsai_reg$group_names)

Table 2: BSAI log-log TAC~ABC regressions with separate other-stock terms (lag 1, OY = 1).

	n	Intercept	log(ABC)	Year (centered)	R2	log(ABC Pollock)	log(ABC Yellowfin sole)	log(ABC Pacific cod)	log(ABC Atka mackerel)	log(ABC Northern rock sole)	log(ABC Other)
Pollock	38	10.218	0.364	−0.005	0.782	NA	−0.016	0.034	0.025	−0.067	−0.079
Yellowfin sole	38	11.507	0.758	0.002	0.866	−0.615	NA	0.166	0.014	−0.132	−0.050
Pacific cod	38	2.805	0.700	0.000	0.862	0.036	−0.164	NA	0.071	0.151	−0.027
Atka mackerel	38	−7.435	0.543	0.008	0.709	0.375	−0.008	−0.363	NA	0.562	0.356
Northern rock sole	38	34.044	0.366	0.013	0.607	−0.772	−0.343	0.028	−0.207	NA	−0.787
Other	38	33.128	−0.175	0.021	0.728	−0.782	−0.202	−0.527	−0.203	0.286	NA

plot_regression_predictions(bsai_reg$pred)

Figure 14: BSAI observed vs predicted TAC from log-linear multiple regressions (lag 1, OY = 1).

GOA

For GOA, we use the top five species by total TAC from the Area == "Total" series (Pollock, Arrowtooth Flounder, Pacific cod, Shallow-water Flatfish, Sablefish), plus a pooled “Other” category.

render_regression_table(goa_reg$results, goa_reg$group_names)

Table 3: GOA log-log TAC~ABC regressions with separate other-stock terms (lag 1, OY = 1, Area = Total).

	n	Intercept	log(ABC)	Year (centered)	R2	log(ABC Pollock)	log(ABC Arrowtooth Flounder)	log(ABC Pacific cod)	log(ABC Shallow-water Flatfish)	log(ABC Sablefish)	log(ABC Other)
Pollock	35	1.131	0.926	0.003	0.987	NA	−0.082	0.051	0.016	−0.004	0.000
Arrowtooth Flounder	35	9.902	−0.454	0.029	0.899	0.221	NA	0.120	−0.302	−0.169	0.639
Pacific cod	35	−3.621	1.052	−0.011	0.982	0.016	0.048	NA	−0.109	0.122	0.178
Shallow-water Flatfish	35	14.347	−0.383	0.037	0.937	0.112	0.069	−0.124	NA	0.018	−0.097
Sablefish	35	2.541	0.691	−0.011	0.928	−0.015	−0.069	−0.046	−0.022	NA	0.188
Other	35	8.114	0.297	0.013	0.958	0.178	−0.086	−0.036	−0.102	0.026	NA

plot_regression_predictions(goa_reg$pred)

Figure 15: GOA observed vs predicted TAC from log-linear multiple regressions (lag 1, OY = 1, Area = Total).

GOA OY-cap rescaling alternative

As a GOA-specific alternative, we considered the weighted OY-cap rescaling method described in the separate GOA OY-cap note. The method starts with projected stock-specific allocations \(C_i\) and, when the total exceeds a GOA-wide cap \(OY\), rescales them to

\[ C'_i = C_i \times r^{(1 / (w_i m))}, \qquad r = \frac{OY}{\sum_i C_i}, \]

where \(w_i > 0\) are stock-specific weights and \(m\) is a common multiplier chosen so that the adjusted allocations sum to the cap. Higher weights imply less reduction, so the method can preserve relatively prioritized stocks while still satisfying the aggregate constraint.

For illustration, we applied the method to the latest GOA Area == "Total" ABC vector and evaluated three lower cap levels: 95%, 90%, and 85% of the observed aggregate ABC. Because this repository does not include ex-vessel value series or explicit policy weights, we use each stock’s historical mean TAC/ABC ratio as a simple proxy weight. That makes this a sensitivity analysis rather than a recommended operational rule, but it is still useful for seeing how an explicit cap-rescaling mechanism differs from the regression and DSEM approaches above.

# Rescale projected catches to satisfy an aggregate OY cap while preserving
# relative priority through stock-specific weights.
rescale_oy_cap <- function(catch, weights, oy_cap, tol = 1e-8, max_iter = 200) {
  catch <- as.numeric(catch)
  weights <- as.numeric(weights)

  total_catch <- sum(catch, na.rm = TRUE)
  if (!is.finite(total_catch) || total_catch <= 0) {
    stop("Total projected catch must be positive.")
  }

  if (oy_cap >= total_catch) {
    return(list(multiplier = NA_real_, scaled = catch, binding = FALSE))
  }

  r <- oy_cap / total_catch
  # Measure the current difference between the rescaled total and the cap.
  total_gap <- function(m) {
    sum(catch * r^(1 / (weights * m))) - oy_cap
  }

  lower <- 1e-6
  upper <- 2
  while (total_gap(upper) < 0 && upper < 1e6) {
    upper <- upper * 2
  }

  for (i in seq_len(max_iter)) {
    mid <- (lower + upper) / 2
    gap_mid <- total_gap(mid)
    if (abs(gap_mid) < tol) break
    if (gap_mid > 0) {
      upper <- mid
    } else {
      lower <- mid
    }
  }

  multiplier <- (lower + upper) / 2
  scaled <- catch * r^(1 / (weights * multiplier))

  list(multiplier = multiplier, scaled = scaled, binding = TRUE)
}

goa_oy_latest_year <- goaOY %>%
  filter(lag == 1) %>%
  summarise(Year = max(Year, na.rm = TRUE)) %>%
  pull(Year)

goa_oy_weights <- goaOY %>%
  filter(lag == 1, is.finite(ABC), ABC > 0, is.finite(TAC), TAC >= 0) %>%
  mutate(weight = TAC / ABC) %>%
  filter(is.finite(weight), weight > 0) %>%
  group_by(Species) %>%
  summarise(weight = mean(weight, na.rm = TRUE), .groups = "drop")

goa_oy_base <- goaOY %>%
  filter(lag == 1, Year == goa_oy_latest_year, is.finite(ABC), ABC > 0) %>%
  select(Species, ABC, TAC) %>%
  left_join(goa_oy_weights, by = "Species") %>%
  mutate(weight = if_else(is.finite(weight) & weight > 0, weight, 0.01)) %>%
  arrange(desc(ABC))

goa_oy_total_abc <- sum(goa_oy_base$ABC, na.rm = TRUE)
goa_oy_total_tac <- sum(goa_oy_base$TAC, na.rm = TRUE)

goa_oy_cap_levels <- tibble(
  Scenario = c("Observed ABC (100%)", "95% cap", "90% cap", "85% cap"),
  cap_factor = c(1.00, 0.95, 0.90, 0.85),
  OY_cap = goa_oy_total_abc * cap_factor
)

goa_oy_alloc <- purrr::map_dfr(seq_len(nrow(goa_oy_cap_levels)), function(i) {
  scenario_row <- goa_oy_cap_levels[i, ]
  fit <- rescale_oy_cap(
    catch = goa_oy_base$ABC,
    weights = goa_oy_base$weight,
    oy_cap = scenario_row$OY_cap
  )

  goa_oy_base %>%
    mutate(
      Scenario = scenario_row$Scenario,
      cap_factor = scenario_row$cap_factor,
      OY_cap = scenario_row$OY_cap,
      multiplier = fit$multiplier,
      ABC_rescaled = fit$scaled,
      pct_reduction = 1 - ABC_rescaled / ABC
    )
})

goa_oy_summary <- goa_oy_alloc %>%
  group_by(Scenario, cap_factor, OY_cap) %>%
  summarise(
    multiplier = first(multiplier),
    total_ABC = sum(ABC, na.rm = TRUE),
    total_rescaled = sum(ABC_rescaled, na.rm = TRUE),
    mean_reduction = mean(pct_reduction, na.rm = TRUE),
    pollock_reduction = pct_reduction[Species == "Pollock"][1],
    pacific_cod_reduction = pct_reduction[Species == "Pacific cod"][1],
    arrowtooth_reduction = pct_reduction[Species == "Arrowtooth Flounder"][1],
    .groups = "drop"
  )

goa_oy_top10 <- goa_oy_base %>%
  slice_head(n = 10) %>%
  pull(Species)

The baseline GOA ABC total in 2024 is 599,784 t, compared with aggregate TAC of 518,320 t. Under the weighting rule used here, the cap-rescaling method reduces lower-priority stocks more than high-utilization stocks such as Pollock, Pacific ocean perch, and Rex sole, while still preserving strictly positive allocations for all stocks with positive baseline ABC.

goa_oy_summary %>%
  mutate(
    OY_cap = comma(round(OY_cap)),
    multiplier = if_else(is.na(multiplier), "NA", format(round(multiplier, 3), nsmall = 3)),
    total_ABC = comma(round(total_ABC)),
    total_rescaled = comma(round(total_rescaled)),
    mean_reduction = percent(mean_reduction, accuracy = 0.1),
    pollock_reduction = percent(pollock_reduction, accuracy = 0.1),
    pacific_cod_reduction = percent(pacific_cod_reduction, accuracy = 0.1),
    arrowtooth_reduction = percent(arrowtooth_reduction, accuracy = 0.1)
  ) %>%
  select(
    Scenario,
    OY_cap,
    multiplier,
    total_ABC,
    total_rescaled,
    mean_reduction,
    pollock_reduction,
    pacific_cod_reduction,
    arrowtooth_reduction
  ) %>%
  gt::gt() %>%
  gt::cols_label(
    OY_cap = "OY cap",
    multiplier = "m",
    total_ABC = "Baseline ABC",
    total_rescaled = "Rescaled ABC",
    mean_reduction = "Mean reduction",
    pollock_reduction = "Pollock reduction",
    pacific_cod_reduction = "Pacific cod reduction",
    arrowtooth_reduction = "Arrowtooth reduction"
  )

Table 4: Illustrative GOA OY-cap rescaling scenarios for the latest annual ABC vector.

Scenario	OY cap	m	Baseline ABC	Rescaled ABC	Mean reduction	Pollock reduction	Pacific cod reduction	Arrowtooth reduction
85% cap	509,816	1.540	599,784	509,816	12.4%	10.2%	12.3%	26.1%
90% cap	539,806	1.550	599,784	539,806	8.2%	6.7%	8.1%	17.7%
95% cap	569,795	1.559	599,784	569,795	4.1%	3.3%	4.0%	9.0%
Observed ABC (100%)	599,784	NA	599,784	599,784	0.0%	0.0%	0.0%	0.0%

goa_oy_plot <- goa_oy_alloc %>%
  filter(cap_factor < 1, Species %in% goa_oy_top10) %>%
  mutate(
    Scenario = factor(Scenario, levels = c("95% cap", "90% cap", "85% cap")),
    Species = factor(Species, levels = rev(goa_oy_top10))
  )

ggplot(goa_oy_plot, aes(x = Species, y = pct_reduction, fill = Scenario)) +
  geom_col(position = position_dodge(width = 0.8)) +
  coord_flip() +
  scale_y_continuous(labels = percent_format(accuracy = 1)) +
  scale_fill_brewer(palette = "Set2") +
  labs(x = NULL, y = "Reduction from baseline ABC", fill = NULL) +
  theme(legend.position = "top")

Figure 16: Illustrative GOA OY-cap rescaling for the ten largest stocks in the latest year, under alternative aggregate cap levels.

DAGs for the DSEM

library(dsem)

# Convert arbitrary stock labels into safe column-name fragments.
safe_names <- function(x) {
  x %>%
    gsub("[^A-Za-z0-9]+", "_", ., perl = TRUE) %>%
    gsub("^_|_$", "", ., perl = TRUE)
}

# Return vertically centered node positions for DAG plots.
centered_positions <- function(n) {
  if (n <= 1) {
    0
  } else {
    seq(from = (n - 1) / 2, to = -(n - 1) / 2, length.out = n)
  }
}

# Build a static or dynamic DAG showing ABC-to-TAC links for the selected stocks.
plot_dsem_dag <- function(dep_species, indep_species, dynamic = FALSE) {
  dep_nodes <- paste0("TAC_", safe_names(dep_species))
  indep_nodes <- paste0("ABC_", safe_names(indep_species))

  static_formulas <- lapply(dep_nodes, function(dep_node) {
    reformulate(indep_nodes, response = dep_node)
  })

  x_coords <- setNames(rep(0, length(indep_nodes)), indep_nodes)
  y_coords <- setNames(centered_positions(length(indep_nodes)), indep_nodes)

  if (dynamic) {
    lag_nodes <- paste0(dep_nodes, "_lag")
    dynamic_formulas <- Map(function(dep_node, lag_node) {
      reformulate(c(indep_nodes, lag_node), response = dep_node)
    }, dep_nodes, lag_nodes)

    labels <- c(
      setNames(paste0("ABC ", indep_species, " (t)"), indep_nodes),
      setNames(paste0("TAC ", dep_species, " (t-1)"), lag_nodes),
      setNames(paste0("TAC ", dep_species, " (t)"), dep_nodes)
    )

    x_coords <- c(
      x_coords,
      setNames(rep(1, length(lag_nodes)), lag_nodes),
      setNames(rep(2, length(dep_nodes)), dep_nodes)
    )
    dep_y <- centered_positions(length(dep_nodes))
    y_coords <- c(
      y_coords,
      setNames(dep_y, lag_nodes),
      setNames(dep_y, dep_nodes)
    )

    dag_obj <- do.call(
      ggdag::dagify,
      c(dynamic_formulas, list(labels = labels, coords = list(x = x_coords, y = y_coords)))
    )
  } else {
    labels <- c(
      setNames(paste("ABC:", indep_species), indep_nodes),
      setNames(paste("TAC:", dep_species), dep_nodes)
    )

    x_coords <- c(x_coords, setNames(rep(2, length(dep_nodes)), dep_nodes))
    y_coords <- c(y_coords, setNames(centered_positions(length(dep_nodes)), dep_nodes))

    dag_obj <- do.call(
      ggdag::dagify,
      c(static_formulas, list(labels = labels, coords = list(x = x_coords, y = y_coords)))
    )
  }

  dag_tidy <- ggdag::tidy_dagitty(dag_obj)

  ggplot(dag_tidy, aes(x = x, y = y, xend = xend, yend = yend)) +
    ggdag::geom_dag_edges() +
    ggdag::geom_dag_node() +
    ggdag::geom_dag_label_repel(aes(label = label), color = "black", size = 3) +
    ggdag::theme_dag()
}

# dsem has changed allowed values for gmrf_parameterization across versions.
# Choose a value that works for the installed package (GitHub Actions may differ
# from local).
# Find a `gmrf_parameterization` value accepted by the installed dsem version.
choose_gmrf_param <- function() {
  candidates <- c(
    "full", "project", "gmrf_project", "mvn_project", # seen in newer versions
    "separable", "projection"                           # seen in older versions
  )
  for (x in candidates) {
    ok <- tryCatch({
      dsem::dsem_control(gmrf_parameterization = x, quiet = TRUE)
      TRUE
    }, error = function(e) FALSE)
    if (ok) return(x)
  }
  stop("Could not find a supported gmrf_parameterization value for dsem")
}

# Apply stored centering and scaling constants to new data.
scale_with <- function(x, center, scale) {
  sweep(sweep(x, 2, center, "-"), 2, scale, "/")
}

# Return the last finite value in a numeric vector.
last_finite <- function(x) {
  x <- x[is.finite(x)]
  if (length(x) == 0) NA_real_ else tail(x, 1)
}

# Fit a regional DSEM, compute leave-one-out diagnostics, and generate an
# illustrative forward TAC forecast from observed or supplied ABC scenarios.
fit_dsem_region <- function(
  data,
  dep_species,
  indep_species,
  n_future = 5,
  seed = 1,
  future_abc_raw = NULL,
  future_years = NULL,
  compute_loo = TRUE
) {
  dsem_df <- data %>%
    filter(Species %in% unique(c(dep_species, indep_species))) %>%
    group_by(Year, Species) %>%
    summarise(
      TAC = sum(TAC, na.rm = TRUE),
      ABC = sum(ABC, na.rm = TRUE),
      .groups = "drop"
    ) %>%
    pivot_wider(
      names_from = Species,
      values_from = c(TAC, ABC),
      names_sep = "_"
    ) %>%
    arrange(Year)

  names(dsem_df) <- safe_names(names(dsem_df))

  dep_safe <- safe_names(dep_species)
  indep_safe <- safe_names(indep_species)
  y_cols <- paste0("TAC_", dep_safe)
  x_cols <- paste0("ABC_", indep_safe)

  ts_raw <- dsem_df[, c(y_cols, x_cols)]
  ts_scaled <- scale(ts_raw)
  ts_center <- attr(ts_scaled, "scaled:center")
  ts_scale <- attr(ts_scaled, "scaled:scale")
  ts_center[is.na(ts_center)] <- 0
  ts_scale[is.na(ts_scale) | ts_scale == 0] <- 1

  tsdata <- ts(ts_scaled, start = min(dsem_df$Year))

  ar_lines <- paste0(y_cols, " -> ", y_cols, ", 1, ar_", dep_safe)
  abc_lines <- unlist(lapply(y_cols, function(y) {
    paste0(x_cols, " -> ", y, ", 0, b_", gsub("ABC_", "", x_cols), "_to_", gsub("TAC_", "", y))
  }))
  sem <- paste(c(ar_lines, abc_lines), collapse = "\n")

  gmrf_param <- choose_gmrf_param()

  # Try one specific DSEM configuration so callers can loop over fallbacks.
  fit_dsem_once <- function(gmrf, family_opt = NULL) {
    args <- list(
      sem = sem,
      tsdata = tsdata,
      estimate_delta0 = TRUE,
      control = dsem_control(
        quiet = TRUE,
        gmrf_parameterization = gmrf,
        newton_loops = 0,
        extra_convergence_checks = FALSE,
        getsd = TRUE
      )
    )
    if (!is.null(family_opt)) args$family <- family_opt
    do.call(dsem::dsem, args)
  }

  gmrf_candidates <- unique(c(
    gmrf_param,
    "full", "project", "gmrf_project", "mvn_project", "separable", "projection"
  ))

  fit <- NULL
  fit_notes <- character()

  for (g in gmrf_candidates) {
    for (fam in list(NULL, "gaussian")) {
      fam_tag <- if (is.null(fam)) "default" else fam
      res <- tryCatch(
        fit_dsem_once(g, fam),
        error = function(e) e
      )
      if (!inherits(res, "error")) {
        fit <- res
        fit_notes <- c(fit_notes, paste0("fit succeeded with gmrf=", g, ", family=", fam_tag))
        break
      } else {
        fit_notes <- c(
          fit_notes,
          paste0("fit failed with gmrf=", g, ", family=", fam_tag, ": ", conditionMessage(res))
        )
      }
    }
    if (!is.null(fit)) break
  }

  if (is.null(fit)) {
    stop(
      paste(
        "DSEM fit failed for all gmrf/family combinations.",
        paste(fit_notes, collapse = " | ")
      )
    )
  }

  dsem_summary <- summary(fit)
  if (isTRUE(compute_loo)) {
    loo_tbl <- dsem::loo_residuals(fit, what = "loo", track_progress = FALSE) %>%
      as_tibble() %>%
      rename(
        Year = 1,
        Series = 2,
        Observed_scaled = obs,
        Predicted_scaled = est,
        SE_scaled = se
      ) %>%
      mutate(
        Year = as.numeric(Year),
        Series = as.character(Series)
      )
  } else {
    loo_tbl <- tibble(
      Year = numeric(),
      Series = character(),
      Observed_scaled = numeric(),
      Predicted_scaled = numeric(),
      SE_scaled = numeric()
    )
  }

  if (is.null(future_abc_raw)) {
    set.seed(seed)

    last_abc <- vapply(dsem_df[x_cols], last_finite, numeric(1))
    growth <- setNames(rep(-0.005, length(x_cols)), x_cols)
    dep_abc_cols <- intersect(paste0("ABC_", dep_safe), x_cols)
    if (length(dep_abc_cols) >= 1) growth[dep_abc_cols[1]] <- 0.02
    if (length(dep_abc_cols) >= 2) growth[dep_abc_cols[2]] <- 0.01

    noise_sd <- 0.03

    new_abc_raw <- sapply(seq_len(n_future), function(h) {
      eps <- rnorm(length(x_cols), mean = 0, sd = noise_sd)
      pmax(last_abc * (1 + growth) ^ h * (1 + eps), 1)
    })
    new_abc_raw <- t(new_abc_raw)
    colnames(new_abc_raw) <- x_cols
    forecast_years <- max(dsem_df$Year) + seq_len(n_future)
  } else {
    new_abc_raw <- future_abc_raw %>%
      as.data.frame() %>%
      select(all_of(x_cols))
    n_future <- nrow(new_abc_raw)
    forecast_years <- if (is.null(future_years)) {
      max(dsem_df$Year) + seq_len(n_future)
    } else {
      future_years
    }
  }

  new_abc_scaled <- scale_with(new_abc_raw, ts_center[x_cols], ts_scale[x_cols])

  newdata <- as.data.frame(tsdata)
  newdata_future <- as.data.frame(new_abc_scaled)
  newdata_future[, y_cols] <- NA_real_
  newdata <- bind_rows(newdata, newdata_future)
  newdata <- newdata[, colnames(as.data.frame(tsdata))]

  newdata_ts <- ts(as.matrix(newdata), start = start(tsdata), frequency = frequency(tsdata))

  pred_scaled <- predict(fit, newdata = newdata_ts, type = "response")
  colnames(pred_scaled) <- colnames(newdata_ts)
  pred <- sweep(pred_scaled, 2, ts_scale[colnames(pred_scaled)], "*")
  pred <- sweep(pred, 2, ts_center[colnames(pred)], "+")

  pred_tac_future <- tail(pred[, y_cols], n_future)
  forecast_tbl <- as_tibble(new_abc_raw) %>%
    mutate(Year = forecast_years) %>%
    bind_cols(as_tibble(pred_tac_future))

  list(
    fit = fit,
    summary = dsem_summary,
    dsem_df = dsem_df,
    dep_species = dep_species,
    y_cols = y_cols,
    x_cols = x_cols,
    pred = pred,
    ts_center = ts_center,
    ts_scale = ts_scale,
    loo_tbl = loo_tbl,
    n_future = n_future,
    forecast_years = forecast_years,
    forecast_tbl = forecast_tbl,
    fit_notes = fit_notes
  )
}

# Refit the DSEM sequentially through time and score one-step-ahead forecasts.
run_dsem_retrospective <- function(data, dep_species, indep_species, start_year = 2015, seed = 1) {
  full_df <- data %>%
    filter(Species %in% unique(c(dep_species, indep_species))) %>%
    group_by(Year, Species) %>%
    summarise(
      TAC = sum(TAC, na.rm = TRUE),
      ABC = sum(ABC, na.rm = TRUE),
      .groups = "drop"
    ) %>%
    pivot_wider(
      names_from = Species,
      values_from = c(TAC, ABC),
      names_sep = "_"
    ) %>%
    arrange(Year)

  names(full_df) <- safe_names(names(full_df))

  dep_safe <- safe_names(dep_species)
  indep_safe <- safe_names(indep_species)
  y_cols <- paste0("TAC_", dep_safe)
  x_cols <- paste0("ABC_", indep_safe)

  target_years <- full_df %>%
    filter(Year >= start_year) %>%
    pull(Year)

  purrr::map_dfr(seq_along(target_years), function(i) {
    target_year <- target_years[[i]]
    target_row <- full_df %>%
      filter(Year == target_year)

    if (nrow(target_row) != 1 || any(!is.finite(unlist(target_row[, x_cols, drop = FALSE])))) {
      return(tibble(
        Year = target_year,
        Training_end = target_year - 1,
        Species = dep_species,
        Observed = as.numeric(target_row[, y_cols, drop = FALSE]),
        Predicted = NA_real_,
        Error = NA_real_,
        Abs_Error = NA_real_,
        APE = NA_real_,
        Status = "missing_future_abc"
      ))
    }

    fit_i <- tryCatch(
      fit_dsem_region(
        data = data %>% filter(Year < target_year),
        dep_species = dep_species,
        indep_species = indep_species,
        n_future = 1,
        seed = seed + i,
        future_abc_raw = target_row[, x_cols, drop = FALSE],
        future_years = target_year,
        compute_loo = FALSE
      ),
      error = function(e) e
    )

    if (inherits(fit_i, "error")) {
      return(tibble(
        Year = target_year,
        Training_end = target_year - 1,
        Species = dep_species,
        Observed = as.numeric(target_row[, y_cols, drop = FALSE]),
        Predicted = NA_real_,
        Error = NA_real_,
        Abs_Error = NA_real_,
        APE = NA_real_,
        Status = paste0("fit_failed: ", conditionMessage(fit_i))
      ))
    }

    predicted <- fit_i$forecast_tbl %>%
      select(all_of(y_cols))

    tibble(
      Year = target_year,
      Training_end = target_year - 1,
      Species = dep_species,
      Observed = as.numeric(target_row[, y_cols, drop = FALSE]),
      Predicted = as.numeric(predicted[1, , drop = FALSE]),
      Status = "ok"
    ) %>%
      mutate(
        Error = Observed - Predicted,
        Abs_Error = abs(Error),
        APE = if_else(Observed != 0, 100 * Abs_Error / abs(Observed), NA_real_)
      )
  })
}

# Re-label raw DSEM coefficient output into management-facing term groups.
summarize_dsem_results <- function(dsem_summary) {
  dsem_summary %>%
    mutate(
      Term = if_else(lag > 0, paste0(path, " (lag ", lag, ")"), path),
      Group = case_when(
        lag > 0 ~ "Autoregressive (lagged TAC)",
        str_detect(first, "^ABC_") ~ "ABC -> TAC",
        TRUE ~ "Other"
      )
    ) %>%
    arrange(Group, Term) %>%
    select(Group, Term, Estimate, Std_Error, z_value, p_value)
}

# Render the future ABC inputs and forecast TAC outputs as a gt table.
render_dsem_forecast_table <- function(region_fit) {
  region_fit$forecast_tbl %>%
    select(Year, starts_with("ABC_"), starts_with("TAC_")) %>%
    gt::gt() %>%
    gt::fmt_number(columns = starts_with("ABC_"), decimals = 0) %>%
    gt::fmt_number(columns = starts_with("TAC_"), decimals = 0)
}

# Render the grouped DSEM coefficient summary as a gt table.
render_dsem_summary_table <- function(region_fit) {
  summarize_dsem_results(region_fit$summary) %>%
    gt::gt(rowname_col = "Term", groupname_col = "Group") %>%
    gt::cols_label(
      Estimate = "Estimate",
      Std_Error = "Std. Error",
      z_value = "z",
      p_value = "p"
    ) %>%
    gt::fmt_number(columns = c(Estimate, Std_Error), decimals = 3) %>%
    gt::fmt_number(columns = z_value, decimals = 2) %>%
    gt::fmt_number(columns = p_value, decimals = 3)
}

# Back-transform leave-one-out residuals into TAC units for one fitted region.
extract_dsem_loo_residuals <- function(region_fit, region_name) {
  region_fit$loo_tbl %>%
    filter(Series %in% region_fit$y_cols) %>%
    mutate(
      Region = region_name,
      Species = region_fit$dep_species[match(Series, region_fit$y_cols)],
      Observed = Observed_scaled * region_fit$ts_scale[Series] + region_fit$ts_center[Series],
      Predicted = Predicted_scaled * region_fit$ts_scale[Series] + region_fit$ts_center[Series],
      SE = SE_scaled * region_fit$ts_scale[Series],
      Residual = Observed - Predicted,
      Std_Residual = if_else(SE_scaled > 0, (Observed_scaled - Predicted_scaled) / SE_scaled, NA_real_)
    ) %>%
    select(Region, Year, Species, Observed, Predicted, SE, Residual, Std_Residual)
}

# Bind leave-one-out residual diagnostics across all fitted regions.
collect_dsem_loo_residuals <- function(region_fits) {
  purrr::imap_dfr(region_fits, extract_dsem_loo_residuals)
}

# Summarize leave-one-out residual performance metrics by region and species.
render_dsem_residual_comparison_table <- function(region_fits) {
  collect_dsem_loo_residuals(region_fits) %>%
    group_by(Region, Species) %>%
    summarise(
      Bias = mean(Residual, na.rm = TRUE),
      MAE = mean(abs(Residual), na.rm = TRUE),
      RMSE = sqrt(mean(Residual^2, na.rm = TRUE)),
      Mean_abs_z = mean(abs(Std_Residual), na.rm = TRUE),
      Max_abs_z = max(abs(Std_Residual), na.rm = TRUE),
      .groups = "drop"
    ) %>%
    arrange(Region, Species) %>%
    gt::gt(rowname_col = "Species", groupname_col = "Region") %>%
    gt::cols_label(
      Bias = "Bias",
      MAE = "MAE",
      RMSE = "RMSE",
      Mean_abs_z = "Mean |z|",
      Max_abs_z = "Max |z|"
    ) %>%
    gt::fmt_number(columns = c(Bias, MAE, RMSE), decimals = 0) %>%
    gt::fmt_number(columns = c(Mean_abs_z, Max_abs_z), decimals = 2)
}

# Plot standardized leave-one-out residuals through time by species and region.
plot_dsem_loo_residuals <- function(region_fits) {
  plot_df <- collect_dsem_loo_residuals(region_fits)

  ggplot(plot_df, aes(x = Year, y = Std_Residual)) +
    geom_hline(yintercept = 0, color = "grey60", linewidth = 0.4) +
    geom_hline(yintercept = c(-2, 2), color = "firebrick", linetype = "dashed", linewidth = 0.4) +
    geom_line(color = "steelblue", linewidth = 0.8) +
    geom_point(color = "steelblue", size = 1.5, alpha = 0.85) +
    facet_grid(Species ~ Region) +
    labs(x = NULL, y = "Standardized leave-one-out residual")
}

# Combine retrospective forecast results from each region into one data frame.
collect_dsem_retrospective <- function(retro_fits) {
  purrr::imap_dfr(retro_fits, ~ mutate(.x, Region = .y), .id = NULL)
}

# Summarize one-step-ahead retrospective forecast accuracy by region and species.
render_dsem_retrospective_summary_table <- function(retro_fits) {
  collect_dsem_retrospective(retro_fits) %>%
    filter(Status == "ok") %>%
    group_by(Region, Species) %>%
    summarise(
      n_years = n(),
      Bias = mean(Error, na.rm = TRUE),
      MAE = mean(Abs_Error, na.rm = TRUE),
      RMSE = sqrt(mean(Error^2, na.rm = TRUE)),
      MAPE = mean(APE, na.rm = TRUE),
      .groups = "drop"
    ) %>%
    arrange(Region, Species) %>%
    gt::gt(rowname_col = "Species", groupname_col = "Region") %>%
    gt::cols_label(
      n_years = "Years",
      Bias = "Bias",
      MAE = "MAE",
      RMSE = "RMSE",
      MAPE = "MAPE (%)"
    ) %>%
    gt::fmt_number(columns = n_years, decimals = 0) %>%
    gt::fmt_number(columns = c(Bias, MAE, RMSE), decimals = 0) %>%
    gt::fmt_number(columns = MAPE, decimals = 2)
}

# Plot observed versus retrospective one-step-ahead TAC predictions.
plot_dsem_retrospective <- function(retro_fits) {
  plot_df <- collect_dsem_retrospective(retro_fits) %>%
    filter(Status == "ok") %>%
    pivot_longer(c(Observed, Predicted), names_to = "Series", values_to = "TAC")

  ggplot(plot_df, aes(x = Year, y = TAC, color = Series)) +
    geom_line(linewidth = 0.9) +
    geom_point(size = 1.5, alpha = 0.85) +
    facet_grid(Species ~ Region, scales = "free_y") +
    scale_color_manual(values = c(Observed = "grey40", Predicted = "steelblue")) +
    scale_y_continuous(labels = comma) +
    labs(x = NULL, y = "TAC", color = NULL)
}

# Plot in-sample fitted TAC values plus the forward forecast extension.
plot_dsem_predictions <- function(region_fit) {
  pred_years <- c(region_fit$dsem_df$Year, region_fit$forecast_years)
  pred_df <- as_tibble(region_fit$pred[, region_fit$y_cols]) %>%
    mutate(Year = pred_years) %>%
    pivot_longer(-Year, names_to = "Species", values_to = "Predicted")

  obs_df <- region_fit$dsem_df %>%
    select(Year, all_of(region_fit$y_cols)) %>%
    pivot_longer(-Year, names_to = "Species", values_to = "Observed")

  plot_df <- left_join(pred_df, obs_df, by = c("Year", "Species"))

  ggplot(plot_df, aes(x = Year)) +
    geom_line(aes(y = Predicted, color = "Predicted"), linewidth = 1) +
    geom_point(aes(y = Observed, color = "Observed"), size = 1.4, alpha = 0.8) +
    facet_grid(Species ~ ., scales = "free_y") +
    scale_color_manual(values = c(Predicted = "steelblue", Observed = "grey40")) +
    scale_y_continuous(labels = comma) +
    labs(x = NULL, y = "TAC", color = NULL)
}

bsai_dsem_dep <- c("Pollock", "Pacific cod")
bsai_dsem_indep <- mainspp
bsai_dsem <- fit_dsem_region(
  data = dfOY %>% filter(lag == 1, Species %in% mainspp),
  dep_species = bsai_dsem_dep,
  indep_species = bsai_dsem_indep
)

goa_dsem_dep <- c("Pollock", "Pacific cod")
goa_dsem_indep <- goaraw %>%
  filter(lag == 1, OY == 1, Area == "Total") %>%
  group_by(Species) %>%
  summarise(total_TAC = sum(TAC, na.rm = TRUE), .groups = "drop") %>%
  filter(is.finite(total_TAC), total_TAC > 0) %>%
  arrange(desc(total_TAC)) %>%
  slice_head(n = 7) %>%
  pull(Species)

goa_dsem <- fit_dsem_region(
  data = goaOY %>% filter(lag == 1, Species %in% goa_dsem_indep),
  dep_species = goa_dsem_dep,
  indep_species = goa_dsem_indep
)

bsai_dsem_retro <- run_dsem_retrospective(
  data = dfOY %>% filter(lag == 1, Species %in% mainspp),
  dep_species = bsai_dsem_dep,
  indep_species = bsai_dsem_indep,
  start_year = 2015
)

goa_dsem_retro <- run_dsem_retrospective(
  data = goaOY %>% filter(lag == 1, Species %in% goa_dsem_indep),
  dep_species = goa_dsem_dep,
  indep_species = goa_dsem_indep,
  start_year = 2015
)

BSAI

plot_dsem_dag(
  dep_species = bsai_dsem_dep,
  indep_species = setdiff(bsai_dsem_indep, bsai_dsem_dep),
  dynamic = FALSE
)

Figure 17: BSAI contemporaneous DSEM DAG with ABC drivers for Pollock and Pacific cod TACs.

plot_dsem_dag(
  dep_species = bsai_dsem_dep,
  indep_species = setdiff(bsai_dsem_indep, bsai_dsem_dep),
  dynamic = TRUE
)

Figure 18: BSAI dynamic DSEM DAG with AR(1) TAC structure and contemporaneous ABC effects.

GOA

plot_dsem_dag(
  dep_species = goa_dsem_dep,
  indep_species = setdiff(goa_dsem_indep, goa_dsem_dep),
  dynamic = FALSE
)

Figure 19: GOA contemporaneous DSEM DAG with ABC drivers for Pollock and Pacific cod TACs.

plot_dsem_dag(
  dep_species = goa_dsem_dep,
  indep_species = setdiff(goa_dsem_indep, goa_dsem_dep),
  dynamic = TRUE
)

Figure 20: GOA dynamic DSEM DAG with AR(1) TAC structure and contemporaneous ABC effects.

DSEM estimates

In both regions, the DSEM estimates TAC dynamics for Pollock and Pacific cod with two components: (1) AR(1) persistence in each TAC time series and (2) contemporaneous influence of ABC inputs, including both focal-stock ABC and cross-stock ABCs. Models are fit on scaled inputs for stability, then predictions are transformed back to the original TAC units for interpretation. The forecast tables below remain illustrative and use the same simple forward ABC scenario in each region.

We also compare leave-one-out residuals from dsem::loo_residuals() for the fitted TAC series. These diagnostics show whether the dynamic model leaves systematic bias or large standardized errors for Pollock and Pacific cod in either region when each year is predicted from the rest of the time series.

BSAI

For BSAI, the DSEM uses the seven main stocks (Pollock, Yellowfin sole, Pacific cod, Atka mackerel, Northern rock sole, Flathead sole, Pacific ocean perch) as contemporaneous ABC inputs.

render_dsem_forecast_table(bsai_dsem)

Table 5: BSAI example DSEM forecast: assumed future ABCs and predicted TACs (illustrative scenario).

Year	ABC_Pollock	ABC_Yellowfin_sole	ABC_Pacific_cod	ABC_Atka_mackerel	ABC_Northern_rock_sole	ABC_Flathead_sole	ABC_Pacific_ocean_perch	TAC_Pollock	TAC_Pacific_cod
2027	2,142,230	267,768	152,126	96,297	158,990	85,114	36,927	1,375,625	133,466
2028	2,276,252	269,546	156,154	95,587	158,479	85,207	33,807	1,395,708	130,103
2029	2,348,122	263,289	159,097	93,559	159,703	87,930	37,025	1,406,115	130,085
2030	2,371,258	262,913	151,171	92,211	154,823	85,557	34,270	1,422,667	126,298
2031	2,329,333	264,287	168,992	89,797	156,103	85,391	34,199	1,422,630	127,743

render_dsem_summary_table(bsai_dsem)

Table 6: Summary of BSAI DSEM results.

	Estimate	Std. Error	z	p
ABC -> TAC
ABC_Atka_mackerel -> TAC_Pacific_cod	0.096	0.105	0.91	0.362
ABC_Atka_mackerel -> TAC_Pollock	0.074	0.105	0.70	0.485
ABC_Atka_mackerel <-> ABC_Atka_mackerel	0.932	0.104	8.94	0.000
ABC_Flathead_sole -> TAC_Pacific_cod	−0.102	0.111	−0.91	0.361
ABC_Flathead_sole -> TAC_Pollock	−0.074	0.115	−0.65	0.519
ABC_Flathead_sole <-> ABC_Flathead_sole	0.965	0.108	8.94	0.000
ABC_Northern_rock_sole -> TAC_Pacific_cod	0.193	0.134	1.44	0.150
ABC_Northern_rock_sole -> TAC_Pollock	−0.279	0.138	−2.02	0.043
ABC_Northern_rock_sole <-> ABC_Northern_rock_sole	0.935	0.104	8.95	0.000
ABC_Pacific_cod -> TAC_Pacific_cod	0.379	0.116	3.28	0.001
ABC_Pacific_cod -> TAC_Pollock	0.037	0.100	0.37	0.709
ABC_Pacific_cod <-> ABC_Pacific_cod	0.908	0.102	8.95	0.000
ABC_Pacific_ocean_perch -> TAC_Pacific_cod	0.237	0.125	1.90	0.057
ABC_Pacific_ocean_perch -> TAC_Pollock	−0.258	0.129	−2.00	0.045
ABC_Pacific_ocean_perch <-> ABC_Pacific_ocean_perch	0.985	0.110	8.95	0.000
ABC_Pollock -> TAC_Pacific_cod	−0.297	0.118	−2.52	0.012
ABC_Pollock -> TAC_Pollock	0.583	0.132	4.41	0.000
ABC_Pollock <-> ABC_Pollock	0.911	0.102	8.95	0.000
ABC_Yellowfin_sole -> TAC_Pacific_cod	−0.451	0.116	−3.89	0.000
ABC_Yellowfin_sole -> TAC_Pollock	−0.107	0.116	−0.92	0.356
ABC_Yellowfin_sole <-> ABC_Yellowfin_sole	0.893	0.100	8.95	0.000
Autoregressive (lagged TAC)
TAC_Pacific_cod -> TAC_Pacific_cod (lag 1)	0.595	0.106	5.59	0.000
TAC_Pollock -> TAC_Pollock (lag 1)	0.398	0.094	4.24	0.000
Other
TAC_Pacific_cod <-> TAC_Pacific_cod	0.489	0.055	8.93	0.000
TAC_Pollock <-> TAC_Pollock	0.503	0.056	8.94	0.000

plot_dsem_predictions(bsai_dsem)

Figure 21: Observed vs predicted BSAI TACs from the DSEM (back-transformed to original units).

GOA

For GOA, the DSEM uses the top seven GOA species by total TAC from the Area == "Total" series (Pollock, Arrowtooth Flounder, Pacific cod, Shallow-water Flatfish, Sablefish, Pacific Ocean Perch, Flathead Sole) as contemporaneous ABC inputs.

render_dsem_forecast_table(goa_dsem)

Table 7: GOA example DSEM forecast: assumed future ABCs and predicted TACs (illustrative scenario).

Year	ABC_Pollock	ABC_Arrowtooth_Flounder	ABC_Pacific_cod	ABC_Shallow_water_Flatfish	ABC_Sablefish	ABC_Pacific_Ocean_Perch	ABC_Flathead_Sole	TAC_Pollock	TAC_Pacific_cod
2025	200,656	119,306	31,778	57,933	47,374	38,548	40,890	198,563	32,057
2026	213,209	120,099	32,619	57,506	47,222	38,590	37,435	207,959	33,902
2027	219,941	117,311	33,234	56,286	47,587	39,823	40,998	216,448	34,420
2028	222,108	117,143	31,578	55,475	46,132	38,749	37,947	218,546	33,204
2029	218,181	117,756	35,301	54,023	46,514	38,673	37,869	215,595	34,875

render_dsem_summary_table(goa_dsem)

Table 8: Summary of GOA DSEM results.

	Estimate	Std. Error	z	p
ABC -> TAC
ABC_Arrowtooth_Flounder -> TAC_Pacific_cod	−0.234	NaN	NaN	NaN
ABC_Arrowtooth_Flounder -> TAC_Pollock	−0.166	NaN	NaN	NaN
ABC_Arrowtooth_Flounder <-> ABC_Arrowtooth_Flounder	1.229	NaN	NaN	NaN
ABC_Flathead_Sole -> TAC_Pacific_cod	0.065	NaN	NaN	NaN
ABC_Flathead_Sole -> TAC_Pollock	0.066	NaN	NaN	NaN
ABC_Flathead_Sole <-> ABC_Flathead_Sole	1.013	NaN	NaN	NaN
ABC_Pacific_Ocean_Perch -> TAC_Pacific_cod	−0.450	NaN	NaN	NaN
ABC_Pacific_Ocean_Perch -> TAC_Pollock	−0.088	NaN	NaN	NaN
ABC_Pacific_Ocean_Perch <-> ABC_Pacific_Ocean_Perch	0.940	NaN	NaN	NaN
ABC_Pacific_cod -> TAC_Pacific_cod	0.833	NaN	NaN	NaN
ABC_Pacific_cod -> TAC_Pollock	0.069	NaN	NaN	NaN
ABC_Pacific_cod <-> ABC_Pacific_cod	0.931	NaN	NaN	NaN
ABC_Pollock -> TAC_Pacific_cod	0.074	NaN	NaN	NaN
ABC_Pollock -> TAC_Pollock	0.829	NaN	NaN	NaN
ABC_Pollock <-> ABC_Pollock	0.948	NaN	NaN	NaN
ABC_Sablefish -> TAC_Pacific_cod	0.173	NaN	NaN	NaN
ABC_Sablefish -> TAC_Pollock	−0.009	NaN	NaN	NaN
ABC_Sablefish <-> ABC_Sablefish	0.996	NaN	NaN	NaN
ABC_Shallow_water_Flatfish -> TAC_Pacific_cod	0.140	NaN	NaN	NaN
ABC_Shallow_water_Flatfish -> TAC_Pollock	0.048	NaN	NaN	NaN
ABC_Shallow_water_Flatfish <-> ABC_Shallow_water_Flatfish	1.022	NaN	NaN	NaN
Autoregressive (lagged TAC)
TAC_Pacific_cod -> TAC_Pacific_cod (lag 1)	0.220	NaN	NaN	NaN
TAC_Pollock -> TAC_Pollock (lag 1)	0.170	NaN	NaN	NaN
Other
TAC_Pacific_cod <-> TAC_Pacific_cod	0.212	NaN	NaN	NaN
TAC_Pollock <-> TAC_Pollock	0.122	NaN	NaN	NaN

plot_dsem_predictions(goa_dsem)

Figure 22: Observed vs predicted GOA TACs from the DSEM (back-transformed to original units).

Leave-one-out residual comparisons

render_dsem_residual_comparison_table(list(BSAI = bsai_dsem, GOA = goa_dsem))

Table 9: Leave-one-out TAC residual summaries for the BSAI and GOA DSEM fits.

	Bias	MAE	RMSE	Mean |z|	Max |z|
BSAI
Pacific cod	−43	12,204	16,734	0.69	3.19
Pollock	65	57,453	73,064	0.76	2.71
GOA
Pacific cod	257	3,495	3,992	0.77	1.71
Pollock	−28	5,606	7,248	0.78	2.97

plot_dsem_loo_residuals(list(BSAI = bsai_dsem, GOA = goa_dsem))

Figure 23: Standardized leave-one-out TAC residuals for the BSAI and GOA DSEM fits.

Retrospective skill test from 2015 onward

For each year from 2015 onward, we refit the DSEM using only data available through the previous year and then forecast that year’s Pollock and Pacific cod TACs using the realized ABC vector for that year. This is a rolling-origin, leave-future-out test of one-step-ahead skill rather than an in-sample fit diagnostic.

render_dsem_retrospective_summary_table(list(BSAI = bsai_dsem_retro, GOA = goa_dsem_retro))

Table 10: Rolling one-step-ahead DSEM skill from 2015 onward, using only prior history to predict each year.

	Years	Bias	MAE	RMSE	MAPE (%)
BSAI
Pacific cod	12	12,225	16,196	20,495	8.97
Pollock	12	−8,727	70,259	80,496	5.24
GOA
Pacific cod	10	7,809	10,552	12,383	64.10
Pollock	10	5,438	12,209	13,748	7.50

plot_dsem_retrospective(list(BSAI = bsai_dsem_retro, GOA = goa_dsem_retro))

Figure 24: Rolling one-step-ahead DSEM forecasts versus observed TAC from 2015 onward.

Discussion

Projection updates under new data and assessment

Multi-year catch advice inevitably embeds a forecasting problem: projections are conditioned on assumptions about future recruitment, selectivity, mortality, and implementation error, all frozen at the time advice is issued (Mid-Atlantic Fishery Management Council 2012; International Council for the Exploration of the Sea (ICES) 2023). The literature above converges on a shared observation: when a new assessment assimilates additional data, the posterior state of the stock can shift abruptly relative to those earlier projections (Sánchez-Maroño, Dolder, and Needle 2021). The impact of updating data and assessment creates a discontinuity and reflects the consequence of evolving data streams and analyses.

Interim management procedure (MP) work reframes the problem by explicitly separating assessment cadence from advice cadence (Huynh et al. 2020; Klibansky et al. 2022). Rather than treating interim updates as ad-hoc corrections, these approaches formalize how limited new information (recent catches, survey indices, environmental indicators) should update catch advice between full assessments. The goal is not to eliminate snap—information shocks are unavoidable—but to bound it, so that revisions are predictable, explainable, and risk-consistent.

This framing is directly relevant to BSAI groundfish, where large stocks (pollock, cod, flatfish) dominate total ABC and where management stability is highly valued. Multi-year ABCs smooth year-to-year variability, but smoothing pushes uncertainty forward in time. When a subsequent assessment revises recruitment trajectories or selectivity patterns, the resulting adjustment can appear disproportionate relative to the modest interim changes that preceded it. The Canadian and U.S. guidance documents emphasize that this tension is intrinsic: stability and responsiveness trade off, and neither can be maximized simultaneously (Krohn et al. 2019; Mid-Atlantic Fishery Management Council 2012; South Atlantic Scientific and Statistical Committee 2022).

The ICES and MSE-oriented literature adds an important synthesis (International Council for the Exploration of the Sea (ICES) 2023). Geromont and Butterworth (Geromont and Butterworth 2015) and Walter et al. (Walter et al. 2023) show that systems relying on simple, pre-tested update rules often experience smaller and more interpretable snaps than systems that defer all learning to infrequent, high-dimensional assessments. In that sense, forecast-to-assimilation snap is best understood as a design property of the advice system, not merely an outcome of assessment error. For BSAI groundfish, this perspective motivates explicit evaluation of how multi-year projections, rollover rules, and interim updates distribute learning over time—shaping not only biological risk, but also the credibility and usability of scientific advice.

Appendix

Function catalog

This appendix documents the custom functions defined in this report. The two safe_names() definitions are intentionally duplicated so the regression and DSEM chunks can run independently during rendering.

function_catalog <- tibble::tribble(
  ~Workflow, ~Function, ~Scope, ~Purpose,
  "Exploratory plots", "lower_points_smooth", "top-level", "Helper for `GGally::ggpairs()` lower panels with points plus a loess smooth.",
  "Regression", "safe_names", "top-level", "Sanitizes species and column names into ASCII-safe identifiers for regression design matrices.",
  "Regression", "fit_tac_abc_regressions", "top-level", "Builds grouped TAC/ABC data, fits stock-specific log-linear regressions, and returns coefficient and prediction objects.",
  "Regression", "render_regression_table", "top-level", "Formats regression coefficients and fit statistics with `gt`.",
  "Regression", "plot_regression_predictions", "top-level", "Plots observed versus predicted TAC trajectories from the regression fits.",
  "GOA OY cap", "rescale_oy_cap", "top-level", "Implements weighted cap rescaling and returns the rescaled catch vector plus the solved multiplier.",
  "GOA OY cap", "total_gap", "nested in `rescale_oy_cap()`", "Objective helper used in the bisection search to match the aggregate OY cap.",
  "DSEM", "safe_names", "top-level", "Sanitizes species and state names used to build DSEM series names and path labels.",
  "DSEM", "centered_positions", "top-level", "Creates vertically centered node positions for static and dynamic DAG plots.",
  "DSEM", "plot_dsem_dag", "top-level", "Builds static or dynamic DAG visualizations of the DSEM structure with `ggdag`.",
  "DSEM", "choose_gmrf_param", "top-level", "Finds a `gmrf_parameterization` value accepted by the installed `dsem` version.",
  "DSEM", "scale_with", "top-level", "Applies stored centering and scaling constants to new predictor matrices.",
  "DSEM", "last_finite", "top-level", "Returns the last finite value in a series, used to seed the forward ABC scenario.",
  "DSEM", "fit_dsem_region", "top-level", "Prepares regional multivariate time series, fits the DSEM, computes forecasts, and caches leave-one-out residual output.",
  "DSEM", "fit_dsem_once", "nested in `fit_dsem_region()`", "Small wrapper that calls `dsem::dsem()` for one GMRF/family combination.",
  "DSEM skill testing", "run_dsem_retrospective", "top-level", "Runs a rolling one-step-ahead retrospective from 2015 onward by fitting each target year only on prior history.",
  "DSEM", "summarize_dsem_results", "top-level", "Re-labels and groups `summary(fit)` output into interpretable DSEM coefficient categories.",
  "DSEM", "render_dsem_forecast_table", "top-level", "Formats the illustrative future ABC/TAC scenario table.",
  "DSEM", "render_dsem_summary_table", "top-level", "Formats the DSEM coefficient summary table with `gt`.",
  "DSEM residuals", "extract_dsem_loo_residuals", "top-level", "Back-transforms TAC leave-one-out predictions and residuals for one regional fit.",
  "DSEM residuals", "collect_dsem_loo_residuals", "top-level", "Combines leave-one-out residual data across multiple regional DSEM fits.",
  "DSEM residuals", "render_dsem_residual_comparison_table", "top-level", "Summarizes bias, MAE, RMSE, and standardized residual magnitudes across regions and species.",
  "DSEM residuals", "plot_dsem_loo_residuals", "top-level", "Plots standardized leave-one-out TAC residuals through time.",
  "DSEM skill testing", "collect_dsem_retrospective", "top-level", "Combines rolling one-step-ahead results across regional retrospective runs.",
  "DSEM skill testing", "render_dsem_retrospective_summary_table", "top-level", "Summarizes rolling retrospective skill using bias, MAE, RMSE, and MAPE.",
  "DSEM skill testing", "plot_dsem_retrospective", "top-level", "Plots observed TAC against rolling one-step-ahead DSEM forecasts from 2015 onward.",
  "DSEM", "plot_dsem_predictions", "top-level", "Plots observed and predicted TAC values from fitted and forecast DSEM output."
)

function_catalog %>%
  gt::gt(rowname_col = "Function", groupname_col = "Workflow") %>%
  gt::cols_label(
    Scope = "Scope",
    Purpose = "Purpose"
  ) %>%
  gt::cols_width(
    Scope ~ gt::px(180),
    Purpose ~ gt::px(540)
  )

Table 11: Catalog of custom functions defined in the report.

	Scope	Purpose
Exploratory plots
lower_points_smooth	top-level	Helper for `GGally::ggpairs()` lower panels with points plus a loess smooth.
Regression
safe_names	top-level	Sanitizes species and column names into ASCII-safe identifiers for regression design matrices.
fit_tac_abc_regressions	top-level	Builds grouped TAC/ABC data, fits stock-specific log-linear regressions, and returns coefficient and prediction objects.
render_regression_table	top-level	Formats regression coefficients and fit statistics with `gt`.
plot_regression_predictions	top-level	Plots observed versus predicted TAC trajectories from the regression fits.
GOA OY cap
rescale_oy_cap	top-level	Implements weighted cap rescaling and returns the rescaled catch vector plus the solved multiplier.
total_gap	nested in `rescale_oy_cap()`	Objective helper used in the bisection search to match the aggregate OY cap.
DSEM
safe_names	top-level	Sanitizes species and state names used to build DSEM series names and path labels.
centered_positions	top-level	Creates vertically centered node positions for static and dynamic DAG plots.
plot_dsem_dag	top-level	Builds static or dynamic DAG visualizations of the DSEM structure with `ggdag`.
choose_gmrf_param	top-level	Finds a `gmrf_parameterization` value accepted by the installed `dsem` version.
scale_with	top-level	Applies stored centering and scaling constants to new predictor matrices.
last_finite	top-level	Returns the last finite value in a series, used to seed the forward ABC scenario.
fit_dsem_region	top-level	Prepares regional multivariate time series, fits the DSEM, computes forecasts, and caches leave-one-out residual output.
fit_dsem_once	nested in `fit_dsem_region()`	Small wrapper that calls `dsem::dsem()` for one GMRF/family combination.
summarize_dsem_results	top-level	Re-labels and groups `summary(fit)` output into interpretable DSEM coefficient categories.
render_dsem_forecast_table	top-level	Formats the illustrative future ABC/TAC scenario table.
render_dsem_summary_table	top-level	Formats the DSEM coefficient summary table with `gt`.
plot_dsem_predictions	top-level	Plots observed and predicted TAC values from fitted and forecast DSEM output.
DSEM skill testing
run_dsem_retrospective	top-level	Runs a rolling one-step-ahead retrospective from 2015 onward by fitting each target year only on prior history.
collect_dsem_retrospective	top-level	Combines rolling one-step-ahead results across regional retrospective runs.
render_dsem_retrospective_summary_table	top-level	Summarizes rolling retrospective skill using bias, MAE, RMSE, and MAPE.
plot_dsem_retrospective	top-level	Plots observed TAC against rolling one-step-ahead DSEM forecasts from 2015 onward.
DSEM residuals
extract_dsem_loo_residuals	top-level	Back-transforms TAC leave-one-out predictions and residuals for one regional fit.
collect_dsem_loo_residuals	top-level	Combines leave-one-out residual data across multiple regional DSEM fits.
render_dsem_residual_comparison_table	top-level	Summarizes bias, MAE, RMSE, and standardized residual magnitudes across regions and species.
plot_dsem_loo_residuals	top-level	Plots standardized leave-one-out TAC residuals through time.

DSEM package specifics

The DSEM portion of the report depends on package-specific behavior in dsem, especially around model specification, version compatibility, prediction, and residual diagnostics. The table below records the main package interfaces and how they are used here.

dsem_package_notes <- tibble::tribble(
  ~Topic, ~Package_interface, ~How_used_here, ~Why_it_matters,
  "Model fitting", "`dsem::dsem()`", "Fits each regional multivariate TAC-ABC time series after converting the data frame to a scaled `ts` object.", "This is the core model fit that produces latent states, coefficient summaries, and predictions.",
  "Version compatibility", "`dsem::dsem_control(gmrf_parameterization = ...)`", "The helper `choose_gmrf_param()` tests several candidate values because accepted names differ across installed package versions.", "Without this compatibility layer, the same report can fail locally or in CI when `dsem` changes its control vocabulary.",
  "SEM specification", "line-based SEM string", "Autoregressive TAC paths and contemporaneous ABC-to-TAC paths are concatenated into the `sem` string passed to `dsem()`.", "The report constructs the model explicitly rather than relying on a canned template, so the path structure is transparent in the source.",
  "State initialization", "`estimate_delta0 = TRUE`", "Enabled for both BSAI and GOA fits inside `fit_dsem_region()`.", "This lets the model estimate initial state offsets rather than forcing a fixed starting value.",
  "Standardization", "base `scale()` plus stored centers/scales", "All TAC and ABC series are fit on a common scaled basis, then predictions and leave-one-out outputs are transformed back to TAC units for display.", "Scaling improves optimizer stability and makes coefficients comparable across series with very different magnitudes.",
  "Prediction", "`predict(fit, newdata = ..., type = \"response\")`", "Used after appending future ABC scenarios with missing TAC values to obtain fitted-plus-forecast TAC paths.", "This provides a single prediction pass covering observed years and the illustrative forecast horizon.",
  "Retrospective skill", "rolling one-step-ahead refits", "For each year from 2015 onward, the report refits the DSEM on data through `year - 1` and predicts that year using the realized ABC vector.", "This is the leave-future-out skill test in the DSEM section and gives a stricter assessment than full-sample fitted values.",
  "Fallback fit strategy", "candidate GMRF and family loop", "The code tries multiple GMRF parameterizations and both default and Gaussian family settings before failing.", "This makes the document more robust to package and platform differences without changing the scientific specification.",
  "Coefficient summaries", "`summary(fit)`", "The raw summary output is regrouped into autoregressive and ABC-to-TAC terms before rendering.", "The package summary is machine-oriented; the report labels it in management-facing terms.",
  "LOO residual diagnostics", "`dsem::loo_residuals(fit, what = \"loo\")`", "The report uses the returned table of observed values, leave-one-out estimates, and standard errors, then back-transforms TAC series and computes standardized residuals.", "This provides a package-native cross-check on fit quality that is more stringent than in-sample residual inspection.",
  "Residual scale handling", "`obs`, `est`, and `se` columns from `loo_residuals()`", "Residual summaries are computed as TAC-unit errors after back-transformation, while standardized residuals remain on the scaled latent-state basis.", "Keeping both scales separates management-scale error magnitude from model-based z-score diagnostics."
)

dsem_package_notes %>%
  gt::gt(rowname_col = "Topic") %>%
  gt::cols_label(
    Package_interface = "Package interface",
    How_used_here = "How used here",
    Why_it_matters = "Why it matters"
  ) %>%
  gt::cols_width(
    Package_interface ~ gt::px(180),
    How_used_here ~ gt::px(300),
    Why_it_matters ~ gt::px(300)
  )

Table 12: DSEM package interfaces and implementation details used in this report.

	Package interface	How used here	Why it matters
Model fitting	`dsem::dsem()`	Fits each regional multivariate TAC-ABC time series after converting the data frame to a scaled `ts` object.	This is the core model fit that produces latent states, coefficient summaries, and predictions.
Version compatibility	`dsem::dsem_control(gmrf_parameterization = ...)`	The helper `choose_gmrf_param()` tests several candidate values because accepted names differ across installed package versions.	Without this compatibility layer, the same report can fail locally or in CI when `dsem` changes its control vocabulary.
SEM specification	line-based SEM string	Autoregressive TAC paths and contemporaneous ABC-to-TAC paths are concatenated into the `sem` string passed to `dsem()`.	The report constructs the model explicitly rather than relying on a canned template, so the path structure is transparent in the source.
State initialization	`estimate_delta0 = TRUE`	Enabled for both BSAI and GOA fits inside `fit_dsem_region()`.	This lets the model estimate initial state offsets rather than forcing a fixed starting value.
Standardization	base `scale()` plus stored centers/scales	All TAC and ABC series are fit on a common scaled basis, then predictions and leave-one-out outputs are transformed back to TAC units for display.	Scaling improves optimizer stability and makes coefficients comparable across series with very different magnitudes.
Prediction	`predict(fit, newdata = ..., type = "response")`	Used after appending future ABC scenarios with missing TAC values to obtain fitted-plus-forecast TAC paths.	This provides a single prediction pass covering observed years and the illustrative forecast horizon.
Retrospective skill	rolling one-step-ahead refits	For each year from 2015 onward, the report refits the DSEM on data through `year - 1` and predicts that year using the realized ABC vector.	This is the leave-future-out skill test in the DSEM section and gives a stricter assessment than full-sample fitted values.
Fallback fit strategy	candidate GMRF and family loop	The code tries multiple GMRF parameterizations and both default and Gaussian family settings before failing.	This makes the document more robust to package and platform differences without changing the scientific specification.
Coefficient summaries	`summary(fit)`	The raw summary output is regrouped into autoregressive and ABC-to-TAC terms before rendering.	The package summary is machine-oriented; the report labels it in management-facing terms.
LOO residual diagnostics	`dsem::loo_residuals(fit, what = "loo")`	The report uses the returned table of observed values, leave-one-out estimates, and standard errors, then back-transforms TAC series and computes standardized residuals.	This provides a package-native cross-check on fit quality that is more stringent than in-sample residual inspection.
Residual scale handling	`obs`, `est`, and `se` columns from `loo_residuals()`	Residual summaries are computed as TAC-unit errors after back-transformation, while standardized residuals remain on the scaled latent-state basis.	Keeping both scales separates management-scale error magnitude from model-based z-score diagnostics.

DSEM object contents used by the report

For the regional DSEM fits, the report retains a small structured list rather than passing the raw dsem fit object through every downstream chunk. The most important retained components are:

• fit: the original fitted dsem object, needed for any additional package-native inspection.
• summary: the coefficient summary returned by summary(fit).
• dsem_df, y_cols, and x_cols: the regional data frame and the TAC/ABC series names used to build the model.
• pred, forecast_years, and forecast_tbl: fitted values plus the illustrative forward ABC scenario and resulting TAC forecasts.
• ts_center and ts_scale: the centering and scaling constants required to back-transform predictions and leave-one-out residuals.
• loo_tbl: the leave-one-out residual table returned by dsem::loo_residuals(..., what = "loo").
• fit_notes: a log of which GMRF/family combinations failed or succeeded during fitting.

References

Asparouhov, Tihomir, Ellen L. Hamaker, and Bengt O. Muthén. 2018. “Dynamic Structural Equation Models.” Structural Equation Modeling 25 (3): 359–88. https://doi.org/10.1080/10705511.2017.1406803.

Faig, Amanda, and Alan Haynie. 2020. The ATTACH Model (Catchfunction Package) Vignette. https://github.com/amandafaig/catchfunction/blob/master/doc/ATTACH_vignette.pdf.

Geromont, H. F., and D. S. Butterworth. 2015. “Complex Assessments or Simple Management Procedures for Efficient Fisheries Management.” ICES Journal of Marine Science 72: 262–74. https://doi.org/10.1093/icesjms/fsu017.

Huynh, Q. C., A. E. Punt, T. R. Carruthers, and C. M. Dichmont. 2020. “The Interim Management Procedure Approach for Assessed Stocks: Responsive Management Advice and Lower Assessment Frequency.” Fish and Fisheries 21 (3): 663–79. https://doi.org/10.1111/faf.12453.

International Council for the Exploration of the Sea (ICES). 2023. “ICES Advice Framework and Advice Sheets.” International Council for the Exploration of the Sea. https://doi.org/10.17895/ices.advice.22116890.

Klibansky, Nikolai, Cassidy Peterson, Kyle Shertzer, Matthew Vincent, and Erik Williams. 2022. “Evaluating Procedures for Updating Catch Advice Between Stock Assessments of Reef Fishes with Management Strategy Evaluation.” NOAA Technical Memorandum NMFS-SEFSC-779. Washington, DC: National Oceanic; Atmospheric Administration, National Marine Fisheries Service, Southeast Fisheries Science Center. https://doi.org/10.25923/7q4k-0b32.

Krohn, M. M., G. Chaput, D. E. Duplisea, N. M. T. Duprey, A. M. Edwards, B. P. Healey, K. L. Howland, B. Lester, M. J. Morgan, and R. F. Tallman. 2019. “Guidelines for Providing Interim-Year Updates and Science Advice for Multi-Year Assessments.” DFO Canadian Science Advisory Secretariat Research Document 2018/035. Ottawa, ON: Fisheries; Oceans Canada. https://publications.gc.ca/collections/collection_2019/mpo-dfo/fs70-5/Fs70-5-2018-035-eng.pdf.

Mid-Atlantic Fishery Management Council. 2012. “Considerations for Setting Multi-Year Acceptable Biological Catch Limits.” Mid-Atlantic Fishery Management Council. https://www.mafmc.org/s/multi-year-ABC-considerations-3-11-12.pdf.

Sánchez-Maroño, S., P. J. Dolder, and C. L. Needle. 2021. “Effects of Advice Lag and Update Frequency on Fisheries Management Performance.” ICES Journal of Marine Science 78: 2495–2507.

South Atlantic Scientific and Statistical Committee. 2022. “Catch Level Projections Workgroup Report.” SEDAR.

Thorson, J. T., A. Andrews, T. Essington, and S. Large. 2024. “Dynamic Structural Equation Models Synthesize Ecosystem Dynamics Constrained by Ecological Mechanisms.” Methods in Ecology and Evolution. https://doi.org/10.1111/2041-210X.14289.

Walter, J. F., A. E. Punt, R. I. C. C. Francis, et al. 2023. “When to Conduct, and When Not to Conduct, Management Strategy Evaluation.” Fish and Fisheries 24: 1–18.

Zellner, Arnold. 1962. “An Efficient Method of Estimating Seemingly Unrelated Regressions and Tests for Aggregation Bias.” Journal of the American Statistical Association 57 (298): 348–68. https://doi.org/10.1080/01621459.1962.10480664.