The model presented in this case study was changed substantially from the operational assessment model that was presented to management in 2022 for Southern New England-Mid Atlantic Yellowtail Flounder (NEFSC, 2022) and should not be considered reflective of the Southern New England-Mid Atlantic Yellowtail Flounder stock. These results are intended to be a demonstration for comparing the results of a FIMS model to the results of an Age-Structured Assessment Program (ASAP) model.
To get the operational model to more closely match what is capable in a FIMS model, the following changes were made:
The following additional changes were made as well though they are not necessary:
Missing years of survey data were filled in using predictions from the simplified ASAP model. The ASAP model was first ran with the missing data treated as missing. Then the missing data was filled in with the expected value of the survey, both in aggregate and for catch at age in place of the missing data.
# read the ASAP rdat files
rdat <- dget(file.path(data_directory, "NEFSC_YT_SIMPLIFIED.RDAT")) # to be used in FIMS, lots of modifications from original
orig <- dget(file.path(data_directory, "NEFSC_YT_ORIGINAL.RDAT")) # where started before modifications for use in FIMS
# set up FIMS data objects
data_4_model <- FIMS::FIMSFrame(get_asap_data(rdat))
# define the dimensions
years <- seq(get_start_year(data_4_model), get_end_year(data_4_model))# clear memory
FIMS::clear()
parameters <- create_default_configurations(data = data_4_model) |>
create_default_parameters(data = data_4_model) |>
tidyr::unnest(data) |>
# ASAP assumes Fmult devs = 0
dplyr::rows_update(
tibble::tibble(
fleet_name = c("survey1", "survey2"),
label = "log_q",
value = log(rdat$initial.guesses$q.year1.init)
),
by = c("fleet_name", "label")
) |>
dplyr::rows_update(
tibble::tibble(
# ASAP can enter either R0 or SB0, use R0 in input file
label = c("log_rzero", "logit_steep"),
value = c(
log(rdat$initial.guesses$SR.inits$SR.scaler.init),
FIMS::logit(0.2, 1.0, rdat$initial.guesses$SR.inits$SR_steepness.init)
),
estimation_type = c("fixed_effects", "constant")
),
by = c("label")
) |>
# Maturity
# NOTE, to match FIMS for a protogynous stock,
# maturity values were obtained by fitting a logistic fcn to the age vector,
# mat.female*prop.female + mat.male*prop.male
# assuming an all female population
dplyr::rows_update(
tibble::tibble(
module_name = "Maturity",
label = c("inflection_point", "slope"),
value = c(1.8, 4)
),
by = c("module_name", "label")
) |>
dplyr::rows_update(
tibble::tibble(
module_name = "Population",
label = "log_init_naa",
age = get_ages(data_4_model),
value = log(rdat$N.age[1, ]),
estimation_type = "fixed_effects"
),
by = c("module_name", "label", "age")
) |>
dplyr::rows_update(
tibble::tibble(
module_name = "Population",
label = rep("log_M", get_n_years(data_4_model) * get_n_ages(data_4_model)),
time = rep(years, each = get_n_ages(data_4_model)),
age = rep(get_ages(data_4_model), get_n_years(data_4_model)),
# M is vector of age1 M X nyrs then age2 M X nyrs
value = log(as.numeric(t(rdat$M.age)))
),
by = c("module_name", "label", "age", "time"),
# TODO: fix this
unmatched = "ignore"
) |>
# Update log_sd for log_devs in the Recruitment module
dplyr::rows_update(
tibble::tibble(
module_name = "Recruitment",
label = "log_sd",
# ASAP takes one value per year but that value is often just repeated
value = mean(rdat$control.parms$recruit.cv)
),
by = c("module_name", "label")
)
# Run the model with optimization
fit <- parameters |>
initialize_fims(data = data_4_model) |>
fit_fims(optimize = TRUE)✔ Starting optimization ...
ℹ Restarting optimizer 3 times to improve gradient.
ℹ Maximum gradient went from 0.06677 to 0.00036 after 3 steps.
✔ Finished optimization
✔ Finished sdreport
ℹ FIMS model version: 0.9.2
ℹ Total run time was 17.87268 seconds
ℹ Number of parameters: fixed_effects=63, random_effects=46, and total=109
ℹ Maximum gradient= 0.00036
ℹ Negative log likelihood (NLL):
• Marginal NLL= 2521.5788
• Total NLL= 2463.04499
ℹ Terminal SB= 843.58881output <- get_estimates(fit)
# eventually change to allow multiple fishing fleets similar to multiple indices - only using 1 fishing fleet for now
# NOTE: FIMS assumes SB calculated at the start of the year, so need to adjust ASAP to do so as well for now, need to make timing of SSB calculation part of FIMS later
# NOTE: for now tricking FIMS into thinking age 0 is age 1, so need to adjust A50 for maturity because FIMS calculations use ages 0-5+ instead of 1-6
sdr <- get_sdreport(fit)
sdr_fixed <- summary(sdr, "fixed")
report <- get_report(fit)convert_output_asap <- function(report, fleet_names) {
years <- seq(report$parms$styr, report$parms$endyr)
catch <- data.frame(
module_name = "Fleet",
module_id = rep(seq(NROW(report$catch.pred)), each = length(years)),
label = "landings_expected",
year_i = rep(seq(NCOL(report$catch.pred)), NROW(report$catch.pred)),
estimate = c(report$catch.pred),
estimation_type = "derived_quantity"
)
index <- data.frame(
module_name = "Fleet",
module_id = rep(seq(report$index.pred) + NROW(report$catch.pred), each = length(years)),
label = "indices_expected",
year_i = unlist(report$index.year.counter),
estimate = unlist(report$index.pred),
estimation_type = "derived_quantity"
)
index_observed <- index |>
dplyr::mutate(
label = "indices_observed",
estimate = unlist(report$index.obs)
)
spawning_biomass <- data.frame(
module_name = "Population",
label = "spawning_biomass",
year_i = seq(years),
estimate = report$SSB * 0.5,
estimation_type = "derived_quantity"
)
expected_recruitment <- data.frame(
module_name = "Population",
label = "expected_recruitment",
year_i = seq(years),
estimate = as.numeric(report$N.age[, 1]),
estimation_type = "derived_quantity"
)
# I don't know the structure of F.report when more than one fleet is present
log_Fmort <- data.frame(
module_name = "Fleet",
module_id = 1,
label = "log_Fmort",
type = "vector",
year_i = seq(years),
estimate = log(report$F.report),
estimation_type = "fixed_effects"
)
log_init_naa <- data.frame(
module_name = "Population",
label = "log_init_naa",
type = "vector",
age_i = rep(
seq(NCOL(report$N.age)),
NROW(report$N.age)
),
year_i = rep(
seq(NROW(report$N.age)),
each = NCOL(report$N.age)
),
estimate = log(c(t(report$N.age))),
estimation_type = "fixed_effects"
)
dplyr::bind_rows(
catch, index, index_observed,
spawning_biomass, expected_recruitment, log_Fmort, log_init_naa
) |>
dplyr::left_join(
data.frame(
year = years,
year_i = seq(years)
),
by = "year_i"
) |>
dplyr::mutate(
uncertainty = NA_real_,
uncertainty_label = "se",
fleet = as.character(factor(module_id, labels = fleet_names)),
era = "time"
)
}
output_asap <- convert_output_asap(rdat, get_fleets(data_4_model))
output_fims <- output |>
dplyr::mutate(
uncertainty_label = "se",
year = year_i + FIMS::get_start_year(data_4_model) - 1,
estimate = estimated,
era = "time"
)
stockplotr::plot_spawning_biomass(
list(
"ASAP" = output_asap,
"FIMS" = output_fims
)
)
stockplotr::plot_recruitment(
list(
"ASAP" = output_asap,
"FIMS" = output_fims
),
facet = "model"
)Ignoring unknown labels:
• colour : "Model"
# Not showing this plot because the hard-coded one is better
# stockplotr::plot_indices(
# dat = list(
# "ASAP" = output_asap,
# "FIMS" = dplyr::mutate(output_fims,
# label = dplyr::if_else(label == "index_expected", "indices_expected", label),
# fleet = as.character(factor(module_id, labels = get_fleets(data_4_model)))
# ) |>
# dplyr::filter(!is.na(observed))
# )
# )
mycols <- c("FIMS" = "blue", "ASAP" = "red", "ASAP_orig" = "darkgreen")
for (i in seq(rdat$parms$nindices)) {
index_results <- data.frame(
survey = i,
year = years,
observed = rdat$index.obs[[i]],
FIMS = report$index_expected[[rdat$parms$nfleet+i]],
ASAP = rdat$index.pred[[i]]
)
if (i==1){
allinds_results <- index_results
}else{
allinds_results <- rbind(allinds_results, index_results)
}
}
#print(allinds_results)
comp_index <- ggplot2::ggplot(allinds_results, ggplot2::aes(x = year, y = observed)) +
ggplot2::geom_point() +
ggplot2::geom_line(ggplot2::aes(x = year, y = FIMS), color = "blue") +
ggplot2::geom_line(ggplot2::aes(x = year, y = ASAP), color = "red") +
ggplot2::facet_wrap(~survey, scales = "free_y", nrow = 2) +
ggplot2::xlab("Year") +
ggplot2::ylab("Index") +
ggplot2::ggtitle("Blue=FIMS, Red=ASAP") +
ggplot2::theme_bw()
print(comp_index)
catch_results <- data.frame(
observed = get_data(data_4_model) |>
dplyr::filter(type == "landings", name == "fishery") |>
dplyr::pull(value),
FIMS = report$landings_expected[[1]],
ASAP = as.numeric(rdat$catch.pred[1,])
)
#print(catch_results)
comp_catch <- ggplot2::ggplot(catch_results, ggplot2::aes(x = years, y = observed)) +
ggplot2::geom_point() +
ggplot2::xlab("Year") +
ggplot2::ylab("Catch (mt)") +
ggplot2::geom_line(ggplot2::aes(x = years, y = FIMS), color = "blue") +
ggplot2::geom_line(ggplot2::aes(x = years, y = ASAP), color = "red") +
ggplot2::ggtitle("Blue=FIMS, Red=ASAP") +
ggplot2::theme_bw()
print(comp_catch)
pop_results <- data.frame(
Year = c(years, max(years)+1, years, years, years, years, max(years)+1, years),
Metric = c(rep("Spawning biomass", 2*get_n_years(data_4_model)+1), rep("Fishing mortality", 2*get_n_years(data_4_model)), rep("Recruitment", 2*get_n_years(data_4_model)+1)),
Model = c(rep("FIMS", get_n_years(data_4_model)+1), rep("ASAP", get_n_years(data_4_model)), rep(c("FIMS", "ASAP"), each = get_n_years(data_4_model)),
rep("FIMS", get_n_years(data_4_model)+1), rep("ASAP", get_n_years(data_4_model))),
Value = c(
report$spawning_biomass[[1]],
rdat$SSB / 2,
output |>
dplyr::filter(label == "log_Fmort", module_id == 1) |>
dplyr::pull(estimated) |>
exp(),
rdat$F.report,
report$expected_recruitment[[1]],
as.numeric(rdat$N.age[,1])
)
)
#print(pop_results)
ggplot2::ggplot(dplyr::filter(pop_results, Year <=2019), ggplot2::aes(x=Year, y=Value, color=Model)) +
ggplot2::geom_line() +
ggplot2::facet_wrap(~Metric, ncol=1, scales = "free_y") +
ggplot2::theme_bw() +
ggplot2::scale_color_manual(values = mycols)
orig_years <- seq(orig$parms$styr, orig$parms$endyr)
orig_pop_results <- data.frame(
Year = rep(orig_years, 3),
Metric = rep(c("Spawning biomass", "Fishing mortality", "Recruitment"), each = length(orig_years)),
Model = "ASAP_orig",
Value = c(orig$SSB / 2, orig$F.report, as.numeric(orig$N.age[,1]))
)
pop_results_3 <- rbind(pop_results, orig_pop_results)
#print(pop_results_3)
comp_FRSSB3 <- ggplot2::ggplot(pop_results_3, ggplot2::aes(x=Year, y=Value, color=Model)) +
ggplot2::geom_line() +
ggplot2::facet_wrap(~Metric, ncol=1, scales = "free_y") +
ggplot2::theme_bw() +
ggplot2::scale_color_manual(values = mycols)
print(comp_FRSSB3)
FIMS_naa_results <- data.frame(
Year = rep(c(years, max(years)+1), each = get_n_ages(data_4_model)),
Age = rep(get_ages(data_4_model), get_n_years(data_4_model)+1),
Metric = "NAA",
Model = "FIMS",
Value = report$numbers_at_age[[1]]
)
ASAP_naa_results <- data.frame(
Year = rep(years, each = get_n_ages(data_4_model)),
Age = rep(get_ages(data_4_model), get_n_years(data_4_model)),
Metric = "NAA",
Model = "ASAP",
Value = as.numeric(t(rdat$N.age))
)
orig_naa_results <- data.frame(
Year = rep(orig_years, each = get_n_ages(data_4_model)),
Age = rep(get_ages(data_4_model), length(orig_years)),
Metric = "NAA",
Model = "ASAP_orig",
Value = as.numeric(t(orig$N.age))
)
naa_results <- rbind(FIMS_naa_results, ASAP_naa_results, orig_naa_results)
#print(naa_results)
comp_naa3 <- ggplot2::ggplot(dplyr::filter(naa_results, Year <= 2019), ggplot2::aes(x=Year, y=Value, color=Model)) +
ggplot2::geom_line() +
ggplot2::facet_wrap(~Age, ncol=1, scales = "free_y") +
ggplot2::ylab("NAA") +
ggplot2::theme_bw() +
ggplot2::scale_color_manual(values = mycols)
print(comp_naa3)
The likelihood components from FIMS and ASAP for the same data are shown in the table below. Note that the ASAP file had to turn on the use likelihood constants option to enable this comparison (this option should not be used when recruitment deviations are estimated).
nll_components <- report[["nll_components"]]
jnlltab <- data.frame(
Component=c("Total","Index","Age Comp", "Rec"),
FIMS = c(
report[["jnll"]],
sum(nll_components[seq(2, length(nll_components), by = 2)]),
sum(nll_components[seq(3, length(nll_components), by = 2)]),
nll_components[1]
),
ASAP = c(
rdat$like$lk.total,
(rdat$like$lk.catch.total + rdat$like$lk.index.fit.total),
(rdat$like$lk.catch.age.comp + rdat$like$lk.index.age.comp),
rdat$like$lk.Recruit.devs
)
)
print(jnlltab) Component FIMS ASAP
1 Total 2463.04499 2410.844
2 Index 1018.74625 1020.548
3 Age Comp 1357.77718 1390.297
4 Rec 86.52156 0.000