AFSC Case Study Gulf of Alaska Walleye Pollock

The setup

Code 
common_name <- "Gulf of Alaska walleye pollock"

# define the dimensions and global variables
years <- 1970:2023
n_years <- length(years)
ages <- 1:10
n_ages <- length(ages)
  • R version: R version 4.5.3 (2026-03-11)
  • TMB version: 1.9.21
  • FIMS commit: 5666d0b
  • Stock name: Gulf of Alaska walleye pollock
  • Region: AFSC
  • Analyst: Cole Monnahan

Simplifications to the original assessment

The model presented in this case study was changed substantially from the operational version and should not be considered reflective of the Gulf of Alaska walleye pollock stock. These results are intended to be a demonstration and nothing more.

To get the operational model to more closely match a FIMS model the following changes were made:

  • Dropped surveys 1, 4, and 5
  • Removed ageing error
  • Removed length compositions
  • Removed bias correction in log-normal index likelihoods
  • Simplified catchability of survey 3 to be constant in time (removed random walk)
  • Updated maturity to be parametric rather than empirical
  • Used constant weight at age for all sources: spawning, fishery, surveys, and biomass calculations. The same matrix was used throughout.
  • Change timing to be Jan 1 for spawning and all surveys
  • Removed prior on catchability for survey 2
  • Removed time-varying fisheries selectivity (constant double logistic)
  • Took off normalization of selectivity
  • Removed age accumulation for fishery age compositions

Setting up the data

Code 
#
## This will fit the models bridging to FIMS (simplifying)
## source("fit_bridge_models.R")
## compare changes to model
pkfitfinal <- readRDS(file.path(data_directory, "pkfitfinal.RDS"))
pkfit0 <- readRDS(file.path(data_directory, "pkfit0.RDS"))
rep0 <- pkfitfinal$rep
parfinal <- pkfitfinal$obj$env$parList()
pkinput0 <- readRDS(file.path(data_directory, "pkinput0.RDS"))
fimsdat <- pkdat0 <- pkinput0$dat
pkinput <- readRDS(file.path(data_directory, "pkinput.RDS"))
data_4_model <- prepare_pollock_data(
  pkfitfinal = pkfitfinal,
  pkfit0 = pkfit0,
  parfinal = parfinal,
  fimsdat = fimsdat,
  pkinput = pkinput,
  years = years,
  n_years = n_years,
  ages = ages,
  n_ages = n_ages
)

Run FIMS model

Code 
FIMS::clear()

estimate_fish_selex <- TRUE
estimate_survey_selex <- TRUE
estimate_q2 <- TRUE
estimate_q3 <- TRUE
estimate_q6 <- TRUE
estimate_F <- TRUE
estimate_recdevs <- TRUE

# Create parameters list and update default values for parameters
default_parameters <- FIMS::create_default_configurations(
  data = data_4_model
) |>
  tidyr::unnest(cols = data) |>
  dplyr::rows_update(
    tibble::tibble(
      module_name = "Selectivity",
      fleet_name = c("fleet1", "survey2", "survey6"),
      module_type = "DoubleLogistic"
    ),
    by = c("module_name", "fleet_name")
  ) |>
  FIMS::create_default_parameters(
    data = data_4_model
  ) |>
  tidyr::unnest(cols = data) |>
  dplyr::rows_update(
    tibble::tibble(
      module_name = "Selectivity",
      fleet_name = "fleet1",
      label = c(
        "inflection_point_asc", "slope_asc",
        "inflection_point_desc", "slope_desc"
      ),
      value = c(
        parfinal$inf1_fsh_mean,
        exp(parfinal$log_slp1_fsh_mean),
        parfinal$inf2_fsh_mean,
        exp(parfinal$log_slp2_fsh_mean)
      ),
      estimation_type = "fixed_effects"
    ),
    by = c("module_name", "fleet_name", "label")
  ) |>
  dplyr::rows_update(
    tibble::tibble(
      module_name = "Selectivity",
      fleet_name = "survey2",
      label = c(
        "inflection_point_asc", "slope_asc",
        "inflection_point_desc", "slope_desc"
      ),
      value = c(
        parfinal$inf1_srv2,
        exp(parfinal$log_slp1_srv2),
        parfinal$inf2_srv2,
        exp(parfinal$log_slp2_srv2)
      ),
      estimation_type = rep(c("fixed_effects", "constant"), each = 2)
    ),
    by = c("module_name", "fleet_name", "label")
  ) |>
  dplyr::rows_update(
    tibble::tibble(
      module_name = "Selectivity",
      fleet_name = "survey6",
      label = c(
        "inflection_point_asc", "slope_asc",
        "inflection_point_desc", "slope_desc"
      ),
      value = c(
        parfinal$inf1_srv6,
        exp(parfinal$log_slp1_srv6),
        parfinal$inf2_srv6,
        exp(parfinal$log_slp2_srv6)
      ),
      estimation_type = rep(c("constant", "fixed_effects"), each = 2)
    ),
    by = c("module_name", "fleet_name", "label")
  ) |>
  dplyr::rows_update(
    tibble::tibble(
      module_name = "Selectivity",
      fleet_name = "survey3",
      label = c("inflection_point", "slope"),
      value = c(
        parfinal$inf1_srv3,
        exp(parfinal$log_slp1_srv3)
      )
    ),
    by = c("module_name", "fleet_name", "label")
  ) |>
  dplyr::rows_update(
    tibble::tibble(
      module_name = "Maturity",
      label = c("inflection_point", "slope"),
      value = c(4.5, 1.5)
    ),
    by = c("module_name", "label")
  ) |>
  dplyr::rows_update(
    tibble::tibble(
      module_name = "Recruitment",
      label = c(
        "log_rzero",
        "logit_steep",
        "log_sd"
      ),
      value = c(
        parfinal$mean_log_recruit + log(1e9),
        FIMS::logit(0.2, 1.0, 0.99999),
        parfinal$sigmaR
      )
    ),
    by = c("module_name", "label")
  ) |>
  dplyr::rows_update(
    tibble::tibble(
      module_name = "Population",
      label = c("log_init_naa"),
      age = ages,
      value = c(log(pkfitfinal$rep$recruit[1]), log(pkfitfinal$rep$initN)) + log(1e9),
      estimation_type = "constant"
    ),
    by = c("module_name", "age", "label")
  ) |>
  dplyr::rows_update(
    tibble::tibble(
      module_name = "Population",
      time = rep(years, each = FIMS::get_n_ages(data_4_model)),
      age = rep(ages, FIMS::get_n_years(data_4_model)),
      label = "log_M",
      value = log(as.numeric(t(matrix(
        rep(pkfitfinal$rep$M, each = FIMS::get_n_years(data_4_model)),
        nrow = FIMS::get_n_years(data_4_model)
      ))))
    ),
    by = c("module_name", "label", "time", "age")
  ) |>
  dplyr::rows_update(
    tibble::tibble(
      module_name = "Recruitment",
      label = "log_devs",
      time = years[-1],
      # The last value of the initial numbers at age is the first
      # recruitment deviation
      value = parfinal$dev_log_recruit[-1],
    ),
    by = c("module_name", "label", "time")
  ) |>
  dplyr::rows_update(
    tibble::tibble(
      module_name = "Fleet",
      fleet_name = "fleet1",
      time = years,
      label = "log_Fmort",
      value = log(pkfitfinal$rep$F)
    ),
    by = c("module_name", "fleet_name", "label", "time")
  ) |>
  dplyr::rows_update(
    tibble::tibble(
      module_name = "Fleet",
      fleet_name = c("survey2", "survey3", "survey6"),
      label = "log_q",
      value = c(parfinal$log_q2_mean, parfinal$log_q3_mean, parfinal$log_q6),
      estimation_type = "constant"
    ),
    by = c("module_name", "fleet_name", "label")
  )

# Put it all together, creating the FIMS model and making the TMB fcn
# Run the model without optimization to help ensure a viable model
test_fit <- default_parameters  |>
  FIMS::initialize_fims(data = data_4_model) |>
  FIMS::fit_fims(optimize = FALSE)
FIMS::get_obj(test_fit)$report()$nll_components |> length()
[1] 9
Code 
FIMS::get_data(data_4_model) |>
  dplyr::filter(name == "fleet1", type == "age_comp") |>
  print(n = 10)
# A tibble: 540 × 7
   type     name     age timing value unit  uncertainty
   <chr>    <chr>  <dbl>  <dbl> <dbl> <chr>       <dbl>
 1 age_comp fleet1     1   1970  -999 ""              1
 2 age_comp fleet1     2   1970  -999 ""              1
 3 age_comp fleet1     3   1970  -999 ""              1
 4 age_comp fleet1     4   1970  -999 ""              1
 5 age_comp fleet1     5   1970  -999 ""              1
 6 age_comp fleet1     6   1970  -999 ""              1
 7 age_comp fleet1     7   1970  -999 ""              1
 8 age_comp fleet1     8   1970  -999 ""              1
 9 age_comp fleet1     9   1970  -999 ""              1
10 age_comp fleet1    10   1970  -999 ""              1
# ℹ 530 more rows
Code 
# Run the  model with optimization
fit <- default_parameters  |>
  FIMS::initialize_fims(data = data_4_model) |>
  FIMS::fit_fims(optimize = TRUE)
✔ Starting optimization ...
ℹ Restarting optimizer 3 times to improve gradient.
ℹ Maximum gradient went from 0.00531 to 0.00048 after 3 steps.
✔ Finished optimization
✔ Finished sdreport
Warning: ! Large condition number detected in Hessian; the matrix may be near singular.
ℹ Condition number of Hessian ("1.518566e+05") exceeds threshold of 1e+05.
ℹ This suggests the model is weakly identified and results may be unreliable.
ℹ Consider simplifying the model, improving data quality, or fixing poorly
  informed parameters.
ℹ The 2 largest standard error values are for parameters:
ℹ 1. "expected_recruitment.53": "8.701804e+09"
ℹ 2. "numbers_at_age.530": "8.701804e+09"
ℹ Standard errors and MLEs may be unreliable.
ℹ FIMS model version: 0.9.2
ℹ Total run time was 30.80789 seconds
ℹ Number of parameters: fixed_effects=66, random_effects=53, and total=119
ℹ Maximum gradient= 0.00048
ℹ Negative log likelihood (NLL):
• Marginal NLL= 3302.808
• Total NLL= 3242.3832
ℹ Terminal SB= 633843606.87812
Code 
## report values for models
rep1 <- FIMS::get_report(test_fit) # FIMS initial values
FIMS::get_max_gradient(fit) # TODO: from Cole, can use TMBhelper::fit_tmb to get val to <1e-10
[1] 0.0004766463
Code 
rep2 <- FIMS::get_report(fit)
output <- get_estimates(fit)

Plotting results

All metrics match quite well between the two models.

Comparisons of biomass, catch, fishing mortality (F), survey indices of abundance, recruitment, and spawning biomass between the Template Model Builder (TMB) model and the optimized FIMS model.

All metrics match quite well between the two models.

Relative error in biomass, catch, fishing mortality (F), survey indices of abundance, recruitment, and spawning biomass between the Template Model Builder (TMB) model and the optimized FIMS model. Relative error is also shown for the pre-optimized FIMS models (FIMS init).

Older ages are overestimated in 2010 by all models.

Proportion at age (paa) by year (panel) for the fishing fleet for the initialized FIMS model (FIMS init), the optimized FIMS model (FIMS opt), the Template Model Builder model (TMB), and observed data (obs).

The best fits are in year 2001.

Proportion at age (paa) by year (panel) for the first survey for the initialized FIMS model (FIMS init), the optimized FIMS model (FIMS opt), the Template Model Builder model (TMB), and observed data (obs).

All models fit the same.

Fits to the indices of abundance (million metric tons) by survey (panel) for the initialized FIMS model (FIMS init), the optimized FIMS model (FIMS opt), the Template Model Builder model (TMB). Observed data with uncertainty measurements are plotted using points.

Comparison table

The likelihood components from the TMB model do not include constants and thus are not directly comparable. To be fixed later. Relative differences between the modified TMB model and FIMS implementation are given in the figure above.

What was your experience using FIMS? What could we do to improve usability?

To do

List any issues that you ran into or found

  • Output more derived quantities like selectivity, maturity, etc.
  • NLL components are not separated by fleet and need to be. So age comp NLL for fleets 1 and 2 need to be separate to make, e.g., the likelihood profile plot above.

What features are most important to add based on this case study?

  • More sophisticated control over selectivity so that ages 1 and 2 can be zeroed out for a double-logistic form, overriding the selectivity curve.
Code 
FIMS::clear()