Pacific Hake

The setup

Code

# Names of required packages
packages <- c("dplyr", "tidyr", "ggplot2", "TMB", "reshape2", "here", "remotes", "lubridate")

# Install packages not yet installed
installed_packages <- packages %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
  install.packages(packages[!installed_packages], repos = "http://cran.us.r-project.org")
}

remotes::install_github("NOAA-FIMS/FIMS")
remotes::install_github("r4ss/r4ss")

# Load packages
invisible(lapply(packages, library, character.only = TRUE))

library(FIMS)

R_version <- version$version.string
TMB_version <- packageDescription("TMB")$Version
FIMS_commit <- substr(packageDescription("FIMS")$GithubSHA1, 1, 7)

source(file.path("R", "utils.R"))

Code

theme_set(theme_bw())

R version: R version 4.4.2 (2024-10-31)
TMB version: 1.9.16
FIMS commit: eca8095
Stock name: Pacific Hake
Region: Joint Technical Committee
Analyst: Kelli F. Johnson

Simplifications to the original assessment

Weight-at-age data

Weight-at-age data was taken from the .ss_new file because this file contains data for all years in the model rather than just since 1975 based on the rules supplied in the input file to the Stock Synthesis model. No simplifications were made for this assessment regarding weight-at-age data because all matrices in the Stock Synthesis model are replicates of each other, e.g., weight at age for fleet two is the same as beginning and middle of the year weight at age in this model.

Maturity

Year-specific maturity ogives for Pacific Hake were translated to a single logistic curve. First, estimates of maturity at age for the maximum day of the year investigated were averaged across years. Subsequently, this ogive was plotted and several values for the inflection point and slope parameters of the logistic curve were used to find a logistic curve that matched the average curve as best as possible. This lead to a curve with an inflection of 2.2 and a slope of 1.8.

Code

maturity_data |>
  dplyr::filter(doy == max(doy)) |>
  dplyr::group_by(age) |>
  dplyr::summarize(empirical = mean(p_mature)) |>
  dplyr::full_join(
    execute_logistic(
      # TODO: un-hard code these age values
      ages,
      hake_maturity_slope,
      hake_maturity_inflection_point
    ) |>
      dplyr::rename(age = x, logit = value),
    by = "age"
  ) |>
  tidyr::pivot_longer(cols = c("empirical", "logit")) |>
  ggplot2::ggplot(ggplot2::aes(age, value, color = name)) +
  ggplot2::geom_line(linewidth = 1.5) +
  ggplot2::xlab("Age") +
  ggplot2::ylab("Probability of being mature") +
  ggplot2::scale_color_brewer(palette = "Dark2", name = "Method")

Selectivity

The logistic and double logistic selectivity curves are currently the only selectivity forms that are available in FIMS. A researcher from the Southwest Fisheries Science Center is working on adding age-specific selectivity, which will be available later this calendar year.

Age-1 index

The relative age-1 index was removed from the data because FIMS currently works on biomass and does not allow for inputs in numbers. This relative index could be included after an agreed upon method to translate numbers of age-1 fish to biomass is established. This could be as simple as multiplying annual numbers of fish by weight of age-1 fish from the weight-at-age estimates from the GLMM.

Data

All data for a FIMS model is contained in one data frame. Here, we call this object fims_data. The creation of this data set is assisted by the R function called get_ss3_data() written by Drs. Ian G. Taylor and Megumi C. Oshima.

Because data-weighting is not available in FIMS, the input sample size for the survey age-composition data was multiplied by 10. Additionally, because non-parametric selectivity is not yet available in FIMS, the terminal year of data was removed from this data set.

Code

inputs <- r4ss::SS_read(base_dir)
inputs[["wtatage"]] <- r4ss::SS_readwtatage(file.path(base_dir, "wtatage.ss_new"))

fims_data <- get_ss3_data(
  ss3_inputs = inputs,
  fleets = seq(inputs[["dat"]][["Nfleets"]]),
  ages = inputs[["dat"]][["agebin_vector"]]
) |>
  dplyr::filter(
    datestart <= as.Date(paste0(inputs[["dat"]][["endyr"]], "-01-01")),
    value != -999,
    # TODO: allow fleet3 after I figure out how to go from numbers to mt
    name != "fleet3"
  ) |>
  dplyr::mutate(
    uncertainty = ifelse(
      test = type == "age" & name == "fleet2",
      uncertainty * 10,
      uncertainty
    )
  ) |>
  dplyr::filter(
    datestart != "2024-01-01"
  ) |>
  FIMS::FIMSFrame()
n_ages_years <- get_n_ages(fims_data) * get_n_years(fims_data)

Code

# Printing the S4 class that is returned from FIMSFrame()
fims_data

tbl_df of class 'FIMSFrame'

with the following 'types': age, landings, weight-at-age, index

# A tibble: 6 × 10
  type  name     age length datestart  dateend    value unit  uncertainty  year
  <chr> <chr>  <dbl>  <dbl> <date>     <date>     <dbl> <chr>       <dbl> <dbl>
1 age   fleet1     1     NA 1966-01-01 1966-12-31  -999 ""             NA  1966
2 age   fleet1     2     NA 1966-01-01 1966-12-31  -999 ""             NA  1966
3 age   fleet1     3     NA 1966-01-01 1966-12-31  -999 ""             NA  1966
4 age   fleet1     4     NA 1966-01-01 1966-12-31  -999 ""             NA  1966
5 age   fleet1     5     NA 1966-01-01 1966-12-31  -999 ""             NA  1966
6 age   fleet1     6     NA 1966-01-01 1966-12-31  -999 ""             NA  1966
additional slots include the following:fleets:
[1] 1 2
n_years:
[1] 58
ages:
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15
n_ages:
[1] 15
lengths:
numeric(0)
n_lengths:
[1] 0
start_year:
[1] 1966
end_year:
[1] 2023

Code

# Summary of the data types available
FIMS::get_data(fims_data) |>
  dplyr::filter(value != -999) |>
  dplyr::group_by(type, name) |>
  dplyr::count()

# A tibble: 5 × 3
# Groups:   type, name [5]
  type          name       n
  <chr>         <chr>  <int>
1 age           fleet1   735
2 age           fleet2   225
3 index         fleet2    15
4 landings      fleet1    58
5 weight-at-age fleet1   870

Model

Parameterization

Currently, only the multinomial distribution is available for fitting composition data in FIMS. Although, a double-logistic curve was explored for the fishery, results below are for the logistic. Survey selectivity was fixed with an inflection point at 2 years and a slope of 1.

Natural mortality can be estimated in FIMS using age- and year-specific values but all values were fixed at the median estimate from Stock Synthesis for this analysis.

Code

fleet1 <- list(
  selectivity = list(form = "LogisticSelectivity"),
  # selectivity = list(form = "DoubleLogisticSelectivity"),
  data_distribution = c(
    Index = "DlnormDistribution",
    AgeComp = "DmultinomDistribution"
  )
)
fleet2 <- list(
  selectivity = list(form = "LogisticSelectivity"),
  data_distribution = c(
    Index = "DlnormDistribution",
    AgeComp = "DmultinomDistribution"
  )
)
fleet3 <- list(
  selectivity = list(form = "DoubleLogisticSelectivity"),
  data_distribution = c(Index = "DlnormDistribution")
)

# Create default parameters
default_parameters <- fims_data |>
  create_default_parameters(
    fleets = list(
      fleet1 = fleet1,
      fleet2 = fleet2
      # TODO: include age-1 index data
      # fleet3 = fleet3
    )
  ) |>
  update_parameters(
    modified_parameters = list(
      maturity = list(
        LogisticMaturity.inflection_point.value =
          hake_maturity_inflection_point,
        LogisticMaturity.inflection_point.estimated = FALSE,
        LogisticMaturity.slope.value = hake_maturity_slope,
        LogisticMaturity.slope.estimated = FALSE
      ),
      population = list(
        # TODO: estimate M
        Population.log_M.value = rep(log(0.235), length = n_ages_years)
      ),
      recruitment = list(
        DnormDistribution.log_sd.value = 1.4
      )
    )
  ) |>
  update_parameters(
    modified_parameters = list(
      fleet2 = list(
        # TODO: Estimate survey selectivity
        LogisticSelectivity.inflection_point.value = 2,
        LogisticSelectivity.inflection_point.estimated = FALSE,
        LogisticSelectivity.slope.value = 1,
        LogisticSelectivity.slope.estimated = FALSE
      )
    )
  )

Time-varying parameters

Parameters can vary with time, though the R interface is not currently fully developed to easily allow this. Hence, the code below.

Code

# Code to create time-varying parameters

default_parameters[["parameters"]][["fleet1"]][["LogisticSelectivity.inflection_point.value"]] <-
  rep(2.49, get_n_years(fims_data))
default_parameters[["parameters"]][["fleet1"]][["LogisticSelectivity.inflection_point.estimated"]] <-
  rep(TRUE, get_n_years(fims_data))

# Code for time-varying Double Logistic

# default_parameters[["parameters"]][["fleet1"]][["DoubleLogisticSelectivity.inflection_point_asc.value"]] <-
#   rep(2.49, get_n_years(fims_data))
# default_parameters[["parameters"]][["fleet1"]][["DoubleLogisticSelectivity.inflection_point_asc.estimated"]] <-
#   rep(TRUE, get_n_years(fims_data))
# default_parameters[["parameters"]][["fleet1"]][["DoubleLogisticSelectivity.inflection_point_desc.value"]] <-
#   rep(10, get_n_years(fims_data))
# default_parameters[["parameters"]][["fleet1"]][["DoubleLogisticSelectivity.inflection_point_desc.estimated"]] <-
#   rep(TRUE, get_n_years(fims_data))

Fit

Code

# Run the model without optimization to help ensure a viable model
test_fit <- default_parameters |>
  initialize_fims(data = fims_data) |>
  fit_fims(optimize = FALSE)
fit <- default_parameters |>
  initialize_fims(data = fims_data) |>
  fit_fims(optimize = TRUE)

✔ Starting optimization ...
ℹ Restarting optimizer 3 times to improve gradient.
ℹ Maximum gradient went from 0.30984 to 0.06327 after 3 steps.
✔ Finished optimization
✔ Finished sdreport
ℹ FIMS model version: 0.3.0.1
ℹ Total run time was 3.09748 seconds
ℹ Number of parameters: total=134, fixed_effects=134, and random_effects=0
ℹ Maximum gradient= 0.06327
ℹ Negative log likelihood (NLL):
• Marginal NLL= 425250.50284
• Total NLL= 425250.50284
ℹ Terminal SB=

Results

Code

get_data(fims_data) |>
  dplyr::filter(type == "landings") |>
  dplyr::mutate(
    year = dplyr::row_number() + get_start_year(fims_data) - 1,
    expected = get_report(fit)$exp_catch[[1]]
  ) |>
  ggplot2::ggplot(
    ggplot2::aes(year, value)
  ) +
  ggplot2::geom_point() +
  ggplot2::geom_line(
    ggplot2::aes(y = expected)
  ) +
  ggplot2::xlab("Year") +
  ggplot2::ylab("Landings (t)")

Code

get_estimates(fit) |>
  dplyr::filter(grepl("SB", name, ignore.case = TRUE)) |>
  dplyr::mutate(
    year = (get_start_year(fims_data)- 1):get_end_year(fims_data)
  ) |>
  ggplot2::ggplot(
    ggplot2::aes(x = year, y = value)
  ) +
  ggplot2::geom_line() +
  ggplot2::xlab("Year") +
  ggplot2::ylab("Spawning biomass (t)")

Code

# Index
get_estimates(fit) |>
  dplyr::filter(grepl("Index", name)) |>
  dplyr::mutate(
    year = rep(1:(dplyr::n() / 2), 2) + get_start_year(fims_data) - 1,
    fleet = rep(c("Fishery estimated", "Estimated"), each = dplyr::n() / 2)
  ) |>
  dplyr::filter(
    fleet == "Estimated"
  ) |>
  ggplot2::ggplot(
    ggplot2::aes(year, value, group = fleet, color = as.factor(fleet))
  ) +
  ggplot2::geom_line() +
  ggplot2::geom_point(
    data = inputs[["dat"]][["CPUE"]] |>
      dplyr::filter(index == 2) |>
      dplyr::mutate(year = year, fleet = "Data"),
    ggplot2::aes(y = obs, x = year)
  ) +
  ggplot2::scale_color_brewer(palette = "Dark2", name = "Method") +
  ggplot2::xlab("Year") +
  ggplot2::ylab("Index of abundance (t)")

Code

# Fishing mortality
get_estimates(fit) |>
  dplyr::filter(grepl("Fmort", name)) |>
  dplyr::mutate(
    year = rep(1:dplyr::n(), 1) + get_start_year(fims_data),
    fleet = rep(1, each = dplyr::n()),
    F = exp(value)
  ) |>
  ggplot2::ggplot(
    ggplot2::aes(year, F, group = fleet)
  ) +
  ggplot2::geom_line() +
  ggplot2::xlab("Year") +
  ggplot2::ylab("Fishing mortality")