South Atlantic Scamp

The setup

  • R version: R version 4.6.1 (2026-06-24)
  • TMB version: 1.9.25
  • FIMS commit: c32d707
  • Stock name: South Atlantic scamp grouper
  • Region: SEFSC
  • Analyst: Kyle Shertzer

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 South Atlantic scamp grouper stock. These results are intended to be a demonstration for comparing the results of a FIMS model to the results of a Beaufort Assessment Model (BAM).

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

The original stock assessment (SEDAR-68OA) was conducted using the Beaufort Assessment Model (BAM). That assessment included options not yet available in FIMS. So, the following simplifications or modifications were made to the BAM configuration to allow more direct comparisons with FIMS output.

  • Set spawning biomass calculations to occur on January 1, rather than the time of peak spawning.
  • Set abundance calculations for matching index data to occur on January 1, rather than mid-year sampling.
  • Dropped the estimation of the variation in size at age, as this parameter is not estimable without length-composition data.
  • Converted spawning biomass to be female only. Because scamp are a protogynous hermaphrodite, the original assessment accounted for males implicitly by computing spawning biomass as the sum of total (male + female) mature biomass. One way to match that accounting would be to assume that the female maturity vector equals that of both sexes combined, and then further assume that the population is 100% female. However, current FIMS does not allow deviation from a 50:50 sex ratio. Thus, for this example, FIMS and BAM use the maturity vector of both sexes combined, but apply that vector only to females and assume a 50:50 sex ratio when computing spawning biomass.
  • Female maturity at age modeled as a logistic function, rather than empirical.
  • Dropped the fishery-dependent growth curve and set mean size at age of the catch equal to the population growth curve.
  • Replaced Dirichlet-multinomial with standard multinomial distribution for fitting age-composition data.
  • Extended age compositions to include ages 1–20+, as modeled in the population. The original assessment modeled ages 1–20+ in the population but fit ages 1–15+ in the age compositions because of many zeros in the 16–20 age range.

The following additional changes were made as well though they are not necessary:

  • Dropped time blocks on selectivity for fleets.
  • Dropped all length-composition data.
  • Converted index data and predictions to occur in weight rather than numbers.
  • Extended estimates of recruitment deviations forward in time to the terminal year. Based on likelihood profiling, values from the last two years were not estimable, and thus were fixed in the original assessment.
  • Extended estimates of recruitment deviations backward in time to the initial year, 1969. The original assessment started estimating recruitment deviations in 1980, when age-composition data become available.
  • Converted all observed and predicted catch to weight (mt). The original assessment used units native to the data collection: commercial in pounds and recreational in numbers.
  • Dropped fishery-dependent indices.

Setting up the data

The data for the FIMS model was created from reading in a simplified version of the BAM model via an .Rdat file.

Code 
grDevices::graphics.off()

# get scamp data and output from simplified version of BAM
rdat_file_name <- file.path(data_directory, "scamp32o.rdat")
sca <- dget(rdat_file_name)

# Set dimensions
styr <- dplyr::first(sca$t.series$year)
endyr <- dplyr::last(sca$t.series$year)
years <- styr:endyr
n_years <- endyr - styr + 1 # the number of years which we have data for
ages <- sca$a.series$age # age vector.
n_ages <- length(ages) # the number of age groups.

# Prepare data; initialize all values with -999 (missing)
# fleet1 is commercial, fleet2 is recreational
fleet1_ac <- fleet2_ac <- survey_ac <- matrix(-999, nrow = n_years, ncol = n_ages)

# CV is arithmetic space and what is used to fit in BAM, which we convert to
# standard deviation and create_default_DlnormDistribution takes the log of it
# for you so we want sqrt(log(1.0 + fleet1_catch_cv^2))
fleet1_catch <- sca$t.series$L.COM.ob
fleet1_catch_cv <- sca$t.series$cv.L.COM
fleet1_catch_sd <- sqrt(log(1.0 + fleet1_catch_cv^2))
fleet2_catch <- sca$t.series$L.REC.ob
fleet2_catch_cv <- sca$t.series$cv.L.REC
fleet2_catch_sd <- sqrt(log(1.0 + fleet2_catch_cv^2))

survey_index <- sca$t.series$U.CVT.ob |>
  tidyr::replace_na(-999)
survey_index_cv <- sca$t.series$cv.U.CVT
survey_index_sd <- sqrt(log(1.0 + survey_index_cv^2)) |>
  tidyr::replace_na(1) |>
  base::replace(n_years - 1, 1) # manually replacing the 2020 CV as a missing value

# COMMENT: These multinomial entries are not whole numbers and many <1. This is not
# technically correct but it does not seem to make a difference based on testing.
fleet1_ac_n <- sca$t.series$acomp.COM.n |> 
  tidyr::replace_na(1)
fleet2_ac_n <- sca$t.series$acomp.REC.n |>
  tidyr::replace_na(1)
survey_ac_n <- sca$t.series$acomp.CVT.n |>
  tidyr::replace_na(1)

fleet1_ac[!is.na(sca$t.series$acomp.COM.n), ] <- sca$comp.mats$acomp.COM.ob 
fleet2_ac[!is.na(sca$t.series$acomp.REC.n), ] <- sca$comp.mats$acomp.REC.ob 
survey_ac[!is.na(sca$t.series$acomp.CVT.n), ] <- sca$comp.mats$acomp.CVT.ob 

## put data into fims friendly form
fleet1_catch_df <- data.frame(
  type = "catch",
  fleet = "fleet1",
  age = NA,
  timing = seq(styr, endyr),
  observed = as.numeric(fleet1_catch),
  unit = "mt",
  uncertainty = paste0(
    "~dlnorm(meanlog = log_catch_expected, sdlog =0.0001",
    # fleet1_catch_sd,
    ")"
  )
)

fleet2_catch_df <- data.frame(
  type = "catch",
  fleet = "fleet2",
  age = NA,
  timing = seq(styr, endyr),
  observed = as.numeric(fleet2_catch),
  unit = "mt",
  uncertainty = paste0(
    "~dlnorm(meanlog = log_catch_expected, sdlog =0.0001",
    # fleet2_catch_sd,
    ")"
  )
)

survey_index_df <- data.frame(
  type = "index",
  fleet = "survey1",
  age = NA,
  timing = seq(styr, endyr),
  observed = as.numeric(survey_index),
  unit = "",
  uncertainty = paste0(
    "~dlnorm(meanlog = log_catch_expected, sdlog =0.01",
    # survey_index_sd,
    ")"
  )
)

fleet1_ac_df <- data.frame(
  type = "age_comp",
  fleet = "fleet1",
  age = rep(seq(1, n_ages), n_years),
  timing = rep(seq(styr, endyr), each = n_ages),
  observed = as.numeric(t(fleet1_ac)),
  unit = "number",
  uncertainty = paste0(
    "~dmultinom(prob = agecomp_proportion, size = ",
    rep(fleet1_ac_n, each = n_ages),
    ")"
  )
)

fleet2_ac_df <- data.frame(
  type = "age_comp",
  fleet = "fleet2",
  age = rep(seq(1, n_ages), n_years),
  timing = rep(seq(styr, endyr), each = n_ages),
  observed = as.numeric(t(fleet2_ac)),
  unit = "number",
  uncertainty = paste0(
    "~dmultinom(prob = agecomp_proportion, size = ",
    rep(fleet2_ac_n, each = n_ages),
    ")"
  )
)

survey_ac_df <- data.frame(
  type = "age_comp",
  fleet = "survey1",
  age = rep(seq(1, n_ages), n_years),
  timing = rep(seq(styr, endyr), each = n_ages),
  observed = as.numeric(t(survey_ac)),
  unit = "number",
  uncertainty = paste0(
    "~dmultinom(prob = agecomp_proportion, size = ",
    rep(survey_ac_n, each = n_ages),
    ")"
  )
)

weight_at_age <- data.frame(
  type = "weight_at_age",
  fleet = NA_character_,
  age = seq(n_ages),
  timing = NA_integer_,
  observed = sca$a.series$wgt.mt,
  unit = "mt",
  uncertainty = NA_character_
)

catch <- rbind(fleet1_catch_df, fleet2_catch_df)
index <- survey_index_df
agecomps <- rbind(fleet1_ac_df, fleet2_ac_df, survey_ac_df)

data_4_model <- FIMS::FIMSFrame(
  rbind(catch, index, agecomps, weight_at_age)
)

Run FIMS model

Empirical weight-at-age rows found. Growth defaults to "EWAA".
✔ Starting optimization ...
ℹ Restarting optimizer 3 times to improve gradient.
ℹ Maximum gradient went from 1387.16397 to 21.91281 after 3 steps.
✔ Finished optimization
Warning: ✖ Optimization failed convergence checks.
ℹ Maximum absolute gradient ("2.191281e+01") is higher than 1. Model does not
  seem converged.
ℹ Skipping sdreport. Consider adjusting control parameters (eval.max, iter.max)
  or model structure.
ℹ Model fit returned for diagnostic purposes only. Results are not reliable.
ℹ FIMS model version: 0.11.0.9000
ℹ Total run time was 1.13396 seconds
ℹ Number of parameters: fixed_effects=52, random_effects=0, and total=52
ℹ Maximum gradient= 21.91281
ℹ Negative log likelihood (NLL):
• Marginal NLL= 5772269184.14019
• Total NLL= 5772269184.14019
ℹ Terminal SB= 593.06944
ℹ Identified fleet names:

ℹ com, rec, and cvt

ℹ Processing parms

ℹ Processing spr.brps

ℹ Processing parm.cons

ℹ Processing parm.tvec

ℹ Processing parm.tvec.cons

ℹ Processing parm.avec

ℹ Processing parm.avec.cons

ℹ Processing sdnr

ℹ Processing like

ℹ Processing N.age

ℹ Processing N.age.mdyr

ℹ Processing N.age.spawn

ℹ Processing B.age

ℹ Processing F.age

ℹ Processing Z.age

ℹ Processing L.age.pred.knum

ℹ Processing L.age.pred.mt

ℹ Processing size.age.fishery

ℹ Processing sel.age

ℹ Processing comp.mats

ℹ Processing t.series

ℹ Processing LD.pr.tseries

ℹ Processing a.series

ℹ Processing eq.series

ℹ Processing pr.series

ℹ Processing CLD.est.mats

ℹ Processing age.error

ℹ Processing proj.t.series

ℹ Processing proj.N.age

✔ Conversion finished!
Code 
######################################################################
# Plot results

# created by running colorspace::sequential_hcl(5, "Viridis")
cols <- c("#4B0055", "#00588B", "#009B95", "#53CC67", "#FDE333")

out.folder <- "figures"
dir.create(out.folder, showWarnings = FALSE)
plot.type <- "png"

get_selectivity_parameter <- function(fit, x) {
  good <- dplyr::filter(
    FIMS::get_estimates(fit),
    module_name == "Selectivity",
    label == x
  )
  out <- dplyr::pull(good, estimated)
  names(out) <- paste(
    dplyr::pull(good, label),
    dplyr::pull(good, parameter_id),
    sep = "_"
  )
  return(out)
}

selex.bam.fleet1 <- 1 / (1 + exp(-sca$parm.cons$selpar_slope_COM2[8] * (ages - sca$parm.cons$selpar_A50_COM2[8])))
selex.fims.fleet1 <- 1 / (1 + exp(-get_selectivity_parameter(fit, "slope")[1] * (ages - get_selectivity_parameter(fit, "inflection_point")[1])))
selex.bam.fleet2 <- 1 / (1 + exp(-sca$parm.cons$selpar_slope1_REC2[8] * (ages - sca$parm.cons$selpar_A50_REC2[8])))
selex.fims.fleet2 <- 1 / (1 + exp(-get_selectivity_parameter(fit, "slope")[2] * (ages - get_selectivity_parameter(fit, "inflection_point")[2])))
selex.bam.survey <- 1 / (1 + exp(-sca$parm.cons$selpar_slope1_CVT[8] * (ages - sca$parm.cons$selpar_A501_CVT[8])))
selex.fims.survey <- 1 / (1 + exp(-get_selectivity_parameter(fit, "slope")[3] * (ages - get_selectivity_parameter(fit, "inflection_point")[3])))

index_results <- data.frame(
  yr = styr:endyr,
  observed = c(
    rep(NA, length(styr:get_start_year(data_4_model)) - 1),
    FIMS::get_data(data_4_model) |>
      dplyr::filter(type == "index", fleet == "survey1") |>
      dplyr::pull(observed)
  ),
  fims.expected = c(rep(NA, length(styr:get_start_year(data_4_model)) - 1), report$index_expected[[3]]),
  bam.expected = sca$t.series$U.CVT.pr
) |>
  dplyr::filter(observed != -999.00)
fleet1_catch_results <- data.frame(
  yr = styr:endyr,
  observed = c(
    rep(NA, length(styr:get_start_year(data_4_model)) - 1),
    FIMS::get_data(data_4_model) |>
      dplyr::filter(type == "catch", fleet == "fleet1") |>
      dplyr::pull(observed)
  ),
  fims.expected = c(rep(NA, length(styr:get_start_year(data_4_model)) - 1), report$catch_expected[[1]]),
  bam.expected = sca$t.series$L.COM.pr
)
fleet2_catch_results <- data.frame(
  yr = styr:endyr,
  observed = c(
    rep(NA, length(styr:get_start_year(data_4_model)) - 1),
    FIMS::get_data(data_4_model) |>
      dplyr::filter(type == "catch", fleet == "fleet2") |>
      dplyr::pull(observed)
  ),
  fims.expected = c(rep(NA, length(styr:get_start_year(data_4_model)) - 1), report$catch_expected[[2]]),
  bam.expected = sca$t.series$L.REC.pr
)

fleet1_F_results <- data.frame(
  yr = styr:endyr,
  fims.F.fleet1 = c(
    rep(NA, length(styr:get_start_year(data_4_model)) - 1), 
    dplyr::filter(FIMS::get_estimates(fit), label == "log_Fmort", module_id == 1) |>
      dplyr::pull(estimated) |>
      exp()
  ),
  bam.F.fleet1 = log(sca$t.series$F.COM)
)
fleet2_F_results <- data.frame(
  yr = styr:endyr,
  fims.F.fleet2 = c(
    rep(NA, length(styr:get_start_year(data_4_model)) - 1),
    dplyr::filter(FIMS::get_estimates(fit), label == "log_Fmort", module_id == 2) |>
      dplyr::pull(estimated) |>
      exp()
  ),
  bam.F.fleet2 = log(sca$t.series$F.REC)
)

# Dropping the last (extra) year from FIMS output, assuming it is a projection yr (not an initialization yr)
fims.naa <- matrix(report$numbers_at_age[[1]], ncol = n_ages, byrow = TRUE)
fims.naa <- fims.naa[-NROW(fims.naa), ]
popn_results <- data.frame(
  yr = styr:endyr,
  fims.ssb = c(
    rep(NA, length(styr:get_start_year(data_4_model)) - 1),
    report$spawning_biomass[[1]][1:get_n_years(data_4_model)]
  ),
  fims.recruits = c(
    rep(NA, length(styr:get_start_year(data_4_model)) - 1),
    report$expected_recruitment[[1]][1:get_n_years(data_4_model)] / 1000
  ),
  fims.biomass = c(
    rep(NA, length(styr:get_start_year(data_4_model)) - 1),
    report$biomass[[1]][1:get_n_years(data_4_model)]
  ),
  fims.abundance = c(
    rep(NA, length(styr:get_start_year(data_4_model)) - 1),
    rowSums(fims.naa) / 1000
  ),
  bam.ssb = sca$t.series$SSB,
  bam.recruits = sca$t.series$recruits / 1000,
  bam.biomass = sca$t.series$B,
  bam.abundance = sca$t.series$N / 1000
)

yr.ind <- 1:n_years

yr.fleet1.ind <- yr.ind[fleet1_ac_n >= 0]
yr.fleet1.ac <- years[yr.fleet1.ind]
fims.fleet1.ncaa <- matrix(report$catch_numbers_at_age[[1]], ncol = n_ages, byrow = TRUE)
fims.fleet1.ncaa <- fims.fleet1.ncaa[yr.fleet1.ind, ]
fims.fleet1.caa <- fims.fleet1.ncaa / rowSums(fims.fleet1.ncaa)
bam.fleet1.caa <- sca$comp.mats$acomp.COM.pr
obs.fleet1.caa <- sca$comp.mats$acomp.COM.ob

yr.fleet2.ind <- yr.ind[fleet2_ac_n >= 0]
yr.fleet2.ac <- years[yr.fleet2.ind]
fims.fleet2.ncaa <- matrix(report$catch_numbers_at_age[[2]], ncol = n_ages, byrow = TRUE)
fims.fleet2.ncaa <- fims.fleet2.ncaa[yr.fleet2.ind, ]
fims.fleet2.caa <- fims.fleet2.ncaa / rowSums(fims.fleet2.ncaa)
bam.fleet2.caa <- sca$comp.mats$acomp.REC.pr
obs.fleet2.caa <- sca$comp.mats$acomp.REC.ob

yr.survey.ind <- yr.ind[survey_ac_n >= 0]
yr.survey.ac <- years[yr.survey.ind]
fims.survey.ncaa <- matrix(report$catch_numbers_at_age[[3]], ncol = n_ages, byrow = TRUE)
fims.survey.ncaa <- fims.survey.ncaa[yr.survey.ind, ]
fims.survey.caa <- fims.survey.ncaa / rowSums(fims.survey.ncaa)
bam.survey.caa <- sca$comp.mats$acomp.CVT.pr
obs.survey.caa <- sca$comp.mats$acomp.CVT.ob
######################################################################
png(filename = paste(out.folder, "/SEFSC_scamp_tseries_fits.", plot.type, sep = ""), width = 8, height = 10, units="in", res=72)
mat <- matrix(1:3, ncol = 1)
layout(mat = mat, widths = rep.int(1, ncol(mat)), heights = rep.int(1, nrow(mat)))
par(las = 1, mar = c(4.1, 4.25, 1.0, 0.5), cex = 1)

plot(index_results$yr, index_results$observed,
  ylim = c(0, max(index_results[, -1])),
  pch = 16, col = cols[1], ylab = "Index", xlab = ""
)
lines(index_results$yr, index_results$bam.expected, lwd = 3, col = cols[2])
lines(index_results$yr, index_results$fims.expected, lwd = 3, col = cols[4])
legend("topright",
  legend = c("observed", "BAM expected", "FIMS expected"),
  pch = c(16, -1, -1), lwd = c(-1, 3, 3), col = c(cols[1], cols[2], cols[4])
)

plot(fleet1_catch_results$yr, fleet1_catch_results$observed,
  ylim = c(0, max(fleet1_catch_results[, -1], na.rm = TRUE)),
  pch = 16, col = cols[1], ylab = "Fleet1 catch (mt)", xlab = ""
)
lines(fleet1_catch_results$yr, fleet1_catch_results$bam.expected, lwd = 3, col = cols[2])
lines(fleet1_catch_results$yr, fleet1_catch_results$fims.expected, lwd = 3, col = cols[4])

plot(fleet2_catch_results$yr, fleet2_catch_results$observed,
  ylim = c(0, max(fleet2_catch_results[, -1], na.rm = TRUE)),
  pch = 16, col = cols[1], ylab = "Fleet2 catch (mt)", xlab = ""
)
lines(fleet2_catch_results$yr, fleet2_catch_results$bam.expected, lwd = 3, col = cols[2])
lines(fleet2_catch_results$yr, fleet2_catch_results$fims.expected, lwd = 3, col = cols[4])

dev.off()
png 
  2
Code 
######################################################################
png(filename = paste(out.folder, "/SEFSC_scamp_tseries_F.", plot.type, sep = ""), width = 8, height = 8, units="in", res=72)
mat <- matrix(1:2, ncol = 1)
layout(mat = mat, widths = rep.int(1, ncol(mat)), heights = rep.int(1, nrow(mat)))
par(las = 1, mar = c(4.1, 4.25, 1.0, 0.5), cex = 1)

plot(fleet1_F_results$yr, fleet1_F_results$bam.F.fleet1,
  ylim = c(0, max(fleet1_F_results[, -1], na.rm = TRUE)),
  type = "l", lwd = 3, col = cols[2], ylab = "Fleet 1 F", xlab = ""
)
lines(fleet1_F_results$yr, fleet1_F_results$fims.F.fleet1, lwd = 3, col = cols[4])
legend("topleft",
  legend = c("BAM predicted", "FIMS predicted"),
  lwd = c(3, 3), col = c(cols[2], cols[4])
)
plot(fleet2_F_results$yr, fleet2_F_results$bam.F.fleet2,
  ylim = c(0, max(fleet2_F_results[, -1], na.rm = TRUE)),
  type = "l", lwd = 3, col = cols[2], ylab = "Fleet 2 F", xlab = ""
)
lines(fleet2_F_results$yr, fleet2_F_results$fims.F.fleet2, lwd = 3, col = cols[4])

dev.off()
png 
  2
Code 
######################################################################
png(filename = paste(out.folder, "/SEFSC_scamp_selex.", plot.type, sep = ""), width = 8, height = 10, units="in", res=72)
mat <- matrix(1:3, ncol = 1)
layout(mat = mat, widths = rep.int(1, ncol(mat)), heights = rep.int(1, nrow(mat)))
par(las = 1, mar = c(4.1, 4.25, 1.0, 0.5), cex = 1)

plot(ages, selex.bam.fleet1, lwd = 3, col = cols[2], type = "l", xlab = "", ylab = "Fleet1 selectivity")
lines(ages, selex.fims.fleet1, lwd = 3, col = cols[4])
legend("bottomright",
  legend = c("BAM predicted", "FIMS predicted"),
  lwd = c(3, 3), col = c(cols[2], cols[4])
)
plot(ages, selex.bam.fleet2, lwd = 3, col = cols[2], type = "l", xlab = "", ylab = "Fleet2 selectivity")
lines(ages, selex.fims.fleet2, lwd = 3, col = cols[4])
plot(ages, selex.bam.survey, lwd = 3, col = cols[2], type = "l", xlab = "Age", ylab = "Survey selectivity")
lines(ages, selex.fims.survey, lwd = 3, col = cols[4])

dev.off()
png 
  2
Code 
######################################################################

png(filename = paste(out.folder, "/SEFSC_scamp_tseries_popn.", plot.type, sep = ""), width = 8, height = 7, units="in", res=72)
mat <- matrix(1:4, ncol = 2)
layout(mat = mat, widths = rep.int(1, ncol(mat)), heights = rep.int(1, nrow(mat)))
par(las = 1, mar = c(4.1, 4.25, 1.0, 0.5), cex = 1)

plot(popn_results$yr, popn_results$bam.ssb,
  ylim = c(0, max(popn_results[, c(2, 6)], na.rm = TRUE)),
  type = "l", lwd = 3, col = cols[2], ylab = "Spawning biomass (mt)", xlab = ""
)
lines(popn_results$yr, popn_results$fims.ssb, lwd = 3, col = cols[4])
legend("topleft",
  legend = c("BAM predicted", "FIMS predicted"),
  lwd = c(3, 3), col = c(cols[2], cols[4])
)

plot(popn_results$yr, popn_results$bam.biomass,
  ylim = c(0, max(popn_results[, c(4, 8)], na.rm = TRUE)),
  type = "l", lwd = 3, col = cols[2], ylab = "Biomass (mt)", xlab = ""
)
lines(popn_results$yr, popn_results$fims.biomass, lwd = 3, col = cols[4])

plot(popn_results$yr, popn_results$bam.recruits,
  ylim = c(0, max(popn_results[, c(3, 7)], na.rm = TRUE)),
  type = "l", lwd = 3, col = cols[2], ylab = "Recruits (1000s)", xlab = ""
)
lines(popn_results$yr, popn_results$fims.recruits, lwd = 3, col = cols[4])

plot(popn_results$yr, popn_results$bam.abundance,
  ylim = c(0, max(popn_results[, c(5, 9)], na.rm = TRUE)),
  type = "l", lwd = 3, col = cols[2], ylab = "Abundance (1000s)", xlab = ""
)
lines(popn_results$yr, popn_results$fims.abundance, lwd = 3, col = cols[4])

dev.off()
png 
  2
Code 
######################################################################
png(filename = paste(out.folder, "/SEFSC_scamp_caa_fleet1.", plot.type, sep = ""), width = 8, height = 11, units="in", res=72)
mat <- matrix(1:18, ncol = 3)
layout(mat = mat, widths = rep.int(1, ncol(mat)), heights = rep.int(1, nrow(mat)))
par(las = 1, mar = c(2.2, 2.7, 0.5, 0.5), cex = 0.75)

for (i in 1:nrow(obs.fleet1.caa))
{
  plot(1:n_ages, obs.fleet1.caa[i, ], col = cols[1], xlab = "", ylab = "", pch = 16)
  lines(1:n_ages, bam.fleet1.caa[i, ], lwd = 3, col = cols[2])
  lines(1:n_ages, fims.fleet1.caa[i, ], lwd = 3, col = cols[4])
  if (i > 1) legend("topright", legend = yr.fleet1.ac[i], cex = 1, bty = "n")
  if (i == 1) {
    legend("topright",
      legend = c("Fleet1 Age Comps", "observed", "BAM expected", "FIMS expected"),
      pch = c(-1, 16, -1, -1), lwd = c(-1, -1, 3, 3), col = c(cols[1], cols[1], cols[2], cols[4]), cex = 0.7
    )
    legend("right", legend = yr.fleet1.ac[i], cex = 1, bty = "n")
  }
}

dev.off()
png 
  2
Code 
######################################################################
png(filename = paste(out.folder, "/SEFSC_scamp_caa_fleet2.", plot.type, sep = ""), width = 8, height = 11, units="in", res=72)

mat <- matrix(1:28, ncol = 4)
layout(mat = mat, widths = rep.int(1, ncol(mat)), heights = rep.int(1, nrow(mat)))
par(las = 1, mar = c(2.2, 2.7, 0.5, 0.5), cex = 0.75)

for (i in 1:nrow(obs.fleet2.caa))
{
  plot(1:n_ages, obs.fleet2.caa[i, ], col = cols[1], xlab = "", ylab = "", pch = 16)
  lines(1:n_ages, bam.fleet2.caa[i, ], lwd = 3, col = cols[2])
  lines(1:n_ages, fims.fleet2.caa[i, ], lwd = 3, col = cols[4])
  if (i > 1) legend("topright", legend = yr.fleet2.ac[i], cex = 1, bty = "n")
  if (i == 1) {
    legend("topright",
      legend = c("Fleet2 Age Comps", "observed", "BAM expected", "FIMS expected"),
      pch = c(-1, 16, -1, -1), lwd = c(-1, -1, 3, 3), col = c(cols[1], cols[1], cols[2], cols[4]), cex = 0.5
    )
    legend("right", legend = yr.fleet2.ac[i], cex = 1, bty = "n")
  }
}

dev.off()
png 
  2
Code 
######################################################################
png(filename = paste(out.folder, "/SEFSC_scamp_caa_survey.", plot.type, sep = ""), width = 8, height = 11, units="in", res=72)

mat <- matrix(1:30, ncol = 5)
layout(mat = mat, widths = rep.int(1, ncol(mat)), heights = rep.int(1, nrow(mat)))
par(las = 1, mar = c(2.2, 2.7, 0.5, 0.5), cex = 0.75)

for (i in 1:nrow(obs.survey.caa))
{
  plot(1:n_ages, obs.survey.caa[i, ], col = cols[1], xlab = "", ylab = "", pch = 16)
  lines(1:n_ages, bam.survey.caa[i, ], lwd = 3, col = cols[2])
  lines(1:n_ages, fims.survey.caa[i, ], lwd = 3, col = cols[4])
  if (i > 1) legend("topright", legend = yr.survey.ac[i], cex = 1, bty = "n")
  if (i == 1) {
    legend("topright",
      legend = c("Survey Age Comps", "observed", "BAM expected", "FIMS expected"),
      pch = c(-1, 16, -1, -1), lwd = c(-1, -1, 3, 3), col = c(cols[1], cols[1], cols[2], cols[4]), cex = 0.5
    )
    legend("right", legend = yr.survey.ac[i], cex = 1, bty = "n")
  }
}

dev.off()
png 
  2
Code 
FIMS::clear()

Comparison figures

Fits-tseries

Fits-CAA-fleet1

Fits-CAA-fleet2

Fits-CAA-survey

Predicted-selex

Predicted-tseries-F

Predicted-tseries-population

Figures with {stockplotr}

Code 
stockplotr::plot_biomass(
  list(
    "bam" = dplyr::filter(output_bam, module_name != "B.age") |>
      dplyr::mutate(module_name = ifelse(module_name == "t.series", "Population", module_name)),
    "fims" = output_fims
  )
)
! Multiple module names found in data. 
ℹ Environment not interactive. Selecting Population.
! biomass_msy not found for fims

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

  • Discards. I chose this assessment as a case study because it did not explicitly model discards, as I knew that capability was not yet available in FIMS. However, most assessments in the Southeast model dead discards.

Acknowledgments

Thanks to Nathan Vaughan and Ian Taylor, who helped troubleshoot earlier versions of the code.