R/get_ref_pts.r

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

#' get_ref_pts
#'
#' Calculations of stock dependent quantities over a range of fishing mortalities (F).
#' This is code based on code from the Beaufort Assessment Model converted from ADMB to R,
#' for both bam functions \code{get_msy} and \code{get_per_recruit_stuff}.
#' @param rdat BAM output rdat file object
#' @param maxF Maximum fishing mortality rate. numeric vector
#' @param nFopt length of F vector for grid search optimization
#' @param nFres length of F vector for filtering grid search results to limit output size
#' @param SPR_props proportions of SPR to compute reference points at. numeric vector
#' @param nages Number of ages in population model. numeric vector
#' @param steep Beverton-Holt steepness parameter used in equilibrium calculations. numeric vector
#' @param R0 Beverton-Holt R0 parameter (virgin recruitment) used in equilibrium calculations. Numbers of fish at first age (often age-0 or age-1). numeric vector
#' @param BiasCor bias correction used in equilibrium calculations. numeric vector
#' @param SR_switch 1 for Beverton-Holt, 2 for Ricker, 3 for none (average recruitment) used in equilibrium calculations. integer
#' @param steep_2 Beverton-Holt steepness parameter used in non-equilibrium calculations. numeric vector
#' @param R0_2 Beverton-Holt R0 parameter (virgin recruitment) used in non-equilibrium calculations. Numbers of fish at first age (often age-0 or age-1). numeric vector
#' @param BiasCor_2 bias correction used in non-equilibrium calculations. numeric vector
#' @param SR_switch_2 1 for Beverton-Holt, 2 for Ricker, 3 for none (average recruitment) used in non-equilibrium calculations. integer
#' @param M Natural mortality rate. numeric vector of length \code{nages}
#' @param sel_wgted_L effort-weighted, recent selectivities toward landings. numeric vector of length \code{nages}
#' @param sel_wgted_D effort-weighted, recent selectivities toward discards. numeric vector of length \code{nages}
#' @param wgt_klb Weight-at-age in thousand pounds. numeric vector of length \code{nages}
#' @param reprod Reproductive potential at-age vector. numeric vector of length \code{nages}
#' @param spawn_time_frac time of year of peak spawning, as a fraction of the year. numeric
#' @param wgt_mt Weight-at-age in metric tons for fish in the population. numeric vector of length \code{nages}
#' @param wgt_klb Weight-at-age in thousand pounds for fish in the population. numeric vector of length \code{nages}
#' @param wgt_klb_L Weight-at-age in thousand pounds for fish in the landings. numeric vector of length \code{nages}
#' @param wgt_klb_D Weight-at-age in thousand pounds for fish in the discards. numeric vector of length \code{nages}
#' @param make_plot Indicate whether or not to plot results. logical
#' @param plot_digits number of significant digits to show in plot. numeric
#'
#' @details This functions computes msy-based reference points with equilibrium methods,
#' and computes SPR-based reference points using non-equilibrium methods. Both sets
#' of calculations use the same function, \code{sr_eq}, to calculate recruitment
#' at each level of F, but separate sets of arguments are passed to each call.
#' By default, \code{SR_switch_2="none"}, setting recruitment to \code{BiasCor_2*R0_2}.
#' To set recruitment equal to a fixed value (e.g. empirical mean recruitment from a set of years),
#' set \code{R0_2} equal to that value, and set \code{BiasCor_2} equal to 1.
#' @returns Invisibly returns a list including data generated by the equilibrium
#' and non-equilibrium calculations (\code{data_eq} and \code{data_neq}), as
#' well as MSY and SPR-based reference points (\code{ref_pts_msy} and \code{ref_pts_SPR}),
#' computed with equilibrium and non-equilibrium methods, respectively.
#' @keywords bam stock assessment fisheries population dynamics
#' @author Erik Williams, Kyle Shertzer, and Nikolai Klibansky
#' @export
#' @examples
#' \dontrun{
#' rdat <- rdat_VermilionSnapper
#' res <- get_ref_pts(rdat,make_plot=TRUE)
#' # Using average recruitment in the non-equilibrium calculations, as in SEDAR 55
#' rec_mean <- mean(rdat$t.series[paste(1976:2013),"recruits"])
#' res2 <- get_ref_pts(rdat,R0_2=rec_mean,BiasCor_2=1,make_plot=TRUE)
#' }

get_ref_pts <-  function(
                     rdat=NULL,
                     maxF=NULL,
                     nFopt=10001,
                     nFres=101,
                     SPR_props = c(0.3,0.4,0.5),
                     nages=NULL,
                     sel_wgted_L=NULL,
                     sel_wgted_D=NULL,
                     M=NULL,
                     reprod=NULL,
                     spawn_time_frac=NULL,
                     R0=NULL,
                     steep=NULL,
                     BiasCor=NULL,
                     SR_switch=NULL,
                     R0_2=NULL,
                     steep_2=NULL,
                     BiasCor_2=NULL,
                     SR_switch_2="none",
                     rec_mean=NULL,
                     wgt_mt=NULL,
                     wgt_klb=NULL,
                     wgt_klb_L=NULL,
                     wgt_klb_D=NULL,
                     make_plot=FALSE,
                     plot_digits=3#,
                     #plot_rp=list("eq"= c("SPR","R","L_klb","D_klb","SSB","B"),
                     #              "neq"=c("SPR","R","L_klb","D_klb","SSB","B"))
                     ){
  if(!is.null(rdat)){
    if(is.null(maxF)){
      maxF <- max(rdat$eq.series$F.eq)
    }
    if(is.null(nFopt)){
      nFopt <- nrow(rdat$eq.series)
    }

    if(is.null(nages)){
      nages <- length(rdat$a.series$age)
    }
    if(is.null(sel_wgted_L)){
      sel_wgted_L <- rdat$sel.age$sel.v.wgted.L
    }
    if(is.null(sel_wgted_D)){
      sel_wgted_D <- rdat$sel.age$sel.v.wgted.D
      if(is.null(sel_wgted_D)){
        sel_wgted_D <- sel_wgted_L*0
      }
    }
    if(is.null(M)){
      M <- rdat$a.series$M
    }
    if(is.null(spawn_time_frac)){
      spawn_time_frac <- rdat$parms$spawn.time
    }
    if(is.null(reprod)){
      reprod <- rdat$a.series$reprod
    }
    if(is.null(SR_switch)){
      SR_switch <- c("BH-steep"="beverton_holt","Ricker-steep"="ricker","Mean"="none")[rdat$info$rec.model]
      if(is.na(SR_switch)){
        warning("value of SR_switch not valid")
      }
    }
    if(is.null(R0)){
      nm_logR0 <- names(rdat$parm.cons)[grepl("^log[\\.\\_]R0",names(rdat$parm.cons))]
      R0 <- exp(rdat$parm.cons[[nm_logR0]][8])
    }
    if(is.null(steep)){
      steep <- rdat$parm.cons$steep[8]
    }
    if(is.null(BiasCor)){
      BiasCor <-  rdat$parms[grepl("biascorr",names(rdat$parms))][[1]]
    }
    if(is.null(SR_switch_2)){
      SR_switch_2 <- c("BH-steep"="beverton_holt","Ricker-steep"="ricker","Mean"="none")[rdat$info$rec.model]
      if(is.na(SR_switch_2)){
        warning("value of SR_switch_2 not valid")
      }
    }
    if(is.null(R0_2)){
      nm_logR0 <- names(rdat$parm.cons)[grepl("^log[\\.\\_]R0",names(rdat$parm.cons))]
      R0_2 <- exp(rdat$parm.cons[[nm_logR0]][8])
    }
    if(is.null(steep_2)){
      steep_2 <- rdat$parm.cons$steep[8]
    }
    if(is.null(BiasCor_2)){
      BiasCor_2 <-  rdat$parms[grepl("biascorr",names(rdat$parms))][[1]]
    }
    if(is.null(wgt_mt)){
      wgt_mt <- rdat$a.series$wgt.mt
    }
    if(is.null(wgt_klb)){
      wgt_klb <- rdat$a.series$wgt.klb
    }
    if(is.null(wgt_klb_L)){
      L_klb_ix <- grep("wgt.wgted.L.klb",names(rdat$a.series))
      if(length(L_klb_ix)==1){
        wgt_klb_L <- rdat$a.series[,L_klb_ix]
      }else{
        wgt_klb_L <- wgt_klb
      }
    }
    if(is.null(wgt_klb_D)){
      D_klb_ix <- grep("wgt.wgted.D.klb",names(rdat$a.series))
      if(length(D_klb_ix)==1){
        wgt_klb_D <- rdat$a.series[,D_klb_ix]
      }else{
        wgt_klb_D <- wgt_klb
      }
    }
  } # end if(!is.null(rdat))

  dzero <- 0.00001

  # Abbreviate stuff
  ni <- nFopt
  na <- nages
  selL <- sel_wgted_L
  selD <- sel_wgted_D
  st <- spawn_time_frac

  # Stuff initialized early in the BAM tpl (e.g. PARAMETER_SECTION)
  # F_msy <- rep(NA,ni) # not needed
  Fvec <- seq(0,maxF,length=ni)
  N_a <- N_a_st <- rep(NA, na) # _st = spawning time
          spr_1 <- SPR_1 <- R_1 <- SSB_1 <- B_1 <- L_klb_1 <- L_knum_1 <- D_klb_1 <- D_knum_1 <- rep(NA,nFopt)
  #F_2 <- spr_2 <- SPR_2 <- R_2 <- SSB_2 <- B_2 <- L_klb_2 <- L_knum_2 <- D_klb_2 <- D_knum_2 <- rep(NA,length(SPR_props))
                                                   L_klb_2a <- L_knum_2a <- D_klb_2a <- D_knum_2a <- rep(NA,nFopt) # per-recruit calcs
   F_2 <- spr_2 <- SPR_2 <- R_2 <- SSB_2 <- B_2 <- L_klb_2 <- L_knum_2 <- D_klb_2 <- D_knum_2 <- rep(NA,nFopt)
  # iter_inc_msy <- maxF/(ni-1) # not needed

  # Usually computed by get_weighted_current()


  #fill in Fs for msy and per-recruit analyses # not needed
  # Fvec[1]=0.0
  # for (ff in 2:ni) {
  #   Fvec[ff] <- Fvec[ff-1] + iter_inc_msy
  #   }

  ##### ref pt calcs 1: equilibrium method (get_msy)
  # compute values as functions of F
  # for(ff=1; ff<=ni; ff++)
  # FOR EACH LEVEL OF F..
  for (ff in 1:ni) {
    F_L_a <- Fvec[ff]*selL
    F_D_a <- Fvec[ff]*selD
    Z_a <- M + F_L_a + F_D_a

    N_a[1] <- 1.0 # per recruit
    for (iage in 2:na){
      N_a[iage] <- N_a[iage-1]*exp(-Z_a[iage-1])
      }
    N_a[na] <- N_a[na]/(1.0-exp(-Z_a[na]))
    N_a_st[1:(na-1)] <- N_a[1:(na-1)]*exp(-Z_a[1:(na-1)]*st)
    N_a_st[na] <- (N_a_st[na-1]*(exp(-(Z_a[na-1]*(1.0-st)+Z_a[na]*st))))/(1.0-exp(-Z_a[na]))

    # Spawners per recruit
    spr_1[ff] <- sum(N_a_st*reprod)

    # Spawning potential ratio
    SPR_1[ff] <- spr_1[ff]/spr_1[1]

    R_1[ff] <- max(dzero,sr_eq(SR_switch=SR_switch, R0=R0, steep=steep, BiasCor=BiasCor, spr_F0=spr_1[1], spr_F=spr_1[ff]))

    # Rescale N_a by recruitment estimate
    N_a <-    N_a*R_1[ff]
    N_a_st <- N_a_st*R_1[ff]

    L_a <- N_a*(F_L_a/Z_a)*(1-exp(-Z_a))
    D_a <- N_a*(F_D_a/Z_a)*(1-exp(-Z_a))

    SSB_1[ff]    <- sum(N_a_st*reprod)
  	B_1[ff]      <- sum(N_a*wgt_mt)
    L_klb_1[ff]  <- sum(L_a*wgt_klb_L)
    L_knum_1[ff] <- sum(L_a)/1000
    D_klb_1[ff]  <- sum(D_a*wgt_klb_D)
    D_knum_1[ff] <- sum(D_a)/1000
  }

  # Compute per-recruit reference points
    msy_klb <- max(L_klb_1)

    ff_msy <- which(L_klb_1==msy_klb)

  ##### ref pt calcs 2: non-equilibrium method (get_per_recruit_stuff)
    # for(i in seq_along(ff_SPR_props)){
    #   ffi <- ff_SPR_props[i]
    #   Fi <- Fvec[ffi]
    #   F_L_a <- Fi*selL
    #   F_D_a <- Fi*selD
    #   Z_a <- M + F_L_a + F_D_a
    # FOR EACH LEVEL OF F..
    for (ff in 1:ni) {

      F_L_a <- Fvec[ff]*selL
      F_D_a <- Fvec[ff]*selD
      Z_a <- M + F_L_a + F_D_a

      N_a[1] <- 1.0 # per recruit
      for (iage in 2:na){
        N_a[iage] <- N_a[iage-1]*exp(-Z_a[iage-1])
      }
      N_a[na] <- N_a[na]/(1.0-exp(-Z_a[na]))
      N_a_st[1:(na-1)] <- N_a[1:(na-1)]*exp(-Z_a[1:(na-1)]*st)
      N_a_st[na] <- (N_a_st[na-1]*(exp(-(Z_a[na-1]*(1.0-st)+Z_a[na]*st))))/(1.0-exp(-Z_a[na]))

      # Spawners per recruit
      spr_2[ff] <- sum(N_a_st*reprod)

      # Spawning potential ratio
      SPR_2[ff] <- spr_2[ff]/spr_2[1]

      L_a_2a <- N_a*(F_L_a/Z_a)*(1-exp(-Z_a)) # per recruit
      D_a_2a <- N_a*(F_D_a/Z_a)*(1-exp(-Z_a)) # per recruit

      R_2[ff] <- max(dzero,sr_eq(SR_switch=SR_switch_2, R0=R0_2, steep=steep_2, BiasCor=BiasCor_2, spr_F0=spr_2[1], spr_F=spr_2[ff]))
      N_a[1] <- R_2[ff] # use mean recruitment
      for (iage in 2:na){
        N_a[iage] <- N_a[iage-1]*exp(-Z_a[iage-1])
      }
      N_a[na] <- N_a[na]/(1.0-exp(-Z_a[na]))
      N_a_st[1:(na-1)] <- N_a[1:(na-1)]*exp(-Z_a[1:(na-1)]*st)
      N_a_st[na] <- (N_a_st[na-1]*(exp(-(Z_a[na-1]*(1.0-st)+Z_a[na]*st))))/(1.0-exp(-Z_a[na]))

      L_a <- N_a*(F_L_a/Z_a)*(1-exp(-Z_a))
      D_a <- N_a*(F_D_a/Z_a)*(1-exp(-Z_a))

      #F_2[ff]      <- Fi
      # spr_2[i]    <- spr_1[ffi] # get the correct value from the msy loop above
      # SPR_2[i]    <- SPR_1[ffi] # get the correct value from the msy loop above
      SSB_2[ff]    <- sum(N_a_st*reprod)
      B_2[ff]      <- sum(N_a*wgt_mt)
      L_klb_2a[ff]  <- sum(L_a_2a*wgt_klb_L)
      L_knum_2a[ff] <- sum(L_a_2a)/1000
      D_klb_2a[ff]  <- sum(D_a_2a*wgt_klb_D)
      D_knum_2a[ff] <- sum(D_a_2a)/1000
      L_klb_2[ff]  <- sum(L_a*wgt_klb_L)
      L_knum_2[ff] <- sum(L_a)/1000
      D_klb_2[ff]  <- sum(D_a*wgt_klb_D)
      D_knum_2[ff] <- sum(D_a)/1000
    }
    ff_SPR_props <- sapply(SPR_props,function(x){which.min(abs(SPR_2-x))})


  # Build data frame of values at F
  data_eq <- data.frame("F"=Fvec,
                        L_klb=L_klb_1, L_knum=L_knum_1,
                        D_klb=D_klb_1, D_knum=D_knum_1,
                        B_mt=B_1, SSB=SSB_1, R=R_1,
                        spr=spr_1, SPR=SPR_1
                        )

  # nm_ref_pts_SPR <- paste0("F",SPR_props*100)
  data_neq <- data.frame(#"ref_pt"=nm_ref_pts_SPR,
                        #"SPR_prop"=SPR_props,
                        #"F"=F_2,
                        "F"=Fvec,
                        L_klb_pr=L_klb_2a, L_knum_pr=L_knum_2a,
                        D_klb_pr=D_klb_2a, D_knum_pr=D_knum_2a,
                        L_klb=L_klb_2, L_knum=L_knum_2,
                        D_klb=D_klb_2, D_knum=D_knum_2,
                        B_mt=B_2, SSB=SSB_2, R=R_2,
                        spr=spr_2, SPR=SPR_2
  )
  # rownames(data_neq) <- nm_ref_pts_SPR

  # Identify msy reference points
  ref_pts_msy <- data_eq[ff_msy,]
  rownames(ref_pts_msy) <- "Fmsy"

  # Identify SPR reference points
  ref_pts_SPR <- data_neq[ff_SPR_props,]
  rownames(ref_pts_SPR) <- paste0("F",SPR_props*100)
  # ref_pts_SPR <- data_neq


  # Subset data_eq to return
  data_eq_res <- data_eq[floor(seq(1,nFopt,length=nFres)),]
  data_neq_res <- data_neq[floor(seq(1,nFopt,length=nFres)),]

  # Warnings
  if(!is.null(rdat)){
    if(rdat$parm.cons$steep[4]<0){
      warning(paste0("steepness was fixed at ",rdat$parm.cons$steep[8]," in this assessment"))
    }
  }
  if(ff_msy==nFopt){
    warning("msy_klb was identified at maxF and probably does not represent the true maximum.")
  }
  if(nFopt%in%ff_SPR_props){
    warning("One of the SPR reference points was identified at maxF. Try a higher maxF to correctly identify this value.")
  }

  if(make_plot){
    plot_rp <- function(...,rp=NA,rp_name=NULL){
      plot(type="l",lwd=2,...)
      if(!is.na(rp)){
        points(ref_pts_msy[,c("F",rp)],type="p",col="blue",pch=16)
        if(is.null(rp_name)){rp_name <- rp}
          rp_labels <- bquote(.(rp_name) == .(signif(ref_pts_msy[,rp],plot_digits)))
                  text(ref_pts_msy[,c("F",rp)],labels=rp_labels,pos=4)
        }
    }
    # Plot equilibrium results
    par(mar=c(2,2,2,1),oma=c(1,1,2,1),mfrow=c(3,2),mgp=c(1.1,0.1,0),tck=0.01)
    with(data_eq_res,{
    plot_rp(F,spr,rp="spr",rp_name = bquote(spr[MSY]))
    plot_rp(F,SPR,rp="SPR",rp_name = bquote(SPR[MSY]))
    plot_rp(F,L_klb, rp="L_klb",rp_name = bquote(MSY))
    plot_rp(F,R,rp="R",rp_name = bquote(R[MSY]))
    plot_rp(F,SSB,rp="SSB",rp_name = bquote(SSB[MSY]))
    plot_rp(F,B_mt,rp="B_mt",rp_name = bquote(B[MSY]))

    mtext("equilibrium", outer=TRUE)
    })

    # Plot non-equilibrium results
    par(mar=c(2,2,2,1),oma=c(1,1,2,1),mfrow=c(3,2),mgp=c(1.1,0.1,0),tck=0.01)
    with(data_neq_res,{
      plot_rp(F,spr)
      points(ref_pts_SPR[,c("F","spr")],type="p",col="blue",pch=16)
      sapply(seq_along(SPR_props),function(x){
        text(ref_pts_SPR[x,c("F","spr")],
             labels=bquote(spr[.(SPR_props[x]*100)] == .(signif(ref_pts_SPR[x,"spr"],plot_digits))),pos=4)})
      plot_rp(F,SPR)
      points(ref_pts_SPR[,c("F","SPR")],type="p",col="blue",pch=16)
      sapply(seq_along(SPR_props),function(x){
        text(ref_pts_SPR[x,c("F","SPR")],
             labels=bquote(F[.(SPR_props[x]*100)] == .(signif(ref_pts_SPR[x,"F"],plot_digits))),pos=4)})
      plot_rp(F,L_klb)
      points(ref_pts_SPR[,c("F","L_klb")],type="p",col="blue",pch=16)
      sapply(seq_along(SPR_props),function(x){
        text(ref_pts_SPR[x,c("F","L_klb")],
             labels=bquote(L[.(SPR_props[x]*100)] == .(signif(ref_pts_SPR[x,"L_klb"],plot_digits))),pos=4)})

      plot_rp(F,R)
      points(ref_pts_SPR[,c("F","R")],type="p",col="blue",pch=16)
      sapply(seq_along(SPR_props),function(x){
        text(ref_pts_SPR[x,c("F","R")],
             labels=bquote(R[.(SPR_props[x]*100)] == .(signif(ref_pts_SPR[x,"R"],plot_digits))),pos=4,srt=45)})

      plot_rp(F,SSB)
      points(ref_pts_SPR[,c("F","SSB")],type="p",col="blue",pch=16)
      sapply(seq_along(SPR_props),function(x){
        text(ref_pts_SPR[x,c("F","SSB")],
             labels=bquote(SSB[.(SPR_props[x]*100)] == .(signif(ref_pts_SPR[x,"SSB"],plot_digits))),pos=4)})

      plot_rp(F,B_mt)
      points(ref_pts_SPR[,c("F","B_mt")],type="p",col="blue",pch=16)
      sapply(seq_along(SPR_props),function(x){
        text(ref_pts_SPR[x,c("F","B_mt")],
             labels=bquote(B[.(SPR_props[x]*100)] == .(signif(ref_pts_SPR[x,"B_mt"],plot_digits))),pos=4)})
      mtext("non-equilibrium", outer=TRUE)
    })

    }

  invisible(list(ref_pts_msy=ref_pts_msy, ref_pts_SPR=ref_pts_SPR, data_eq=data_eq_res, data_neq=data_neq_res))
}