R/vbgfCI_time.R

In rschwamborn/fishboot: Bootstrap-based methods for the study of fish stocks and aquatic populations

#' VBGF plot and CI
#'
#' @param res list of class `lfqBoot`
#' @param CI vector. Confidence interval levels to plot
#' @param agemax numeric. Maximum number of years to project.
#' @param add_legend logical. Should CI and max. density legend be added
#'   (Default: `add_legend = TRUE`).
#' @param add_max_dens_legend logical. Should maximum density line be added
#'   (Default: `add_max_dens_legend = TRUE`).
#' @param xlab Label for x-axis
#' @param ylab Label for y-axis
#' @param perm.col Color for each resample estimate line. See `?par`.
#' @param perm.lwd Line width for each resample estimate line. See `?par`.
#' @param ci.col Color for CI line. See `?par`.
#' @param ci.lty Line type for CI line. See `?par`.
#' @param ci.lwd Line width for CI line. See `?par`.
#' @param maxd.col Color for maximum density line. See `?par`.
#' @param maxd.lty Line type maximum density line. See `?par`.
#' @param maxd.lwd Line width maximum density line. See `?par`.
#'
#' @return plot and list containing: `limCI` - data.frame with CI limits by time,
#'   `inCi` -  data.frame with logical values defining whether bootstrapping samples
#'   are within each of the defined CIs, `density` - the multivariate kernel density
#'   estimates for each sample, and `max_dens` is a list with the VBGF parameter
#'   combination having the maximum density estimate.
#'
#' @examples
#'
#' data(alba_boot)
#' CIinfo <- vbgfCI_time(
#'   res = alba_boot,
#'   agemax = 2, CI = 95,
#'   perm.col = adjustcolor("grey50",0.2)
#' )
#'
#' # plot more CI levels
#' CIinfo <- vbgfCI_time(
#'   res = alba_boot,
#'   agemax = 2, CI = c(95, 50),
#'   ci.lty = 1, ci.lwd = 2, ci.col = c("red", "orange"),
#'   perm.col = adjustcolor("grey50",0.2)
#' )
#'
#' # using output in lfq plot (see ?TropFishR::plot.lfq)
#' library(TropFishR)
#' data(alba)
#' alba <- lfqRestructure(alba, MA = 7)
#' plot(alba, Fname = "rcounts", ylim=c(0,15))
#' for(i in seq(nrow(alba_boot$bootRaw))){
#'   x <- as.list(alba_boot$bootRaw[i,])
#'   tmp <- lfqFitCurves(
#'     lfq = alba, par = x,
#'     col = adjustcolor("grey50",0.1), draw = TRUE, lty=1
#'   )
#' }
#' for(j in seq(CIinfo$limCI)){
#'   usr <- par()$usr
#'   tmax <- diff(range(CIinfo$limCI[[j]]$t))
#'   date_start <- as.Date(usr[1], origin = "1970-01-01")
#'   yeardec_start <- date2yeardec(date_start)
#'   date_end <- as.Date(usr[2], origin = "1970-01-01")
#'   yeardec_end <- date2yeardec(date_end)
#'   years <- (floor(yeardec_start)-tmax):ceiling(yeardec_end)
#'
#'   for(y in years){
#'     time <- y + CIinfo$limCI[[j]]$t
#'     time <- as.numeric(yeardec2date(time))
#'     lines(time, CIinfo$limCI[[j]]$min, lty = 2, col=2+j)
#'     lines(time, CIinfo$limCI[[j]]$max, lty = 2, col=2+j)
#'   }
#'
#'   tmp <- lfqFitCurves(lfq = alba, par = CIinfo$max_dens,
#'     col = 1, draw = TRUE, lty=1, lwd=2
#'   )
#' }
#'
#' @importFrom grDevices adjustcolor blues9 colorRampPalette rgb
#' @importFrom graphics abline box contour hist image
#'   layout legend lines mtext par plot points polygon
#'   rect segments text
#' @importFrom stats cov quantile runif
#'
#' @export
#'
#'
vbgfCI_time <- function(res, CI = 95, agemax = NULL,
  add_legend = TRUE, add_max_dens_legend = TRUE,
  xlab = "Relative time", ylab = "Length",
  perm.col = adjustcolor("grey50",0.1), perm.lwd = 1,
  ci.col = 1, ci.lty = 2, ci.lwd = 1,
  maxd.col = 1, maxd.lty = 1, maxd.lwd = 2
){

  res <- res$bootRaw

  if(is.null(agemax)){agemax <- max(ceiling((1/-res$K)*log(1-((res$Linf*0.95)/res$Linf))))}

  # expand ci line attributes to length of CI (if needed, values are recycled)
  ci.col <- rep_len(ci.col, length(CI))
  ci.lty <- rep_len(ci.lty, length(CI))
  ci.lwd <- rep_len(ci.lwd, length(CI))

  # remove phi' if included in res
  phiLcol <- which(names(res) == "phiL")
  if(length(phiLcol > 0)){
    x <- res[,-phiLcol]
  } else {
    x <- res
  }
  d <- ncol(x)

  # First a fitting round to look for shifts in t_anchor
  age = seq(0, agemax, 0.01)
  Lt0 <- matrix(NaN, nrow=length(age), ncol=nrow(x))
  Lt_minus <- matrix(NaN, nrow=length(age), ncol=nrow(x))
  Lt_plus <- matrix(NaN, nrow=length(age), ncol=nrow(x))
  for(i in seq(ncol(Lt0))){
    par0 <- par_minus <- par_plus <- as.list(x[i,])
    par_minus$t_anchor <- par0$t_anchor-1
    par_plus$t_anchor <- par0$t_anchor+1
    Lt0[,i] <- TropFishR::VBGF(param = par0, t = age)
    Lt_minus[,i] <- TropFishR::VBGF(param = par_minus, t = age)
    Lt_plus[,i] <- TropFishR::VBGF(param = par_plus, t = age)
  }

  # replace negative lengths with NA
  Lt0 <- replace(Lt0, Lt0<0, NaN)
  Lt_minus <- replace(Lt_minus, Lt_minus<0, NaN)
  Lt_plus <- replace(Lt_plus, Lt_plus<0, NaN)

  # determine if a positive or negative shift improves overall covariance for each permutation
  cov0 <- cov(Lt0, use = "pair")
  shift <- 0 * seq(ncol(Lt0))
  for(i in seq(ncol(Lt0))){
    cov_minus <- cov(x = Lt_minus[,i], y = Lt0, use = "pair")
    cov_plus <- cov(x = Lt_plus[,i], y = Lt0, use = "pair")
    shift[i] <- c(0,-1,1)[which.max(c(sum(cov0[i,]), sum(cov_minus), sum(cov_plus)))]
  }

  # fixed predictions
  agenew <- seq(0+min(shift), agemax+max(shift), 0.01)
  Lt <- matrix(NaN, nrow=length(agenew), ncol=nrow(x))
  for(i in seq(ncol(Lt))){
    par0 <- as.list(x[i,])
    par0$t_anchor <- par0$t_anchor + shift[i]
    Lt[,i] <- TropFishR::VBGF(param = par0, t = agenew)
  }

  for(i in seq(ncol(Lt))){
    if(i == 1) {
      plot(agenew, Lt[,i], t="n",
        xlim = c(min(x$t_anchor+shift), max(agenew))+c(-0.1,0),
        ylim = c(0,max(Lt)*1.05),
        xlab = xlab, ylab = ylab,
        xaxs = "i", yaxs = "i"
      )
    }
    lines(agenew, Lt[,i], col = perm.col, lwd = perm.lwd)
  }
  box()

  # multivariate kernel density estimate for max. density
  H <- ks::Hpi(x, nstage = 1)
  fhat <- ks::kde(x = x, H = H, eval.points = x)

  # maximum density
  max_dens <- fhat$eval.points[which.max(fhat$estimate),] # maximum density

  # predict density estimate of original data
  x$estimate <- fhat$estimate


  # determin which resamples are in the CI
  limCI <- vector(mode = "list", length(CI))
  inCI <- vector(mode = "list", length(CI))
  names(inCI) <- names(limCI) <- paste0("CI", CI)
  for(j in seq(limCI)){
    inCI[[j]] <- x$estimate > quantile(x$estimate, probs = (100-CI[j])/100 )
    limCI[[j]] <- data.frame(t = agenew, min = NaN, max = NaN)
    for(i in seq(agenew)){
      limCI[[j]]$min[i] <- min(Lt[i, which(inCI[[j]]) ], na.rm = TRUE)
      limCI[[j]]$max[i] <- max(Lt[i, which(inCI[[j]]) ], na.rm = TRUE)
    }
    lines(max ~ t, limCI[[j]], col = ci.col[j], lwd = ci.lwd[j], lty = ci.lty[j])
    lines(min ~ t, limCI[[j]], col = ci.col[j], lwd = ci.lwd[j], lty = ci.lty[j])
  }
  inCI <- as.data.frame(inCI)

  lines(
    agenew, TropFishR::VBGF(param = as.list(max_dens), t = agenew),
    col = maxd.col, lwd = maxd.lwd, lty = maxd.lty
  )

  if(add_legend){
    legend("bottomright", legend = c(paste0("CI = ", CI, "%"), "Max. Dens."),
      bty = "n", col = c(ci.col, maxd.col), lty = c(ci.lty, maxd.lty),
      lwd = c(ci.lwd, maxd.lwd)
    )
  }

  if(add_max_dens_legend){
    legend("topleft",
      legend = c(paste0(names(max_dens), " = ", sprintf("%.2f", round(max_dens,2) ))),
      bty = "n", title = "Max. Dens. \nparameters:",
      inset = c(0,0.1)
    )
  }

  RES <- list(limCI = limCI, inCI = inCI, density = x$estimate, max_dens = as.list(max_dens))
  return(RES)
}