R/ComputeIndices.R

In nwfsc-assess/nwfscDeltaGLM: Bayesian delta-GLM package for constructing indices of CPUE data

#' Compute index of abundance
#'
#' @param Data The data to be fit
#' @param Model The model object (as a list)
#' @param FileName The name of the file/directory
#' @param maxDims aesthetics argument for plot as dimensions of rows/columns, defaults to 6
#' @param Folder Optional argument for named folder to store results (in current working directory)
#' @param Weights Defaults to "StrataAreas", other option is "Equal"
#' @param StrataTable Table of strata, required
#' @param PlotStrataYearMcmc Whether to plot MCMC chains for each Strata:Year, defaults to \code{TRUE}
#'
#' @return Returns a list of output, by year and by Strata:Year
#' @import R2jags
#' @export
#'
ComputeIndices <- function(Data, Model, FileName, maxDims = 6, Folder = NA,
                           Weights = "StrataAreas", StrataTable, PlotStrataYearMcmc = TRUE) {
  if (is.na(Folder)) Folder <- paste(getwd(), "/", sep = "")

  # Attach stuff -- listed by search()
  attach(Model$BUGSoutput$sims.list, warn.conflicts = FALSE)
  on.exit(detach(Model$BUGSoutput$sims.list))
  modelStructure <- Model$modelStructure
  Dist <- Model$likelihood

  # Compute average AreaSwept - This AreaSwept used to calculate densities, and matters if Catch is a nonlinear function of AreaSwept
  MeanLogAreaSwept <- mean(Data$logeffort)

  # New objects
  Chains <- array(NA, dim = c(nrow(Sdev), length(unique(StrataTable[, "year"])), length(unique(StrataTable[, "strata"])), 2))
  Year <- Strata <- Area <- PosMedian <- PresMedian <- IndexMedian <- IndexMedianWeighted <- PosMean <- PresMean <- IndexMean <- IndexMeanWeighted <- CvMedian <- SdLog <- RawPos <- RawPres <- Raw <- RawWeighted <- RawVar <- RawVarWeighted <- RawCV <- matrix(NA, nrow = length(unique(StrataTable[, "year"])), ncol = length(unique(StrataTable[, "strata"])))
  CvMedianYear <- SdLogYear <- rep(NA, length(unique(StrataTable[, "year"])))

  # Calculate values
  for (YearI in 1:nrow(PosMean)) {
    for (StratI in 1:ncol(PosMean)) {
      # Derived indicators
      StrataYearI <- which(levels(strataYear) == paste(levels(strata)[StratI], ":", levels(year)[YearI], sep = ""))
      AreaI <- which(StrataTable[, "strataYear"] == paste(levels(strata)[StratI], ":", levels(year)[YearI], sep = ""))
      Which <- which(strata == levels(strata)[StratI] & year == levels(year)[YearI])
      # Year, strata, and area
      Year[YearI, StratI] <- levels(year)[YearI]
      Strata[YearI, StratI] <- levels(strata)[StratI]
      if (Weights == "StrataAreas") Area[YearI, StratI] <- StrataTable[AreaI, "Area_Hectares"]
      if (Weights == "Equal") Area[YearI, StratI] <- 1
      # Save Positive catch chain in normal-space and correct for transformation biases
      Chains[, YearI, StratI, 1] <- exp(cMx(Sdev)[, StratI] + cMx(Ydev)[, YearI] + cMx(SYdev)[, StrataYearI] + log(1) * B.pos[, 1] + log(1)^2 * B.pos[, 2]) # wardJAGS uses logeffort offset
      # Lognormal -- Bias correction
      if (Dist == "lognormal") {
        Sigma <- sqrt(log(CV[, 1]^2 + 1))
        Chains[, YearI, StratI, 1] <- Chains[, YearI, StratI, 1] * exp(Sigma^2 / 2)
      }
      # LognormalECE -- Bias correction + incorporate ECE
      if (Dist == "lognormalECE" | Dist == "lognormalECE2") {
        Sigma <- sqrt(log(CV[, 1]^2 + 1))
        Sigma2 <- sqrt(log(CV[, 2]^2 + 1))
        Chains[, YearI, StratI, 1] <- Chains[, YearI, StratI, 1] * exp(Sigma^2 / 2) * p.ece[, 1] + ratio * Chains[, YearI, StratI, 1] * exp(Sigma2^2 / 2) * p.ece[, 2]
      }
      # GammaECE -- incorporate ECE
      if (Dist == "gammaECE" | Dist == "gammaECE2") {
        Chains[, YearI, StratI, 1] <- Chains[, YearI, StratI, 1] * p.ece[, 1] + ratio * Chains[, YearI, StratI, 1] * p.ece[, 2]
      }
      # Don't make mean-unbiased for unobserved vessel
      # if(modelStructure$VesselYear.positiveTows=="random") Chains[,YearI,StratI,1] = Chains[,YearI,StratI,1] * exp(sigmaVY[,1]^2/2)
      # Save Presence/Absence chain in normal-space and correct for transformation biases
      Chains[, YearI, StratI, 2] <- plogis(cMx(pSdev)[, StratI] + cMx(pYdev)[, YearI] + cMx(pSYdev)[, StrataYearI] + (1) * B.zero[, 1] + (1)^2 * B.zero[, 2]) # wardJAGS predicts the probability of 0 catch, and uses offset as effort, not logeffort
      # Don't make mean-unbiased for unobserved vessel (this was done incorrectly anyway)
      # if(modelStructure$VesselYear.zeroTows=="random") Chains[,YearI,StratI,2] = mean(plogis(rnorm(n=dim(Chains)[1]*1e3, mean=qlogis(Chains[,YearI,StratI,2]), sd=sigmaVY[,2])))
      # Index (median)
      PosMedian[YearI, StratI] <- median(Chains[, YearI, StratI, 1]) # / 2e4 # Convert kilograms to metric tons
      PresMedian[YearI, StratI] <- median(Chains[, YearI, StratI, 2])
      IndexMedian[YearI, StratI] <- median(Chains[, YearI, StratI, 1] * Chains[, YearI, StratI, 2])
      IndexMedianWeighted[YearI, StratI] <- IndexMedian[YearI, StratI] * Area[YearI, StratI]
      # Index (mean)
      PosMean[YearI, StratI] <- mean(Chains[, YearI, StratI, 1]) # / 2e4 # Convert kilograms to metric tons
      PresMean[YearI, StratI] <- mean(Chains[, YearI, StratI, 2])
      IndexMean[YearI, StratI] <- mean(Chains[, YearI, StratI, 1] * Chains[, YearI, StratI, 2])
      IndexMeanWeighted[YearI, StratI] <- IndexMean[YearI, StratI] * Area[YearI, StratI]
      # CV of median (from J. Wallace "Survey.Biomass.GlmmBUGS.ver.3.00.r)
      Temp <- Area[YearI, StratI] * Chains[, YearI, StratI, 1] * Chains[, YearI, StratI, 2]
      CvMedian[YearI, StratI] <- sqrt(var(Temp)) / median(Temp)
      # SD of log of index (from J. Wallace "Survey.Biomass.GlmmBUGS.ver.3.00.r")
      SdLog[YearI, StratI] <- sd(log(Temp))
      # Raw
      RawPos[YearI, StratI] <- mean(ifelse(Data[Which, "HAUL_WT_KG"] > 0, Data[Which, "HAUL_WT_KG"] / Data[Which, "effort"], NA), na.rm = TRUE)
      RawPres[YearI, StratI] <- mean(Data[Which, "HAUL_WT_KG"] > 0)
      Raw[YearI, StratI] <- RawPos[YearI, StratI] * RawPres[YearI, StratI]
      RawWeighted[YearI, StratI] <- Raw[YearI, StratI] * Area[YearI, StratI]
      RawVar[YearI, StratI] <- var(Data[Which, "HAUL_WT_KG"] / Data[Which, "effort"], na.rm = TRUE) / length(Which)
      RawVarWeighted[YearI, StratI] <- RawVar[YearI, StratI] * Area[YearI, StratI]^2
      RawCV[YearI, StratI] <- sqrt(RawVarWeighted[YearI, StratI]) / RawWeighted[YearI, StratI]
    }
    # CV of median (from J. Wallace "Survey.Biomass.GlmmBUGS.ver.3.00.r")
    # Temp = Area[YearI,StratI] * rowSums( cMx(Chains[, YearI, ,1]) * cMx(Chains[,YearI,,2]) )  # By using StratI outside of the Strata loop only the final area estimate is obained.
    Temp <- (cMx(Chains[, YearI, , 1]) * cMx(Chains[, YearI, , 2])) %*% array(Area[YearI, ]) # array(Area[YearI,], c(3,1))
    CvMedianYear[YearI] <- sqrt(var(Temp)) / median(Temp)
    # SD of log of index (from J. Wallace "Survey.Biomass.GlmmBUGS.ver.3.00.r")
    SdLogYear[YearI] <- sd(log(Temp))
  } # 1085-115

  # Plot MCMC chains for each Strata:Year
  maxDims2 <- maxDims^2
  if (PlotStrataYearMcmc == TRUE) {
    for (Type in c("Trace", "ACF")) {
      for (ChainI in 1:2) {
        Nstrat <- length(unique(strataYear[nonZeros]))
        Nplots <- ceiling(Nstrat / maxDims2)
        for (PlotI in 1:Nplots) {
          ParSet <- (PlotI - 1) * maxDims^2 + 1:min(Nstrat - (PlotI - 1) * maxDims2, maxDims2)
          Ncol <- ceiling(sqrt(length(ParSet)))
          Nrow <- ceiling(length(ParSet) / Ncol)
          jpeg(paste(Folder, "/", FileName, "Chain_", c("Positive", "Presence")[ChainI], "_", Type, "_by_StrataYear_", PlotI, ".jpg", sep = ""), width = Ncol * 1.5, height = Nrow * 1.5, res = 150, units = "in")
          par(mfrow = c(Nrow, Ncol), mar = c(2, 2, 2, 0), mgp = c(1.25, 0.25, 0), tck = -0.02)
          for (StrataYearI in ParSet) {
            StrataI <- strata[which(strataYear == unique(strataYear)[StrataYearI])[1]]
            YearI <- year[which(strataYear == unique(strataYear)[StrataYearI])[1]]
            Mat <- matrix(Chains[, YearI, StrataI, ChainI], ncol = Model$mcmc.control$chains, byrow = FALSE)
            if (Type == "Trace") {
              matplot(Mat, type = "l", lty = "solid", col = rainbow(Model$mcmc.control$chains, alpha = 0.4), main = paste(YearI, StrataI), xaxt = "n", xlab = "", ylab = "")
            }
            if (Type == "ACF") {
              Acf <- apply(Mat, MARGIN = 2, FUN = function(Vec) {
                acf(Vec, plot = FALSE)$acf
              })
              matplot(Acf, type = "h", lty = "solid", col = rainbow(Model$mcmc.control$chains, alpha = 0.4), main = paste(YearI, StrataI), xaxt = "n", xlab = "", ylab = "", ylim = c(0, 1), lwd = 2)
            }
          }
          dev.off()
        }
      }
    }
  }

  # Compile into matrices
  Results1 <- data.frame(Year = as.vector(Year), Strata = as.vector(Strata), Raw = as.vector(RawWeighted), RawCV = as.vector(RawCV), IndexMedian = as.vector(IndexMedianWeighted), IndexMean = as.vector(IndexMeanWeighted), CvMedian = as.vector(CvMedian), SdLog = as.vector(SdLog), Area = as.vector(Area), PosMedian = as.vector(PosMedian), PresMedian = as.vector(PresMedian), PosMean = as.vector(PosMean), PresMean = as.vector(PresMean), RawPos = as.vector(RawPos), RawPres = as.vector(RawPres), RawSD = ifelse(as.vector(RawVarWeighted) == 0, NA, as.vector(sqrt(RawVarWeighted))))
  Results2 <- data.frame(Year = Year[, 1], Raw = rowSums(RawWeighted, na.rm = TRUE), RawCV = sqrt(rowSums(RawVarWeighted, na.rm = TRUE)) / rowSums(RawWeighted, na.rm = TRUE), IndexMedian = rowSums(IndexMedianWeighted, na.rm = TRUE), IndexMean = rowSums(IndexMeanWeighted, na.rm = TRUE), CvMedian = CvMedianYear, SdLog = SdLogYear)

  # Write and print output
  write.csv(Results1, file = paste(Folder, "/", FileName, "ResultsByYearAndStrata.csv", sep = ""))
  write.csv(Results2, file = paste(Folder, "/", FileName, "ResultsByYear.csv", sep = ""))

  # Return output
  Return <- list(byYearAndStrata = Results1, byYear = Results2)
  return(Return)
}