dbsra(year = NULL, catch = NULL, catchCV = NULL,
catargs = list(dist = "none", low = 0, up = Inf, unit = "MT"),
agemat = NULL, k = list(low = 0, up = NULL, tol = 0.01, permax = 1000),
btk = list(dist = "unif", low = 0, up = 1, mean = 0, sd = 0, refyr = NULL),
fmsym = list(dist = "unif", low = 0, up = 1, mean = 0, sd = 0),
bmsyk = list(dist = "unif", low = 0, up = 1, mean = 0, sd = 0),
M = list(dist = "unif", low = 0, up = 1, mean = 0, sd = 0), nsims = 10000,
catchout = 0, grout = 1,
graphs = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14),
grargs = list(lwd = 1, 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 specificed 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.
k (carrying capacity). low and up are
the lower and upper bounds of the minimization routine and tol is the tolerance level
for minimization. A simple equation ((btk)-(b[refyr]/k))^2 is used as
the objective function. R function optimize is used to find k. btk is described
below. permax is the absolute percent difference between the maximum biomass estimate
and k that should not be exceeded in the rule for accepting a run (see details).
refyr).
dist is the statistical distribution name from which to sample btk.
low and up are the lower and upper bounds of btk
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 constant. refyr is the selected terminal year and can range from the first year
to the year after the last year of catch (t+1).
dist is the statistical distribution name from which
to sample Fmsy/M. low and up are the lower and upper bounds of Fmsy/M in
the selected distribution. mean and sd are the mean and standard deviation for selected
distributions. Valid distributions are the same as in btk. If dist= "none", the mean is used as
a fixed constant.
dist is the statistical distribution name from which
to sample Bmsy/k. low and up are the lower and upper bounds of Bmsy/k in
the selected distribution. mean and sd are the mean and standard deviation for selected
distributions. Valid distributions are the same as in btk. If dist= "none", the mean is used as
a fixed constant.
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 btk. If dist= "none", the mean is used as
a fixed constant. M is used to determine exploitation rate (Umsy) at MSY.
setwd before running function to direct .tif files
to a specific directory. Each name of each file is automatically determined.k values, 3 = histogram of plausible Bmsy values,
4 = histogram of plausible MSY values, 5 = histogram of plausible Fmsy values, 6 = histogram of Umsy values,
7 = histogram of plausible Cmsy , 8 = histogram of Bmsy from plausible M, 9 = histogram of plausible Bt/k values,
10 = histogram of plausible Fmsy/M values, 11 = histogram of plausible Bmsy/k values and 12 = histogram of
plausible biomasses in termyr, 13 = line plots of accepted and rejected biomass trajectores
with median and 2.5th and 97.5th percentiles (in red) and 14 = stacked histograms of accepted and
rejected values for each input parameter and resulting estimates and if grout=2,
.tif files are saved with "AR" suffix. Any combination of graphs can be
selected within c(). Default is all.
lwd is the line width for graph = 1 and 13,
nclasses is the nclass argument for the histogram plots (graphs 2-12,14), 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 median value;
llty and llwd are the line type and line width of the 2.5
ulwd are the line type and line width of the 97.5
tiff help file in R.dlproj.dlproj.catchmsy object for use in function dlproj.nsims). This file is loaded to plot
graph 13 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-dbsra.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 delay-difference model is used to propogate biomass: B[t+1]<-B[t]+P[Bt-a]-C[t]
where B[t] is biomass in year t, P[Bt-a] is latent annual production based
on parental biomass agemat years earlier and C[t] is the catch in year
t. Biomass in the first year is assumed equal to k.
If Bmsy/k>=0.5, then P[t] is calculated as
P[t]<-g*MSY*(B[t-agemat]/k)-g*MSY*(B[t-agemat]/k)^n
where MSY is k*Bmsy/k*Umsy, n is solved iteratively using the equation, Bmsy/k=n^(1/(1-n)), and g is (n^(n/(n-1)))/(n-1). Fmsy is calculated as Fmsy=Fmsy/M*M and Umsy is calculated as (Fmsy/(Fmsy+M))*(1-exp(-Fmsy-M)).
If Bsmy/k < 0.5, Bjoin is calculated based on linear rules:
If Bmsy/k<0.3, bjoin="0.5*Bmsy/k*k" if="" bmsy="" k="">0.3 and Bmsy/k<0.5, bjoin="(0.75*Bmsy/k-0.075)*k" if="" any="" b[t-a] P[Bt-agematt where c =(1-n)*g*MSY*Bjoin^(n-2)*K^(-n)
Biomass at MSY is calculated as: Bmsy=(Bmsy/k)*k The overfishing limit (OFL) is Umsy*B[termyr]. The rule for accepting a run is:
if(min(B)>0 && max(B)<=k &&="" <="" p=""> (objective function minimum<=tol^2) &&="" abs(((max(b)-k)="" k)*100)<="permax If using the R Gui (not Rstudio), run graphics.off()
windows(width=10, height=12,record=TRUE)
.SavedPlots <- NULL before running the dbsra function to recall plots.length(year)+1 biomass estimates are made for each run.
Dick, E. J. and A. D. MacCall. 2011. Depletion-based stock reduction analysis: a catch-based method for determining sustainable yield for data-poor fish stocks. Fisheries Research 110: 331-341.
catchmsy dlproj ## Not run:
# data(cowcod)
# dbsra(year =cowcod$year, catch = cowcod$catch, catchCV = NULL,
# catargs = list(dist="none",low=0,up=Inf,unit="MT"),
# agemat=11, k = list(low=100,up=15000,tol=0.01,permax=1000),
# btk = list(dist="beta",low=0.01,up=0.99,mean=0.1,sd=0.1,refyr=2009),
# fmsym = list(dist="lnorm",low=0.1,up=2,mean=-0.223,sd=0.2),
# bmsyk = list(dist="beta",low=0.05,up=0.95,mean=0.4,sd=0.05),
# M = list(dist="lnorm",low=0.001,up=1,mean=-2.90,sd=0.4),
# nsims = 10000)
# ## End(Not run)
Run the code above in your browser using DataLab