This function estimates MSY following Martell and Froese(2012).
catchmsy(year = NULL, catch = NULL, catchCV = NULL,
catargs = list(dist = "none", low = 0, up = Inf, unit = "MT"),
l0 = list(low = 0, up = 1, step = 0), lt = list(low = 0, up = 1,
refyr = NULL),
sigv = 0, k = list(dist = "unif", low = 0, up = 1, mean = 0, sd = 0),
r = list(dist = "unif", low = 0, up = 1, mean = 0, sd = 0),
M = list(dist = "unif", low = 0.2, up = 0.2, mean = 0, sd = 0),
nsims = 10000, catchout = 0, grout = 1,
graphs = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11),
grargs = list(lwd = 1, pch = 16, cex = 1, nclasses = 20, mains = " ",
cex.main = 1,
cex.axis = 1, cex.lab = 1),
pstats = list(ol = 1, mlty = 1, mlwd = 1.5, llty = 3, llwd = 1, ulty = 3,
ulwd = 1),
grtif = list(zoom = 4, width = 11, height = 13, pointsize = 10))vector containing the time series of numeric year labels.
vector containing the time series of catch data (in weight). Missing values are not allowed.
vector containing the time series of coefficients of variation associated with catch if resampling of catch is desired; otherwise, catchCV = NULL.
list arguments associated with resampling of catch. dist is the specification
of the resampling distribution to use ("none" = no resampling, "unif"=uniform, "norm" = normal,
and "lnorm" =log-normal). If "lnorm" is selected, catch is log transformed and standard deviation
on the log scale is calculated from the specified CVs using the relationship sdlog=sqrt(log(CV^2+1)).
low and up are the lower and upper limit of distribution (if truncation is desired).
unit is the weight unit of catch (used in graph labels; default="MT"). If "unif", the
catch must be incorporated in low and up arguments. For instance, if the
lower limit to sample is the value of catch, specify low=catch. If some maximum
above catch will be the upper limit, specify up=50*catch. The limits for each year will
be applied to catch internally.
list arguments for the relative biomass in year 1. low and up are the lower
and upper bounds of the starting value of relative biomass (in relation to k) in year 1. step
is the step increment to examine. If step=0, then l0 is randomly selected from a
uniform distribution using the lower and upper starting values. If step>0, then step increments
are used (in this case, the number of simulations (nsims) are used for each increment).
list arguments for the depletion level in the selected reference year (refyr).
low and up are the lower and upper bounds of depletion level in refyr.
refyr can range from the first year to the year after the last year of catch (t+1).
standard deviation of the log-normal random process error. signv = 0 for no
process error.
list arguments for the carrying capacity. dist is the statistical distribution name
from which to sample k. low and up are the lower and upper bounds of k
in the selected distribution.
mean and sd are the mean and standard deviation for selected distributions. The
following are valid distributions: "none", "unif" - uniform, "norm" - normal, "lnorm" - log-normal,
"gamma" - gamma, and "beta" - beta distributions. "unif" requires non-missing values for low
and up. "norm", "lnorm", "gamma" and "beta", require non-missing values for
low,up, mean and sd. If "lnorm" is used, mean and sd
must be on the natural log scale (keep low and up on the original scale). If
dist = "none", the mean is used as a fixed value.
list arguments for the intrinsic growth rate. dist is the statistical distribution name
from which to sample r. low and up are the lower and upper bounds of r
in the selected distribution. mean and sd are the mean and standard deviation for
selected distributions. Valid distributions are the same as in k. If dist = "none",
the mean is used as a fixed value.
list arguments for natural mortality. dist is the statistical distribution name from
which to sample M. low and up are the lower and upper bounds of M in the
selected distribution. mean and sd are the mean and standard deviation for selected
distributions. Valid distributions are the same as in k. M is used to determine exploitation
rate (Umsy) at MSY. If dist = "none", the mean is used as a fixed value.
number of Monte Carlos samples.
If resampling catch, save catch trajectories from each Monte Carlos simulation
- 0 = No (default), 1 = Yes.
numeric argument specifying whether graphs should be printed to console only (1) or to
both the console and TIF graph files (2).Use setwd before running function to direct .tif files
to a specific directory. Each name of each file is automatically determined.
vector specifying which graphs should be produced. 1 = line plot of observed catch versus
year,2 = point plot of plausible k versus r values, 3 = histogram of plausible r values,
4 = histogram of plausible k values, 5 = histogram of M values,
6 = histogram of MSY from plausible values of l0,k,r, and Bmsy/k, 7 = histogram of Bmsy from plausible
values of l0,k,r, and Bmsy/k, 8 = histogram of Fmsy from plausible values of l0,k,r, and Bmsy/k, 9 =
histogram of Umsy values from Fmsy and M, 10 = histogram of overfishing limit (OFL) in last year+1 values
from Umsys, and 11 = line plots of accepted and rejected biomass trajectores with median and 2.5th and 97.5th
percentiles (in red). Any combinations of graphs can be selected within c(). Default is all.
list control arguments for plotting functions. lwd is the line width for graph = 1 and 11,
pch and cex are the symbol character and character expansion value used in graph = 2,
nclasses is the nclass argument for the histogram plots (graphs 3-11), mains and
cex.main are the titles and character expansion values for the graphs, cex.axis is the
character expansion value(s) for the x and y-axis tick labels and cex.lab is the character
expansion value(s) for the x and y-axis labels. Single values of nclasses,mains,
cex.main,cex.axis, cex.lab are applied to all graphs. To change arguments for
specific graphs, enclose arguments within c() in order of the number specified in graphs.
list control arguments for plotting the mean and 95
and management quantities on respective graphs. ol = 0, do not overlay values on plots, 1 =
overlay values on plots. mlty and mlwd are the line type and line width of the mean value;
llty and llwd are the line type and line wdith of the 2.5
ulwd are the line type and line width of the 97.5
list arguments for the .TIF graph files. See tiff help file in R.
dataframe containing the initial values for each explored parameter.
dataframe containing the mean, median, 2.5th and 97.5 plausible (likelihood=1) parameters.
dataframe containing the mean, median, 2.5th and 97.5 of the management quantities (i.e., MSY, Bmsy, etc.) for the plausible parameters (likelihood=1)
dataframe containing the values of l0, k, r, Bmsy/k, M and associated management quantities for all (likelihood=0 and likelihood=1) random draws.
value of the last year of catch data + 1 for use in function dlproj.
designates the output object as a catchmsy object for use in function dlproj.
The biomass estimates from each simulation are not stored in memory but are automatically saved to a .csv file named "Biotraj-cmsy.csv". Yearly values for each simulation are stored across columns. The first column holds the likelihood values for each simulation (1= accepted, 0 = rejected). The number of rows equals the number of simulations (nsims). This file is loaded to plot graph 11 and it must be present in the default or setwd() directory. When catchout=1, catch values randomly selected are saved to a .csv file named "Catchtraj-cmsy.csv". Yearly values for each simulation are stored across columns. The first column holds the likelihood values (1= accepted, 0 = rejected). The number of rows equals the number of simulations (nsims). Use setwd() before running the function to change the directory where .csv files are stored.
The method of Martell and Froese (2012) is used to produce estimates of MSY where only catch and information on resilience is known.
The Schaefer production model is
B[t+1]<-B[t]+r*B[t]*(1-B[t]/k)-catch[t]
where B is biomass in year t, r is the instrince rate of increase, k is the carrying
capacity and catch is the catch in year t. Biomass is the first year is calculated by
B[1]=k*l0. For sigv>0, the production equation is multiplied by exp(rnorm(1,0,sigv))
if process error is desired.
The maximum sustainable yield (MSY) is calculated as
MSY=r*k/4
Biomass at MSY is calculated as
Bmsy=k/2
Fishing mortality at MSY is calculated as
Fmsy=r/2
Exploitation rate at MSY is calculated as
Umsy=(Fmsy/(Fmsy+M))*(1-exp(-Fmsy-M))
The overfishing limit in last year+1 is calculated as
OFL=B[last year +1]*Umsy
length(year)+1 biomass estimates are made for each run.
If using the R Gui (not Rstudio), run
graphics.off() windows(width=10, height=12,record=TRUE) .SavedPlots <- NULL
before running the catchmsy function to recall plots.
Martell, S. and R. Froese. 2012. A simple method for estimating MSY from catch and resilience. Fish and Fisheries 14:504-514.
# NOT RUN {
# }
# NOT RUN {
data(lingcod)
outpt<-catchmsy(year=lingcod$year,
catch=lingcod$catch,catchCV=NULL,
catargs=list(dist="none",low=0,up=Inf,unit="MT"),
l0=list(low=0.8,up=0.8,step=0),
lt=list(low=0.01,up=0.25,refyr=2002),sigv=0,
k=list(dist="unif",low=4333,up=433300,mean=0,sd=0),
r=list(dist="unif",low=0.015,up=0.1,mean=0,sd=0),
M=list(dist="unif",low=0.18,up=0.18,mean=0.00,sd=0.00),
nsims=30000)
# }
Run the code above in your browser using DataLab