In natureanalytics-ca/BycatchEstimator: Generic Model-Based Bycatch Estimation Procedure

Although the model was originally developed for bycatch estimation, it can also produce a CPUE index of abundance. For this exercise we will turn off bycatch estimation. We will use the LLSIM blue marlin data. The following models will show how to:

1. Compare model types.

2. Add random effects

3. Add more variables.

# devtools::install_github("ebabcock/BycatchEstimator")  #Make sure library is installed
library(BycatchEstimator)
library(MuMIn)
library(tidyverse)

Model fitting and model selection can be slow. So, for the purposes of this demonstration, we will limit the data to a 15 year period so it will run faster.

obsdatTest<-filter(LLSIM_BUM_Example_observer,Year %in% 2001:2015)
logdatTest<-filter(LLSIM_BUM_Example_logbook,Year %in% 2001:2015)

1. Compare model types.

To set up the model for bycatch estimation only, use the following inputs in bycatchSetup. The formula in indexModel gives the variables to be kept separate in the index. This will always include Year. If you need separate indices by fleet or area, include them here. Note that the best model in between simpleModel and complexModel will be used to generate the index, and this may have different variables. Generally, the index will be calculated at the mean value of any numerical predictor variables, and at the most common level of any factors. If a variable, such as area, is included in indexModel but not in the selected best model according to the information criteria, the areas will all have the same index.

setupObj<-bycatchSetup(
 modelTry = c("TMBdelta-Lognormal","TMBdelta-Gamma","TMBnbinom2","TMBtweedie"), #model types to compare
 obsdat = obsdatTest,
  logdat = NULL,
  yearVar = "Year",
  obsEffort = "hooks",
  logEffort = "hooks",
  factorNames = c("Year","season","area"),
  EstimateIndex = TRUE,
  EstimateBycatch = FALSE,
  logNum=NULL,
  sampleUnit = "trips",
  complexModel = formula(y~Year+season+area+fleet +area:season), #upper range of models to compare with information criteria,
  simpleModel = formula(y~Year), #lower range of models to compare
  indexModel = formula(y~Year+area), #Variables to include in index.
  designMethods ="None",
  baseDir = getwd(),
  runName = "LLSIMBUMCPUE1",
  runDescription = "LLSIm BUM CPUE",
  common = "Blue marlin",
  sp = "Makaira nigricans",
  obsCatch = "BUM",
  catchUnit = "number",
  catchType = "catch"
)

bycatchFit doesn't need any changes.

bycatchFit(
  setupObj,
  selectCriteria = "BIC",
  DoCrossValidation = TRUE,
  DredgeCrossValidation = FALSE,
  ResidualTest = FALSE,
  CIval = 0.05,
  VarCalc = "DeltaMethod",
  useParallel = TRUE,  #Make this FALSE if parallel processing doesn't work. 
  nSims = 100,
  baseDir = getwd(),
)

Now, looking at the output .pdf, there will be no total bycatch estimates, but there will be abundance indices plotted, plus and minus a standard error. The data are also in a .csv file.

nbinom2Index<-read.csv("Output LLSIMBUMCPUE1/Blue marlin catch/Blue marlincatchTMBnbinom2Index.csv")
ggplot(nbinom2Index,aes(x=Year,y=Index,ymin=ymin,ymax=ymax))+
  geom_line()+
  geom_ribbon(alpha=0.3)+
  facet_wrap(area~.)

We can see that all the models have similar diagnostics, but imply somewhat different trends.

1. Add random effects

CPUE standardization models commonly include random effects, for example to handle interactions between Year and other variables. It is also worthwhile to include fishing vessel as a random effect when using sets as a sample unit, to deal with the clumping in the data. Let's try Year:area as a random effect. Of course, random effects can also be used for bycatch estimation, but this is less commonly done. If random effects are included, they will be automatically included in all models. randomEffects is the list of random effects to include in the model. For delta models, randomEffects2 can be used to give a different set of randomEffects for the two model components.

setupObj<-bycatchSetup(
 modelTry = c("TMBdelta-Lognormal","TMBnbinom1","TMBdelta-Gamma","TMBnbinom2","TMBtweedie"), #model types to compare
 obsdat = obsdatTest,
  logdat = NULL,
  yearVar = "Year",
  obsEffort = "hooks",
  logEffort = "hooks",
  factorNames = c("Year","season","area"),
  randomEffects = c("Year:area"),
  randomEffects2 = c("Year:area"),
  EstimateIndex = TRUE,
  EstimateBycatch = FALSE,
  logNum=NULL,
  sampleUnit = "trips",
  complexModel = formula(y~Year+season+area+fleet +area:season), #upper range of models to compare with information criteria, fixed effects only. 
  simpleModel = formula(y~Year), #lower range of models to compare, fixed effects only
  indexModel = formula(y~Year+area), #Variables to include in index.
  designMethods ="None",
  baseDir = getwd(),
  runName = "LLSIMBUMCPUE2",
  runDescription = "LLSIm BUM CPUE with Year:area",
  common = "Blue marlin",
  sp = "Makaira nigricans",
  obsCatch = "BUM",
  catchUnit = "number",
  catchType = "catch"
)

bycatchFit(
  setupObj,
  selectCriteria = "BIC",
  DoCrossValidation = TRUE,
  DredgeCrossValidation = FALSE,
  ResidualTest = FALSE,
  CIval = 0.05,
  VarCalc = "DeltaMethod",
  useParallel = TRUE,  #Make this FALSE if parallel processing doesn't work.
  nSims = 100,
  baseDir = getwd(),
)

1. Add more variables.

For CPUE standardization, it may be desirable to add more variables, such as gear variables and environmental variables. Numerical variables will be set to their mean for the prediction. Here we add two numerical variables (hbf and habBUM) and more factors (light).

setupObj<-bycatchSetup(
 modelTry = c("TMBdelta-Lognormal","TMBnbinom1","TMBdelta-Gamma","TMBnbinom2","TMBtweedie"), #model types to compare
 obsdat = obsdatTest,
  logdat = NULL,
  yearVar = "Year",
  obsEffort = "hooks",
  logEffort = "hooks",
  factorNames = c("Year","season","area","light"),
  randomEffects = NULL,
  randomEffects2 = NULL,
  EstimateIndex = TRUE,
  EstimateBycatch = FALSE,
  logNum=NULL,
  sampleUnit = "trips",
  complexModel = formula(y~Year+season+area+fleet +light+hbf+habBUM), #upper range of models to compare with information criteria, fixed effects only. 
  simpleModel = formula(y~Year), #lower range of models to compare, fixed effects only
  indexModel = formula(y~Year+area), #Variables to include in index.
  designMethods ="None",
  baseDir = getwd(),
  runName = "LLSIMBUMCPUE3",
  runDescription = "LLSIm BUM CPUE more variables",
  common = "Blue marlin",
  sp = "Makaira nigricans",
  obsCatch = "BUM",
  catchUnit = "number",
  catchType = "catch"
)

bycatchFit(
  setupObj,
  selectCriteria = "BIC",
  DoCrossValidation = TRUE,
  DredgeCrossValidation = FALSE,
  ResidualTest = FALSE,
  CIval = 0.05,
  VarCalc = "DeltaMethod",
  useParallel = TRUE,  #Make this FALSE if parallel processing doesn't work. 
  nSims = 100,
  baseDir = getwd(),
)

In this case the trends differ somewhat, perhaps in part because the IC selected different predictor variables for different models.

natureanalytics-ca/BycatchEstimator documentation built on June 2, 2025, 12:38 p.m.