[feat]: Adds index intervals, follows PR #465

R/plot_index.R

@@ -1,21 +1,28 @@
 #' 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
-#'
-#' Options: 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, but 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
+#' @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
 #' @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
+#' (i.e. select names of fleets to zoom into the plot). Default = NULL
 #' @returns A plot showing the expected and predicted index.
+#' @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.
 #'
-#' @details The input is from an assessment model output file
+#' 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
 #' containing associated caption and alternative text for the figure.
@@ -31,6 +38,14 @@
 #'   interactive = FALSE
 #' )
 #'
+#' # 95% intervals around the observed index, on a log scale
+#' plot_index(
+#'   dat = stockplotr::example_data,
+#'   interval = 0.95,
+#'   log = TRUE,
+#'   interactive = FALSE
+#' )
+#'
 plot_index <- function(
   dat,
   unit_label = NULL,
@@ -40,6 +55,8 @@ plot_index <- function(
   interactive = TRUE,
   module = NULL,
   focus = NULL,
+  interval = 0.95,
+  log = FALSE,
   make_rda = FALSE,
   figures_dir = getwd(),
   ...
@@ -65,18 +82,35 @@ plot_index <- function(
     interactive = interactive,
     module = module
   )
-  if (!is.null(unit_label) && !is.na(unique(prepared_data$unit))) unit_label <- unique(prepared_data$unit)
-  if (is.na(unit_label) || is.null(unit_label)) {
-    cli::cli_alert_warning("No unit label found in the data. Please provide a unit label using the 'unit_label' argument.")
-    unit_label <- ""
+
+  # There's 3 cases handled here
+  # plot_index(dat) (unit_label = NULL): Uses prepared_data$unit if present, otherwise warns / falls back to ""
+  # plot_index(dat, unit_label = "fish/hr"): Skips the is.null() block and keeps "fish/hr"
+  # plot_index(dat, unit_label = ""): Skips the is.null() block and intentionally leaves the y-axis without units ("Estimated Index")
+  if (is.null(unit_label)) {
+    data_units <- unique(stats::na.omit(prepared_data$unit))
+    if (length(data_units) == 1 && nzchar(data_units)) {
+      unit_label <- data_units
+    } else {
+      cli::cli_alert_warning(
+        "No unit label found in the data. Please provide a unit label using the 'unit_label' argument."
+      )
+      unit_label <- ""
+    }
   }
+
   # Set cpue unit label for plot
   u_units <- ifelse(
     is.null(unit_label) | unit_label == "",
     glue::glue("Estimated Index"),
     glue::glue("Estimated Index ({unit_label})")
   )
 
+  # Observation-error intervals for the observed index
+  if (!is.null(interval)) {
+    prepared_data <- index_observed_interval(prepared_data, interval)
+  }
+
   # Subset data if focus
   if (!is.null(focus)) {
     prepared_data <- prepared_data |>
@@ -126,6 +160,31 @@ plot_index <- function(
         minor.ticks = TRUE
       )
     )
+  if (!is.null(interval)) {
+    obs <- prepared_data[
+      grepl("index_observed", prepared_data$label) & !is.na(prepared_data$estimate_lower),
+    ]
+    # Draw the intervals beneath the observed and predicted index
+    intervals <- ggplot2::geom_linerange(
+      data = obs,
+      ggplot2::aes(
+        x = .data$year,
+        ymin = .data$estimate_lower,
+        ymax = .data$estimate_upper,
+        color = .data$group_var
+      ),
+      alpha = 0.6,
+      show.legend = FALSE
+    )
+    plt$layers <- c(list(intervals), plt$layers)
+  }
+  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)")
+  }
   # Overwrite facets from base plot_timeseries bc scales need to be free
   if (!is.null(facet) && "fleet" %in% colnames(prepared_data)) {
     facet <- paste("~", paste(facet, collapse = " + "))
@@ -157,7 +216,7 @@ plot_index <- function(
       cpue.min,
       cpue.max,
       cpue.units,
-      dir=figures_dir
+      dir = figures_dir
     )
 
     # Add key quantities to captions/alt text
@@ -167,7 +226,7 @@ plot_index <- function(
       cpue.min,
       cpue.max,
       cpue.units,
-      dir=figures_dir
+      dir = figures_dir
     )
 
     create_rda(
@@ -179,10 +238,30 @@ plot_index <- function(
       dir = figures_dir
     )
   }
-  
+
   if (exists("selected_module")) {
     rm(selected_module, envir = .GlobalEnv)
   }
-  
+
   plt
 }
+
+# Lognormal intervals around the observed index from its log-scale standard
+# error (uncertainty labeled "se", "sd", or "log_sd") or its coefficient of
+# variation ("cv"); other rows get no interval
+index_observed_interval <- function(dat, interval) {
+  if (!is.numeric(interval) || length(interval) != 1 || interval <= 0 || interval >= 1) {
+    cli::cli_abort("{.arg interval} must be a number between 0 and 1 (e.g., 0.95).")
+  }
+  z <- stats::qnorm(1 - (1 - interval) / 2)
+  type <- tolower(as.character(dat$uncertainty_label))
+  sigma <- dplyr::case_when(
+    type %in% c("se", "sd", "log_sd", "logsd", "stddev") ~ dat$uncertainty,
+    type == "cv" ~ sqrt(log(1 + dat$uncertainty^2)),
+    TRUE ~ NA_real_
+  )
+  observed <- grepl("index_observed", dat$label) & !is.na(sigma)
+  dat$estimate_lower <- ifelse(observed, dat$estimate * exp(-z * sigma), NA_real_)
+  dat$estimate_upper <- ifelse(observed, dat$estimate * exp(z * sigma), NA_real_)
+  dat
+}

tests/testthat/test-plot_index.R

@@ -30,3 +30,64 @@ test_that("plot_index generates plots without errors", {
 #   # expect that both figures dir and the index_figure.rda file exist
 #   check_and_remove_tmp_files("index_figure.rda")
 # })
+
+test_that("plot_index draws lognormal intervals around the observed index", {
+  p <- plot_index(stockplotr::example_data, interval = 0.95, interactive = FALSE)
+  # The intervals are the first layer, beneath the points and lines
+  expect_s3_class(p$layers[[1]]$geom, "GeomLinerange")
+  ints <- p$layers[[1]]$data
+  obs <- stockplotr::example_data[
+    stockplotr::example_data$label %in% "index_observed",
+  ]
+  obs <- obs[order(obs$fleet, obs$year), ]
+  ints <- ints[order(ints$fleet, ints$year), ]
+  expect_equal(nrow(ints), nrow(obs))
+  expect_equal(ints$estimate_lower, obs$estimate * exp(-stats::qnorm(0.975) * obs$uncertainty))
+  expect_equal(ints$estimate_upper, obs$estimate * exp(stats::qnorm(0.975) * obs$uncertainty))
+  expect_no_error(ggplot2::ggplot_build(p))
+
+  # 95% intervals by default, and none with interval = NULL
+  p_default <- plot_index(stockplotr::example_data, interactive = FALSE)
+  expect_equal(p_default$layers[[1]]$data, p$layers[[1]]$data)
+  p <- plot_index(stockplotr::example_data, interval = NULL, interactive = FALSE)
+  expect_false(any(sapply(p$layers, function(l) inherits(l$geom, "GeomLinerange"))))
+
+  expect_error(plot_index(stockplotr::example_data, interval = 95, interactive = FALSE), "between 0 and 1")
+})
+
+test_that("plot_index intervals use log-scale SDs and CVs", {
+  dat <- data.frame(
+    label = "index_observed", year = 1:3, fleet = "a",
+    estimate = c(10, 20, 30), uncertainty = c(0.2, 0.2, NA),
+    uncertainty_label = c("cv", "log_sd", "cv")
+  )
+  out <- index_observed_interval(dat, interval = 0.9)
+  z <- stats::qnorm(0.95)
+  expect_equal(out$estimate_lower[1], 10 * exp(-z * sqrt(log(1 + 0.2^2))))
+  expect_equal(out$estimate_upper[2], 20 * exp(z * 0.2))
+  # No uncertainty, no interval
+  expect_true(is.na(out$estimate_lower[3]))
+
+  # Rceattle (log_sd) and BAM (cv) output
+  for (dat in list(stockplotr::pollock, stockplotr::vermilion_snapper)) {
+    p <- suppressMessages(plot_index(dat, interval = 0.95, interactive = FALSE))
+    expect_gt(nrow(p$layers[[1]]$data), 0)
+    expect_no_error(ggplot2::ggplot_build(p))
+  }
+})
+
+test_that("plot_index can use a log scale", {
+  p <- plot_index(stockplotr::example_data, interval = 0.95, log = TRUE, interactive = FALSE)
+  expect_equal(p$scales$get_scales("y")$trans$name, "log-10")
+  expect_equal(p$labels$y, "Estimated Index (log scale)")
+  expect_no_error(ggplot2::ggplot_build(p))
+})
+
+test_that("plot_index intervals work for FIMS output", {
+  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, interval = 0.95, interactive = FALSE))
+  ints <- p$layers[[1]]$data
+  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)))
+})