Introduction

Reproducing the single stock advice (SSA) is a quality control step done by each WGMIXFISH case study before proceeding with mixed fishery modelling scenarios. Mixed fishery modelling is conducted with a different operating model than the SSA (e.g. FCube and FLBEIA), which is deterministic and thus will deviate from the stochastic forecasts used by many assessment forecasts.

Currently, assessors are asked to provide assessment summaries in the form of FLStock objects. These provide a starting point for setting up mixed fishery forecasts, but they are also use to see if we are able to approximate the SSA using a deterministic forecast in FLR. Specifically, we use information from the headline advice of each stock for forecast years. This includes assumptions used for the intermediate year (e.g. status quo F, TAC uptake) and advice year fishing targets (e.g. Fmsy). Additionally, assumptions for future biological parameters must also be considered (e.g. mean weights, maturity, recruitment).

During this process, we hope to identify and understand any deviations from SSA that might exist. In the past, larger deviations have been attributable to both simple reporting errors as well as inconsistencies between specified assumptions and those of the forecast.

The aim of this document is to demonstrate the procedure employed by WGMIXFISH so that assessors can help catch issues at an earlier stage in the advice calendar year, potentially reducing the frequency of re-openings.

In the past, assessors have been asked to deliver two FLStock objects to WGMIXFISH: 1. containing input data, and 2. containing all assessment estimated data. To improve consistency between SSA and WGMIXFISH forecasts, we would now ask that the FLStock containing estimated data also include the forecast years. This will ensure even greater consistency in the mixed fishery forecasts regarding biological parameter assumptions and allow for further quality control checks.

The following section outlines the basic procedure in setting up an FLR-based deterministic forecast. After that, an additional section highlights typical problems encountered, as relating to mis-specified forecast assumptions. Finally, a specialized approach is presented for Nephrops, whose advice is based on harvest rates and fixed stock dynamics, requiring a slightly altered methodology to forecast the stock.

The procedure

1. Start with an FLStock object with all assessment-estimated values in historical years
2. Extend FLStock into the short-term forecast years using FLasher::stf, and create a and, if necessary, manually adjust forecast year slot values to match assessment as closely as possible
3. Create FLasher::fwdControl object to define the targets and limits of the short-term forecast years
4. Run deterministic forecast with FLasher::fwd
5. Compare resulting stock parameters with those reported by the assessment

To demonstrate the procedure, the following example uses the assessment oand forecast information from North Sea saithe (pok.27.3a46) as assessed in 2024.

data("wgnsskStocks")

stkInp <- wgnsskStocks[["POK_input"]]
stkEst <- wgnsskStocks[["POK"]]

L <- FLStocks(list(input = stkInp, estimated = stkEst))
plot(L) + 
  aes(linetype = stock) +
  scale_color_manual(values = c(8,1)) + 
  scale_linetype_manual(values = c(1,2)) +
  theme_bw()

# forecast year definitions
yrAssess <- 2023 # final year of assessment data
yrNow <- 2024 # intermediate year
yrTAC <- 2025 # advice year
yrTACp1 <- 2026 # advice year +1 (needed to get SSB at end of yrTAC)

# extend FLStock object
stkProj <- stf(object = stkEst, nyears = 3, wts.nyears = 3, 
  fbar.nyears = 3, f.rescale = TRUE, disc.nyears = 3)

# stock-recruitment model (manual input within a geometric mean model)
srPar <- FLPar(c(94047, 94047, 94047), 
  dimnames = list(params="a", year = c(yrNow, yrTAC, yrTACp1), iter = 1))
srMod <- FLSR(model = "geomean", params = srPar)

# View the extended FLStock
df <- as.data.frame(stkProj)
df <- subset(df, slot %in% c("landings.wt", "discards.wt", "catch.wt", "m", "mat", "harvest") & year > (yrAssess-20))
df$forecast <- df$year %in% c(yrNow, yrTAC, yrTACp1)

ggplot(df) + aes(x = year, y = data, group = age, color = forecast) +
  facet_wrap(~slot, scales = "free_y") +
  geom_line(show.legend = F) + 
  scale_color_manual(values = c(8,1)) + 
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) + 
  theme_bw()

# stf control (Fsq, followed by 2 years at Fmsy)
Fmsy <- 0.316
ctrl <- fwdControl( 
  data.frame(
    year = c(yrNow, yrTAC, yrTACp1),
    value = c(1, Fmsy, Fmsy),
    quant = c("f"),
    relYear = c(yrAssess, NA, NA)                               
  )
)

# projection
stkProj <- fwd(object = stkProj, control = ctrl, sr = srMod)

# plot
L <- FLStocks(list(assessment = stkEst, forecast = stkProj[,ac(yrAssess:yrTACp1)]))
plot(L) + 
  # aes(linetype = stock) +
  scale_color_manual(values = c(8,1)) + 
  # scale_linetype_manual(values = c(2,1)) +
  theme_bw()

TACNow <- 73815
Fmsy <- 0.316

ctrlAlt <- fwdControl( 
  data.frame(
    year = c(yrNow, yrTAC, yrTACp1),
    value = c(TACNow, Fmsy, Fmsy),
    quant = c("catch", "f", "f")
  )
)

stkProjAlt <- fwd(object = stkProj, control = ctrlAlt, sr = srMod)

df <- data.frame(year = yrAssess:yrTACp1, 
  catch = c(catch(stkProjAlt)[, ac(yrAssess:yrTACp1)]),
  fbar = c(fbar(stkProjAlt)[, ac(yrAssess:yrTACp1)])
)

kable(df, digits = 3)

# Reported output from single stock headline advice
stfRef <- data.frame(
  model = "NSSK",
  year = 2023:2026,
  catch = c(65865, 78251, 79071, NA),
  landings = c(62691, 74968, 75880, NA),
  fbar = c(0.32, 0.32, 0.316, NA),
  ssb = c(161756, 185632, 195899, 197298)
)

stfDet <- data.frame(
  model = "FLR",
  year = ac(yrAssess:yrTACp1),
  catch = c(catch(stkProj[,ac(yrAssess:yrTACp1)])),
  landings = c(landings(stkProj[,ac(yrAssess:yrTACp1)])),
  fbar = c(fbar(stkProj[,ac(yrAssess:yrTACp1)])),
  ssb = c(ssb(stkProj[,ac(yrAssess:yrTACp1)]))
)  

df <- merge(stfRef, stfDet, all = T)
df <- pivot_longer(df, cols = c(catch, landings, fbar, ssb), 
  names_to = "variable", values_to = "value")
df <- df |>
  filter(
    (variable %in% c("catch", "landings", "fbar") & year <= yrTAC) |
    (variable %in% c("ssb") & year <= yrTACp1))
    
df2 <- pivot_wider(df, names_from = model, values_from = value)
df2$percErr <- round((df2$FLR - df2$NSSK)/df2$NSSK * 100, 1)

ggplot(df) + aes(x = year, y = value, group = model, color = model, shape = model) +
  facet_wrap(~variable, scales = "free_y") +
  geom_line() +
  geom_point(size = 3, stroke = 1) +
  scale_shape_discrete(solid = F) +
  coord_cartesian(ylim = c(0, NA)) +
  theme_bw()

kable(df2, digits = 3)

ggplot(df2) + aes(x = year, y = percErr) + 
  facet_wrap(~variable) +
  geom_col(fill = 4, color = 4) +
  geom_hline(yintercept = 0, linetype = 1) +
  geom_hline(yintercept = c(-10,10), linetype = 3) + 
  theme_bw()

Typical problems

Where deviations are larger than the +/-10% threshold, one might start investigation as to whether future slot values have been correctly defined. State-space models may model future conditions in a more sophisticated way than a simple average of the recent past, such as time-varying processes (e.g. random walk or auto-regressive), which may be used for a host of values (e.g. individual weights, maturity ogive, natural mortality, selectivity, etc.). These may need to be input manually into the forecast years of the slots.

Note that while time-varying selectivity may be largely reproduced in a deterministic forecast, the inclusion of this process in the subsequent mixed fishery model is currently impossible due to the splitting of catches across fleets and métiers. Nevertheless, it’s inclusion in the reproduce the advice quality control procedure is encouraged to successfully identify remaining issues.

The following example comes from North Sea haddock (had.27.46a20) as assessed in 2024, where manual input of forecast values for weight-at-age slots was required.

First, the procedure is conducted without these adjustments in order to illustrate the poor reproduction of advice.

# load haddock
stkEst <- wgnsskStocks[["HAD"]]

# trim to only historical years
minyear <- range(stkEst)["minyear"]
stkEst <- stkEst[,ac(minyear:yrAssess)]

# simple extension of FLStock
stkProj <- stf(object = stkEst, nyears = 3, wts.nyears = 3, 
  fbar.nyears = 3, f.rescale = TRUE, disc.nyears = 3)

# control
Fmsy <- 0.174
ctrl <- fwdControl( 
  data.frame(
    year = c(yrNow, yrTAC, yrTACp1),
    value = c(1, Fmsy, Fmsy),
    quant = c("f"),
    relYear = c(yrAssess, NA, NA)                               
  )
)

# stock-recruitment model
srPar <- FLPar(c(1788860, 1788860, 1788860), 
  dimnames = list(params="a", year = c(yrNow, yrTAC, yrTACp1), iter = 1))
srMod <- FLSR(model = "geomean", params = srPar)

# projection
stkProj1 <- fwd(object = stkProj, control = ctrl, sr = srMod)

# plot
L <- FLStocks(list(assessment = stkEst, forecast = stkProj1[,ac(yrAssess:yrTACp1)]))
plot(L) + 
  scale_color_manual(values = c(8,1)) + 
  theme_bw()

# Reported output from single stock headline advice
stfRef <- data.frame(
  model = "NSSK",
  year = 2023:2026,
  catch = c(60979, 66309, 112435, NA),
  landings = c(37702, 53277, 93661, NA),
  fbar = c(0.084, 0.084, 0.174, NA),
  ssb = c(580727, 606378, 535682, 420140)
)

stfDet1 <- data.frame(
  model = "FLR",
  year = ac(yrAssess:yrTACp1),
  catch = c(catch(stkProj1[,ac(yrAssess:yrTACp1)])),
  landings = c(landings(stkProj1[,ac(yrAssess:yrTACp1)])),
  fbar = c(fbar(stkProj1[,ac(yrAssess:yrTACp1)])),
  ssb = c(ssb(stkProj1[,ac(yrAssess:yrTACp1)]))
)  

df1 <- merge(stfRef, stfDet1, all = T)
df1 <- pivot_longer(df1, cols = c(catch, landings, fbar, ssb), 
  names_to = "variable", values_to = "value")
df1 <- df1 |>
  filter(
    (variable %in% c("catch", "landings", "fbar") & year <= yrTAC) |
    (variable %in% c("ssb") & year <= yrTACp1))
df1 <- pivot_wider(df1, names_from = model, values_from = value)
df1$percErr <- round((df1$FLR - df1$NSSK)/df1$NSSK * 100, 1)

ggplot(df1) + aes(x = year, y = percErr) + 
  facet_wrap(~variable) +
  geom_col(fill = 4, color = 4) +
  geom_hline(yintercept = 0, linetype = 1) +
  geom_hline(yintercept = c(-10,10), linetype = 3) + 
  theme_bw()

Deviations above the threshold are seen for all summary parameters except F, which is matched more or less perfectly due to it’s input as a controlling parameter (in fwdControl object).

Now we adjust the weight-at-age values assumed by the SSA, and compare again.

stkProj2 <- stkProj

# replace individual weights
data("haddockWts")
catch.wt(stkProj2)[, ac(yrNow:yrTACp1)][] <- t(haddockWts$catch) # catch weight-at-age
landings.wt(stkProj2)[, ac(yrNow:yrTACp1)][] <- t(haddockWts$landings) # landings weight-at-age
discards.wt(stkProj2)[, ac(yrNow:yrTACp1)][] <- t(haddockWts$discards) # discards weight-at-age
stock.wt(stkProj2)[, ac(yrNow:yrTACp1)][] <- t(haddockWts$stock) # catcstockh weight-at-age

stkProj2 <- fwd(object = stkProj2, control = ctrl, sr = srMod)

stfDet2 <- data.frame(
  model = "FLR",
  year = ac(yrAssess:yrTACp1),
  catch = c(catch(stkProj2[,ac(yrAssess:yrTACp1)])),
  landings = c(landings(stkProj2[,ac(yrAssess:yrTACp1)])),
  fbar = c(fbar(stkProj2[,ac(yrAssess:yrTACp1)])),
  ssb = c(ssb(stkProj2[,ac(yrAssess:yrTACp1)]))
)

df2 <- merge(stfRef, stfDet2, all = T)
df2 <- pivot_longer(df2, cols = c(catch, landings, fbar, ssb), 
  names_to = "variable", values_to = "value")
df2 <- df2 |>
  filter(
    (variable %in% c("catch", "landings", "fbar") & year <= yrTAC) |
    (variable %in% c("ssb") & year <= yrTACp1))
df2 <- pivot_wider(df2, names_from = model, values_from = value)
df2$percErr <- round((df2$FLR - df2$NSSK)/df2$NSSK * 100, 1)
df2$stf <- "adjusted"

df <- df1
df$stf <- "original"
df <- rbind(df, df2)
df$stf <- factor(df$stf, levels = c("original", "adjusted"))

ggplot(df) + aes(x = year, y = percErr, fill = stf, colour = stf) + 
  facet_wrap(~variable) +
  geom_col(position = position_dodge(), width = 0.75) +
  scale_fill_manual(values = c(4,5)) + 
  scale_color_manual(values = c(4,5)) + 
  geom_hline(yintercept = 0, linetype = 1) +
  geom_hline(yintercept = c(-10,10), linetype = 3) + 
  theme_bw()

Residual differences are likely attributable to the use of the intermediate year IBTS survey index in the forecast, which influences the starting numbers in 2024. This information is not currently used in the FLR forecast, nor in the the mixed fishery model.

Special case: fixed dynamics and harvest rate advice (Nephrops)

Nephrops stocks are generally treated as a fixed dynamic stock, except in a few rare cases where surplus production models have been fit (e.g. SPiCT). Some Nephrops functional units (FUs) also contain abundance estimates via an underwater TV-survey, while others rely on catch trends to provide advice.

WGMIXFISH includes all Nephrops FUs in mixed fishery scenarios, and thus requires an FLStock for each. Simply translating the Nephrops assessment into an FLStock provides a good quality control check for the accounting done during the assessment.

In cases where fixed dynamics are assumed and abundance estimates are available, harvest rates (hr, in percent) can be used in place of fishing mortality (F). This is specified in the @harvest slot through the units (units(harvest(stkEst)) <- "hr"). A short-term forecast can test whether an advised harvest rate results in the same TAC as reported.

As advised harvest rates for Nephrops are based on dead removals (landings + dead discards), particular care should be given to only include only dead discards in the discards.n slot.

Several functions exist to calculate total weight based on the product of numbers and individual weights (computeCatch, computeLandings, computeDiscards), providing a further test of the assessment calculations. Harvest rates can also be recalculated using computeHarvest.

The FLasher::fwd function does not currently handle harvest rates, and thus quality control tests for the advice year are calculated manually. Total removals are calculated based on the advised harvest rate multiplied by the stock size (catch.n = stock.n * hr).

The following example for FU 6 (nep.fu.6) shows a few tests that ensure the quality control of the FLStock. A TV-survey estimate is used to provide the assumed stock size in the advice year. In this case, the discards slot required adjustment to reflect dead discards only. The resulting comparison showed good agreement with reported values, with minor differences likely due to rounded values in the advice document.

stkEst <- wgnsskStocks[["NEP6"]]
stkEst <- window(stkEst, start = 2001) # trim data, as earlier years only contain landings
stkEst@landings <- computeLandings(stkEst)

# replace discards totals with dead discards only
stkEst@discards[,ac(2001:2023)] <- c(2034, 676, 608, 523, 608, 893, 367, 141, 
  392, 171, 209, 293, 383, 169, 162, 231, 170, 166, 385, 264, 356, 289, 281)
stkEst@discards.n <- stkEst@discards / stkEst@discards.wt

# re-compute catch
stkEst@catch.n <- stkEst@landings.n + stkEst@discards.n
stkEst@catch <- stkEst@landings + stkEst@discards
stkEst@catch.wt[] <- stkEst@catch / stkEst@catch.n

# re-compute harvest rate (hr)
stkEst@harvest <- computeHarvest(stkEst)

# double-check dead discard rate
# [email protected] / [email protected]

# for clarity, units can be redefined
units(stkEst)[c("catch", "landings", "discards")] <- "tonnes"
units(stkEst)[c("catch.n", "landings.n", "discards.n")] <- "millions"
units(stkEst)[c("catch.wt", "landings.wt", "discards.wt")] <- "g"

# extend FLStock
stkProj <- stf(object = stkEst, nyears = 3, wts.nyears = 3, 
  fbar.nyears = 3, f.rescale = TRUE, disc.nyears = 3)

# view mean weights
# [email protected][,ac(yrNow:yrTACp1)] # correct
# [email protected][,ac(yrNow:yrTACp1)] # correct

# view dead discard rates (forecast years are a placeholder; ratio of catch)
# [email protected][,ac(yrNow:yrTACp1)] # discard rate in forecast years
# [email protected][,ac(yrNow:yrTACp1)] # landings rate (1 - discard rate)

# add assumed stock size in forecast years
stockN <- 760  # stock size (million inds)
Fmsy <- 0.0812 # Fmsy harvest rate
Btrigger <- 858 # Trigger stock size
hr <- ifelse(stockN < Btrigger, Fmsy * (stockN/Btrigger), Fmsy) # hcr-based hr

# input stock size in advice year(s)
stkProj@stock.n[,ac(yrNow:yrTACp1)] <- stockN

# compute total removals (catch.n), followed by splits into landings.n and 
# discards.n. 
# then, compute total catch, landings, discards and harvest rate
stkProj@catch.n[,ac(yrNow:yrTACp1)] <- stkProj@stock.n[,ac(yrNow:yrTACp1)] * hr
stkProj@landings.n[,ac(yrNow:yrTACp1)] <- stkProj@catch.n[,ac(yrNow:yrTACp1)] * 
  stkProj@landings.n[,ac(yrNow:yrTACp1)]
stkProj@discards.n[,ac(yrNow:yrTACp1)] <- stkProj@catch.n[,ac(yrNow:yrTACp1)] * 
  stkProj@discards.n[,ac(yrNow:yrTACp1)]

stkProj@catch <- computeCatch(stkProj)
stkProj@landings <- computeLandings(stkProj)
stkProj@discards <- computeDiscards(stkProj)
stkProj@harvest <- computeHarvest(stkProj)

# compare output to reference
stfRef <- data.frame(
  model = "NSSK",
  year = ac(yrTAC),
  catch = c(1204),
  landings = c(1056),
  discards = c(148),
  fbar = c(0.072)
)

stfDet <- data.frame(
  model = "FLR",
  year = ac(yrTAC),
  catch = c(catch(stkProj[,ac(yrTAC)])),
  landings = c(landings(stkProj[,ac(yrTAC)])),
  discards = c(discards(stkProj[,ac(yrTAC)])),
  fbar = c(fbar(stkProj[,ac(yrTAC)]))
)

df <- merge(stfRef, stfDet, all = T)
df <- pivot_longer(df, cols = c(catch, landings, discards, fbar), 
  names_to = "variable", values_to = "value")
df <- pivot_wider(df, names_from = model, values_from = value)
df$percErr <- round((df$FLR - df$NSSK)/df$NSSK * 100, 1)

ggplot(df) + aes(x = year, y = percErr) + 
  facet_wrap(~variable) +
  geom_col(fill = 4, color = 4) +
  geom_hline(yintercept = 0, linetype = 1) +
  geom_hline(yintercept = c(-10,10), linetype = 3) + 
  theme_bw()

Additional information

Software Versions

• R version 4.4.3 (2025-02-28)
• FLCore: 2.6.20.9329
• FLasher: 0.7.1.9225
• ggplot2: 3.5.1
• ggplotFL: 2.7.0.9139
• mixfishtools: 0.5.0
• tidyr: 1.3.1
• knitr: 1.3.1
• Compiled: 2025-Mar-12