Depletion-Based Stock Reduction Analysis (DB-SRA) is an assessment
method for data-limited fisheries in cases where approximate catches are
known from the start of exploitation, but there are no indices of
abundance or informative composition data. DB-SRA requires an
approximate natural mortality rate (*M*) and age at maturity. A
production function is specified based on the relative location of
maximum productivity to carrying capacity, and the relationship of
*F*_{MSY} to *M*. Given a probability distribution
for the depletion level (i.e. *B*/*B*_{0}) near
the end of the time series, the method derives a probability
distribution for unfished biomass *B*_{0}, and
distributions for management reference points. The method was published
by Dick and MacCall
*(Fisheries
Research 110:331-341)*.
*Arnold
and Heppell* have evaluated DB-SRA using retrospective
analysis.

DB-SRA is best suited for situations in which there is a time-series of historical catches, ideally from the start of the fishery. Also, DB-SRA requires priors on multiple parameters. Some of these can be inferred from the literature but the prior on recent depletion (biomass relative to carrying capacity) drives the results to a substantial extent. Senstivity to this prior should be ideally be examined.

This package is run from R and installed from CRAN. R can be installed from https://cran.r-project.org/.

The recommended way to run R is using RStudio. After you have installed R, then install RStudio. It can be found at https://www.rstudio.com/products/rstudio/download/.

There are several implementations of DB-SRA. The supported version of
DB-SRA is available from the
*fishmethods*
package within R. Another version can be found in the
*Data-Limited Methods
Toolkit*. To use DB-SRA, you first need to install the
*fishmethods* package.

A reference manual for *fishmethods* and therefore
*dbsra* can be found at “https://cran.r-project.org/web/packages/fishmethods/fishmethods.pdf”.

If you are using this code for the first time, then you need to change eval=FALSE to eval=TRUE.

```
Check_Install.packages <- function(pkg){
new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if (length(new.pkg))
install.packages(new.pkg,dependencies = TRUE) else print(paste0("'",pkg,"' has been installed already!"))
sapply(pkg, require, character.only = TRUE)
}
Check_Install.packages("fishmethods")
```

`## [1] "'fishmethods' has been installed already!"`

`## Loading required package: fishmethods`

```
## fishmethods
## TRUE
```

NB Once you have installed *fishmethods*, you don’t need to run
the install.packages code chunk again.

`library(fishmethods)`

```
Example <- read.csv("Catches.csv")
# Read in the simulated annual catch data set
Example <- Example[-c(1:11),c(1,3)]
# We have removed the first 11 rows of data as these are all prefishery and zero, also removed the fleet column as only one fleet so not needed
names(Example) <- c("year","catch")
# this makes sure the headings are correct for later use
print(head(Example))
```

```
## year catch
## 12 1971 98.9982
## 13 1972 199.1610
## 14 1973 991.5380
## 15 1974 988.8430
## 16 1975 1978.3500
## 17 1976 1999.4000
```

Our data set was created using Stock Synthesis based on an age at 50%
maturity of roughly 10-11 years and a rate of natural mortality of
0.2yr^{-1} for females and 0.25yr^{-1} for males. The
catch data are in metric tons.

DB-SRA is run using the function *dbsra*. There are many
inputs to this function. These can be divided into general (required)
inputs, those related to the priors and those related to the other
aspects of the assessment.

*general (required) inputs:***year**: a vector of the years for which historical catches are available.**catch**: the vector of the historical catches corresponding to the vector*year*.**catchCV**: a vector of CVs to capture uncertainty in catches (*catchCV=NULL*implies that catches are known exactly and therefore*catchCV*is not required).

*priors and constraints:***catargs**: specifies how catch uncertainty is quantified using a list*list(dist,low,up,unit)*. The options for the probability distribution for the historical catches are:**norm**for normal catches,**lnorm**for log-normally catches, and**unif**for uniformally distributed catches. The**unit**argument specifies the unit for reporting. For example*catargs=list(dist=“unif”,low=catch*would specify that the historical catches are uniform between 90% and 110% of the values specified in*0.9,up=catch*1.1,unit=“MT”)**catch**argument.**maxn**: (default**25**) the model computes the Pella-Tomlinson shape parameter and this parameter places an upper bound on it.**k**: specifies the bounds on carrying capacity using a list*list(low,up,tol,permax)*. The computed carrying capacities are constrained to be between the*up*and*low*arguments. Only simulations that match the generated and computed final depletion within the*tol*argument (default**0.01**) are included in the final sample.**b1k**: specifies the prior for the relative depletion level (**B1K**) in the first year using a list*list(dist,low,up,mean,sd)*. The options for the selected distribution*dist*argument are: (a)*none*(B1K equals the mean), (b)*unif*(B1K is drawn between the*low*and*up*parameters), (c)*norm*,*lnorm*,*gamma*, and*beta*(the distribution for B1K is defined by the parameters*mean*and*sd*within bounds defined by*low*and*up*.**btk**: specifies the prior for the relative depletion level in a reference year using a list*list(dist,low,up,mean,sd,refyr)*. The parameters are specified as for**b1k**.**M**: specifies the prior for the rate for natural mortality (*M*) using a list*list(dist,low,up,mean,sd)*. The parameters are specified as for**b1k**.**fmsym**: specifies the prior for the rate for*F*_{MSY}/*M*using a list*list(dist,low,up,mean,sd)*. The parameters are specified as for**b1k**.**bmsyK**: specifies the prior for the rate for*B*_{MSY}/*K*using a list*list(dist,low,up,mean,sd)*. The parameters are specified as for**b1k**.

*other inputs***agemat**: specifies the age (in years) that animals enter into the reproductive biomass.**nsims**: number of Monte Carlo samples. This number should ideally be 10,000 or larger.**grout**: send the output the console (1) or to both console and a TIF files (2). The file names are automatically generated if a TIF file option is selected.**graphs**: the graphs that can be produced depend on this argument:- 1: line plot of observed catch versus year,
- 2: histogram of
*K*values - 3: histogram of
*B*_{MSY}values - 4: histogram of
*MSY*values - 5: histogram of
*F*_{MSY}values, - 6: histogram of
*U*_{MSY}(exploitation rate at*MSY*) values, - 7: histogram of overfishing limit (OFL) in last year+1 values
- 8: histogram of
*M*values - 9: histogram of
*B*_{t}/*K*values - 10: histogram of
*F*_{MSY}/*M*values - 11: histogram of
*B*_{MSY}/*K*values - 12: histogram in biomass in the reference year
- 13: line plots of accepted and rejected biomass trajectories with median and 2.5th and 97.5th percentiles (in red).
- 14: stacked histograms of accepted and rejected values for each input parameter and resulting estimates

The model is run by calling *dbsra*. The results are the
posteriors for the parameters (*K*, *B*_{MSY},
*MSY*, *F*_{MSY}, *U*_{MSY}, OFL,
*B*_{MSY}, *B*_{t}/*K*,
*F*_{MSY}/*M*,
*B*_{MSY}/*K*, biomass in the last year). The OFL
is the catch for the first projection year based on applying
*F*_{MSY} to current biomass. Note that the examples here
are based on 100 simulations. This is insufficient for an actual
assessment.

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). Thus an error “None of the runs had a likelihood equal to 1” means that none of the runs were accepted.

There are 14 graphs. The default is to produce all of them. A combination of graphs can be selected using the c() statement e.g. c(1,13) or c(1:13).

Below we run the model close to the example settings. This is a reasonably uninformed model (as would often be the case) and shows how DBSRA works in accepting and rejecting various parameter permutations.

```
par(mfrow=c(2,2))
dbsraOut <- dbsra(year = Example$year,
catch = Example$catch,
catchCV = NULL,
catargs = list(dist="none",low=0,up=Inf,unit="MT"),
agemat=4,
k = list(low=100,up=500000,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.3,sd=0.1,
refyr=max(Example$year)+1),
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=-1.609,sd=0.4),
graph=c(1:14),
nsims = 100,grout=2)
```

The output from the function includes a summary of the specified
priors (object *Initial*), the selected parameters (object
*Parameters*) and the values for the management reference points
(object *Estimates*).

`print(dbsraOut$Initial)`

```
## Distr Lower Upper Mean SD
## Fmsy/M lnorm 0.1 2 -0.223 0.2
## Br/K beta 0.01 0.99 0.3 0.1
## Bmsy/K beta 0.05 0.95 0.4 0.05
## M lnorm 0.001 1 -1.609 0.4
## refyr <NA> 2002 <NA> <NA> <NA>
## B1/K none 0.01 0.99 1 0.1
```

`print(dbsraOut$Parameters)`

```
## Mean (ll=1) Median (ll=1) 2.5% (ll=1) 97.5% (ll=1) min (ll=1) max (ll=1)
## Fmsy/M 0.752 0.7432 0.5458 1.0090 0.53002146 1.0231011
## Bt/K 0.343 0.3434 0.2036 0.5303 0.19863864 0.5644684
## Bmsy/K 0.388 0.3920 0.3037 0.4369 0.25346124 0.4758752
## M 0.168 0.1624 0.0861 0.3477 0.07081843 0.4886967
## B1/K 1.000 1.0000 1.0000 1.0000 1.00000000 1.0000000
```

`print(dbsraOut$Estimates)`

```
## Mean (ll=1) Median (ll=1) 2.5% (ll=1) 97.5% (ll=1) min (ll=1)
## MSY 2782.432 2727.4521 1840.9497 3774.0821 1.749838e+03
## Bmsy 27703.817 27374.7302 17700.6517 36471.7207 1.658075e+04
## Fmsy 0.125 0.1163 0.0590 0.3060 5.447684e-02
## Umsy 0.106 0.1019 0.0554 0.2115 5.040933e-02
## OFL 2567.882 2363.7366 1142.6355 4619.6790 1.097507e+03
## Brefyr 24515.234 22959.4675 14181.2751 40499.8337 1.321638e+04
## K 71635.047 70786.7114 46861.5149 88343.7228 4.633733e+04
## max (ll=1)
## MSY 3.830901e+03
## Bmsy 3.824943e+04
## Fmsy 3.155663e-01
## Umsy 2.310451e-01
## OFL 5.110425e+03
## Brefyr 4.336656e+04
## K 9.785325e+04
```

The biomass trajectory is saved to a file *Biotraj-dbsra.csv*
in the working directory. The first column indicates whether a
simulation (in rows) is accepted (indicated as a “1”) or rejected (“0”),
which should be ignored. These results could therefore be plotted in
your own way. We use the example plot of the biomass trajectories here
as this is quite a useful plot and easily modified.

```
par(mfrow=c(1,1))
# Read in the biomass output, remove the rejected trajectories, and add headers
Years <- c(Example$year,Example$year[length(Example$year)]+1)
BioT <- read.csv("Biotraj-dbsra.csv",head=F)
BioT <- BioT[BioT[,1]==1,-1]
colnames(BioT) <- Years
Nyear <- length(Years)
BiomassSummary <- matrix(0,ncol=5,nrow=Nyear)
for (Iyear in 1:Nyear)
BiomassSummary[Iyear,] <- quantile(BioT[,Iyear],
prob=c(0.05,0.25,0.5,0.75,0.95))
ymax <- max(BiomassSummary)
plot(Years,BiomassSummary[,3],xlab="Year",
ylab="Biomass",ylim=c(0,ymax*1.1),
yaxs="i",type="l",xaxs="i")
polygon(c(Years,rev(Years)),c(BiomassSummary[,1],
rev(BiomassSummary[,5])),col="gray50")
polygon(c(Years,rev(Years)),c(BiomassSummary[,2],
rev(BiomassSummary[,4])),col="gray90")
lines(Years,BiomassSummary[,3],col='red',lwd=2)
```

Stock status (biomass relative to carrying capacity) can be evaluated
from the *B*_{t}/*K* output using
*dbsraOut$Parameters*. It is possible to conduct projections
based on the results of DB-SRA using the *dlproj* function. The
key inputs for this function are:

**dlobj**: the output object from*dbsra*.**projyears**: the number of projection years**projtype**: 0 - median MSY; 1 - mean MSY; 2 - user-specified catch**projcatch**: the vector of projected catches

`ProjRes <- dlproj(dlobj=dbsraOut,projyears=5,projtype=2,projcatch=c(1400,1400,1000,1000,100))`

`print(summary(t(ProjRes$ProjBio)))`

```
## 2002 2003 2004 2005
## Min. :13216 Min. :13670 Min. :14231 Min. :15270
## 1st Qu.:17970 1st Qu.:19488 1st Qu.:20106 1st Qu.:21204
## Median :22959 Median :24354 Median :25869 Median :27562
## Mean :24515 Mean :25659 Mean :26862 Mean :28503
## 3rd Qu.:28043 3rd Qu.:29123 3rd Qu.:30353 3rd Qu.:32048
## Max. :43367 Max. :45414 Max. :47427 Max. :49783
## 2006 2007
## Min. :16335 Min. :18318
## 1st Qu.:22390 1st Qu.:24680
## Median :29650 Median :32708
## Mean :30157 Mean :32714
## 3rd Qu.:34182 3rd Qu.:37345
## Max. :52076 Max. :55192
```

It is possible to apply control rules other than the Overfishing
Limit (OFL) control rule to the results from DB-SRA. For example, the
following code computes a posterior for the catch limit based on a
control rule that sets the exploitation rate to *E*_{MSY}
when stock size is above 35% of *K* and reduces it linearly to
zero at 20% of *K*. Note that this code does not assume that
*refyr* is the last year+1.

```
# Extract the biomass for last year+1 and the parameters that are selected
CurrentBiomass <- BioT[,Nyear]
Parameters <- dbsraOut$Values[dbsraOut$Values[,1]==1,-1]
Nvectors <- length(CurrentBiomass)
# Find the multiplier for the catch limits
Multiplier <- rep(1,Nvectors)
for (Isim in 1:Nvectors)
{
if (CurrentBiomass[Isim] < 0.2*Parameters$K[Isim])
{
Multiplier[Isim] <- 0
}
else
{
if (CurrentBiomass[Isim] < 0.35*Parameters$K[Isim])
Multiplier[Isim] <- (CurrentBiomass[Isim]/Parameters$K[Isim]-0.2)/(0.35-0.2)
}
}
CatchLimits <- CurrentBiomass*Parameters$Umsy*Multiplier
par(mfrow=c(2,2))
hist(CurrentBiomass/Parameters$K[Isim],xlab="Depletion",main="")
hist(Multiplier,xlab="Multiplier",main="")
hist(CatchLimits,xlab="Catch Limit",main="")
```

Out of interest, how well do we do if we use (close to) the actual
values from the simulation? The run below sets the parameters close to
the actual values, although with a uniform distribution around the
values, not as a beta or lognormal distribution, which would make the
model very informative. The exception is the value for *K* for
which we use an uninformative prior.

```
par(mfrow=c(2,2))
dbsraOutNew <- dbsra(year = Example$year,
catch = Example$catch,
catchCV = NULL,
catargs = list(dist="none",low=0,up=Inf,unit="MT"),
agemat=10,
k = list(dist="unif",low=20000,up=250000,mean=21000, sd=0.05),
b1k = list(dist="unif",low=0.95,up=0.99, mean=0.98, sd=0.05),
btk = list(dist="unif",low=0.20,up=0.30,mean=0.27,sd=0.04,
refyr=max(Example$year)+1),
fmsym = list(dist="unif",low=0.65,up=0.85, mean=0.8, sd=0.05),
bmsyk = list(dist="unif",low=0.2,up=0.3,mean=0.261,sd=0.05),
M = list(dist="unif",low=0.2,up=0.3, mean=0.25,sd=0.05),
graph=c(1:14),
nsims = 100,grout=1)
```

```
## Warning in dbsra(year = Example$year, catch = Example$catch, catchCV = NULL, :
## <3 runs were rejected!
```