#' Plot indices of abundance and associated quantities.
#'
#' Plot indices of abundance with or without model fit as well as other diagnostic
#' plots such as observed vs. expected index and plots related to time-varying
#' catchability (if present).
#'
#'
#' @template replist
#' @param subplots vector controlling which subplots to create
#' Numbering of subplots is as follows, where subplot 9 (comparison of all indices) is
#' provided first:
#' \itemize{
#' \item 1 index data by fleet
#' \item 2 index data with fit by fleet
#' \item 3 observed vs expected index values with smoother
#' \item 4 index data by fleet on a log scale (lognormal error only)
#' \item 5 index data with fit by fleet on a log scale (lognormal error only)
#' \item 6 log(observed) vs log(expected) with smoother (lognormal error only)
#' \item 7 time series of time-varying catchability (only if actually time-varying)
#' \item 8 catchability vs. vulnerable biomass (if catchability is not constant)
#' \item 9 comparison of all indices
#' \item 10 index residuals based on total uncertainty
#' \item 11 index residuals based on input uncertainty (not currently provided)
#' \item 12 index deviations (independent of index uncertainty)
#' }
#'
#' @template plot
#' @template print
#' @template fleets
#' @template fleetnames
#' @param smooth add smoothed line to plots of observed vs. expected sample
#' sizes
#' @param add add to existing plot (not yet implemented)
#' @param datplot make plot of data only?
#' @template labels
#' @param fleetcols vector of colors for all fleets (including those
#' with no index data)
#' @param col1 vector of colors for points in each season for time series plot.
#' Default is red for single season models and a rainbow using the
#' rich.colors.short function for multiple seasons.
#' @param col2 vector of colors for points in each season for obs. vs. exp.
#' plot. Default is blue for single season models and a rainbow using the
#' rich.colors.short function for multiple seasons.
#' @param col3 color of line showing expected index in time series plot.
#' Default is blue.
#' @param col4 color of smoother shown in obs. vs. exp. plots. Default is red.
#' @param pch1 single value or vector of plotting characters (pch parameter)
#' for time-series plots of index fit. Default=21.
#' @param pch2 single value or vector of plotting characters (pch parameter)
#' for sample size plots of index fit. Default=16.
#' @param cex character expansion factor for points showing observed values.
#' Default=1.
#' @param bg Background color for points with pch=21.
#' @param legend add a legend to seasonal colors (only for seasonal models)
#' @template legendloc
#' @param seasnames optional vector of names for each season to replace
#' defaults if a legend is used
#' @template pwidth
#' @template pheight
#' @template punits
#' @template res
#' @template ptsize
#' @template cex.main
#' @template mainTitle
#' @template plotdir
#' @param minyr First year to show in plot (for zooming in on a subset of
#' values)
#' @param maxyr Last year to show in plot (for zooming in on a subset of
#' values)
#' @param maximum_ymax_ratio Maximum allowed value for ymax (specified
#' as ratio of y), which overrides any
#' value of ymax that is greater (default = Inf)
#' @param show_input_uncertainty Switch controlling whether to add thicker
#' uncertainty interval lines indicating the input uncertainty relative to
#' the total uncertainty which may result from estimating a parameter for
#' extra standard deviations. This is only added for the plots with index
#' fit included (the data-only plots only show the input uncertainty).
#' @template verbose
#' @param \dots Extra arguments to pass to calls to `plot`
#' @author Ian Stewart, Ian Taylor, James Thorson
#' @export
#' @seealso [SS_plots()], [SS_output()]
SSplotIndices <-
function(replist,
subplots = c(1:10, 12), # IGT 2021/4/15: not sure why 11 is skipped
plot = TRUE, print = FALSE,
fleets = "all", fleetnames = "default",
smooth = TRUE, add = FALSE, datplot = TRUE,
labels = c(
"Year", # 1
"Index", # 2
"Observed index", # 3
"Expected index", # 4
"Log index", # 5
"Log observed index", # 6
"Log expected index", # 7
"Standardized index", # 8
"Catchability (Q)", # 9
"Time-varying catchability", # 10
"Vulnerable biomass", # 11
"Catchability vs. vulnerable biomass", # 12
"Residual", # 13
"Deviation"
), # 14
fleetcols = NULL,
col1 = "default", col2 = "default", col3 = "blue", col4 = "red",
pch1 = 21, pch2 = 16, cex = 1, bg = "white",
legend = TRUE, legendloc = "topright", seasnames = NULL,
pwidth = 6.5, pheight = 5.0, punits = "in", res = 300, ptsize = 10, cex.main = 1,
mainTitle = FALSE, plotdir = "default", minyr = NULL, maxyr = NULL,
maximum_ymax_ratio = Inf, show_input_uncertainty = TRUE, verbose = TRUE, ...) {
# get some quantities from replist
cpue <- replist[["cpue"]]
SS_versionNumeric <- replist[["SS_versionNumeric"]]
# confirm that some CPUE values are present
if (is.null(dim(cpue))) {
message("skipping index plots: no index data in this model")
return()
}
# table to store information on each plot
plotinfo <- NULL
# define a bunch of internal functions
index.fn <- function(addexpected = TRUE, log = FALSE, ...) {
# plot of time series of observed values with fit (if requested)
# don't do anything for log-scale plot if normal error structure is used
if (error == -1 & log == TRUE) {
return()
}
# function to get uncertainty intervals around points
# (with or without extra standard error included)
get_intervals <- function(total = TRUE) {
if (total) {
colname <- "SE"
} else {
colname <- "SE_input"
}
if (error == 0) {
if (!log) {
lower <- qlnorm(.025,
meanlog = log(y[include]),
sdlog = cpueuse[[colname]][include]
)
upper <- qlnorm(.975,
meanlog = log(y[include]),
sdlog = cpueuse[[colname]][include]
)
} else {
lower <- qnorm(.025,
mean = log(y[include]),
sd = cpueuse[[colname]][include]
)
upper <- qnorm(.975,
mean = log(y[include]),
sd = cpueuse[[colname]][include]
)
}
}
# normal error interval
if (error == -1) {
lower <- qnorm(.025, mean = y[include], sd = cpueuse[[colname]][include])
upper <- qnorm(.975, mean = y[include], sd = cpueuse[[colname]][include])
}
# T-distribution interval
if (error > 0) {
lower <- log(y[include]) + qt(.025, df = error) * cpueuse[[colname]][include]
upper <- log(y[include]) + qt(.975, df = error) * cpueuse[[colname]][include]
if (!log) {
lower <- exp(lower)
upper <- exp(upper)
}
}
return(data.frame(lower = lower, upper = upper))
}
# get total uncertainty
total_intervals <- get_intervals(total = TRUE)
lower_total <- total_intervals[["lower"]]
upper_total <- total_intervals[["upper"]]
# get input uncertainty
# SE_input column was first available in SS version 3.30.15 but is calculated
# elsewhere in this function for models that didn't report it
input_intervals <- get_intervals(total = FALSE)
lower_input <- input_intervals[["lower"]]
upper_input <- input_intervals[["upper"]]
if (max(upper_total) == Inf) {
warning(
"Removing upper interval on indices with infinite upper quantile values.\n",
"Check the uncertainty inputs for the indices."
)
upper_total[upper_total == Inf] <- 100 * max(cpueuse[["Obs"]][upper_total == Inf])
}
# plot title
main <- paste0(labels[2], Fleet)
if (log) {
main <- paste0(labels[5], Fleet)
}
# no title
if (!mainTitle) {
main <- ""
}
xlim <- c(max(minyr, min(x)), min(maxyr, max(x)))
if (legend & length(colvec1) > 1) {
xlim[2] <- xlim[2] + 0.25 * diff(xlim)
}
if (!add) {
# get range for expected values
zrange <- NULL
if (addexpected) {
zrange <- range(z, na.rm = TRUE)
}
logzrange <- range(log(z))
# y-limits for non-log plot
if (!log) {
# ylim for standard scale (if lognormal)
if (error != -1) {
ylim <- c(0, 1.05 * min(
max(upper_total, zrange, na.rm = TRUE),
max(maximum_ymax_ratio * y, na.rm = TRUE)
))
} else {
ylim <- 1.05 * c(
min(lower_total, zrange, na.rm = TRUE),
min(
max(upper_total, zrange, na.rm = TRUE),
max(maximum_ymax_ratio * y, na.rm = TRUE)
)
)
}
}
if (log) {
# ylim for log scale plot
ylim <- range(c(lower_total, upper_total), na.rm = TRUE)
}
plot(
x = x[include], y = y[include], type = "n", xlab = labels[1],
ylab = ifelse(!log, labels[2], labels[5]),
main = main, cex.main = cex.main, xlim = xlim,
ylim = ylim,
yaxs = ifelse(log, "r", "i"),
...
)
# add line at 0 if it's not a log-scale plot
# and the axes include zero
if (!log & min(ylim) < 0) {
abline(h = 0, lty = 3)
}
}
# set bounds for arrows at total uncertainty
lower <- lower_total
upper <- upper_total
if (addexpected) {
# show thicker lines behind final lines for input uncertainty
# only in plot with expected value as well
if (show_input_uncertainty & !all(lower_input == lower_total, na.rm = TRUE)) {
segments(x[include], lower_input,
x[include], upper_input,
col = colvec1[s], lwd = 3, lend = 1
)
}
} else {
# change bounds for arrows at input uncertainty for plots without fit
# if values are available
lower <- lower_input
upper <- upper_input
}
# add intervals
arrows(
x0 = x[include], y0 = lower,
x1 = x[include], y1 = upper,
length = 0.03, angle = 90, code = 3, col = colvec1[s]
)
# add points and expected values on standard scale
if (!log) {
points(
x = x[include], y = y[include],
pch = pch1, cex = cex, bg = bg, col = colvec1[s]
)
if (addexpected) {
lines(x, z, lwd = 2, col = col3)
if (length(x) == 1) {
points(x, z, pch = 23, col = col3)
}
}
} else {
# add points and expected values on log scale
points(
x = x[include], y = log(y[include]),
pch = pch1, cex = cex, bg = bg, col = colvec1[s]
)
if (addexpected) {
lines(x, log(z), lwd = 2, col = col3)
if (length(x) == 1) {
points(x, log(z), pch = 23, col = col3)
}
}
}
if (legend & length(colvec1) > 1) {
legend(x = legendloc, legend = seasnames, pch = pch1, col = colvec1, cex = cex)
}
}
index_resids.fn <- function(option = 1, ...) {
# plot of time series of residuals
# options
# 1: residuals based on total SE
# 2: residuals based on input SE
# 3: deviations (independent of index variability)
# choose y value and y-axis label
if (option == 1) { # residuals based on total SE
ylab <- labels[13]
y <- (log(cpueuse[["Obs"]]) - log(cpueuse[["Exp"]])) / cpueuse[["SE"]]
}
if (error == 0 & option == 2) { # residuals based on input SE
ylab <- labels[13]
# manually calculating residual based on SE_input
y <- (log(cpueuse[["Obs"]]) - log(cpueuse[["Exp"]])) / cpueuse[["SE_input"]]
}
if (option == 3) { # deviations
ylab <- labels[14]
# Dev should be equal to log(Obs/Exp)
y <- cpueuse[["Dev"]]
}
# plot title
main <- paste(ylab, Fleet)
# no plot title
if (!mainTitle) {
main <- ""
}
# xlim (maybe reduced by inputs minyr and maxyr)
xlim <- c(max(minyr, min(x)), min(maxyr, max(x)))
if (legend & length(colvec1) > 1) {
xlim[2] <- xlim[2] + 0.25 * diff(xlim)
}
# ylim is symetrical around 0
ylim <- c(-1.05, 1.05) * max(abs(y[include]))
if (!add) {
plot(
x = x[include], y = y[include], type = "n",
xlab = labels[1], xlim = xlim,
ylab = ylab, ylim = ylim, yaxs = "i",
main = main, cex.main = cex.main,
...
)
}
# add points
points(
x = x[include], y = y[include],
pch = pch1, cex = cex,
bg = adjustcolor(colvec1[s], alpha.f = 0.7),
col = adjustcolor(colvec1[s], alpha.f = 0.7)
)
# add line at 0
abline(h = 0, lty = 3)
# add legend if more than one color (indicating season) was used
if (legend & length(colvec1) > 1) {
legend(
x = legendloc, legend = seasnames, pch = pch1,
pt.bg = colvec1, col = colvec1, cex = cex
)
}
}
obs_vs_exp.fn <- function(log = FALSE, ...) {
# plot of observed vs. expected with smoother
# plot title
main <- paste(labels[2], Fleet, sep = " ")
# no title
if (!mainTitle) {
main <- ""
}
if (!add) {
if (!log) {
# standard plot
plot(y[include], z[include],
type = "n",
xlab = labels[3], ylab = labels[4],
main = main, cex.main = cex.main,
ylim = c(0, 1.05 * max(z)), xlim = c(0, 1.05 * max(y)),
xaxs = "i", yaxs = "i", ...
)
} else {
# log-scale plot doesn't specificy y limits
plot(log(y[include]), log(z[include]),
type = "n",
xlab = labels[6], ylab = labels[7],
main = main, cex.main = cex.main
)
}
}
if (!log) {
points(y[include], z[include], col = colvec2[s], pch = pch2, cex = cex)
} else {
points(log(y[include]), log(z[include]),
col = colvec2[s], pch = pch2, cex = cex
)
}
abline(a = 0, b = 1, lty = 3)
if (smooth && npoints > 6 && diff(range(y)) > 0) {
if (!log) {
psmooth <- loess(z[include] ~ y[include], degree = 1)
lines(psmooth[["x"]][order(psmooth[["x"]])], psmooth[["fitted"]][order(psmooth[["x"]])],
lwd = 1.2, col = col4, lty = "dashed"
)
} else {
psmooth <- loess(log(z[include]) ~ log(y[include]), degree = 1)
lines(psmooth[["x"]][order(psmooth[["x"]])], psmooth[["fitted"]][order(psmooth[["x"]])],
lwd = 1.2, col = col4, lty = "dashed"
)
}
}
if (legend & length(colvec2) > 1) {
legend(x = legendloc, legend = seasnames, pch = pch2, col = colvec2, cex = cex)
}
}
timevarying_q.fn <- function() {
# plot of time-varying catchability (if present)
main <- paste(labels[10], Fleet, sep = " ")
if (!mainTitle) main <- ""
q <- cpueuse[["Calc_Q"]]
if (!add) {
plot(x, q,
type = "o", xlab = labels[1], main = main,
cex.main = cex.main, ylab = labels[9],
col = colvec2[1], pch = pch2
)
}
}
q_vs_vuln_bio.fn <- function() {
# plot of time-varying catchability (if present)
main <- paste(labels[12], Fleet, sep = " ")
if (!mainTitle) main <- ""
v <- cpueuse[["Vuln_bio"]]
q1 <- cpueuse[["Calc_Q"]]
q2 <- cpueuse[["Eff_Q"]]
if (all(q1 == q2)) ylab <- labels[9] else ylab <- "Effective catchability"
if (!add) {
plot(v, q2,
type = "o", xlab = labels[11], main = main,
cex.main = cex.main, ylab = ylab,
col = colvec2[1], pch = pch2
)
}
}
# check for super periods
if (length(grep("supr_per", cpue[["Supr_Per"]]))) {
warning(
"Some indices have superperiods. Values will be plotted\n",
"in year/season associated with data in report file."
)
cpue <- cpue[!is.na(cpue[["Dev"]]), ]
}
FleetNames <- replist[["FleetNames"]]
nfleets <- replist[["nfleets"]]
nseasons <- replist[["nseasons"]]
# find any extra SD parameters
parameters <- replist[["parameters"]]
Q_extraSD_info <- parameters[grep("Q_extraSD", parameters[["Label"]]), ]
# calculate how many of these parameters there are
nSDpars <- nrow(Q_extraSD_info)
if (nSDpars > 0) {
# parse the parameter label to get the fleet number
Q_extraSD_info[["Fleet"]] <- NA
for (ipar in 1:nSDpars) {
if (SS_versionNumeric >= 3.3) {
# parsing label with ending like "(2)" assuming only one set of parentheses
num <- strsplit(Q_extraSD_info[["Label"]][ipar], split = "[()]", fixed = FALSE)[[1]][2]
} else {
num <- strsplit(substring(Q_extraSD_info[["Label"]][ipar], nchar("Q_extraSD_") + 1),
split = "_", fixed = TRUE
)[[1]][1]
}
Q_extraSD_info[["Fleet"]][ipar] <- as.numeric(num)
}
# NOTE: important columns in Q_extraSD_info to use below are $Value and $Fleet
}
if (nseasons > 1) {
# if seasons, put CPUE at season midpoint
cpue[["YrSeas"]] <- cpue[["Yr"]] + (cpue[["Seas"]] - 0.5) / nseasons
} else {
# if no seasons, put at integer year value
cpue[["YrSeas"]] <- cpue[["Yr"]]
}
if (plotdir == "default") plotdir <- replist[["inputs"]][["dir"]]
if (fleetnames[1] == "default") {
fleetnames <- FleetNames
}
if (fleets[1] == "all") {
fleets <- 1:nfleets
} else {
if (length(intersect(fleets, 1:nfleets)) != length(fleets)) {
return("Input 'fleets' should be 'all' or a vector of values between 1 and nfleets.")
}
}
# subset fleets as requested
fleetvec <- intersect(fleets, unique(as.numeric(cpue[["Fleet"]])))
# empty data.frame to store data for comparison among indices
allcpue <- data.frame()
# keep track of whether any indices with negative observations is excluded
any_negative <- FALSE
# loop over fleets
for (ifleet in fleetvec) {
# use fancy colors only if the individual index spans more than one season
usecol <- FALSE
if (length(unique(cpue[["Seas"]][cpue[["Fleet"]] == ifleet])) > 1) {
usecol <- TRUE
}
# turn off use of legend if there's never more than 1 season per index
if (!usecol) {
legend <- FALSE
}
if (col1[1] == "default") {
colvec1 <- "black"
if (usecol & nseasons == 4) {
colvec1 <- c("blue4", "green3", "orange2", "red3")
}
if (usecol & !nseasons %in% c(1, 4)) {
colvec1 <- rich.colors.short(nseasons)
}
} else {
colvec1 <- col1
# if user provides single value (or vector of length less than nseasons)
# make sure it's adequate to cover all seasons
if (length(colvec1) < nseasons) {
colvec1 <- rep(col1, nseasons)
}
}
if (col2[1] == "default") {
colvec2 <- "blue"
if (usecol & nseasons == 4) {
colvec2 <- c("blue4", "green3", "orange2", "red3")
}
if (usecol & !nseasons %in% c(1, 4)) {
colvec2 <- rich.colors.short(nseasons)
}
} else {
colvec2 <- col2
# if user provides single value (or vector of length less than nseasons)
# make sure it's adequate to cover all seasons
if (length(colvec1) < nseasons) {
colvec1 <- rep(col1, nseasons)
}
}
if (is.null(seasnames)) seasnames <- paste("Season", 1:nseasons, sep = "")
Fleet <- fleetnames[ifleet]
error <- replist[["survey_error"]][ifleet]
if (error == 0) {
error_caption <- "lognormal error"
}
if (error == -1) {
error_caption <- "normal error"
}
if (error == 1) {
error_caption <- paste0(
"T-distributed error with ", error,
" degree of freedom"
)
}
if (error > 1) {
error_caption <- paste0(
"T-distributed error with ", error,
" degrees of freedom"
)
}
cpueuse <- cpue[cpue[["Fleet"]] == ifleet, ]
cpueuse <- cpueuse[order(cpueuse[["YrSeas"]]), ]
# look for time-vary
time <- diff(range(cpueuse[["Calc_Q"]])) > 0
# look for time-varying effective Q
time2 <- diff(range(cpueuse[["Eff_Q"]])) > 0
# Teresa's model had NA values in Eff_Q for unknown reasons
# line below will allow model to play on
if (is.na(time2)) {
time2 <- FALSE
}
# if "SE_input" column not available, look for extra SD and
# calculate input SD (if different from final value)
if (!"SE_input" %in% names(cpue)) {
if (exists("Q_extraSD_info") && ifleet %in% Q_extraSD_info[["Fleet"]]) {
# input uncertainty is final value minus extra SD parameter (if present)
cpueuse[["SE_input"]] <- cpueuse[["SE"]] - Q_extraSD_info[["Value"]][Q_extraSD_info[["Fleet"]] == ifleet]
} else {
cpueuse[["SE_input"]] <- cpueuse[["SE"]]
}
}
# use short variable names for often-used quantities
x <- cpueuse[["YrSeas"]]
y <- cpueuse[["Obs"]]
z <- cpueuse[["Exp"]]
npoints <- length(z)
include <- !is.na(cpueuse[["Like"]])
if (any(include)) {
if (usecol) {
s <- cpueuse[["Seas"]][which(include)]
} else {
s <- 1 # only use colorvector if more than 1 season
}
if (datplot) {
# add index data to data frame which is used to compare all indices
if (min(cpueuse[["Obs"]] >= 0)) {
cpueuse[["Index"]] <- rep(ifleet, length(cpueuse[["YrSeas"]]))
cpueuse[["stdvalue"]] <- cpueuse[["Obs"]] / mean(cpueuse[["Obs"]])
tempcpue <- cbind(cpueuse[["Index"]], cpueuse[["YrSeas"]], cpueuse[["Obs"]], cpueuse[["stdvalue"]])
colnames(tempcpue) <- c("Index", "year", "value", "stdvalue")
allcpue <- rbind(allcpue, tempcpue)
} else {
if (verbose & 9 %in% subplots & datplot) {
message(
"Excluding fleet ", ifleet,
" from index comparison figure because it has negative values"
)
}
any_negative <- TRUE
}
}
addlegend <- function(pch, colvec) {
names <- paste(seasnames, "observations")
}
if (plot) {
if (1 %in% subplots & datplot) {
index.fn(addexpected = FALSE)
}
if (2 %in% subplots) {
index.fn()
}
if (3 %in% subplots) {
obs_vs_exp.fn()
}
}
if (print) {
if (1 %in% subplots & datplot) {
file <- paste0("index1_cpuedata_", Fleet, ".png")
caption <- paste0(
"Index data for ", Fleet, ". ",
"Lines indicate 95% uncertainty interval around index values ",
"based on the model assumption of ", error_caption, ". ",
"Thicker lines (if present) indicate input uncertainty before addition of ",
"estimated additional uncertainty parameter."
)
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
index.fn(addexpected = FALSE)
dev.off()
}
if (2 %in% subplots) {
file <- paste0("index2_cpuefit_", Fleet, ".png")
caption <- paste0(
"Fit to index data for ", Fleet, ". ",
"Lines indicate 95% uncertainty interval around index values ",
"based on the model assumption of ", error_caption, ". ",
"Thicker lines (if present) indicate input uncertainty before addition of ",
"estimated additional uncertainty parameter."
)
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
index.fn()
dev.off()
}
if (3 %in% subplots) {
file <- paste0("index3_obs_vs_exp_", Fleet, ".png")
caption <- paste("Observed vs. expected index values with smoother for", Fleet)
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
obs_vs_exp.fn()
dev.off()
}
}
# same plots again in log space
# check for lognormal error
if (error != -1) {
# plot subplots 4-6 to current device
if (plot) {
if (4 %in% subplots & datplot) {
index.fn(log = TRUE, addexpected = FALSE)
}
if (5 %in% subplots) {
index.fn(log = TRUE)
}
if (6 %in% subplots) {
obs_vs_exp.fn(log = TRUE)
}
}
# print subplots 4-6 to PNG files
if (print) {
if (4 %in% subplots & datplot) {
file <- paste0("index4_logcpuedata_", Fleet, ".png")
caption <- paste0(
"Log index data for ", Fleet, ". ",
"Lines indicate 95% uncertainty interval around index values ",
"based on the model assumption of ", error_caption, ". ",
"Thicker lines (if present) indicate input uncertainty before addition of ",
"estimated additional uncertainty parameter."
)
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
index.fn(log = TRUE, addexpected = FALSE)
dev.off()
}
if (5 %in% subplots) {
file <- paste0("index5_logcpuefit_", Fleet, ".png")
caption <- paste0(
"Fit to log index data on log scale for ", Fleet, ". ",
"Lines indicate 95% uncertainty interval around index values ",
"based on the model assumption of ", error_caption, ". ",
"Thicker lines (if present) indicate input uncertainty before addition of ",
"estimated additional uncertainty parameter."
)
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
index.fn(log = TRUE)
dev.off()
}
if (6 %in% subplots) {
file <- paste0("index6_log_obs_vs_exp_", Fleet, ".png")
caption <- paste("log(observed) vs. log(expected) index values with smoother for", Fleet)
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
obs_vs_exp.fn(log = TRUE)
dev.off()
}
}
} # end plots that require lognormal error
# plots 7 and 8 related to time-varying catchability
if (plot) {
if (7 %in% subplots & time) {
timevarying_q.fn()
}
if (8 %in% subplots & time2) {
q_vs_vuln_bio.fn()
}
} # end plot to graphics device
if (print) {
if (7 %in% subplots & time) {
file <- paste0("index7_timevarying_q_", Fleet, ".png")
caption <- paste("Timeseries of catchability for", Fleet)
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
timevarying_q.fn()
dev.off()
}
if (8 %in% subplots & time2) {
file <- paste0("index8_q_vs_vuln_bio_", Fleet, ".png")
caption <-
paste0(
"Catchability vs. vulnerable biomass for fleet ", Fleet, "<br> \n",
"This plot should illustrate curvature of nonlinear catchability relationship<br> \n",
"or reveal patterns associated with random-walk catchability."
)
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
q_vs_vuln_bio.fn()
dev.off()
}
} # end print to PNG
# residual/deviation plots
if (plot) {
if (10 %in% subplots & all(cpueuse[["Obs"]] >= 0)) {
index_resids.fn(option = 1)
}
if (11 %in% subplots) {
index_resids.fn(option = 2)
}
if (12 %in% subplots) {
index_resids.fn(option = 3)
}
}
if (print) {
#### residuals based on total uncertainty
if (10 %in% subplots & all(cpueuse[["Obs"]] >= 0)) {
file <- paste0("index10_resids_SE_total_", Fleet, ".png")
caption <- paste0("Residuals of fit to index for ", Fleet, ".")
if (error == 0) {
caption <- paste0(
caption,
"<br>Values are (log(Obs) - log(Exp))/SE ",
"where SE is the total standard error including any ",
"estimated additional uncertainty."
)
} else {
caption <- paste0(
caption,
"<br>Values are based on the total standard error ",
"including any estimated additional uncertainty."
)
}
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
index_resids.fn(option = 1)
dev.off()
}
#### residuals based on input uncertainty
if (11 %in% subplots &
show_input_uncertainty &&
any(!is.null(cpueuse[["SE_input"]][include])) &&
any(cpueuse[["SE_input"]] > cpueuse[["SE"]])) {
file <- paste0("index11_resids_SE_input_", Fleet, ".png")
caption <- paste0("Residuals for fit to index for ", Fleet, ".")
if (error == 0) {
caption <- paste0(
caption,
"<br>Values are (log(Obs) - log(Exp))/SE_input ",
"where SE_input is the input standard error",
"excluding any estimated additional uncertainty."
)
} else {
caption <- paste0(
caption,
"<br>Values are based on the input standard error ",
"excluding any estimated additional uncertainty."
)
}
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
index_resids.fn(option = 2)
dev.off()
}
#### simple deviation plot
if (12 %in% subplots) {
file <- paste0("index12_resids_SE_total_", Fleet, ".png")
caption <- paste0("Deviations for fit to index for ", Fleet, ".")
if (error != -1) {
# lognormal or T-distributed error
caption <- paste0(
caption,
"<br>Values are log(Obs) - log(Exp) ",
"and thus independent of index uncertainty."
)
}
if (error == -1) {
# normal error
caption <- paste0(
caption,
"<br>Values are Obs - Exp ",
"and thus independent of index uncertainty."
)
}
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
index_resids.fn(option = 3)
dev.off()
}
} # end if(print)
} # end check for any values to include
} # end loop over fleets
### standardized plot of all CPUE indices
if (datplot == TRUE & nrow(allcpue) > 0) {
all_index.fn <- function() {
main <- "All index plot"
if (!mainTitle) {
main <- ""
}
xlim <- c(
min(allcpue[["year"]], na.rm = TRUE) - 1,
max(allcpue[["year"]], na.rm = TRUE) + 1
)
# change year range if requested
xlim[1] <- max(xlim[1], minyr)
xlim[2] <- min(xlim[2], maxyr)
# set y limits
ylim <- c(0, 1.05 * max(allcpue[["stdvalue"]], na.rm = TRUE))
# set colors
if (!is.null(fleetcols) & length(fleetcols) >= nfleets) {
usecols <- fleetcols
} else {
usecols <- rich.colors.short(max(allcpue[["Index"]], na.rm = TRUE), alpha = 0.7)
if (max(allcpue[["Index"]], na.rm = TRUE) >= 2) {
usecols <- rich.colors.short(max(allcpue[["Index"]], na.rm = TRUE) + 1,
alpha = 0.7
)[-1]
}
}
# make empty plot
if (!add) {
plot(0,
type = "n", xlab = labels[1], main = main, cex.main = cex.main,
col = usecols[1], ylab = labels[8], xlim = xlim, ylim = ylim,
yaxs = "i"
)
}
# add points and lines for each fleet
for (ifleet in fleetvec) {
lines(
x = allcpue[["year"]][allcpue[["Index"]] == ifleet],
y = allcpue[["stdvalue"]][allcpue[["Index"]] == ifleet],
col = adjustcolor(usecols[ifleet], alpha.f = 0.7),
lwd = 2
)
points(
x = allcpue[["year"]][allcpue[["Index"]] == ifleet],
y = allcpue[["stdvalue"]][allcpue[["Index"]] == ifleet],
pch = pch1,
bg = adjustcolor(usecols[ifleet], alpha.f = 0.7),
col = gray(0, alpha = 0.7),
cex = cex
)
}
legend(legendloc,
legend = fleetnames[fleetvec],
ncol = 2,
bty = "n",
pch = pch1,
col = gray(0, alpha = 0.7),
pt.bg = usecols[fleetvec]
)
} # end all_index.fn
if (plot & (9 %in% subplots)) {
all_index.fn()
}
if (print & (9 %in% subplots)) {
file <- paste0("index9_standcpueall", ".png")
caption <- paste(
"Standardized indices overlaid.",
"Each index is rescaled to have mean observation = 1.0."
)
if (any_negative) {
caption <- paste(
caption,
"Indices with negative observations have been excluded."
)
}
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
all_index.fn()
dev.off()
}
} # end datplot
if (!is.null(plotinfo)) plotinfo[["category"]] <- "Index"
return(invisible(plotinfo))
} # end function
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.