| Title: | Functions Related to ICES Advice |
|---|---|
| Description: | A collection of functions that facilitate computational steps related to advice for fisheries management, according to ICES guidelines. These include methods for calculating reference points and model diagnostics. |
| Authors: | Arni Magnusson [aut], Colin Millar [aut, cre], Anne Cooper [aut], Ruben Verkempynck [ctb], Raphael Girardin [ctb] |
| Maintainer: | Colin Millar <[email protected]> |
| License: | GPL-3 |
| Version: | 2.1.1 |
| Built: | 2024-11-17 02:53:02 UTC |
| Source: | https://github.com/ices-tools-prod/icesAdvice |
A collection of functions that facilitate computational steps related to advice for fisheries management, according to ICES guidelines. These include methods for calculating reference points and model diagnostics.
Calculate ICES advice:
DLS3.2 |
DLS method 3.2 |
icesRound |
rounding method |
Calculate PA reference points and sigma:
Bpa |
from Blim |
Fpa |
from Flim |
sigmaCI |
from confidence interval |
sigmaPA |
from PA reference points |
Other calculations:
agesFbar |
suitable age range for Fbar |
mohn |
Mohn's rho retrospective diagnosis |
Read and write files:
read.dls |
read DLS3.2 results from file |
write.dls |
write DLS3.2 results to file
|
Example tables:
gss |
Greater silver smelt catch at age |
shake |
Southern hake retro |
Arni Magnusson, Colin Millar, and Anne Cooper.
ICES advice: https://ices.dk/advice
Calculate a suitable age range for average fishing mortality (Fbar).
agesFbar(cn, probs = c(0.1, 0.9), plot = FALSE, main = NULL, xlab = NULL, ylab = NULL, col = NULL)agesFbar(cn, probs = c(0.1, 0.9), plot = FALSE, main = NULL, xlab = NULL, ylab = NULL, col = NULL)
cn |
a vector, matrix, table, or data frame containing the catch at age
as integers or decimal values. Ages are indicated in names if
|
probs |
a vector of probabilities, e.g. |
plot |
whether to indicate the cut points on a bar plot. |
main |
passed to |
xlab |
passed to |
ylab |
passed to |
col |
passed to |
Vector of two values, the lower and upper age. These can be rounded to construct an Fbar age range.
The calculations are more complicated than traditional quantiles, since the ages for Fbar are discrete and not continuous. For example, with a continuous variable, Fbar[5-6] has width one, but for a discrete variable it has width two.
The algorithm cuts off exactly probs[1] and 1-probs[2] from the
left and right tails of the distribution and returns the corresponding values
that can be rounded with a simple round(agesFbar(x)) to get a discrete
age range that covers close to diff(probs) of the catch.
With plot = TRUE, the white bars are outside the Fbar age range and
the dark gray bars are inside it. The bars affected by the lower and upper
cut points are drawn in a shade between white and dark gray, providing a
visual cue whether that age should be in the Fbar age range or not.
The calculated limits are somewhat unintuitive for interpretation, but they
are statistically precise and unbiased. A lower limit of 6.8 means that 80%
of age six was cut off the left tail, so rounding the lower age of Fbar to 7
appropriate. An upper limit of 14.1 means that all of age 14 and 10% of age
15 is inside the range, so rounding the upper age of Fbar to 14 is
appropriate. See plot = TRUE example below.
Arni Magnusson.
gss is an example catch-at-age table.
quantile is the base function to calculate quantiles.
icesAdvice-package gives an overview of the package.
agesFbar(gss) agesFbar(gss, plot=TRUE) agesFbar(gss, plot=TRUE, main="Greater silver smelt in 5b and 6a") # 50% coverage instead of 80% agesFbar(gss, c(0.25, 0.75)) # Include only the last 20 years agesFbar(tail(gss, 20)) # Tests cn1 <- setNames(c(100,400,400,100), 3:6) cn2 <- setNames(c(99,401,401,99), 3:6) cn3 <- setNames(c(101,399,399,101), 3:6) cn4 <- setNames(c(500,400,50,50), 3:6) cn5 <- setNames(c(50,50,400,500), 3:6) agesFbar(cn1) agesFbar(cn2) agesFbar(cn3) agesFbar(cn4) agesFbar(cn5)agesFbar(gss) agesFbar(gss, plot=TRUE) agesFbar(gss, plot=TRUE, main="Greater silver smelt in 5b and 6a") # 50% coverage instead of 80% agesFbar(gss, c(0.25, 0.75)) # Include only the last 20 years agesFbar(tail(gss, 20)) # Tests cn1 <- setNames(c(100,400,400,100), 3:6) cn2 <- setNames(c(99,401,401,99), 3:6) cn3 <- setNames(c(101,399,399,101), 3:6) cn4 <- setNames(c(500,400,50,50), 3:6) cn5 <- setNames(c(50,50,400,500), 3:6) agesFbar(cn1) agesFbar(cn2) agesFbar(cn3) agesFbar(cn4) agesFbar(cn5)
Calculate the value of Bpa from Blim and sigmaB.
Bpa(Blim, sigmaB)Bpa(Blim, sigmaB)
Blim |
the value of the Blim reference point. |
sigmaB |
the estimation uncertainty in B (standard error of logSSB in the terminal year). |
Value of Bpa.
The purpose of PA reference points is to apply a precautionary approach in fisheries management.
By comparing the current B to Bpa, one can answer the question: are we at least 95% sure that B is above Blim, given the estimation uncertainty?
The ICES (2017) technical guidelines define Bpa as:
Arni Magnusson.
ICES (2017). ICES fisheries management reference points for category 1 and 2 stocks. ICES Advice Technical Guidelines 12.4.3.1. doi:10.17895/ices.pub.3036.
Fpa calculates that reference point from Flim and sigmaF.
sigmaPA calculates the implicit sigma from PA reference points.
icesAdvice-package gives an overview of the package.
Bpa(100, 0.15)Bpa(100, 0.15)
Apply ICES method 3.2 to calculate catch advice for data-limited stocks (DLS).
DLS3.2(lastadvice, index, len = c(3, 2), buffer = FALSE, i1, i2)DLS3.2(lastadvice, index, len = c(3, 2), buffer = FALSE, i1, i2)
lastadvice |
last catch advice given for this stock. |
index |
stock size index. |
len |
two integers, indicating the desired lengths of reference vectors. |
buffer |
whether to apply a -20% precautionary buffer. |
i1 |
included for backward compatibility, use |
i2 |
included for backward compatibility, use |
This function compares the average values of two reference vectors i1
and i2. In the simplest case, only lastadvice and index
are required to calculate the advice.
The default value of len = c(3, 2) produces vectors i1 and
i2 of lengths 3 and 2,
i1 = (I[n-4], I[n-3], I[n-2])
i2 = (I[n-1], I[n])
where I is a stock size index of length n.
Other vector lengths can be used, such as len = c(5, 2) to get
i1 = (I[n-6], I[n-5], I[n-4], I[n-3], I[n-2])
i2 = (I[n-1], I[n])
Finally, a -20% precautionary buffer can be applied at the end of all calculations.
See the ICES (2012) guidance report for details.
A list containing the resulting advice and other elements showing
intermediate steps in the calculations.
Anne Cooper and Arni Magnusson.
ICES (2012). ICES DLS guidance report: ICES implementation of advice for data-limited stocks in 2012 in its 2012 advice. ICES CM 2012/ACOM:68. doi:10.17895/ices.pub.5322.
read.dls and write.dls read and write DLS3.2
results to file.
icesAdvice-package gives an overview of the package.
# Three hypothetical surveys survey <- data.frame(year=2001:2010, randu[1:10,]) DLS3.2(1000, survey$x) DLS3.2(1000, survey$y) DLS3.2(1000, survey$y, len=c(5,2)) DLS3.2(1000, survey$z) DLS3.2(1000, survey$z, buffer=TRUE) # Plot output <- DLS3.2(1000, survey$y) plot(y~year, survey, ylab="index", type="b", lty=3) segments(2006, output$i1bar, 2008, lwd=2) segments(2009, output$i2bar, 2010, lwd=2)# Three hypothetical surveys survey <- data.frame(year=2001:2010, randu[1:10,]) DLS3.2(1000, survey$x) DLS3.2(1000, survey$y) DLS3.2(1000, survey$y, len=c(5,2)) DLS3.2(1000, survey$z) DLS3.2(1000, survey$z, buffer=TRUE) # Plot output <- DLS3.2(1000, survey$y) plot(y~year, survey, ylab="index", type="b", lty=3) segments(2006, output$i1bar, 2008, lwd=2) segments(2009, output$i2bar, 2010, lwd=2)
Calculate the value of Fpa from Flim and sigmaF.
Fpa(Flim, sigmaF)Fpa(Flim, sigmaF)
Flim |
the value of the Flim reference point. |
sigmaF |
the estimation uncertainty in F (standard error of logF in the terminal year). |
Value of Fpa.
The purpose of PA reference points is to apply a precautionary approach in fisheries management.
By comparing the current F to Fpa, one can answer the question: are we at least 95% sure that F is below Flim, given the estimation uncertainty?
The ICES (2017) technical guidelines define Fpa as:
The Fpa function can also be used to calculate reference points based
on harvest rate: Hpa from Hlim and sigmaH.
Arni Magnusson.
ICES (2017). ICES fisheries management reference points for category 1 and 2 stocks. ICES Advice Technical Guidelines 12.4.3.1. doi:10.17895/ices.pub.3036.
Bpa calculates that reference point from Blim and sigmaB.
sigmaPA calculates the implicit sigma from PA reference points.
icesAdvice-package gives an overview of the package.
Fpa(0.90, 0.15)Fpa(0.90, 0.15)
Catch at age of greater silver smelt in areas 5b and 6a, in crosstab format.
gssgss
Data frame containing 18 columns:
Year |
year |
5 |
number of five-year-olds in the catch (thousands) |
6 |
number of six-year-olds in the catch (thousands) |
...
|
|
21 |
number of twenty-one-year-olds in the catch (thousands) |
The data were presented and analyzed at the 2020 ICES Benchmark Workshop on Greater Silver Smelt.
agesFbar calculates a suitable age range for Fbar.
icesAdvice-package gives an overview of the package.
gss agesFbar(gss)gss agesFbar(gss)
Round values according to the ICES Advice Technical Guidelines.
icesRound(x, percent = FALSE, sign = percent, na = "")icesRound(x, percent = FALSE, sign = percent, na = "")
x |
the values to round. |
percent |
whether to format values with a percent suffix. |
sign |
whether to format values with a sign prefix. |
na |
what to return when x is |
Rounded values as a noquote string vector, retaining trailing zeros.
This function implements the following ICES rounding method:
| i) | Round to two significant figures when the first non-zero digit is 2 or larger. |
| ii) | Round to three significant figures when the first non-zero digit is 1. |
As indicated in the ICES (2017) technical guidelines, this rounding method
should not be applied to biomass, catch, or number of individuals. For those
quantities, use the normal round function instead.
Colin Millar and Arni Magnusson.
ICES (2017). Rounding rules to be applied in ICES advice. ICES Advice Technical Guidelines 16.5.4. doi:10.17895/ices.pub.3038.
signif rounds values to a specified number of significant
digits.
icesAdvice-package gives an overview of the package.
icesRound(0.123456) icesRound(0.2468) ## Formatted string or numeric icesRound(1.0) as.numeric(icesRound(1.0)) ## Percent, sign, NA icesRound(33.33, percent = TRUE) icesRound(33.33, sign = TRUE) icesRound(c(1, NA, 3)) icesRound(c(1, NA, 3), na = NA) ## Example from the ICES Technical Guidelines Actual <- c(0.35776, 0.34665, 0.202, 0.12665, 0.001567, 0.002567, 0.013415, 0.02315, 1.168, 2.15678) Rounded <- icesRound(Actual) print(data.frame(Actual = as.character(Actual), Rounded), row.names = FALSE) ## Continued example from Guidelines, now rounding percentages Actual <- c(9.546, 10.546, 23.445, -1.482, -9.09, 0.51, 130.11, 584) Rounded <- icesRound(Actual, percent = TRUE) print(data.frame(Actual = as.character(Actual), Rounded), row.names = FALSE)icesRound(0.123456) icesRound(0.2468) ## Formatted string or numeric icesRound(1.0) as.numeric(icesRound(1.0)) ## Percent, sign, NA icesRound(33.33, percent = TRUE) icesRound(33.33, sign = TRUE) icesRound(c(1, NA, 3)) icesRound(c(1, NA, 3), na = NA) ## Example from the ICES Technical Guidelines Actual <- c(0.35776, 0.34665, 0.202, 0.12665, 0.001567, 0.002567, 0.013415, 0.02315, 1.168, 2.15678) Rounded <- icesRound(Actual) print(data.frame(Actual = as.character(Actual), Rounded), row.names = FALSE) ## Continued example from Guidelines, now rounding percentages Actual <- c(9.546, 10.546, 23.445, -1.482, -9.09, 0.51, 130.11, 584) Rounded <- icesRound(Actual, percent = TRUE) print(data.frame(Actual = as.character(Actual), Rounded), row.names = FALSE)
Calculate Mohn's rho, the average relative bias of retrospective estimates.
mohn(x, peels = 5, details = FALSE, plot = FALSE, ...)mohn(x, peels = 5, details = FALSE, plot = FALSE, ...)
x |
a matrix or data frame containing retrospective estimates in columns, with years as row names. |
peels |
the number of retrospective peels to use in the calculation of
rho, or |
details |
whether to return the intermediate calculations of relative bias. |
plot |
whether to plot the retrospective trajectories. |
... |
passed to |
The default value peels = 5 is based on the ICES (2018) guidelines.
The basic plot = TRUE functionality is intended to quickly visualize
the calculation of Mohn's rho. To produce a fully formatted plot, bypass the
mohn function and plot the x data directly.
Mohn's rho, along with intermediate calculations if details = TRUE.
Relative bias is defined as
and Mohn's rho is the average relative bias:
See Mohn (1999), Brooks and Legault (2016), ICES (2018), and
mohn(shake, details=TRUE) for details.
Arni Magnusson, with a contribution from Ruben Verkempynck.
Brooks, E.N. and Legault, C.M. (2016). Retrospective forecasting — evaluating performance of stock projections in New England groundfish stocks. Canadian Journal of Fisheries and Aquatic Sciences, 73, 935–950. doi:10.1139/cjfas-2015-0163.
ICES (2018). Guidelines for calculating Mohn's rho: Retrospective bias in assessment. Draft document version 7 (2018-04-03), available at the Expert Groups area on the ICES Sharepoint.
ICES (2020). Workshop on Catch Forecast from Biased Assessments (WKFORBIAS; outputs from 2019 meeting). ICES Scientific Reports 2(28). doi:10.17895/ices.pub.5997.
Mohn, R. (1999). The retrospective problem in sequential population analysis: An investigation using cod fishery and simulated data. ICES Journal of Marine Science, 56, 473–488. doi:10.1006/jmsc.1999.0481.
shake is a retrospective example table.
icesAdvice-package gives an overview of the package.
mohn(shake) mohn(shake, details=TRUE) mohn(shake, plot=TRUE) mohn(shake, peels=3, plot=TRUE, col="black", ylim=0:1, yaxs="i") lines(as.numeric(rownames(shake)), shake$base, lwd=3) ## Plot last 10 years x <- rbind(matrix(1,28,6,dimnames=list(1981:2008,names(shake))), shake) mohn(tail(x, 10), plot=TRUE, lwd=2, main="main")mohn(shake) mohn(shake, details=TRUE) mohn(shake, plot=TRUE) mohn(shake, peels=3, plot=TRUE, col="black", ylim=0:1, yaxs="i") lines(as.numeric(rownames(shake)), shake$base, lwd=3) ## Plot last 10 years x <- rbind(matrix(1,28,6,dimnames=list(1981:2008,names(shake))), shake) mohn(tail(x, 10), plot=TRUE, lwd=2, main="main")
Read results from the DLS3.2 advisory method from a file into a list.
read.dls(file)read.dls(file)
file |
a filename. |
A list containing advice and other elements showing intermediate steps
in the calculations.
Arni Magnusson.
write.dls writes DLS3.2 results to a file.
DLS3.2 can be used to calculate catch advice for data-limited
stocks (DLS).
icesAdvice-package gives an overview of the package.
## Not run: survey <- data.frame(year=2001:2010, randu[1:10,]) dls <- icesAdvice::DLS3.2(1000, survey$y) write.dls(dls, "dls.txt") read.dls("dls.txt") file.remove("dls.txt") ## End(Not run)## Not run: survey <- data.frame(year=2001:2010, randu[1:10,]) dls <- icesAdvice::DLS3.2(1000, survey$y) write.dls(dls, "dls.txt") read.dls("dls.txt") file.remove("dls.txt") ## End(Not run)
Retrospective estimates of Southern hake fishing mortality.
shakeshake
Data frame containing 6 columns:
base |
base model estimates |
-1 |
1st retro peel |
-2 |
2nd retro peel |
-3 |
3rd retro peel |
-4 |
4th retro peel |
-5 |
5th retro peel |
This dataset is an example from the ICES (2018) Advice Technical Guidelines on quantifying and reporting retrospective bias.
ICES (2018). Guidelines for calculating Mohn's rho: Retrospective bias in assessment. Draft document version 7 (2018-04-03), available at the Expert Groups area on the ICES Sharepoint.
mohn calculates Mohn's rho.
icesAdvice-package gives an overview of the package.
shake mohn(shake)shake mohn(shake)
Calculate the implicit sigma that was used to construct a confidence interval.
sigmaCI(lo, hi, log = TRUE, level = 0.95)sigmaCI(lo, hi, log = TRUE, level = 0.95)
lo |
the lower confidence bound. |
hi |
the upper confidence bound. |
log |
whether the confidence interval is lognormal. |
level |
the confidence level. |
Implicit value of sigma.
The purpose of PA reference points is to apply a precautionary approach in fisheries management.
This function is useful for reviewing PA reference points, when the report provides a CI but not the value of sigma.
Arni Magnusson.
sigmaPA calculates the implicit sigma from PA reference points.
icesAdvice-package gives an overview of the package.
sigmaCI(100, 200)sigmaCI(100, 200)
Calculate the implicit sigma that was used to calculate PA reference points from limit reference points (Xpa from Xlim).
sigmaPA(lim, pa)sigmaPA(lim, pa)
lim |
the value of the limit reference point, e.g., Blim or Flim. |
pa |
the value of the PA reference point, e.g., Bpa or Fpa. |
The order of the parameters does not matter, so sigmaPA(Fpa, Flim) and
sigmaPA(Flim, Fpa) are equivalent.
Implicit value of sigma.
The purpose of PA reference points is to apply a precautionary approach in fisheries management.
This function is useful for reviewing PA reference points, when the advice sheet provides the value of Xlim and Xpa but not the value of sigma.
The inference is based on the following relationships:
Arni Magnusson.
sigmaCI calculates the implicit sigma from a confidence
interval.
Bpa and Fpa calculate those reference points from
the limit reference points, based on a given sigma.
icesAdvice-package gives an overview of the package.
sigmaPA(100, 120)sigmaPA(100, 120)
Write results from the DLS3.2 advisory method to a file.
write.dls(x, file = "")write.dls(x, file = "")
x |
a list generated by the |
file |
a filename. |
The resulting text file has Dos line endings (CRLF).
Arni Magnusson.
read.dls reads DLS3.2 results from a file back into R.
DLS3.2 can be used to calculate catch advice for data-limited
stocks (DLS).
icesAdvice-package gives an overview of the package.
## Not run: survey <- data.frame(year=2001:2010, randu[1:10,]) dls <- icesAdvice::DLS3.2(1000, survey$y) write.dls(dls, "dls.txt") read.dls("dls.txt") file.remove("dls.txt") ## End(Not run)## Not run: survey <- data.frame(year=2001:2010, randu[1:10,]) dls <- icesAdvice::DLS3.2(1000, survey$y) write.dls(dls, "dls.txt") read.dls("dls.txt") file.remove("dls.txt") ## End(Not run)