dbsra: Depletion-Based Stock Reduction Analysis

In fishmethods: Fishery Science Methods and Models

Depletion-Based Stock Reduction Analysis

Description

This function estimates MSY from catch following Dick and MAcCall (2011).

Usage

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))

Arguments

year	vector containing the time series of numeric year labels.
catch	vector containing the time series of catch data (in weight). Missing values are not allowed.
catchCV	vector containing the time series of coefficients of variation associated with catch if resampling of catch is desired; otherwise, catchCV = NULL.
catargs	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.
agemat	median age at entry to the reproductive biomass.
maxn	the maximum limit of the Pella-Tomlinson shape parameter that should not be exceeded in the rule for accepting a run.
k	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).
b1k	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.
btk	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).
fmsym	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.
bmsyk	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.
M	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.
nsims	number of Monte Carlos samples.
catchout	if catch is resampled, output the time series from every MC sample to a .csv file. 0 = no (default), 1 = yes.
grout	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.
graphs	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.
grargs	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.
pstats	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
grtif	list arguments for the .TIF graph files. See tiff help file in R.

Details

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.

Value

Initial	dataframe containing the descriptive statistics for each explored parameter.
Parameters	dataframe containing the mean, median, 2.5th and 97.5 of the plausible (accepted: likelihood(ll)=1) parameters.
Estimates	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)
Values	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	agemat for use in function dlproj.
end1yr	value of the last year of catch data + 1 for use in function dlproj.
type	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.

Note

The random distribution function was adapted from Nadarajah, S and S. Kotz. 2006. R programs for computing truncated distributions. Journal of Statistical Software 16, code snippet 2.

Author(s)

Gary A. Nelson, Massachusetts Division of Marine Fisheries [e-mail]

References

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.

See Also

catchmsy dlproj

Examples

 ## 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.4,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)

fishmethods documentation built on April 27, 2023, 9:10 a.m.