R/StockRec.plots.R

Defines functions StockRec.plots

Documented in StockRec.plots

#' Stock-recruitment plots
#'
#' The function \code{StockRec.plots} generates plots of stock vs. recruitment in
#' linear and logarithmic scales.
#'
#' @param x an R list with output from the assessment models.
#' @param DataName string used in plot titles.  Defaults to argument \code{x}.
#' @param draft modifies plots for use in a report.  When \code{FALSE} main titles
#' are omitted.
#' @param graphics.type a vector of graphics file types to which graphics are saved.
#' When \code{NULL}, no plots are saved.
#' @param use.color plots are made in grayscale when \code{FALSE}
#' @param start.drop Number of years at the start of the data to be omitted from
#' plots, as when a model includes an initialization period.
#' @param units.ssb A text string (e.g. \code{"tons"} or \code{"10^9 eggs"}) used
#' in labeling plots of spawning stock.
#' @param units.rec A text string (e.g. \code{"million fish"}) for labeling plots
#' of recruitment
#' @param rec.model Specifies the type of recruitment function to draw.  Valid
#' values are: \code{"BH", "BH-steep", "Ricker", "Ricker-steep", "Mean"}.
#' @param draw.model If \code{TRUE}, a function curve is drawn on plots of stock
#' and recruitment.
#' @param draw.lowess If \code{TRUE}, a lowess smooth is drawn on plots of stock
#' and recruitment.
#' @param draw.time If \code{TRUE}, stock-recruitment points are connected to indicate
#' progression of time.
#' @param year.pos An integer (= 1, 2, 3, or 4) defining the position of text relative to points.
#' The text indicates first and last years in the time series, and applies only if draw.time=TRUE. A
#' value of year.pos=0 turns off the text feature.
#'
#' @return Graphics
#'
#' @author M Prager
#' @author E Williams
#' @author K Shertzer
#' @author R Cheshire
#' @author K Purcell
#'
#'
#' @examples \donttest{
#' StockRec.plots(gag)
#' }
#' @export
StockRec.plots <- function(x, DataName = deparse(substitute(x)), draft = TRUE,
    graphics.type = NULL, use.color = TRUE, start.drop = 0,
    units.ssb = x$info$units.ssb, units.rec = x$info$units.rec,
    rec.model = x$info$rec.model, draw.model = TRUE, draw.lowess = FALSE,
    draw.time = TRUE, year.pos=1)
#=================================================================================
{   ### Check for required data:
    if (! ("t.series" %in% names(x)))
    {   Errmsg <- paste("Data frame 't.series' not found in data object:",
            deparse(substitute(x)))
        warning(Errmsg, immediate. = TRUE)
        return(invisible(-1))
    }
    if (! ("parms" %in% names(x)))
    {   Errmsg <- paste("List 'parms' not found in data object:", deparse(substitute(x)))
        warning(Errmsg, immediate. = TRUE)
        return(invisible(-1))
    }
    if (! ("rec.lag" %in% names(x$parms)))
    {   Errmsg <- paste("Value 'parms$rec.lag' not found in data object:",
            deparse(substitute(x)))
        warning(Errmsg, immediate. = TRUE)
        return(invisible(-1))
    }
#=================================================================================
    ### Get data, limits, titles, labels, etc.
#=================================================================================
    ts <- x$t.series
    rec.lag <- x$parms$rec.lag
    draw.bias <- FALSE

    ### Get starting and stopping indices for S and R in S-R plots
    sndx = as.integer(rep(0,2))  # initialize
    sndx[1] <- start.drop + 1
    sndx[2] <- nrow(ts) - rec.lag
    if (sndx[1] >= sndx[2])
    {   Errmsg <- "Length of S-R data series is <= 1"
        warning(Errmsg, immediate. = TRUE)
        return(invisible(-1))
    }
    # Expand index variables to full-length
    sndx <- sndx[1]:sndx[2]
    rndx <- sndx + x$parms$rec.lag

    ### Set up plotting-related stuff:
    plot.options = FGGetOptions()
    PlotTitle <- ""
    savepar <- FGSetPar(draft)
    if (! is.null(graphics.type))
    {   write.graphs <- TRUE
        GraphicsDirName <- paste(DataName, "-figs/SR", sep="")
    }
    else write.graphs <- FALSE

    # Get colors for graphical objects:
    if (use.color) parlist <- plot.options$color
    else parlist <- plot.options$bw
    clr.points <- parlist$clr.line
    clr.lowess <- parlist$clr.lightline

    ### Set plotting limits, etc., for first plot:
    ssb.max <- max(ts$SSB[sndx], na.rm = TRUE)
    r.max   <- max(ts$recruits[rndx], na.rm = TRUE)
    lab.x   <- FGMakeLabel("Spawning stock", units.ssb)
    lab.y   <- FGMakeLabel("Recruitment", units.rec)
    rec.sim <- as.numeric(NA)   # Needs a definition for use in log plot
    lim.x = range(0, 1.1 * ssb.max)
    lim.y = range(0, 1.1 * r.max)
#=================================================================================
    ### If recruitment curve will be drawn, generate simulated
    ### data & adjust Y limits of plot:
#=================================================================================
    if (is.null(rec.model))
    {   draw.model <- FALSE
        rec.model <- "none"
        curve.OK <- TRUE
    }
    # Make sure curve type is recognized, and if so, get parameters:
    if (draw.model)
    {   stock.sim <- seq(from = 0.001 * ssb.max, to = 1.1 * ssb.max, length = 200)
        curve.OK <- FALSE
        # Beverton-Holt curve (standard)
        if (rec.model == "BH")
        {   alpha <- x$parms$BH.alpha
            beta <-  x$parms$BH.beta
            bias.corr <- x$parms$BH.biascorr
            rec.sim <- alpha * stock.sim / (1.0 + beta * stock.sim)
            curve.OK <- TRUE
        }
        # Beverton-Holt curve (steepness)
        if (rec.model == "BH-steep")
        {   R0 <- x$parms$BH.R0
            h <- x$parms$BH.steep
            phi0 <- x$parms$BH.Phi0
            bias.corr <- x$parms$BH.biascorr
            rec.sim <- (0.8 * R0 * h * stock.sim) / (0.2 * phi0 * R0 * (1.0 - h) +
                (h - 0.2) * stock.sim)
            curve.OK <- TRUE
        }
        # Ricker curve (standard)
        if (rec.model == "Ricker")
        {   alpha <- x$parms$Ricker.alpha
            beta <-  x$parms$Ricker.beta
            bias.corr <- x$parms$Ricker.biascorr
            rec.sim <- alpha * stock.sim * exp(-beta * stock.sim)
            curve.OK <- TRUE
        }
        # Ricker curve (steepness)
        if (rec.model == "Ricker-steep")
        {   R0 <- x$parms$Ricker.R0
            h <- x$parms$Ricker.steep
            phi0 <- x$parms$Ricker.Phi0
            bias.corr <- x$parms$Ricker.biascorr
            rec.sim <-  stock.sim / phi0 * exp(h * (1.0 - stock.sim / (phi0 * R0)))
            curve.OK <- TRUE
        }
        # Average recruitment model
        if (rec.model == "Mean")
        {   R0 <- x$parms$Mean.R0
        bias.corr <- x$parms$Mean.biascorr
        rec.sim <-  rep(R0, length=length(stock.sim))
        curve.OK <- TRUE
        }
        if (! curve.OK)
        {   warning("Bad value of argument rec.model\n")
            draw.model <- FALSE
            rec.model <- "none"
            bias.corr <- NULL
        }
    }
    if (draw.model)
    {   if (is.numeric(bias.corr) && bias.corr != 1.0) draw.bias <- TRUE
        # Adjust plotting limits to include the curves:
        {   if (draw.bias) lim.y <- range(0, max(bias.corr * rec.sim, ts$recruits[rndx],
                na.rm = TRUE))
            else lim.y <- range(0, max(rec.sim, ts$recruits[rndx], na.rm = TRUE))
        }
    }
#=================================================================================
    ### Now make the plots:
#=================================================================================
    ### Plot of stock vs. recruitment -- arithmetic scale
    if(draft) PlotTitle <- FGMakeTitle("Stock-recruitment", DataName)
    FGTimePlot(x = ts$SSB[sndx], y = ts$recruits[rndx], y2 = NULL,
        lab.x = lab.x, lab.y = lab.y, FGtype = "circles", main = PlotTitle,
        use.color = use.color, xlim = lim.x, ylim = lim.y)
if (draw.time){
    lines(ts$SSB[sndx], ts$recruits[rndx], col = clr.points, lty=3)
    if (year.pos>0){
        text (x=ts$SSB[sndx[1]], y=ts$recruits[rndx[1]],
              labels=ts$year[rndx[1]], pos=year.pos, cex=0.85, offset=0.35)
        text (x=ts$SSB[sndx[length(sndx)]], y=ts$recruits[rndx[length(rndx)]],
              labels=ts$year[rndx[length(rndx)]], pos=year.pos, cex=0.85, offset=0.35)
    }
}
redraw <- FALSE
    if(draw.lowess)
    {   lines(lowess(ts$SSB[sndx], ts$recruits[rndx],f = 0.55),
            col = clr.lowess, lwd = 2)
        redraw <- TRUE
    }
    if (draw.model && curve.OK)
    {   lines(stock.sim, rec.sim, lwd = 2)
        if (draw.bias)
        {   lines(stock.sim, bias.corr * rec.sim, lwd = 2, lty = "dashed")
            leg.text <- c(rec.model, "Expected")
            legend("topright", legend = leg.text, lty = c("solid", "dashed"),
                lwd = 2, inset = 0.01, bg = "white")
            redraw <- TRUE
        }
    }
    # Replot to reveal any points under legend:
    if (redraw) {points(ts$SSB[sndx], ts$recruits[rndx], col = clr.points)}

    # Save plots to file:
    if (write.graphs) FGSavePlot(GraphicsDirName, DataName,
        GraphName = paste("SR.",rec.model, sep = ""), graphics.type)
    #==================================================================
#     ##### Plot of stock & recruitment with function -- log scale
#     if(draft) PlotTitle <- FGMakeTitle("Stock-recruitment (log R)", DataName)
#     {   if (draw.bias)
#             lim.y <- range(ts$recruits[rndx], bias.corr * rec.sim, na.rm = TRUE)
#         else
#             lim.y <- range(ts$recruits[rndx], rec.sim, na.rm = TRUE)
#     }
#     lim.y[1] = max(lim.y[1], min(ts$recruits[rndx], na.rm = TRUE) / 5.0)
#     FGTimePlot(x = ts$SSB[sndx], y = ts$recruits[rndx],
#         y2 = NULL, lab.x = lab.x, lab.y = lab.y, FGtype = "circles",
#         main = PlotTitle, use.color = use.color, xlim = lim.x,  ylim = lim.y,
#         log = "y")
#     if(draw.lowess)
#     {   lines(lowess(ts$SSB[sndx], ts$recruits[rndx],f = 0.60),
#             col = clr.lowess, lwd = 2)
#     }
#     if (draw.model && curve.OK)
#     {   lines(stock.sim, rec.sim, lwd = 2)
#         if (draw.bias)
#         {   lines(stock.sim, bias.corr * rec.sim, lwd = 2, lty = "dashed")
#             legend("bottomright", legend = leg.text, lty = c("solid", "dashed"),
#                 lwd = 2, inset = 0.01, bg = "white")
#         }
#     }
#     # Replot to reveal any points under legend:
#     if (redraw) points(ts$SSB[sndx], ts$recruits[rndx], col = clr.points)
#     # Save plots to file:
#     if (write.graphs) FGSavePlot(GraphicsDirName, DataName,
#         GraphName = paste("SR.", rec.model, ".log", sep = ""), graphics.type)
#=================================================================================
### Plot of stock vs. log(recruitment/spawners)
if(draft) PlotTitle <- FGMakeTitle("Stock v log(R/S)", DataName)
logRS=log(ts$recruits[rndx]/ts$SSB[sndx])
lab.y2   <- "log(recruits/spawner)"

if (draw.model && curve.OK) {
    logRS.sim=log(rec.sim/stock.sim);
    lim.y <- range(logRS, logRS.sim, na.rm = TRUE)
} else {lim.y <- range(logRS, na.rm = TRUE)}

if (any(logRS <= 0))
{   Errmsg <- "Warning: attempted to take log of a non-positive R/S"
    warning(Errmsg, immediate. = TRUE)
    return(invisible(-1))
}


lim.y[1] = 0.9*lim.y[1]; lim.y[2] = 1*lim.y[2]

FGTimePlot(x = ts$SSB[sndx], y = logRS, y2 = NULL,
           lab.x = lab.x, lab.y = lab.y2, FGtype = "circles", main = PlotTitle,
           use.color = use.color, xlim = lim.x, ylim = lim.y)
if (draw.time){
  lines(ts$SSB[sndx], logRS, col = clr.points, lty=3)
  if (year.pos>0){
    text (x=ts$SSB[sndx[1]], y=logRS[1],
          labels=ts$year[rndx[1]], pos=year.pos, cex=0.85, offset=0.35)
    text (x=ts$SSB[sndx[length(sndx)]], y=logRS[length(logRS)],
          labels=ts$year[rndx[length(rndx)]], pos=year.pos, cex=0.85, offset=0.35)
  }
}

redraw <- FALSE
if(draw.lowess)
{   lines(lowess(ts$SSB[sndx], logRS,f = 0.55),
          col = clr.lowess, lwd = 2)
    redraw = TRUE
}

if (draw.model && curve.OK)
  {lines(stock.sim, logRS.sim, lwd = 2)}
# Replot to reveal any points under legend:
if (redraw) {points(ts$SSB[sndx], logRS, col = clr.points)}
# Save plots to file:
if (write.graphs) FGSavePlot(GraphicsDirName, DataName,
                             GraphName = paste("SlogRS.",rec.model, sep = ""), graphics.type)


#==================================================================================
   par(savepar)    # reset graphics device
   return(invisible(0))
} # end function definition
rcheshire/FishGraph documentation built on Feb. 23, 2024, 11:27 a.m.