Installation

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

## Warning: package 'fishmethods' was built under R version 4.0.3

## fishmethods 
##        TRUE

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

library(fishmethods)

Reading in the data

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.

Description of the model settings

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=catch0.9,up=catch1.1,unit=“MT”) would specify that the historical catches are uniform between 90% and 110% of the values specified in 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

Running the model

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.786        0.7496      0.5304       1.0891 0.51134260  1.2452057
## Bt/K         0.354        0.3459      0.1815       0.5602 0.12823316  0.6728200
## Bmsy/K       0.386        0.3877      0.2901       0.4633 0.26068355  0.4785907
## M            0.185        0.1675      0.0882       0.3971 0.08245986  0.4855176
## 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       2954.447     2948.7304   2046.6456    3852.3308 2.022274e+03
## Bmsy     26067.545    26056.3429  16769.7159   36269.9088 1.555801e+04
## Fmsy         0.143        0.1244      0.0793       0.3209 6.744497e-02
## Umsy         0.119        0.1079      0.0722       0.2257 6.248115e-02
## OFL       2816.503     2714.0153   1165.1159    6119.5642 7.248675e+02
## Brefyr   24055.045    22775.4699  11096.4857   43115.6530 9.001708e+03
## K        67944.077    66470.1197  45886.4352   93217.1327 4.564667e+04
##          max (ll=1)
## MSY    4.026953e+03
## Bmsy   3.730537e+04
## Fmsy   3.408840e-01
## Umsy   2.413262e-01
## OFL    6.605754e+03
## Brefyr 4.608397e+04
## K      1.050514e+05

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)

Interpreting results in a management context

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.   : 9002   Min.   : 8937   Min.   : 8950   Min.   : 9413  
##  1st Qu.:18637   1st Qu.:19500   1st Qu.:20522   1st Qu.:22365  
##  Median :22775   Median :24217   Median :25756   Median :27954  
##  Mean   :24055   Mean   :25337   Mean   :26689   Mean   :28481  
##  3rd Qu.:28486   3rd Qu.:29875   3rd Qu.:31275   3rd Qu.:33256  
##  Max.   :46084   Max.   :46999   Max.   :48973   Max.   :51530  
##       2006            2007      
##  Min.   : 9863   Min.   :11185  
##  1st Qu.:24625   1st Qu.:27423  
##  Median :29730   Median :32646  
##  Mean   :30282   Mean   :32982  
##  3rd Qu.:35130   3rd Qu.:37851  
##  Max.   :53884   Max.   :56907

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

Comparing to the simulated values

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!

Let’s compare the two sets of parameter values:

Compare <- as.data.frame(cbind(rownames(dbsraOut$Estimates),dbsraOut$Estimates[,"Median (ll=1)"],dbsraOutNew$Estimates[,"Median (ll=1)"]))
names(Compare) <- c("Estimate", "Median Uninformed", "Median Informed")
print(Compare)

##   Estimate Median Uninformed Median Informed
## 1      MSY         2948.7304       3164.0267
## 2     Bmsy        26056.3429      20624.3116
## 3     Fmsy            0.1244           0.189
## 4     Umsy            0.1079          0.1521
## 5      OFL         2714.0153       3171.5216
## 6   Brefyr        22775.4699      21106.1792
## 7        K        66470.1197      85104.8896

Key references

DB-SRA Methodology

Dick, E.J. and A.D. MacCall. 2011. Depletion-Based Stock Reduction Analysis: A catch-based method for determining sustainable yields for data-poor fish stocks. Fisheries Research 110:331-341. https://doi.org/10.1016/j.fishres.2011.05.007.

DB-SRA Evaluation

Arnold, L.M. and S.S. Heppell, 2015. Testing the robustness of data-poor assessment methods to uncertainty in catch and biology: a retrospective approach. ICES Journal of Marine Science 72:243-250. https://doi.org/10.1093/icesjms/fsu077.

Carruthers, T.R., Punt, A.E., Walters, C.J., MacCall, A., McAllister, M.K., Dick, E.J. and J. Cope. 2015. Evaluating methods for setting catch limits in data-limited fisheries. Fisheries Research 153:48-68. https://doi.org/10.1016/j.fishres.2013.12.014.

Wetzel, C.R. and A.E. Punt. 2011. Evaluating the performance of data-moderate and catch-only assessment methods for U.S. west coast groundfish. Fisheries Research 171:170-187. https://doi.org/10.1016/j.fishres.2015.06.005.

DB-SRA Applications

Dick, E.J. and A.D. MacCall. 2010. Estimates of sustainable yield for 50 data-poor stocks in the Pacific Coast groundfish fishery management plan. Technical memorandum. Southwest fisheries Science Centre, Santa Cruz, CA. National Marine Fisheries Service, National Oceanic and Atmospheric Administration of the U.S. Department of Commerce NOAA-TM-NMFS-SWFSC-460.

Cope, J., Dick, E.J., MacCall, A., Monk, M., Soper, B. and C. Wetzel. 2013. Data-moderate stock assessments for brown, China, copper, sharpchin, stripetail and yellowtail rockfishes and English and rex soles in 2013.

DBSRA Installation

17 February, 2021

Depletion-Based Stock Stock Reduction Analysis (DB-SRA)