--- title: "Fitting SPiCT with priors derived from an FLBRP" author: "Laurence Kell" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: true fig_width: 7 fig_height: 4.5 vignette: > %\VignetteIndexEntry{Fitting SPiCT with priors derived from an FLBRP} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, knitr, include=FALSE} knitr::opts_chunk$set(echo=TRUE, message=FALSE, warning=FALSE, fig.align="center", comment="#>") pkgs=c("FLCore", "FLBRP", "FLRebuild", "icesdata", "spict", "mpb", "ggplot2") ok =all(vapply(pkgs, requireNamespace, logical(1), quietly=TRUE)) knitr::opts_chunk$set(eval=ok) ``` # Introduction 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 1. take an `FLStock` (North Sea plaice from `icesdata`); 2. fit a stock-recruitment relationship and build an `FLBRP` (`FLRebuild::eql`); 3. convert the `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; 4. build the SPiCT input from the `FLStock` catch and a biomass index, and fit; 5. check convergence and diagnostics; 6. compare the SPiCT estimates with the `FLBRP` equivalents; 7. coerce the fit to a `biodyn` object and project recovery at zero catch from $0.3B_{MSY}$. ```{r, pkgs} library(FLCore) library(FLBRP) library(mpb) library(spict) library(ggplot2) ``` `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. # Stock ```{r, stock} data(icesdata, package="icesdata") stk=icesdata[["ple.27.420"]] stk=window(stk, end=2020) summary(stk) ``` # Equilibrium reference points `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. ```{r, eql} 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. ```{r, quant} quant="eb" # or "ssb" ``` # Priors from the FLBRP 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)}$. ```{r, pellat} cv=c(r=0.3, K=0.5, n=0.3, psi=0.3) pr=spictPriors(stk, eq, quant, cv) pt =pr$pars psi=pr$psi idx=switch(quant, eb=FLRebuild::ebiomass(stk), ssb=ssb(stk)) pt pr$table ``` # SPiCT input 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. ```{r, inp} 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. ```{r, plotData, fig.height=4} plot(FLQuants(Catch=catch(stk)/scale, Index=idx/scale))+ labs(x="Year", y="Thousand tonnes") ``` # Fit ```{r, fit} fit=fit.spict(inp) fit=calc.process.resid(fit) fit ``` ## Convergence ```{r, convergence} 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). ```{r, plotTS, fig.height=6.5} plotSpictTS(fit, scale=scale) ``` ```{r, plotKobe, fig.height=5} plotSpictKobe(fit) ``` ```{r, plotProd, fig.height=4.5} plotSpictProduction(fit, scale=scale) ``` ```{r, plotFitIdx, fig.height=4} plotSpictFit(fit, scale=scale) ``` ## Prior and posterior ```{r, priors, fig.height=5} plotSpictPriors(fit) ``` ## Residuals One-step-ahead residuals should show no bias, autocorrelation or departure from normality. ```{r, osa, fig.height=4} fit=calc.osa.resid(fit) plotSpictRuns(fit) ``` ```{r, osaACF, fig.height=3.5} plotSpictACF(fit) ``` ```{r, osaQQ, fig.height=3.5} plotSpictQQ(fit) ``` ## Retrospective analysis Five one-year peels; each refit takes a few seconds. ```{r, retro, fig.height=5} fit=retro(fit, nretroyear=5) plotSpictRetro(fit) ``` ```{r, hindcast, fig.height=4} plotSpictHindcast(fit) ``` # SPiCT compared with the FLBRP `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. ```{r, gp} gp=function(nm, fit, i=1) spictPar(fit, nm, i) ``` 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. ```{r, compare} 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), ".")) ``` ```{r, fbar} c(FbarMSY=c(rf["msy", "harvest"])) ``` 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. ```{r, q} gp("logq", fit) ``` # Comparison figures 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. ```{r, spAnn} 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") ``` ## Trajectories 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. ```{r, trajectories, fig.height=6} 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") ``` ## Production functions 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. ```{r, production, fig.height=5} 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 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. ```{r, processError, fig.width=10, fig.height=8} 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"])) ``` # Coercing to biodyn `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. ```{r, spict2biodyn} bd=spict2biodyn(fit, scale) params(bd) mpb::refpts(bd) ``` `mpb::refpts` gives the same $F_{MSY}$, $B_{MSY}$ and MSY as the deterministic SPiCT reference points (to rounding). ```{r, plotBiodyn, fig.height=5} plot(bd) ``` # Recovery from 0.3BMSY at zero catch 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. ```{r, recover} 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])) tRecover ``` ```{r, plotRecover} ggplot(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. # Notes * The index here is an assessment output, so this is a consistency check between the age-structured and biomass dynamic models rather than an independent assessment. Survey indices can be added as further elements of `inp$obsI` and `inp$timeI`. * Setting `quant="ssb"` gives a much lower $n$ (a $B_{MSY}/K$ around 0.25 for plaice), since recruitment compensation acts on spawning biomass. * The FLBRP depends on the stock-recruitment relationship and its steepness prior in `FLRebuild::eql`; different choices give different priors. # Session information ```{r, session} sessionInfo() ``` # References 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.