#' Make plot of likelihood contributions by fleet
#'
#' This style of plot was officially named a "Piner Plot" at the
#' CAPAM Selectivity Workshop, La Jolla March 2013. This is in
#' honor of Kevin Piner's contributions to interpreting likelihood
#' profiles. He's surely not the first person to make such a plot
#' but the name seems to have stuck.
#' @param summaryoutput List created by the function
#' [SSsummarize()].
#' @template plot
#' @template print
#' @param component Which likelihood component to plot. Default is "Length_like".
#' @param main Title for plot. Should match component.
#' @param models Optional subset of the models described in
#' `summaryoutput`. Either "all" or a vector of numbers indicating
#' columns in summary tables.
#' @template fleets
#' @template fleetnames
#' @param profile.string Character string used to find parameter over which the
#' profile was conducted. If `exact=FALSE`, this can be a substring of
#' one of the SS parameter labels found in the Report.sso file.
#' For instance, the default input 'R0'
#' matches the parameter 'SR_LN(R0)'. If `exact=TRUE`, then
#' profile.string needs to be an exact match to the parameter label.
#' @param profile.label Label for x-axis describing the parameter over which
#' the profile was conducted.
#' @param exact Should the `profile.string` have to match the parameter
#' label exactly, or is a substring OK.
#' @param ylab Label for y-axis. Default is "Change in -log-likelihood".
#' @param col Optional vector of colors for each line.
#' @param pch Optional vector of plot characters for the points.
#' @param lty Line total for the likelihood components.
#' @param lty.total Line type for the total likelihood.
#' @template lwd
#' @param lwd.total Line width for the total likelihood.
#' @param cex Character expansion for the points representing the likelihood
#' components.
#' @param cex.total Character expansion for the points representing the total
#' likelihood.
#' @param xlim Range for x-axis. Change in likelihood is calculated relative to
#' values within this range.
#' @param ymax Maximum y-value. Default is 10\% greater than largest value
#' plotted.
#' @param xaxs The style of axis interval calculation to be used for the x-axis
#' (see ?par for more info)
#' @param yaxs The style of axis interval calculation to be used for the y-axis
#' (see ?par for more info).
#' @param type Line type (see ?plot for more info).
#' @template legend
#' @template legendloc
#' @template pwidth
#' @template pheight
#' @template punits
#' @template res
#' @template ptsize
#' @template cex.main
#' @template plotdir
#' @param add_cutoff Add dashed line at ~1.92 to indicate 95% confidence interval
#' based on common cutoff of half of chi-squared of p=.95 with 1 degree of
#' freedom: `0.5*qchisq(p=cutoff_prob, df=1)`. The probability value
#' can be adjusted using the `cutoff_prob` below.
#' @param cutoff_prob Probability associated with `add_cutoff` above.
#' @template verbose
#' @param fleetgroups Optional character vector, with length equal to
#' the number of declared fleets, where fleets with the same value are
#' aggregated
#' @param likelihood_type choice of "raw" or "raw_times_lambda" (the default)
#' determines whether or not likelihoods plotted are adjusted by lambdas
#' (likelihood weights)
#' @param minfraction Minimum change in likelihood (over range considered) as a
#' fraction of change in total likelihood for a component to be included in the
#' figure.
#' @references Kevin Piner says that he's not the originator of this idea so
#' Athol Whitten is going to add a reference here.
#' @author Ian G. Taylor, Kevin R. Piner, James T. Thorson
#' @export
#' @family profile functions
PinerPlot <-
function(summaryoutput,
plot = TRUE, print = FALSE,
component = "Length_like",
main = "Changes in length-composition likelihoods by fleet",
models = "all",
fleets = "all",
fleetnames = "default",
profile.string = "R0",
profile.label = expression(log(italic(R)[0])),
exact = FALSE,
ylab = "Change in -log-likelihood",
col = "default",
pch = "default",
lty = 1, lty.total = 1,
lwd = 2, lwd.total = 3,
cex = 1, cex.total = 1.5,
xlim = "default",
ymax = "default",
xaxs = "r", yaxs = "r",
type = "o",
legend = TRUE, legendloc = "topright",
pwidth = 6.5, pheight = 5.0, punits = "in", res = 300, ptsize = 10, cex.main = 1,
plotdir = NULL,
add_cutoff = FALSE,
cutoff_prob = 0.95,
verbose = TRUE,
fleetgroups = NULL,
likelihood_type = "raw_times_lambda",
minfraction = 0.01) {
# this function is very similar to SSplotProfile, but shows fleet-specific likelihoods
# for a single components rather than multiple components aggregated across fleets
if (print & is.null(plotdir)) {
stop("to print PNG files, you must supply a directory as 'plotdir'")
}
# get stuff from summary output into shorter variable names
n <- summaryoutput[["n"]]
lbf <- summaryoutput[["likelihoods_by_fleet"]]
lbtg <- summaryoutput[["likelihoods_by_tag_group"]]
if (is.null(lbf)) {
stop(
"Input 'summaryoutput' needs to be a list output from SSsummarize\n",
"and have an element named 'likelihoods_by_fleet'."
)
}
# count of fleets
nfleets <- ncol(lbf) - 3
pars <- summaryoutput[["pars"]]
# names of fleets
FleetNames <- summaryoutput[["FleetNames"]][[1]]
# stop if lengths don't match
if (length(FleetNames) != nfleets) {
stop(
"problem with FleetNames: length!= ", nfleets, "\n",
paste(FleetNames, collapse = "\n")
)
}
# stop if component input isn't found in table
component_options <- c(
unique(lbf[["Label"]][-grep("_lambda", lbf[["Label"]])]),
unique(lbtg[["Label"]][-grep("_lambda", lbtg[["Label"]])])
)
if (!component %in% component_options) {
stop(
"input 'component' needs to be one of the following\n",
paste(" ", component_options, "\n")
)
}
if (fleetnames[1] == "default") fleetnames <- FleetNames # note lower-case value is the one used below (either equal to vector from replist, or input by user)
# check number of models to be plotted
if (models[1] == "all") {
models <- 1:n
} else {
if (!all(models %in% 1:n)) {
stop("Input 'models' should be a vector of values from 1 to n=", n, " (for your inputs).\n")
}
}
# check number of fleets to be plotted
if (fleets[1] == "all") {
fleets <- 1:nfleets
} else {
if (!all(fleets %in% 1:nfleets)) {
stop("Input 'fleets' should be a vector of values from 1 to nfleets=", nfleets, " (for your inputs).\n")
}
}
# find the parameter that the profile was over
if (exact) {
parnumber <- match(profile.string, pars[["Label"]])
} else {
parnumber <- grep(profile.string, pars[["Label"]])
}
if (length(parnumber) <= 0) {
stop("No parameters matching profile.string='", profile.string, "'", sep = "")
}
parlabel <- pars[["Label"]][parnumber]
if (length(parlabel) > 1) {
stop("Multiple parameters matching profile.string='", profile.string, "':\n",
paste(parlabel, collapse = ", "),
"\nYou may need to use 'exact=TRUE'.",
sep = ""
)
}
parvec <- as.numeric(pars[pars[["Label"]] == parlabel, models])
message(
"Parameter matching profile.string = '", profile.string, "': '",
parlabel,
"\nParameter values (after subsetting based on input 'models'): ",
paste0(parvec, collase = ", ")
)
if (xlim[1] == "default") xlim <- range(parvec)
# rearange likelihoods to be in columns by type
if (likelihood_type == "raw") prof.table <- lbf[which(lbf[["model"]] %in% models & lbf[["Label"]] == component), ]
if (likelihood_type == "raw_times_lambda") {
prof.table <- lbf[which(lbf[["model"]] %in% models & lbf[["Label"]] == component), ]
prof.table[, -c(1:3)] <- prof.table[, -c(1:3)] * lbf[which(lbf[["model"]] %in% models & lbf[["Label"]] == component) - 1, ][, -c(1:3)]
}
# Aggregate by input fleetgroups (a character vector, where two fleets with the same value are aggregated)
if (!is.null(fleetgroups)) {
if (length(fleetgroups) != nfleets) stop("fleetgroups, if specified, must have length equal to the number of declared fleets")
FleetNames <- unique(fleetgroups)
prof.table_new <- data.frame(matrix(
nrow = nrow(prof.table),
ncol = 3 + length(unique(fleetgroups)),
dimnames = list(
rownames(prof.table),
c(
colnames(prof.table)[1:3],
unique(fleetgroups)
)
)
))
prof.table_new[, 1:3] <- prof.table[, 1:3]
for (rowI in 1:nrow(prof.table)) {
prof.table_new[rowI, -c(1:3)] <- tapply(as.numeric(prof.table[rowI, -c(1:3)]),
FUN = sum,
INDEX = as.numeric(factor(fleetgroups,
levels = unique(fleetgroups)
))
)
}
prof.table <- prof.table_new
nfleets <- ncol(prof.table) - 3
}
# subtract minimum value from each likelihood component (over requested parameter range)
subset <- parvec >= xlim[1] & parvec <= xlim[2]
for (icol in 3:ncol(prof.table)) {
prof.table[, icol] <- prof.table[, icol] -
min(prof.table[subset, icol], na.rm = TRUE)
}
# remove columns that have change less than minfraction change relative to total
column.max <- apply(data.frame(prof.table[, -c(1:3)]), 2, max, na.rm = TRUE)
change.fraction <- column.max / max(prof.table[, 3], na.rm = TRUE)
include <- change.fraction >= minfraction
message(
"Fleet-specific likelihoods showing max change as fraction of total change.\n",
"To change which components are included, change input 'minfraction'.\n",
paste0(utils::capture.output(print(data.frame(
frac_change = round(change.fraction, 4),
include = include
))), collapse = "\n")
)
# subset values and reorder values
# Note: first 3 columns are "model", "Label", and "ALL", and
# are excluded from subsetting process
# a future option to exclude the "ALL" column is possible if requested
prof.table <- prof.table[order(parvec), ]
prof.table <- prof.table[, c(1:3, 3 + intersect(
(1:nfleets)[fleets],
(1:nfleets)[include]
))]
nfleets <- ncol(prof.table) - 3
# replace column names with fleetnames unless "fleetgroup" is used
if (is.null(fleetgroups)) {
for (icol in 4:ncol(prof.table)) {
if (names(prof.table)[icol] %in% FleetNames) {
names(prof.table)[icol] <- fleetnames[which(FleetNames == names(prof.table)[icol])]
}
if (names(prof.table)[icol] %in% paste("X", FleetNames, sep = "")) {
names(prof.table)[icol] <- fleetnames[which(paste("X", FleetNames, sep = "") == names(prof.table)[icol])]
}
}
}
# set default y-limits
if (ymax == "default") ymax <- 1.1 * max(prof.table[subset, -(1:2)], na.rm = TRUE)
ylim <- c(0, ymax)
parvec <- parvec[order(parvec)]
# default colors and plot characters
nlines <- ncol(prof.table) - 2
if (col[1] == "default") col <- rich.colors.short(nlines)
if (pch[1] == "default") pch <- 1:nlines
lwd <- c(lwd.total, rep(lwd, nlines - 1))
cex <- c(cex.total, rep(cex, nlines - 1))
lty <- c(lty.total, rep(lty, nlines - 1))
# return(prof.table)
# make plot
plotprofile <- function() {
plot(0,
type = "n", xlim = xlim, ylim = ylim, xlab = profile.label, ylab = ylab,
yaxs = yaxs, xaxs = xaxs, main = main
)
abline(h = 0, col = "grey")
# optionally add horizontal line at ~1.92 (or other value depending
# on chosen probability)
if (add_cutoff) {
abline(h = 0.5 * qchisq(p = cutoff_prob, df = 1), lty = 2)
}
matplot(parvec, prof.table[, -(1:2)],
type = type,
pch = pch, col = col,
cex = cex, lty = lty, lwd = lwd, add = TRUE
)
if (legend) {
legend(legendloc,
bty = "n", legend = names(prof.table)[-(1:2)],
lwd = lwd, pt.cex = cex, lty = lty, pch = pch, col = col
)
}
box()
}
if (plot) {
plotprofile()
}
if (print) {
save_png(
plotinfo = NULL,
file = "profile_plot_likelihood.png",
plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize
)
plotprofile()
dev.off()
}
out <- data.frame(parvec = parvec, prof.table)
names(out)[1] <- parlabel
return(invisible(out))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.