Getting started

This vignette explains basic and more advanced functions of the spict package. The package is installed from gihtub using the devtools package:

devtools::install_github("DTUAqua/spict/spict")

installs the stable version of spict. When loading the package the user is greeted and the installed version is shown:

library(spict)
   Loading required package: TMB
   Welcome to spict_v1.3.8

The printed version follows the format ver@SHA, where ver is the spict version number and SHA is a unique github commit on github. The content of this vignette pertains to the version printed above that can be found here.

A specific version of spict can be installed using the following:

devtools::install_github("DTUAqua/spict/spict", ref = "v1.3.8")

Loading built-in example data

The package contains the catch and index data analysed in Polacheck, Hilborn, and Punt (1993) that can be loaded using

data(pol)

Data on three stocks are contained in this dataset: South Atlantic albacore, northern Namibian hake, and New Zealand rock lobster. Here focus will be on the South Atlantic albacore data. This dataset contains the following

pol$albacore
   $obsC
    [1] 15.9 25.7 28.5 23.7 25.0 33.3 28.2 19.7 17.5 19.3 21.6 23.1 22.5 22.5 23.6
   [16] 29.1 14.4 13.2 28.4 34.6 37.5 25.9 25.3
   
   $timeC
    [1] 1967 1968 1969 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981
   [16] 1982 1983 1984 1985 1986 1987 1988 1989
   
   $obsI
    [1] 61.89 78.98 55.59 44.61 56.89 38.27 33.84 36.13 41.95 36.63 36.33 38.82
   [13] 34.32 37.64 34.01 32.16 26.88 36.61 30.07 30.75 23.36 22.36 21.91
   
   $timeI
    [1] 1967 1968 1969 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981
   [16] 1982 1983 1984 1985 1986 1987 1988 1989

Note that data are structured as a list containing the entries obsC (catch observations), timeC (time of catch observations), obsI (index observations), and timeI (time of index observations). If times are not specified it is assumed that the first observation is observed at time 1 and then sequentially onward with a time step of one year. It is therefore recommended to always specify observation times.

Each catch observation relates to a time interval. This is specified using dtc. If dtc is left unspecified (as is the case here) each catch observation is assumed to cover the time interval until the next catch observation. For this example with annual catches dtc therefore is

inp <- check.inp(pol$albacore)
inp$dtc
    [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

It is important to specify dtc if the default assumption is not fulfilled.

Note that the index observations for all data sets analysed in Polacheck, Hilborn, and Punt (1993) represent commercial CPUE data and the respective times (inp$timeI) should thus correspond to the middle of the year, e.g. 1988.50, rather than the beginning of the year, e.g. 1988.00. The ‘incorrect’ times of the commercial CPUE observations were kept only for reasons of consistency regarding previous studies using models which were not implemented in ‘continuous time’. We emphasise the importance of correct times of the index observations, which should most accurately correspond to the timing of the survey or to the middle of the year in case of commercial CPUE.

Plotting data

The data can be plotted using the command

plotspict.data(pol$albacore)

Note that the number of catch and index observations are given in the respective plot headers. Furthermore, the color of individual points shows when the observation was made and the corresponding colors are shown in the color legend in the top right corner. For illustrative purposes let’s try shifting the biomass index data a bit

inpshift <- pol$albacore
inpshift$timeC <- inpshift$timeC
inpshift$timeI <- inpshift$timeI + 0.8
plotspict.data(inpshift)

Now the colours show that the index takes place in autumn.

Advanced data plotting

There is also a more advanced function for plotting data, which at the same time does some basic model fitting (linear regression) and shows the results

plotspict.ci(pol$albacore)

The two top plots come from plotspict.data, with the dashed horizontal line representing a guess of MSY. This guess comes from a linear regression between the index and the catch divided by the index (middle row, left). This regression is expected to have a negative slope. A similar plot can be made showing catch versus catch/index (middle row, right) to approximately find the optimal effort (or effort proxy). The proportional increase in the index as a function of catch (bottom row, right) should show primarily positive increases in index at low catches and vice versa. Positive increases in index at large catches could indicate model violations. In the current plot these are not seen.

Fitting the model

The model is fitted to data by running

res <- fit.spict(pol$albacore)

Here the call to fit.spict is wrapped in the system.time command to check the time spent on the calculations. This is obviously not required, but done here to show that fitting the model only takes a few seconds. The result of the model fit is stored in res, which can either be plotted using plot or summarised using summary.

The results are returned as a list that contains output as well as input. The content of this list is

names(res)
    [1] "value"           "sd"              "cov"             "par.fixed"      
    [5] "cov.fixed"       "pdHess"          "gradient.fixed"  "par.random"     
    [9] "diag.cov.random" "env"             "inp"             "obj"            
   [13] "opt"             "pl"              "Cp"              "report"         
   [17] "computing.time"

Many of these variables are generated by TMB::sdreport(). In addition to these spict includes the list of input values (inp), the object used for fitting (obj), the result from the optimiser (opt), the time spent on fitting the model (computing.time), and more less useful variables.

Interpreting summary of results

The results are summarised using

summary(res)

    [1] "Convergence: 0  MSG: relative convergence (4)"                  
    [2] "Objective function at optimum: 2.0654937"                       
    [3] "Euler time step (years):  1/16 or 0.0625"                       
    [4] "Nobs C: 23,  Nobs I1: 23"                                       
    [5] ""                                                               
    [6] "Priors"                                                         
    [7] "     logn  ~  dnorm[log(2), 2^2]"                               
    [8] " logalpha  ~  dnorm[log(1), 2^2]"                               
    [9] "  logbeta  ~  dnorm[log(1), 2^2]"                               
   [10] ""                                                               
   [11] "Model parameter estimates w 95% CI "                            
   [12] "           estimate       cilow       ciupp    log.est  "       
   [13] " alpha    8.5381049   1.2233146  59.5915697  2.1445391  "       
   [14] " beta     0.1212590   0.0180695   0.8137344 -2.1098262  "       
   [15] " r        0.2556015   0.1010611   0.6464616 -1.3641355  "       
   [16] " rc       0.7435358   0.1445758   3.8239150 -0.2963384  "       
   [17] " rold     0.8180032   0.0019102 350.2949907 -0.2008891  "       
   [18] " m       22.5827677  17.0682739  29.8789088  3.1171871  "       
   [19] " K      201.4754010 138.1203391 293.8910914  5.3056673  "       
   [20] " q        0.3512548   0.1942710   0.6350919 -1.0462434  "       
   [21] " n        0.6875299   0.0636729   7.4238412 -0.3746500  "       
   [22] " sdb      0.0128136   0.0018407   0.0891983 -4.3572484  "       
   [23] " sdf      0.3673760   0.2673623   0.5048024 -1.0013694  "       
   [24] " sdi      0.1094038   0.0808977   0.1479547 -2.2127093  "       
   [25] " sdc      0.0445477   0.0073372   0.2704703 -3.1111956  "       
   [26] " "                                                              
   [27] "Deterministic reference points (Drp)"                           
   [28] "         estimate      cilow      ciupp    log.est  "           
   [29] " Bmsyd 60.7442667 15.4035002 239.547239  4.1066727  "           
   [30] " Fmsyd  0.3717679  0.0722879   1.911957 -0.9894856  "           
   [31] " MSYd  22.5827677 17.0682739  29.878909  3.1171871  "           
   [32] "Stochastic reference points (Srp)"                              
   [33] "       estimate     cilow      ciupp    log.est  rel.diff.Drp  "
   [34] " Bmsys 60.73662 15.403659 239.484437  4.1065468 -1.259605e-04  "
   [35] " Fmsys  0.37178  0.072281   1.912265 -0.9894529  3.257870e-05  "
   [36] " MSYs  22.58066 17.062739  29.883028  3.1170939 -9.324958e-05  "
   [37] ""                                                               
   [38] "States w 95% CI (inp$msytype: s)"                               
   [39] "                  estimate      cilow       ciupp    log.est  " 
   [40] " B_1989.94      56.6971159 30.1359483 106.6687174  4.0377233  " 
   [41] " F_1989.94       0.4464499  0.2145076   0.9291863 -0.8064282  " 
   [42] " B_1989.94/Bmsy  0.9334915  0.2961564   2.9423855 -0.0688234  " 
   [43] " F_1989.94/Fmsy  1.2008441  0.2864363   5.0343707  0.1830247  " 
   [44] ""                                                               
   [45] "Predictions w 95% CI (inp$msytype: s)"                          
   [46] "                prediction      cilow      ciupp    log.est  "  
   [47] " B_1991.00      54.3059314 27.8530542 105.881896  3.9946335  "  
   [48] " F_1991.00       0.4464501  0.1573061   1.267069 -0.8064276  "  
   [49] " B_1991.00/Bmsy  0.8941218  0.2498623   3.199578 -0.1119133  "  
   [50] " F_1991.00/Fmsy  1.2008447  0.2390672   6.031894  0.1830252  "  
   [51] " Catch_1990.00  24.7359915 15.3329650  39.905477  3.2082593  "  
   [52] " E(B_inf)       49.9856298         NA         NA  3.9117356  "

Here is a line-by-line description of the summary:

• Line 1: Convergence of the model fit, which has code 0 if the fit was succesful. If this is not the case convergence was not obtained and reported results should not be used. In case of non-convergence results will still be reported to aid diagnosis of the problem.
• Line 2: Objective function value at the optimum. The objective function is the likelihood function if priors are not used and the posterior density function if priors are used.
• Line 3: The Euler time step used in the calculation.
• Line 4: Number of observations for the time series used.
• Line 6-9: Summary of the priors used in the fit. The priors shown here are the default priors that are applied when priors are unspecified. These are relatively uninformative and are applied because most data-limited situations do not allow simultaneous estimation of all noise parameters and logn. The default priors can be disabled (see the section on priors).
• Line 11-25: Summary of the parameter estimates and their 95% CIs. These can be extracted as a data frame with sumspict.parest(res).
• Line 27-31: Estimates of deterministic reference points with 95% CIs. These are the reference points one would derive if stochasticity were ignored. Can be extracted with sumspict.drefpoints(res).
• Line 32-36: Estimates of stochastic reference points with 95% CIs. These are the reference points of the stochastic model. The column ´rel.diff.Drp´ shows the relative difference when compared to the deterministic reference points. The information can be extracted with sumspict.srefpoints(res).
• Line 38-43: State estimates in the final year where data were available. The states of the model are biomass (B) and fishing mortality (F) with the year of the estimates appended. The year is shown as a decimal number as estimates within year are possible. Both absolute (B and F) and relative estimates (B/Bmsy and F/Fmsy) are shown. The relative estimates are calculated using the type of reference points given by msytype (line 38), where s is stochastic and d is deterministic. Here msytype is ´s´. This information can be extracted using sumspict.states(res).
• Line 45-52: Predictions of absolute and relative biomass and fishing mortality at the time indicated by inp$timepredi, here 1990 (line 47-50). In addition, predicted catch at the time indicated by inp$timepredc (line 51). Finally, the equilibrium biomass, indicated by E(B_inf), if current conditions remain constant. There predictions or forecasts are calculated under the fishing scenario given by inp$ffac. See the section on forecasting for more information. The prediction summary can be extracted using sumspict.predictions(res).

Interpreting plots of results

spict comes with several plotting abilities. The basic plotting of the results is done using the generic function plot that produces a multipanel plot with the most important outputs.

plot(res)

Some general comments can be made regarding the style and colours of these plots:

• Estimates (biomass, fishing mortality, catch, production) are shown using blue lines.
• 95% CIs of absolute quantities are shown using dashed blue lines.
• 95% CIs of relative biomass and fishing mortality are shown using shaded blue regions.
• Estimates of reference points (B_MSY, F_MSY, MSY) are shown using black lines.
• 95% CIs of reference points are shown using grey shaded regions.
• The end of the data range is shown using a vertical grey line.
• Predictions beyond the data range are shown using dotted blue lines.
• Data are shown using points colored by season. Different index series use different point characters (not shown here).

The individual plots can be plotted separately using the plotspict.* family of plotting functions; all functions are summarised in Table 1 and their common arguments that control their look in Table 2:

We will now look at them one at a time. The top left is the plot of absolute biomass

plotspict.biomass(res)

Note that this plot has a y-axis on the right side related to the relative biomass (B_t/B_MSY). The shaded 95% CI region relates to this axis, while the dashed blue lines relate to the left y-axis indicating absolute levels. The dashed lines and the shaded region are shown on the same plot to make it easier to assess whether the relative or absolute levels are most accurately estimated. Here, the absolute are more accurate than the relative. Later, we will see examples of the opposite. The horizontal black line is the estimate of B_MSY with 95% CI shown as a grey region.

The plot of the relative biomass is produced using

plotspict.bbmsy(res)

This plot contains much of the same information as given by plotspict.biomass, but without the information about absolute biomass and without the 95% CI around the B_MSY reference point.

The plots of fishing mortality follow the same principles

plotspict.f(res, main='', qlegend=FALSE, rel.axes=FALSE, rel.ci=FALSE)
plotspict.ffmsy(res, main='', qlegend=FALSE)

The estimate of F_MSY is shown with a horizontal black line with 95% CI shown as a grey region (left plot). The 95% CI of F_MSY is very wide in this case. As shown here it is quite straightforward to remove the information about relative levels from the plot of absolute fishing mortality. Furthermore, the argument main='' removes the heading and qlegend=FALSE removes the colour legend for data points.

The plot of the catch is produced using

plotspict.catch(res)

This plot shows estimated catches (blue line) versus observed catches (points) with the estimate of MSY plotted as a horizontal black line with its 95% CI given by the grey region.

A phase plot (or kobe plot) of fishing mortality versus biomass is plotted using

plotspict.fb(res, ylim=c(0, 1.3), xlim=c(0, 300))

The plot shows the development of biomass and fishing mortality since the initial year (here 1967) indicated with a circle until the terminal year (here 1990) indicated with a square. The yellow diamond indicates the mean biomass over a long period if the current (1990) fishing pressure remains. This point can be interpreted as the fished equilibrium and is denoted E(B_∞) in the legend as a statistical way of expressing the expectation of the biomass as t → ∞. As the current fishing mortality is close to F_MSY the expected long term biomass is close to B_MSY.

A vertical dashed red line at B_t = 0 indicates the biomass level below which the stock has crashed. The grey shaded banana-shaped area indicates the 95% confidence region of the pair F_MSY, B_MSY. This region is important to visualise jointly as the two reference points are highly (negatively) correlated.

Residuals and diagnostics

Before proceeding with the results for an actual assessment it is very important that the model residuals are checked and possible model deficiencies identified. Residuals can be calculated using calc.osa.resid(). OSA stands for one-step-ahead, which are the proper residuals for state-space models. More information about OSA residuals is contained in Pedersen and Berg (2017). To calculate and plot residuals and diagnostics do

res <- calc.osa.resid(res)
plotspict.diagnostic(res)

The first column of the plot contains information related to catch data and the second column contains information related to the index data. The rows contain

1. Log of the input data series.
2. OSA residuals with the p-value of a test for bias (i.e. that the mean of the residuals is different from zero) in the plot header. If the header is green the test was not significant, otherwise the header would be red.
3. Empirical autocorrelation of the residuals. Two tests for significant autocorrelation is performed. A Ljung-Box simultaneous test of multiple lags (here 4) with p-value shown in the header, and tests for individual lags shown by dashed horizontal lines in the plot. Here no violation is identified.
4. Tests for normality of the residuals both as a QQ-plot and with a Shapiro test with p-value shown in the plot header.

This data did not have any significant violations of the assumptions, which increases confidence in the results. For a discussion of possible violations and remedies the reader is referred to Pedersen and Berg (2017).

res <- calc.process.resid(res)
plotspict.diagnostic.process(res)

Extracting parameter estimates

To extract an estimated quantity, here logBmsy use

get.par('logBmsy', res)
                 ll      est       ul        sd        cv
   logBmsy 2.734605 4.106547 5.478488 0.6999831 0.1704554

This returns a vector with ll being the lower 95% limit of the CI, est being the estimated value, ul being the upper 95% limit of the CI, sd being the standard deviation of the estimate, and cv being the coefficient of variation of the estimate. The estimated quantity can also be returned on the natural scale (as opposed to log scale) by running

get.par('logBmsy', res, exp=TRUE)
                 ll      est       ul        sd        cv
   logBmsy 15.40366 60.73662 239.4844 0.6999831 0.7951589

This essentially takes the exponential of ll, est and ul of the values in log, while sd is unchanged as it is the standard deviation of the quantity on the scale that it is estimated (here log). When transforming using exp=TRUE the $CV = \sqrt{e^{\sigma^2}-1}$. Most parameters are log-transformed under estimation and should therefore be extracted using exp=TRUE.

For a standard fit (not using robust observation error, seasonality etc.), the quantities that can be extracted using this method are

list.quantities(res)
    [1] "Bm"                 "Bmsy"               "Bmsy2"             
    [4] "BmsyB0"             "Bmsyd"              "Bmsys"             
    [7] "Cp"                 "Emsy"               "Emsy2"             
   [10] "Fmsy"               "Fmsyd"              "Fmsys"             
   [13] "K"                  "MSY"                "MSYd"              
   [16] "MSYs"               "gamma"              "isdb2"             
   [19] "isdc2"              "isde2"              "isdf2"             
   [22] "isdi2"              "logB"               "logBBmsy"          
   [25] "logBl"              "logBlBmsy"          "logBlK"            
   [28] "logBm"              "logBmBmsy"          "logBmK"            
   [31] "logBmsy"            "logBmsyPluslogFmsy" "logBmsyd"          
   [34] "logBmsys"           "logBp"              "logBpBmsy"         
   [37] "logBpK"             "logCp"              "logCpred"          
   [40] "logEmsy"            "logEmsy2"           "logEp"             
   [43] "logF"               "logFFmsy"           "logFFmsynotS"      
   [46] "logFl"              "logFlFmsy"          "logFlFmsynotS"     
   [49] "logFm"              "logFmFmsy"          "logFmFmsynotS"     
   [52] "logFmsy"            "logFmsyd"           "logFmsys"          
   [55] "logFnotS"           "logFp"              "logFpFmsy"         
   [58] "logFpFmsynotS"      "logFs"              "logIp"             
   [61] "logIpred"           "logK"               "logMSY"            
   [64] "logMSYd"            "logMSYs"            "logalpha"          
   [67] "logbeta"            "logbkfrac"          "logm"              
   [70] "logn"               "logq"               "logq2"             
   [73] "logr"               "logrc"              "logrold"           
   [76] "logsdb"             "logsdc"             "logsdf"            
   [79] "logsdi"             "m"                  "p"                 
   [82] "q"                  "r"                  "rc"                
   [85] "rold"               "sdb"                "sdc"               
   [88] "sde"                "sdf"                "sdi"               
   [91] "seasonsplinefine"

These should be relatively self-explanatory when knowing that reference points ending with s are stochastic and those ending with d are deterministic, quantities ending with p are predictions and quantities ending with l are estimates in the final year. If a quantity is available both on natural and log scale, it is preferred to transform the quantity from log as most quantities are estimated on the log scale.

res$cov.fixed
                   logm         logK         logq         logn       logsdb
   logm    0.0204039157 -0.005706842  0.026834631 -0.156739841  0.032239968
   logK   -0.0057068418  0.037105158 -0.038075151 -0.011341936 -0.021784067
   logq    0.0268346313 -0.038075151  0.091311499 -0.227633882  0.040958484
   logn   -0.1567398408 -0.011341936 -0.227633882  1.473734413 -0.190011439
   logsdb  0.0322399684 -0.021784067  0.040958484 -0.190011439  0.980090469
   logsdf  0.0014253331 -0.002407442  0.003270286 -0.006648007 -0.002144168
   logsdi -0.0007603661 -0.002521063  0.002848536  0.007158184  0.010535657
   logsdc  0.0015534798  0.001300449 -0.008392720  0.001998025  0.137919960
                 logsdf        logsdi       logsdc
   logm    0.0014253331 -0.0007603661  0.001553480
   logK   -0.0024074422 -0.0025210626  0.001300449
   logq    0.0032702864  0.0028485361 -0.008392720
   logn   -0.0066480073  0.0071581839  0.001998025
   logsdb -0.0021441680  0.0105356567  0.137919960
   logsdf  0.0262881625 -0.0002784344 -0.035159250
   logsdi -0.0002784344  0.0237200081  0.005885894
   logsdc -0.0351592502  0.0058858943  0.846809016

cov2cor(res$cov.fixed)
                 logm         logK        logq         logn      logsdb
   logm    1.00000000 -0.207406292  0.62169323 -0.903884687  0.22798421
   logK   -0.20740629  1.000000000 -0.65412652 -0.048502088 -0.11423225
   logq    0.62169323 -0.654126516  1.00000000 -0.620532535  0.13691406
   logn   -0.90388469 -0.048502088 -0.62053254  1.000000000 -0.15810189
   logsdb  0.22798421 -0.114232253  0.13691406 -0.158101887  1.00000000
   logsdf  0.06154312 -0.077083002  0.06674872 -0.033775498 -0.01335813
   logsdi -0.03456277 -0.084978502  0.06120707  0.038285632  0.06909890
   logsdc  0.01181833  0.007336409 -0.03018195  0.001788539  0.15139143
               logsdf      logsdi       logsdc
   logm    0.06154312 -0.03456277  0.011818330
   logK   -0.07708300 -0.08497850  0.007336409
   logq    0.06674872  0.06120707 -0.030181947
   logn   -0.03377550  0.03828563  0.001788539
   logsdb -0.01335813  0.06909890  0.151391434
   logsdf  1.00000000 -0.01115027 -0.235649625
   logsdi -0.01115027  1.00000000  0.041530036
   logsdc -0.23564963  0.04153004  1.000000000

cov2cor(get.cov(res, 'logBmsy', 'logFmsy'))
              [,1]       [,2]
   [1,]  1.0000000 -0.9982507
   [2,] -0.9982507  1.0000000

Comparison plots

The function plotspict.compare() can be used to make comparison plots of one or more spict runs. The function is able to show biomass and fishing mortality (absolute and relative), catch, and production curve. You can use the varname argument to select which of those plots to produce, default is all of them (c("B", "F", "C", "BBmsy", "FFmsy", "P")).

As an example we show below the comparison of two runs with and without logbkfrac prior.

data(pol)
inp <- pol$hake
rep1 <- fit.spict(inp)
inp$priors$logbkfrac <- c(log(0.3), 0.1, 1)
rep2 <- fit.spict(inp)
plotspict.compare(list(def = rep1, bkprior = rep2))

Retrospective plots

Retrospective plots are sometimes used to evaluate the robustness of the model fit to the introduction of new data, i.e. to check whether the fit changes substantially when new data becomes available. Such calculations and plotting thereof can be performed using retro() as shown here

## res <- fit.spict(pol$albacore)
res <- retro(res, nretroyear = 4, mc.cores = 1)
plotspict.retro(res)
        FFmsy      BBmsy 
    0.4668503 -0.2441031

By default retro creates 5 scenarios with catch and index time series which are shortened by the 1 to 5 last observations. The number of scenarios and thus observations which are removed can be changed with the argument nretroyear in the function retro. The graphs show the different peels with different colors. For the relative biomass and fishing mortality, the Mohn’s rho are shown inside the corresponding plots. Mohn’s rho for other quantities can be calculated using the function mohns_rho.

## res <- fit.spict(pol$albacore)
## res <- retro(res)
mohns_rho(res, what = c("FFmsy", "BBmsy"))
        FFmsy      BBmsy 
    0.4668503 -0.2441031

The function plotspict.retro.fixed makes a plot of the paramter estimates, with corresponding 95% confidence intervals for each of the retrospective runs.

## res <- fit.spict(pol$albacore)
## res <- retro(res)
plotspict.retro.fixed(res)

Estimation using two or more biomass indices

The estimation can be done using more than one biomass index, for example when scientific surveys are performed more than once every year or when there are both commercial and survey CPUE time-series available. The following example emulates a situation where a long but noisy first quarter index series and a shorter and less noisy second quarter index series are available with different catchabilities

set.seed(123)
inp <- list(timeC=pol$albacore$timeC, obsC=pol$albacore$obsC)
inp$timeI <- list(pol$albacore$timeI, pol$albacore$timeI[10:23]+0.25)
inp$obsI <- list()
inp$obsI[[1]] <- pol$albacore$obsI * exp(rnorm(23, sd=0.1)) # Index 1
inp$obsI[[2]] <- 10*pol$albacore$obsI[10:23] # Index 2
res <- fit.spict(inp)
sumspict.parest(res)
              estimate        cilow       ciupp    log.est
   alpha1   5.54953281   0.94907993  32.4496530  1.7137137
   alpha2   3.45424303   0.58225750  20.4922992  1.2396033
   beta     0.12057032   0.01913001   0.7599159 -2.1155222
   r        0.18601585   0.08916105   0.3880831 -1.6819234
   rc       1.20960483   0.19388362   7.5465058  0.1902937
   rold     0.26864004   0.05297139   1.3623857 -1.3143829
   m       24.30175989  17.98207990  32.8424485  3.1905488
   K      220.56430343 148.87830749 326.7676317  5.3961893
   q1       0.35403668   0.22557324   0.5556598 -1.0383548
   q2       3.60405904   2.32439321   5.5882290  1.2820607
   n        0.30756467   0.03135919   3.0165330 -1.1790699
   sdb      0.02006319   0.00372284   0.1081248 -3.9088686
   sdf      0.36879350   0.26948920   0.5046905 -0.9975184
   sdi1     0.11134131   0.08158825   0.1519445 -2.1951549
   sdi2     0.06930312   0.04587601   0.1046936 -2.6692653
   sdc      0.04446555   0.00764514   0.2586198 -3.1130406
plotspict.biomass(res)

The model estimates seperate observation noises and finds that the first index (sdi1) is more noisy than the second (sdi2). It is furthermore estimated that the catchabilities are different by a factor 10 (q1 versus q2). The biomass plot shows both indices, with circles indicating the first index and squares indicating the second index (the two series can also be distinguished by their colours). It is possible to force the model to assume that two or more index time-series have the same level of observation noise (CV). For example, to assume that sdi1 equals sdi2 one must add inp$mapsdi <- c(1,1) before calling fit.spict(inp). The length of mapsdi should equal the number of indices. In case of 3 index series one could for example use inp$mapsdi <- c(1, 1, 2) to have series 1 and 2 share sdi and have a separate sdi for series 3.

Using effort data instead of commercial CPUE

It is possible to use effort data directly in the model instead of calculating commercial CPUE and inputting this as an index. It is beyond the scope of this vignette to discuss all problems associated with indices based on commercial CPUEs, however it is intuitively clear that using the same information twice (catch as catch and catch in catch/effort) induces a correlation, which the model does not account for. These problems are easily avoided by putting catch and effort seperately

inpeff <- list(timeC=pol$albacore$timeC, obsC=pol$albacore$obsC,
               timeE=pol$albacore$timeC, obsE=pol$albacore$obsC/pol$albacore$obsI)
repeff <- fit.spict(inpeff)
sumspict.parest(repeff)
              estimate        cilow        ciupp    log.est
   beta     0.07385348   0.01464766   0.37236911 -2.6056722
   r        0.23822940   0.10657013   0.53254366 -1.4345212
   rc       1.14399151   0.21452927   6.10041032  0.1345235
   rold     0.40826824   0.03851526   4.32771228 -0.8958309
   m       24.16376089  18.62126768  31.35593936  3.1848540
   K      189.53177582 132.49614433 271.11954260  5.2445567
   qf       0.41117178   0.25562848   0.66135915 -0.8887442
   n        0.41648805   0.04180369   4.14944936 -0.8758975
   sdb      0.01430671   0.00236760   0.08645135 -4.2470266
   sdf      0.37625010   0.27938937   0.50669121 -0.9775012
   sde      0.09820098   0.06985386   0.13805153 -2.3207391
   sdc      0.02778738   0.00568278   0.13587327 -3.5831734
par(mfrow=c(2, 2))
plotspict.bbmsy(repeff)
plotspict.ffmsy(repeff, qlegend=FALSE)
plotspict.catch(repeff, qlegend=FALSE)
plotspict.fb(repeff)

Here the model runs without an index of biomass and instead uses effort as an index of fishing mortality Note that index observations are missing from the biomass plot, but effort observations are present in the plot of fishing mortality. Note also that q is missing from the summary of parameter estimates and instead qf is present, which is the commercial catchability.

Overall for this data set the results in terms of stock status etc. do not change much, and this will probably often be the case, however using effort data directly instead of commercial CPUE is cleaner and avoids inputting the same data twice.

Scaling the uncertainty of individual data points

It is not always appropriate to assume that the observation noise of a data series is constant in time. Knowledge that certain data points are more uncertain than others can be implemented using stdevfacC, stdevfacI, and stdevfacE, which are vectors containing factors that are multiplied onto the standard deviation of the data points of the corresponding observation vectors. An example where the first 10 years of the biomass index are considered uncertain relative to the remaining time series and therefore are scaled by a factor 5.

inp <- pol$albacore
res1 <- fit.spict(inp)
inp$stdevfacC <- rep(1, length(inp$obsC))
inp$stdevfacC[1:10] <- 5
res2 <- fit.spict(inp)
par(mfrow=c(2, 1))
plotspict.catch(res1, main='No scaling')
plotspict.catch(res2, main='With scaling', qlegend=FALSE)

From the plot it is noted that the scaling factor widens the 95% CIs of the initial ten years of catch data, while narrowing the 95% CIs of the remaining years.

Simulating data

The package has built-in functionality for simulating data, which is useful for testing.

inp <- check.inp(pol$albacore)
sim <- sim.spict(inp)
   The stabilise option was used for fitting / is specified in the input list (inp$stabilise). This activates six very vague priors. Acknowledging these priors when simulating is not yet implemented. 
   Additional priors were used for fitting / are specified in the input list (inp$priors): logn, logalpha, logbeta. Acknowledging these priors when simulating is not yet implemented.
plotspict.data(sim)

inp <- list(ini=list(logK=log(100), logm=log(10), logq=log(1)))
sim <- sim.spict(inp, nobs=50)
   The stabilise option was used for fitting / is specified in the input list (inp$stabilise). This activates six very vague priors. Acknowledging these priors when simulating is not yet implemented. 
   Additional priors were used for fitting / are specified in the input list (inp$priors): logn, logalpha, logbeta. Acknowledging these priors when simulating is not yet implemented.
plotspict.data(sim)

set.seed(31415926)
inp <- list(ini=list(logK=log(100), logm=log(10), logq=log(1),
                     logbkfrac=log(1), logF0=log(0.3), logsdc=log(0.1),
                     logsdf=log(0.3)))
sim <- sim.spict(inp, nobs=30)
   The stabilise option was used for fitting / is specified in the input list (inp$stabilise). This activates six very vague priors. Acknowledging these priors when simulating is not yet implemented. 
   Additional priors were used for fitting / are specified in the input list (inp$priors): logn, logalpha, logbeta. Acknowledging these priors when simulating is not yet implemented.
res <- fit.spict(sim)
sumspict.parest(res)
             estimate  true       cilow       ciupp true.in.ci     log.est
   alpha   1.04607706  -9.0  0.31605055   3.4623487         -9  0.04504703
   beta    0.13191759  -9.0  0.02269952   0.7666352         -9 -2.02557789
   r       1.07672249  -9.0  0.35062210   3.3064982         -9  0.07392170
   rc      0.63474119  -9.0  0.34468649   1.1688778         -9 -0.45453794
   rold    0.45001541  -9.0  0.22253267   0.9100411         -9 -0.79847345
   m      14.48076842  10.0 10.46879070  20.0302652          0  2.67282145
   K      76.02606873 100.0 46.11361200 125.3417999          1  4.33107629
   q       1.19617690   1.0  0.85014237   1.6830583          1  0.17913055
   n       3.39263471   2.0  1.36227805   8.4490610          1  1.22160682
   sdb     0.19525936   0.2  0.08858085   0.4304115          1 -1.63342653
   sdf     0.34363005   0.3  0.25198877   0.4685987          1 -1.06818963
   sdi     0.20425634   0.2  0.12167035   0.3428991          1 -1.58837950
   sdc     0.04533085   0.1  0.00814495   0.2522895          1 -3.09376752
par(mfrow=c(2, 2))
plotspict.biomass(res)
plotspict.f(res, qlegend=FALSE)
plotspict.catch(res, qlegend=FALSE)
plotspict.fb(res)

set.seed(1234)
inp <- list(nseasons=4, splineorder=3)
inp$timeC <- seq(0, 30-1/inp$nseasons, by=1/inp$nseasons)
inp$timeI <- seq(0, 30-1/inp$nseasons, by=1/inp$nseasons)
inp$ini <- list(logK=log(100), logm=log(20), logq=log(1),
                logbkfrac=log(1), logsdf=log(0.4), logF0=log(0.5),
                logphi=log(c(0.05, 0.1, 1.8)))
seasonsim <- sim.spict(inp)
   The stabilise option was used for fitting / is specified in the input list (inp$stabilise). This activates six very vague priors. Acknowledging these priors when simulating is not yet implemented. 
   Additional priors were used for fitting / are specified in the input list (inp$priors): logn, logalpha, logbeta. Acknowledging these priors when simulating is not yet implemented.
plotspict.data(seasonsim)

set.seed(432)
inp <- list(nseasons=4, seasontype=2)
inp$timeC <- seq(0, 30-1/inp$nseasons, by=1/inp$nseasons)
inp$timeI <- seq(0, 30-1/inp$nseasons, by=1/inp$nseasons)
inp$ini <- list(logK=log(100), logm=log(20), logq=log(1),
                logbkfrac=log(1), logsdf=log(0.4), logF0=log(0.5))
seasonsim2 <- sim.spict(inp)
   The stabilise option was used for fitting / is specified in the input list (inp$stabilise). This activates six very vague priors. Acknowledging these priors when simulating is not yet implemented. 
   Additional priors were used for fitting / are specified in the input list (inp$priors): logn, logalpha, logbeta. Acknowledging these priors when simulating is not yet implemented.
plotspict.data(seasonsim2)

Estimation using quarterly data

Catch information available in sub-annual aggregations, e.g. quarterly catch, can be used to estimate the seasonal pattern of the fishing mortality. The user can choose between two types of seasonality by setting seasontype to 1, 2 or 3:

1. using cyclic B-splines.
2. using coupled stochastic differential equations (SDEs).
3. combination of a cyclic spline (fixed seasonal pattern) and mean-zero autoregressive process

Technical description of the first two season types is found in Pedersen and Berg (2017). The third is described by the following equation:

Here, an example of a spline-based model fitted to quarterly data simulated in section is shown

seasonres <- fit.spict(seasonsim)
plotspict.biomass(seasonres)
plotspict.f(seasonres, qlegend=FALSE)
plotspict.season(seasonres)

The model is able to estimate the seasonal variation in fishing mortality as seen both in the plot of F and in the plot of the estimated spline, where blue is the estimated spline, orange is the true spline, and green is the spline if time were truly continuous (it is discretised with the Euler steps shown by the blue line).

To fit the coupled SDE model run

seasonres2 <- fit.spict(seasonsim2)
sumspict.parest(seasonres2)
             estimate  true       cilow       ciupp true.in.ci    log.est
   alpha    1.1638022  -9.0  0.67490896   2.0068419         -9  0.1516924
   beta     0.4961820  -9.0  0.30347470   0.8112589         -9 -0.7008125
   r        0.8526430  -9.0  0.26693957   2.7234630         -9 -0.1594144
   rc       0.8468346  -9.0  0.56021633   1.2800926         -9 -0.1662499
   rold     0.8411047  -9.0  0.42324337   1.6715140         -9 -0.1730391
   m       22.6501793  20.0 18.08395681  28.3693789          1  3.1201678
   K      106.7057696 100.0 61.53788893 185.0261924          1  4.6700752
   q        1.1111248   1.0  0.81097416   1.5223647          1  0.1053729
   n        2.0137179   2.0  0.84869990   4.7779668          1  0.6999827
   sdb      0.1619837   0.2  0.10231809   0.2564425          1 -1.8202597
   sdu      0.1259182   0.1  0.07102503   0.2232368          1 -2.0721226
   sdf      0.3565349   0.4  0.25806035   0.4925868          1 -1.0313233
   sdi      0.1885170   0.2  0.15969214   0.2225448          1 -1.6685673
   sdc      0.1769062   0.2  0.13548624   0.2309887          1 -1.7321358
   lambda   0.1032635   0.1  0.02267910   0.4701841          1 -2.2704712
plotspict.biomass(seasonres2)
plotspict.f(seasonres2, qlegend=FALSE)

Two parameters related to the coupled SDEs are estimated (sdu and lambda) as evident from the summary of estimated parameters. In the plot of fishing mortality it is noted that the amplitude of the seasonal pattern varies over time. This is a property of the coupled SDE model, which is not possible to obtain with the spline based seasonal model. The spline based model has a fixed amplitude and phases, which will lead to biased estimates and autocorrelation in residuals if in reality the seasonal pattern shifts a bit. This is illustrated by fitting a spline based model to data generated with a coupled SDE model

inp2 <- list(obsC=seasonsim2$obsC, obsI=seasonsim2$obsI,
             timeC=seasonsim2$timeC, timeI=seasonsim2$timeI,
             seasontype=1, true=seasonsim2$true)
rep2 <- fit.spict(inp2)
rep2 <- calc.osa.resid(rep2)
plotspict.diagnostic(rep2)

From the diagnostics it is clear that autocorrelation is present in the catch residuals.

Setting initial parameter values

Initial parameter values used as starting guess of the optimiser can be set using inp$ini. For example, to specify the initial value of logK set

inp <- pol$albacore
inp$ini$logK <- log(100)

This procedure generalises to all other model parameters. If initial values are not specified they are set to default values. To see the default initial value of a parameter, here logK, run

inp <- check.inp(pol$albacore)
inp$ini$logK
   [1] 5.010635

This can also be done posterior to fitting the model by printing res$inp$ini$logK.

set.seed(123)
check.ini(pol$albacore, ntrials=4)
   Checking sensitivity of fit to initial parameter values...
    Trial 1 ... model fitted!
    Trial 2 ... model fitted!
    Trial 3 ... model fitted!
    Trial 4 ... model fitted!
   $propchng
            logm  logK  logq  logn logsdb logsdf logsdi logsdc
   Trial 1 -1.41  0.26 -0.12 -2.75  -1.26   1.30  -0.08  -1.12
   Trial 2  0.34 -0.04  0.62  0.34  -0.51  -0.21   1.14  -1.14
   Trial 3 -1.69 -0.42 -0.23 -3.26  -1.11  -0.55  -0.40  -1.41
   Trial 4  1.03  0.19  0.06 -0.68   0.60   1.01  -1.32  -1.15
   
   $inimat
           Distance  logn logK logm  logq logsdb logsdf logsdi logsdc
   Basevec     0.00  0.69 5.01 3.41 -0.64  -1.61  -1.61  -1.61  -1.61
   Trial 1     4.22 -0.29 6.34 2.99  1.12   0.42  -3.70  -1.48   0.20
   Trial 2     3.48  0.93 4.81 5.51 -0.86  -0.79  -1.27  -3.44   0.23
   Trial 3     4.52 -0.48 2.90 2.61  1.45   0.18  -0.72  -0.96   0.67
   Trial 4     3.64  1.41 5.97 3.61 -0.21  -2.58  -3.23   0.52   0.24
   
   $resmat
           Distance     m      K    q    n  sdb  sdf  sdi  sdc
   Basevec        0 22.58 201.48 0.35 0.69 0.01 0.37 0.11 0.04
   Trial 1        0 22.58 201.48 0.35 0.69 0.01 0.37 0.11 0.04
   Trial 2        0 22.58 201.48 0.35 0.69 0.01 0.37 0.11 0.04
   Trial 3        0 22.58 201.48 0.35 0.69 0.01 0.37 0.11 0.04
   Trial 4        0 22.58 201.48 0.35 0.69 0.01 0.37 0.11 0.04
   $obsC
    [1] 15.9 25.7 28.5 23.7 25.0 33.3 28.2 19.7 17.5 19.3 21.6 23.1 22.5 22.5 23.6
   [16] 29.1 14.4 13.2 28.4 34.6 37.5 25.9 25.3
   
   $timeC
    [1] 1967 1968 1969 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981
   [16] 1982 1983 1984 1985 1986 1987 1988 1989
   
   $obsI
    [1] 61.89 78.98 55.59 44.61 56.89 38.27 33.84 36.13 41.95 36.63 36.33 38.82
   [13] 34.32 37.64 34.01 32.16 26.88 36.61 30.07 30.75 23.36 22.36 21.91
   
   $timeI
    [1] 1967 1968 1969 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981
   [16] 1982 1983 1984 1985 1986 1987 1988 1989
   
   $check.ini
   $check.ini$propchng
            logm  logK  logq  logn logsdb logsdf logsdi logsdc
   Trial 1 -1.41  0.26 -0.12 -2.75  -1.26   1.30  -0.08  -1.12
   Trial 2  0.34 -0.04  0.62  0.34  -0.51  -0.21   1.14  -1.14
   Trial 3 -1.69 -0.42 -0.23 -3.26  -1.11  -0.55  -0.40  -1.41
   Trial 4  1.03  0.19  0.06 -0.68   0.60   1.01  -1.32  -1.15
   
   $check.ini$inimat
           Distance  logn logK logm  logq logsdb logsdf logsdi logsdc
   Basevec     0.00  0.69 5.01 3.41 -0.64  -1.61  -1.61  -1.61  -1.61
   Trial 1     4.22 -0.29 6.34 2.99  1.12   0.42  -3.70  -1.48   0.20
   Trial 2     3.48  0.93 4.81 5.51 -0.86  -0.79  -1.27  -3.44   0.23
   Trial 3     4.52 -0.48 2.90 2.61  1.45   0.18  -0.72  -0.96   0.67
   Trial 4     3.64  1.41 5.97 3.61 -0.21  -2.58  -3.23   0.52   0.24
   
   $check.ini$resmat
           Distance     m      K    q    n  sdb  sdf  sdi  sdc
   Basevec        0 22.58 201.48 0.35 0.69 0.01 0.37 0.11 0.04
   Trial 1        0 22.58 201.48 0.35 0.69 0.01 0.37 0.11 0.04
   Trial 2        0 22.58 201.48 0.35 0.69 0.01 0.37 0.11 0.04
   Trial 3        0 22.58 201.48 0.35 0.69 0.01 0.37 0.11 0.04
   Trial 4        0 22.58 201.48 0.35 0.69 0.01 0.37 0.11 0.04

Phases and how to fix parameters

The package has the ability to estimate parameters in phases. Users familiar with AD model builder will know that this means that some parameters are held constant in phase 1, some are then released and estimated in phase 2, more are released in phase 3 etc. until all parameters are estimated. Per default all parameters are estimated in phase 1. As an example the standard deviation on the biomass process, logsdb, is estimated in phase 2:

inp <- pol$albacore
inp$phases$logsdb <- 2
res <- fit.spict(inp)
   Estimating - phase 1 
   Estimating - phase 2

Phases can also be used to fix parameters to their initial value by setting the phase to -1. For example

inp <- pol$albacore
inp$phases$logsdb <- -1
inp$ini$logsdb <- log(0.1)
res <- fit.spict(inp)
summary(res)
   Convergence: 0  MSG: relative convergence (4)
   Objective function at optimum: 5.8647408
   Euler time step (years):  1/16 or 0.0625
   Nobs C: 23,  Nobs I1: 23
   
   Priors
        logn  ~  dnorm[log(2), 2^2]
    logalpha  ~  dnorm[log(1), 2^2]
     logbeta  ~  dnorm[log(1), 2^2]
   
   Fixed parameters
        fixed.value  
    sdb         0.1  
   
   Model parameter estimates w 95% CI 
              estimate      cilow       ciupp    log.est  
    alpha    1.0503613  0.6791573   1.6244527  0.0491342  
    beta     0.1713918  0.0256782   1.1439720 -1.7638031  
    r        0.3471988  0.1465267   0.8226963 -1.0578579  
    rc       1.7791120  0.2676182  11.8274426  0.5761143  
    rold     0.5694637  0.0689506   4.7032055 -0.5630603  
    m       27.4846393 18.9164856  39.9337073  3.3136273  
    K      144.5708271 82.7920750 252.4483661  4.9737695  
    q        0.4966500  0.2541128   0.9706762 -0.6998696  
    n        0.3903057  0.0392963   3.8766587 -0.9408250  
    sdf      0.3738950  0.2594031   0.5389198 -0.9837802  
    sdi      0.1050361  0.0679157   0.1624453 -2.2534509  
    sdc      0.0640825  0.0113806   0.3608403 -2.7475832  
    
   Deterministic reference points (Drp)
           estimate      cilow      ciupp    log.est  
    Bmsyd 30.897032  6.2840156 151.913464  3.4306601  
    Fmsyd  0.889556  0.1338091   5.913721 -0.1170328  
    MSYd  27.484639 18.9164856  39.933707  3.3136273  
   Stochastic reference points (Srp)
            estimate      cilow      ciupp    log.est  rel.diff.Drp  
    Bmsys 30.8170236  6.3482980 149.597411  3.4280673 -0.0025962352  
    Fmsys  0.8901021  0.1349423   5.871261 -0.1164191  0.0006135517  
    MSYs  27.4303399 18.8038684  40.014296  3.3116497 -0.0019795390  
   
   States w 95% CI (inp$msytype: s)
                     estimate      cilow     ciupp    log.est  
    B_1989.94      44.7251083 21.9038498 91.323458  3.8005351  
    F_1989.94       0.5737431  0.2575356  1.278197 -0.5555735  
    B_1989.94/Bmsy  1.4513117  0.3383892  6.224506  0.3724678  
    F_1989.94/Fmsy  0.6445813  0.1047014  3.968286 -0.4391543  
   
   Predictions w 95% CI (inp$msytype: s)
                   prediction      cilow     ciupp    log.est  
    B_1991.00      45.2794307 20.9341261 97.937064  3.8128529  
    F_1991.00       0.5737434  0.1907893  1.725367 -0.5555730  
    B_1991.00/Bmsy  1.4692993  0.3184910  6.778340  0.3847856  
    F_1991.00/Fmsy  0.6445816  0.0900523  4.613825 -0.4391539  
    Catch_1990.00  25.8407665 15.9510708 41.862093  3.2519533  
    E(B_inf)       45.9908947         NA        NA  3.8284434

Priors

SPiCT is a generalisation of previous surplus production models in the sense that stochastic noise is included in both observation and state processes of both fishing and biomass. Estimating all model parameters is only possible if data contain sufficient information, which may not be the case for short time series or time series with limited contrast. The basic data requirements of the model are limited to only catch and biomass index time series. More information may be available, which can be used to improve the model fit. This is particularly advantageous if the model is not able to converge with only catch and index time series. Additional information can then be included in the fit via prior distributions for model parameters.

inp <- pol$albacore
inp$priors$logn <- c(1, 1, 0)
inp$priors$logalpha <- c(1, 1, 0)
inp$priors$logbeta <- c(1, 1, 0)
fit.spict(inp)
   Convergence: 0  MSG: relative convergence (4)
   Objective function at optimum: 5.0598269
   Euler time step (years):  1/16 or 0.0625
   Nobs C: 23,  Nobs I1: 23
   
   No priors are used
   
   Model parameter estimates w 95% CI 
              estimate       cilow        ciupp    log.est  
    alpha   39.0515852   0.0402403 3.789800e+04  3.6648835  
    beta     0.0245072   0.0000265 2.268976e+01 -3.7087885  
    r        0.1955746   0.0313382 1.220535e+00 -1.6318133  
    rc       1.2604846   0.0097773 1.625003e+02  0.2314962  
    rold     0.2835716   0.0024030 3.346293e+01 -1.2602908  
    m       24.3112585  12.0983558 4.885270e+01  3.1909396  
    K      210.4517437 123.9798143 3.572351e+02  5.3492564  
    q        0.3842439   0.2070728 7.130026e-01 -0.9564777  
    n        0.3103166   0.0004171 2.308856e+02 -1.1701624  
    sdb      0.0028039   0.0000029 2.728242e+00 -5.8767274  
    sdf      0.3809117   0.2827823 5.130932e-01 -0.9651878  
    sdi      0.1094986   0.0812454 1.475768e-01 -2.2118440  
    sdc      0.0093351   0.0000103 8.474456e+00 -4.6739763  
    
   Deterministic reference points (Drp)
            estimate      cilow      ciupp   log.est  
    Bmsyd 38.5744642  0.5970162 2492.37666  3.652591  
    Fmsyd  0.6302423  0.0048887   81.25016 -0.461651  
    MSYd  24.3112585 12.0983558   48.85270  3.190940  
   Stochastic reference points (Srp)
            estimate      cilow      ciupp    log.est  rel.diff.Drp  
    Bmsys 38.5743443  0.5970258 2492.32112  3.6525874 -3.108375e-06  
    Fmsys  0.6302434  0.0048887   81.24967 -0.4616493  1.767620e-06  
    MSYs  24.3112241 12.0982041   48.85317  3.1909381 -1.414970e-06  
   
   States w 95% CI (inp$msytype: s)
                     estimate      cilow      ciupp    log.est  
    B_1989.94      52.6709906 29.3008926 94.6808445  3.9640648  
    F_1989.94       0.4858021  0.2429986  0.9712142 -0.7219539  
    B_1989.94/Bmsy  1.3654410  0.0268855 69.3470419  0.3114774  
    F_1989.94/Fmsy  0.7708167  0.0077214 76.9497275 -0.2603047  
   
   Predictions w 95% CI (inp$msytype: s)
                   prediction      cilow     ciupp    log.est  
    B_1991.00      51.2819174 28.0681883 93.694507  3.9373382  
    F_1991.00       0.4858023  0.1724966  1.368165 -0.7219535  
    B_1991.00/Bmsy  1.3294307  0.0231032 76.499549  0.2847508  
    F_1991.00/Fmsy  0.7708171  0.0072436 82.025622 -0.2603042  
    Catch_1990.00  25.2158724 15.5409044 40.913978  3.2274737  
    E(B_inf)       49.5039157         NA        NA  3.9020518

list.possible.priors()
    [1] "logn"      "logalpha"  "logbeta"   "logr"      "logK"      "logm"     
    [7] "logq"      "logqf"     "logbkfrac" "logB"      "logF"      "logBBmsy" 
   [13] "logFFmsy"  "logsdb"    "logsdf"    "logsdi"    "logsde"    "logsdc"   
   [19] "logsdm"    "logpsi"    "mu"        "BmsyB0"    "logngamma"

inp <- pol$albacore
inp$priors$logK <- c(log(300), 2, 1)
fit.spict(inp)
   Convergence: 0  MSG: relative convergence (4)
   Objective function at optimum: 3.6972089
   Euler time step (years):  1/16 or 0.0625
   Nobs C: 23,  Nobs I1: 23
   
   Priors
        logK  ~  dnorm[log(300), 2^2]
        logn  ~  dnorm[log(2), 2^2]
    logalpha  ~  dnorm[log(1), 2^2]
     logbeta  ~  dnorm[log(1), 2^2]
   
   Model parameter estimates w 95% CI 
              estimate       cilow       ciupp    log.est  
    alpha    8.5541387   1.2276474  59.6044821  2.1464152  
    beta     0.1213066   0.0180843   0.8137052 -2.1094339  
    r        0.2543932   0.0999840   0.6472628 -1.3688740  
    rc       0.7408800   0.1423292   3.8565754 -0.2999166  
    rold     0.8120644   0.0018385 358.6840858 -0.2081756  
    m       22.5701062  17.0284024  29.9152957  3.1166263  
    K      202.2161203 138.7005283 294.8176176  5.3093370  
    q        0.3499363   0.1930357   0.6343668 -1.0500041  
    n        0.6867327   0.0623534   7.5633733 -0.3758102  
    sdb      0.0127865   0.0018399   0.0888582 -4.3593675  
    sdf      0.3672884   0.2672952   0.5046883 -1.0016080  
    sdi      0.1093773   0.0808917   0.1478939 -2.2129522  
    sdc      0.0445545   0.0073421   0.2703737 -3.1110419  
    
   Deterministic reference points (Drp)
          estimate      cilow      ciupp    log.est  
    Bmsyd 60.92783 15.2916983 242.759230  4.1096901  
    Fmsyd  0.37044  0.0711646   1.928288 -0.9930638  
    MSYd  22.57011 17.0284024  29.915296  3.1166263  
   Stochastic reference points (Srp)
            estimate      cilow      ciupp   log.est  rel.diff.Drp  
    Bmsys 60.9201702 15.2918949 242.695046  4.109564 -1.257930e-04  
    Fmsys  0.3704521  0.0711577   1.928599 -0.993031  3.268228e-05  
    MSYs  22.5680073 17.0228324  29.919519  3.116533 -9.300394e-05  
   
   States w 95% CI (inp$msytype: s)
                     estimate      cilow       ciupp    log.est  
    B_1989.94      56.9259673 30.1898050 107.3397379  4.0417516  
    F_1989.94       0.4446516  0.2131880   0.9274211 -0.8104642  
    B_1989.94/Bmsy  0.9344355  0.2946969   2.9629416 -0.0678127  
    F_1989.94/Fmsy  1.2002944  0.2841527   5.0701835  0.1825668  
   
   Predictions w 95% CI (inp$msytype: s)
                   prediction      cilow      ciupp    log.est  
    B_1991.00      54.5233252 27.9289532 106.441261  3.9986286  
    F_1991.00       0.4446518  0.1564571   1.263702 -0.8104637  
    B_1991.00/Bmsy  0.8949963  0.2487066   3.220736 -0.1109357  
    F_1991.00/Fmsy  1.2002950  0.2373802   6.069202  0.1825674  
    Catch_1990.00  24.7355104 15.3287800  39.914819  3.2082399  
    E(B_inf)       50.1633799         NA         NA  3.9152853

inp <- pol$albacore
inp$priors$logB <- c(log(80), 0.1, 1, 1980)
par(mfrow=c(1, 2), mar=c(5, 4.1, 3, 4))
plotspict.biomass(fit.spict(pol$albacore), ylim=c(0, 500))
plotspict.biomass(fit.spict(inp), qlegend=FALSE, ylim=c(0, 500))

set.seed(123)
inp <- list(timeC=pol$albacore$timeC, obsC=pol$albacore$obsC)
inp$timeI <- list(pol$albacore$timeI, pol$albacore$timeI[10:23]+0.25)
inp$obsI <- list()
inp$obsI[[1]] <- pol$albacore$obsI * exp(rnorm(23, sd=0.1)) # Index 1
inp$obsI[[2]] <- 10*pol$albacore$obsI[10:23] # Index 2
inp$priors$logsdi <- list(c(1, 1, 0),          # No prior for index 1
                          c(log(0.1), 0.2, 1)) # Set a prior on index 2
res <- fit.spict(inp)
sumspict.priors(res)
   
   Priors
        logn  ~  dnorm[log(2), 2^2]
    logalpha  ~  dnorm[log(1), 2^2]
     logbeta  ~  dnorm[log(1), 2^2]
     logsdi2  ~  dnorm[log(0.1), 0.2^2]

inp <- pol$albacore
inp$priors$logn <- c(log(2), 1e-3)
inp$priors$logalpha <- c(log(1), 1e-3)
inp$priors$logbeta <- c(log(1), 1e-3)
fit.spict(inp)
   Convergence: 0  MSG: relative convergence (4)
   Objective function at optimum: -13.3777234
   Euler time step (years):  1/16 or 0.0625
   Nobs C: 23,  Nobs I1: 23
   
   Priors
        logn  ~  dnorm[log(2), 0.001^2] (fixed)
    logalpha  ~  dnorm[log(1), 0.001^2] (fixed)
     logbeta  ~  dnorm[log(1), 0.001^2] (fixed)
   
   Model parameter estimates w 95% CI 
              estimate      cilow       ciupp    log.est  
    alpha    1.0000039  0.9980459   1.0019658  0.0000039  
    beta     0.9999976  0.9980396   1.0019595 -0.0000024  
    r        0.5046213  0.1875615   1.3576486 -0.6839470  
    rc       0.5046222  0.1875615   1.3576536 -0.6839452  
    rold     0.5046231  0.1875608   1.3576639 -0.6839434  
    m       22.0086911 17.0614719  28.3904275  3.0914374  
    K      174.4568977 74.7565781 407.1241615  5.1616777  
    q        0.3549422  0.1279007   0.9850140 -1.0358004  
    n        1.9999964  1.9960803   2.0039201  0.6931454  
    sdb      0.0966006  0.0704060   0.1325409 -2.3371705  
    sdf      0.2102260  0.1562934   0.2827692 -1.5595724  
    sdi      0.0966010  0.0704064   0.1325411 -2.3371666  
    sdc      0.2102255  0.1562931   0.2827683 -1.5595748  
    
   Deterministic reference points (Drp)
            estimate      cilow       ciupp   log.est  
    Bmsyd 87.2283874 37.3782180 203.5621807  4.468530  
    Fmsyd  0.2523111  0.0937808   0.6788268 -1.377092  
    MSYd  22.0086911 17.0614719  28.3904275  3.091437  
   Stochastic reference points (Srp)
            estimate      cilow       ciupp   log.est rel.diff.Drp  
    Bmsys 86.1721719 37.1712866 199.7682591  4.456347 -0.012257037  
    Fmsys  0.2500268  0.0921878   0.6781096 -1.386187 -0.009136251  
    MSYs  21.5429415 16.5473970  28.0466063  3.070048 -0.021619589  
   
   States w 95% CI (inp$msytype: s)
                     estimate      cilow      ciupp    log.est  
    B_1989.94      54.8231823 17.4451855 172.287152  4.0041131  
    F_1989.94       0.4313749  0.1433961   1.297694 -0.8407778  
    B_1989.94/Bmsy  0.6362052  0.3800251   1.065080 -0.4522342  
    F_1989.94/Fmsy  1.7253144  0.9443708   3.152056  0.5454093  
   
   Predictions w 95% CI (inp$msytype: s)
                   prediction      cilow      ciupp    log.est  
    B_1991.00      50.1747639 14.1661714 177.712585  3.9155122  
    F_1991.00       0.4313751  0.1324968   1.404446 -0.8407773  
    B_1991.00/Bmsy  0.5822618  0.2912706   1.163965 -0.5408351  
    F_1991.00/Fmsy  1.7253153  0.8254232   3.606287  0.5454098  
    Catch_1990.00  22.6023879 14.8442486  34.415210  3.1180556  
    E(B_inf)       20.7048450         NA         NA  3.0303677

Robust estimation (reducing influence of extreme observations)

The presence of extreme observations may inflate estimates of observation noise and increase the general uncertainty of the fit. To reduce this effect it is possible to apply a robust estimation scheme, which is less sensitive to extreme observations. An example with an extreme observation in the catch series is

inp <- pol$albacore
inp$obsC[10] <- 3*inp$obsC[10]
res1 <- fit.spict(inp)
inp$robflagc <- 1
res2 <- fit.spict(inp)
sumspict.parest(res2)
              estimate        cilow        ciupp    log.est
   alpha    8.01383384   1.14068244  56.30097453  2.0811693
   beta     0.13903418   0.02002496   0.96532028 -1.9730355
   r        0.25685421   0.10125244   0.65158021 -1.3592466
   rc       0.73334983   0.13979381   3.84710861 -0.3101324
   rold     0.85759781   0.00140626 522.99835534 -0.1536200
   m       22.57344297  16.94368700  30.07375710  3.1167741
   K      202.06196124 137.06165230 297.88810726  5.3085744
   q        0.34766877   0.18846324   0.64136419 -1.0565050
   n        0.70049573   0.06409210   7.65608050 -0.3559670
   sdb      0.01371149   0.00196492   0.09568068 -4.2895209
   sdf      0.37086361   0.26568724   0.51767565 -0.9919209
   sdi      0.10988163   0.08106946   0.14893367 -2.2083516
   sdc      0.05156272   0.00843522   0.31519203 -2.9649564
   pp       0.95304961   0.72887564   0.99351802  3.0105755
   robfac  20.83563432   2.67338471 236.12369134  2.9874800
par(mfrow=c(1, 2))
plotspict.catch(res1, main='Regular fit')
plotspict.catch(res2, qlegend=FALSE, main='Robust fit')

It is evident from the plot that the presence of the extreme catch observation generally inflates the uncertainty of the estimated catches, while the robust fit is less sensitive. Robust estimation can be applied to index and effort data using robflagi and robflage respectively.

Robust estimation is implemented using a mixture of light-tailed and a heavy-tailed Gaussian distribution as described in Pedersen and Berg (2017). This entails two additional parameters (pp and robfac) that require estimation. This may not always be possible given the increased model complexity. In such cases these parameters should be fixed by setting their phases to -1.

Forecasting and management scenarios

All quantities estimated by SPiCT can be forecasted into the future. By default, the forecast period starts one time step after the last observation, which can be defined by the last index observation or by the end of the last catch or effort interval. Forecasting catch, fishing mortality, and biomass is a crucial step in fisheries management that allows determining future yields and risk of overfishing under current or alternative management scenarios. The following provides details to the forecasting and management functionality of SPiCT and includes examples based on the South Atlantic albacore data set of Polacheck, Hilborn, and Punt (1993), which contains catch and commercial CPUE data from 1967 to 1990.

inp <- pol$albacore
inp$maninterval <- c(1990, 1991)
inp$maneval <- 1991
inp$ffac <- 0.75
rep <- fit.spict(inp)

sumspict.predictions(rep)
                  prediction      cilow      ciupp     log.est
   B_1991.00      59.6456548 32.5515458 109.291404  4.08842130
   F_1991.00       0.3348376  0.1179796   0.950302 -1.09410958
   B_1991.00/Bmsy  0.9820378  0.2791407   3.454882 -0.01812545
   F_1991.00/Fmsy  0.9006337  0.1793004   4.523921 -0.10465670
   Catch_1990.00  19.4447465 11.7527439  32.171055  2.96757693
   E(B_inf)       67.1855141         NA         NA  4.20745766

plot2(rep)

inp <- pol$albacore
rep <- fit.spict(inp)
rep <- manage(rep)
   Selected scenario(s):  currentCatch, currentF, Fmsy, noF, reduceF25, increaseF25, msyHockeyStick, ices

sumspict.manage(rep)
   SPiCT timeline:
                                                     
         Observations              Management        
       1967.00 - 1990.00        1990.00 - 1991.00    
    |-----------------------| ----------------------|
   
   Management evaluation: 1991.00
   
   Predicted catch for management period and states at management evaluation time:
   
                               C B/Bmsy F/Fmsy
   1. Keep current catch    25.2   0.89   1.23
   2. Keep current F        24.7   0.89   1.20
   3. Fish at Fmsy          21.3   0.95   1.00
   4. No fishing             0.0   1.30   0.00
   5. Reduce F by 25%       19.4   0.98   0.90
   6. Increase F by 25%     29.5   0.81   1.50
   7. MSY hockey-stick rule 21.3   0.95   1.00
   8. ICES advice rule      19.3   0.98   0.89

plot2(rep)

plotspict.hcr(rep)

man.timeline(inp)
   SPiCT timeline:
                                                     
         Observations              Management        
       1967.00 - 1990.00        1990.00 - 1991.00    
    |-----------------------| ----------------------|
   
   Management evaluation: 1991.00

man.timeline(rep)
   SPiCT timeline:
                                                     
         Observations              Management        
       1967.00 - 1990.00        1990.00 - 1991.00    
    |-----------------------| ----------------------|
   
   Management evaluation: 1991.00

inp$maninterval <- c(1991,1992)
man.timeline(inp)
   SPiCT timeline:
                                                                             
         Observations             Intermediate             Management        
       1967.00 - 1990.00        1990.00 - 1991.00       1991.00 - 1992.00    
    |-----------------------| ----------------------| ----------------------|
   
   Management evaluation: 1992.00

repIntPer <- manage(rep, scenarios = c(1,2),
                    maninterval = c(1991,1992), maneval = 1992)
   Selected scenario(s):  currentCatch, currentF
plotspict.catch(repIntPer)

repIntPerC <- manage(rep, scenarios=c(8), maninterval = c(1991,1992), intermediatePeriodCatch = 20)
   Selected scenario(s):  ices
par(mfrow=c(1,2))
plotspict.biomass(repIntPerC)
plotspict.catch(repIntPerC)

repIntPer <- add.man.scenario(repIntPer, ffac = 1.5)
sumspict.manage(repIntPer)
   SPiCT timeline:
                                                                             
         Observations             Intermediate             Management        
       1967.00 - 1990.00        1990.00 - 1991.00       1991.00 - 1992.00    
    |-----------------------| ----------------------| ----------------------|
   
   Management evaluation: 1992.00
   
   Predicted catch for management period and states at management evaluation time:
   
                                 C B/Bmsy F/Fmsy
   1. Keep current catch (@$) 25.2   0.85   1.28
   2. Keep current F (@$)     23.9   0.87   1.20
   3. customScenario_1        33.9   0.74   1.80
   (@) This scenario assumes another management period. Thus, the estimates might not be comparable to the other scenarios.
   ($) This scenario assumes another management evaluation time. Thus, the estimates might not be comparable to the other scenarios.
plotspict.f(repIntPer)

repIntPer <- add.man.scenario(repIntPer, cfac = 0.64,
                              scenarioTitle = "reduced_catch")
names(repIntPer$man)
   [1] "currentCatch"     "currentF"         "customScenario_1" "reduced_catch"

repIntPer <- add.man.scenario(repIntPer, scenarioTitle = "ices",
                              breakpointB = 0.5,
                              fractiles = list(catch=0.35, bbmsy=0.35, ffmsy=0.35))
sumspict.manage(repIntPer)
   SPiCT timeline:
                                                                             
         Observations             Intermediate             Management        
       1967.00 - 1990.00        1990.00 - 1991.00       1991.00 - 1992.00    
    |-----------------------| ----------------------| ----------------------|
   
   Management evaluation: 1992.00
   
   Predicted catch for management period and states at management evaluation time:
   
                                 C B/Bmsy F/Fmsy
   1. Keep current catch (@$) 25.2   0.85   1.28
   2. Keep current F (@$)     23.9   0.87   1.20
   3. customScenario_1        33.9   0.74   1.80
   4. reduced_catch           16.2   1.04   0.73
   5. ICES advice rule        15.0   1.06   0.67
   (@) This scenario assumes another management period. Thus, the estimates might not be comparable to the other scenarios.
   ($) This scenario assumes another management evaluation time. Thus, the estimates might not be comparable to the other scenarios.

repVarPer <- manage(rep, scenarios = c(1,4), maninterval = c(1990.5,1991))
   Selected scenario(s):  currentCatch, noF

repVarPer <- add.man.scenario(repVarPer, maninterval = c(1991, 1993),
                              breakpointB = 0.5)
plotspict.f(repVarPer)

sumspict.manage(repVarPer)
   SPiCT timeline:
                                                                             
         Observations             Intermediate             Management        
       1967.00 - 1990.00        1990.00 - 1990.50       1990.50 - 1991.00    
    |-----------------------| ----------------------| ----------------------|
   
   Management evaluation: 1991.00
   
   Predicted catch for management period and states at management evaluation time:
   
                                C B/Bmsy F/Fmsy
   1. Keep current catch (@) 12.6   0.89   1.24
   2. No fishing (@)          0.0   1.10   0.00
   3. customScenario_1 (@)   41.7   0.89   1.00
   (@) This scenario assumes another management period. Thus, the estimates might not be comparable to the other scenarios.

names(rep$man)
   [1] "currentCatch"   "currentF"       "Fmsy"           "noF"           
   [5] "reduceF25"      "increaseF25"    "msyHockeyStick" "ices"
rep <- man.select(rep, scenarios = c(1,3))
   Selected scenario(s):  currentCatch, Fmsy

man.tac(rep)
   $currentCatch
   [1] 25.23699
   
   $Fmsy
   [1] 21.25489

get.TAC(rep, maninterval = c(1990.25, 1991),
        safeguardB = list(limitB = 0.3), verbose = TRUE)
   Note that the specified management interval: 1990.25 - 1991.00 differs from existing scenarios in rep$man.
   [1] 15.73659

## selection by index
length(repIntPer$man)
   [1] 5
repSelect1 <- man.select(repIntPer, scenarios = c(3,1))
   Selected scenario(s):  currentCatch, customScenario_1

## selection by scenario name
names(repIntPer$man)
   [1] "currentCatch"     "currentF"         "customScenario_1" "reduced_catch"   
   [5] "ices"
repSelect2 <- man.select(repIntPer, scenarios = c("currentCatch"))
   Selected scenario(s):  currentCatch

## selected object as standard spictcls object
repSelect3 <- man.select(repIntPer, scenarios = c("currentCatch"), spictcls = TRUE)

## Plot objects with selected scenarios
opar <- par(mfrow=c(3,1))
plotspict.catch(repSelect1)
plotspict.catch(repSelect2)
plotspict.catch(repSelect3)
par(opar)

rep <- man.select(rep, scenarios = NULL)
   All scenarios removed!

catchunit - Define unit of catch observations

This will print the unit of the catches on relevant plots.

Example: inp$catchunit <- "'000 t".

dteuler and eulertype - Temporal discretisation and time step

To solve the continuous-time system an Euler discretisation scheme is used. This requires a time step to be specified (dteuler). The smaller the time step the more accurate the approximation to the continuous-time solution, however with the cost of increased memory requirements and computing time. The default value of dteuler is 1/16, which seems sufficiently fine for most cases, and perhaps too fine for some cases. When fitting quarterly data and species with fast growth it is important to have a small dteuler. The influence of dteuler on the results can be checked by using different values and comparing resulting model estimates. If dteuler <- 1 the model essentially becomes a discrete-time model with one Euler step per year.

There are two possible temporal discretisation schemes which can be set to either eulertype = 'hard' (default) or eulertype = 'soft'. If eulertype = 'hard' then time is discretised into intervals of length dteuler. Observations are then assigned to these intervals. For annual and quarterly data dteuler = 1/16 is appropriate, however if fitting to monthly data dteuler should be changed to e.g. 1/24. If eulertype = 'soft' (careful, this feature has not been thoroughly tested), then time is discretised into intervals of length dteuler and additional time points corresponding to the times of observation are added to the discretisation. This feature is particularly useful if observations (most likely index series) are observed at odd times during the year. The model then estimates values of biomass and fishing mortality at the exact time of the observation instead of assigning the observation to an interval.

msytype - Stochastic and deterministic reference points

As default the stochastic reference points are reported and used for calculation of relative levels of biomass and fishing mortality. It is, however, possible to use the deterministic reference points by setting inp$msytype <- 'd'.

do.sd.report - Perform SD report calculations

The sdreport step calculates the uncertainty of all quantities that are reported in addition to the model parameters. For long time series and with small dteuler this step may have high memory requirements and a substantial computing time. Thus, if one is only interested in the point estimates of the model parameters it is advisable to set do.sd.report <- 0 to increase speed.

reportall - Report all derived quantities

If uncertainties of some quantities (such as reference points) are required, but uncertainty on state variables (biomass and fishing mortality) are not needed, then reportall <- 0 can be used to increase speed.

optim.method - Report all derived quantities

Parameter estimation is per default performed using R’s nlminb() optimiser. Alternatively it is possible to use optim by setting inp$optim.method <- 'optim'.