knitr::opts_chunk$set(eval = TRUE) # set to FALSE locally if you don't want to build when developing
The msy
R package is a collection of methods to estimate equilibrium reference points for fish stocks
You are welcome to:
This document is as much as the msy
-package itself a living document.
The origin of this package is from an intial coding by John Simmonds which was restructured by Colin Millar into an R-package with additional development, including coding the Buckland method. Einar Hjörleifsson compartmentalized the structure of the code as well as providing the output of the analysis in a more structured format. Further contributions have been made by Carmen Fernandez, Max Cardinale, Martin Pastoors and Niels Hintzen.
The package is currently being maintianed by the ICES secretariat.
The developmental repository for the msy
package is located on github, more specifically on github.com/ices-tools-prod/msy
The easiest way to install the msy
package is to use the function install_github
in the devtools
package. Since the input into the msy is in the form of a FLStock
object it is advised to obtain the latest release of FLCore
from the same site.
The Rtools.exe
software is needed for building packages under Microsoft Windows.
Run the following lines to install the latest versions of msy
and FLCore
.
# install.packages("devtools") install.packages("FLCore", repo = "http://flr-project.org/R") library(devtools) install_github("ices-tools-prod/msy")
devtools::install("D:/projects/git/ices-tools-prod/msy", quiet = TRUE) set.seed(305934095)
The above is equivalent to install.packages
and hence need only to be performed once. However, since the msy
package is currently under development (including bug-fixing) one may expect more frequent code updating in the package than what one may be familiar with for packages on CRAN
. Once installed the msy
packages is simply loaded via the familiar:
library(msy)
Besides functions the package comes with the following data:
FLStock
object of the Celtic Sea codFLStock
object of the Eastern Baltic codFLStock
object of the Icelandic codFLStock
object of the North Sea codFLStock
object of the Western Baltic codFLStock
object of the West of Scotland codFLStock
object of the Faeroe saitheFLStock
object of the Icelandic saitheFLStock
object of the Kattegat soleThese are all stored in the icesStocks list object.
The current version of the msy
implements two methods that go under the working names eqsim
and plotMSY
. Only usage of functions for the eqsim
approaches are described in the following sections.
eqsim
is a stochastic equilibrium software that may be used to explore MSY reference points. Productivity parameters (i.e. year vectors for natural mortality, weights-at-age and maturities) as well as selection are re-sampled at random from user specified range of years from the assessment. Fixing these parameters to an average over specified years can also be set by the user. Recruitments are re-sampled from their predictive distribution. Uncertainty in the stock-recruitment model is taken into account by applying model averaging using smooth AIC weights (Buckland et al. 1997). In addition assessment errors can be emulated by applying a user-specified error (CV and autocorrelation) to the intended target fishing mortality.
The current version of eqsim
only takes FLStock
objects as inputs.
In the following subsections we will simulate the north sea cod stock into the future under some basic assumptions. For the simulations we need to choose which years we will use to generate noise in the quantities: weight at age, maturity at age, natural mortality at age, and selection pattern. We also need to choose a set of Fbar values to simulate over in order estimate F reference points.
The eqsim approach consists of three components:
This can be done in one go with the following code:
FIT <- eqsr_fit(icesStocks$saiNS, nsamp = 1000, models = c("Ricker", "Segreg") ) SIM <- eqsim_run(FIT, bio.years = c(2004, 2013), sel.years = c(2004, 2013), Fcv = 0.24, Fphi = 0.42, Blim = 106000, Bpa = 200000, Fscan = seq(0, 1.2, len = 40), verbose = FALSE )
The stock recruitment function can be plotted by:
eqsr_plot(FIT,n=2e4)
Summary of the key results can be obtained by:
SIM$Refs
Summary plots conditioned on maximizing catch are obtained by:
eqsim_plot(SIM,catch=TRUE)
Summary plots of Fmsy range (WKMSYREF3) are obtained by:
eqsim_plot_range(SIM, type="mean") eqsim_plot_range(SIM, type="median") eqsim_plot_range(SIM, type="ssb")
Model fitting is done by maximum likelihood using the nlminb
optimiser in R. By refitting to non-parametric bootstrap resamples of the stock and recruit pairs, samples from the approximate joint distribution of the model parameters can be made. This is done by invoking the eqrs_fit
function. The function first sets up the stock and recruit pairs based on the information in the FLStock
object and removes any incomplete pairs, before dispatching on the model fitting / averaging algorithm chosen. Currently only a bootstrap based model averaging method called smooth AIC is implemented fully. The details can be found in eqrs_Buckland
function. The algorithm implemented is:
This process provides a robust way to average over several models, as long as the bootstrap resampling procedure provides an adequate approximation to the empirical distribution of the stock and recruit pairs.
The arguments to the fitting function are
args(eqsr_fit)
Here:
stk
is an FLStock
objectnsamp
is the number of simulations to run (often referred to as iterations)models
is the models to average over (any of the combination of these can be supplied, including only a single model)method
the method used (only Buckland as of now)id.sr
placeholder if one wants to name the fitremove.years
is used to remove years from the fitdelta
and nburn
are related to an MCMC based fitting procedure (not implemented)The results from the fitting process are returned to the user as a list:
str(FIT, 2, give.attr=FALSE)
where:
sr.sto
is the the (joint) stochastic distribution of the estimated model and parameters. The number of rows of the data frame is equivalent to the value supplied to nsamp
in the eqsr_fit
function.sr.det
is the conventional determinimstic predictive estimate. The n
indicates the number of times a particular function is drawn in the stochastic sample and the prop
the proportion, given nsamp
.pRec
contains the fitted parameters to the observed datastk
retains the original FLStock
object passed to the function.rby
(results by year) contains a summary of the ssb and rec data used in the fitting as well as other stock summary information used later down the lineThe contribution of each of the models can be obtained by printing out the sr.det
:
FIT$sr.det
Here the a, b and cv are the estimated parameters from the deterministic fit for each model. The n
and prop
is a summary of the number and proportion that each model contributes to the final fit.
Again to obtain a plot one simply calls:
eqsr_plot(FIT,n=2e4)
The n
supplied to the eqsr_plot
stands here for the number of stochastic recruitment points desired to include in the plot. The various black dashed lines represent the best fit of the different recruitment models and the yellow and blue lines the median and 5% and 95% percentiles of the distributions of the stochastic recruits drawn from the models. The input data are represented by red points.
An alternative to the base
plot is a ggplot2
version (with too many fancy colours :-):
eqsr_plot(FIT,n=2e4,ggPlot=TRUE,Scale=1e3)
Here the model fits are represented in different colours with the yellow lines indicating the 5%, 50% and 95% percentiles of the stochastic recruitment distribution. The input data are represented by text indicating year class. The weight of each model in the final stochastic recruitment draw is indicated as a proportion in the legends and by different colours for the stochastic draws.
Simulating forward is done using the eqsim_run
function. The function takes as input the output from the eqsr_fit
function. Simulations are run independently for each sample from the distribution of model and parameters. This is done for a range of $F_{advisory}$ values. For example if we scanned over 10 values of $F_{advisory}$ and had taken 2000 samples from the stock-recruit relationship then 20000 simulations would be run in total. These simulations are run for 200 years (default, specified with Nrun
), and the last 50 years are retained to calculate summaries, like the proportion of times the stock crashes at a given $F_{advisory}$. It is important to note that each simulation is conditioned on a single stock recruit relationship with fixed parameters and cv.
Error is introduced within the simulations by generating process error about the constant stock-recruit fit, and by using variation in maturity, natural mortality, weight at age and selection estimates. Note that if there is no variability in these quantities in the stock object then no variability will be taken in to the simulations. The user can also specify using average values for these parameters.
The arguments to the simulation function are:
args(eqsim_run)
where:
fit
is the output list from eqsr_fit
bio.years
is the start and end year from which to generate noise in maturity, M and weights.bio.const
is a flag indicating if the average maturity, M and weights over the specified years should be used (TRUE
) or not (FALSE
).sel.years
is the start and end year from which to generated noise in the selection at agesel.const
is a flag indicating if the average selection over the specified years should be used (TRUE
) or not (FALSE
).Fscan
is the range of $F_{advisory}$ values to scan overBtrigger
is the location of a modifier of a HCR upon which $F_{advisory}$ becomes linearily reduced. If Btrigger
is 0 (default) this is equivalent to a constant F-rule.Fcv
The assessment error of fishing mortality in the advisory year.Fphi
The autocorrelation of fishing mortality in assessment errorSSBcv
The assessment error in the spawning stock in the advisory yearrhologRec
Use or not use autocorrelation in recruitment residuals.Blim
$B_{lim}$Bpa
$B_{pa}$Nrun
is the number of years to simulate forward (fixed for now is that the last 50 years from those are used for summarising equilibrium conditions)verbose
controls if progress bar is displayed during the simulationextreme.trim
A numerical vector of length 2 containing the lower and upper percentiles. If specified, recruitement values outside this range are trimmed (ignored).The results from the simulation process are returned to the user as a list
str(SIM, 2, give.attr = FALSE)
where
ibya
(input by year and age) contains the biological and fisheries input data.rby
(results by year) contains the stock summary data.rbp
(results by probability) contains the 0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975 percentiles of the simulations of SSB, catch, landings and recruitment for each Fscan value.Blim
$B_{lim}$ input valueBpa
$B_{pa}$ input valueRefs
Calculated reference pointspProfile
The probability profiles for a given target F for $B_{lim}$, $B_{pa}$ and $F_{msy}$ (both for catch and landings).refs_interval
The $F_{msy}$ interval and 5% percentile risk of going below $B_{lim}$ (per WKFMSYREF3)The calculation associated with the $F_{msy}$ range values can be accessed simply by:
t(SIM$refs_interval)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.