Title: | Estimation of Equilibrium Reference Points for Fsisheries |
---|---|
Description: | Methods to estimate equilibrium reference points for fisheries data. Currently data must be converted into FLStock objects of the FLR (Fisheries Library in R) style, defined in the R package FLCore. |
Authors: | John Simmonds [aut], Einar Hjorleifsson [aut], Carmen Fernandez [ctb], Colin Millar [aut, cre] |
Maintainer: | Colin Millar <[email protected]> |
License: | GPL (>= 2) |
Version: | 0.1.19 |
Built: | 2024-10-23 02:58:27 UTC |
Source: | https://github.com/ices-tools-prod/msy |
Methods to estimate equilibrium reference points for fisheries data. Currently data must be converted into FLStock objects of the FLR (Fisheries Library in R) style, defined in the R package FLCore
Model fitting and simulation:
eqsr_fit |
fitting stock recruit models to data |
eqsim_run |
simulation am 'equilibrium' population state |
Plotting:
eqsr_plot |
plot stock recuitment fit |
eqsim_plot |
plot summary of simulation showing reference points |
eqsim_plot_range |
plot summary of MSY ranges reference points |
Example data:
icesStocks |
A list of various stocks |
John Simmonds, Einar Hjorleifsson, Carmen Fernandez and Colin Millar.
ICES (2015) Report of the Workshop to consider F MSY ranges for stocks in ICES categories 1 and 2 in Western Waters (WKMSYREF4). 01 WKMSYREF4 Report.pdf
ICES (2017) ICES fisheries management reference points for category 1 and 2 stocks. DOI: 10.17895/ices.pub.3036
Buckland, S.T., K.P. Burnham & N.H. Augustin (1997). Model selection: An integral part of inference. Biometrics 53, 603-618. DOI: 10.2307/2533961
To explore the code of the package see the GitHub repo: ices-tools-prod/msy
stock recruitment function
Bevholt(ab, ssb)
Bevholt(ab, ssb)
ab |
the model parameters |
ssb |
a vector of ssb |
log recruitment according to model
stock recruitment function
bevholt2(ab, ssb)
bevholt2(ab, ssb)
ab |
the model parameters |
ssb |
a vector of ssb |
log recruitment according to model
XXX
eqsim_ggplot(sim, Scale = 1, plotit = TRUE)
eqsim_ggplot(sim, Scale = 1, plotit = TRUE)
sim |
An object returned from the function eqsim_run |
Scale |
A value, the scaling on the yaxis |
plotit |
Boolean, if TRUE (default) returns a plot |
Einar Hjorleifsson [email protected]
XXX
eqsim_plot(sim, ymax.multiplier = 1.2, catch = TRUE)
eqsim_plot(sim, ymax.multiplier = 1.2, catch = TRUE)
sim |
An object returned from the function eqsim_run |
ymax.multiplier |
A value that acts as a multiplier of the maximum observed variable being plotted. E.g. 1.2 means that for each of the three panels a, b and c the ymax is set to 1.2 of the maximum observed recruitment, spawning stock biomass and yield (catch or landings, depending on user input. |
catch |
Boolean, if TRUE (default) returns a plot based on catch. If false returns a plot based on landings. |
Einar Hjorleifsson [email protected]
XXX
eqsim_plot_range(sim, interval = 0.95, type = "median")
eqsim_plot_range(sim, interval = 0.95, type = "median")
sim |
XXX |
interval |
XXX |
type |
XXX |
Simulate a fish stock forward in time given biological parameters, fishery parameters and advice parameters.
eqsim_run( fit, bio.years = c(-5, -1) + FLCore::dims(fit$stk)$maxyear, bio.const = FALSE, sel.years = c(-5, -1) + FLCore::dims(fit$stk)$maxyear, sel.const = FALSE, Fscan = seq(0, 2, len = 40), Fcv = 0, Fphi = 0, SSBcv = 0, rhologRec = TRUE, Blim, Bpa, recruitment.trim = c(3, -3), Btrigger = 0, Nrun = 200, process.error = TRUE, verbose = TRUE, extreme.trim = c(0, 1), R.initial = mean(fit$rby$rec), keep.sims = FALSE )
eqsim_run( fit, bio.years = c(-5, -1) + FLCore::dims(fit$stk)$maxyear, bio.const = FALSE, sel.years = c(-5, -1) + FLCore::dims(fit$stk)$maxyear, sel.const = FALSE, Fscan = seq(0, 2, len = 40), Fcv = 0, Fphi = 0, SSBcv = 0, rhologRec = TRUE, Blim, Bpa, recruitment.trim = c(3, -3), Btrigger = 0, Nrun = 200, process.error = TRUE, verbose = TRUE, extreme.trim = c(0, 1), R.initial = mean(fit$rby$rec), keep.sims = FALSE )
fit |
A list returned from the function fitModels |
bio.years |
The years to sample maturity, weights and M from, given as a vector of length 2, i.e. c(2010, 2015) select from the years 2010 to 2015 inclusive. |
bio.const |
A flag (default FALSE), if TRUE mean of the biological values from the years selected are used |
sel.years |
The years to sample the selection patterns from, given as a vector of length 2, i.e. c(2010, 2015) select from the years 2010 to 2015 inclusive. |
sel.const |
A flag (default FALSE), if TRUE mean of the selection patterns from the years selected are used |
Fscan |
F values to scan over, i.e. seq(0, 2, by = 0.05) |
Fcv |
Assessment error in the advisory year |
Fphi |
Autocorrelation in assessment error in the advisory year |
SSBcv |
Spawning stock biomass error in the advisory year |
rhologRec |
A flag for recruitment autocorrelation, default (TRUE), or a vector of numeric values specifcying the autocorrelation parameter for the residuals for each SR model. |
Blim |
SSB limit reference point |
Bpa |
SSB precuationary reference point |
recruitment.trim |
A numeric vector with two log-value clipping the extreme recruitment values from a continuous lognormal distribution. The values must be set as c("high","low"). |
Btrigger |
If other than 0 (default) the target F applied is reduced by SSB/Btrigger. This is the "ICES Advice Rule". |
Nrun |
The number of years to run in total (the last 50 years from that will be retained to compute equilibrium values from) |
process.error |
Use stochastic recruitment or mean recruitment? TRUE (default) uses the predictive distribution of recruitment, model estimate of recruitment + simulated observation error. FALSE uses model prediction of recruitment with no observation error. |
verbose |
Flag, if TRUE (default) indication of the progress of the simulation is provided in the console. Useful to turn to FALSE when knitting documents. |
extreme.trim |
a pair of quantiles (low, high) which are used to trim
the equilibrium catch values, across simulations within
an F scenario, when calculating the mean catch and
landings for that F scenario. These mean values
calculated accross simulations within an F scenario
are used to find which F scenario gave the maximum catch.
|
R.initial |
Initial recruitment for the simulations. This is common accross all simulations. Default = mean of all recruitments in the series. |
keep.sims |
Flag, if TRUE returns a matrix of population tragectories for each value of F in Fscan (see examples). |
Details of the steps required to evaluate reference points are given in ICES (2017). WHile, details of the calculation of MSY ranges is given in ICES (2015).
A list containing the results from the forward simulation and the reference points calculated from it.
ICES (2015) Report of the Workshop to consider F MSY ranges for stocks in ICES categories 1 and 2 in Western Waters (WKMSYREF4). 01 WKMSYREF4 Report.pdf
ICES (2017) ICES fisheries management reference points for category 1 and 2 stocks. DOI: 10.17895/ices.pub.3036
eqsr_fit
fits multiple stock recruitment models to a data set.
eqsr_plot
plots the results from eqsr_fit.
eqsim_plot
summary plot of the forward simulation showing estimates
of various reference points.
eqsim_plot_range
summary plots of the forward simulation showing
the estimates of MSY ranges (ICES, 2015)
msy-package
gives an overview of the package.
## Not run: data(icesStocks) 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) ) # extract tragectories ssbsim <- SIM$rbya$ssb years <- SIM$rbya$simyears models <- SIM$rbya$srmodels$model Ftarget <- SIM$rbya$Ftarget Fval <- which(Ftarget == 0) Fval <- which(Ftarget > .3)[1] x <- ssbsim[Fval,,] df <- data.frame(year = 1:nrow(x), ssb = c(x), sim = rep(1:ncol(x), each = nrow(x)), model = rep(models, each = nrow(x))) xyplot(ssb ~ year | model, groups = sim, data = df, type = "l", col = grey(0.5, alpha = 0.5)) fit <- density(x[x>1e-3], from = 0) plot(fit$x,fit$y*mean(x>1e-3),col="red", type = "l") lines(x = 0, y = mean(x<=1e-3), type = "h", lwd = 3) ## End(Not run)
## Not run: data(icesStocks) 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) ) # extract tragectories ssbsim <- SIM$rbya$ssb years <- SIM$rbya$simyears models <- SIM$rbya$srmodels$model Ftarget <- SIM$rbya$Ftarget Fval <- which(Ftarget == 0) Fval <- which(Ftarget > .3)[1] x <- ssbsim[Fval,,] df <- data.frame(year = 1:nrow(x), ssb = c(x), sim = rep(1:ncol(x), each = nrow(x)), model = rep(models, each = nrow(x))) xyplot(ssb ~ year | model, groups = sim, data = df, type = "l", col = grey(0.5, alpha = 0.5)) fit <- density(x[x>1e-3], from = 0) plot(fit$x,fit$y*mean(x>1e-3),col="red", type = "l") lines(x = 0, y = mean(x<=1e-3), type = "h", lwd = 3) ## End(Not run)
Fits one or more stock recruitment relationship to data containted in an FLStock object. If more than one stock recruit relationship is provided, the models are weighted based on smooth AIC weighting (See Buckland et al., 1997).
eqsr_fit( stk, nsamp = 1000, models = c("Ricker", "Segreg", "Bevholt"), id.sr = FLCore::name(stk), remove.years = NULL, rshift = 0 )
eqsr_fit( stk, nsamp = 1000, models = c("Ricker", "Segreg", "Bevholt"), id.sr = FLCore::name(stk), remove.years = NULL, rshift = 0 )
stk |
FLStock object |
nsamp |
Number of samples (iterations) to take from the stock recruitment fit (default is 1000). If 0 (zero) then only the fits to the data are returned and no simulations are made. |
models |
A character vector containing stock recruitment models to use in the model averaging. User can set any combination of "Ricker", "Segreg", "Bevholt", "Smooth_hockey". |
id.sr |
A character vector specifying an id or name for the stock recruitment fit being run. The default is to use the slot "name" in the stk parameter is provided |
remove.years |
A vector specifying the years to remove from the model fitting. |
rshift |
lag ssb by aditional years (default = 0). As an example, for some herring stocks, age 1 (1 winter ring) fish were spawned 2 years previously, in this case, rshift = 1. |
A list containing the following objects:
'sr.sto' data.frame containing the alpha (a), beta (b), cv and model names. The number of rows correspond to the value set of 'nsamp' in the function call.
'sr.det' The parameters in the stock recruitment model corresponding to the "best fit" of any given model.
'stk' An FLStock object, same as provided as input by the user.
'rby' A data.frame containing the recruitment (rec), spawning stock biomass (ssb) and year used in the fitting of the data.
'id.sr' A string containing run name (taken from the 'id.sr' argument)
Buckland, S.T., K.P. Burnham & N.H. Augustin (1997). Model selection: An integral part of inference. Biometrics 53, 603-618. DOI: 10.2307/2533961
eqsr_plot
plots a simulation of predictive recruitment
from the fit, and shows a summary of the contributions of each stock
recruitment model to the model average fit.
data(icesStocks) FIT <- eqsr_fit(icesStocks$saiNS, nsamp = 0, models = c("Ricker", "Segreg")) # summary of individual fits FIT$sr.det eqsr_plot(FIT) # fit a bounded segmented regression Segreg_bounded <- function(ab, ssb) { ab$b <- min_ssb + ab$b Segreg(ab, ssb) } min_ssb <- min(FLCore::ssb(icesStocks$saiNS)) FIT <- eqsr_fit(icesStocks$saiNS, nsamp = 0, models = c("Segreg", "Segreg_bounded")) # summary of individual fits FIT$sr.det FIT$sr.det$b[2] + min_ssb eqsr_plot(FIT) ## Not run: FIT <- eqsr_fit(icesStocks$saiNS, nsamp = 2000, models = c("Segreg", "Segreg_bounded")) # summary of individual fits FIT$sr.det eqsr_plot(FIT) ## End(Not run)
data(icesStocks) FIT <- eqsr_fit(icesStocks$saiNS, nsamp = 0, models = c("Ricker", "Segreg")) # summary of individual fits FIT$sr.det eqsr_plot(FIT) # fit a bounded segmented regression Segreg_bounded <- function(ab, ssb) { ab$b <- min_ssb + ab$b Segreg(ab, ssb) } min_ssb <- min(FLCore::ssb(icesStocks$saiNS)) FIT <- eqsr_fit(icesStocks$saiNS, nsamp = 0, models = c("Segreg", "Segreg_bounded")) # summary of individual fits FIT$sr.det FIT$sr.det$b[2] + min_ssb eqsr_plot(FIT) ## Not run: FIT <- eqsr_fit(icesStocks$saiNS, nsamp = 2000, models = c("Segreg", "Segreg_bounded")) # summary of individual fits FIT$sr.det eqsr_plot(FIT) ## End(Not run)
Plot Simulated Predictive Distribution of Recruitment
eqsr_plot( fit, n = 20000, x.mult = 1.1, y.mult = 1.4, ggPlot = FALSE, Scale = 1 )
eqsr_plot( fit, n = 20000, x.mult = 1.1, y.mult = 1.4, ggPlot = FALSE, Scale = 1 )
fit |
an fitted stock recruit model returned from |
n |
Number of random recruitment draws to plot |
x.mult |
max value for the y axis (ssb) as a multiplier of maximum observed ssb |
y.mult |
max value for the x axis (rec) as a multiplier of maismum observed rec |
ggPlot |
Flag, if FALSE (default) plot using base graphics, if TRUE do a ggplot |
Scale |
Numeric value for scaling varibles in plot. |
NULL produces a plot
eqsr_fit
Fits several stock recruitment models to a data set
and calculates the proportion contribution of each model based on a bootstrap
model averaging procedure.
## Not run: data(icesStocks) FIT <- eqsr_fit(icesStocks$saiNS, nsamp = 1000, models = c("Ricker", "Segreg")) eqsr_plot(FIT, n = 20000) # Scale argument only available for ggPlot = TRUE eqsr_plot(FIT, n = 20000, ggPlot = TRUE, Scale = 1000) ## End(Not run)
## Not run: data(icesStocks) FIT <- eqsr_fit(icesStocks$saiNS, nsamp = 1000, models = c("Ricker", "Segreg")) eqsr_plot(FIT, n = 20000) # Scale argument only available for ggPlot = TRUE eqsr_plot(FIT, n = 20000, ggPlot = TRUE, Scale = 1000) ## End(Not run)
A list FLStock object for various stocks
icesStocks
icesStocks
a list
NA
NA
Starting values for hockey stick models
initial(model, data)
initial(model, data)
model |
XXX |
data |
XXX |
vector of starting values
the log-likelihood of the rectuit function
llik(param, data, model, logpar = FALSE)
llik(param, data, model, logpar = FALSE)
param |
the model parameters |
data |
the rec and ssb data |
model |
the stock recruit model to use |
logpar |
are the parameters on the log scale |
the log-likelihood
Non exported function, plots the prgress of an iterative procedure using "[=> ]", "[==> ]", etc.
loader(p)
loader(p)
p |
Value |
stock recruitment function
Ricker(ab, ssb)
Ricker(ab, ssb)
ab |
the model parameters |
ssb |
a vector of ssb |
log recruitment according to model
stock recruitment function
Segreg(ab, ssb)
Segreg(ab, ssb)
ab |
the model parameters |
ssb |
a vector of ssb |
log recruitment according to model
stock recruitment function
segreg2(ab, ssb)
segreg2(ab, ssb)
ab |
the model parameters |
ssb |
a vector of ssb |
log recruitment according to model
stock recruitment function
Smooth_hockey(ab, ssb, gamma = 0.1)
Smooth_hockey(ab, ssb, gamma = 0.1)
ab |
the model parameters |
ssb |
a vector of ssb |
gamma |
a smoother parameter |
log recruitment according to model