Package 'spict'

Description

Check whether ACF of residuals is significant in any lags.

Usage

acf.signf(resid, lag.max = 4, return.p = FALSE)

Arguments

Details

This corresponds to plotting the ACF using acf() and checking whether any lags has an acf value above the CI limit.

Value

Vector of TRUE and FALSE indicating whether significant lags were present. If return.p is TRUE then p-values are returned instead.

Description

Add catch unit to label

Usage

add.catchunit(lab, cu)

Arguments

Value

Label with added catch unit

Description

Add a legend explaining colors of points (vertical orientation)

Usage

add.col.legend()

Value

Nothing.

Description

Add a legend explaining colors of points (horizontal orientation)

Usage

add.col.legend.hor()

Value

Nothing.

Description

Define management scenario

Usage

add.man.scenario(
  rep,
  scenarioTitle = "",
  maninterval = NULL,
  maneval = NULL,
  ffac = NULL,
  fabs = NULL,
  cfac = NULL,
  cabs = NULL,
  fractiles = list(catch = 0.5, bbmsy = 0.5, ffmsy = 0.5),
  breakpointB = 0,
  safeguardB = list(limitB = 0, prob = 0.95),
  intermediatePeriodCatch = NULL,
  intermediatePeriodCatchSDFac = 1,
  intermediatePeriodCatchList = NULL,
  ctol = 0.001,
  evalBreakpointB = 0,
  verbose = TRUE,
  dbg = 0,
  mancheck = TRUE
)

get.TAC(
  rep,
  scenarioTitle = "",
  maninterval = NULL,
  maneval = NULL,
  ffac = NULL,
  fabs = NULL,
  cfac = NULL,
  cabs = NULL,
  fractiles = list(catch = 0.5, bbmsy = 0.5, ffmsy = 0.5),
  breakpointB = 0,
  safeguardB = list(limitB = 0, prob = 0.95),
  intermediatePeriodCatch = NULL,
  intermediatePeriodCatchSDFac = 1,
  intermediatePeriodCatchList = NULL,
  ctol = 0.001,
  evalBreakpointB = 0,
  verbose = TRUE,
  dbg = 0,
  mancheck = TRUE
)

make.man.inp(
  rep,
  scenarioTitle = "",
  maninterval = NULL,
  maneval = NULL,
  ffac = NULL,
  fabs = NULL,
  cfac = NULL,
  cabs = NULL,
  fractiles = list(catch = 0.5, bbmsy = 0.5, ffmsy = 0.5),
  breakpointB = 0,
  safeguardB = list(limitB = 0, prob = 0.95),
  intermediatePeriodCatch = NULL,
  intermediatePeriodCatchSDFac = 1,
  intermediatePeriodCatchList = NULL,
  ctol = 0.001,
  evalBreakpointB = 0,
  verbose = TRUE,
  dbg = 0,
  mancheck = TRUE
)

Arguments

Details

Default management scenario

The default management scenario is fish at $F_{MSY}$. This is when ffac, cfac, fabs, cabs are all NULL, and breakpointB and safeguardB$limitB are 0. In practice ffac is set equal to $F_{MSY}/F_m$.

Catch scenarios

Management scenarios can be defined based on a desired catch during the management period. Common examples include scenarios like "increase catch by 25%", "keep current catch", or "zero catch". The catch can be relative to the predicted "previous catch", using the multiplier cfac, or in absolute terms using cabs catch in the same units as the input data. By default, the respective previous catch corresponds to that part of the previous year which corresponds to the management interval. For example, if the management period is $[1991, 1992[$, the whole catch from the year $[1990, 1991[$ is being used. If the management period is $[1991.5, 1991.75[$, the same interval from the previous year $[1990.5, 1990.75[$ is being used. If the management period spans several years, e.g. $[1991, 1993[$, the whole catch from the previous year $[1990, 1991[$ is being used two times.

Harvest Control Rules (HCRs)

The combination of the arguments "fractiles", "breakpointB", and "safeguardB" allow the specification of a number of different harvest control rules:

• MSY hockey-stick rule: Fishing at F_MSY above a certain biomass reference level (here defined as a fraction of B_MSY with breakpointB). Below the reference level, fishing is reduced linearly to 0 as suggested in ICES (2017).

• MSY (hockey-stick) rule with additional precautionary buffer: As long as the probability of the predicted biomass relative to a reference biomass level (e.g. 0.3 B_MSY, defined by safeguardB$limitB) is smaller or equal to a specified risk aversion probability (e.g. 95%, defined by safeguardB$prob), fishing at F_MSY or following the hockey-stick rule (if breakpoint != 0), otherwise reduce fishing mortality to meet specified risk aversion probability (safeguardB$prob) as introduced in ICES (2018).

• By ICES (2019) recommended MSY hockey-stick rule with 35th percentiles: Fishing at 35th percentile of F_MSY above the 35th percentile of 0.5 $B/B_{MSY}$ (breakpointB = 0.5) and 35th percentile of linearly reduced F_MSY below the 35th percentile of 0.5 $B/B_{MSY}$. TAC corresponds to 35th percentile of predicted catch. Rule is applied with fractiles = list(catch=0.35, bbmsy=0.35, ffmsy=0.35), breakpointB = 0.5, and safeguardB = list(limitB = 0, prob = 0.95).

Fractiles

By default, the median (fractile of 0.5) is used for the stock status ($B/B_{MSY}$, $F/F_{MSY}$) and predicted catch distribution. A more precautionary approach is to used fractiles lower than the median (0.5) to account for the estimated uncertainty. The arguments of the 'fractiles' are:

• catch - Fractile of the predicted catch distribution

• bbmsy - Fractile of the $B/B_{MSY}$ distribution

• ffmsy - Fractile of the $F/F_{MSY}$ distribution

Note that the fractile for the $F/F_{MSY}$ distribution is 1 minus the fractile specified. As the current fishing mortality is divided by the value of this distribution $F{y+1} = \frac{F_y}{F_y/F_{MSY}}$, a lower percentile of the $F/F_{MSY}$ distribution is more conservative than a larger one. This allows a consistent setting of fractiles among the different quantities.

Biomass safeguard

The argument list "safeguardB" includes:

• limitB - Reference level for the evaluation of the predicted biomass defined as fraction of $B/B_{MSY}$. By default (safeguardB$limitB == 0) the PA buffer is not used. Theoretically, any value smaller than 1 is meaningful, but an ICES recommended value would be 30% safeguardB$limitB = 0.3 (ICES, 2018).

• prob - Risk aversion probability of the predicted biomass relative to specified reference level (safeguardB$limitB) for all rules with PA buffer (safeguardB$limitB != 0). Default: 0.95 as recommended by ICES (2018).

Intermediate period assumptions

Dependent on the start of the management period (e.g. advice year), there might be a time lag between the last observation and the start of the management period, often referred to as the intermediate period. If this is the case, an assumption about the catch during intermediate time period (e.g. assessment year) has to be made. Two meaningful assumptions are:

• 1: The catch in the intermediate period is based on the fishing mortality which is extrapolated from the previous year. This is the default assumption;

• 2: The catch in the intermediate period is directly specified. This could for example be the TAC recommended in the previous year. The catch can be specified by means of the argument intermediatePeriodCatch. Be aware that this catch might correspond to several years or a fraction of a year depending on the time between the last observation and the start of the management period. The function man.timeline can help visualising the default or specified intermediate period in your data. The argument intermediatePeriodCatchSDFac allows to specify the factor with which to multiply the standard deviation of the catch ($\sigma_C$) with. It is thus a measure of the certainty around the catch in the intermediate period. The argument intermediatePeriodCatchList allows to define a list with catches and their intervals. It is a list with the elements 'obsC', 'timeC', 'dtc' and the optional element 'stdevfacC' (which is equal to 1 if not provided).

make.man.inp

Internal function that creates the required input list for the specific HCR.

Value

add.man.scenario returns the input object rep with the specified HCR added to the man list. get.TAC returns the total allowable catch (TAC) based on the specified scenario. make.man.inp returns the updated inp list based on specified HCR.

References

ICES. 2017. Report of the Workshop on the Development of the ICES approach to providing MSY advice for category 3 and 4 stocks (WKMSYCat34), 6-10 March 2017, Copenhagen, Denmark. ICES CM 2017/ ACOM:47. 53 pp.

ICES. 2018. Report of the Eighth Workshop on the Development of Quantitative Assessment Methodologies based on LIFE-history traits, exploitation characteristics, and other relevant parameters for data-limited stocks (WKLIFE VIII), 8-12 October 2018, Lisbon, Portugal. ICES CM 2018/ACOM:40. 172 pp.

ICES.2019. Ninth Workshop on the Development of Quantitative Assessment Methodologies based on LIFE-history traits, exploitation characteristics, and other relevant parameters for data-limited stocks (WKLIFE IX). ICES Scientific Reports. 1:77. 131 pp.http://doi.org/10.17895/ices.pub.5550

ICES 2020. Report of the Ninth Workshop on the Development of Quantitative Assessment Methodologies based on LIFE-history traits, exploitation characteristics, and other relevant parameters for data-limited stocks (WKLIFE X), ICES Scientific Reports. 2:98. 72 pp. http://doi.org/10.17895/ices.pub.5985

Examples

data(pol)
rep <- fit.spict(pol$albacore)

## Fishing at Fmsy
rep <- add.man.scenario(rep)

## MSY hockey-stick rule
rep <- add.man.scenario(rep, breakpointB = 0.5)

## ICES (2019) recommended HCR
rep <- add.man.scenario(rep, fractiles = list(catch=0.35, bbmsy=0.35, ffmsy=0.35), breakpointB=0.5)

## Get the TAC for the ICES (2020) recommended HCR (as used in WKMSYSPICT)
rep <- add.man.scenario(rep, fractiles = list(catch=0.35), breakpointB = c(0.3, 0.5))

## Now `rep` includes 3 management scenarios

## Get the TAC when fishing mortality is equal to Fmsy
get.TAC(rep)

## Get TAC for the MSY hockey-stick rule (only using Btrigger)
get.TAC(rep, breakpointB = 0.5)

## Get TAC for the MSY hockey-stick rule (with Btrigger and Blim)
get.TAC(rep, breakpointB = c(0.3, 0.5))

## Get the TAC for the ICES (2019) recommended HCR
get.TAC(rep, fractiles = list(catch=0.35, bbmsy=0.35, ffmsy=0.35), breakpointB=0.5)

## Get the TAC for the ICES (2020) recommended HCR (as used in WKMSYSPICT)
get.TAC(rep, fractiles = list(catch=0.35), breakpointB = c(0.3, 0.5))

Description

Add lines to plot indicating result of management scenarios.

Usage

add.manlines(
  rep,
  par,
  par2 = NULL,
  index.shift = 0,
  plot.legend = TRUE,
  verbose = TRUE,
  ...
)

Arguments

Value

Nothing

Description

Convert from quarterly (or other sub-annual) data to annual means, sums or a custom function.

Usage

annual(intime, vec, type = "mean")

Arguments

Value

A list containing the annual means $annvec and a corresponding time vector $anntime.

Description

Draw a line with arrow heads.

Usage

arrow.line(
  x,
  y,
  length = 0.25,
  angle = 30,
  code = 2,
  col = par("fg"),
  lty = par("lty"),
  lwd = par("lwd"),
  ...
)

Arguments

Details

Add to an existing plot a continuous line with arrow heads showing the direction between each data point

Value

Nothing, but an arrow line is added to the current plot.

Description

Calculates the Bmsy/K ratio

Usage

calc.bmsyk(rep)

Arguments

Value

Bmsy/K

Description

Calculate E(Binfinity), i.e. the fished equilibrium.

Usage

calc.EBinf(K, n, Fl, Fmsy, sdb2)

Arguments

Details

If a seasonal pattern in F is imposed the annual average F is used for calculating the expectation. Max() is used to avoid negative values.

Value

E(Binf).

Description

Calculate gamma from n

Usage

calc.gamma(n)

Arguments

Description

Calculates influence statistics of observations.

Usage

calc.influence(rep, mc.cores = 1)

Arguments

Details

TBA

Value

A list equal to the input with the added key "infl" containing influence statistics.

Description

Calculate Mean Absolute Scaled Error (MASE)

Usage

calc.mase(rep, verbose = TRUE)

Arguments

Details

This function calculates the Mean Absolute Scaled Error (MASE) for each index time series by hindcasting. Thus, the application of this method requires a fitted spict object with the results of the hindcasting analysis hindcast.

The smaller MASE, the higher the predictive power of the spict model regarding the prediction of index observations. In contrast, a MASE above 1 suggests that the naive prediction of the index observations assuming the preceding index observations have a higher predictive power than the spict model. Note, however, that the absolute MASE value depends on multiple factors such as the number of peels, assumed priors, etc.

Note that a difference in the timing of the index observations of less than a month are considered acceptable for the estimation of the naive prediction residuals and no warning is printed. If the variability exceeds a month, the predictions are still calculated, but a warning is printed.

Value

A data frame with estimate the MASE and the number of runs used for the estimation for each index.

References

Carvalho, F., Winker, H., Courtney, D., Kapur, M., Kell, L., Cardinale, M., Schirripa, M., Kitakado, T., Yemane, D., Piner, K.R. Maunder, M.N., Taylor, I., Wetzel, C.R., Doering, K., Johnsonm, K.F., Methot, R. D. (2021). A cookbook for using model diagnostics in integrated stock assessments. Fisheries Research, 240, 105959.

Kell, L. T., Kimoto, A., & Kitakado, T. (2016). Evaluation of the prediction skill of stock assessment using hindcasting. Fisheries research, 183, 119-127.

Kell, L. T., Sharma, R., Kitakado, T., Winker, H., Mosqueira, I., Cardinale, M., & Fu, D. (2021). Validation of stock assessment methods: is it me or my model talking?. ICES Journal of Marine Science, 78(6), 2244-2255.

Winker, H., Carvalho, F., & Kapur, M. (2018). JABBA: just another Bayesian biomass assessment. Fisheries Research, 204, 275-288.

Examples

data(pol)
inp <- pol$albacore
rep <- fit.spict(inp)
rep <- hindcast(rep)
calc.mase(rep)

Description

Calculates the order of magnitude for the relative reference levels B/Bmsy and F/Fmsy

Usage

calc.om(rep, CI = 0.95)

Arguments

Details

The lower, upper values and the CI range are based on the 95% confidence interval (CI; default).

Value

Matrix containing the order of magnitude for B/Bmsy and F/Fmsy.

Description

Calculate one-step-ahead residuals.

Usage

calc.osa.resid(rep)

Arguments

Details

In TMB one-step-ahead residuals are calculated by sequentially including one data point at a time while keeping the model parameters fixed at their ML estimates. The calculated residuals are tested for independence, bias, and normality.

Value

An updated result report, which contains one-step-ahead residuals stored in $osarC and $osarI.

Examples

data(pol)
rep <- fit.spict(pol$albacore)
rep <- calc.osa.resid(rep)
plotspict.osar(rep)

Description

Calculate process residuals

Usage

calc.process.resid(rep, dt = NULL)

Arguments

Details

Calculates process residuals for biomass and fishing mortality process.

Value

Data frame with year, process residuals for biomass and fishing mortality in the columns.

Examples

data(pol)
inp <- pol$albacore
rep <- fit.spict(inp)
rep <- calc.process.resid(rep)
plotspict.diagnostic.process(rep)

Description

Calculate Total Allowable Catch (TAC)

Usage

calc.tac(rep, inp = NULL, fractileCatch = 0.5, exp = TRUE)

Arguments

Value

Total Allowable Catch (TAC)

Description

Check catch list

Usage

check.catchList(catchList, sdfac = 1)

Arguments

Details

Internal function that checks if catchList is complete.

Value

Checked catchList

Description

Check sensitivity of fit to initial parameter values

Usage

check.ini(input, ntrials = 10, verbose = TRUE, numdigits = 2)

Arguments

Value

List containing results of sensitivity check and associated initial values.

Description

Check list of input variables

Usage

check.inp(inp, verbose = TRUE, mancheck = TRUE)

Arguments

Details

Fills in default values if missing.

Required inputs:

• "inp$obsC" Vector of catch observations.

• "inp$obsI and/or inp$obsE" List containing vectors of index observations and/or a vector of effort information.

Optional inputs:

- Data

• "inp$timeC" Vector of catch times. Default: even time steps starting at 1.

• "inp$timeI" List containing vectors of index times. Default: even time steps starting at 1.

• "inp$timeE" Vector of effort times. Default: even time steps starting at 1.

• "inp$dtc" Time interval for catches, e.g. for annual catches inp$dtc=1, for quarterly catches inp$dtc=0.25. Can be given as a scalar, which is then used for all catch observations. Can also be given as a vector specifying the catch interval of each catch observation. Default: min(diff(inp$timeC)).

• "inp$dte" Time interval for effort observations. For annual effort inp$dte=1, for quarterly effort inp$dte=0.25. Default: min(diff(inp$timeE)).

• "inp$nseasons" Number of within-year seasons in data. If inp$nseasons > 1 then a seasonal pattern is used in F. Valid values of inp$nseasons are 1, 2 or 4. Default: number of unique within-year time points present in data.

• "start.in.first.data.point" Logical. If TRUE (default), modelling time starts at the first available data point, otherwise it starts in the beginning of that year.

- Initial parameter values

• "inp$ini$logn" Pella-Tomlinson exponent determining shape of production function. Default: log(2) corresponding to the Schaefer formulation.

• "inp$ini$logm" Initial value for logm (log maximum sustainable yield). Default: log(mean(catch)).

• "inp$ini$logK" Initial value for logK (log carrying capacity). Default: log(4*max(catch)).

• "inp$ini$logq" Initial value for logq (log catchability of index). Default: log(max(index)/K).

• "inp$ini$logsdb" Initial value for logsdb (log standard deviation of biomass process). Default: log(0.2).

• "inp$ini$logsdf" Initial value for logsdf (log standard deviation of fishing mortality process). Default: log(0.2).

• "inp$ini$logsdi" Initial value for logsdi (log standard deviation of index observation error). Default: log(0.2).

• "inp$ini$logsdc" Initial value for logsdc (log standard deviation of catch observation error). Default: log(0.2).

• "inp$ini$phi" Vector for cyclic B spline representing within-year seasonal variation. Default: rep(1, inp$nseasons).

• "inp$ini$logsdu" Initial value for logsdu (log standard deviation of log U, the state of the coupled SDE representation of seasonality). Default: log(0.1).

• "inp$ini$loglambda" Initial value for loglambda (log damping parameter of the coupled SDE representation of seasonality). Default: log(0.1).

- Initial values for unobserved states estimated as random effects

• "inp$ini$logF" Log fishing mortality. Default: log(0.2*r), with r derived from m and K.

• "inp$ini$logB" Log biomass. Default: log(0.5*K).

• "inp$ini$logU" Log U, the state of the coupled SDE representation of seasonality. Default: log(1).

- Priors

Priors on model parameters are assumed generally assumed Gaussian and specified in a vector of length 2: c(log(mean), stdev in log domain, useflag [optional]). NOTE: if specifying a prior for a value in a temporal vector e.g. logB, then a fourth element is required specifying the year the prior should be applied. log(mean): log of the mean of the prior distribution. stdev in log: standard deviation of the prior distribution in log domain. useflag: if 1 then the prior is used, if 0 it is not used. Default is 1. To list parameters to which priors can be applied run list.possible.priors(). Example: intrinsic growth rate of 0.8 inp$priors$logr <- c(log(0.8), 0.1) inp$priors$logr <- c(log(0.8), 0.1, 1) # This includes the optional useflag Example: Biomass prior of 200 in 1985 inp$priors$logB <- c(log(200), 0.2, 1985) inp$priors$logB <- c(log(200), 0.2, 1, 1985) # This includes the optional useflag

- Settings/Options/Preferences

• "inp$maninterval" Start and end time of management period. Default: One year interval starting at the beginning of the new year after the last observation. Example: inp$maninterval <- c(2020.25,2021.25)

• "inp$maneval" Time for the estimation of predicted model states (biomass and fishing mortality), which can be used to evaluate the implications of management scenarios. Default: At the end of the management interval inp$maninterval[2]. Example: inp$maneval <- 2021.25

• "inp$timepredc" Deprecated: Predict accummulated catch in the interval starting at $timepredc and $dtpredc into the future. Default depends on inp$maninterval.

• "inp$dtpredc" Deprecated: Length of catch prediction interval in years. Default depends on inp$maninterval.

• "inp$timepredi" Deprecated: Predict index until this time. Default depends on inp$maneval.

• "inp$manstart" Deprecated: Start of the management period. Updated argument inp$maninterval. Default depends on inp$maninterval.

• "inp$do.sd.report" Flag indicating whether SD report (uncertainty of derived quantities) should be calculated. For small values of inp$dteuler this may require a lot of memory. Default: TRUE.

• "inp$reportall" Flag indicating whether quantities derived from state vectors (e.g. B/Bmsy, F/Fmsy etc.) should be calculated by SD report. For small values of inp$dteuler (< 1/32) reporting all may have to be set to FALSE for sdreport to run. Additionally, if only reference points of parameter estimates are of interest one can set to FALSE to gain a speed-up. Default: TRUE.

• "inp$reportmode" Integer between 0 and 2 determining which objects will be adreported. Default: 0 = all quantities are adreported. Example: inp$reportmode <- 1

• "inp$reportRel" Flag indicating whether mean 1 standardized states (i.e. B/mean(B), F/mean(F) etc.) should be calculated by SD report. Default: FALSE.

• "inp$robflagc" Flag indicating whether robust estimation should be used for catches (either 0 or 1). Default: 0.

• "inp$robflagi" Vector of flags indicating whether robust estimation should be used for indices (either 0 or 1). Default: 0.

• "inp$ffac" Management scenario represented by a factor to multiply F with when calculating the F of the next time step. ffac=0.8 means a 20% reduction in F over the next year. The factor is only used when predicting beyond the data set. Default: 1 (0% reduction).

• "inp$dteuler" Length of Euler time step in years. Default: 1/16 year.

• "inp$phases" Phases can be used to fix/free parameters and estimate in different stages or phases. To fix e.g. logr at inp$ini$logr set inp$phases$logr <- -1. To free logalpha and estimate in phase 1 set inp$phases$logalpha <- 1.

• "inp$osar.method" Method to use in TMB's oneStepPredict function. Valid methods include: "oneStepGaussianOffMode", "fullGaussian", "oneStepGeneric", "oneStepGaussian", "cdf". See TMB help for more information. Default: "none" (i.e. don't run this).

• "inp$osar.trace" If TRUE print OSAR calculation progress to screen. Default: FALSE.

• "inp$osar.parallel" If TRUE parallelise OSAR calculation for speed-up. Default: FALSE.

• "inp$catchunit" Specify unit of catches to be used in plotting legends. Default: ”.

• "inp$stdevfacC" Factors to multiply the observation error standard deviation of each individual catch observation. Can be used if some observations are more uncertain than others. Must be same length as observation vector. Default: 1.

• "inp$stdevfacI" Factors to multiply the observation error standard deviation of each individual index observation. Can be used if some observations are more uncertain than others. A list with vectors of same length as observation vectors. Default: 1.

• "inp$stdevfacE" Factors to multiply the observation error standard deviation of each individual effort observation. Can be used if some observations are more uncertain than others. A list with vectors of same length as observation vectors. Default: 1.

• "inp$mapsdi" Vector of length equal to the number of index series specifying which indices that should use the same sdi. For example: in case of 3 index series use inp$mapsdi <- c(1, 1, 2) to have series 1 and 2 share sdi and have a separate sdi for series 3. Default: 1:nindex, where nindex is number of index series.

• "inp$seasontype" If set to 1 use the spline-based representation of seasonality. If set to 2 use the oscillatory SDE system (this is more unstable and difficult to fit, but also more flexible).

• "inp$sim.random.effects"Should random effects (logB, logF, etc.) be simulated (default) or the same random effects be used (as specified in inp$ini or in an fitted spict object)?

• "inp$sim.fit"Should the estimated parameters from the last fit of a fitted spict object be used for simulation (env$last.par, default) or the inital values (specified in inp$ini)?. Note, that this only works if a fitted spict object is provided as input to sim.spict.

Value

An updated list of input variables checked for consistency and with defaults added.

Examples

data(pol)
(inp <- check.inp(pol$albacore))

Description

Check the consistency of management scenarios in rep

Usage

check.man(
  rep,
  maninterval = NULL,
  maneval = NULL,
  verbose = TRUE,
  reportmode0 = TRUE
)

Arguments

Details

Internal function that checks if the fitted spict objects in rep$man have a consistent management interval.

Value

TRUE/FALSE

Description

Checks and corrects management time to be within model time

Usage

check.man.time(
  x,
  maninterval = NULL,
  maneval = NULL,
  verbose = TRUE,
  printTimeline = FALSE,
  mancheck = TRUE
)

Arguments

Value

Updated input list or fitted spict object dependent on type of input.

Examples

data(pol)
inp <- check.inp(pol$albacore)
rep <- fit.spict(inp)

## with an input list
check.man.time(inp)

## with an output list
check.man.time(rep)

Description

Check rep list

Usage

check.rep(rep, reportmode0 = TRUE)

Arguments

Details

Internal function that checks if rep is fitted spict object.

Value

Nothing

Description

Extract hindcast info from a fitted spict object

Usage

extract.hindcast.info(rep, CI = 0.95, verbose = TRUE)

Arguments

Details

Note that a difference in the timing of the index observations of less than a month are considered acceptable for the estimation of the naive prediction residuals and no warning is printed. If the variability exceeds a month, the predictions are still calculated, but a warning is printed.

Value

A list with hindcast information.

Description

Extracts relevant statistics from the estimation of a simulated data set.

Usage

extract.simstats(rep, inp = NULL, exp = NULL, parnames = NULL)

Arguments

Details

TBA

Value

A list containing the relevant statistics.

Examples

data(pol)
repin <- fit.spict(pol$albacore)
sim <- sim.spict(repin)
rep <- fit.spict(sim)
extract.simstats(rep)

Description

Format date

Usage

fd(d, dec = 2)

Arguments

Value

Correctly formatted date.

Description

Fits aspic to the data contained in the input file

Usage

fit.aspic(
  input,
  do.boot = FALSE,
  nboot = NULL,
  ciperc = NULL,
  verbose = FALSE,
  filebase = "tmp",
  savefile = NULL
)

Arguments

Details

Only works on Linux. This furthermore requires that wine is installed and that aspic7 is installed and available to the PATH.

Value

List containing aspic results.

Description

Fit the Meyer & Millar model using rjags

Usage

fit.jags(
  inp,
  fn,
  n.iter = 10000,
  n.chains = 1,
  burnin = round(n.iter/2),
  thin = 1000
)

Arguments

Value

The raw output of rjags::coda.samples.

Description

Fit the model of Meyer & Millar (1999)

Usage

fit.meyermillar(mminp)

Arguments

Details

Same input structure as for fit.spict(). Fitting the model of Meyer & Millar requires the packages rjags and coda. It furthermore requires that priors are specified for K, r, q, sigma2 (process error variance) and tau2 (observation error variance). Following Meyer & Millar (1999) the priors are:

• "K" log-normal.

• "r" log-normal.

• "q" inverse-gamma.

• "tau2" inverse-gamma.

• "sigma2" inverse-gamma.

See example for how to specify priors.

Value

List containing results

Examples

priors <- list()
priors$K <- c(5.042905, 3.76)
priors$r <- c(-1.38, 3.845)
priors$iq <- c(0.001, 0.0012)
priors$itau2 <- c(1.709, 0.00861342)
priors$isigma2 <- c(3.785518, 0.0102232)
priors$logPini <- -0.223
data(pol)
inp <- pol$albacore
inp$meyermillar$n.iter <- 10000
inp$meyermillar$burnin <- 1000
inp$meyermillar$thin <- 10
inp$meyermillar$n.chains <- 1
inp$meyermillar$priors <- priors
res <- fit.meyermillar(inp)
summary(res$jags)

Description

Fit a continuous-time surplus production model to data.

Usage

fit.spict(inp, verbose = TRUE, dbg = 0)

Arguments

Details

Fits the model using the TMB package and returns a result report containing estimates of model parameters, random effects (biomass and fishing mortality), reference points (Fmsy, Bmsy, MSY) including uncertainties given as standard deviations.

Model parameters using the formulation of Fletcher (1978):

• "logn" Parameter determining the shape of the production curve as in the generalised form of Pella & Tomlinson (1969).

• "logm" Log of maximum sustainable yield.

• "logK" Log of carrying capacity.

• "logq" Log of catchability vector.

• "logsdb" Log of standard deviation of biomass process error.

• "logsdf" Log of standard deviation of fishing mortality process error.

• "logsdi" Log of standard deviation of index observation error.

• "logsdc" Log of standard deviation of catch observation error.

Unobserved states estimated as random effects:

• "logB" Log of the biomass process given by the stochastic differential equation: dB_t = r*B_t*(1-(B_t/K)^n)*dt + sdb*dW_t, where dW_t is Brownian motion.

• "logF" Log of the fishing mortality process given by: dlog(F_t) = f(t, sdf), where the function f depends on the choice of seasonal model.

Other parameters (which are only needed in certain cases):

• "logphi" Log of parameters used to specify the cyclic B spline representing seasonal variation. Used when inp$nseasons > 1 and inp$seasontype = 1.

• "logU" Log of the state of the coupled SDE system used to represent seasonal variation, i.e. when inp$nseasons > 1 and inp$seasontype = 2.

• "loglambda" Log of damping parameter when using the coupled SDE system to represent seasonal variation, i.e. when inp$nseasons > 1 and inp$seasontype = 2.

• "logsdu" Log of standard deviation of process error of U_t (the state of the coupled SDE system) used to represent seasonal variation, i.e. when inp$nseasons > 1 and inp$seasontype = 2.

• "logsde" Log of standard deviation of observation error of effort data. Only used if effort data is part of input.

• "logp1robfac" Log plus one of the coefficient to the standard deviation of the observation error when using a mixture distribution robust toward outliers, i.e. when either inp$robflag = 1 and/or inp$robflagi = 1.

• "logitpp"Logit of the proportion of narrow distribution when using a mixture distribution robust toward outliers, i.e. when either inp$robflag = 1 and/or inp$robflagi = 1.

Parameters that can be derived from model parameters:

• "logr" Log of intrinsic growth rate (r = 4m/K).

• "logalpha" Proportionality factor for the observation noise of the indices and the biomass process noise: sdi = exp(logalpha)*sdb. (normally set to logalpha=0)

• "logbeta" Proportionality factor for the observation noise of the catches and the fishing mortality process noise: sdc = exp(logbeta)*sdf. (this is often difficult to estimate and can result in divergence of the optimisation. Normally set to logbeta=0)

• "logBmsy" Log of the equilibrium biomass (Bmsy) when fished at Fmsy.

• "logFmsy" Log of the fishing mortality (Fmsy) leading to the maximum sustainable yield.

• "MSY" The yield when the biomass is at Bmsy and the fishing mortality is at Fmsy, i.e. the maximum sustainable yield.

The above parameter values can be extracted from the fit.spict() results using get.par().

Model assumptions

1. The intrinsic growth rate (r) represents a combination of natural mortality, growth, and recruitment.

2. The biomass B_t refers to the exploitable part of the stock. Estimates in absolute numbers (K, Bmsy, etc.) should be interpreted in light of this.

3. The stock is closed to migration.

4. Age and size-distribution are stable in time.

5. Constant catchability of the gear used to gather information for the biomass index.

Value

A result report containing estimates of model parameters, random effects (biomass and fishing mortality), reference points (Fmsy, Bmsy, MSY) including uncertainties given as standard deviations.

Examples

data(pol)
rep <- fit.spict(pol$albacore)
Bmsy <- get.par('logBmsy', rep, exp=TRUE)
summary(rep)
plot(rep)

Description

Calculate AIC from a rep list.

Usage

get.AIC(rep)

Arguments

Value

AIC

Description

Find observations of catch and index that overlap

Usage

get.catchindexoverlap(inp)

Arguments

Value

List containing overlapping catch (y) and index (z) observations and their time vectors.

Description

Get column names for data.frames.

Usage

get.colnms()

Value

Vector containing column names of data frames.

Description

Get covariance matrix of two reported quantities not of fixed model parameters. Covariance of fixed model parameters can be found in rep$cov.fixed.

Usage

get.cov(rep, parname1, parname2, cor = FALSE)

Arguments

Value

Covariance matrix of specified parameters.

Description

Calculate E(Binfinity) the fished equilibrium.

Usage

get.EBinf(rep)

Arguments

Details

If a seasonal pattern in F is imposed the annual average F is used for calculating the expectation.

Value

E(Binf).

Description

Estimate fishing mortality factor minimising probability of specified model variable hitting a specified reference level under given fishing mortality

Usage

get.ffac(
  rep,
  var = "logBpBmsy",
  ref = 1,
  problevel = 0.95,
  reportmode = 1,
  verbose = TRUE
)

Arguments

Value

Optimised Fishing mortality for P(Bp<Blim)

Description

Estimate catch for management period based on last catch observations

Usage

get.manC(rep, inp)

Arguments

Details

Internal function that estimates the catch in the management period based on the catch observations in the last year. Only catch observations in the last year are considered. If the management period is longer than a year the catches of the last year are raised. If the management period is shorter than a year, but only annual catches are available, the respective fraction of the last annual catch observation is used. If both the management period and the catch observations are subannual, the subannual catches of the respective 'seasons' of the year corresponding to the 'season' of the management period are used. Be aware that the estimated catch might correspond to a different season(s) than the management period both are subannual and some catch observations are missing.

Value

Table with catches for management period based on last observed catches and corresponding times.

Description

Get limts of any parameter considering all spict objects in rep$man

Usage

get.manlimits(rep, par, CI = 0.95)

Arguments

Value

plotting limits for all reps in rep$man

Description

Get spict object in rep$man with longest time series

Usage

get.manmax(rep)

Arguments

Value

rep in rep$man which hast the longest time series

Description

Get mfrow from the number of plots to be plotted

Usage

get.mfrow(n)

Arguments

Value

Nothing

Description

If multiple growth rates (r) are used (e.g. for a seasonal model), return specified reference point for all instances of r.

Usage

get.msyvec(inp, msy)

Arguments

Value

A list containing reference point estimates with upper and lower CI bounds.

Description

Get number of active priors

Usage

get.no.active.priors(inp)

Arguments

Value

number of active priors

Description

Get order of printed quantities.

Usage

get.order()

Value

Vector containing indices of printed quantities.

Description

Check whether ACF of catch and index residuals is significant in any lags.

Usage

get.osar.pvals(rep)

Arguments

Value

Vector of p-values of length equal to the number of data series.

Description

Extract parameters from a result report as generated by fit.spict.

Usage

get.par(
  parname,
  rep = rep,
  exp = FALSE,
  random = FALSE,
  fixed = FALSE,
  CI = 0.95
)

list.quantities(rep)

Arguments

Details

get.par is a helper function for extracting the value and uncertainty of a specific model parameter, random effect or derived quantity. list.quantities gives the names of all quantities.

Value

get.par returns a matrix with four columns containing respectively: 1) the lower 95% confidence limit; 2) the parameter estimate; 3) the upper 95% confidence limit; 4) the parameter standard deviation in the domain it was estimated (log or non-log). 'list.quantities' returns a vector with the names of all estimated parameters and derived quantities.

Examples

## Run the South Atlantic albacore assessment
data(pol)
rep <- fit.spict(pol$albacore)

## See all quantitites that can be extracted
list.quantities(rep)

## Extract the Bmsy reference point
Bmsy <- get.par('logBmsy', rep, exp=TRUE)

## Extract the exploitable biomass estimates
Best <- get.par('logB', rep, exp=TRUE)

## Extract the estimated carrying capacity
K <- get.par('logK', rep, exp=TRUE)

Description

Get the values of the seasonal spline for F.

Usage

get.spline(logphi, order, dtfine = 1/100)

Arguments

Value

Spline values at the points between 0 and 1 with dtfine as time step.

Description

Get version of spict including git sha1 version if available.

Usage

get.version(pkg = "spict")

Arguments

Value

Package version

Description

Use a simple linear regression to guess m (MSY).

Usage

guess.m(inp, all.return = FALSE)

Arguments

Details

Equations 9.1.7 and 9.1.8 on page 284 of FAO's tropical assessment book are used to guess MSY.

Value

The guess on MSY.

Description

Conduct hindcasting analysis

Usage

hindcast(
  rep,
  npeels = 7,
  reduce.output.size = TRUE,
  mc.cores = 1,
  peel.dtc = FALSE
)

Arguments

Details

This method creates a number of subsets (or peels) specified with the argument npeels by sequentially omitting all observations from the most recent time step (by default corresponding to a year). Then, spict is fitted to each subset while excluding all index observations in the most recent year with data. The resulting fitted spict objects are attached to the base spict object as a list element labeled hindcast and resulted by this method.

The prediction of the excluded index observations relative to the actual excluded index can then be compared to a 'naive' prediction using the preceding index observation. This allows to estimate the Mean Absolute Scaled Error (MASE) for each index calc.mase.

The predicted indices can be visualised with the function plotspict.hindcast.

Value

A spictcls list with the added element hindcast containing the results of the hindcasting analysis. Use plotspict.hindcast to plot these results.

References

Carvalho, F., Winker, H., Courtney, D., Kapur, M., Kell, L., Cardinale, M., Schirripa, M., Kitakado, T., Yemane, D., Piner, K.R. Maunder, M.N., Taylor, I., Wetzel, C.R., Doering, K., Johnsonm, K.F., Methot, R. D. (2021). A cookbook for using model diagnostics in integrated stock assessments. Fisheries Research, 240, 105959.

Kell, L. T., Kimoto, A., & Kitakado, T. (2016). Evaluation of the prediction skill of stock assessment using hindcasting. Fisheries research, 183, 119-127.

Kell, L. T., Sharma, R., Kitakado, T., Winker, H., Mosqueira, I., Cardinale, M., & Fu, D. (2021). Validation of stock assessment methods: is it me or my model talking?. ICES Journal of Marine Science, 78(6), 2244-2255.

Winker, H., Carvalho, F., & Kapur, M. (2018). JABBA: just another Bayesian biomass assessment. Fisheries Research, 204, 275-288.

Examples

data(pol)
inp <- pol$albacore
rep <- fit.spict(inp)
rep <- hindcast(rep, npeels = 5)
plotspict.hindcast(rep)

Description

Inverse logit transform.

Usage

invlogit(a)

Arguments

Value

Inverse logit.

Description

Inverse log "plus one" transform

Usage

invlogp1(a)

Arguments

Details

If a = log(b-1), then the inverse transform is b = 1 + exp(a). Useful for values with lower bound at 1.

Value

Inverse logp1.

Description

Generate latex code for including a figure.

Usage

latex.figure(figfile, reportfile, caption = "")

Arguments

Value

Nothing.

Description

Create profile likelihood

Usage

likprof.spict(input, verbose = FALSE, mc.cores = 1)

Arguments

Details

The "likprof" list must containg the following keys:

• "pars" A character vector of length equal 1 or 2 containing the name(s) of the parameters to calculate the profile likelihood for.

• "parrange" A vector containing the parameter range(s) to profile over: parrange = c(min(par1), max(par1), min(par2), max(par2)).

Optional:

• "nogridpoints" Number of grid points to evaluate the profile likelihood for each parameter. Default: 9. Note: with two parameters the calculation time increases quadratically when increasing the number of gridpoints.

Value

The output is the input with the likelihood profile information added to the likprof key of either inp or rep$inp.

Examples

data(pol)
inp <- pol$albacore
inp$likprof <- list()
inp$likprof$pars <- 'logK'
inp$likprof$parrange <- c(log(80), log(400))
inp$likprof$nogridpoints <- 15
rep <- fit.spict(inp)
rep <- likprof.spict(rep)
plotspict.likprof(rep, logpar=TRUE)

Description

List parameters to which priors can be added

Usage

list.possible.priors()

Value

Prints parameters to which priors can be added.

Description

Create data list used as input to TMB::MakeADFun.

Usage

make.datin(inp, dbg = 0)

Arguments

Value

List to be used as data input to TMB::MakeADFun.

Description

Make fcon vector

Usage

make.fconvec(inp, fcon)

Arguments

Value

Input list containing fconvec

Description

Make ffac vector

Usage

make.ffacvec(inp, ffac)

Arguments

Value

Input list containing ffacvec

Description

Create TMB obj using TMB::MakeADFun and squelch screen printing.

Usage

make.obj(datin, pl, inp, phase = 1)

Arguments

Value

List to be used as data input to TMB.

Description

Creates a pdf file containing the summary output and result plots

Usage

make.report(
  rep,
  reporttitle = "",
  reportfile = "report.tex",
  summaryoutfile = "summaryout.txt",
  keep.figurefiles = FALSE,
  keep.txtfiles = FALSE,
  keep.texfiles = FALSE
)

Arguments

Details

This function probably requires that you are running linux and that you have latex functions installed (pdflatex).

Value

Nothing.

Description

Calculate confidence ellipsis for reference points.

Usage

make.rpellipse(rep)

Arguments

Details

Calculates the confidence ellipsis of logBmsy and logFmsy (last if multiple)

Value

A matrix with two columns containing the x and y coordinates of the ellipsis.

Description

Make a spline design matrix

Usage

make.splinemat(nseasons, order, dtfine = 1/100)

Arguments

Value

Spline design matrix.

Description

Load color of management scenarios.

Usage

man.cols()

Value

Color vector

Description

Select management scenarios

Usage

man.select(rep, scenarios = "all", spictcls = FALSE, verbose = TRUE)

Arguments

Value

A fitted spict object wit selected management scenarios in preferred order in rep$man. This function can also be used to select a specific scenarios in rep$man as the new main spictcls object. By setting the argument spictcls to TRUE, management related catch observations are removed and the retaped spict object of class 'spictcls' is returned, comparable to the object returned by fit.spict. This only works if one scenario is selected (length(scnearios) == 1).

Examples

data(pol)
rep <- fit.spict(pol$albacore)
rep <- manage(rep, c(2,4,6))

## based on names
names(rep$man)
rep1 <- man.select(rep, c("currentF","noF"))

## based on indices
length(rep$man)
rep2 <- man.select(rep, c(1,3))

## select specific scenario as new spictcls object
rep3 <- man.select(rep, 1, spictcls = TRUE)

Description

Get the TAC for the management scenarios

Usage

man.tac(rep, fractileCatch = 0.5, exp = TRUE, verbose = TRUE)

Arguments

Value

rep wit selected management scenarios

Examples

data(pol)
rep <- fit.spict(pol$albacore)
rep <- manage(rep, c(3,4,5))

## Median of predicted catch distributions
man.tac(rep)

## 30th percentile of catch distributions
man.tac(rep, fractileCatch = 0.3)

Description

Print a schematic to the console visualising the management timeline

Usage

man.timeline(x, verbose = TRUE, obsonly = FALSE)

Arguments

Value

Nothing

Examples

data(pol)
inp <- check.inp(pol$albacore)
inp$maninterval <- c(1991,1992)
rep <- fit.spict(inp)

## based on an input list
man.timeline(inp)

## based on an output list
man.timeline(rep)

Description

Calculate predictions under 8 default management scenarios

Usage

manage(
  rep,
  scenarios = "all",
  maninterval = NULL,
  maneval = NULL,
  intermediatePeriodCatch = NULL,
  intermediatePeriodCatchSDFac = 1,
  intermediatePeriodCatchList = NULL,
  verbose = TRUE,
  dbg = 0
)

Arguments

Details

The 8 default scenarios are:

• "1""currentCatch": Keep the catch of the current year (i.e. the last observed catch).

• "2""currentF": Keep the F of the current year.

• "3""Fmsy": Fish at Fmsy i.e. F=Fmsy.

• "4""noF": No fishing, reduce to 1% of current F.

• "5""reduceF25": Reduce F by X%. Default X = 25.

• "6""increaseF25": Increase F by X%. Default X = 25.

• "7""msyHockeyStick": Use ICES MSY hockey-stick advice rule (ICES, 2017).

• "8""ices": Use ICES MSY 35th hockey-stick advice rule (ICES, 2019).

Scenario 7 implements the ICES MSY advice rule for stocks that are assessed using spict (ICES 2017). MSY B_trigger is set equal to B_MSY / 2. Then fishing mortality in the short forecast is calculated as:

F(y+1) = F(y) * min 1, median[B(y+1) / MSY B_trigger] / median[F(y)/F_MSY

Scenario 8 is similar to scenario 7, but includes assessment uncertainty in the predictions by using the 35th percentile of the distributions of the predicted catch, B/B_MSY and F/F_MSY.

Dependent on the start of the management period (e.g. advice year), there might be a time lag between the last observation and the start of the management period, often referred to as the intermediate period. If this is the case, an assumption about the catch during intermediate time period (e.g. assessment year) has to be made. Two meaningful assumptions are:

• 1:The catch in the intermediate period is based on the fishing mortality which is extrapolated from the previous year. This is the default assumption;

• 2:The catch in the intermediate period is directly specified. This could for example be the TAC recommended in the previous year. The catch can be specified by means of the argument intermediatePeriodCatch. Be aware that this catch might correspond to several years or a fraction of a year depending on the time between the last observation and the start of the management period. The function man.timeline can help visualising the default or specified intermediate period in your data. The argument intermediatePeriodCatchSDFac allows to specify the factor with which to multiply the standard deviation of the catch ($\sigma_C$) with. It is thus a measure of the certainty around the catch in the intermediate period. The argument intermediatePeriodCatchList allows to define a list with catches and their intervals. It is a list with the elements 'obsC', 'timeC', 'dtc' and the optional element 'stdevfacC' (which is equal to 1 if not provided).

Value

List containing results of management calculations.

References

ICES. 2017. Report of the Workshop on the Development of the ICES approach to providing MSY advice for category 3 and 4 stocks (WKMSYCat34), 6-10 March 2017, Copenhagen, Denmark. ICES CM 2017/ACOM:47. 53 pp.

Examples

data(pol)
rep <- fit.spict(pol$albacore)
repman <- manage(rep, c(2,4,8))
sumspict.manage(repman) # To print projections

Description

Convert mean and variance to shape and rate of gamma distribution

Usage

meanvar2shaperate(mean, var)

Arguments

Value

Vector containing shape and rate parameters.

Description

Convert mode and fractile to shape and rate in Gamma distribution (only for mode>0, i.e. shape >= 1)

Usage

modefrac2shaperate(mode, xf, f = 0.9)

Arguments

Value

Vector containing shape and rate parameters.

Description

Calculate Mohn's rho for different estimates

Usage

mohns_rho(rep, what = c("FFmsy", "BBmsy"), annualfunc = mean)

Arguments

Details

A function that calculates Mohn's rho for selected estimated quantities. The function allows the user to define the method of aggrgating from the subannual time steps (1/dteuler) into annual values; the default is to take the mean.

Value

A named vector with the Monh's rho value for each quantity.

Examples

data(pol)
inp <- pol$albacore
rep <- fit.spict(inp)
rep <- retro(rep, nretroyear = 4)
mohns_rho(rep)

Description

Plot osar acf

Usage

osar.acf.plot(res, lag.max, pval, ylab)

Arguments

Value

Nothing.

Description

Plot osar qq

Usage

osar.qq.plot(res, pval)

Arguments

Value

Nothing.

Description

Plot model points colored depending on the quarter to which they belong.

Usage

## S3 method for class 'col'
plot(
  time,
  obs,
  obsx = NULL,
  pch = 1,
  add = FALSE,
  typ = "p",
  do.line = TRUE,
  add.legend = FALSE,
  add.vline.at = NULL,
  ...
)

Arguments

Value

Nothing.

Description

Plot summarising spict results.

Usage

## S3 method for class 'spictcls'
plot(x, stamp = get.version(), verbose = TRUE, CI = 0.95, ...)

Arguments

Details

Create a plot containing the following:

• 1. Estimated biomass using plotspict.biomass().

• 2. Estimated fishing mortality using plotspict.f().

• 3. Observed versus predicted catches using plotspict.catch().

• 4. Estimated biomass relative to Bmsy using plotspict.bbmsy().

• 5. Estimated fishing mortality relative to Fmsy using plotspict.ffmsy().

• 6. Estimated F versus estimated B using plotspict.fb().

• 7. Observed versus theoretical production using plotspict.production().

Optional plots included if relevant:

• Estimated seasonal spline using plotspict.season().

• Calculated time-constant using plotspict.tc().

• First prior and corresponding posterior distribution using plotspict.priors().

• One-step-ahead residuals of catches using plotspict.osar().

• One-step-ahead residuals of catches using plotspict.osar().

If no management scenarios are included in rep$man, the grey vertical line corresponds to the time of the last observation. If management scenarios are included in rep$man, the prediction and confidence intervals of the base scenario (rep) are omitted and instead the projections of the different management scenarios are drawn in different colours. Dotted lines of the management scenarios reflect the intermediate period, while solid lines reflect the management period. Additionally, two vertical lines correspond to the start and end of the management period.

Be aware that potential catch intervals of more a year, e.g. biennial assessment so that the intermediate period spans two years, or management period spans two years, are equally split up into annual intervals.

Be aware of the fact that the catches represent intervals, where the length of the interval is indicated by dtc, e.g. with $dtc = 1, C(1990) = [1990,1990[$. In the plot the catches (and vertical lines) correspond to the beginning of the catch interval. It might thus seem as if the time of the vertical lines and the management interval would not align.

Value

Nothing.

Examples

data(pol)
rep <- fit.spict(pol$albacore)
plot(rep)

Description

Plot summarising spict results (alternative plot composition)

Usage

plot2(rep, stamp = get.version(), verbose = TRUE, CI = 0.95, ...)

Arguments

Details

Create a plot containing the following:

• 1. Estimated biomass relative to Bmsy using plotspict.bbmsy().

• 2. Estimated fishing mortality relative to Fmsy using plotspict.ffmsy().

• 3. Observed versus predicted catches using plotspict.catch().

• 4. Estimated F versus estimated B using plotspict.fb().

If no management scenarios are included in rep$man, the grey vertical line corresponds to the time of the last observation. If management scenarios are included in rep$man, the prediction and confidence intervals of the base scenario (rep) are omitted and instead the projections of the different management scenarios are drawn in different colours. Dotted lines of the management scenarios reflect the intermediate period, while solid lines reflect the management period. Additionally, two vertical lines correspond to the start and end of the management period.

Be aware that potential catch intervals of more a year, e.g. biennial assessment so that the intermediate period spans two years, or management period spans two years, are equally split up into annual intervals.

Be aware of the fact that the catches represent intervals, where the length of the interval is indicated by dtc, e.g. with $dtc = 1, C(1990) = [1990,1990[$. In the plot the catches (and vertical lines) correspond to the beginning of the catch interval. It might thus seem as if the time of the vertical lines and the management interval would not align.

Value

Nothing.

Examples

data(pol)
rep <- fit.spict(pol$albacore)
plot2(rep)

Description

Plot priors of Meyer & Millar model

Usage

## S3 method for class 'priors'
plot(nm, priorsin, add = TRUE, ...)

Arguments

Value

Nothing.

Description

Plot estimated B/Bmsy.

Usage

plotspict.bbmsy(
  rep,
  logax = FALSE,
  main = "Relative biomass",
  ylim = NULL,
  plot.obs = TRUE,
  qlegend = TRUE,
  lineat = 1,
  xlab = "Time",
  stamp = get.version(),
  verbose = TRUE,
  CI = 0.95
)

Arguments

Details

Plots estimated B/Bmsy.

If no management scenarios are included in rep$man, the grey vertical line corresponds to the time of the last observation. If management scenarios are included in rep$man, the prediction and confidence intervals of the base scenario (rep) are omitted and instead the projections of the different management scenarios are drawn in different colours. Dotted lines of the management scenarios reflect the intermediate period, while solid lines reflect the management period. Additionally, two vertical lines correspond to the start and end of the management period.

Value

Nothing.

Examples

data(pol)
rep <- fit.spict(pol$albacore)
plotspict.bbmsy(rep)

Description

Plot estimated biomass.

Usage

plotspict.biomass(
  rep,
  logax = FALSE,
  main = "Absolute biomass",
  ylim = NULL,
  plot.obs = TRUE,
  qlegend = TRUE,
  xlab = "Time",
  ylab = NULL,
  rel.axes = TRUE,
  rel.ci = TRUE,
  stamp = get.version(),
  verbose = TRUE,
  CI = 0.95
)

Arguments

Details

Plots estimated biomass, Bmsy with confidence limits.

If no management scenarios are included in rep$man, the grey vertical line corresponds to the time of the last observation. If management scenarios are included in rep$man, the prediction and confidence intervals of the base scenario (rep) are omitted and instead the projections of the different management scenarios are drawn in different colours. Dotted lines of the management scenarios reflect the intermediate period, while solid lines reflect the management period. Additionally, two vertical lines correspond to the start and end of the management period.

Value

Nothing.

Examples

data(pol)
rep <- fit.spict(pol$albacore)
plotspict.biomass(rep)

Description

Plot the expected biomass trend

Usage

plotspict.btrend(rep, CI = 0.95)

Arguments

Value

Nothing.

Description

Plot observed catch and predictions.

Usage

plotspict.catch(
  rep,
  main = "Catch",
  ylim = NULL,
  qlegend = TRUE,
  lcol = "blue",
  xlab = "Time",
  ylab = NULL,
  stamp = get.version(),
  verbose = TRUE,
  CI = 0.95
)

Arguments

Details

Plots observed catch and predictions using the current F and Fmsy. The plot also contains the equilibrium catch if the current F is maintained. If no management scenarios are included in rep$man, the grey vertical line corresponds to the time of the last observation.

If management scenarios are included in rep$man, the prediction and confidence intervals of the base scenario (rep) are omitted and instead the projections of the different management scenarios are drawn in different colours. Generally, dotted lines of the management scenarios reflect the intermediate period, while solid lines reflect the management period. The catch of management period which are longer than 1 year are split up equally into annual intervals. Two vertical lines correspond to the start and end of the management period, respectively. However, there are special cases in which there is only one or no vertical line drawn, the catch trajectories are missing completely, or the line of the catch trajectory is solid even in the intermediate period. These cases and their implications on the annual catch plot are described in the following:

• If the management period is shorter than a year, no catch trajectories are drawn and there is only one vertical line indicating the start of the assessment period.

• If the management timeline differs between the scenarios in rep$man, no vertical lines are drawn as they would be at different times for each scenario.

• If the management period cannot be split equally into annual intervals, e.g. because it is 1.5 years long, the uneven remains are not displayed, in this example only the catch representative of one year is displayed. Additionally, the second vertical line indicating the end of the management period is omitted.

• If the intermediate period is shorter or longer than a year, e.g. 0.5 or 1.25 years, the lines of the management period start at the time of the last observation, because the catch in the intermediate period cannot be aggregated and displayed correctly. Additionally, the first vertical line indicating the start of the management period is omitted.

All catches in SPiCT represent intervals, where the length of the interval is indicated by dtc, e.g. with $dtc = 1, C(1990) = [1990,1990[$. In the plot the catches (and vertical lines) correspond to the beginning of the catch interval. It might thus seem as if the time of the vertical lines and the management interval would not align.

Value

Nothing.

Examples

data(pol)
rep <- fit.spict(pol$albacore)
plotspict.catch(rep)

Description

Plot catch and index data.

Usage

plotspict.ci(inp, stamp = get.version())

Arguments

Value

Nothing

Description

Compare different spict fits

Usage

plotspict.compare(
  rep,
  ...,
  varname = c("B", "F", "C", "BBmsy", "FFmsy", "P"),
  exp = TRUE,
  CI = 0.95,
  plot.unc = TRUE,
  col = c("dodgerblue2", "darkorange1", "forestgreen", "goldenrod1", "purple2",
    "firebrick3", "skyblue4", "darkgreen", "salmon3", "brown2"),
  asp = 2,
  plot.legend = TRUE,
  stamp = get.version()
)

Arguments

Details

This function plots the results of spict fits in a single plot

Examples

data(pol)
inp <- pol$albacore
inp$dteuler <- 1/8
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))

Description

Compare one variable of different spict fits

Usage

plotspict.compare.one(
  rep,
  ...,
  varname = c("B", "F", "C", "BBmsy", "FFmsy", "P"),
  exp = TRUE,
  CI = 0.95,
  plot.unc = 1,
  col = c("dodgerblue2", "darkorange1", "forestgreen", "goldenrod1", "purple2",
    "firebrick3", "skyblue4", "darkgreen", "salmon3", "brown2"),
  xlab = NULL,
  ylab = NULL,
  main = NULL,
  plot.legend = TRUE,
  stamp = get.version()
)

Arguments

Details

This function plots the results of the hindcasting cross validation analysis for each index.

Examples

data(pol)
inp <- pol$albacore
inp$dteuler <- 1/8
rep1 <- fit.spict(inp)
inp$priors$logbkfrac <- c(log(0.3), 0.1, 1)
rep2 <- fit.spict(inp)
## Not run: 
plotspict.compare.one(rep1, rep2, varname = "F")

## End(Not run)

Description

Plot input data

Usage

plotspict.data(
  inpin,
  MSY = NULL,
  one.index = NULL,
  qlegend = TRUE,
  stamp = get.version()
)

Arguments

Value

Nothing

Description

Plot model diagnostic (data, residuals, and more)

Usage

plotspict.diagnostic(
  rep,
  lag.max = 4,
  qlegend = TRUE,
  plot.data = TRUE,
  mfcol = FALSE,
  stamp = get.version()
)

Arguments

Value

Nothing.

Examples

data(pol)
rep <- fit.spict(pol$albacore)
rep <- calc.osa.resid(rep)
plotspict.diagnostic(rep)

Description

Plot model diagnostics regarding processes and process residuals

Usage

plotspict.diagnostic.process(
  rep,
  lag.max = 4,
  qlegend = TRUE,
  plot.data = TRUE,
  mfcol = FALSE,
  add.loess = FALSE,
  span = 0.75,
  stamp = get.version()
)

Arguments

Value

Invisible NULL

See Also

process.resid

Examples

data(pol)
inp <- pol$albacore
rep <- fit.spict(inp)
rep <- calc.process.resid(rep)
plotspict.diagnostic.process(rep)

Description

Plot estimated fishing mortality.

Usage

plotspict.f(
  rep,
  logax = FALSE,
  main = "Absolute fishing mortality",
  ylim = NULL,
  plot.obs = TRUE,
  qlegend = TRUE,
  xlab = "Time",
  ylab = NULL,
  rel.axes = TRUE,
  rel.ci = TRUE,
  stamp = get.version(),
  verbose = TRUE,
  CI = 0.95
)

Arguments

Details

Plots estimated fishing mortality with Fmsy and associated confidence interval.

If no management scenarios are included in rep$man, the grey vertical line corresponds to the time of the last observation. If management scenarios are included in rep$man, the prediction and confidence intervals of the base scenario (rep) are omitted and instead the projections of the different management scenarios are drawn in different colours. Dotted lines of the management scenarios reflect the intermediate period, while solid lines reflect the management period. Additionally, two vertical lines correspond to the start and end of the management period.

Value

Nothing.

Examples

data(pol)
rep <- fit.spict(pol$albacore)
plotspict.f(rep)

Description

Plot fishing mortality versus biomass.

Usage

plotspict.fb(
  rep,
  logax = FALSE,
  plot.legend = TRUE,
  man.legend = TRUE,
  ext = TRUE,
  rel.axes = FALSE,
  xlim = NULL,
  ylim = NULL,
  labpos = c(1, 1),
  xlabel = NULL,
  stamp = get.version(),
  verbose = TRUE,
  CI = 0.95
)

Arguments

Details

Plots estimated fishing mortality as a function of biomass together with reference points and the prediction for next year given a constant F. The equilibrium biomass for F fixed to the current value is also plotted.

The predicted trajectory (or trajectories of different management scenarios) are only plotted for annnual data.

Value

Nothing.

Examples

data(pol)
rep <- fit.spict(pol$albacore)
plotspict.fb(rep)

Description

Plot estimated relative fishing mortality.

Usage

plotspict.ffmsy(
  rep,
  logax = FALSE,
  main = "Relative fishing mortality",
  ylim = NULL,
  plot.obs = TRUE,
  qlegend = TRUE,
  lineat = 1,
  xlab = "Time",
  stamp = get.version(),
  verbose = TRUE,
  CI = 0.95
)

Arguments

Details

Plots estimated fishing mortality with Fmsy and associated confidence interval.

If no management scenarios are included in rep$man, the grey vertical line corresponds to the time of the last observation. If management scenarios are included in rep$man, the prediction and confidence intervals of the base scenario (rep) are omitted and instead the projections of the different management scenarios are drawn in different colours. Dotted lines of the management scenarios reflect the intermediate period, while solid lines reflect the management period. Additionally, two vertical lines correspond to the start and end of the management period.

Value

Nothing.

Examples

data(pol)
rep <- fit.spict(pol$albacore)
plotspict.ffmsy(rep)

Description

Plot estimated time-varying growth

Usage

plotspict.growth(
  rep,
  logax = FALSE,
  main = "Time-varying growth",
  ylim = NULL,
  xlim = NULL,
  xlab = "Time",
  plot.ci = TRUE,
  stamp = get.version(),
  CI = 0.95
)

Arguments

Details

Plots estimated time-varying growth

Value

Nothing.

Description

A comparison of all management scenarios in rep$man. For each scenario, the fishing mortality in the management period is shown against the relative biomass at the point of the harvest control rule (HCR) evaluation.

Usage

plotspict.hcr(rep, xlim = c(0, 3), CI = 0.95)

Arguments

Value

Invisible NULL

See Also

manage

Examples

rep <- fit.spict(pol$hake)
rep <- manage(rep)
plotspict.hcr(rep)

Description

Hindcasting plot for indices

Usage

plotspict.hindcast(
  rep,
  add.mase = TRUE,
  CI = 0.95,
  verbose = TRUE,
  xlim = NULL,
  ylim = NULL,
  xlab = "Year",
  ylab = NA,
  plot.log = TRUE,
  mfrow = NULL,
  mar = c(2, 2, 3, 1) + 0.1,
  oma = c(3, 3, 1, 1),
  legend.title = NULL,
  legend.pos = "topright",
  legend.ncol = 1,
  asp = 2,
  stamp = get.version()
)

Arguments

Details

This function plots the results of the hindcasting cross validation analysis for each index.

Value

MASE or Invisible NULL

References

Carvalho, F., Winker, H., Courtney, D., Kapur, M., Kell, L., Cardinale, M., Schirripa, M., Kitakado, T., Yemane, D., Piner, K.R. Maunder, M.N., Taylor, I., Wetzel, C.R., Doering, K., Johnsonm, K.F., Methot, R. D. (2021). A cookbook for using model diagnostics in integrated stock assessments. Fisheries Research, 240, 105959.

Kell, L. T., Kimoto, A., & Kitakado, T. (2016). Evaluation of the prediction skill of stock assessment using hindcasting. Fisheries research, 183, 119-127.

Kell, L. T., Sharma, R., Kitakado, T., Winker, H., Mosqueira, I., Cardinale, M., & Fu, D. (2021). Validation of stock assessment methods: is it me or my model talking?. ICES Journal of Marine Science, 78(6), 2244-2255.

Winker, H., Carvalho, F., & Kapur, M. (2018). JABBA: just another Bayesian biomass assessment. Fisheries Research, 204, 275-288.

See Also

hindcast

Examples

data(pol)
inp <- pol$albacore
rep <- fit.spict(inp)
rep <- hindcast(rep, npeels = 5)
plotspict.hindcast(rep)

Description

Plots influence statistics of observations.

Usage

plotspict.infl(rep, stamp = get.version())

Arguments

Details

TBA

Value

Nothing.

Description

Plots summary of influence statistics of observations.

Usage

plotspict.inflsum(rep, stamp = get.version())

Arguments

Details

TBA

Value

Nothing.

Description

Plots result of likelihood profiling.

Usage

plotspict.likprof(input, logpar = FALSE, stamp = get.version(), CI = 0.95)

Arguments

Details

TBA

Value

Nothing but shows a plot.

Description

Plot one-step-ahead residuals

Usage

plotspict.osar(rep, collapse.I = TRUE, qlegend = TRUE)

Arguments

Details

Plots observed versus predicted catches.

Value

Nothing.

Examples

data(pol)
rep <- fit.spict(pol$albacore)
rep <- calc.osa.resid(rep)
plotspict.osar(rep)

Description

Plot priors and posterior distribution.

Usage

plotspict.priors(
  rep,
  do.plot = NULL,
  stamp = get.version(),
  automfrow = TRUE,
  CI = 0.95
)

Arguments

Value

Nothing

Description

Plot theoretical production curve and estimates.

Usage

plotspict.production(
  rep,
  n.plotyears = 40,
  main = "Production curve",
  stamp = get.version(),
  CI = 0.95
)

Arguments

Details

Plots the theoretical production curve (production as a function of biomass) as calculated from the estimated model parameters. Overlaid is the estimated production/biomass trajectory.

Value

Nothing.

Examples

data(pol)
rep <- fit.spict(pol$albacore)
plotspict.production(rep)

Description

Plot results of retrospective analysis

Usage

plotspict.retro(rep, stamp = get.version(), add.mohn = TRUE, CI = 0.95)

plotspict.retro.fixed(rep, CI = 0.95)

Arguments

Value

Ivisible NULL. If add.mohn is TRUE, plotspict.retro returns the Mohn's rho for B/Bmsy and F/Fmsy.

Note

The retrospective runs that did not converge are excluded from the plots and from the calculation of Mohn's rho. A message is displayed in such a case.

Description

Plot the mean F cycle

Usage

plotspict.season(rep, stamp = get.version(), CI = 0.95)

Arguments

Details

If seasonal data are available the seasonal cycle in the fishing mortality can be estimated. This function plots this mean F cycle.

Value

Nothing.

Description

Plot time constant.

Usage

plotspict.tc(rep, main = "Time to Bmsy", stamp = get.version(), CI = 0.95)

Arguments

Details

Plots the time required for the biomass to reach a certain proportion of Bmsy. The time required to reach 95% of Bmsy is highlighted.

Value

Nothing.

Examples

data(pol)
rep <- fit.spict(pol$albacore)
plotspict.tc(rep)

Description

Fisheries data included in Polacheck et al. (1993).

Usage

data(pol)

Format

Data are lists containing data and initial values for estimation formatted to be used as an input to fit.spict().

Details

Fisheries data for South Atlantic albacore, northern Namibian hake, and New Zealand rock lobster.

Source

Polacheck et al. (1993), Canadian Journal of Fisheries and Aquatic Science, vol 50, pp. 2597-2607.

Examples

data(pol)
rep <- fit.spict(inp=pol$albacore)
rep <- fit.spict(inp=pol$hake)
rep <- fit.spict(inp=pol$lobster)

Description

Helper function for sim.spict().

Usage

## S3 method for class 'b'
predict(B0, F0, gamma, m, K, n, dt, sdb, btype)

Arguments

Value

Predicted biomass at the end of dt.

Description

Helper function for sim.spict().

Usage

## S3 method for class 'logf'
predict(logF0, dt, sdf, efforttype)

Arguments

Value

Predicted F at the end of dt.

Description

Helper function for sim.spict().

Usage

## S3 method for class 'logmre'
predict(logmre0, dt, sdm, psi, logm)

Arguments

Value

Predicted mre at the end of dt.

Description

Output a summary of a fit.spict() run.

Usage

## S3 method for class 'spictcls'
print(x, ...)

Arguments

Value

Nothing.

Description

Estimate deviation between targeted and realised probability of specified model variable hitting a specified reference level under given fishing mortality

Usage

probdev(
  ffac = 1,
  rep,
  var = "logBpBmsy",
  ref = 1,
  problevel = 0.95,
  reportmode = 1,
  getFrac = FALSE,
  verbose = FALSE
)

Arguments

Value

Returns deviation between targeted and realised probability of hitting specified reference levels under given fishing mortality

Description

Prune a fitted spict object to the core elements

Usage

prune.baserun(rep)

Arguments

Value

Fitted spict core object.

Description

Adds the x-axis to influence plots

Usage

put.xax(rep)

Arguments

Details

TBA

Value

Nothing.

Description

Reads ASPIC input file.

Usage

read.aspic(filename)

Arguments

Details

Reads an input file following the ASPIC 7 format described in the ASPIC manual (found here http://www.mhprager.com/aspic.html).

Value

A list of input variables that can be used as input to fit.spict().

Examples

## Not run: 
filename <- 'YFT-SSE.a7inp' # or some other ASPIC 7 input file
inp <- read.aspic(filename)
rep <- fit.spict(inp)
summary(rep)
plot(rep)

## End(Not run)

Description

Reads the parameter estimates of an Aspic result file.

Usage

read.aspic.res(filename)

Arguments

Details

TBA

Value

Vector containing the parameter estimates.

Description

Draw CI around a reference point using polygon

Usage

refpointci(t, ll, ul, cicol = "ivory2")

Arguments

Value

Spline design matrix.

Description

Helper function for calc.osar.resid that calculates residual statistics.

Usage

res.diagn(resid, id, name = "")

Arguments

Value

List containing residual statistics in 'diagn', shapiro output in 'shapiro', and bias output in 'bias'.

Description

Retape a fitted spict model based on an updated input list

Usage

retape.spict(rep, inp, verbose = FALSE, dbg = 0, mancheck = TRUE)

Arguments

Details

This function reevaluates derived variables of a fitted spict model after updating model settings, such as the prediction horizon or fishing mortality during the predicted time period.

Value

Retaped spict model

Description

Conduct retrospective analysis

Usage

retro(rep, nretroyear = 5, reduce_output_size = TRUE, mc.cores = 1)

Arguments

Details

A retrospective analysis consists of estimating the model with later data points removed sequentially one year at a time.

Value

A rep list with the added key retro containing the results of the retrospective analysis. Use plotspict.retro() to plot these results.

Examples

data(pol)
inp <- pol$albacore
rep <- fit.spict(inp)
rep <- retro(rep, nretroyear=6)
plotspict.retro(rep)

Description

Load season colors.

Usage

season.cols(modin)