# Depletion-Based Stock Reduction Analysis This function estimates MSY from catch following Dick and MAcCall (2011). ```r 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]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. ## Returns - **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. ## 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. ## Author(s) Gary A. Nelson, Massachusetts Division of Marine Fisheries gary.nelson@mass.gov ## 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. ## See Also `catchmsy` `dlproj` ## Examples ```r ## 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) ```