Nothing
#' Plot tagging data and fits
#'
#' Plot observed and expected tag recaptures in aggregate and by tag group.
#'
#'
#' @template replist
#' @param subplots vector controlling which subplots to create
#' @param latency period of tag mixing to exclude from plots (in future could
#' be included in SS output)
#' @param taggroups which tag groups to include in the plots. Default=NULL
#' causes all groups to be included.
#' @param rows number or rows of panels for regular plots
#' @param cols number or columns of panels for regular plots
#' @param tagrows number or rows of panels for multi-panel plots
#' @param tagcols number or columns of panels for multi-panel plots
#' @param plot plot to active plot device?
#' @param print print to PNG files?
#' @param pntscalar maximum bubble size for balloon plots; each plot scaled
#' independently based on this maximum size and the values plotted. Often some
#' plots look better with one value and others with a larger or smaller value.
#' Default=2.6
#' @param minnbubble minimum number of years below which blank years will be
#' added to bubble plots to avoid cropping
#' @param pwidth default width of plots printed to files in units of
#' `punits`. Default=7.
#' @param pheight default height width of plots printed to files in units of
#' `punits`. Default=7.
#' @param punits units for `pwidth` and `pheight`. Can be "px"
#' (pixels), "in" (inches), "cm" or "mm". Default="in".
#' @param ptsize point size for plotted text in plots printed to files (see
#' help("png") in R for details). Default=12.
#' @template res
#' @param cex.main character expansion parameter for plot titles
#' @param col1 color for bubbles
#' @param col2 color for lines with expected values
#' @param col3 shading color for observations within latency period
#' @param col4 shading color for observations after latency period
#' @param labels vector of labels for plots (titles and axis labels)
#' @param plotdir directory where PNG files will be written. by default it will
#' be the directory where the model was run.
#' @param verbose return updates of function progress to the R GUI?
#' @author Andre E. Punt, Ian G. Taylor, Ashleigh J. Novak
#' @import ggplot2
#' @export
#' @seealso [SS_plots()], [SS_output()]
SSplotTags <-
function(replist = replist, subplots = 1:10, latency = NULL, taggroups = NULL,
rows = 1, cols = 1,
tagrows = 3, tagcols = 3,
plot = TRUE, print = FALSE,
pntscalar = 2.6, minnbubble = 8,
pwidth = 6.5, pheight = 5.0, punits = "in", ptsize = 10, res = 300, cex.main = 1,
col1 = rgb(0, 0, 1, .7), col2 = "red", col3 = "grey95", col4 = "grey70",
labels = c(
"Year", # 1
"Frequency", # 2
"Tag Group", # 3
"Fit to tag recaptures by tag group", # 4
"Post-latency tag recaptures aggregated across tag groups", # 5
"Observed tag recaptures by year and tag group", # 6
"Residuals for post-latency tag recaptures: (obs-exp)/sqrt(exp)", # 7
"Observed and expected post-latency tag recaptures by year and tag group", # 8
"Summarized observed and expected numbers of recaptures by fleet", # 9
"Pearson residuals by tag group"
), # 10
plotdir = "default",
verbose = TRUE) {
# table to store information on each plot
plotinfo <- NULL
if (plotdir == "default") plotdir <- replist[["inputs"]][["dir"]]
tagdbase2 <- replist[["tagdbase2"]]
if (is.null(tagdbase2) || nrow(tagdbase2) == 0) {
if (verbose) {
message("skipping tag plots because there's no tagging data")
}
} else {
# filter tag groups if requested
if (!is.null(taggroups)) {
tagdbase2 <- tagdbase2[tagdbase2[["Repl."]] %in% taggroups, ]
message(
"Filtered tag groups for plotting based on input vector taggroups\n",
"Plots will show", length(unique(tagdbase2[["Repl."]])),
"out of", length(unique(replist[["tagdbase2"]][["Repl."]])),
"total included in the model."
)
}
# calculations needed for printing to multiple PNG files
grouprange <- unique(tagdbase2[["Repl."]])
ngroups <- length(unique(tagdbase2[["Repl."]]))
npages <- ceiling(ngroups / (tagrows * tagcols))
nseasons <- replist[["nseasons"]]
width <- 0.5 / nseasons
tagreportrates <- replist[["tagreportrates"]]
tagrelease <- replist[["tagrelease"]]
tagsalive <- replist[["tagsalive"]]
tagtotrecap <- replist[["tagtotrecap"]]
if (is.null(latency)) {
latency <- replist[["tagfirstperiod"]]
}
tagfun1 <- function() {
# obs & exp recaps by tag group
par(
mfcol = c(tagrows, tagcols), mar = c(2.5, 2.5, 2, 1),
cex.main = cex.main, oma = c(2, 2, 2, 0)
)
mu <- tagdbase2[["Exp"]] # mean expected number of recaptures
# Get overdispersion parameter from model output
parameters <- replist[["parameters"]]
overdispersion <- parameters %>%
dplyr::filter(stringr::str_detect(.data[["Label"]], "TG_overdispersion_")) %>%
dplyr::select(.data[["Value"]]) # grabs the overdispersion parms for each Tag Group
tau <- overdispersion[["Value"]][1]
k <- c(1:length(tagdbase2[["Exp"]]))
for (i in 1:length(tagdbase2[["Exp"]])) {
k[i] <- mu[i] / (tau - 1) # variance
}
i <- c(1:length(tagdbase2[["Exp"]]))
CI_down <- stats::qnbinom(c(0.975), size = k[i], mu = mu[i])
CI_up <- stats::qnbinom(c(0.025), size = k[i], mu = mu[i])
new_tagdbase2 <- cbind(tagdbase2, CI_up, CI_down)
new_tagdbase2 <- new_tagdbase2 %>%
dplyr::mutate(
CI_down = ifelse(is.nan(CI_down), NA, CI_down),
CI_up = ifelse(is.nan(CI_up), NA, CI_up)
)
new_tagdbase2[["title"]] <- paste("TG_", as.character(new_tagdbase2[["Repl."]]), sep = "")
newfigure1 <- ggplot(new_tagdbase2, aes(x = .data[["Yr"]], y = .data[["Obs"]])) +
geom_bar(aes(fill = as.factor(.data[["Seas"]])), position = "dodge", stat = "identity", alpha = 0.6) +
geom_point(aes(x = .data[["Yr"]], y = .data[["Exp"]], fill = as.factor(.data[["Seas"]])), position = position_dodge(0.9)) +
facet_wrap(forcats::fct_inorder(as.factor(new_tagdbase2[["title"]])), scales = "free") +
geom_errorbar(aes(ymin = CI_down, ymax = CI_up, color = as.factor(.data[["Seas"]])), position = position_dodge(0.9), width = 0.25) +
scale_color_viridis_d() +
scale_fill_viridis_d() +
theme(strip.background = element_blank()) +
labs(x = "Year", y = "Frequency", fill = "Season", color = "Season")
print(newfigure1)
# restore default single panel settings
par(mfcol = c(rows, cols), mar = c(5, 5, 4, 2) + .1, oma = rep(0, 4))
}
# new system which takes latency value as input
tgroups <- sort(unique(tagdbase2[["Repl."]]))
x <- NULL
for (igroup in tgroups) {
# subset results for only 1 tag group
temp <- tagdbase2[tagdbase2[["Repl."]] == igroup, ]
# remove the first rows corresponding to the latency period
temp <- temp[-(1:latency), ]
x <- rbind(x, temp)
}
# obs vs exp tag recaptures by year aggregated across group
tagobs <- aggregate(x[["Obs"]], by = list(x[["Yr.S"]], x[["Repl."]]), FUN = sum, na.rm = TRUE)
tagexp <- aggregate(x[["Exp"]], by = list(x[["Yr.S"]], x[["Repl."]]), FUN = sum, na.rm = TRUE)
Recaps <- data.frame(
Yr.S = tagobs[, 1], Group = tagobs[, 2],
Obs = tagobs[, 3], Exp = tagexp[, 3]
)
xlim <- range(Recaps[["Yr.S"]])
xx2 <- aggregate(Recaps[["Obs"]], by = list(Recaps[["Yr.S"]]), FUN = sum, na.rm = TRUE)
xx3 <- aggregate(Recaps[["Exp"]], by = list(Recaps[["Yr.S"]]), FUN = sum, na.rm = TRUE)
RecAg <- data.frame(Yr.S = xx2[, 1], Obs = xx2[, 2], Exp = xx3[, 2])
tagfun2 <- function() {
# obs vs exp tag recaptures by year aggregated across group
plot(0,
xlim = xlim + c(-0.5, 0.5),
ylim = c(0, max(RecAg[["Obs"]], RecAg[["Exp"]]) * 1.05),
type = "n", xaxs = "i", yaxs = "i",
xlab = labels[1], ylab = labels[2], main = labels[5],
cex.main = cex.main
)
for (iy in 1:nrow(RecAg)) {
xx <- c(
RecAg[["Yr.S"]][iy] - width, RecAg[["Yr.S"]][iy] - width,
RecAg[["Yr.S"]][iy] + width, RecAg[["Yr.S"]][iy] + width
)
yy <- c(0, RecAg[["Obs"]][iy], RecAg[["Obs"]][iy], 0)
polygon(xx, yy, col = col4)
}
lines(RecAg[["Yr.S"]], RecAg[["Exp"]], type = "o", pch = 16, lty = 1, lwd = 2)
}
Recaps[["Pearson"]] <- (Recaps[["Obs"]] - Recaps[["Exp"]]) / sqrt(Recaps[["Exp"]])
Recaps[["Pearson"]][Recaps[["Exp"]] == 0] <- NA
tagfun3 <- function() {
# bubble plot of observed recapture data
plottitle <- labels[6]
bubble3(
x = Recaps[["Yr.S"]], y = Recaps[["Group"]], z = Recaps[["Obs"]],
xlab = labels[1], ylab = labels[3], col = col1,
las = 1, main = plottitle, cex.main = cex.main, maxsize = pntscalar,
allopen = FALSE, minnbubble = minnbubble
)
}
tagfun4 <- function() {
# bubble plot of residuals
plottitle <- labels[7]
bubble3(
x = Recaps[["Yr.S"]], y = Recaps[["Group"]], z = Recaps[["Pearson"]],
xlab = labels[1], ylab = labels[3], col = col1,
las = 1, main = plottitle, cex.main = cex.main, maxsize = pntscalar,
allopen = FALSE, minnbubble = minnbubble
)
}
tagfun5 <- function() {
# line plot by year and group
plottitle <- labels[8]
plot(0,
type = "n", xlim = range(Recaps[["Yr.S"]]),
ylim = range(Recaps[["Group"]]) + c(0, 1), xlab = labels[1], ylab = labels[3],
main = plottitle, cex.main = cex.main
)
rescale <- .9 * min(ngroups - 1, 5) / max(Recaps[["Obs"]], Recaps[["Exp"]])
for (igroup in sort(unique(Recaps[["Group"]]))) {
lines(Recaps[["Yr.S"]][Recaps[["Group"]] == igroup],
igroup + 0 * Recaps[["Obs"]][Recaps[["Group"]] == igroup],
col = "grey", lty = 3
)
points(Recaps[["Yr.S"]][Recaps[["Group"]] == igroup],
igroup + rescale * Recaps[["Obs"]][Recaps[["Group"]] == igroup],
type = "o", pch = 16, cex = .5
)
lines(Recaps[["Yr.S"]][Recaps[["Group"]] == igroup],
igroup + rescale * Recaps[["Exp"]][Recaps[["Group"]] == igroup],
col = col2, lty = "42", lwd = 2
)
}
legend("topleft",
bty = "n", lty = c("91", "42"),
pch = c(16, NA), pt.cex = c(.5, NA),
col = c(1, 2), lwd = c(1, 2), legend = c("Observed", "Expected")
)
}
tagfun6 <- function() {
# a function to plot tag parameters after transformation
# into reporting rate and tag loss quantities
par(mfrow = c(2, 2))
# first plot is reporting rate parameters
barplot(
height = tagreportrates[["Init_Reporting"]],
names.arg = tagreportrates[["Fleet"]], ylim = c(0, 1), yaxs = "i",
ylab = "Reporting rate", xlab = "Fleet number",
main = "Initial reporting rate"
)
box()
# second plot shows any decay in reporting rate over time
matplot(0:5, exp((0:5) %*% t(tagreportrates[["Report_Decay"]])),
type = "l", lwd = 3, lty = 1, col = rich.colors.short(nrow(tagreportrates)),
ylim = c(0, 1.05), yaxs = "i",
ylab = "Reporting rate", xlab = "Time at liberty (years)",
main = "Reporting rate decay"
)
# third plot shows initial tag loss
barplot(
height = tagrelease[["Init_Loss"]],
names.arg = tagrelease[["Fleet"]], ylim = c(0, 1), yaxs = "i",
ylab = "Initial tag loss", xlab = "Tag group",
main = "Initial tag loss\n(fraction of tags lost at time of tagging)"
)
box()
# fourth plot shows chronic tag loss
barplot(
height = tagrelease[["Chron_Loss"]],
names.arg = tagrelease[["Fleet"]], ylim = c(0, 1), yaxs = "i",
ylab = "Chronic tag loss", xlab = "Tag group",
main = "Chronic tag loss\n(fraction of tags lost per year)"
)
box()
# restore default single panel settings
par(mfcol = c(rows, cols), mar = c(5, 5, 4, 2) + .1, oma = rep(0, 4))
}
tagfun7 <- function() {
# a function to plot the "tags alive" matrix
xvals <- as.numeric(substring(names(tagsalive)[-1], 7))
matplot(xvals, t(tagsalive[, -1]),
type = "l", lwd = 3,
col = rich.colors.short(nrow(tagsalive)),
xlab = "Period at liberty",
ylab = "Estimated number of alive tagged fish",
main = "'Tags alive' by tag group"
)
abline(h = 0, col = "grey")
}
tagfun8 <- function() {
# a function to plot the "total recaptures" matrix
xvals <- as.numeric(substring(names(tagtotrecap)[-1], 7))
matplot(xvals, t(tagtotrecap[, -1]),
type = "l", lwd = 3,
col = rich.colors.short(nrow(tagtotrecap)),
xlab = "Period at liberty",
ylab = "Estimated number of recaptures",
main = "'Total recaptures' by tag group"
)
abline(h = 0, col = "grey")
}
tagdata <- replist[["tagdbase1"]]
# need to get the max number of fleets you have so you can rep over that.
max_num_fleets <- max(tagdata[["Fleet"]])
# make a new column to bind to your other frame that contains
# the breakdown of obs + exp recaps by fleet
expected_by_fleets <- as.data.frame(rep(tagdbase2[["Exp"]], each = max_num_fleets))
names(expected_by_fleets)[1] <- "Expected" # rename column
new_tagdata <- cbind(tagdata, expected_by_fleets) # bind new column to tagdata dataframe
fleet_numbers <- new_tagdata %>%
dplyr::mutate(
Numbers_Obs = round(.data[["Obs"]] * .data[["Nsamp_adj"]]),
Numbers_Exp = round(.data[["Exp"]] * .data[["Expected"]])
) # generate exp. and obs. recaptures by fleet
fleet_numbers2 <- fleet_numbers %>%
dplyr::group_by(.data[["Fleet"]], .data[["Yr"]]) %>%
dplyr::summarize(sum_exp = sum(.data[["Numbers_Exp"]]), sum_obs = sum(.data[["Numbers_Obs"]]))
fleet_numbers2[["fleet_title"]] <- paste("Fleet_", as.character(fleet_numbers2[["Fleet"]]), sep = "")
mycols <- rep("black", max_num_fleets) # set colors for plotting expected values
myfill <- rep("gray35", max_num_fleets)
tagfun9 <- function() {
# summarized observed and expected numbers of recaptures by fleet
fleet_plot2 <- ggplot(fleet_numbers2, aes(x = .data[["Yr"]], y = .data[["sum_obs"]])) +
geom_bar(aes(fill = as.factor(.data[["Fleet"]])), stat = "identity", alpha = 0.5) +
geom_line(aes(y = .data[["sum_exp"]], color = as.factor(.data[["Fleet"]])), linetype = 1, size = 1) +
facet_wrap(forcats::fct_inorder(as.factor(fleet_numbers2[["fleet_title"]])), scales = "free") +
scale_fill_manual(values = myfill) +
scale_color_manual(values = mycols) +
labs(x = "Year", y = "Number of recaptures", subtitle = "Observed (bar) and expected (line)") +
theme(strip.background = element_blank(), legend.position = "none")
print(fleet_plot2)
}
# Calculate Pearson Residuals by fleet
fleetnumbers_PRs <- fleet_numbers %>%
dplyr::mutate(Pearson = (.data[["Numbers_Obs"]] - .data[["Numbers_Exp"]]) / sqrt(.data[["Numbers_Exp"]]))
fleetnumbers_PRs[["Pearson"]][is.nan(fleetnumbers_PRs[["Pearson"]])] <- NA
fleetnumbers_PRs[["Pearson"]][!is.finite(fleetnumbers_PRs[["Pearson"]])] <- NA
# Adjust names
fleetnumbers_PRs[["TGTitle"]] <- paste("TG_", as.character(fleetnumbers_PRs[["Repl."]]), sep = "")
# Set limits and breaks of residuals for plotting
limits <- as.numeric(round(range(fleetnumbers_PRs[["Pearson"]], na.rm = TRUE)))
breaks <- seq(limits[1], limits[2], by = 3)
legend_size <- c(abs(breaks))
tagfun10 <- function() {
# pearson residuals by tag group and fleets
pearson_plot1 <- ggplot(fleetnumbers_PRs, aes(x = .data[["Yr"]], y = forcats::fct_inorder(as.factor(.data[["Fleet"]])))) +
geom_point(aes(size = .data[["Pearson"]], color = .data[["Pearson"]]), alpha = 0.75, na.rm = TRUE) +
scale_color_continuous(name = "Pearson\nResiduals", limits = limits, breaks = breaks, type = "viridis") +
scale_size_continuous(name = "Pearson\nResiduals", limits = limits, breaks = breaks) +
facet_wrap(forcats::fct_inorder(as.factor(fleetnumbers_PRs[["TGTitle"]])), scales = "free_x") +
xlab("Year") +
ylab("Fleets") +
guides(color = guide_legend(), size = guide_legend(override.aes = list(size = legend_size)))
print(pearson_plot1)
}
# make plots
if (plot) {
if (1 %in% subplots) tagfun1()
if (2 %in% subplots) tagfun2()
if (3 %in% subplots) tagfun3()
if (4 %in% subplots) tagfun4()
if (5 %in% subplots) tagfun5()
if (6 %in% subplots) tagfun6()
if (7 %in% subplots) tagfun7()
if (8 %in% subplots) tagfun8()
if (9 %in% subplots) tagfun9()
if (10 %in% subplots) tagfun10()
}
# send to files if requested
if (print) {
if (1 %in% subplots) {
file <- "tags_by_group.png"
caption <- "Expected and observed recaptures by tag group"
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
tagfun1()
dev.off()
}
if (2 %in% subplots) {
file <- "tags_aggregated.png"
caption <- labels[5]
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
tagfun2()
dev.off()
}
if (3 %in% subplots) {
file <- "tags_data_bubbleplot.png"
caption <- labels[6]
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
tagfun3()
dev.off()
}
if (4 %in% subplots) {
file <- "tags_residuals.png"
caption <- labels[7]
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
tagfun4()
dev.off()
}
if (5 %in% subplots) {
file <- "tags_lines.png"
caption <- labels[8]
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
tagfun5()
dev.off()
}
if (6 %in% subplots) {
file <- "tags_parameters.png"
caption <- "Tag-related parameters"
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
tagfun6()
dev.off()
}
if (7 %in% subplots) {
file <- "tags_alive.png"
caption <- "'Tags alive' by tag group"
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
tagfun7()
dev.off()
}
if (8 %in% subplots) {
file <- "tags_total_recaptures.png"
caption <- "Total tag recaptures"
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
tagfun8()
dev.off()
}
if (9 %in% subplots) {
file <- "summarized_recaptures_fleet.png"
caption <- "summarized observed and expected numbers of recaptures by fleet"
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
tagfun9()
dev.off()
}
if (10 %in% subplots) {
file <- "pearson_residuals_taggroup.png"
caption <- "Pearson residuals by tag group"
plotinfo <- save_png(
plotinfo = plotinfo, file = file, plotdir = plotdir, pwidth = pwidth,
pheight = pheight, punits = punits, res = res, ptsize = ptsize,
caption = caption
)
tagfun10()
dev.off()
}
}
flush.console()
} # end if data
if (!is.null(plotinfo)) plotinfo[["category"]] <- "Tag"
return(invisible(plotinfo))
} # end SSplotTags
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.