This vignette shows how to go from an age-structured stock assessment to a surplus production model fitted with SPiCT (Pedersen and Berg, 2017), using the equilibrium reference points of the age-structured model as priors. The steps are
FLStock (North Sea plaice from
icesdata);FLBRP
(FLRebuild::eql);FLBRP reference points (MSY, \(B_{MSY}\), \(K\)) into Pella-Tomlinson parameters (\(r\), \(n\)) and an initial depletion \(B_0/K\), which are used as SPiCT
priors;FLStock catch and a
biomass index, and fit;FLBRP
equivalents;biodyn object and project recovery
at zero catch from \(0.3B_{MSY}\).FLRebuild and icesdata mask several
FLCore generics (e.g. fwd, msy,
bmsy), so they are not attached and are called with
::. For the same reason FLBRP::refpts and
mpb::refpts are used explicitly.
FLRebuild::eql fits a Beverton-Holt stock-recruitment
relationship (steepness parameterisation, with a prior on steepness) and
returns an FLBRP. FLRebuild::refptsEB adds
exploitable biomass (eb, i.e. biomass weighted by
selectivity) to the standard reference points.
eq=FLRebuild::eql(stk, model="bevholtSV")
rf=FLRebuild::refptsEB(eq)
rf[c("virgin", "msy"), c("harvest", "yield", "ssb", "biomass", "eb")]A biomass dynamic model has a single biomass, so we have to choose
which age-structured biomass it represents. Catch is taken from the
exploitable biomass, so exploitable biomass is the natural choice;
spawning stock biomass ("ssb") can be used by changing
quant below.
For the Pella-Tomlinson model, surplus production is
\[P(B)=\frac{r}{p}B\left(1-\left(\frac{B}{K}\right)^p\right),\]
with \(B_{MSY}/K=(1+p)^{-1/p}\),
\(MSY=rB_{MSY}/(1+p)\) and \(F_{MSY}=r/(1+p)\). SPiCT uses \(n=p+1\) and the same definition of \(r\) (\(r=mn^{n/(n-1)}/K\) in the SPiCT source,
where \(m\) is MSY), so
biodyn, FLRebuild and SPiCT share \(r\) and \(K\), and \(n=p+1\).
spictPriors uses FLRebuild::pellatParams,
which solves for \(p\) from \(B_{MSY}/K\) and then for \(r\) from MSY and \(B_{MSY}\). Initial depletion is the first
year biomass divided by \(K\). The
priors are log-normal; SPiCT takes the mean and standard deviation on
the log scale, so the CVs are converted with \(\sigma=\sqrt{\log(1+CV^2)}\).
SPiCT is fitted to catch and to the biomass index from the
FLStock. The index is an output of the age-structured
assessment rather than a survey, so it has almost no observation error;
without some information the index standard deviation (sdi)
collapses towards zero. A weak prior centred on 0.1 (log-scale SD 0.5)
keeps it away from zero.
Values are rescaled to thousand tonnes, which helps the optimiser, and the Euler time step is set to a quarter of a year to keep the runtime short.
scale=1000
yrs =an(dimnames(catch(stk))$year)
inp=spictInp(stk, eq, quant, scale=scale, dteuler=1/4, cv=cv, sdi=c(0.1, 0.5))
inp$priors[c("logr", "logK", "logn", "logbkfrac", "logsdi")]Only the priors are set; initial values (inp$ini) are
left to SPiCT. Setting the initial values to the prior means led to a
poorer local optimum in testing.
c(convergence=fit$opt$convergence,
pdHess =fit$pdHess,
finiteSE =all(is.finite(sqrt(diag(fit$cov.fixed)))))convergence should be 0 and the Hessian positive
definite.
The SPiCT figures use the ggplotFL-based plotSpict*
functions in mpb, which convert the fit to a common assessment summary
(spict2cmp) so the same plots can be used for other
assessment methods (see the assessment-comparison
vignette).
One-step-ahead residuals should show no bias, autocorrelation or departure from normality.
spictPar extracts the estimate and CV of a parameter
from a SPiCT fit on the natural scale. The deterministic reference
points (Bmsyd, Fmsyd) are used, since these
correspond to the equilibrium reference points of the
FLBRP; the stochastic versions (Bmsy,
Fmsy) also depend on the process error.
The FLBRP harvest rate at MSY is \(MSY/B_{MSY}\) for the chosen biomass, which
is comparable with the SPiCT \(F_{MSY}\); the age-based \(F_{bar}\) at MSY is shown separately.
Biomass and catch are in thousand tonnes.
tl =match(max(yrs), fit$inp$time)
bmsy=c(rf["msy", quant])
hrTerm=c(catch(stk)[, ac(max(yrs))])/c(idx[, ac(max(yrs))])
flbrp=c(r =c(pt["r"]),
K =c(pt["k"])/scale,
n =c(pt["p"])+1,
BmsyK =bmsy/c(pt["k"]),
MSY =c(rf["msy", "yield"])/scale,
Bmsy =bmsy/scale,
Fmsy =c(rf["msy", "yield"])/bmsy,
BBmsy =c(idx[, ac(max(yrs))])/bmsy,
FFmsy =hrTerm/(c(rf["msy", "yield"])/bmsy))
spct=rbind(r =gp("logr", fit),
K =gp("logK", fit),
n =gp("logn", fit),
BmsyK=c(gp("logBmsyd", fit)[1]/gp("logK", fit)[1], NA),
MSY =gp("logm", fit),
Bmsy =gp("logBmsyd", fit),
Fmsy =gp("logFmsyd", fit),
BBmsy=gp("logBBmsy", fit, tl),
FFmsy=gp("logFFmsy", fit, tl))
cmp=data.frame(FLBRP =flbrp,
SPiCT =spct[names(flbrp), 1],
SPiCT.cv =spct[names(flbrp), 2],
ratio =spct[names(flbrp), 1]/flbrp)
knitr::kable(cmp, digits=3,
caption=paste0("FLBRP (", quant, ") and SPiCT estimates, ",
"terminal year ", max(yrs), "."))SPiCT estimates a catchability \(q\)
for the index, even though the index is itself a biomass. If \(q<1\) SPiCT biomass is larger than the
FLStock biomass, so \(B/B_{MSY}\) can differ even when the
production functions are similar. Fixing \(q=1\) (e.g. a tight prior on
logq) forces the two biomass scales to agree.
Three figures compare the age-structured stock (FLStock
and FLBRP) with SPiCT: the time series, the production
functions and the process error. All biomass and catch values are in
thousand tonnes.
spAnn extracts a SPiCT quantity at the start of each
year, with its 95% confidence interval.
ann =match(yrs, round(fit$inp$time, 6))
spAnn=function(nm) {
x=get.par(nm, fit, exp=TRUE)[ann, , drop=FALSE]
data.frame(year=yrs, lower=x[, "ll"], est=x[, "est"], upper=x[, "ul"])}
C =c(catch(stk))/scale
B =c(idx)/scale
nY =length(yrs)
fmsyBR=c(rf["msy", "yield"])/bmsy
clrs =c(FLStock="grey40", FLBRP="grey40", "Pella-Tomlinson (FLBRP)"="black",
SPiCT="red")Biomass and fishing mortality on the absolute scale (top) and
relative to each model’s own reference points (bottom). For the
FLStock, fishing mortality is the harvest rate (catch
divided by the chosen biomass) and \(F_{MSY}\) is \(MSY/B_{MSY}\) from the FLBRP;
SPiCT \(F\) is an instantaneous rate,
so the two are only approximately comparable on the absolute scale.
Shaded areas are the SPiCT 95% confidence intervals.
pnl=c("Biomass (kt)", "F or harvest rate", "B/BMSY", "F/FMSY")
fls=data.frame(year=yrs, lower=NA, upper=NA, source="FLStock")
trj=rbind(cbind(fls, est=B, panel=pnl[1]),
cbind(fls, est=C/B, panel=pnl[2]),
cbind(fls, est=B*scale/bmsy, panel=pnl[3]),
cbind(fls, est=C/B/fmsyBR, panel=pnl[4]),
cbind(spAnn("logB"), source="SPiCT", panel=pnl[1]),
cbind(spAnn("logFs"), source="SPiCT", panel=pnl[2]),
cbind(spAnn("logBBmsy"), source="SPiCT", panel=pnl[3]),
cbind(spAnn("logFFmsy"), source="SPiCT", panel=pnl[4]))
trj$panel=factor(trj$panel, levels=pnl)
ggplot(trj, aes(year, est, colour=source))+
geom_ribbon(aes(ymin=lower, ymax=upper, fill=source),
data=subset(trj, source=="SPiCT"), alpha=0.2, colour=NA,
show.legend=FALSE)+
geom_line()+
geom_hline(aes(yintercept=1), data=data.frame(panel=factor(pnl[3:4], levels=pnl)),
linetype=2)+
facet_wrap(~panel, scales="free_y")+
scale_colour_manual(values=clrs)+scale_fill_manual(values=clrs, guide="none")+
expand_limits(y=0)+
labs(x="Year", y="", colour="")+
theme_bw()+theme(legend.position="bottom")The FLBRP equilibrium curve is obtained by evaluating
the FLBRP over a range of \(F_{bar}\) and plotting equilibrium yield
against the chosen biomass. It is compared with the Pella-Tomlinson
curve fitted to it (FLRebuild::pellatParams, used as the
SPiCT prior) and with the SPiCT surplus production curve \(P(B)=\gamma m\left(B/K-(B/K)^n\right)\),
\(\gamma=n^{n/(n-1)}/(n-1)\). Points
are the annual surplus production \(B_{t+1}-B_t+C_t\) calculated from the
FLStock biomass and from the SPiCT biomass. Large symbols
mark (\(B_{MSY}\), MSY) for each
model.
eqp=eq
fbar(eqp)=FLQuant(seq(0, 5*c(rf["msy", "harvest"]), length.out=101))
eqB=switch(quant, eb=FLRebuild::ebiomass(eqp), ssb=ssb(eqp))
n=gp("logn", fit)[1]
K=gp("logK", fit)[1]
m=gp("logm", fit)[1]
pSP=function(B) n^(n/(n-1))/(n-1)*m*(B/K-(B/K)^n)
pBR=function(B, r=c(pt["r"]), k=c(pt["k"])/scale, p=c(pt["p"])) r/p*B*(1-(B/k)^p)
bGrid=seq(0, max(K, c(pt["k"])/scale), length.out=201)
crv=rbind(data.frame(biomass=c(eqB)/scale, yield=c(catch(eqp))/scale, source="FLBRP"),
data.frame(biomass=bGrid, yield=pBR(bGrid), source="Pella-Tomlinson (FLBRP)"),
data.frame(biomass=bGrid, yield=pSP(bGrid), source="SPiCT"))
crv=subset(crv, yield>=0)
spB=spAnn("logB")$est
obs=rbind(data.frame(biomass=B[-nY], yield=diff(B)+C[-nY], source="FLStock"),
data.frame(biomass=spB[-nY], yield=diff(spB)+C[-nY], source="SPiCT"))
ref=data.frame(biomass=c(bmsy/scale, gp("logBmsyd", fit)[1]),
yield =c(c(rf["msy", "yield"])/scale, m),
source =c("FLBRP", "SPiCT"))
ggplot(crv, aes(biomass, yield, colour=source))+
geom_point(data=obs, alpha=0.5, size=1)+
geom_line(aes(linetype=source))+
geom_point(data=ref, size=4, shape=18)+
geom_hline(yintercept=0, colour="grey70")+
scale_colour_manual(values=clrs)+
scale_linetype_manual(values=c(FLBRP=1, "Pella-Tomlinson (FLBRP)"=2, SPiCT=1),
guide="none")+
labs(x=paste0("Biomass (", quant, ", kt)"), y="Surplus production (kt)",
colour="", caption="Dashed: Pella-Tomlinson fitted to the FLBRP (SPiCT prior)")+
theme_bw()+theme(legend.position="bottom")Process error is summarised with the four-panel diagnostic used at the ICES SPiCT benchmark (WKBMSYSPiCT4): residuals over time, residuals against biomass, their autocorrelation, and their distribution. Following Walters et al. (2008), the realised surplus production \(SP_t=B_{t+1}-B_t+C_t\) is compared with the production function \(P(B_t)\). If \(P(B)\) is correctly specified the residuals are white noise around zero; trends, dependence on biomass or autocorrelation show that process error is absorbing misspecification or non-stationary productivity.
For SPiCT the residuals are its own process residuals
(calc.process.resid, computed after the fit). These are
standardised Euler-step residuals averaged within each year, so they are
multiplied by \(\sigma_B/\sqrt{\Delta
t}\) to give annual log deviations, comparable with the estimated
sdb (\(\sigma_B\)). For
the age-structured stock the residual is the log deviation of the
FLStock biomass from a one-year Pella-Tomlinson step with
the FLBRP-derived parameters, \(\log\left(B_{t+1}/(B_t+P(B_t)-C_t)\right)\).
Residuals are indexed by the start year \(t\) of each interval.
spictPrd returns the SPiCT production curve and annual
residuals, prdResid the residuals for a biomass and catch
series and a production function, and plotPE draws the four
panels.
sdb=gp("logsdb", fit)[1]
peSP=spictPrd(fit)$traj
peBR=prdResid(B, C, yrs, pBR)
pe =rbind(cbind(peBR[, c("year", "b", "rsdl")], source="FLStock"),
cbind(peSP[, c("year", "b", "rsdl")], source="SPiCT"))
peSD=sapply(split(pe$rsdl, pe$source), sd, na.rm=TRUE)
peAC=sapply(split(pe$rsdl, pe$source), function(x) acf(x, plot=FALSE, na.action=na.pass)$acf[2])
plotPE(pe, clrs[c("FLStock", "SPiCT")])+
patchwork::plot_annotation(subtitle=sprintf(
"SD: FLStock %.3f, SPiCT %.3f; SPiCT sdb %.3f. Lag-1 ACF: FLStock %.2f, SPiCT %.2f",
peSD["FLStock"], peSD["SPiCT"], sdb, peAC["FLStock"], peAC["SPiCT"]))spict2biodyn copies \(r\), \(K\), \(p=n-1\) and \(B_0/K\) into the biodyn
parameters, and the catch, index and SPiCT annual biomass into the
corresponding slots, rescaling back to tonnes.
mpb::refpts gives the same \(F_{MSY}\), \(B_{MSY}\) and MSY as the deterministic
SPiCT reference points (to rounding).
To link the production function to rebuilding time, both the SPiCT
and the FLBRP-derived Pella-Tomlinson parameters are projected from
\(0.3B_{MSY}\) with no catch.
biodyn is a discrete annual model whereas SPiCT is in
continuous time, so recovery times for the same parameters can differ
slightly.
rcv=rbind(cbind(model="SPiCT", recoverBD(params(bd))),
cbind(model="FLBRP", recoverBD(FLPar(r=c(pt["r"]), k=c(pt["k"]), p=c(pt["p"])))))
tRecover=sapply(split(rcv, rcv$model), function(x) min(x$year[x$BBmsy>=1]))
tRecoverggplot(rcv, aes(year, BBmsy, colour=model))+
geom_line()+
geom_hline(yintercept=1, linetype=2)+
geom_vline(xintercept=tRecover, linetype=3, colour=c("grey40", "red"))+
scale_colour_manual(values=c(FLBRP="grey40", SPiCT="red"))+
labs(x="Years since start of recovery", y=expression(B/B[MSY]), colour="")+
coord_cartesian(xlim=c(0, 30))+
theme_bw()+theme(legend.position="bottom")Time to recover depends on both \(r\) and the shape. With no catch the per-capita growth rate is \((r/p)(1-(B/K)^p)\), which tends to \(r/p\) as \(B\rightarrow 0\); for \(p<1\) (i.e. \(B_{MSY}/K<0.5\)) this is much larger than \(r\). A lower \(r\), or a larger \(p\), therefore means a longer rebuilding time.
inp$obsI and inp$timeI.quant="ssb" gives a much lower \(n\) (a \(B_{MSY}/K\) around 0.25 for plaice), since
recruitment compensation acts on spawning biomass.FLRebuild::eql; different choices give
different priors.Pedersen, M. W. and Berg, C. W. (2017). A stochastic surplus production model in continuous time. Fish and Fisheries, 18: 226-243.
Walters, C. J., Hilborn, R. and Christensen, V. (2008). Surplus production dynamics in declining and recovering fish populations. Canadian Journal of Fisheries and Aquatic Sciences, 65: 2536-2551.