R/SSplotMCMC_ExtraSelex.R

In r4ss/r4ss: R Code for Stock Synthesis

#' Plot uncertainty around chosen selectivity ogive from MCMC.
#'
#' Plot uncertainty in selectivity from an MCMC output for whichever fleet/year
#' was chosen in the optional extra "more stddev reporting"
#'
#'
#' @param post A data frame containing either derived_posteriors.sso or a good
#' subset of it. This can be an element of the list created by the the
#' [SSgetMCMC()] function.
#' @param add TRUE/FALSE option to add results to an existing plot.
#' @param nsexes Number of sexes in the model (should match model values but is
#' only used in the title).
#' @param shift Optional adjustment to the x values to avoid overlap of
#' intervals when overplotting on an existing plot.
#' @param fleetname Optional input to make the title better. Default will be
#' something like "Fleet 1", using the numbering from the model.
#' @param col Color for points and lines.
#' @author Ian Taylor
#' @export
SSplotMCMC_ExtraSelex <- function(
  post,
  add = FALSE,
  nsexes = 1,
  shift = 0,
  fleetname = "default",
  col = "blue"
) {
  # post is a data.frame containing either derived_posteriors.sso or a good subset of it
  #      it can be an element of the list created by the the SSgetMCMC function
  # add will add to existing plot
  # nsexes should match model values but is only used in the title
  # shift will shift the x values by for overplotting on an existing plot
  # fleetname will make title better
  # col will color points and lines

  cols <- grep("Selex_std", names(post))
  if (length(cols) == 0) {
    stop("no columns in posteriors include text 'Selex_std'")
  } else {
    sel <- post[, cols]
    names <- names(post)[cols]
    splitnames <- strsplit(names, "_")
    namesDF <- as.data.frame(matrix(
      unlist(strsplit(names, "_")),
      ncol = 6,
      byrow = T
    ))
    i <- as.numeric(as.character(namesDF[["V3"]]))[1]
    m <- as.character(namesDF[["V4"]])[1]
    agelen <- as.character(namesDF[["V5"]])[1]
    bin <- sort(unique(as.numeric(as.character(namesDF[["V6"]])))) + shift
    quants <- apply(sel, 2, quantile, probs = c(0.025, 0.5, 0.975))

    xlab <- "Age (years)"
    if (fleetname == "default") {
      fleetname <- paste("Fleet", i)
    }

    if (m == "Fem" & nsexes == 1) {
      sextitle3 <- ""
    }
    if (m == "Fem" & nsexes == 2) {
      sextitle3 <- "females"
    }
    if (m == "Mal") {
      sextitle3 <- "males"
    }
    main <- paste("Uncertainty in selectivity for", fleetname, sextitle3)
    no0 <- as.numeric(quants[3, ] - quants[1, ]) != 0
    upper <- quants[3, ]
    lower <- quants[1, ]
    if (!add) {
      matplot(
        bin,
        t(quants),
        lty = c(3, 1, 3),
        lwd = c(1, 3, 1),
        xlab = xlab,
        ylim = c(0, 1),
        main = main,
        ylab = "Selectivity",
        type = "n",
        xlim = c(0, max(bin))
      )
    }

    lines(bin, quants[2, ], lty = 1, col = col, lwd = 1, type = "o")
    arrows(
      x0 = bin[no0],
      y0 = quants[1, no0],
      x1 = bin[no0],
      y1 = quants[3, no0],
      length = 0.01,
      angle = 90,
      code = 3,
      col = col
    )
    abline(h = 0, col = "grey")
  }
}