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))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.
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).
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).
signv = 0 for no
process error.
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.
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.
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.
catch, save catch trajectories from each Monte Carlos simulation
- 0 = No (default), 1 = Yes.setwd before running function to direct .tif files
to a specific directory. Each name of each file is automatically determined.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.
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.
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
tiff help file in R.dlproj.catchmsy object for use in function dlproj.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 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.
dbsra dlproj ## 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)
# ## End(Not run)
Run the code above in your browser using DataLab