Mean Length Estimation Using CPUE and Size Frequency Data

June 2026 · Seattle, WA

---

1 Introduction

This document describes the processing of Japanese longline size-frequency and CPUE data to estimate catch size composition using both nominal and model-derived CPUE. The model-based index is derived from a generalized additive model (GAM).

2 Data Input and Cleaning

As a first step, files are read in and labelled consistently.

knitr::opts_chunk$set(
  echo = TRUE,
  message = FALSE,
  warning = FALSE
)
suppressPackageStartupMessages({
  library(tidyverse)
  library(here)
})

clean_names <- function(x) {
  nm <- names(x)
  nm <- gsub("[^[:alnum:]]+", "_", tolower(nm))
  nm <- gsub("^_+|_+$", "", nm)
  names(x) <- make.unique(nm, sep = "_")
  x
}

size_file_candidates <- c(
  Sys.getenv("JP_SIZE_FILE", unset = NA_character_),
  here("doc", "JP_Size_Age4plus.csv"),
  here("scratch", "JP_Size_Age4plus.csv"),
  here("JP_Size_Age4plus.csv"),
  "JP_Size_Age4plus.csv"
)
size_file <- size_file_candidates[!is.na(size_file_candidates) & file.exists(size_file_candidates)][1]

model_file_candidates <- c(
  Sys.getenv("JP_MODEL_EST_FILE", unset = NA_character_),
  here("doc", "MeanPredictGAM2022_YMLL.csv"),
  here("scratch", "MeanPredictGAM2022_YMLL.csv"),
  here("MeanPredictGAM2022_YMLL.csv"),
  "MeanPredictGAM2022_YMLL.csv"
)
model_file <- model_file_candidates[!is.na(model_file_candidates) & file.exists(model_file_candidates)][1]

if (is.na(size_file)) {
  stop(
    "Could not find JP_Size_Age4plus.csv. Set JP_SIZE_FILE or copy data to ",
    here("doc", "JP_Size_Age4plus.csv"),
    call. = FALSE
  )
}

if (is.na(model_file)) {
  stop(
    "Could not find MeanPredictGAM2022_YMLL.csv. Set JP_MODEL_EST_FILE or copy data to ",
    here("doc", "MeanPredictGAM2022_YMLL.csv"),
    call. = FALSE
  )
}

df1 <- read_csv(size_file, show_col_types = FALSE) |> clean_names()
names(df1) <- c("year", "month", "lat5", "lon5", "length", "freq", "prec", "hooks")

df2 <- read_csv(model_file, show_col_types = FALSE) |>
  clean_names() |>
  rename(value = mean_pred1)
latest_data_year <- min(max(df1$year, na.rm = TRUE), max(df2$year, na.rm = TRUE))

assessment_length_freq_file <- here(
  "ESC31", "csv_2026", "lf_assessment.csv"
)
assessment_length_freq <- read_csv(
  assessment_length_freq_file,
  show_col_types = FALSE
)
latest_assessment_year <- max(assessment_length_freq$Year, na.rm = TRUE)

The Japanese size file JP_Size_Age4plus.csv runs through 2025. The GAM CPUE file MeanPredictGAM2022_YMLL.csv also runs through 2025, so the joined GAM-weighted estimates run through 2025. The assessment input ESC31/csv_2026/lf_assessment.csv runs through 2025.

2.1 Spatial Adjustments and Filtering

The data are filtered to include only records after 1968. Spatial coordinates are adjusted by shifting latitude and longitude by -2.5 and +2.5 degrees, respectively:

df1 <- df1 |>
  filter(year > 1968) |>
  mutate(lat5 = lat5 - 2.5, lon5 = lon5 + 2.5)

df2 <- df2 |>
  filter(year > 1968) |>
  mutate(jap_cpue = value)

2.2 Normalized Length Frequencies by Cell

Duplicate records for the same year, month, spatial cell, and length are first collapsed by summing their frequencies. The collapsed frequencies are then normalized within each spatiotemporal cell:

duplicate_cell_length_keys <- df1 |>
  count(year, month, lat5, lon5, length, name = "records") |>
  filter(.data$records > 1L)

cell_metadata_conflicts <- df1 |>
  summarise(
    hook_values = n_distinct(.data$hooks),
    precision_values = n_distinct(.data$prec),
    .by = c(year, month, lat5, lon5)
  ) |>
  filter(.data$hook_values != 1L | .data$precision_values != 1L)
if (nrow(cell_metadata_conflicts)) {
  stop(
    "Hooks and precision must each be constant within a spatiotemporal cell.",
    call. = FALSE
  )
}

df1.1 <- df1 |>
  summarise(
    freq = sum(.data$freq),
    hooks = first(.data$hooks),
    prec = first(.data$prec),
    .by = c(year, month, lat5, lon5, length)
  ) |>
  mutate(
    cell_frequency = sum(.data$freq),
    prop = .data$freq / .data$cell_frequency,
    .by = c(year, month, lat5, lon5)
  )

stopifnot(
  all(abs(
    df1.1 |>
      summarise(total = sum(.data$prop), .by = c(year, month, lat5, lon5)) |>
      pull(.data$total) -
      1
  ) < 1e-12)
)

\[ \text{prop}_{l} = \frac{\text{freq}_{l}}{\sum_{l} \text{freq}_{l}} \]

2.3 Merge with CPUE Data

dfj <- df1.1 |>
  inner_join(
    df2,
    by = join_by(year, month, lat5, lon5),
    relationship = "many-to-one"
  )

3 Size composition from CPUE-weighted Frequencies

We compute length-frequency estimates for each year-month-area cell (i), weighted by:

• Nominal CPUE:

\[ u_i = \frac{\text{catch}_i}{\text{hooks}_i} \]

(by year, month, and 5x5 degree spatial block)

• Model-based CPUE: from the GAM model output we have the predicted CPUE by year, month, 5x5 cell to get \(\text{CPUE}_i\). See jap_cpue code in previous section.

Then for each length bin:

\[ \bar{lf}^{\text{nom}}_l = \sum_i u_i \cdot \text{prop}_{il}, \quad \]

\[ \bar{lf}^{\text{gam}}_l = \sum_i \text{CPUE}_i \cdot \text{prop}_{il} \]

And annual proportions:

\[ p_l^{\text{nom}} = \frac{lf_l^{\text{nom}}}{\sum_l lf_l^{\text{nom}}}, \quad p_l^{\text{gam}} = \frac{lf_l^{\text{gam}}}{\sum_l lf_l^{\text{gam}}} \]

Mean lengths (for a given year; index dropped for clarity):

\[ \bar{L}_{\text{nom}} = \sum_l l \cdot p_l^{\text{nom}}, \quad \bar{L}_{\text{gam}} = \sum_l l \cdot p_l^{\text{gam}} \]

As a diagnostic to compare with the GAM CPUE, we can also compute mean lengths from the assessment model’s length frequency data. This is done by summing the product of length and proportion for each year. Comparing these values shows that the mean lengths from the GAM CPUE data and the nominal CCSBT data are similar, but not identical. The GAM CPUE is a model-based estimate that may differ from the nominal CPUE due to smoothing and other adjustments in the GAM (Figure 1). The GAM-weighted and nominal compositions differ more substantially from the catch-at-length data used elsewhere in the CCSBT assessment. The GAM-weighted composition is therefore the more internally consistent choice for estimating the selectivity associated with the Japanese CPUE index.

mnlen_cpue <- dfj |>
  mutate(
    u = sum(.data$freq) / first(.data$hooks),
    .by = c(year, month, lat5, lon5)
  ) |>
  summarise(
    lf_u_nom = sum(.data$u * .data$prop),
    lf_u_gam = sum(.data$jap_cpue * .data$prop),
    .by = c(year, length)
  ) |>
  mutate(
    p_u_nom = .data$lf_u_nom / sum(.data$lf_u_nom),
    p_u_gam = .data$lf_u_gam / sum(.data$lf_u_gam),
    .by = year
  ) |>
  summarise(
    mean_len_u_nom = sum(.data$length * .data$p_u_nom),
    mean_len_u_gam = sum(.data$length * .data$p_u_gam),
    .by = year
  ) |>
  pivot_longer(
    cols = starts_with("mean_len"),
    names_to = "type",
    values_to = "Length"
  ) |>
  mutate(
    type = recode(
      .data$type,
      mean_len_u_nom = "Nominal CPUE",
      mean_len_u_gam = "GAM CPUE"
    )
  )

mnlen_LL1 <- assessment_length_freq |>
  filter(Year > 1968, Fishery == 1) |>
  pivot_longer(
    cols = -c(Fishery, Year, N),
    names_to = "len",
    values_to = "proportion"
  ) |>
  filter(proportion > 0) |>
  mutate(year = Year, len = as.numeric(len)) |>
  summarise(
    type = "Assessment LL1",
    Length = sum(.data$len * .data$proportion),
    .by = year
  )

bind_rows(mnlen_cpue, mnlen_LL1) |>
  ggplot(aes(x = .data$year, y = .data$Length, color = .data$type)) +
  geom_point() +
  geom_line() +
  labs(x = "Year", y = "Mean length (cm)", color = NULL) +
  ggthemes::theme_few()

Figure 1: Annual mean SBT length from the accepted assessment LL1, nominal CPUE-weighted, and GAM CPUE-weighted length compositions.

3.1 Length Frequency by GAM CPUE

To compute and optionally visualize full length frequency distributions weighted by GAM-predicted CPUE:

Length-frequency proportions by year are computed as:

\[ \text{prop}^{\text{gam}}_{l,\text{year}} = \frac{\sum_i \text{CPUE}_i \cdot \text{prop}_{il}}{\sum_l \sum_i \text{CPUE}_i \cdot \text{prop}_{il}} \]

Notes:

• floor(length / 2) * 2 groups lengths 40 and 41 into bin 40, 42 and 43 into bin 42, and so on.

• prop_gam is the proportion in each 2-cm bin within each year:

\[ \text{prop}^{\text{gam}}_{b, \text{year}} = \frac{\text{freq}^{\text{gam}}_b}{\sum_b \text{freq}^{\text{gam}}_b} \]

• Use length_bin for plotting or comparison with model outputs.

3.2 Format to bring into the assessment as a new set of length-composition data

We write a matrix with years in rows and length bins in columns. For comparisons, we show the length frequency (in proportions) for the GAM CPUE data (Figure 2) and the accepted assessment Fishery 1 (LL1) length-composition input (Figure 3). The latter is a contextual catch-composition comparison, not the frozen Fishery 7 CPUE composition used for the duplicate-key audit below.

# str(lf_gam)
sbt_length_bins <- as.character(seq(32, 250, by = 2))
lf_gam_sbt_wide <- lf_gam |>
  mutate(length_bin = pmin(pmax(length_bin, 32), 250)) |>
  group_by(year, length_bin) |>
  summarise(prop_gam = sum(prop_gam), .groups = "drop") |>
  pivot_wider(
    names_from = length_bin,
    values_from = prop_gam,
    values_fill = 0
  )

accepted_cpue_lf <- assessment_length_freq |>
  filter(.data$Fishery == 7L) |>
  arrange(.data$Year)
if (!nrow(accepted_cpue_lf)) {
  stop("The accepted assessment input has no Fishery 7 CPUE compositions.")
}

lf_gam_sbt_fishery_corrected <- lf_gam_sbt_wide |>
  rename(Year = year) |>
  mutate(Fishery = 7L, .before = Year) |>
  left_join(
    accepted_cpue_lf |>
      select(.data$Year, .data$N),
    by = "Year",
    relationship = "one-to-one"
  ) |>
  relocate(N, .after = Year)

missing_bins <- setdiff(
  sbt_length_bins,
  names(lf_gam_sbt_fishery_corrected)
)
lf_gam_sbt_fishery_corrected[missing_bins] <- 0
lf_gam_sbt_fishery_corrected <- lf_gam_sbt_fishery_corrected |>
  select(Fishery, Year, N, all_of(sbt_length_bins)) |>
  arrange(Fishery, Year)

accepted_cpue_lf <- accepted_cpue_lf |>
  select(Fishery, Year, N, all_of(sbt_length_bins))
stopifnot(identical(
  lf_gam_sbt_fishery_corrected$Year,
  accepted_cpue_lf$Year
))
composition_difference <- as.matrix(
  lf_gam_sbt_fishery_corrected[sbt_length_bins]
) - as.matrix(accepted_cpue_lf[sbt_length_bins])
max_composition_difference <- max(abs(composition_difference))
material_difference_years <- accepted_cpue_lf$Year[
  apply(abs(composition_difference) > 1e-6, 1, any)
]

# Publish the frozen, accepted ESC31 input by preserving its original CSV
# fields. The corrected duplicate-key aggregation above is intentionally
# diagnostic until its one affected year is accepted as a production-data
# change.
assessment_length_freq_lines <- readLines(
  assessment_length_freq_file,
  warn = FALSE
)
accepted_cpue_lf_lines <- c(
  assessment_length_freq_lines[[1L]],
  assessment_length_freq_lines[
    -1L
  ][startsWith(assessment_length_freq_lines[-1L], "7,")]
)
stopifnot(length(accepted_cpue_lf_lines) == nrow(accepted_cpue_lf) + 1L)
accepted_cpue_lf_file <- here(
  "doc", "lf_gam_cpue_sbt_fishery_length_freq.csv"
)
writeLines(
  accepted_cpue_lf_lines,
  accepted_cpue_lf_file,
  useBytes = TRUE
)
published_cpue_lf <- read_csv(
  accepted_cpue_lf_file,
  show_col_types = FALSE
)
stopifnot(isTRUE(all.equal(
  accepted_cpue_lf,
  published_cpue_lf,
  tolerance = 0,
  check.attributes = FALSE
)))

flen <- 80
llen <- 190
library(ggridges)
p1 <- ggplot(
  lf_gam,
  aes(
    x = .data$length_bin,
    y = as.factor(.data$year),
    height = .data$prop_gam
  )
) +
  geom_density_ridges(
    stat = "identity",
    scale = 4,
    alpha = 0.7,
    fill = "salmon",
    color = "black"
  ) +
  ggthemes::theme_few() +
  labs(x = "Length (2-cm bin)", y = "Year") +
  scale_x_continuous(
    limits = c(flen, llen),
    breaks = seq(flen, llen, 6)
  ) +
  scale_y_discrete(limits = rev(levels(as.factor(lf_gam$year))))
p1

Figure 2: Length frequency matrix from GAM CPUE data.

The corrected GAM CPUE matrix in Figure 2 collapses 23 duplicate cell-length keys. Relative to the frozen accepted ESC31 Fishery 7 input, the correction exceeds \(10^{-6}\) in 1996 only, with a maximum absolute bin-proportion difference of 0.000529306. The downloadable lf_gam_cpue_sbt_fishery_length_freq.csv therefore remains an exact copy of the accepted assessment input; the deduplicated result is a documented pre-sign-off data decision and is not silently substituted into the fitted model.

lf_ass <- assessment_length_freq |>
  filter(Year > 1968, Fishery == 1) |>
  pivot_longer(
    cols = -c(Fishery, Year, N),
    names_to = "length_bin",
    values_to = "proportion"
  )
p2 <- ggplot(
  lf_ass,
  aes(
    x = as.numeric(.data$length_bin),
    y = as.factor(.data$Year),
    height = .data$proportion
  )
) +
  geom_density_ridges(
    stat = "identity",
    scale = 4,
    alpha = 0.7,
    fill = "salmon",
    color = "black"
  ) +
  ggthemes::theme_few() +
  labs(x = "Length (2-cm bin)", y = "Year") +
  scale_x_continuous(
    limits = c(flen, llen),
    breaks = seq(flen, llen, 6)
  ) +
  scale_y_discrete(limits = rev(levels(as.factor(lf_ass$Year))))
p2

Figure 3: Length frequency matrix from the accepted assessment Fishery 1 (LL1) input.

4 Summary

This workflow merges catch-at-length observations with spatially resolved GAM predictions to estimate trends in average catch size. The final product includes corrected mean-length diagnostics and full length-frequency distributions weighted by GAM CPUE. The accepted Fishery 7 composition remains frozen pending the explicit decision on the small 1996 duplicate-key correction described above.

---