catch_curve: Fit catch curves to age composition data

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

Fit catch curves to age composition data

Description

Simple geometric mean function using only is.finite values of log(x).

Usage

catch_curve(
  M,
  age_min = NULL,
  age_max = NULL,
  nobs_min = 3,
  age_min_method = "peak_overall"
)

Arguments

M	matrix of proportions at age where rows are years (catch years or cohort years) and columns are ages. Rows should sum to 1.
age_min	minimum age to use in catch curves. If a value is not supplied, see methods used in age_min_method
age_max	maximum age to use in catch curves. If the value is not supplied, the default is to use the second oldest age class, since the oldest age class is often a plus group.
nobs_min	the minimum number of usable observations (i.e. catch > 0) required before a catch curve should be calculated
age_min_method	character. allows you to set the method that should be used to calculate the minimum age used in the catch curves.

Author(s)

Nikolai Klibansky

Examples

## Not run: 
# 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 <- 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")


## End(Not run)

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