Code
common_name <- "Petrale sole"NWFSC Case Study Petrale Sole
The setup
- R version: R version 4.6.1 (2026-06-24)
- TMB version: 1.9.25
- FIMS commit: c32d707
- Stock name: Petrale sole
- Region: NWFSC
- Analyst: Ian G. Taylor
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 Petrale sole stock. These results are intended to be a demonstration and nothing more.
To get the operational model to more closely match a FIMS model a simpler SS3 model was developed using code in https://github.com/NOAA-FIMS/case-studies/blob/petrale-with-new-r4ss-converter/content/R/NWFSC-petrale_simplify_assessment.R
These changes included * Remove data * Remove male lengths and ages * Remove discard fractions and discard comps * Simplify selectivity * Remove length-based retention functions * Convert to age-based logistic (from asymptotic length-based double-normal with time blocks) * Use time-invariant population-level female mean weight at age from the SS3 weight-at-age file for all fleets * Use time-invariant maturity x fecundity at age from the SS3 weight-at-age file
Setting up the data
Code
# convert SS3 data into FIMS format using function in r4ss and convert to FIMSFrame object
data_4_model <- r4ss::ss3_data_to_fims(
ss3_dir = "https://raw.githubusercontent.com/pfmc-assessments/petrale/main/models/2023.a050.003_FIMS_case-study_wtatage/",
ss_new = TRUE,
fleets = c(-4, 1, 2, 4)
) |>
FIMS::FIMSFrame()
# read SS3 output files from modified petrale sole assessment on GitHub
rep <- r4ss::SS_output(
"https://raw.githubusercontent.com/pfmc-assessments/petrale/main/models/2023.a050.003_FIMS_case-study_wtatage",
printstats = FALSE,
verbose = FALSE,
covar = FALSE
)
# # commented out pending fix to https://github.com/nmfs-ost/stockplotr/pull/317
# # alternative output processing using stockplotr
# output_petrale <- stockplotr::convert_output(
# "https://raw.githubusercontent.com/pfmc-assessments/petrale/main/models/2023.a050.003_FIMS_case-study_wtatage/Report.sso",
# model = "ss3"
# )Run FIMS model
Code
# lots of code below copied PIFSC-opakapaka.qmd
years <- (data_4_model |> FIMS::get_start_year()):(data_4_model |>
FIMS::get_end_year())
ages <- data_4_model |> FIMS::get_ages()
# Create default parameters
default_parameters <- FIMS::setup_default_parameters(data = data_4_model)Empirical weight-at-age rows found. Growth defaults to "EWAA".Code
# Note: the last value of the initial numbers at age is the first
# recruitment deviation
recdev_years <- default_parameters |>
dplyr::filter(module_name == "Recruitment" & label == "log_devs") |>
dplyr::pull(timing)
recdevs <- rep$recruitpars |> # $recruitpars is a subset of $parameters with Yr as an additional column
dplyr::filter(Yr %in% recdev_years) |>
dplyr::pull(Value)
# initial numbers-at-age from SS3 output
init_naa_ss3 <- 1000 *
rep$natage |>
dplyr::filter(Time == min(years), Sex == 1) |> # numbers at age in FIMS start year, not SS3 model
dplyr::select(paste(ages)) |>
as.numeric() # TODO: figure out if *1000 is needed
# alternative initial numbers-at-age based on exponential decay using female M
init_naa <- exp(
rep$parameters |> dplyr::filter(grepl("R0", Label)) |> dplyr::pull(Value)
) *
1000 *
exp(
-(ages - 1) *
rep$parameters |>
dplyr::filter(grepl("NatM.*Fem", Label)) |>
dplyr::pull(Value)
)
# init_naa calculation for plus group is based on sum of infinite series,
# which is equivalent to dividing the initial number at the last age by M
init_naa[length(ages)] <- init_naa[length(ages)] /
rep$parameters |>
dplyr::filter(grepl("NatM.*Fem", Label)) |>
dplyr::pull(Value) # sum of infinite series
# ratio is close to 1, suggesting that approximation is reasonable
# range(init_naa / init_naa_ss3)
# # [1] 1.037574 1.228785
# function to get estimated parameter value from SS3 model
get_ss3_estimate <- function(string) {
rep$parameters |>
dplyr::filter(grepl(
pattern = string,
Label,
fixed = TRUE
)) |>
dplyr::pull(Value)
}
parameters <- default_parameters |>
# dplyr::rows_update(
# tibble::tibble(
# module_name = "Maturity",
# label = c("inflection_point", "slope"),
# value = c(7.0, 1.0), # maturity based on the derived age-based pattern from original length-based SS3 model, arbitrary guess of slope
# ),
# by = c("module_name", "label")
# ) |>
# dplyr::rows_update(
# tibble::tibble(
# module_name = "Selectivity",
# fleet = "North",
# label = c("inflection_point", "slope"),
# value = c(10, 2),
# estimation_type = "constant"
# ),
# by = c("module_name", "fleet", "label")
# ) |>
# dplyr::rows_update(
# tibble::tibble(
# module_name = "Selectivity",
# fleet = "South",
# label = c("inflection_point", "slope"),
# value = c(10, 2),
# estimation_type = "constant"
# ),
# by = c("module_name", "fleet", "label")
# ) |>
# dplyr::rows_update(
# tibble::tibble(
# module_name = "Selectivity",
# fleet = "WCGBTS",
# label = c("inflection_point", "slope"),
# value = c(6, 2),
# estimation_type = "constant"
# ),
# by = c("module_name", "fleet", "label")
# ) |>
# set initial values for stock-recruitment parameters based on SS3 output
dplyr::rows_update(
tibble::tibble(
module_name = "Recruitment",
label = c("log_rzero", "logit_steep", "log_sd"),
value = c(
log(1000) + # adding log(1000) to account for R0 units in 1000s in SS3
get_ss3_estimate("SR_LN(R0)"),
FIMS::logit(
0.2,
1.0,
get_ss3_estimate("SR_BH_steep")
),
get_ss3_estimate("SR_sigmaR")
)
),
by = c("module_name", "label")
) |>
# set initial values for recruitment deviations based on SS3 output,
# starting with the second year of the SS3 model since the first recruitment deviation
# corresponds to the initial numbers at age
dplyr::rows_update(
tibble::tibble(
module_name = "Recruitment",
label = "log_devs",
timing = recdev_years,
value = recdevs,
),
by = c("module_name", "label", "timing")
) |>
# set initial value for log_q based on SS3 output
dplyr::rows_update(
tibble::tibble(
fleet = "WCGBTS",
label = "log_q",
value = get_ss3_estimate("LnQ_base_WCGBTS(4)")
),
by = c("fleet", "label")
) |>
# set initial numbers at age based on SS3 output calculated above
dplyr::rows_update(
tibble::tibble(
module_name = "Population",
label = "log_init_naa",
age = FIMS::get_ages(data_4_model),
value = log(init_naa),
estimation_type = "fixed_effects"
),
by = c("module_name", "label", "age")
) |>
# set initial M based on SS3 output
dplyr::rows_update(
tibble::tibble(
module_name = "Population",
label = "log_M",
value = log(get_ss3_estimate("NatM_uniform_Fem_GP_1")),
estimation_type = "constant"
),
by = c("module_name", "label")
)
# confirm that all columns except "value" and "estimation_type" are identical between parameters and default_parameters
if (!identical(
parameters |>
dplyr::select(-value, -estimation_type),
default_parameters |>
dplyr::select(-value, -estimation_type)
)) {
cli::cli_alert_danger("Columns other than 'value' and 'estimation_type' are not identical between parameters and default_parameters")
}
# figure out if any rows have different estimate_type
cli::cli_alert_info("Rows with different estimation_type between default_parameters parameters:")ℹ Rows with different estimation_type between default_parameters parameters:Code
diff_estimation_type <- parameters |>
dplyr::filter(estimation_type != default_parameters$estimation_type) |>
dplyr::select(module_name, fleet, label, estimation_type)
if (nrow(diff_estimation_type) == 0) {
cli::cli_alert_success("No rows with different estimation_type.")
} else {
print(diff_estimation_type, n = Inf, width = Inf)
}✔ No rows with different estimation_type.Code
# parameters which have been updated from defaults
cli::cli_alert_info("Rows with updated initial values from default_parameters to parameters:")ℹ Rows with updated initial values from default_parameters to parameters:Code
updated_parameter_counts <- parameters |>
dplyr::filter(value != default_parameters$value) |>
dplyr::select(module_name, fleet, label, value) |>
dplyr::count(module_name, label)
if (nrow(updated_parameter_counts) == 0) {
cli::cli_alert_success("No rows with updated initial values.")
} else {
print(updated_parameter_counts, n = Inf, width = Inf)
}Code
# Run the model without optimization to help ensure a viable model
test_fit <- parameters |>
FIMS::initialize_fims(data = data_4_model) |>
FIMS::fit_fims(optimize = FALSE)
# Run the model with optimization
fit <- parameters |>
FIMS::initialize_fims(data = data_4_model) |>
FIMS::fit_fims(optimize = TRUE, get_sd = FALSE)✔ Starting optimization ...
ℹ Restarting optimizer 3 times to improve gradient.
ℹ Maximum gradient went from 0.02246 to 0.00283 after 3 steps.
✔ Finished optimization
ℹ FIMS model version: 0.11.0.9000
ℹ Total run time was 58.86712 minutes
ℹ Number of parameters: fixed_effects=321, random_effects=146, and total=467
ℹ Maximum gradient= 0.00283
ℹ Negative log likelihood (NLL):
• Marginal NLL= 14103.62215
• Total NLL= 13964.60611
ℹ Terminal SB= 10153.87401Code
fit_est <- fit # backup in case we want to set fit <- test_fitPlotting results
TODO: Add plotting code for the FIMS model results
Code
FIMS::clear()