JABBA requires the installation of R and JAGS and the following R packages that can be directly installed within R

Data

JABBA requires a minimum of two input comma-separated value files (.csv) in the form of catch and abundance indices. The Catch input file contains the time series of year and catch by weight, aggregated across fleets for the entire fishery. Missing catch years or catch values are not allowed. JABBA is formulated to accommodate abundance indices from multiple sources (i.e., fleets) in a single cpue file, which contains all considered abundance indices. The first column of the cpue input is year, which must match the range of years provided in the Catch file. In contrast to the Catch input, missing abundance index values are allowed, such that different abundance indices may correspond to smaller portions of the catch time series. Optionally, an additional se input can be passed onto JABBA, containing standard error estimates associated with the abundance indices on a log scale. The se input is a third file, structurally identical to the cpue input. Alternatively, this feature can be used to apply different weighting to individual abundance indices by assigning varying coefficients of variation (CV) to each time series. If such weighting is implemented, it is advised that the CV chosen for each indexed year approximates the observed standard error on the log scale, such that the data weights are congruent with expectations as to how well the model should fit these data.

JABBA provides provides the option to use a single averaged CPUE index instead of the individual abundance indices (see 2.5.1. State-Space model for averaging of abundance indices in Winker et al, 2018). This feature can be activated by setting meanCPUE = TRUE.

rm(list = ls()) # clear memory

all_data <- read.csv(file = here::here("StockA1_catch_cpueSE.csv"), fileEncoding="UTF-8-BOM") 

#inspect
head(all_data)
##   Year Species Total_catch  cpue_sc   cpue_SE
## 1 1992 StockA1      829.20 1.798236        NA
## 2 1993 StockA1      567.97 1.075917 0.5905036
## 3 1994 StockA1      402.71 2.685084 0.5794874
## 4 1995 StockA1      366.71 1.767176 0.5776929
## 5 1996 StockA1      447.09 1.871632 0.5793944
## 6 1997 StockA1      685.99 2.576474 0.5789248
#which species do we have?
unique(all_data$Species)
## [1] "StockA1"

#StockA1

SETUP AND RUN

Setup model parameters

For surplus production models we basically need carrying capacity K and population growth rate r Since K and r are never known, a common approach is to set priors for K at 3-10 times the maximum catch in the catch time series, whereas r is assumed to be lognormally distributed with mean r of 0.2. But you can set it differently. Let’s also make projections into the future assuming that catch in the next 10 years will stay the same as the last year’s catch, be 50% smaller and 50% larger

# K range will be set 3-10 times of maximum catch 
max(catch$catch)
## [1] 1012.672
Kmin <- as.numeric(3*max(catch$catch))
Kmax <- as.numeric(10*max(catch$catch))

#last year's catch
lastcatch <- round((1*as.numeric(catch[nrow(catch),][2])),0)
#Set the TAC_range
TAC_range <- c(0.5*lastcatch, lastcatch, 1.5*lastcatch)

#when to start projections with the new TAC range? Usually from the year after the last year
TAC_year <- catch$year[nrow(catch)] +1

setup: Schaeffer

Schaefer with CPUE SE

scenario_name = "Schaeffer_both"

input_both_sch <- build_jabba(
  catch = catch,
  cpue = cpue,
  se = se,
  assessment = species,
  scenario = scenario_name,
  model.type = c("Schaefer", "Fox", "Pella", "Pella_m") [1],  # with [1] we select the first option
  add.catch.CV = TRUE,
  catch.cv = 0.1,
  catch.error = c("random", "under") [1],          #"under" mean directional underreporting
  Plim = 0.25,                                    #limit of biomass where recruitment can become impaired Plim = Blim/K
  r.dist = c("lnorm", "range")[1],
  r.prior = c(0.4, 0.2),             #IMPORTANT: What is the expected mean R and lognormal errors? Default is 0.2 and 0.25
  K.dist = c("lnorm", "range")[2],
  K.prior = c(Kmin,Kmax),
  psi.dist = c("lnorm", "beta")[1],
  psi.prior = c(0.6, 0.2),
  b.prior = c(FALSE, 0.3, NA, c("bk", "bbmsy", "ffmsy")[1]),
  BmsyK = 0.4,
  shape.CV = 0.3,
  sets.q = 1:(ncol(cpue) - 1),
  sigma.est = TRUE,
  sets.var = 1:(ncol(cpue) - 1),
  #fixed.obsE = 0.01,
  fixed.obsE = ifelse(is.null(se), 0.1, 0.001),
  sigma.proc = TRUE,
  proc.dev.all = TRUE,
  igamma = c(4, 0.01),
  projection = TRUE,
  TACs = TAC_range, # we set a range in here
  #TACs = c(350,1000), # we set a range in here
  TACint = NULL,
  imp.yr = TAC_year,        #when does TAC start
  pyrs = 10,            #how many years to project     
  P_bound = c(0.02, 1.3),
  sigmaobs_bound = 1,
  sigmaproc_bound = 0.2,
  q_bounds = c(10^-30, 1000),
  K_bounds = c(0.01, 10^10),
  harvest.label = c("Hmsy", "Fmsy")[2],
  catch.metric = "(t)"
)

setup: Fox

Fox with CPUE SE

setup: Pella

Pella with CPUE SE

RUN

We will run Schaefer, Fox and Pella scenarios. They make different assumptions about density dependence. We run them with identical parameters and priors, but we expect to get different posterior estimates, because they assume a slightly different production curve

Run: Schaeffer

# creating a directory to save the output files
dir.create(species) 


#note this will save an .RData object 

run_both_sch <- fit_jabba(
  input_both_sch,
  ni = 30000,
  nt = 5,
  nb = 5000,
  nc = 2,
  init.values = FALSE,
  init.K = NULL,
  init.r = NULL,
  init.q = NULL,
  peels = NULL,
  save.all = FALSE,
  save.trj = TRUE,
  save.jabba = TRUE,  #will save the RData object
  save.csvs = TRUE,
  output.dir = species,  ## make sure you create a folder in your working directory to save outputs
  quickmcmc = TRUE
)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 92
##    Unobserved stochastic nodes: 128
##    Total graph size: 2350
## 
## Initializing model

Run: Fox

Run: Pella

MAKE OUTPUTS

all output plots into directory

You can explore these plots from the working directory

jabba_plots(run_both_sch, output.dir = species)
## 
## ><> jbplot_catch()  <><
## 
## ><> jbplot_catcherror()  <><
## 
## ><> jbplot_cpue() - fits to CPUE <><
## 
## ><> jbplot_cpue() - fits to CPUE <><
## 
## ><> jbplot_mcmc() - mcmc chains  <><
## 
## ><> jbplot_ppist() - prior and posterior distributions  <><
## 
## ><> jbplot_procdev() - Process error diviations on log(biomass)  <><
## 
## ><> jbplot_trj() - B trajectory  <><
## 
## ><> jbplot_trj() - F trajectory  <><
## 
## ><> jbplot_trj() - BBmsy trajectory  <><
## 
## ><> jbplot_trj() - FFmsy trajectory  <><
## 
## ><> jbplot_trj() - BB0 trajectory  <><
## 
## ><> jbplot_spphase() - JABBA Surplus Production Phase Plot  <><
## 
## ><> jbplot_cpue() - fits to CPUE <><
## 
## ><> jbplot_runstest()   <><
## 
## ><> jbplot_summary()
## 
## ><> jbplot_kobe() - Stock Status Plot  <><

## png 
##   2
#jabba_plots(run_both_fox, output.dir = species)
#jabba_plots(run_both_pella, output.dir = species)

summary plot in pdf

If you want to generate your summary plot

#load the data file, but it is always saved as name "jabba". Chose which file you want: Schaeffer, Fox or Pella

 load(file = paste0(species,"/",species,"_Schaeffer_both_jabba.rdata"))
#load(file = paste0(species,"/",species,"_Fox_both_jabba.rdata"))
# load(file = paste0(species,"/",species,"_Pella_both_jabba.rdata"))

 #rename it, to not get confused
 model_out <- jabba

 #make a name for a summary pdf file 
pdf(file = paste0(species,"/",species,"_summary.pdf"), width = 11, height = 12)

# plot various outputs 

par(mfrow=c(3,2),mar = c(5, 5, 4, 3))
jbplot_trj(model_out,type="BBmsy",add=T)
## 
## ><> jbplot_trj() - BBmsy trajectory  <><
jbplot_trj(model_out,type="F",add=T)
## 
## ><> jbplot_trj() - F trajectory  <><
jbplot_catch(model_out, add = T)
## 
## ><> jbplot_catch()  <><
jbplot_spphase(model_out,add=T)
## 
## ><> jbplot_spphase() - JABBA Surplus Production Phase Plot  <><
jbplot_kobe(model_out, add = T)
## 
## ><> jbplot_kobe() - Stock Status Plot  <><
dev.off()
## png 
##   2