In jabbamodel/ss3diags: What the Package Does (One Line, Title Case)
library(knitr)
hook_output = knit_hooks$get('output')
knit_hooks$set(output = function(x, options) {
# this hook is used only when the linewidth option is not NULL
if (!is.null(n <- options$linewidth)) {
x = knitr:::split_lines(x)
# any lines wider than n should be wrapped
if (any(nchar(x) > n)) x = strwrap(x, width = n)
x = paste(x, collapse = '\n')
}
hook_output(x, options)
})
knitr::opts_chunk$set(collapse = TRUE, comment = " " ,fig.align = 'center', cache=FALSE,tidy.opts=list(width.cutoff=80), tidy=TRUE)
Getting started {#s1}
This vignette introduces you to the ss3diags R package, which accompanies the paper “A cookbook for using model diagnostics in integrated stock assessments” by Carvalho, Winker et al. (2021).
The ss3diags comprises a explains set of functions for applying advanced model diagnostics for Stock Synthesis models. with the ss3diags R package for model. The ss3diag R package accompanies the paper “A cookbook for using model diagnostics in integrated stock assessments” by Carvalho, Winker et al. (2021). The ss3diags package builds on the widely used R package r4ss (Taylor et al. 2021), which is designed to support the use of the Stock Synthesis software modeling framework (Methot and Wetzel, 2013).
This vignette is divided into four sections. Section 1 consists of installing ss3diags and loading the data from the two case studies presented in the cookbook. Section 2 describes the plotting of various model diagnostics as described in the Cookbook. Section 3 provides a detailed explanation on how to assess model uncertainty using ss3diags. In Section 4 we provide a series “cookbook recipes” on how to implement selected model diagnostics in Stock Synthesis.
Installation
Both ss3diags and r4ss can be installed from gihtub using library(devtools):
installed.packages("devtools")
devtools::install_github("JABBAmodel/r4ss")
devtools::install_github("JABBAmodel/ss3diags")
library(r4ss)
library(ss3diags)
Loading built-in example data
The package contains two fully worked examples of Stock Synthesis assessments as presented in the Cookbook.
The Pacif hake (Merluccius productus) base case model
data("pac.hke")
Individual list objects generated using R package r4ss:
ss3phk: output from Stock Synthesis as read by r4ss::SS_output()
dir.path = "C:/Users/henni/Dropbox/ss3diags_demo/PacificHake"
ss3pke = r4ss::SS_output(file.path(dir.path,"Reference_Run"))
retro.phk: list of retrospective runs with r4ss:SS_doRetro() as read by r4ss::SSgetoutput()
The retrospective runs, produced with r4ss:SS_doRetro(), were located in the subfolders /Retro_Reference_Run and named "retro0","retro-1",..., "retro-7" and so on. To load those at once we use the r4ss::SSgetoutput(). Eight retrospective runs were conducted, where "retro0" corresponds to the "Reference_Run" and "retro-1" to "retro-7" are "peels". To assign the model names, specify start.retro = 0 and end.retro = 7 below.
start.retro = 0
end.retro = 7
retro.runs = "Retro_Reference_Run"
# load models
retro.phk <- r4ss::SSgetoutput(dirvec=file.path(dir.path,
retro.runs,paste0("retro",start.retro:-end.retro)))
The list object retro.phk[[1]] ("retro0") corresponds to the reference run ss3phk
aspm.phk: Comprises of two runs the "Reference_Run" and the "ASPM", which can be loaded together using r4ss::SSgetoutput() and then summarized with r4ss::SSsummarize()
asem.phk <- r4ss::SSgetoutput(dirvec=
file.path(dir.path,paste0("Reference_Run","APSM")))
asem.phk <- r4ss::SSsummarize(asem.phk)
North Atlantic shortfin mako (Isurus oxyrinchus)
data("natl.sma")
Individual list objects generated using R package r4ss:
-
ss3sma: output from Stock Synthesis as read by r4ss::SS_output()
-
retro.sma: list of retrospective runs with r4ss:SS_doRetro() as read by r4ss::SSgetoutput()
-
aspm.sma: list of Stock Synthesis referece and aspm created with r4ss::SSsummarize()
Model Diagnostics with ss3diags {#s2}
Residual diagnostics
The plotting options are kept mainly to those provided by r4ss.
Like with r4ss, if, for example, `SSplotRunstest() called with no further specifications several windows will open, which depends on the number abundance indices.
SSplotRunstest(ss3sma)
To visiualize the runs test for multiple indices, it is recommended to use of the function sspar(). The option plot add=TRUE in sspar() facilitates setting the graphic parameters so that they are suitable for ss3diags plots. The option add=TRUE prevents the plotting functions from over-writing sspar().
sspar(mfrow=c(3,2),plot.cex=0.8)
rt = SSplotRunstest(ss3sma,add=T,verbose=F)
It is also possible to select the indices that should be plotted. For example, we excluded CPUE2 as it was not fitted (zero weight to the likelihood). This creates space to add the joint-residual plot SSplotJABBAres to summarize all selected indices.
sspar(mfrow=c(3,2),plot.cex=0.8)
rt= SSplotRunstest(ss3sma,add=T,indexselect = c(1,3:6),legendcex = 0.8,verbose=F)
jr = SSplotJABBAres(ss3sma,add=T,indexselect = c(1,3:6),legendcex = 0.55,verbose=F)
The default for SSplotRunstest() and SSplotJABBAres() is plot the residual runs for the abundance indices. However, it is also possible to plot the composition data by specifying subplots="len" (or "age")
sspar(mfrow=c(3,2),plot.cex=0.8)
rt= SSplotRunstest(ss3sma,add=T,legendcex = 0.8,subplot="len",verbose=F)
jr = SSplotJABBAres(ss3sma,add=T,
legendcex = 0.55,legendloc="bottomright",subplot="len",verbose=F)
\newpage
To facilitate automated processing, results from several diagnostic tests can also be called without plotting
rti= SSrunstest(ss3sma,quant="cpue",verbose=F)
rtl= SSrunstest(ss3sma,quant="len",verbose=F)
rbind(rti,rtl)
The Pacific hake assessment provides an example of fits to age composition instead of length composition data, which can visualized by specifying subplots="age"
sspar(mfrow=c(2,2),plot.cex=0.8)
rti = SSplotRunstest(ss3phk,add=T,legendcex = 0.8,subplot="cpue",verbose=F)
rta = SSplotRunstest(ss3phk,add=T,legendcex = 0.8,subplot="age",verbose=F)
jra = SSplotJABBAres(ss3phk,add=T,legendcex = 0.7,subplot="age",verbose=F)
Retrospective and Forecast bias
In Stock Synthesis, retrospective analysis can be routinely using r4ss:SS_doRetro() (see Section 4.1). ss3diags provides the function SSplotRetro() to visualize the retrospective patterns of SBB and F and compute the associated Mohn's rho value (i.e. retrospective bias). This would require first loading the retrospective runs (Section 1.2]), which are already inbuilt into ss3diags in this case. The next step is to summarize the list of retrospective runs using r4ss::SSsummarize().
retroI.phk <- r4ss::SSsummarize(retro.phk,verbose=F)
We use notation "retroI" because r4ss::SSsummarize() summarizes the modeled quantities and abundance indices, but not length or age composition data. Using retroI.phk it is possible to produce some basic retrospective plots.
sspar(mfrow=c(1,2),plot.cex=0.8)
rb = SSplotRetro(retroI.phk,add=T,forecast = F,legend = F,verbose=F)
rf = SSplotRetro(retroI.phk,add=T,subplots="F", ylim=c(0,0.4),
forecast = F,legendloc="topleft",legendcex = 0.8,verbose=F)
An intuitive extension of the retrospective analysis is to assess potential forecast bias by adding the additional step of forward projecting quantities, such as SSB, over the truncated years. In Stock Synthesis the forecasts are automatically done when using r4ss:SS_doRetro().The forecasts are based on the settings specified in ‘forecast.ss’, which are also evoked when conducting future projections with the same model. The observed catches are used for the retrospective forecasts. Retrospective forecasts with Stock Synthesis are therefore only a matter of visualization, which can be done by setting the SSplotRetro() option forecast=TRUE.
sspar(mfrow=c(1,2),plot.cex=0.8)
rb = SSplotRetro(retroI.phk,add=T,forecast = T,legend = F,verbose=F,xmin=2000)
rf = SSplotRetro(retroI.phk,add=T,subplots="F", ylim=c(0,0.4),
forecast = T,legendloc="topleft",legendcex = 0.8,verbose=F,xmin=2000)
\newpage
The statistics from the retrospective analysis with forecasting, mohn's rho and forecast bias, can be called without plotting using the function SShcbias()
SShcbias(retroI.phk,quant="SSB",verbose=F)
SShcbias(retroI.phk,quant="F",verbose=F)
Hindcast Cross-Validation and prediction skill
Implementing the Hindcast Cross-Validation (HCxval) diagnostic in Stock Synthesis requires the same model outputs generated by r4ss:SS_doRetro() as described in Section 3.1. Therefore, no additional step is needed for HCxval if conducted in conjunction with retrospective analysis.
As a robust measure of prediction skill, we implemented the mean absolute
scaled error (MASE). In brief, the MASE score scales
the mean absolute error (MAE) of forecasts (i.e., prediction residuals) to
MAE of a naïve in-sample prediction, which is realized in the form of a simple 'persistence algorithm', i.e. tomorrow’s weather will be the same as today’s (see Eq. 3, p.5 in Carvalho and Winker et al. 2021). A MASE score > 1 indicates that the average model forecasts are worse than a random walk. Conversely, a MASE score of 0.5 indicates that the model forecasts twice as accurately as a naïve baseline prediction; thus, the model has prediction skill.
HCxval is implemented using function SSplotHCxval(), which produces the novel HCxval diagnostic plot and computes the MASE scores for CPUE indices, mean lengths or mean ages that have observations falling within the hindcast evaluation period.
Plotting HCxval for abundance indices requires the same step of summarizing the list of retrospective runs as for the retrospective analysis, which therefore only needs be done once. Below is a summary of the retrospective runs for shortfin mako.
retroI.sma <- r4ss::SSsummarize(retro.sma,verbose=F)
sspar(mfrow=c(3,2),plot.cex=0.8)
hci = SSplotHCxval(retroI.sma,add=T,verbose=F,ylimAdj = 1.3,legendcex = 0.7)
\newpage
The forecast length- and age-composition are located in the Stock Synthesis report.sso as "ghost files". To extract and summarize the composition data in the form of observed and expected mean lengths and age ss3diags provides the function SSretroComps().
retroC.sma = SSretroComps(retro.sma)
sspar(mfrow=c(1,2),plot.cex=0.8)
hcl = SSplotHCxval(retroC.sma,subplots="len",add=T,verbose=F,
ylimAdj = 1.3,legendcex = 0.7,indexselect = c(1,2))
The Figure above provides some additional, so called adjusted MASE values, in brackets. This gets invoked in cases where the inter-annual variation in the observed values is very small (default MAE < 0.1 for naive predictions log(y[t+1])-log(y[t])).
The reasoning is that prediction residuals must be already very accurate to fall below this threshold. The adjusted MASE essential keep the naive prediction MAE denominator of the MASE to a maximum. Below we show the effect of changing adjustment threshold from the default MAE.base.adj = 0.1
mase1 = SSmase(retroC.sma,quant="len",MAE.base.adj = 0.1,indexselect = c(1:2))
mase1
to a larger value MAE.base.adj = 0.15
SSmase(retroC.sma,quant="len",MAE.base.adj = 0.15,indexselect = c(1:2))
where MASE is the ratio of the mean absolute error of the prediction residuals MAE.PR to the residuals of the naive predictions MAE.base
mase1$MAE.PR/mase1$MAE.base
mase1$MASE
and MASE.adj
mase1$MAE.PR/pmax(mase1$MAE.base,0.1)
mase1$MASE.adj
Applying HCxval for composition data requires correctly specifying the composition data type fitted in the model. For example, age composition data need to be specified as "age" in SSplotHCxval and SSmase, as shown below for the Pacific hake model.
retroC.phk = SSretroComps(retro.phk) # summarize comps
sspar(mfrow=c(1,2),plot.cex=0.8)
hcl = SSplotHCxval(retroC.phk,subplots="age",add=T,
verbose=F,ylimAdj = 1.3,legendcex = 0.7,indexselect = c(1,2))
SSmase(retroC.phk,quants="age")
\newpage
Model uncertainty {#s3}
The management advice frameworks increasingly require translating the estimated uncertainty about the stock status into probabilistic statements (Kell et al. 2016). A classical example is the Kobe framework used in tuna Regional Fisheries Management Organisations (tRFMOs) around the world. The key quantities of interest are typically the ratios $SSB/SSB_{MSY}$ and $F/F_{MSY}$. It is reasonably straight forward in Stock Synthesis to approximate uncertainty of individual quantities (e.g. $SSB$) from the asymptotic standard errors (SE) derived from the Hessian matrix using the delta method. However, the joint distribution of $SSB/SSB_{MSY}$ and $F/F_{MSY}$
requires to adequately account for the covariance structure between these two derived quantities. Joint distributions were typically constructed using bootstrap or Markov Chain Monte-Carlo (MCMC) methods. However, these methods can be computationally intense and time-consuming in integrated assessments.
As an alternative, ss3diags implements a rapid delta-Multivariate lognormal approximation with SSdeltaMVLN() to generate joint error distributions for $SSB/SSB_{ref}$ and $F/F_{ref}$,where the $ref$ may refer to $MSY$, but also other reference points (e.g., $SSB_{40}$ and $F_{40}$). In Stock Synthesis, these ratios are determined by the derived quantities Bratio and F, where either can take the form of ratios (e.g. $F/F_{ref}$) or absolute value (e.g. absF) depending on settings in the starter.ss file.
Let Bratio be $u = SSB/SSB_{ref}$, F be $v = F$, and $w = F_{ref}$ be the F reference point of interest (e.g. $F_{MSY}$), with $x = \log(u)$, $y = \log(v)$ and $z = \log(w)$, then the variance-covariance matrix $VCM$ has the form:
$$
VCM_{x,y,z} =
\begin{pmatrix}
\sigma^2_{x} & cov_{x,y} & cov_{x,y} \
cov_{x,y} & \sigma^2_{y} & cov_{y,z} \
cov_{x,z} & cov_{y,z} & \sigma^2_{z}
\end{pmatrix}
$$
where, e.g., $\sigma^2_{x}$ is the variance of $x$ and $cov_{x,y}$ is the covariance of $x$ and $y$. Deriving those requires conducting a few normal to lognormal transformations. First, the variances are approximated as:
$$
\sigma^2_{x} = \log\left(1+\left(\frac{SE_u}{u}\right)^2\right)
$$
where $SE_{u}$, $SE_{v}$ and $SE_{z}$ are the asymptotic standard error estimates for $u = SSB/SSB_{ref}$, $v = F$ and $z = F_{ref}$.
The corresponding covariance for $x$ and $y$, can then be approximated on the log-scale by:
$$
COV_{x,y} = \log \left( {1+\rho_{u,v} \sqrt{\sigma^2_{x}\sigma^2_{y}}} \right)
$$
where $rho_{u,v}$ denotes the correlation of $u$ and $v$.
To generate a joint distribution of $\tilde{u}$ = $SSB/SSB_{ref}$, $\tilde{v}$ = $F$ and $\tilde{z}$ = $F_{ref}$, a multivariate random generator is used, which is available in the R package ‘mvtnorm’, to obtain a large number (e.g. nsim = 10,000) iterations, such that
$$
JD(\tilde{u},\tilde{v},\tilde{w}) = \exp(MVN(\mu_{x,y,z},VCM_{x,y,z}))
$$
so that
$$
\tilde{SSB}/\tilde{SSB}_{{MSY}} = \tilde{u}
$$
and
$$
\tilde{F}/\tilde{F}_{{MSY}} = \tilde{v}/\tilde{w}
$$
The reference points depend on the settings in the starter.ss file that determine the derived quantities Bratio and Fvalue.
In the first example, we consider the ss3sma Stock Synthesis model, which was run with starter.ss settings that are common to produce target Kobe plot estimates of $SSB/SSB_{MSY}$ and $F/F_{MSY}$ in tRFMO assessments:
2 # Depletion basis: 2=rel SPBmsy; 3=rel X*SPB_styr; 4=rel X*SPB_endyr
2 # F_report_basis: 0=raw_F_report; 1=F/Fspr; 2=F/Fmsy ; 3=F/Fbtgt
To generate a joint MVLN of $SSB/SSB_{MSY}$ and $F/F_{MSY}$, the default options can be used.
```r$ and $F/F_{MSY}$ for North Atlantic shortfin mako Stock model."}
mvln = SSdeltaMVLN(ss3sma,run="SMA")
We provide the function `SSsettingsBratioF(ss3sma)` to the `starter.ss` settings:
```r
SSsettingsBratioF(ss3sma)
This function is also inbuilt in the SSdeltaMVLN() to prevent misleading results. The SSdeltaMVLN() output include the maximum likelihood estimates (mles) and the MVLN monte-Carlo distributions $kb of $SSB/SSB_{MSY}$, $F/F_{MSY}$ and $F$. Note the additional quantities $SSB$ and $Rec$ are generated independently from lognormal distributions for practical reasons. These can be plotted by
```r$, $F/F_{MSY}$, $SSB$, $F$ ,Recruitment and Catch trajectories for the North Atlantic shortfin mako SS3 model"}
sspar(mfrow=c(3,2),plot.cex = 0.7)
SSplotEnsemble(mvln$kb,ylabs=mvln$labels,add=T,verbose=F)
The `SSdeltaMVLN()` provides the option to set alternative `Fref` values, but this is only possible for the recommended `starter.ss` option 0 for `F_report_basis`. For option 2, `SSdeltaMVLN()` prompts an error if `Fref` is changed.
By comparison, the Pacific Hake base case model `ss3phk` is run with settings that are common in NOAA assessments, with `Bratio` set $SSB/SSB_{0}$ and `F` is typically kept at absolute quantity.
`1 # Depletion basis: 2=rel SPBmsy; 3=rel X*SPB_styr; 4=rel X*SPB_endyr`
`0 # F_report_basis: 0=raw_F_report; 1=F/Fspr; 2=F/Fmsy ; 3=F/Fbtgt`
The management quantities in this case are $SSB/SSB_{40}$ and $F/F_{spr40}$, where the target of 40\% is specified in the `forecast.ss` file.
`0.4 # SPR target (e.g. 0.40)`
`0.4 # Biomass target (e.g. 0.40)`
```r$ and $F/F_{SPR40}$ for North Atlantic shortfin mako Stock model."}
mvln = SSdeltaMVLN(ss3phk,run="Pac.Hake", Fref="SPR",plot=TRUE)
```r$, $F/F_{SPR40}$, $SSB$, $F$, Recruitment and Catch trajectories for the Pacific Hake SS3 model"}
sspar(mfrow=c(3,2),plot.cex = 0.7)
SSplotEnsemble(mvln$kb,ylabs=mvln$labels,add=T,verbose=F)
<br>
It is important to note that MVLN approximation differs notably from MCMC posterior for this model as documented in [Stewart et al. (2013)](https://www.sciencedirect.com/science/article/abs/pii/S0165783612002081) and [Taylor et al. (2021)](https://www.sciencedirect.com/science/article/abs/pii/S0165783621000527). Such differences may be more likely to occur in cases where key parameters such as steepness $h$ or natural $M$ are estimated using informative priors, which can result in left skewed (non-lognormal) distributions of the benchmarks $F_{ref}$ and $B_{ref}$.
In other instances, mismatches between the`SSdeltaMVLN` and MCMC may also be caused by the latter's poor performance due to poor regularization [(Monnahan et al., 2019)](https://academic.oup.com/icesjms/article/76/6/1477/5475859?login=true).
To facilitate a comparison between the `SSdeltaMLVN()` and MCMC outputs, we provide the function `SSdiagsMCMC()`, which is illustrated on the example of the Stock Synthesis model for the ICES Gulf of Bothian Herring stock [(ICES, 2021)](https://www.ices.dk/sites/pub/Publication%20Reports/Stock%20Annexes/2021/her.27.3031_SA.pdf).
`SSdiagsMCMC()` requires loading both the `report.sso` and MCMC output in the `posterior.sso` file, where the MCMC was in this case run in the subfolder of the assessment file `/mcmc`.
```r
# Alread loaded to ss3diags
ss3her = SS_output("gob_her")
mcmc.her = SSreadMCMC("gob_her/mcmc")
The options and output of SSdiagsMCMC() are largely identical to SSdeltaMVLN. In this case, the starter.ss is the same as for ss3phk, only we use Fref = "Btgt" for illustration with $SB_{50}$ and $F_{SB50}$ for this illustration.
mvln = SSdeltaMVLN(ss3her,Fref="Btgt",plot=F,run="mvln")
mcmc = SSdiagsMCMC(mcmc.her,ss3her,Fref="Btgt",plot=F,run="mcmc")
Comparing delta-MVLN with MCMC simply requires combining the $kb outputs by, e.g.,
```r$ and $F/F_{50}$ for the ICES Gulf of Bothia Herring SS3 model"}
sspar(mfrow=c(1,1),plot.cex = 0.8)
SSplotKobe(rbind(mvln$kb,mcmc$kb),joint=F,
xlab=mvln$labels[1],ylab=mvln$labels[2],fill=F)
```r$, $F/F_{SB40}$, $SSB$ and $F$ for the ICES Gulf of Bothia Herring SS3 model"}
sspar(mfrow=c(2,2),plot.cex = 0.7)
SSplotEnsemble(rbind(mvln$kb,mcmc$kb),ylabs=mvln$labels,add=T,subplots = c("stock", "harvest", "SSB", "F"),verbose=F)
This works equally for joining a model ensemble.
\newpage
Cookbook Recipies {#s4}
Retrospectives with hindcasts {#r1}
Retrospective analysis can be run for Stock Synthesis using the function r4ss::SS_doRetro() available in r4ss. This setup of the retrospective analysis has the advantage that forecasts are conducted automatically given the catch. This makes it possible to apply retrospective forecasting and hindcast cross-validations of observations based on the same output.
library(r4ss)
Below is a step-by-step cookbook recipe for retrospective analysis in Stock Synthesis.
Step1: Identify restrospective period
Specify the range pf peels that will then determine the end.yr.vec of runs in r4ss::SS_doRetro()
start.retro <- 0 # end year of reference year
end.retro <- 7 # number of years for retrospective e.g.,
Step 2: Identify the base directory
Specify the path directory that holds the folder with the base case run. In this case the Pacific Hake folder with a model folder `Reference_Run"
dirname.base = "C:/Users/henni/Dropbox/ss3diags_demo/PacificHake"
run = "Reference_Run"
model.run <- file.path(dirname.base,run)
model.run
Step 3: DAT and CONTROL files
Specify the names of the data and control files.
Note these files are named differently from the DATA.ss and CONTROL.ss. In this case
DAT = "phk.dat"
CTL = "phk.ctl"
The names of the DAT and CONTROL are declared on the top of the 'starter.ss`, e.g.
#C Hake starter file
phk.dat
phk.ctl
Step 4: Create a subdirectory for the Retrospectives
There are several ways to organize the retrospective output structure.
First create a new subfolder for the retrospective runs output
dir.retro <- paste0(dirname.base,'/Retro_',run)
dir.create(path=dir.retro, showWarnings = F)
Also create a subdirectory for the retrospective model folders
dir.create(path=file.path(dir.retro,"retros"), showWarnings = F)
then copy model run files to the new retrospective folder
file.copy(file.path(model.run,"starter.ss_new"), file.path(dir.retro,"starter.ss"))
file.copy(file.path(model.run,"control.ss_new"),
file.path(dir.retro, CTL))
file.copy(file.path(model.run,"data.ss_new"),
file.path(dir.retro, DAT))
file.copy(file.path(model.run,"forecast.ss"),
file.path(dir.retro, "forecast.ss"))
file.copy(file.path(model.run,"SS.exe"),
file.path(dir.retro, "SS.exe"))
# Automatically ignored for models without wtatage.ss
file.copy(file.path(model.run,"wtatage.ss"),
file.path(dir.retro, "wtatage.ss"))
Step 5: Modify Starter.ss file
Modifying the Starter File helps to speed up model runs
starter <- readLines(paste0(dir.retro, "/starter.ss"))
#[8] "2 # run display detail (0,1,2)"
linen <- grep("# run display detail", starter)
starter[linen] <- paste0( 1 , " # run display detail (0,1,2)" )
# write modified starter.ss
write(starter, file.path(dir.retro, "starter.ss"))
Step 6: Execute retrospective runs
Run the retrospective analyses using r4SS function r4ss::SS_doRetro
Ideally, the runs should be done with Hessian to evaluate the retrospective trajectories with respect to the confidence interval coverage of the reference model.
r4ss::SS_doRetro(masterdir=dir.retro, oldsubdir="",
newsubdir="", years=start.retro:-end.retro)
However, for larger models it might be desirable to shorten run times, by not inverting the hessian, using the option extras = "-nohess" (much faster)
r4ss::SS_doRetro(masterdir=dir.retro,oldsubdir="",
newsubdir="", years=start.retro:-end.retro,extras = "-nohess")
Step 7: Read SS_doRetro() output
retro.phk <- r4ss::SSgetoutput(dirvec=file.path(dir.retro,
paste0("retro",start.retro:-end.retro)))
It is often useful to save the retro model runs as .rdata file for further processing with ss3diags, considering that reading the models with r4ss::SSgetoutput() can be quite time-consuming for more complex models.
save(retro.phk, file=file.path(dir.retro,"retro.phk.rdata"))
Step 8: Check
library(ss3diags)
check.retro =r4ss::SSsummarize(retro.phk)
sspar(mfrow=c(1,1))
SSplotRetro(check.retro,forecast = T,add=T,legendcex = 0.7,legendloc = "topleft")
R0 profiling {#r2}
ASPM diagnostic {#r3}
Jittering {#r4}
jabbamodel/ss3diags documentation built on Oct. 8, 2024, 11:49 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(knitr) hook_output = knit_hooks$get('output') knit_hooks$set(output = function(x, options) { # this hook is used only when the linewidth option is not NULL if (!is.null(n <- options$linewidth)) { x = knitr:::split_lines(x) # any lines wider than n should be wrapped if (any(nchar(x) > n)) x = strwrap(x, width = n) x = paste(x, collapse = '\n') } hook_output(x, options) })
knitr::opts_chunk$set(collapse = TRUE, comment = " " ,fig.align = 'center', cache=FALSE,tidy.opts=list(width.cutoff=80), tidy=TRUE)
Getting started {#s1}
This vignette introduces you to the ss3diags R package, which accompanies the paper “A cookbook for using model diagnostics in integrated stock assessments” by Carvalho, Winker et al. (2021).
The ss3diags comprises a explains set of functions for applying advanced model diagnostics for Stock Synthesis models. with the ss3diags R package for model. The ss3diag R package accompanies the paper “A cookbook for using model diagnostics in integrated stock assessments” by Carvalho, Winker et al. (2021). The ss3diags package builds on the widely used R package r4ss (Taylor et al. 2021), which is designed to support the use of the Stock Synthesis software modeling framework (Methot and Wetzel, 2013).
This vignette is divided into four sections. Section 1 consists of installing ss3diags and loading the data from the two case studies presented in the cookbook. Section 2 describes the plotting of various model diagnostics as described in the Cookbook. Section 3 provides a detailed explanation on how to assess model uncertainty using ss3diags. In Section 4 we provide a series “cookbook recipes” on how to implement selected model diagnostics in Stock Synthesis.
Installation
Both ss3diags and r4ss can be installed from gihtub using library(devtools):
installed.packages("devtools") devtools::install_github("JABBAmodel/r4ss") devtools::install_github("JABBAmodel/ss3diags")
library(r4ss) library(ss3diags)
Loading built-in example data
The package contains two fully worked examples of Stock Synthesis assessments as presented in the Cookbook.
The Pacif hake (Merluccius productus) base case model
data("pac.hke")
Individual list objects generated using R package r4ss:
ss3phk: output from Stock Synthesis as read byr4ss::SS_output()
dir.path = "C:/Users/henni/Dropbox/ss3diags_demo/PacificHake" ss3pke = r4ss::SS_output(file.path(dir.path,"Reference_Run"))
retro.phk: list of retrospective runs withr4ss:SS_doRetro()as read byr4ss::SSgetoutput()
The retrospective runs, produced with r4ss:SS_doRetro(), were located in the subfolders /Retro_Reference_Run and named "retro0","retro-1",..., "retro-7" and so on. To load those at once we use the r4ss::SSgetoutput(). Eight retrospective runs were conducted, where "retro0" corresponds to the "Reference_Run" and "retro-1" to "retro-7" are "peels". To assign the model names, specify start.retro = 0 and end.retro = 7 below.
start.retro = 0 end.retro = 7 retro.runs = "Retro_Reference_Run" # load models retro.phk <- r4ss::SSgetoutput(dirvec=file.path(dir.path, retro.runs,paste0("retro",start.retro:-end.retro)))
The list object retro.phk[[1]] ("retro0") corresponds to the reference run ss3phk
aspm.phk: Comprises of two runs the "Reference_Run" and the "ASPM", which can be loaded together usingr4ss::SSgetoutput()and then summarized withr4ss::SSsummarize()
asem.phk <- r4ss::SSgetoutput(dirvec= file.path(dir.path,paste0("Reference_Run","APSM"))) asem.phk <- r4ss::SSsummarize(asem.phk)
North Atlantic shortfin mako (Isurus oxyrinchus)
data("natl.sma")
Individual list objects generated using R package r4ss:
-
ss3sma: output from Stock Synthesis as read byr4ss::SS_output() -
retro.sma: list of retrospective runs withr4ss:SS_doRetro()as read byr4ss::SSgetoutput() -
aspm.sma: list of Stock Synthesis referece and aspm created withr4ss::SSsummarize()
Model Diagnostics with ss3diags {#s2}
Residual diagnostics
The plotting options are kept mainly to those provided by r4ss. Like with r4ss, if, for example, `SSplotRunstest() called with no further specifications several windows will open, which depends on the number abundance indices.
SSplotRunstest(ss3sma)
To visiualize the runs test for multiple indices, it is recommended to use of the function sspar(). The option plot add=TRUE in sspar() facilitates setting the graphic parameters so that they are suitable for ss3diags plots. The option add=TRUE prevents the plotting functions from over-writing sspar().
sspar(mfrow=c(3,2),plot.cex=0.8) rt = SSplotRunstest(ss3sma,add=T,verbose=F)
It is also possible to select the indices that should be plotted. For example, we excluded CPUE2 as it was not fitted (zero weight to the likelihood). This creates space to add the joint-residual plot SSplotJABBAres to summarize all selected indices.
sspar(mfrow=c(3,2),plot.cex=0.8) rt= SSplotRunstest(ss3sma,add=T,indexselect = c(1,3:6),legendcex = 0.8,verbose=F) jr = SSplotJABBAres(ss3sma,add=T,indexselect = c(1,3:6),legendcex = 0.55,verbose=F)
The default for SSplotRunstest() and SSplotJABBAres() is plot the residual runs for the abundance indices. However, it is also possible to plot the composition data by specifying subplots="len" (or "age")
sspar(mfrow=c(3,2),plot.cex=0.8) rt= SSplotRunstest(ss3sma,add=T,legendcex = 0.8,subplot="len",verbose=F) jr = SSplotJABBAres(ss3sma,add=T, legendcex = 0.55,legendloc="bottomright",subplot="len",verbose=F)
\newpage
To facilitate automated processing, results from several diagnostic tests can also be called without plotting
rti= SSrunstest(ss3sma,quant="cpue",verbose=F) rtl= SSrunstest(ss3sma,quant="len",verbose=F) rbind(rti,rtl)
The Pacific hake assessment provides an example of fits to age composition instead of length composition data, which can visualized by specifying subplots="age"
sspar(mfrow=c(2,2),plot.cex=0.8) rti = SSplotRunstest(ss3phk,add=T,legendcex = 0.8,subplot="cpue",verbose=F) rta = SSplotRunstest(ss3phk,add=T,legendcex = 0.8,subplot="age",verbose=F) jra = SSplotJABBAres(ss3phk,add=T,legendcex = 0.7,subplot="age",verbose=F)
Retrospective and Forecast bias
In Stock Synthesis, retrospective analysis can be routinely using r4ss:SS_doRetro() (see Section 4.1). ss3diags provides the function SSplotRetro() to visualize the retrospective patterns of SBB and F and compute the associated Mohn's rho value (i.e. retrospective bias). This would require first loading the retrospective runs (Section 1.2]), which are already inbuilt into ss3diags in this case. The next step is to summarize the list of retrospective runs using r4ss::SSsummarize().
retroI.phk <- r4ss::SSsummarize(retro.phk,verbose=F)
We use notation "retroI" because r4ss::SSsummarize() summarizes the modeled quantities and abundance indices, but not length or age composition data. Using retroI.phk it is possible to produce some basic retrospective plots.
sspar(mfrow=c(1,2),plot.cex=0.8) rb = SSplotRetro(retroI.phk,add=T,forecast = F,legend = F,verbose=F) rf = SSplotRetro(retroI.phk,add=T,subplots="F", ylim=c(0,0.4), forecast = F,legendloc="topleft",legendcex = 0.8,verbose=F)
An intuitive extension of the retrospective analysis is to assess potential forecast bias by adding the additional step of forward projecting quantities, such as SSB, over the truncated years. In Stock Synthesis the forecasts are automatically done when using r4ss:SS_doRetro().The forecasts are based on the settings specified in ‘forecast.ss’, which are also evoked when conducting future projections with the same model. The observed catches are used for the retrospective forecasts. Retrospective forecasts with Stock Synthesis are therefore only a matter of visualization, which can be done by setting the SSplotRetro() option forecast=TRUE.
sspar(mfrow=c(1,2),plot.cex=0.8) rb = SSplotRetro(retroI.phk,add=T,forecast = T,legend = F,verbose=F,xmin=2000) rf = SSplotRetro(retroI.phk,add=T,subplots="F", ylim=c(0,0.4), forecast = T,legendloc="topleft",legendcex = 0.8,verbose=F,xmin=2000)
\newpage
The statistics from the retrospective analysis with forecasting, mohn's rho and forecast bias, can be called without plotting using the function SShcbias()
SShcbias(retroI.phk,quant="SSB",verbose=F) SShcbias(retroI.phk,quant="F",verbose=F)
Hindcast Cross-Validation and prediction skill
Implementing the Hindcast Cross-Validation (HCxval) diagnostic in Stock Synthesis requires the same model outputs generated by r4ss:SS_doRetro() as described in Section 3.1. Therefore, no additional step is needed for HCxval if conducted in conjunction with retrospective analysis.
As a robust measure of prediction skill, we implemented the mean absolute
scaled error (MASE). In brief, the MASE score scales
the mean absolute error (MAE) of forecasts (i.e., prediction residuals) to
MAE of a naïve in-sample prediction, which is realized in the form of a simple 'persistence algorithm', i.e. tomorrow’s weather will be the same as today’s (see Eq. 3, p.5 in Carvalho and Winker et al. 2021). A MASE score > 1 indicates that the average model forecasts are worse than a random walk. Conversely, a MASE score of 0.5 indicates that the model forecasts twice as accurately as a naïve baseline prediction; thus, the model has prediction skill.
HCxval is implemented using function SSplotHCxval(), which produces the novel HCxval diagnostic plot and computes the MASE scores for CPUE indices, mean lengths or mean ages that have observations falling within the hindcast evaluation period.
Plotting HCxval for abundance indices requires the same step of summarizing the list of retrospective runs as for the retrospective analysis, which therefore only needs be done once. Below is a summary of the retrospective runs for shortfin mako.
retroI.sma <- r4ss::SSsummarize(retro.sma,verbose=F)
sspar(mfrow=c(3,2),plot.cex=0.8) hci = SSplotHCxval(retroI.sma,add=T,verbose=F,ylimAdj = 1.3,legendcex = 0.7)
\newpage
The forecast length- and age-composition are located in the Stock Synthesis report.sso as "ghost files". To extract and summarize the composition data in the form of observed and expected mean lengths and age ss3diags provides the function SSretroComps().
retroC.sma = SSretroComps(retro.sma)
sspar(mfrow=c(1,2),plot.cex=0.8) hcl = SSplotHCxval(retroC.sma,subplots="len",add=T,verbose=F, ylimAdj = 1.3,legendcex = 0.7,indexselect = c(1,2))
The Figure above provides some additional, so called adjusted MASE values, in brackets. This gets invoked in cases where the inter-annual variation in the observed values is very small (default MAE < 0.1 for naive predictions log(y[t+1])-log(y[t])).
The reasoning is that prediction residuals must be already very accurate to fall below this threshold. The adjusted MASE essential keep the naive prediction MAE denominator of the MASE to a maximum. Below we show the effect of changing adjustment threshold from the default MAE.base.adj = 0.1
mase1 = SSmase(retroC.sma,quant="len",MAE.base.adj = 0.1,indexselect = c(1:2)) mase1
to a larger value MAE.base.adj = 0.15
SSmase(retroC.sma,quant="len",MAE.base.adj = 0.15,indexselect = c(1:2))
where MASE is the ratio of the mean absolute error of the prediction residuals MAE.PR to the residuals of the naive predictions MAE.base
mase1$MAE.PR/mase1$MAE.base mase1$MASE
and MASE.adj
mase1$MAE.PR/pmax(mase1$MAE.base,0.1) mase1$MASE.adj
Applying HCxval for composition data requires correctly specifying the composition data type fitted in the model. For example, age composition data need to be specified as "age" in SSplotHCxval and SSmase, as shown below for the Pacific hake model.
retroC.phk = SSretroComps(retro.phk) # summarize comps sspar(mfrow=c(1,2),plot.cex=0.8) hcl = SSplotHCxval(retroC.phk,subplots="age",add=T, verbose=F,ylimAdj = 1.3,legendcex = 0.7,indexselect = c(1,2)) SSmase(retroC.phk,quants="age")
\newpage
Model uncertainty {#s3}
The management advice frameworks increasingly require translating the estimated uncertainty about the stock status into probabilistic statements (Kell et al. 2016). A classical example is the Kobe framework used in tuna Regional Fisheries Management Organisations (tRFMOs) around the world. The key quantities of interest are typically the ratios $SSB/SSB_{MSY}$ and $F/F_{MSY}$. It is reasonably straight forward in Stock Synthesis to approximate uncertainty of individual quantities (e.g. $SSB$) from the asymptotic standard errors (SE) derived from the Hessian matrix using the delta method. However, the joint distribution of $SSB/SSB_{MSY}$ and $F/F_{MSY}$
requires to adequately account for the covariance structure between these two derived quantities. Joint distributions were typically constructed using bootstrap or Markov Chain Monte-Carlo (MCMC) methods. However, these methods can be computationally intense and time-consuming in integrated assessments.
As an alternative, ss3diags implements a rapid delta-Multivariate lognormal approximation with SSdeltaMVLN() to generate joint error distributions for $SSB/SSB_{ref}$ and $F/F_{ref}$,where the $ref$ may refer to $MSY$, but also other reference points (e.g., $SSB_{40}$ and $F_{40}$). In Stock Synthesis, these ratios are determined by the derived quantities Bratio and F, where either can take the form of ratios (e.g. $F/F_{ref}$) or absolute value (e.g. absF) depending on settings in the starter.ss file.
Let Bratio be $u = SSB/SSB_{ref}$, F be $v = F$, and $w = F_{ref}$ be the F reference point of interest (e.g. $F_{MSY}$), with $x = \log(u)$, $y = \log(v)$ and $z = \log(w)$, then the variance-covariance matrix $VCM$ has the form:
$$ VCM_{x,y,z} = \begin{pmatrix} \sigma^2_{x} & cov_{x,y} & cov_{x,y} \ cov_{x,y} & \sigma^2_{y} & cov_{y,z} \ cov_{x,z} & cov_{y,z} & \sigma^2_{z} \end{pmatrix} $$
where, e.g., $\sigma^2_{x}$ is the variance of $x$ and $cov_{x,y}$ is the covariance of $x$ and $y$. Deriving those requires conducting a few normal to lognormal transformations. First, the variances are approximated as:
$$
\sigma^2_{x} = \log\left(1+\left(\frac{SE_u}{u}\right)^2\right)
$$
where $SE_{u}$, $SE_{v}$ and $SE_{z}$ are the asymptotic standard error estimates for $u = SSB/SSB_{ref}$, $v = F$ and $z = F_{ref}$.
The corresponding covariance for $x$ and $y$, can then be approximated on the log-scale by:
$$ COV_{x,y} = \log \left( {1+\rho_{u,v} \sqrt{\sigma^2_{x}\sigma^2_{y}}} \right) $$
where $rho_{u,v}$ denotes the correlation of $u$ and $v$.
To generate a joint distribution of $\tilde{u}$ = $SSB/SSB_{ref}$, $\tilde{v}$ = $F$ and $\tilde{z}$ = $F_{ref}$, a multivariate random generator is used, which is available in the R package ‘mvtnorm’, to obtain a large number (e.g. nsim = 10,000) iterations, such that
$$ JD(\tilde{u},\tilde{v},\tilde{w}) = \exp(MVN(\mu_{x,y,z},VCM_{x,y,z})) $$ so that
$$
\tilde{SSB}/\tilde{SSB}_{{MSY}} = \tilde{u}
$$
and
$$
\tilde{F}/\tilde{F}_{{MSY}} = \tilde{v}/\tilde{w}
$$
The reference points depend on the settings in the starter.ss file that determine the derived quantities Bratio and Fvalue.
In the first example, we consider the ss3sma Stock Synthesis model, which was run with starter.ss settings that are common to produce target Kobe plot estimates of $SSB/SSB_{MSY}$ and $F/F_{MSY}$ in tRFMO assessments:
2 # Depletion basis: 2=rel SPBmsy; 3=rel X*SPB_styr; 4=rel X*SPB_endyr
2 # F_report_basis: 0=raw_F_report; 1=F/Fspr; 2=F/Fmsy ; 3=F/Fbtgt
To generate a joint MVLN of $SSB/SSB_{MSY}$ and $F/F_{MSY}$, the default options can be used.
```r$ and $F/F_{MSY}$ for North Atlantic shortfin mako Stock model."}
mvln = SSdeltaMVLN(ss3sma,run="SMA")
We provide the function `SSsettingsBratioF(ss3sma)` to the `starter.ss` settings: ```r SSsettingsBratioF(ss3sma)
This function is also inbuilt in the SSdeltaMVLN() to prevent misleading results. The SSdeltaMVLN() output include the maximum likelihood estimates (mles) and the MVLN monte-Carlo distributions $kb of $SSB/SSB_{MSY}$, $F/F_{MSY}$ and $F$. Note the additional quantities $SSB$ and $Rec$ are generated independently from lognormal distributions for practical reasons. These can be plotted by
```r$, $F/F_{MSY}$, $SSB$, $F$ ,Recruitment and Catch trajectories for the North Atlantic shortfin mako SS3 model"} sspar(mfrow=c(3,2),plot.cex = 0.7) SSplotEnsemble(mvln$kb,ylabs=mvln$labels,add=T,verbose=F)
The `SSdeltaMVLN()` provides the option to set alternative `Fref` values, but this is only possible for the recommended `starter.ss` option 0 for `F_report_basis`. For option 2, `SSdeltaMVLN()` prompts an error if `Fref` is changed. By comparison, the Pacific Hake base case model `ss3phk` is run with settings that are common in NOAA assessments, with `Bratio` set $SSB/SSB_{0}$ and `F` is typically kept at absolute quantity. `1 # Depletion basis: 2=rel SPBmsy; 3=rel X*SPB_styr; 4=rel X*SPB_endyr` `0 # F_report_basis: 0=raw_F_report; 1=F/Fspr; 2=F/Fmsy ; 3=F/Fbtgt` The management quantities in this case are $SSB/SSB_{40}$ and $F/F_{spr40}$, where the target of 40\% is specified in the `forecast.ss` file. `0.4 # SPR target (e.g. 0.40)` `0.4 # Biomass target (e.g. 0.40)` ```r$ and $F/F_{SPR40}$ for North Atlantic shortfin mako Stock model."} mvln = SSdeltaMVLN(ss3phk,run="Pac.Hake", Fref="SPR",plot=TRUE)
```r$, $F/F_{SPR40}$, $SSB$, $F$, Recruitment and Catch trajectories for the Pacific Hake SS3 model"} sspar(mfrow=c(3,2),plot.cex = 0.7) SSplotEnsemble(mvln$kb,ylabs=mvln$labels,add=T,verbose=F)
<br>
It is important to note that MVLN approximation differs notably from MCMC posterior for this model as documented in [Stewart et al. (2013)](https://www.sciencedirect.com/science/article/abs/pii/S0165783612002081) and [Taylor et al. (2021)](https://www.sciencedirect.com/science/article/abs/pii/S0165783621000527). Such differences may be more likely to occur in cases where key parameters such as steepness $h$ or natural $M$ are estimated using informative priors, which can result in left skewed (non-lognormal) distributions of the benchmarks $F_{ref}$ and $B_{ref}$.
In other instances, mismatches between the`SSdeltaMVLN` and MCMC may also be caused by the latter's poor performance due to poor regularization [(Monnahan et al., 2019)](https://academic.oup.com/icesjms/article/76/6/1477/5475859?login=true).
To facilitate a comparison between the `SSdeltaMLVN()` and MCMC outputs, we provide the function `SSdiagsMCMC()`, which is illustrated on the example of the Stock Synthesis model for the ICES Gulf of Bothian Herring stock [(ICES, 2021)](https://www.ices.dk/sites/pub/Publication%20Reports/Stock%20Annexes/2021/her.27.3031_SA.pdf).
`SSdiagsMCMC()` requires loading both the `report.sso` and MCMC output in the `posterior.sso` file, where the MCMC was in this case run in the subfolder of the assessment file `/mcmc`.
```r
# Alread loaded to ss3diags
ss3her = SS_output("gob_her")
mcmc.her = SSreadMCMC("gob_her/mcmc")
The options and output of SSdiagsMCMC() are largely identical to SSdeltaMVLN. In this case, the starter.ss is the same as for ss3phk, only we use Fref = "Btgt" for illustration with $SB_{50}$ and $F_{SB50}$ for this illustration.
mvln = SSdeltaMVLN(ss3her,Fref="Btgt",plot=F,run="mvln") mcmc = SSdiagsMCMC(mcmc.her,ss3her,Fref="Btgt",plot=F,run="mcmc")
Comparing delta-MVLN with MCMC simply requires combining the $kb outputs by, e.g.,
```r$ and $F/F_{50}$ for the ICES Gulf of Bothia Herring SS3 model"} sspar(mfrow=c(1,1),plot.cex = 0.8) SSplotKobe(rbind(mvln$kb,mcmc$kb),joint=F, xlab=mvln$labels[1],ylab=mvln$labels[2],fill=F)
```r$, $F/F_{SB40}$, $SSB$ and $F$ for the ICES Gulf of Bothia Herring SS3 model"}
sspar(mfrow=c(2,2),plot.cex = 0.7)
SSplotEnsemble(rbind(mvln$kb,mcmc$kb),ylabs=mvln$labels,add=T,subplots = c("stock", "harvest", "SSB", "F"),verbose=F)
This works equally for joining a model ensemble.
\newpage
Cookbook Recipies {#s4}
Retrospectives with hindcasts {#r1}
Retrospective analysis can be run for Stock Synthesis using the function r4ss::SS_doRetro() available in r4ss. This setup of the retrospective analysis has the advantage that forecasts are conducted automatically given the catch. This makes it possible to apply retrospective forecasting and hindcast cross-validations of observations based on the same output.
library(r4ss)
Below is a step-by-step cookbook recipe for retrospective analysis in Stock Synthesis.
Step1: Identify restrospective period
Specify the range pf peels that will then determine the end.yr.vec of runs in r4ss::SS_doRetro()
start.retro <- 0 # end year of reference year end.retro <- 7 # number of years for retrospective e.g.,
Step 2: Identify the base directory
Specify the path directory that holds the folder with the base case run. In this case the Pacific Hake folder with a model folder `Reference_Run"
dirname.base = "C:/Users/henni/Dropbox/ss3diags_demo/PacificHake" run = "Reference_Run" model.run <- file.path(dirname.base,run) model.run
Step 3: DAT and CONTROL files
Specify the names of the data and control files.
Note these files are named differently from the DATA.ss and CONTROL.ss. In this case
DAT = "phk.dat" CTL = "phk.ctl"
The names of the DAT and CONTROL are declared on the top of the 'starter.ss`, e.g.
#C Hake starter file
phk.dat
phk.ctl
Step 4: Create a subdirectory for the Retrospectives
There are several ways to organize the retrospective output structure.
First create a new subfolder for the retrospective runs output
dir.retro <- paste0(dirname.base,'/Retro_',run) dir.create(path=dir.retro, showWarnings = F)
Also create a subdirectory for the retrospective model folders
dir.create(path=file.path(dir.retro,"retros"), showWarnings = F)
then copy model run files to the new retrospective folder
file.copy(file.path(model.run,"starter.ss_new"), file.path(dir.retro,"starter.ss")) file.copy(file.path(model.run,"control.ss_new"), file.path(dir.retro, CTL)) file.copy(file.path(model.run,"data.ss_new"), file.path(dir.retro, DAT)) file.copy(file.path(model.run,"forecast.ss"), file.path(dir.retro, "forecast.ss")) file.copy(file.path(model.run,"SS.exe"), file.path(dir.retro, "SS.exe")) # Automatically ignored for models without wtatage.ss file.copy(file.path(model.run,"wtatage.ss"), file.path(dir.retro, "wtatage.ss"))
Step 5: Modify Starter.ss file
Modifying the Starter File helps to speed up model runs
starter <- readLines(paste0(dir.retro, "/starter.ss")) #[8] "2 # run display detail (0,1,2)" linen <- grep("# run display detail", starter) starter[linen] <- paste0( 1 , " # run display detail (0,1,2)" ) # write modified starter.ss write(starter, file.path(dir.retro, "starter.ss"))
Step 6: Execute retrospective runs
Run the retrospective analyses using r4SS function r4ss::SS_doRetro
Ideally, the runs should be done with Hessian to evaluate the retrospective trajectories with respect to the confidence interval coverage of the reference model.
r4ss::SS_doRetro(masterdir=dir.retro, oldsubdir="", newsubdir="", years=start.retro:-end.retro)
However, for larger models it might be desirable to shorten run times, by not inverting the hessian, using the option extras = "-nohess" (much faster)
r4ss::SS_doRetro(masterdir=dir.retro,oldsubdir="", newsubdir="", years=start.retro:-end.retro,extras = "-nohess")
Step 7: Read SS_doRetro() output
retro.phk <- r4ss::SSgetoutput(dirvec=file.path(dir.retro, paste0("retro",start.retro:-end.retro)))
It is often useful to save the retro model runs as .rdata file for further processing with ss3diags, considering that reading the models with r4ss::SSgetoutput() can be quite time-consuming for more complex models.
save(retro.phk, file=file.path(dir.retro,"retro.phk.rdata"))
Step 8: Check
library(ss3diags) check.retro =r4ss::SSsummarize(retro.phk) sspar(mfrow=c(1,1)) SSplotRetro(check.retro,forecast = T,add=T,legendcex = 0.7,legendloc = "topleft")
R0 profiling {#r2}
ASPM diagnostic {#r3}
Jittering {#r4}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.