diff --git a/R/plot_index.R b/R/plot_index.R index 3111531a..e788f93b 100644 --- a/R/plot_index.R +++ b/R/plot_index.R @@ -1,16 +1,32 @@ #' Plot Index of Abundance #' #' @inheritParams plot_spawning_biomass -#' @param facet Character vector. Column name or names used for faceting. "fleet" -#' is always added on to any faceting selections. Default: NULL, but may include, -#' but is not limited to: "fleet", "year", "area", NULL -#' @param unit_label String. Index of abundance/CPUE units. Default: NULL +#' @param facet Character vector. Column name or names used for faceting. "fleet" is always added on to any faceting selections +#' +#' Default: NULL +#' +#' Options: May include, but is not limited to: "fleet", "year", "area", NULL +#' @param unit_label String. Index of abundance/CPUE units +#' +#' Default: NULL #' @param interval Number. Coverage of intervals around the observed index #' based on its observation error (e.g., 0.95), drawn as vertical lines. -#' NULL draws no intervals, default = 0.95 +#' NULL draws no intervals. +#' +#' Default: 0.95 #' @param log Logical. TRUE/FALSE; plot the index on a log scale. Indices #' with lognormal errors are often easier to judge on a log scale, where the -#' intervals are symmetric. Default = FALSE +#' intervals are symmetric. +#' +#' Default: FALSE +#' @param relative Logical. TRUE/FALSE; divide each fleet's observed and +#' predicted index (and intervals) by the mean of its observed index, and plot +#' all fleets on one panel, colored by fleet. Indices with conflicting trends +#' (data conflicts) are a common cause of poor fits and are easiest to see +#' this way (`r4ss::SSplotIndices()` subplot 9). Any `group` other than +#' "fleet" is added to `facet`. +#' +#' Default: FALSE #' @param focus Character vector. Column name or names used to filter data. #' (i.e. select names of fleets to zoom into the plot). Default = NULL #' @returns A plot showing the expected and predicted index. @@ -22,6 +38,14 @@ #' \eqn{\sigma = \sqrt{\log(1 + CV^2)}}. Observations without an #' uncertainty have no interval. #' +#' @details Intervals assume lognormal observation errors, as in Stock +#' Synthesis, FIMS, BAM, and Rceattle: \eqn{o \exp(\pm z \cdot \sigma)}, +#' where \eqn{o} is the observed index and \eqn{\sigma} is its log-scale +#' standard error. Uncertainty labeled "se", "sd", or "log_sd" is used as +#' \eqn{\sigma}, and a coefficient of variation ("cv") is converted to +#' \eqn{\sigma = \sqrt{\log(1 + CV^2)}}. Observations without an +#' uncertainty have no interval. +#' #' The input is from an assessment model output file #' translated to a standardized output (\link[stockplotr]{convert_output}). #' There are options to return a `ggplot2` object or export an .rda object @@ -45,6 +69,12 @@ #' log = TRUE, #' interactive = FALSE #' ) +#' # All indices relative to their means on one panel +#' plot_index( +#' dat = stockplotr::example_data, +#' relative = TRUE, +#' interactive = FALSE +#' ) #' plot_index <- function( dat, @@ -57,13 +87,22 @@ plot_index <- function( focus = NULL, interval = 0.95, log = FALSE, + relative = FALSE, make_rda = FALSE, figures_dir = getwd(), ... ) { check_data(dat) - if (is.null(facet) || "fleet" %notin% facet) { + if (relative) { + # All fleets on one panel, distinguished by color + if (!is.null(group) && group != "fleet") { + facet <- c(facet, group) + } + group <- "fleet" + facet <- setdiff(facet, "fleet") + if (length(facet) == 0) facet <- NULL + } else if (is.null(facet) || "fleet" %notin% facet) { facet <- c(facet, "fleet") } # Filter data @@ -110,6 +149,9 @@ plot_index <- function( if (!is.null(interval)) { prepared_data <- index_observed_interval(prepared_data, interval) } + if (relative) { + prepared_data <- index_relative(prepared_data) + } # Subset data if focus if (!is.null(focus)) { @@ -126,14 +168,24 @@ plot_index <- function( group <- processed_data[[2]] facet <- processed_data[[3]] - # move label to grouping and set grouping into facet - if (!is.null(group)) { + # move label to grouping and set grouping into facet (for relative + # indices, the group is shown by color on one panel) + if (!is.null(group) && !relative) { facet <- c(facet, group) } # process_data() drops columns with a single value (e.g., an index from one # fleet), so facet only by the columns that remain facet <- intersect(facet, colnames(prepared_data)) if (length(facet) == 0) facet <- NULL + # Relative indices are on one panel; columns with one value for each fleet + # (e.g., the month of each survey) are shown by the fleet colors already + if (relative && !is.null(facet) && "fleet" %in% colnames(prepared_data)) { + nested <- vapply(facet, function(v) { + all(tapply(prepared_data[[v]], prepared_data$fleet, function(x) length(unique(x)) == 1)) + }, logical(1)) + facet <- facet[!nested] + if (length(facet) == 0) facet <- NULL + } # transform prep data so group_var = label # prepared_data <- prepared_data |> # dplyr::mutate( @@ -178,12 +230,23 @@ plot_index <- function( ) plt$layers <- c(list(intervals), plt$layers) } + if (relative) { + plt$layers <- c( + list(ggplot2::geom_hline(yintercept = 1, color = "gray70")), + plt$layers + ) + plt <- plt + ggplot2::labs(y = "Index relative to its mean") + # The model legend only has an entry for each model + if (length(unique(prepared_data$model)) == 1) { + plt <- plt + ggplot2::guides(shape = "none", linetype = "none") + } + } if (log) { # plot_obsvpred() expands the y-axis to include zero, which has no log plt$layers <- Filter(function(l) !inherits(l$geom, "GeomBlank"), plt$layers) plt <- suppressMessages( plt + ggplot2::scale_y_log10(labels = scales::label_comma()) - ) + ggplot2::labs(y = "Estimated Index (log scale)") + ) + ggplot2::labs(y = paste(if (relative) "Index relative to its mean" else "Estimated Index", "(log scale)")) } # Overwrite facets from base plot_timeseries bc scales need to be free if (!is.null(facet) && "fleet" %in% colnames(prepared_data)) { @@ -265,3 +328,20 @@ index_observed_interval <- function(dat, interval) { dat$estimate_upper <- ifelse(observed, dat$estimate * exp(z * sigma), NA_real_) dat } + +# Divide the observed and predicted index (and intervals) of each fleet and +# model by the mean of its observed index +index_relative <- function(dat) { + keys <- intersect(c("model", "fleet"), names(dat)) + dat |> + dplyr::group_by(dplyr::across(dplyr::all_of(keys))) |> + dplyr::mutate( + index_mean = mean(.data$estimate[grepl("index_observed", .data$label)], na.rm = TRUE), + dplyr::across( + dplyr::any_of(c("estimate", "estimate_lower", "estimate_upper")), + ~ .x / .data$index_mean + ) + ) |> + dplyr::ungroup() |> + dplyr::select(-"index_mean") +} diff --git a/man/plot_index.Rd b/man/plot_index.Rd index 27152f2d..20003694 100644 --- a/man/plot_index.Rd +++ b/man/plot_index.Rd @@ -14,6 +14,7 @@ plot_index( focus = NULL, interval = 0.95, log = FALSE, + relative = FALSE, make_rda = FALSE, figures_dir = getwd(), ... @@ -26,7 +27,9 @@ returned from \link[stockplotr]{convert_output}. If inputting a list of tibbles, the first tibble's reference point defined in \code{ref_line} is used to plot a reference line or calculate relative spawning biomass.} -\item{unit_label}{String. Index of abundance/CPUE units. Default: NULL} +\item{unit_label}{String. Index of abundance/CPUE units + +Default: NULL} \item{group}{String. Single column that groups the data. @@ -35,9 +38,11 @@ Set group = "none" to summarize data over all indexing values. Default: NULL Options: Including, but not limited to: "year", "area", "fleet", "sex", "none", NULL} -\item{facet}{Character vector. Column name or names used for faceting. "fleet" -is always added on to any faceting selections. Default: NULL, but may include, -but is not limited to: "fleet", "year", "area", NULL} +\item{facet}{Character vector. Column name or names used for faceting. "fleet" is always added on to any faceting selections + +Default: NULL + +Options: May include, but is not limited to: "fleet", "year", "area", NULL} \item{interactive}{Logical. TRUE/FALSE; indicate whether the environment is interactive. @@ -56,11 +61,24 @@ module_name in the console. @seealso \code{\link[=filter_data]{filter_data()}}} \item{interval}{Number. Coverage of intervals around the observed index based on its observation error (e.g., 0.95), drawn as vertical lines. -NULL draws no intervals, default = 0.95} +NULL draws no intervals. + +Default: 0.95} \item{log}{Logical. TRUE/FALSE; plot the index on a log scale. Indices with lognormal errors are often easier to judge on a log scale, where the -intervals are symmetric. Default = FALSE} +intervals are symmetric. + +Default: FALSE} + +\item{relative}{Logical. TRUE/FALSE; divide each fleet's observed and +predicted index (and intervals) by the mean of its observed index, and plot +all fleets on one panel, colored by fleet. Indices with conflicting trends +(data conflicts) are a common cause of poor fits and are easiest to see +this way (\code{r4ss::SSplotIndices()} subplot 9). Any \code{group} other than +"fleet" is added to \code{facet}. + +Default: FALSE} \item{make_rda}{Logical. TRUE/FALSE; indicate whether to save the object and make an automated caption and alternative text in the form of an \code{rda} object. If TRUE, @@ -91,6 +109,14 @@ standard error. Uncertainty labeled "se", "sd", or "log_sd" is used as \eqn{\sigma = \sqrt{\log(1 + CV^2)}}. Observations without an uncertainty have no interval. +Intervals assume lognormal observation errors, as in Stock +Synthesis, FIMS, BAM, and Rceattle: \eqn{o \exp(\pm z \cdot \sigma)}, +where \eqn{o} is the observed index and \eqn{\sigma} is its log-scale +standard error. Uncertainty labeled "se", "sd", or "log_sd" is used as +\eqn{\sigma}, and a coefficient of variation ("cv") is converted to +\eqn{\sigma = \sqrt{\log(1 + CV^2)}}. Observations without an +uncertainty have no interval. + The input is from an assessment model output file translated to a standardized output (\link[stockplotr]{convert_output}). There are options to return a \code{ggplot2} object or export an .rda object @@ -110,6 +136,12 @@ plot_index( log = TRUE, interactive = FALSE ) +# All indices relative to their means on one panel +plot_index( + dat = stockplotr::example_data, + relative = TRUE, + interactive = FALSE +) } \seealso{ diff --git a/tests/testthat/test-plot_index.R b/tests/testthat/test-plot_index.R index 4d7e2299..2ae88b18 100644 --- a/tests/testthat/test-plot_index.R +++ b/tests/testthat/test-plot_index.R @@ -91,3 +91,47 @@ test_that("plot_index intervals work for FIMS output", { rows <- est[est$label == "index_expected" & !is.na(est$observed), ] expect_equal(ints$estimate_lower, rows$observed * exp(-stats::qnorm(0.975) * exp(rows$log_sd))) }) + +test_that("plot_index can plot all indices relative to their means", { + p <- plot_index(stockplotr::example_data, relative = TRUE, interval = 0.95, interactive = FALSE) + b <- ggplot2::ggplot_build(p) + # One panel, with fleets shown by color + expect_equal(length(unique(b$layout$layout$PANEL)), 1) + expect_equal(p$labels$y, "Index relative to its mean") + expect_s3_class(p$layers[[1]]$geom, "GeomHline") + + # Observed values average one for each fleet, and the predicted index and + # intervals are scaled by the same mean + pts <- p$layers[[3]]$data + obs <- pts[pts$label == "index_observed", ] + expect_equal(as.vector(tapply(obs$estimate, obs$fleet, mean)), c(1, 1)) + raw <- stockplotr::example_data[stockplotr::example_data$label %in% c("index_observed", "index_predicted"), ] + tri <- raw[raw$fleet == "Triennial", ] + tri_mean <- mean(tri$estimate[tri$label == "index_observed"]) + pred <- p$layers[[4]]$data + expect_equal( + sort(pred$estimate[pred$fleet == "Triennial"]), + sort(tri$estimate[tri$label == "index_predicted"] / tri_mean) + ) + # Intervals keep their width relative to each observation + ints <- p$layers[[2]]$data + ints0 <- plot_index(stockplotr::example_data, interval = 0.95, interactive = FALSE)$layers[[1]]$data + ints <- ints[order(ints$fleet, ints$year), ] + ints0 <- ints0[order(ints0$fleet, ints0$year), ] + expect_equal(ints$estimate_lower / ints$estimate, ints0$estimate_lower / ints0$estimate) + expect_no_error(b) + + # Survey months that are nested within fleets do not split the panel + p <- suppressMessages(plot_index(stockplotr::pollock, relative = TRUE, interactive = FALSE)) + expect_equal(length(unique(ggplot2::ggplot_build(p)$layout$layout$PANEL)), 1) + + # Works on a log scale, with one fleet, and with FIMS output + p <- plot_index(stockplotr::example_data, relative = TRUE, log = TRUE, interactive = FALSE) + expect_equal(p$labels$y, "Index relative to its mean (log scale)") + expect_no_error(ggplot2::ggplot_build(p)) + expect_no_error(ggplot2::ggplot_build(suppressMessages(plot_index(stockplotr::ss3_simple, relative = TRUE, interactive = FALSE)))) + est <- readRDS(test_path("fixtures", "fims_output", "fims_estimates.rds")) + fims <- suppressMessages(convert_output(est, model = "fims", start_year = 2001)) + p <- suppressMessages(plot_index(fims, relative = TRUE, interactive = FALSE)) + expect_equal(mean(p$layers[[2]]$data$estimate[p$layers[[2]]$data$label == "index_observed"]), 1) +})