[feat]: Add relative option to plot_index, follows #468

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")
+}

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)
+})