R/Sinclair_Z_functions.r

In cmlegault/ASAPplots: Creates Standard Plots and PDFs for ASAP3

# Estimate Z from catch curves using 4 year moving window with cohort-specific coefficients
# based on Sinclair 2001 ICES J Mar Sci 58: 1-10

# matrix has year in rows and age in columns
calc_Sinclair_Z <- function(mat){
  res <- list()
  mat <- as.matrix(mat)
  if (all(is.na(mat))){
    res$error <- TRUE
    return(res)
  }
  year <- as.numeric(rownames(mat))
  age <- as.numeric(colnames(mat))
  ny <- length(year)
  est.Sinclair.Z <- matrix(NA, nrow=(ny-3), ncol=3)
  plot.year <- rep(NA, ny-3)
  for (i in 1:(ny-3)){
    submat <- mat[year %in% year[i:(i+3)],]
    data <- reshape2::melt(submat)
    colnames(data) <- c("Year","Age","Value")
    data$cohort <- data$Year - data$Age
    data$Value[data$Value == 0] <- NA
    data <- data[!is.na(data$Value) ,]
    data$lnVal <- log(data$Value)
    can.calc <- FALSE
    if (length(data[,1]) >= 2){
      if (max(table(data$cohort)) >= 3){
        can.calc <- TRUE
      }
    }
    if (can.calc == TRUE){
      my.lm <- lm(data$lnVal ~ as.factor(data$cohort) + data$Age)
      data$pred <- predict(my.lm)
      data$resid <- residuals(my.lm)
      res[[i]] <- data
      est.Sinclair.Z[i,1] <- -1 * my.lm$coefficients[names(my.lm$coefficients) == "data$Age"]
      est.Sinclair.Z[i,2:3] <- -1 * rev(confint(my.lm, "data$Age", level=0.90))
    }
    else{  
      res[[i]] <- data
      est.Sinclair.Z[i,] <- rep(NA, 3)
    }
    plot.year[i] <- year[i] + 1.5
  }
  colnames(est.Sinclair.Z) <- c("Sinclair_Z","low90%","high90%")
  res$est.Sinclair.Z <- est.Sinclair.Z
  res$plot.year <- plot.year
  res$error <- FALSE
  return(res)
}

#-------------------------------------------------------------------------------------
# function to plot the results from calc_Sinclair_Z
# res is a list generated by the function calc_Sinclair_Z
plot_Sinclair_Z <- function(res,est.Z.name,regress.name,resid.name,my.title,save.plots,od,plotf){
  # first the plot of estimated Z with CI over time
  est.Z <- res$est.Sinclair.Z
  if (!all(is.na(est.Z))){
    par(mfcol=c(1,1), oma=rep(1,4), mar=rep(4,4))
    plot.year <- res$plot.year
    Hmisc::errbar(plot.year,est.Z[,1],est.Z[,3],est.Z[,2],ylim=c(0,max(est.Z, na.rm=T)),xlab="Year",ylab="Sinclair Z")
    if(!is.na(my.title)) title(main=my.title, outer=F)
    if(save.plots) savePlot(paste0(od,est.Z.name,".",plotf), type=plotf)
  }
  
  # plot the individual regressions
  ny <- length(res$plot.year)
  np.row <- 3             # number of plot rows
  np.col <- 3             # number of plot columns
  np.p <- np.row * np.col # number of plots per page
  par(mfcol=c(np.row, np.col))
  my.col <- 1:100
  my.col[7] <- "blue"
  icount <- 0
  for (i in 1:ny){
    data <- res[[i]]
    if (!is.na(res$est.Sinclair.Z[i,1])){
      icount <- icount + 1
      n.coh <- length(unique(data$cohort))
      i.coh <- (data$cohort-min(data$cohort, na.rm=T)+1)
      plot(data$Age,data$lnVal,pch=i.coh,col=my.col[i.coh],xlab="Age",ylab="ln Val")
      title(main=paste0("Years ",min(data$Year)," to ",max(data$Year),"\nZ = ",round(res$est.Sinclair.Z[i,1],3)), outer=F)
      for (j in 1:n.coh){
        subcoh <- data[i.coh == j,]
        if (length(subcoh$lnVal) >= 2){
          lines(subcoh$Age,subcoh$pred,col=my.col[j],lty=j) 
        }
      }
      if(save.plots == T & (icount %% np.p) == 0){
        savePlot(paste0(od,regress.name,"_",i,".",plotf), type=plotf)  
      }
    }
    if (save.plots == TRUE & i == ny) savePlot(paste0(od,regress.name,"_",i,".",plotf), type=plotf)
  }
  par(mfcol=c(1,1))
  
  # plot the distribution of residuals by age
  if (!all(is.na(est.Z))){
    resids <- matrix(NA, nrow=0, ncol=2)
    colnames(resids) <- c("Age","Resid")
    for (i in 1:ny){
      if (!is.na(res$est.Sinclair.Z[i,1])){
        data <- res[[i]]
        resids <- rbind(resids,cbind(data$Age,data$resid))
      }  
    }
    boxplot(resids[,2] ~ resids[,1],xlab="Age",ylab="Residual")
    abline(h=0, col="red", lty=2)
    if(!is.na(my.title)) title(main=my.title, outer=F)
    if(save.plots) savePlot(paste0(od,resid.name,".",plotf), type=plotf)
  }
  else{
    resids <- NA
  }
  
  return(resids)
}
#-------------------------------------------------------------------------------------
# function to remove ages from matrix of index catch at age, used by get_Sinclair_Z
drop_partially_selected_ages <- function(mat,sel,sel.start.age,sel.end.age,sel.cutoff){
  sel2 <- sel[sel.start.age:sel.end.age]
  age.col <- 1:length(sel2)
  drop.col <- age.col[sel2 < sel.cutoff]
  if (length(drop.col) > 0) mat[,drop.col] <- NA
  return(mat) 
}
#-------------------------------------------------------------------------------
#' Get Sinclair Z
#' 
#' Set of functions to prepare index matrices and compute and plot Sinclair Z. Estimate Z from catch curves using 4 year moving window with cohort-specific coefficients. Based on Sinclair 2001 ICES J Mar Sci 58: 1-10. Note hard wired to use only ages that are at least 90 percent selected.
#' @param asap name of the variable that read in the asap.rdat file
#' @param a1 list file produced by grab.aux.files function
#' @param index.names names of indices
#' @param save.plots save individual plots
#' @param od output directory for plots and csv files 
#' @param plotf type of plot to save
#' @export

get_Sinclair_Z <- function(asap,a1,index.names,save.plots,od,plotf){
  asap.name <- a1$asap.name
  index.mats <- ConvertSurveyToAtAge(asap)
  nages <- asap$parms$nages
  Sinclair.Z <- NA
  for (j in 1:length(index.mats$ob)){
    if (asap$control.parms$index.age.comp.flag[j] == 1){  # used age composition for the index
      mat <- index.mats$ob[[j]]
      sel.start.age <- asap$control.parms$index.sel.start.age[j]
      sel.end.age <- asap$control.parms$index.sel.end.age[j]
      mat <- drop_partially_selected_ages(mat,asap$index.sel[j,],sel.start.age,sel.end.age,0.90)
      if (sel.end.age == nages) mat[,length(mat[1,])] <- NA # drop oldest age because plus group in ASAP
      Sinclair.Z <- calc_Sinclair_Z(mat)
      if (Sinclair.Z$error == FALSE){
        est.Z.name <- paste0("Sinclair_Z_estimates_index_",j)
        regress.name <- paste0("Sinclair_Z_regress_index_",j)
        resid.name <- paste0("Sinclair_Z_resids_index_",j)
        resids <- plot_Sinclair_Z(Sinclair.Z,est.Z.name,regress.name,resid.name,index.names[j],
                                  save.plots,od,plotf)
        Sinclair.Z.table <- cbind(Sinclair.Z$plot.year,Sinclair.Z$est.Sinclair.Z)
        write.csv(Sinclair.Z.table, file=paste0(od,"Sinclair_Z_index_",j,"_",asap.name,".csv"), row.names=F)
      }
    }
  }
  return(Sinclair.Z)
}