This function estimates MSY from catch following Dick and MAcCall (2011).
dbsra(year = NULL, catch = NULL, catchCV = NULL,
catargs = list(dist = "none", low = 0, up = Inf, unit = "MT"),
agemat = NULL, maxn=25, k = list(low = 0, up = NULL, tol = 0.01, permax = 1000),
b1k = list(dist = "unif", low = 0, up = 1, mean = 0, sd = 0),
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, 15),
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))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 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.
median age at entry to the reproductive biomass.
the maximum limit of the Pella-Tomlinson shape parameter that should not be exceeded in the rule for accepting a run.
list arguments for estimation of 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).
list arguments for B1/K, the relative depletive level in the first year.
dist is the statistical distribution name from which to sample b1k.
low and up are the lower and upper bounds of b1k
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.
list arguments for Bt/K, the relative depletive level in a specific reference year (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).
list arguments for Fmsy/M. 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.
list arguments for Bmsy/k. 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.
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 btk. If dist= "none", the mean is used as
a fixed constant. M is used to determine exploitation rate (Umsy) at MSY.
number of Monte Carlos samples.
if catch is resampled, output the time series from every MC sample to a .csv file. 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 = histogram of plausible (accepted) 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.
list control arguments for plotting functions. 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.
list control arguments for plotting the median and 2.5
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 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
list arguments for the .TIF graph files. See tiff help file in R.
dataframe containing the descriptive statistics for each explored parameter.
dataframe containing the mean, median, 2.5th and 97.5 of the plausible (accepted: likelihood(ll)=1) parameters.
dataframe containing the mean, median, 2.5th and 97.5 of the management quantities (i.e., MSY, Bmsy, etc.) from the plausible parameters (likelihood=1)
dataframe containing the values of likelihood, k, Bt/k, Bmsy/k, M and associated management quantities for all (likelihood=0 and likelihood=1) random draws.
agemat for use in function dlproj.
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-dbsra.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 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 method of Dick and MAcCall (2011) is used to produce estimates of MSY where only catch and information on resilience and current relative depletion is known.
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]<Bjoin, then the Schaefer model is used to calculated P:
P[Bt-agematt<Bjoin]<-B[t-agemat]*(P(Bjoin)/Bjoin+c(B[t-agemat]-Bjoin))
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].
length(year)+1 biomass estimates are made for each run.
The rule for accepting a run is: if(min(B)>0 && max(B)<=k &&
(objective function minimum<=tol^2) && abs(((max(B)-k)/k)*100)<=permax && n<=maxn
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.
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.
# NOT RUN {
# }
# 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),
b1k = list(dist="none",low=0.01,up=0.99,mean=1,sd=0.1),
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)
# }
Run the code above in your browser using DataLab