In nikolaifish/bamExtras: Functions and Data in Support of Stock Assessments with the Beaufort Assessment Model (BAM)

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 7,
  warning = FALSE,
  errors = FALSE,
  messages = FALSE
)

Description

bamExtras is an R package that contains functions and data to support stock assessments of fish stocks in the Southeast US Atlantic using the Beaufort Assessment Model (BAM). It provides support and organization for data and functions used to run stock assessments and diagnostics with BAM. Currently it contains all of the most recent versions of BAM for 14 stocks and can be used to modify and run BAM, if ADMB is installed. It also contains various support and plotting functions

Example usage

Load the package

First, install and load the package.

# Install and load package
#devtools::install_github("nikolaifish/bamExtras")
library(bamExtras)
# Load dependencies
library(stringr)

Data

bamExtras contains BAM input files: dat (data), tpl (template model code), cxx (convert C++ output to R objects), and rdat (R output) files for the most recent models. The dat, tpl, and cxx files are stores as character vectors with each item representing a line of the original file. Other files like "admb2r.cpp" and "cleanup.bat" are also stored in bamExtras for running BAM.

# see all data objects in bamExtras
data(package="bamExtras")$result[,"Item"]

Since different naming conventions are sometimes used in BAM, the function standardize_rdat() can be used to modify names to follow common naming conventions. Common naming conventions are important for some functions in bamExtras and allow for general analysis.

rdat_BS <- standardize_rdat(rdat_BlackSeaBass)
rdat_RP <- standardize_rdat(rdat_RedPorgy)

We can simply use the stored data to plot results from an assessment

par(mfrow=c(3,2),mar=c(2,2,1,1),mgp=c(1,0.2,0),tck=-0.01)
with(rdat_BS$t.series,{
     plot(year,N,type="o")
     plot(year,recruits,type="o")  
     plot(year,SSB,type="o")
     plot(year,F.full,type="o")
     plot(year,SSB.msst,type="o")
     plot(year,F.Fmsy,type="o")
  })

We can plot results from multiple assessments

spcol <- c("BS"="black","RP"="red")

par(mfrow=c(3,2),mar=c(2,2,1,1),mgp=c(1,0.2,0),tck=-0.01)
ylim <- range(c(0,rdat_BS$t.series$N,rdat_RP$t.series$N),na.rm=TRUE)
with(rdat_BS$t.series,plot(year,N,type="o",col=spcol["BS"],ylim=ylim))
with(rdat_RP$t.series,points(year,N,type="o",col=spcol["RP"]))

legend("top",legend=c("Black Sea Bass","Red Porgy"),
       pch=1,col=spcol,lty=1)

ylim <- range(c(0,rdat_BS$t.series$recruits,rdat_RP$t.series$recruits),na.rm=TRUE)
with(rdat_BS$t.series,plot(year,recruits,type="o",col=spcol["BS"],ylim=ylim))
with(rdat_RP$t.series,points(year,recruits,type="o",col=spcol["RP"]))

ylim <- range(c(0,rdat_BS$t.series$SSB,rdat_RP$t.series$SSB),na.rm=TRUE)
with(rdat_BS$t.series,plot(year,SSB,type="o",col=spcol["BS"],ylim=ylim))
with(rdat_RP$t.series,points(year,SSB,type="o",col=spcol["RP"]))

ylim <- range(c(0,rdat_BS$t.series$F.full,rdat_RP$t.series$F.full),na.rm=TRUE)
with(rdat_BS$t.series,plot(year,F.full,type="o",col=spcol["BS"],ylim=ylim))
with(rdat_RP$t.series,points(year,F.full,type="o",col=spcol["RP"]))

ylim <- range(c(0,rdat_BS$t.series$SSB.msst,rdat_RP$t.series$SSB.msst),na.rm=TRUE)
with(rdat_BS$t.series,plot(year,SSB.msst,type="o",col=spcol["BS"],ylim=ylim))
with(rdat_RP$t.series,points(year,SSB.msst,type="o",col=spcol["RP"]))

ylim <- range(c(0,rdat_BS$t.series$F.Fmsy,rdat_RP$t.series$F.Fmsy),na.rm=TRUE)
with(rdat_BS$t.series,plot(year,F.Fmsy,type="o",col="black",ylim=ylim))
with(rdat_RP$t.series,points(year,F.Fmsy,type="o",col=spcol["RP"]))

Functions

One of core functions in this package is bam2R. is a helpful function compares corresponding dat and tpl files from an assessment, and builds an R list L_init of all the init objects initialized in the DATA_SECTION of the tpl. It then combines L_init and the BAM files dat, tpl, and cxx into an output list.

bam_RP <- bam2r(dat_obj=dat_RedPorgy,tpl_obj=tpl_RedPorgy,cxx_obj=cxx_RedPorgy)

Drawing upon data objects stored in the package, bam2R can also be used in a simpler way with the CommonName argument

bam_RP <- bam2r("RedPorgy")

Running BAM

If ADMB is installed, we can then rerun the stock assessment. The code in this section is provided by commented out for the vignette.

# Run a bam model and assign rdat output to object
#rdat_RP <- run_bam(bam=bam_RP,fileName="RP")

And bam2r lets us modify the data inputs. Then we can rerun BAM.

# Modify data input from a BAM model and incorporate it back into the dat file object
# L_init2 <- bam_RP$L_init
# L_init2$set_steep[c(1,5)] <- paste(0.6) # Change steepness
# L_init2$set_steep[4] <- paste(-abs(as.numeric(L_init2$set_steep[4]))) # Fix steepness
# bam_RP2 <- bam2r("RedPorgy",L_init_user=L_init2)
# rdat_RP2 <- run_bam(bam=bam_RP2,fileName="RP2")

Then we can compare the two models

# spcol <- c("RP"="black","RP2"="red")
# 
# par(mfrow=c(3,2),mar=c(2,2,1,1),mgp=c(1,0.2,0),tck=-0.01)
# ylim <- range(c(0,rdat_RP$t.series$N,rdat_RP2$t.series$N),na.rm=TRUE)
# with(rdat_RP$t.series,plot(year,N,type="o",col=spcol["RP"],ylim=ylim))
# with(rdat_RP2$t.series,points(year,N,type="o",col=spcol["RP2"]))
# 
# legend("top",legend=c("Red Porgy","Red Porgy2"),
#        pch=1,col=spcol,lty=1)
# 
# ylim <- range(c(0,rdat_RP$t.series$recruits,rdat_RP2$t.series$recruits),na.rm=TRUE)
# with(rdat_RP$t.series,plot(year,recruits,type="o",col=spcol["RP"],ylim=ylim))
# with(rdat_RP2$t.series,points(year,recruits,type="o",col=spcol["RP2"]))
# 
# ylim <- range(c(0,rdat_RP$t.series$SSB,rdat_RP2$t.series$SSB),na.rm=TRUE)
# with(rdat_RP$t.series,plot(year,SSB,type="o",col=spcol["RP"],ylim=ylim))
# with(rdat_RP2$t.series,points(year,SSB,type="o",col=spcol["RP2"]))
# 
# ylim <- range(c(0,rdat_RP$t.series$F.full,rdat_RP2$t.series$F.full),na.rm=TRUE)
# with(rdat_RP$t.series,plot(year,F.full,type="o",col=spcol["RP"],ylim=ylim))
# with(rdat_RP2$t.series,points(year,F.full,type="o",col=spcol["RP2"]))
# 
# ylim <- range(c(0,rdat_RP$t.series$SSB.msst,rdat_RP2$t.series$SSB.msst),na.rm=TRUE)
# with(rdat_RP$t.series,plot(year,SSB.msst,type="o",col=spcol["RP"],ylim=ylim))
# with(rdat_RP2$t.series,points(year,SSB.msst,type="o",col=spcol["RP2"]))
# 
# ylim <- range(c(0,rdat_RP$t.series$F.Fmsy,rdat_RP2$t.series$F.Fmsy),na.rm=TRUE)
# with(rdat_RP$t.series,plot(year,F.Fmsy,type="o",col="black",ylim=ylim))
# with(rdat_RP2$t.series,points(year,F.Fmsy,type="o",col=spcol["RP2"]))

Conduct a catch curve analysis

# read in black ses bass age comps used in the recent stock assessment
rdat <- rdat_BlackSeaBass
comp <- rdat$comp.mats # List of composition matrices
acomp <- comp[grepl("^acomp.*.ob$",names(comp))] # list of observed age composition matrices
comp_i <- acomp$acomp.Mcvt.ob
cc_i <- catch_curve(comp_i)
# Plot comps and add catch curves (note the log-tranformation of the composition data)
comp_plot(log(comp_i),cc=cc_i,fillComp = FALSE,ylab= "log(proportion)",xlab="age",title="black sea bass commercial handline catch curves")
# Plot Z estimates over time for multiple sets of age compositions
par(mfrow=c(1,1),mgp=c(1,.2,0),mar=c(2,2,1,1),oma=c(0,0,0,0),tck=0.01)
M <- rdat$parms$M.msst
cc <- lapply(acomp,catch_curve)
abb <- gsub("^(acomp.)|(.ob)$","",names(cc)) # abbreviations of age composition fleets
cols <- rainbow(length(cc))
xlim <- range(unlist(lapply(cc,function(x){as.numeric(x$year)})),na.rm=TRUE)
ylim <- range(c(0,unlist(lapply(cc,function(x){-as.numeric(x$slope)}))),na.rm=TRUE)
for(i in 1:length(cc)){
  xi <- cc[[i]]
  year <- as.numeric(xi$year)
  Z <- -xi$slope
  if(i==1){
    plot(year,Z,xlim=xlim,ylim=ylim,type="o",col=cols[i])
  }else{
    points(year,Z,type="o",col=cols[i])
  }
}
# Add reference line for natural mortality
abline(h=M,lty=2,lwd=2)
text(par("usr")[1]+0.05*diff(par("usr")[1:2]),M,labels="M",pos=3)
legend("bottomleft",legend=abb,col=cols,pch=1,lty=1,bty="n")

nikolaifish/bamExtras documentation built on July 21, 2023, 8:26 a.m.