Introduction

Diagnostics for different assessment methods are usually produced by method-specific code (e.g. spict::plotspict.*, JABBA::jbplot_*, ss3diags), which makes it hard to compare methods and ties the plots to package internals. Here, extraction is separated from plotting:

1. each assessment is converted to a common summary, a cmpAssess, holding
   • ts, an FLQuants of annual biomass (stock), harvest, catch, \(B/B_{MSY}\) (bbmsy) and \(F/F_{MSY}\) (ffmsy), with lower and upper FLQuants (or iters) for the uncertainty;
   • refpts, an FLPar with bmsy, fmsy, msy, k, r and n;
   • index, a data.frame of observations and fits (series, year, obs, hat, lower, upper), resid and process residuals (series, year, residual), the production function curve (stock, yield), and prior/posterior densities or draws;
   • peels, a named list of cmpAssess from retrospective fits;
2. the cmp* functions plot one assessment or a named list of them, coloured by method.

Mohn’s rho, the hindcast and MASE, and the runs test are calculated from the common summary, so no method-specific code is needed. The converters are flbrp2cmp (an FLStock and FLBRP), spict2cmp (using only the exported SPiCT API), biodyn2cmp and jabba2cmp.

North Sea plaice is used to compare

• the age-structured assessment (ICES) and its equilibrium reference points (FLBRP), using exploitable biomass and the harvest rate (catch/biomass);
• SPiCT, fitted to catch and exploitable biomass, with priors from the FLBRP-derived Pella-Tomlinson parameters;
• JABBA, with priors from the FLBRP (FLRebuild::runJABBA).

library(FLCore)
library(FLBRP)
library(ggplotFL)
library(ggplot2)
library(spict)
library(mpb)

Assessments

data(icesdata, package="icesdata")
stk=window(icesdata[["ple.27.420"]], end=2020)

eq =FLRebuild::eql(stk, model="bevholtSV")
rf =FLRebuild::refptsEB(eq)
pt =FLRebuild::pellatParams(FLPar(msy =rf["msy",   "yield"],
                                  bmsy=rf["msy",   "eb"],
                                  k   =rf["virgin","eb"]))
idx=FLRebuild::ebiomass(stk)
psi=c(idx[, 1])/c(pt["k"])

asFLBRP=flbrp2cmp(stk, eq, name="FLBRP", biomass="eb", harvest="hr")
asFLBRP

cv2sd=function(x) sqrt(log(1+x^2))
scale=1000
yrs  =an(dimnames(catch(stk))$year)

inp=list(obsC=c(catch(stk))/scale, timeC=yrs,
         obsI=c(idx)/scale,        timeI=yrs,
         dteuler=1/4)
inp$priors=list(
  logr     =c(log(c(pt["r"])),       cv2sd(0.3), 1),
  logK     =c(log(c(pt["k"])/scale), cv2sd(0.5), 1),
  logn     =c(log(c(pt["p"])+1),     cv2sd(0.3), 1),
  logbkfrac=c(log(psi),              cv2sd(0.3), 1),
  logsdi   =c(log(0.1),              0.5,        1))
inp=check.inp(inp, verbose=FALSE)

fit=fit.spict(inp)
fit=calc.osa.resid(fit)
fit=calc.process.resid(fit)
fit=retro(fit, nretroyear=5)

asSPiCT=spict2cmp(fit, name="SPiCT", scale=scale, process="generic")

jb=FLRebuild::runJABBA(stk, method="flbrp", quick=TRUE, nc=2)
asJABBA=jabba2cmp(jb, name="JABBA")

asmt=list(FLBRP=asFLBRP, SPiCT=asSPiCT)
if (exists("asJABBA")) asmt$JABBA=asJABBA
names(asmt)

Time series

cmpTS(asmt)

cmpStatus(asmt)

Since ts is an FLQuants, ggplotFL can be used directly, e.g. to compare \(B/B_{MSY}\) across methods

bb=FLQuants(lapply(asmt, function(x) x$ts$bbmsy))

ggplot(bb, aes(year, data, colour=qname)) +
  geom_line() +
  geom_hline(yintercept=1, linetype=2) +
  labs(x="Year", y=expression(B/B[MSY]), colour=NULL) +
  theme_minimal()

or plot a single method with plot

plot(asmt$SPiCT$ts[c("stock", "harvest", "catch")])

Kobe phase plot

cmpKobe(asmt)

Production functions

The FLBRP curve is the equilibrium yield against equilibrium exploitable biomass; for the biomass dynamic models it is the Pella-Tomlinson production function. Points are the realised surplus production \(B_{t+1}-B_t+C_t\) of each method.

cmpProduction(asmt)

Fits to the data

cmpFit(asmt)

Residuals

Observation residuals (one-step-ahead for SPiCT, \(\log(obs/hat)\) otherwise), with runs tests.

cmpRuns(asmt)

cmpACF(asmt)

cmpQQ(asmt)

cmpProcess(asmt)

cmpACF(asmt, "process")

Retrospective analysis

Mohn’s rho is calculated by mohnsRho from the peels, for methods that have them.

cmpRetro(asmt)

Hindcast cross-validation

One-year-ahead predictions of the index from each peel, with the MASE (cmpMASE) relative to a naive (last observation) prediction.

cmpHindcast(asmt, series="Index 1")

hc=cmpHindcastDF(asmt)
do.call(rbind, lapply(split(hc, hc[, c("method", "series")], drop=TRUE), function(x)
  data.frame(method=x$method[1], series=x$series[1],
             MASE=cmpMASE(x$obs, x$pred, x$naive))))

Priors and posteriors

Priors are dashed and posteriors solid. For SPiCT the posterior is the normal approximation on the log scale; for JABBA it is the MCMC sample (with quick=TRUE, so a bimodal posterior may indicate chains that have not mixed).

cmpPriors(asmt, parameters=c("r", "k", "n", "psi"))