Preamble

This document is as much as the msy-package itself still in developement.

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.

At this moment there is no person responsible for the maintenance of the package, including the above mentioned names.

Installation

The developmental repository for the msy package is located on github, more specifically on github.com/wgmg/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")
library(devtools)
install_github("wgmg/msy")
install_github("flr/FLCore")

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:

codCS: FLStock object of the Celtic Sea cod
codEB: FLStock object of the Eastern Baltic cod
codIS: FLStock object of the Icelandic cod
codNS: FLStock object of the North Sea cod
codWB: FLStock object of the Western Baltic cod
codWS: FLStock object of the West of Scotland cod
saiFO: FLStock object of the Faeroe saithe
saiIS: FLStock object of the Icelandic saithe
solKA: FLStock object of the Kattegat sole

These 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

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.

A quick start

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:

Estimate the stock recruitment relationship
Simulate a stock to equilibrium and continue simulating for some years
Calculate reference points from the simulated stock at equilibrium (last 50 years of the runs are used)

This can be done in one go with the following code:

FIT <- eqsr_fit(icesStocks$codNS,
                nsamp = 1000, 
                models = c("ricker", "smooth_hockey", "bevholt"))
SIM <- eqsim_run(FIT, Fcv=0.25, Fphi=0.30,
                 Blim=70000,Bpa=150000,
                 Fscan = seq(0,1.2,len=40),
                 verbose=FALSE,extreme.trim=c(0.05,0.95))

The stock recruitment function can be plotted by:

eqsr_plot(FIT,n=2e4)

plot of chunk unnamed-chunk-2

Summary of the key results can be obtained by:

SIM$Refs

               F05       F10       F50 medianMSY   meanMSY  FCrash05
catF     7.165e-01 7.373e-01 7.972e-01 2.769e-01 3.077e-01 7.692e-01
lanF            NA        NA        NA 2.154e-01 2.462e-01        NA
catch    3.232e+05 2.990e+05 1.285e+05 4.275e+05 4.259e+05 2.301e+05
landings        NA        NA        NA 3.618e+05 3.598e+05        NA
catB     2.149e+05 1.888e+05 7.037e+04 1.045e+06 9.175e+05 1.331e+05
lanB            NA        NA        NA 1.383e+06 1.198e+06        NA
          FCrash50
catF        0.8615
lanF            NA
catch    1820.9747
landings        NA
catB      896.1669
lanB            NA

Summary plots conditioned on maximizing catch are obtained by:

eqsim_plot(SIM,catch=TRUE)

plot of chunk unnamed-chunk-4

Summary plots of Fmsy range (WKMSYREF3) are obtained by:

eqsim_plot_range(SIM)

plot of chunk unnamed-chunk-5

The recruitment model

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:

Take a resample with replacement from the stock and recruit pairs
Fit every stock-recruit model under consideration and store the AIC of each
Retain the parameter estimates from the best model
Repeat

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)

function (stk, nsamp = 5000, models = c("ricker", "segreg", "bevholt"), 
    method = "Buckland", id.sr = NULL, remove.years = NULL, delta = 1.3, 
    nburn = 10000) 
NULL

Here:

stk is an FLStock object
nsamp 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 fit
remove.years is used to remove years from the fit
delta 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)

List of 6
 $ sr.sto:'data.frame': 1000 obs. of  4 variables:
  ..$ a    : num [1:1000] 3.2 3 6.13 5.63 3.13 ...
  ..$ b.b  : num [1:1000] 1.60e+05 1.83e+05 5.19e-07 2.89e-07 1.64e+05 ...
  ..$ cv   : num [1:1000] 0.412 0.453 0.415 0.467 0.501 ...
  ..$ model: chr [1:1000] "smooth_hockey" "smooth_hockey" "bevholt" "ricker" ...
 $ sr.det:'data.frame': 3 obs. of  6 variables:
  ..$ a    : num [1:3] 6.3 3.04 6.25
  ..$ b    : num [1:3] 8.54e-07 1.77e+05 8.57e-07
  ..$ cv   : num [1:3] 0.482 0.463 0.482
  ..$ model: Factor w/ 3 levels "ricker","smooth_hockey",..: 1 2 3
  ..$ n    : int [1:3] 77 901 22
  ..$ prop : num [1:3] 0.077 0.901 0.022
 $ pRec  : num [1:1000, 1:49] 139906 131183 132419 122157 136642 ...
 $ stk   :Formal class 'FLStock' [package "FLCore"] with 20 slots
 $ rby   :'data.frame': 49 obs. of  6 variables:
  ..$ year    : int [1:49] 1963 1964 1965 1966 1967 1968 1969 1970 1971 1972 ...
  ..$ rec     : num [1:49] 845768 1067681 1375049 1274418 654744 ...
  ..$ ssb     : num [1:49] 153638 165830 205112 228117 252019 ...
  ..$ fbar    : num [1:49] 0.478 0.509 0.541 0.556 0.609 ...
  ..$ landings: num [1:49] 109409 117035 167733 207369 235070 ...
  ..$ catch   : num [1:49] 128686 130740 210237 259416 276387 ...
 $ id.sr : chr "North Sea cod"

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 data
stk 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 line
id.rs` is the user specified id

The contribution of each of the models can be obtained by printing out the sr.det:

FIT$sr.det

      a         b     cv         model   n  prop
1 6.296 8.541e-07 0.4820        ricker  77 0.077
2 3.043 1.771e+05 0.4626 smooth_hockey 901 0.901
3 6.253 8.574e-07 0.4825       bevholt  22 0.022

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)

plot of chunk unnamed-chunk-9

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)

plot of chunk unnamed-chunk-10

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.

The simulation

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)

function (fit, bio.years = c(2008, 2012), bio.const = FALSE, 
    sel.years = c(2008, 2012), sel.const = FALSE, Fscan = seq(0, 
        1, len = 20), Fcv = 0, Fphi = 0, SSBcv = 0, rhologRec = FALSE, 
    Blim, Bpa, recruitment.trim = c(3, -3), Btrigger = 0, Nrun = 200, 
    process.error = TRUE, verbose = TRUE, extreme.trim) 
NULL

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 age
sel.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 over
Btrigger 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 error
SSBcv The assessment error in the spawning stock in the advisory year
rhologRec 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)
`process.error} allows the simulations to be run using the predictive distribition of recruitment or the mean recruitment
verbose controls if progress bar is displayed during the simulation
extreme.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)

List of 10
 $ ibya         :List of 7
  ..$ Mat  : num [1:7, 1:5] 0.01 0.05 0.23 0.62 0.86 1 1 0.01 0.05 0.23 ...
  ..$ M    : num [1:7, 1:5] 1.039 0.696 0.487 0.232 0.2 ...
  ..$ Fprop: Named num [1:7] 0 0 0 0 0 0 0
  ..$ Mprop: Named num [1:7] 0 0 0 0 0 0 0
  ..$ west : num [1:7, 1:5] 0.33 0.904 1.971 3.834 5.692 ...
  ..$ weca : num [1:7, 1:5] 0.33 0.904 1.971 3.834 5.692 ...
  ..$ sel  : num [1:7, 1:5] 0.246 0.795 1.087 1.118 1.142 ...
 $ rbya         :List of 1
  ..$ ferr: num [1:40, 1:50, 1:1000] 0.174 0.174 0.174 0.174 0.174 ...
 $ rby          :'data.frame':  49 obs. of  6 variables:
  ..$ year    : int [1:49] 1963 1964 1965 1966 1967 1968 1969 1970 1971 1972 ...
  ..$ rec     : num [1:49] 845768 1067681 1375049 1274418 654744 ...
  ..$ ssb     : num [1:49] 153638 165830 205112 228117 252019 ...
  ..$ fbar    : num [1:49] 0.478 0.509 0.541 0.556 0.609 ...
  ..$ landings: num [1:49] 109409 117035 167733 207369 235070 ...
  ..$ catch   : num [1:49] 128686 130740 210237 259416 276387 ...
 $ rbp          :'data.frame':  160 obs. of  10 variables:
  ..$ Ftarget : num [1:160] 0 0.0308 0.0615 0.0923 0.1231 ...
  ..$ variable: chr [1:160] "Recruitment" "Recruitment" "Recruitment" "Recruitment" ...
  ..$ p025    : num [1:160] 427445 434452 440630 444548 446988 ...
  ..$ p05     : num [1:160] 507608 513592 519459 522581 524475 ...
  ..$ p25     : num [1:160] 822383 826736 830537 832592 834781 ...
  ..$ p50     : num [1:160] 1158240 1163868 1168463 1172209 1176053 ...
  ..$ p75     : num [1:160] 1674294 1686354 1697002 1706038 1714254 ...
  ..$ p95     : num [1:160] 4505939 4815581 5166433 5545947 6050406 ...
  ..$ p975    : num [1:160] 31719521 35322635 39346841 42457322 46962108 ...
  ..$ Mean    : num [1:160] NA NA NA NA NA NA NA NA NA NA ...
 $ Blim         : num 70000
 $ Bpa          : num 150000
 $ Refs         : num [1:6, 1:7] 7.17e-01 NA 3.23e+05 NA 2.15e+05 ...
 $ pProfile     :'data.frame':  1104 obs. of  3 variables:
  ..$ Ftarget : num [1:1104] 0.169 0.17 0.17 0.171 0.172 ...
  ..$ value   : num [1:1104] 1.28e-05 3.20e-05 6.04e-05 1.01e-04 1.59e-04 ...
  ..$ variable: Factor w/ 4 levels "pFmsyCatch","pFmsyLandings",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ id.sim       : chr "North Sea cod"
 $ refs_interval:'data.frame':  1 obs. of  8 variables:
  ..$ FmsyMedianC     : num 0.277
  ..$ FmsylowerMedianC: num 0.169
  ..$ FmsyupperMedianC: num 0.452
  ..$ FmsyMedianL     : num 0.217
  ..$ FmsylowerMedianL: num 0.145
  ..$ FmsyupperMedianL: num 0.332
  ..$ F5percRiskBlim  : num 0.717
  ..$ Btrigger        : num 0

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 value
Bpa \(B_{pa}\) input value
Refs Calculated reference points
pProfile 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)

##                    [,1]
## FmsyMedianC      0.2774
## FmsylowerMedianC 0.1688
## FmsyupperMedianC 0.4523
## FmsyMedianL      0.2171
## FmsylowerMedianL 0.1447
## FmsyupperMedianL 0.3317
## F5percRiskBlim   0.7165
## Btrigger         0.0000

Using the eqsim routines in the msy package

John Simmond, Colin Millar and Einar Hjörleifsson

2015-03-14