Introduction

SPiCT (Pedersen and Berg, 2017) provides base graphics plots (plot.spictcls, plotspict.*). This vignette shows the equivalent plots built with ggplot2 and ggplotFL, in the same style as the diagnostics used for the other assessment models (time series with confidence intervals, Kobe phase plots, production functions, runs tests on residuals, retrospective and hindcast cross-validation, and priors against posteriors).

The plotSpict* functions are thin wrappers: the fit is converted to the method-agnostic cmpAssess summary with spict2cmp (which uses only the exported SPiCT API) and plotted with the generic cmp* functions, which can also overlay other assessment methods (see the assessment-comparison vignette).

The example is North Sea plaice, fitted as in the spict vignette, with priors for \(r\), \(K\), \(n\) and the initial depletion derived from the FLBRP reference points and exploitable biomass as the index.

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

library(mpb)

Fit

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"])

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.process.resid(fit)
fit=calc.osa.resid(fit)
fit$opt$convergence

The retrospective (spict::retro) refits the model to five peels; the same peels are used for the hindcast.

fit=retro(fit, nretroyear=5)

Time series

spictTS extracts annual estimates and 95% confidence intervals as a data.frame, and plotSpictTS plots them. Dashed lines are the deterministic reference points (\(B_{MSY}\), \(F_{MSY}\), MSY); scale converts back to tonnes.

plotSpictTS(fit, scale=scale)

The same plot from the common summary, which is what plotSpictTS does internally:

sp=spict2cmp(fit, scale=scale)
sp

Stock status relative to MSY reference points.

plotSpictStatus(fit)

As FLQuants

spict2FLQuants converts the fit to FLQuants. With iter the uncertainty is represented by log-normal draws, so that ggplotFL::plot shows the median and quantile ribbons as for any other FLR model.

flq=spict2FLQuants(fit, iter=500)
plot(flq)

Without iter the point estimates can be combined with other FLQuants, e.g. the exploitable biomass from the age-structured assessment.

bsp=FLQuant(c(spict2FLQuants(fit, "B")[["B"]])*scale, dimnames=dimnames(idx))
cmp=FLQuants(SPiCT=bsp, Assessment=idx)

ggplot(cmp, aes(year, data, colour=qname)) +
  geom_line() +
  labs(x="Year", y="Exploitable biomass (t)", colour=NULL) +
  theme_minimal()

Kobe phase plot

The terminal year is shown as a large point with its 95% confidence intervals.

plotSpictKobe(fit)

Production function

The Pella-Tomlinson production function implied by the SPiCT estimates, with the realised surplus production \(B_{t+1}-B_t+C_t\) in each year.

plotSpictProduction(fit)

Model fit

Observed (points) and fitted catch and index, with 95% confidence intervals.

plotSpictFit(fit)

Residuals

One-step-ahead residuals

The 3-sigma band is green if the runs test is passed and red if it fails; points outside the band are outlined in red.

plotSpictRuns(fit)

plotSpictACF(fit)

plotSpictQQ(fit)

Process residuals

plotSpictProcess(fit)

plotSpictACF(fit, "process")

Retrospective analysis

plotSpictRetro(fit)

Hindcast cross-validation

Each peel predicts the index one year ahead as \(qB_{T+1}\), where \(T\) is the peel’s last year; MASE is calculated by cmpMASE from these predictions. Unlike spict::hindcast, catches are also removed from the peeled years.

plotSpictHindcast(fit)

Priors and posteriors

plotSpictPriors(fit)

References

Pedersen, M.W. and Berg, C.W. (2017). A stochastic surplus production model in continuous time. Fish and Fisheries, 18(2), 226-243.

ggplot2 and ggplotFL plots for SPiCT fits