## Warning: package 'knitr' was built under R version 4.4.3

Introduction

Biomass-dynamic models such as SPiCT represent uncertainty in population dynamics through a process-error term. If the production function and state equation are correctly specified, and all relevant covariates are included, these process residuals should be independent and identically distributed around zero (i.e. white noise), with approximately constant variance. In practice, however, the process-error term often compensates for structural misspecification and non-stationarity, for example time-varying productivity, changing selectivity or catchability, e.g. in stocks which are not the main target.

Here we propose a set of diagnostics for process error using an example data set. The aim is not to formally test ypotheses, but to provide a visual summary to highlights when a model is explain model misspecification via process-error focussing on four aspects namely:

  1. Randomness and variance over time
  2. Trends and shifts including potential bias at high or low biomass.
  3. Magnitude and distribution of residuals.
  4. Autocorrelation temporal structure beyond white noise.

Framework

Use is made of existing diagnostics used for model fits, such as runs tests and residual diagnostics to help identify non-stationarity, the impact of environmental drivers, and the balance between stochasticity and model misspecification. If process residuals appear random with approximately constant variance, this suggests a correctly specified production function with genuine stochastic variation, conditional on the chosen model structure and fixed parameters. Positive or negative trends in residuals may indicate production-curve misspecification or structural changes, for example an incorrect shape parameter, time-varying productivity, changes in discarding or landings, or regime shifts. Consistently large residuals relative to the assumed process variance may reflect high biological variability, such as recruitment-driven dynamics, but may also signal structural model inadequacy or a mis-specified variance structure. Significant temporal autocorrelation may point to environmental forcing, omitted covariates, or serially correlated implementation error rather than white-noise process error. But we live in a world full of red noise.

Data and R Implementation

The following R code loads precomputed process-error, e.g. output from SPiCT, JABBA or SS3, and runs the set of diagnostic plots for a selected stock.

library(spict)
library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)
load("P:/rfmo/ices/WKBMSYSPiCT4/fle.27,21-23/finalFit.RData")

fit = calc.process.resid(fit)
getSpictPrd<-function(fit, xMax = 1.3, nX = 200) {

  ## -- Parameters for Pella–Tomlinson curve --
  n   = get.par("logn",  fit, exp = TRUE)[1, "est"]
  msy = get.par("MSY",   fit, exp =!TRUE)[1, "est"]
  K   = get.par("logK",  fit, exp = TRUE)[1, "est"]

  x = seq(0.01, xMax, length.out = nX)    # relative biomass B/K

  gamma = (1 / (n - 1)) * n^(n / (n - 1))
  Pfun  = gamma * msy * (x - x^n)

  prodFun = data.frame(
    bRel = x,
    b    = x * K,
    sp   = Pfun)

  ## -- Realised production trajectory: P_t = B_{t+1} - B_t + C_t --
  # Biomass on SPiCT grid
  Bfull = get.par("logB", fit, exp = TRUE)[, "est"]
  tFull = fit$inp$time

  # Keep integer years
  sel = (tFull - floor(tFull)) == 0
  B   = Bfull[sel]
  t   = tFull[sel]

  # Catches aligned to those intervals
  C = fit$inp$obsC
  C = C[seq_len(length(B) - 1)]

  P = B[-1] - B[-length(B)] + C

  prodTraj = data.frame(
    year = tail(t, -1),          # end of interval
    sp   = P,
    bRel = B[-length(B)] / K,    # starting biomass of interval, relative to K
    b    = B[-length(B)])

  names(fit$process.resid)=c("year","bRsdl","fRsdl")

  list(pf=prodFun, traj=merge(prodTraj,fit$process.resid))}

pf=getSpictPrd(fit)

ggplot(pf$traj)+
  geom_line(aes(b,sp), data=pf$pf)+
  geom_path( aes(b,sp))+
  geom_point(aes(b,sp,col=year))+
  geom_text(aes(label=year, x = b, y = sp), size=3, data=subset(pf$traj,min(year)==year)) +
  geom_point(aes(x = b, y = sp),            data=subset(pf$traj,min(year)==year)) +
  geom_text(aes(label=year, x = b, y = sp), col="blue",  size=3, data=subset(pf$traj,max(year)==year)) +
  geom_point(aes(x = b, y = sp), col="blue",        data=subset(pf$traj,max(year)==year))

library(dplyr)
library(ggplot2)
library(grid)

library(plyr)
windowOverlap = function(dat, yearCol = "year", step = 10, overlap=5) {
  
  years = dat[[yearCol]]
  yMin  = min(years, na.rm = TRUE)
  yMax  = max(years, na.rm = TRUE)
  
  centres = seq(
    from = ceiling(yMin / step) * step,
    to   = floor(yMax / step) * step,
    by   = step
  )
  
  outList = lapply(centres, function(y0) {
    from = y0 - overlap
    to   = y0 + 9+overlap
    
    sel = years >= from & years <= to
    w   = dat[sel, ]
    w$centre = y0
    w$from   = from
    w$to     = to
    w
  })
  
  names(outList) = as.character(centres)
  
  out = ldply(outList, .id = ".id")
  out}


trajWin=windowOverlap(pf$traj,overlap=2)

trajArrows =
  trajWin %>%
  arrange(.id, year) %>%
  group_by(.id) %>%
  mutate(bEnd  = lead(b),
         spEnd = lead(sp)) %>%
  ungroup() %>%
  filter(!is.na(bEnd), !is.na(spEnd))

ggplot(trajWin) +
  facet_wrap(~ .id) +
  geom_line(data = pf$pf, aes(x = b, y = sp)) +
  #geom_path(aes(x = b, y = sp),col="grey85") +
  geom_segment(data = trajArrows,
               aes(x = b, y = sp,col=year-from,
                   xend = bEnd, yend = spEnd),linewidth = 0.2,
               arrow = arrow(length = unit(0.15, "cm"))) +
  geom_text(aes(label=year, x = b, y = sp), size=3, data=subset(trajWin,from==year)) +
  geom_point(aes(x = b, y = sp),            data=subset(trajWin,  from==year)) +
  geom_text(aes(label=year, x = b, y = sp), col="blue",  size=3, data=subset(trajWin,to==year)) +
  geom_point(aes(x = b, y = sp), col="blue",        data=subset(trajWin,  to==year))

Production Curves

The surplus production function \(P(B)\) (e.g. Schaefer, Pella–Tomlinson, Fletcher) is the expected \(B_{t+1}-B_t + C_t\) as a function of \(B_t\) under equilibrium assumptions. The realised surplus production (\(B_{t+1}-B_t + C_t\)) in year \(t\) is plotted against \(B_t\) and overlayed on the equilibrium production curve. The residual \(SP_t - P(B_t)\) is the process error for that year (on the biomass scale), assuming observation error in biomass is negligible or already accounted for. This allows a check that process error is small, stationary, and biomass‑dependent as assumed.

Potential patterns

  1. Good fit / stationary process error
    • Points should be scattered symmetrically around the curve across the range of biomasses, with no trend over time or with biomass.
    • Residuals should look have a mean of zero and be white noise; supporting the assumption of stationary productivity and correctly specified \(P(B)\).
  2. Systematic bias with biomass (shape mis‑specification)
    • Points at low \(B\) all lie above the curve while points at high \(B\) lie below (or vice versa).
    • Indicates the chosen production curve is the wrong shape (e.g.wrong n in Pella–Tomlinson); may need a different functional form or to allow time‑varying productivity.
  3. Cycles indicating non‑stationarity

    • When points connected in time order, there may be “hooks” or cyclic patterns:

    • Upward or clockwise hooks: Periods where production increases, biomass then increases, then production drops and biomass drops, forming loops around the curve. These can indicate productivity regime shifts or transient dynamics not captured by a stationary \(P(B)\).

    • Downward or counter‑clockwise hooks: Higher production for a given biomass during declines than during recoveries (i.e. reduced productivity on the way back up), driven by age‑structure or ecosystem changes.

Interpretation is that for the same \(B_t\), \(SP_t\) differs systematically by period, i.e. time‑varying productivity or unmodelled processes.

  1. Trend in residuals (time‑varying productivity)

    • If residuals become increasingly negative (points consistently below the curve in recent years), it suggests productivity has declined relative to the assumed stationary curve (e.g. environmental change, depensation, changes in vital rates).
    • A run of positive residuals implies higher‑than‑assumed productivity (potential rebuilding regimes, ecosystem changes that favour the stock).

  2. Heteroscedasticity / episodes of large process error

    • Clusters of large positive or negative deviations in particular years or periods indicate episodes where the model is missing something (e.g. unrecorded mortality, strong year classes, changes in selectivity / catchability).
    • Check whether the estimated process deviations are random or structurally patterned.

Diagnostics

  • Model structure: If the cloud of points looks like noisy draws around the curve, the process equation is structurally adequate; if not, consider alternative \(P(B)\), time‑varying productivity, or adding explicit covariates (e.g. environment).
  • Check process vs observation error balance: Large, structured deviations in the phase plot may indicate that what is attributed to process error is actually mis‑specified observation error (e.g. biased CPUE) or mis‑specified selectivity/catch; revisiting the observation model might reduce the apparent process error.
  • Regime identification: Colouring points by period or regime and seeing distinct clusters above/below the curve is a potential emprirical check for non‑stationary surplus production.

Example interpretation

If the \(B_{t+1}-B_t+C_t\) vs \(B_t\) phase plot shows early‑period points centred on the curve, but recent years consistently below the curve at similar biomass levels, with a counter‑clockwise loop, then this may be evidence that productivity has declined in recent decades and that recoveries are occurring under lower surplus production than declines at comparable biomass; this could argue for time‑varying productivity or for more conservative advice.

Four-Panel Diagnostic Plot

A single figure with four panels is used to summarise process error: i.e. 1) residuals vs time, 2) residuals vs fitted biomass, 3) autocorrelation function (ACF), and 4) a histogram with normal curve for comparison.

1. Residuals vs time

p_time = ggplot(pf$traj, aes(x=year, y=bRsdl))  +
  geom_hline(yintercept = 0, colour = "grey60") +
  geom_line( colour = "steelblue") +
  geom_point(colour = "steelblue") +
  labs(x = "Year", y = "Process residual",
       title = "Over time") +
  theme_bw()

2. Residuals vs fitted biomass

p_fit = ggplot(pf$traj, aes(x=b, y=bRsdl)) +
  geom_hline(yintercept = 0, colour = "grey60") +
  geom_point(alpha = 0.6, colour = "steelblue") +
  labs(x = "Fitted biomass (or state)", y = "Process residual",
       title = "Residuals vs Biomass") +
  theme_bw()

3. ACF (pre-compute then plot with ggplot)

acf_obj = acf(pf$traj$bRsdl, plot = FALSE)
acf_df = with(acf_obj, data.frame(lag = lag, acf = acf))[-1, ]
crit = qnorm(0.975) / sqrt(nrow(pf$traj))
p_acf = ggplot(acf_df, aes(x = lag, y = acf)) +
  geom_hline(yintercept = 0, colour = "grey60") +
  geom_hline(yintercept = c(-crit, crit), linetype = "dashed",
             colour = "grey60") +
  geom_col(fill = "steelblue") +
  labs(x = "Lag", y = "ACF",
       title = "Autocorrelation") +
  theme_bw()

4. Histogram + normal overlay

p_hist = ggplot(pf$traj, aes(x=bRsdl)) +
  geom_histogram(aes(y = ..density..), bins = 20,
                 fill = "steelblue", colour = "white", alpha = 0.7) +
  stat_function(fun = dnorm,
                args = list(mean = mean(pf$bRsdl), sd = sd(pf$bRsdl)),
                colour = "red", linewidth = 0.8) +
  labs(x = "Residual", y = "Density",
       title = "Distribution") +
  theme_bw()
(p_time | p_fit) /
  (p_acf | p_hist)

Interpretation

Randomness and Variance

The top-left panel (“Over time”) shows process residuals plotted against year. For a well-specified residuals would be expected to fluctuate around zero with no obvious trend and roughly constant variance. However, residuals are mostly positive after about 2000, with an upward shift in the mean. The spread also appears slightly larger in recent years. Together this suggest that the process-error term is not stationary and probably reflects systematic model bias rather than purely stochastic fluctuations.

Magnitude and Distribution

The bottom-right panel (“Distribution”) shows the empirical distribution of residuals with a normal density overlay. The shape is roughly Gaussian but centred above zero. Residual magnitudes are not particulary extreme, i.e. relative to the assumed process standard deviation, but the mean of 0.2 points to bias in the production curve or to a missing process, rather than to occasional large shocks superimposed on an otherwise correct specified model.

Autocorrelation

The bottom-left panel (“Autocorrelation”) is the ACF of the process residuals. Several low-lag bars exceed the approximate 95% bounds, and there is a pattern of positive autocorrelation at short lags followed by negative values at longer lags. This is inconsistent with a white-noise process and suggests slow, temporally structured deviations from the model, for example an environmental regime, time-varying productivity, or changes in catchability.

Summary

Together, these plots suggest that the process-error term in this example is accounting for model misspecification or non-stationarity rather than representing pure stochastic noise. This is because

  • Residuals exhibit a shift in mean and variance over time.
  • The model tends to under-predict biomass at higher levels, consistent with a biased production curve.
  • The residual distribution is centred above zero, even though individual deviations are not huge.
  • Autocorrelation is present, indicating red rather than the white noise assumed.

Process-error appears not to be purely stochastic and so confidence intervals or fractiles from the fitted model will be biased.

Possible next steps are to

  1. Explore time-varying productivity or selectivity, or alternative production-function shapes.
  2. Consider including covariates or explicit regime-shift structures if there is supporting evidence.
  3. Evaluate HCR performance in an operating-model/MSE framework or an ensemble that reproduces the process-error properties, rather than basing advice solely on a singlr fitted model.
  4. Give up

The same diagnostic framework can be applied systematically across stocks to identify cases where process error behaves as desired and those where it is clearly signalling deeper structural issues. These diagnostics can be applied to any model with an implicit or explicit production function.