diff --git a/NAMESPACE b/NAMESPACE
index f7fd3c70..6f409ca8 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -22,6 +22,7 @@ export(plot_natural_mortality)
 export(plot_obsvpred)
 export(plot_recruitment)
 export(plot_recruitment_deviations)
+export(plot_runs_test)
 export(plot_selectivity)
 export(plot_spawning_biomass)
 export(plot_spr)
diff --git a/R/conout_fims.R b/R/conout_fims.R
new file mode 100644
index 00000000..490000b1
--- /dev/null
+++ b/R/conout_fims.R
@@ -0,0 +1,143 @@
+#' Convert FIMS output to the standard output format
+#'
+#' Reshape the estimates from a FIMS model (`FIMS::get_estimates()`) into the
+#' standardized long format returned by [convert_output()].
+#'
+#' @param file A tibble returned by `FIMS::get_estimates()`, a `FIMSFit`
+#' object, or the path to an .rds file containing either.
+#' @param out_new The empty standard output data frame from
+#' [convert_output()].
+#' @param start_year Number. First year of the model. FIMS reports years as
+#' indices (`year_i`, 1 = first model year), which are converted to years as
+#' `year_i + start_year - 1`. NULL keeps the indices.
+#'
+#' @returns A data frame with the columns of `out_new`.
+#'
+#' @details Each data point (index, landings, and age and length compositions)
+#' is one row in FIMS, with the observed value in `observed` and the expected
+#' value in `estimated`. These are split into two rows, labeled
+#' "*_observed" and "*_predicted" as for other models (e.g., "index_observed"
+#' and "index_predicted"). Observed rows hold the observation error in
+#' `uncertainty` and the negative log-likelihood of the observation in
+#' `likelihood`. FIMS reports the observation error as `log_sd`, the log of
+#' its standard deviation. For lognormal data this standard deviation
+#' (\eqn{\sigma}) is on the log scale, so it is converted to the coefficient
+#' of variation of the observation, \eqn{\sqrt{\exp(\sigma^2) - 1}}, and
+#' labeled "cv"; for other data `exp(log_sd)` is labeled "se". Observations
+#' coded as missing (-999) are not observations and are not included.
+#' Compositions are kept as counts, as reported by FIMS.
+#'
+#' All other estimates keep one row each, with the value in `estimate` and
+#' its standard error in `uncertainty`. FIMS codes estimates that could not
+#' be calculated as -999; these are NA. Rows in years after the last year of
+#' data (e.g., spawning biomass at the start of the year after the model
+#' ends) have `era` "fore"; other rows with a year have `era` "time".
+#'
+#' @noRd
+conout_fims <- function(file, out_new, start_year = NULL) {
+  dat <- read_fims_estimates(file)
+  names(dat) <- tolower(names(dat))
+  if (!all(c("label", "estimated") %in% names(dat))) {
+    cli::cli_abort(c(
+      "{.arg file} is not FIMS output.",
+      "i" = "Supply the output of {.code FIMS::get_estimates()} or a {.cls FIMSFit} object."
+    ))
+  }
+  col <- function(name) {
+    if (name %in% names(dat)) dat[[name]] else rep(NA, nrow(dat))
+  }
+  year_i <- suppressWarnings(as.numeric(col("year_i")))
+  if (is.null(start_year)) {
+    cli::cli_alert_info(
+      "FIMS reports years as indices (1 = first model year); use {.arg start_year} to convert them to years."
+    )
+    year <- year_i
+  } else {
+    year <- year_i + start_year - 1
+  }
+  observed <- suppressWarnings(as.numeric(col("observed")))
+  observed[observed %in% -999] <- NA
+  # The last year of data; later years are after the model ends
+  data_years <- year_i[!is.na(observed)]
+  last_year <- if (length(data_years) > 0) max(data_years) else max(year_i, na.rm = TRUE)
+  era <- dplyr::case_when(
+    is.na(year_i) ~ NA_character_,
+    year_i > last_year ~ "fore",
+    TRUE ~ "time"
+  )
+  uncertainty <- dplyr::na_if(suppressWarnings(as.numeric(col("uncertainty"))), -999)
+  base <- data.frame(
+    module_name = as.character(col("module_name")),
+    label = as.character(dat$label),
+    era = era,
+    year = year,
+    initial = suppressWarnings(as.numeric(col("input"))),
+    estimate = dplyr::na_if(suppressWarnings(as.numeric(dat$estimated)), -999),
+    uncertainty = uncertainty,
+    uncertainty_label = ifelse(is.na(uncertainty), NA_character_, "se"),
+    likelihood = NA_real_,
+    fleet = as.character(col("fleet")),
+    age = suppressWarnings(as.numeric(col("age"))),
+    len_bins = suppressWarnings(as.numeric(col("length"))),
+    stringsAsFactors = FALSE
+  )
+
+  # Data points: split into observed and predicted rows
+  stream <- dplyr::case_when(
+    grepl("^index_?expected$", base$label, ignore.case = TRUE) ~ "index",
+    grepl("^(catch|landings)_?expected$", base$label, ignore.case = TRUE) ~ "landings",
+    grepl("^age_?comp(osition)?_?expected$", base$label, ignore.case = TRUE) ~ "age_composition",
+    grepl("^length_?comp(osition)?_?expected$", base$label, ignore.case = TRUE) ~ "length_composition",
+    TRUE ~ NA_character_
+  )
+  # Fleets with data for a stream; expected values for other fleets (e.g.,
+  # the expected index of a fishery without an index) keep their FIMS label
+  has_data <- stats::ave(!is.na(observed), stream, base$fleet, FUN = any)
+  is_data <- !is.na(stream) & has_data %in% TRUE
+  obs_rows <- base[is_data & !is.na(observed), , drop = FALSE]
+  if (nrow(obs_rows) > 0) {
+    is_obs <- is_data & !is.na(observed)
+    sigma <- exp(suppressWarnings(as.numeric(col("log_sd")))[is_obs])
+    # A lognormal sigma is on the log scale; report the CV of the observation
+    lognormal <- grepl("^(d?lnorm|log_?normal)$", col("distribution")[is_obs], ignore.case = TRUE)
+    obs_rows$label <- paste0(stream[is_obs], "_observed")
+    obs_rows$estimate <- observed[is_obs]
+    obs_rows$uncertainty <- ifelse(lognormal, sqrt(exp(sigma^2) - 1), sigma)
+    obs_rows$uncertainty_label <- dplyr::case_when(
+      is.na(sigma) ~ NA_character_,
+      lognormal ~ "cv",
+      TRUE ~ "se"
+    )
+    obs_rows$initial <- NA_real_
+    # Stock Synthesis reports the negative log-likelihood of each observation
+    obs_rows$likelihood <- -suppressWarnings(as.numeric(col("likelihood")))[is_obs]
+  }
+  pred_rows <- base[is_data, , drop = FALSE]
+  pred_rows$label <- paste0(stream[is_data], "_predicted")
+  out <- rbind(base[!is_data, , drop = FALSE], obs_rows, pred_rows)
+
+  missing <- setdiff(names(out_new), names(out))
+  out[missing] <- NA
+  out[, names(out_new)]
+}
+
+# The estimates of a FIMS model from get_estimates() output, a FIMSFit
+# object, or an .rds file with either
+read_fims_estimates <- function(file) {
+  if (is.character(file) && length(file) == 1) {
+    file <- readRDS(file)
+  }
+  if (inherits(file, "FIMSFit")) {
+    get_estimates <- tryCatch(
+      getExportedValue("FIMS", "get_estimates"),
+      error = function(e) {
+        cli::cli_abort("The {.pkg FIMS} package is needed to read a {.cls FIMSFit} object.")
+      }
+    )
+    file <- get_estimates(file)
+  }
+  if (!is.data.frame(file)) {
+    cli::cli_abort("{.arg file} must be FIMS output from {.code FIMS::get_estimates()} or a {.cls FIMSFit} object.")
+  }
+  as.data.frame(file)
+}
diff --git a/R/convert_output.R b/R/convert_output.R
index 62bda5b8..d3b2ac99 100644
--- a/R/convert_output.R
+++ b/R/convert_output.R
@@ -3,7 +3,8 @@
 #' Format stock assessment output file to a standardized tibble.
 #'
 #' @param file Path. Assessment model output file path. Currently accepts SS3
-#' `r4ss::ss_output()` or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS model object.
+#' `r4ss::ss_output()` or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS output
+#' (the tibble returned by `FIMS::get_estimates()`, a `FIMSFit` object, or an .rds file with either).
 #' @param model String. Assessment model used in evaluation
 #'
 #' Default: NULL
@@ -20,6 +21,11 @@
 #' in the file path provided.
 #'
 #' Default: NULL
+#' @param start_year Number. First year of the model, for FIMS output only.
+#' FIMS reports years as indices (1 = first model year), which are converted
+#' to years when `start_year` is given.
+#'
+#' Default: NULL (years are left as indices)
 #'
 #' @author Samantha Bredeck
 #'
@@ -30,6 +36,14 @@
 #' version of a model output file. Converted data frame is always returned.
 #' It will also be saved if save_dir is not NULL.
 #'
+#' For FIMS, each data point (index, landings, and age and length
+#' compositions) is split into an observed row (e.g., "index_observed") and a
+#' predicted row (e.g., "index_predicted"), as for the other models. Observed
+#' rows hold the observation error (a log-scale standard error for lognormal
+#' data) in `uncertainty` and the negative log-likelihood in `likelihood`.
+#' Other estimates are relabeled with `fims_var_names.csv` (e.g.,
+#' "expected_recruitment" becomes "recruitment").
+#'
 #'
 #' @export
 #'
@@ -46,7 +60,8 @@ convert_output <- function(
   file,
   model = NULL,
   fleet_names = NULL,
-  save_dir = NULL
+  save_dir = NULL,
+  start_year = NULL
 ) {
   
   if (is.character(file) && grepl("Report.sso", file)) {
@@ -163,8 +178,12 @@ convert_output <- function(
         cli::cli_abort("Unknown file type. Please indicate model.")
       )
     }
-  } else {
+  } else if (is.null(model)) {
     model <- switch(class(file)[1],
+      "FIMSFit" = {
+        cli::cli_alert_info("Processing FIMS output file...")
+        "fims"
+      },
       "fims" = {
         cli::cli_alert_info("Processing FIMS output file...")
         "fims"
@@ -210,17 +229,7 @@ convert_output <- function(
     cli::cli_abort("JABBA output not currently compatible.")
     #### FIMS ####
   } else if (model == "fims") {
-    if (is.character(file) && length(file) == 1) {
-      fims_output <- readRDS(file)
-    } else {
-      fims_output <- file
-    }
-    # TEMPORARY -- add in era to df
-    if ("era" %notin% colnames(fims_output)) {
-      fims_output$era <- NA
-    }
-    fims_output[setdiff(tolower(names(out_new)), tolower(names(fims_output)))] <- NA
-    out_new <- fims_output
+    out_new <- conout_fims(file = file, out_new = out_new, start_year = start_year)
     #### Rceattle ####
   } else if (model == "rceattle") {
    out_new <- conout_rceattle(file = file, out_new = out_new, fleet_names = fleet_names)
diff --git a/R/plot_index.R b/R/plot_index.R
index 4acb0a13..03552cbc 100644
--- a/R/plot_index.R
+++ b/R/plot_index.R
@@ -59,7 +59,9 @@ plot_index <- function(
   # Filter data
   prepared_data <- filter_data(
     dat,
-    label_name = "index",
+    # Only the observed and predicted index are plotted; other labels that
+    # contain "index" (e.g., catchability or index numbers at age) are not
+    label_name = "^index_(observed|predicted)$",
     era = NULL,
     geom = "line",
     # ifelse guarantees the code doesn't miss grouping when label has > 1 value
@@ -90,6 +92,10 @@ plot_index <- function(
   if (!is.null(group)) {
     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
   # transform prep data so group_var = label
   # prepared_data <- prepared_data |>
   #   dplyr::mutate(
@@ -117,7 +123,7 @@ plot_index <- function(
       )
     )
   # Overwrite facets from base plot_timeseries bc scales need to be free
-  if ("fleet" %in% colnames(prepared_data)) {
+  if (!is.null(facet) && "fleet" %in% colnames(prepared_data)) {
     facet <- paste("~", paste(facet, collapse = " + "))
     facet_formula <- stats::reformulate(facet)
     plt <- plt + ggplot2::facet_wrap(facet_formula, scales = "free", labeller = function(labels) {
diff --git a/R/plot_runs_test.R b/R/plot_runs_test.R
new file mode 100644
index 00000000..8fcdf442
--- /dev/null
+++ b/R/plot_runs_test.R
@@ -0,0 +1,348 @@
+#' Plot residual runs tests for indices of abundance
+#'
+#' Plot the residuals of the fit to each index over time with three-sigma
+#' control limits and a runs test for non-random residuals (Carvalho et al.
+#' 2017), as in `JABBA` and `ss3diags` (`SSplotRunstest()`).
+#'
+#' @inheritParams plot_index
+#' @param dat A tibble from [convert_output()] with labels "index_observed"
+#' and "index_predicted".
+#' @param fleet Character vector. Fleet(s) to plot. Default: NULL (all fleets)
+#' @param ncol Number. Number of columns of panels. Default: NULL
+#'
+#' @returns A `ggplot2` object.
+#'
+#' @details Residuals are \eqn{\log(o) - \log(p)} for observed \eqn{o} and
+#' predicted \eqn{p}, assuming lognormal errors. The one-sided runs test
+#' checks whether residual signs are random over time (normal
+#' approximation, as in `ss3diags`). The shaded band spans the three-sigma
+#' limits, with sigma from the average moving range; it is green when the
+#' test passes (p >= 0.05) and red when it fails. Residuals outside the
+#' limits are filled red.
+#'
+#' No observations are dropped silently: a fleet observed in more than one
+#' area, season, subseason, month, or sex becomes separate series (e.g.,
+#' "Survey (season 2)"), observations without a positive observed and
+#' predicted value are reported in a message, and output from more than one
+#' model or with values that cannot be told apart is an error. Years with no
+#' observed index (BAM reports every year) are not observations.
+#'
+#' @references Carvalho, F., et al. 2017. Can diagnostic tests help identify
+#' model misspecification in integrated stock assessments? Fisheries Research
+#' 192: 28-40.
+#'
+#' @seealso [plot_index()], [convert_output()]
+#' @export
+#'
+#' @examples
+#' plot_runs_test(stockplotr::example_data, interactive = FALSE)
+plot_runs_test <- function(
+  dat,
+  fleet = NULL,
+  module = NULL,
+  interactive = TRUE,
+  ncol = NULL,
+  make_rda = FALSE,
+  figures_dir = getwd()
+) {
+  check_data(dat)
+  resid <- index_residuals(dat, fleet = fleet, module = module, interactive = interactive)
+
+  ctl <- resid |>
+    dplyr::arrange(.data$series, .data$year) |>
+    dplyr::group_by(.data$series) |>
+    dplyr::summarise(
+      lcl = residual_control(.data$residual)$lcl,
+      ucl = residual_control(.data$residual)$ucl,
+      p = residual_control(.data$residual)$p,
+      xmin = min(.data$year) - 0.5,
+      xmax = max(.data$year) + 0.5,
+      .groups = "drop"
+    ) |>
+    dplyr::mutate(
+      result = factor(
+        dplyr::case_when(
+          is.na(.data$p) ~ "Too few residuals",
+          .data$p < 0.05 ~ "Failed (p < 0.05)",
+          TRUE ~ "Passed"
+        ),
+        levels = c("Passed", "Failed (p < 0.05)", "Too few residuals")
+      ),
+      panel = paste0(
+        wrap_label(.data$series), "\n",
+        dplyr::if_else(
+          is.na(.data$p),
+          "p = NA",
+          paste0("p = ", format.pval(round(.data$p, 3), eps = 0.001))
+        )
+      )
+    )
+  resid <- resid |>
+    dplyr::left_join(ctl[, c("series", "lcl", "ucl", "panel")], by = "series") |>
+    dplyr::mutate(
+      outside = !is.na(.data$lcl) & (.data$residual < .data$lcl | .data$residual > .data$ucl)
+    )
+
+  plt <- ggplot2::ggplot(resid, ggplot2::aes(x = .data$year, y = .data$residual)) +
+    ggplot2::geom_rect(
+      data = ctl,
+      ggplot2::aes(
+        xmin = .data$xmin, xmax = .data$xmax,
+        ymin = .data$lcl, ymax = .data$ucl,
+        fill = .data$result
+      ),
+      inherit.aes = FALSE,
+      alpha = 0.35,
+      na.rm = TRUE
+    ) +
+    ggplot2::geom_hline(yintercept = 0, linetype = "dashed", color = "gray40") +
+    ggplot2::geom_segment(ggplot2::aes(xend = .data$year, yend = 0), color = "gray40") +
+    ggplot2::geom_point(fill = "white", shape = 21, size = 2, stroke = 0.8) +
+    ggplot2::geom_point(data = resid[resid$outside, ], color = "#B2182B", size = 2) +
+    ggplot2::scale_fill_manual(
+      values = c(
+        "Passed" = "#4DAC26",
+        "Failed (p < 0.05)" = "#D7191C",
+        "Too few residuals" = "gray80"
+      ),
+      drop = TRUE
+    ) +
+    ggplot2::labs(x = "Year", y = "Residual (log observed / predicted)", fill = "Runs test") +
+    theme_noaa() +
+    ggplot2::theme(legend.position = "bottom") +
+    ggplot2::facet_wrap(ggplot2::vars(.data$panel), ncol = ncol, scales = "free")
+
+  if (make_rda) {
+    fit_export_rda(
+      plt,
+      topic_label = "runs_test",
+      tags = runs_test_tags(resid, ctl),
+      kqs = runs_test_kqs(resid, ctl),
+      dat = dat,
+      figures_dir = figures_dir
+    )
+  }
+  plt
+}
+
+#------------------------------------------------------------------------------
+# Index residuals
+#------------------------------------------------------------------------------
+
+# Columns that distinguish observations of a fleet in the same year
+index_series_cols <- c("area", "season", "subseason", "month", "sex")
+
+# Observed and predicted index with log residuals, one row per series and
+# year. A series is a fleet, or a fleet in one area, season, subseason,
+# month, or sex when it has observations in more than one of these.
+index_residuals <- function(dat, fleet = NULL, module = NULL, interactive = TRUE) {
+  idx <- filter_data(
+    dat,
+    label_name = "^index_(observed|predicted)$",
+    module = module,
+    geom = "point",
+    interactive = interactive
+  )
+  models <- unique(stats::na.omit(as.character(idx$model)))
+  if (length(models) > 1) {
+    cli::cli_abort(c(
+      "{.arg dat} has output from {length(models)} models.",
+      "i" = "Diagnostic plots show one model at a time; supply the output of one model."
+    ))
+  }
+  if (!is.null(fleet)) {
+    idx <- idx[idx$fleet %in% fleet, , drop = FALSE]
+    if (nrow(idx) == 0) {
+      cli::cli_abort("No index found for fleet {.val {fleet}}.")
+    }
+  }
+  ids <- intersect(index_series_cols, names(idx))
+  type <- tolower(as.character(idx$uncertainty_label))
+  long <- tibble::tibble(
+    fleet = as.character(idx$fleet),
+    year = as.numeric(idx$year),
+    kind = sub("^index_", "", idx$label),
+    value = as.numeric(idx$estimate),
+    se = dplyr::case_when(
+      type %in% c("se", "sd", "log_sd", "logsd", "stddev") ~ as.numeric(idx$uncertainty),
+      type == "cv" ~ sqrt(log(1 + as.numeric(idx$uncertainty)^2)),
+      TRUE ~ NA_real_
+    )
+  )
+  for (col in ids) long[[col]] <- as.character(idx[[col]])
+  long$fleet[is.na(long$fleet)] <- "All"
+  # Rows that are exact copies are combined; other duplicates are an error
+  long <- dplyr::distinct(long)
+  keys <- c("fleet", "year", "kind", ids)
+  dups <- long |>
+    dplyr::group_by(dplyr::across(dplyr::all_of(keys))) |>
+    dplyr::filter(dplyr::n() > 1) |>
+    dplyr::ungroup() |>
+    dplyr::distinct(.data$fleet, .data$year)
+  if (nrow(dups) > 0) {
+    examples <- utils::head(paste(dups$fleet, dups$year), 3)
+    cli::cli_abort(c(
+      "{.arg dat} has more than one index value for {nrow(dups)} fleet{?s} and year{?s} that cannot be told apart (e.g., {examples}).",
+      "i" = "Filter {.arg dat} to one value per fleet and year, or add a column (e.g., {.field season}) that distinguishes them."
+    ))
+  }
+  wide <- tidyr::pivot_wider(
+    long,
+    id_cols = dplyr::all_of(c("fleet", "year", ids)),
+    names_from = "kind",
+    values_from = c("value", "se")
+  )
+  if (!"value_predicted" %in% names(wide)) wide$value_predicted <- NA_real_
+  if (!"se_observed" %in% names(wide)) wide$se_observed <- NA_real_
+  out <- tibble::tibble(
+    fleet = wide$fleet,
+    year = wide$year,
+    observed = wide$value_observed,
+    predicted = wide$value_predicted,
+    se = wide$se_observed
+  )
+  # Name series after the columns that vary within a fleet
+  extra <- rep("", nrow(out))
+  for (col in ids) {
+    vals <- wide[[col]]
+    varies <- tapply(vals, out$fleet, function(v) length(unique(stats::na.omit(v))) > 1)
+    use <- out$fleet %in% names(varies)[varies %in% TRUE] & !is.na(vals)
+    extra[use] <- paste0(extra[use], ifelse(extra[use] == "", "", ", "), col, " ", vals[use])
+  }
+  out$series <- ifelse(extra == "", out$fleet, paste0(out$fleet, " (", extra, ")"))
+  out <- out[!is.na(out$observed), , drop = FALSE]
+  out$residual <- suppressWarnings(log(out$observed) - log(out$predicted))
+  out$residual[!is.finite(out$residual)] <- NA_real_
+  no_resid <- is.na(out$residual)
+  if (any(no_resid)) {
+    cli::cli_alert_info(
+      "{sum(no_resid)} observation{?s} without a positive observed and predicted index {?has/have} no residual and {?is/are} not shown."
+    )
+  }
+  out <- out[!no_resid, , drop = FALSE]
+  if (nrow(out) == 0) {
+    cli::cli_abort("No index residuals found in {.arg dat}.")
+  }
+  dplyr::arrange(out, .data$series, .data$year)
+}
+
+# Wrap long labels for panel titles at spaces or underscores (fleet names
+# such as "Survey_1_shelikof_acoustic" often have no spaces), keeping every
+# character of the label
+wrap_label <- function(x, width = 28) {
+  vapply(x, function(lab) {
+    if (nchar(lab) <= width) {
+      return(lab)
+    }
+    # Split after each space or underscore, keeping the separator
+    pieces <- regmatches(lab, gregexpr("[^ _]+[ _]?", lab))[[1]]
+    lines <- character(0)
+    current <- ""
+    for (piece in pieces) {
+      if (nchar(current) > 0 && nchar(current) + nchar(piece) > width) {
+        lines <- c(lines, trimws(current, "right"))
+        current <- piece
+      } else {
+        current <- paste0(current, piece)
+      }
+    }
+    paste(c(lines, trimws(current, "right")), collapse = "\n")
+  }, character(1), USE.NAMES = FALSE)
+}
+
+#------------------------------------------------------------------------------
+# Runs test and control limits
+#------------------------------------------------------------------------------
+
+# Three-sigma control limits and runs test for a series of residuals
+# (Carvalho et al. 2017; Winker et al. 2018). The standard deviation is
+# estimated from the average moving range (Montgomery 2009), after removing
+# moving ranges greater than the upper limit for moving ranges.
+residual_control <- function(x) {
+  x <- x[!is.na(x)]
+  if (length(x) < 2) {
+    return(list(lcl = NA_real_, ucl = NA_real_, p = residual_runs_test(x)))
+  }
+  mr <- abs(diff(x))
+  amr <- mean(mr)
+  mr <- mr[mr < 3.267 * amr]
+  stdev <- if (length(mr) > 0) mean(mr) / 1.128 else NA_real_
+  list(
+    lcl = -3 * stdev,
+    ucl = 3 * stdev,
+    p = residual_runs_test(x)
+  )
+}
+
+# One-sided Wald-Wolfowitz runs test on the signs of x (in order). Returns
+# the p-value for fewer runs than expected (positive autocorrelation), using
+# the normal approximation.
+residual_runs_test <- function(x) {
+  s <- sign(x[x != 0 & !is.na(x)])
+  n1 <- sum(s > 0)
+  n2 <- sum(s < 0)
+  n <- n1 + n2
+  if (n < 3) {
+    return(NA_real_)
+  }
+  # A single run: probability that all n signs are the same by chance
+  if (n1 == 0 || n2 == 0) {
+    return(2 * 0.5^n)
+  }
+  runs <- 1 + sum(s[-1] != s[-length(s)])
+  mu <- 2 * n1 * n2 / n + 1
+  v <- 2 * n1 * n2 * (2 * n1 * n2 - n) / (n^2 * (n - 1))
+  if (v <= 0) {
+    return(NA_real_)
+  }
+  stats::pnorm((runs - mu) / sqrt(v))
+}
+
+#------------------------------------------------------------------------------
+# Caption and alternative text tags
+#------------------------------------------------------------------------------
+
+# `resid` has one row per residual and `ctl` one row per series
+runs_test_tags <- function(resid, ctl) {
+  passed <- ctl$series[!is.na(ctl$p) & ctl$p >= 0.05]
+  failed <- ctl$series[!is.na(ctl$p) & ctl$p < 0.05]
+  outside <- resid[resid$outside, , drop = FALSE]
+  outside_text <- paste0(outside$series, " in ", outside$year)
+  if (length(outside_text) > 5) {
+    outside_text <- c(outside_text[1:5], paste(length(outside_text) - 5, "more"))
+  }
+  tags <- list(
+    panel.desc = fit_panel_text(ctl$series, what = "index"),
+    start.year = min(resid$year),
+    end.year = max(resid$year),
+    summary = paste0(
+      "The runs test passed for ", length(passed), " of ", nrow(ctl), " indices",
+      if (length(failed) > 0) paste0(" and failed for ", fit_list_text(failed)), ". ",
+      if (nrow(outside) == 0) {
+        "No residuals are outside the control limits."
+      } else {
+        paste0(
+          if (nrow(outside) == 1) "One residual is" else paste(nrow(outside), "residuals are"),
+          " outside the control limits: ", fit_list_text(outside_text), "."
+        )
+      }
+    )
+  )
+  stats::setNames(tags, paste0("runs.test.", names(tags)))
+}
+
+# Key quantities recorded in key_quantities.csv (runs.test.*)
+runs_test_kqs <- function(resid, ctl) {
+  tested <- ctl[!is.na(ctl$p), , drop = FALSE]
+  failed <- tested$series[tested$p < 0.05]
+  kqs <- list(
+    start.year = min(resid$year),
+    end.year = max(resid$year),
+    n.indices = nrow(ctl),
+    n.passed = sum(tested$p >= 0.05),
+    failed = if (length(failed) == 0) "none" else paste(failed, collapse = "; "),
+    p.values = paste0(ctl$series, ": ", ifelse(is.na(ctl$p), "NA", format(round(ctl$p, 3), nsmall = 3)), collapse = "; "),
+    n.outside = sum(resid$outside)
+  )
+  stats::setNames(kqs, paste0("runs.test.", names(kqs)))
+}
diff --git a/R/utils.R b/R/utils.R
index 4688f4fe..551f8eaa 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -1,6 +1,6 @@
-# @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+############################################################################
 # General utility functions
-# @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
+############################################################################
 
 # Check end_year isn't past current year for non-projections plots
 # make year character if not null
diff --git a/R/utils_fit_captions.R b/R/utils_fit_captions.R
new file mode 100644
index 00000000..307543f5
--- /dev/null
+++ b/R/utils_fit_captions.R
@@ -0,0 +1,114 @@
+# Captions, alternative text, and key quantities for the fit diagnostic
+# plots. Each plot has a row in captions_alt_text_template.csv whose
+# placeholders start with the plot's topic label (e.g., runs.test.summary),
+# and rows in key_quantity_template.csv for the quantities it records (e.g.,
+# runs.test.n.passed). The plot describes what it shows as tags, which
+# fit_export_rda() inserts with insert_kqs(), and its results as key
+# quantities, which fit_export_rda() records with export_kqs().
+
+# Export key quantities, insert tags into captions_alt_text.csv, and export
+# the figure as an rda
+fit_export_rda <- function(plt, topic_label, tags, dat, figures_dir, kqs = list()) {
+  if (length(kqs) > 0) {
+    fit_add_kq_rows(names(kqs), dir = figures_dir)
+    fit_call_with_names(export_kqs, lapply(kqs, fit_tag_value), dir = figures_dir)
+  }
+  fit_add_caption_rows(topic_label, dir = figures_dir)
+  fit_call_with_names(insert_kqs, lapply(tags, fit_tag_value), dir = figures_dir)
+  create_rda(
+    object = plt,
+    topic_label = topic_label,
+    fig_or_table = "figure",
+    dat = dat,
+    dir = figures_dir
+  )
+  invisible(tags)
+}
+
+# Call export_kqs() or insert_kqs(), which take the names of the quantities
+# from the names of their arguments, with one variable per value
+fit_call_with_names <- function(fun, values, dir) {
+  env <- list2env(values, parent = environment())
+  do.call(fun, c(rlang::syms(names(values)), list(dir = dir)), envir = env)
+}
+
+# Add the template rows for key quantities to an existing key_quantities.csv
+# that was made before they existed (e.g., by an older version of
+# stockplotr), so that their values are recorded
+fit_add_kq_rows <- function(kq_names, dir = getwd()) {
+  path <- fs::path(dir, "key_quantities.csv")
+  if (!file.exists(path)) {
+    return(invisible(NULL))
+  }
+  kqs <- utils::read.csv(path)
+  missing <- setdiff(kq_names, kqs$key_quantity)
+  if (length(missing) == 0) {
+    return(invisible(NULL))
+  }
+  template <- utils::read.csv(
+    system.file("resources", "key_quantity_template.csv", package = "stockplotr")
+  )
+  utils::write.csv(
+    rbind(kqs, template[template$key_quantity %in% missing, , drop = FALSE]),
+    path,
+    row.names = FALSE
+  )
+  invisible(NULL)
+}
+
+# Add the template rows for a topic to an existing captions_alt_text.csv
+# that was made before the topic existed (e.g., by an older version of
+# stockplotr), so that its caption and alternative text are filled in
+fit_add_caption_rows <- function(topic_label, dir = getwd()) {
+  path <- fs::path(dir, "captions_alt_text.csv")
+  if (!file.exists(path)) {
+    return(invisible(NULL))
+  }
+  caps <- utils::read.csv(path)
+  if (any(caps$label == topic_label & caps$type == "figure")) {
+    return(invisible(NULL))
+  }
+  template <- utils::read.csv(
+    system.file("resources", "captions_alt_text_template.csv", package = "stockplotr")
+  )
+  new_rows <- template[template$label == topic_label & template$type == "figure", , drop = FALSE]
+  utils::write.csv(rbind(caps, new_rows), path, row.names = FALSE)
+  invisible(NULL)
+}
+
+# One string per tag. Remove characters with special meaning in the
+# replacement strings of stringr::str_replace_all().
+fit_tag_value <- function(x) {
+  x <- paste(x, collapse = " ")
+  gsub("[$\\\\]", "", x)
+}
+
+#------------------------------------------------------------------------------
+# Text helpers
+#------------------------------------------------------------------------------
+
+# "A", "A and B", or "A, B, and C"
+fit_list_text <- function(x) {
+  x <- unique(as.character(x))
+  n <- length(x)
+  if (n == 0) {
+    return("")
+  }
+  if (n <= 2) {
+    return(paste(x, collapse = " and "))
+  }
+  paste0(paste(x[-n], collapse = ", "), ", and ", x[n])
+}
+
+# Description of the panels in a faceted plot
+fit_panel_text <- function(panels, what = "fleet") {
+  panels <- unique(panels)
+  if (length(panels) == 1) {
+    paste0("A single panel for ", panels, ".")
+  } else {
+    paste0(
+      "A set of ", length(panels), " panels, one for each ", what,
+      " (", fit_list_text(panels), ")."
+    )
+  }
+}
diff --git a/inst/WORDLIST b/inst/WORDLIST
index c725269d..a9a58810 100644
--- a/inst/WORDLIST
+++ b/inst/WORDLIST
@@ -3,6 +3,7 @@ BMSY
 Barbeaux
 Birthseason
 CPUE
+Carvalho
 Dorn
 FIMS
 FLIMIT
@@ -29,6 +30,7 @@ SPR
 Spawners
 StockSmart
 Subseason
+al
 alttext
 asar
 assmt
@@ -44,6 +46,7 @@ dat
 dataframe's
 diagonostic
 dplyr
+et
 facetting
 fims
 flextable
@@ -56,6 +59,7 @@ is’
 kableExtra
 kobe
 len
+misspecification
 msy
 natively
 noaa
@@ -64,6 +68,7 @@ pkgdown
 preformatted
 rceattle
 rda
+rds
 rda's
 rdat
 rowid
diff --git a/inst/resources/captions_alt_text_template.csv b/inst/resources/captions_alt_text_template.csv
index e5cfabab..5ea9e619 100644
--- a/inst/resources/captions_alt_text_template.csv
+++ b/inst/resources/captions_alt_text_template.csv
@@ -48,3 +48,4 @@ catchability,table,Catchability over time for fleet catchability.fleet.,
 sensitivity.runs,table,"Forecasted catch, biomass, spawning biomass, recruitment, fishing mortality, and stock status as projected from the base model configuration.",
 catch_comp,figure,Fishery age composition (caa.start.year-caa.end.year). The area of the circle is proportional to the catch. Diagonal lines indicated the top 5\\% strongest year classes.,"Bubble plot showing estimated total catch at age. The x axis shows the year, which spans from caa.start.year to caa.end.year. The y axis shows age, which spans from caa.age.min to caa.age.max. The size of the bubbles range from tot.catch.min to tot.catch.max."
 age_to_length_conversion,figure,"Age-to-length conversion shown as the probability distribution of length bins at each age. Density curves show the proportion of fish in each length bin at a given age, with wider curves indicating greater variability in length at age.","A set of density plots showing the age-to-length conversion for species, with one panel for each morph. The x axis shows age in years and the y axis shows length bins (in length.units). At each age, a vertical density curve shows the distribution of lengths. Mean length increases from about length.at.min.age length.units at age min.age to about length.at.max.age length.units at age max.age, after which the distributions stabilize, indicating that growth slows as fish age. "
+runs_test,figure,"Residuals of the fit to indices of abundance over time, with three-sigma control limits (shaded) and a runs test for non-random residuals. The shading is green when the residuals pass the runs test (p >= 0.05) and red when they fail, and points filled in red are outside the control limits. Residuals are the log of the observed index divided by the predicted index.","runs.test.panel.desc The x axis shows year, from runs.test.start.year to runs.test.end.year, and the y axis shows the residual, with a vertical line from zero to each residual. runs.test.summary"
diff --git a/inst/resources/fims_var_names.csv b/inst/resources/fims_var_names.csv
index 3d295826..7b3c2cd1 100644
--- a/inst/resources/fims_var_names.csv
+++ b/inst/resources/fims_var_names.csv
@@ -1,46 +1,59 @@
-module_name,label,alt_label
-Fleet,age_to_length_conversion,
-Fleet,agecomp_expected,age_composition_expected
-Fleet,agecomp_proportion,age_composition_proportion
-Fleet,index_expected,
-Fleet,index_naa,index
-Fleet,index_nal,index
-Fleet,index_numbers,
-Fleet,index_waa,index_weight
-Fleet,index_weight,
-Fleet,landings_expected,
-Fleet,landings_naa,landings
-Fleet,landings_nal,landings
-Fleet,landings_numbers,landings_numbers
-Fleet,landings_waa,landings_weight
-Fleet,landings_weight,landings_weight
-Fleet,lengthcomp_expected,length_composition_expected
-Fleet,lengthcomp_proportion,length_composition_proportion
-Fleet,log_Fmort,log_fishing_mortality
-Fleet,log_q,log_catchability
-Population,Biomass,biomass
-Population,NAA,abundance
-Population,Recruitment,recruitment
-Population,SSB,spawning_biomass
-Population,log_M,log_natural_mortality
-Population,log_init_naa,log_abundance_init
-maturity,inflection_point,
-maturity,slope,
-recruitment,log_devs,log_recruitment_deviations
-recruitment,log_rzero,log_recruitment_unfished
-recruitment,logit_steep,
-selectivity,inflection_point,
-selectivity,slope,
-,AgeCompositionExpected,age_composition_expected
-,AgeCompositionProportion,age_composition_proportion
-,FMort,fishing_mortality
-,IndexExpected,index_expected
-,IndexNumberAtAge,index_numbers
-,IndexNumberAtLength,index_numbers
-,LandingsExpected,landings_expected
-,LandingsNumberAtAge,landings_numbers
-,LandingsNumberAtLength,landings_numbers
-,LengthCompositionExpected,length_composition_expected
-,LengthCompositionProportion,length_composition_proportion
-,LogRecDev,log_recruitment_deviations
-,Q,catchability
+module_name,label,alt_label,unit
+Fleet,age_to_length_conversion,,
+Fleet,agecomp_expected,age_composition_expected,
+Fleet,agecomp_proportion,age_composition_proportion,
+Fleet,index_expected,,
+Fleet,index_naa,index,numbers
+Fleet,index_nal,index,numbers
+Fleet,index_numbers,,numbers
+Fleet,index_waa,index_weight,mt
+Fleet,index_weight,,mt
+Fleet,landings_expected,,
+Fleet,landings_naa,landings,numbers
+Fleet,landings_nal,landings,numbers
+Fleet,landings_numbers,landings_numbers,numbers
+Fleet,landings_waa,landings_weight,mt
+Fleet,landings_weight,landings_weight,mt
+Fleet,lengthcomp_expected,length_composition_expected,
+Fleet,lengthcomp_proportion,length_composition_proportion,
+Fleet,log_Fmort,log_fishing_mortality,
+Fleet,log_q,log_catchability,
+Population,Biomass,biomass,mt
+Population,NAA,abundance,numbers
+Population,Recruitment,recruitment,numbers
+Population,SSB,spawning_biomass,mt
+Population,log_M,log_natural_mortality,
+Population,log_init_naa,log_abundance_init,
+maturity,inflection_point,,
+maturity,slope,,
+recruitment,log_devs,log_recruitment_deviations,
+recruitment,log_rzero,log_recruitment_unfished,
+recruitment,logit_steep,,
+selectivity,inflection_point,,
+selectivity,slope,,
+,AgeCompositionExpected,age_composition_expected,
+,AgeCompositionProportion,age_composition_proportion,
+,FMort,fishing_mortality,
+,IndexExpected,index_expected,
+,IndexNumberAtAge,index_numbers,
+,IndexNumberAtLength,index_numbers,
+,LandingsExpected,landings_expected,
+,LandingsNumberAtAge,landings_numbers,
+,LandingsNumberAtLength,landings_numbers,
+,LengthCompositionExpected,length_composition_expected,
+,LengthCompositionProportion,length_composition_proportion,
+,LogRecDev,log_recruitment_deviations,
+,Q,catchability,
+Population,expected_recruitment,recruitment,numbers
+Population,numbers_at_age,abundance,numbers
+Population,unfished_spawning_biomass,spawning_biomass_unfished,mt
+Population,unfished_biomass,biomass_unfished,mt
+Population,mortality_M,natural_mortality,
+Recruitment,log_devs,log_recruitment_deviations,
+Recruitment,log_rzero,log_recruitment_unfished,
+Recruitment,log_sd,log_recruitment_sd,
+Recruitment,logit_steep,logit_steepness,
+Population,biomass,,mt
+Population,spawning_biomass,,mt
+Fleet,landings_observed,,mt
+Fleet,landings_predicted,,mt
diff --git a/inst/resources/key_quantity_template.csv b/inst/resources/key_quantity_template.csv
index 04fb14ae..c6d201ae 100644
--- a/inst/resources/key_quantity_template.csv
+++ b/inst/resources/key_quantity_template.csv
@@ -112,6 +112,13 @@ rel.recruitment.max,,maximum relative recruitment (R/R0),,,,
 rel.recruitment.min,,minimum relative recruitment (R/R0),,,,
 rel.ssb.max,,maximum relative spawning biomass,,,,
 rel.ssb.min,,minimum relative spawning biomass,,,,
+runs.test.end.year,,final year of the index residuals in the runs test plot,,,,
+runs.test.failed,,indices whose residuals fail the runs test (p < 0.05),,,,
+runs.test.n.indices,,number of indices in the runs test plot,,,,
+runs.test.n.outside,,number of index residuals outside the three-sigma control limits,,,,
+runs.test.n.passed,,number of indices whose residuals pass the runs test (p >= 0.05),,,,
+runs.test.p.values,,p-value of the runs test for each index,,,,
+runs.test.start.year,,start year of the index residuals in the runs test plot,,,,
 selectivity.end.year,,final year for selectivity estimates,,,,
 selectivity.start.year,,start year for selectivity estimates,,,,
 spr.start.year,,start year of SPR timeseries,,,,
diff --git a/man/conout_bam.Rd b/man/conout_bam.Rd
index 7b3ee972..8db1ca02 100644
--- a/man/conout_bam.Rd
+++ b/man/conout_bam.Rd
@@ -8,7 +8,8 @@ conout_bam(file, out_new, fleet_names = NULL)
 }
 \arguments{
 \item{file}{Path. Assessment model output file path. Currently accepts SS3
-\code{r4ss::ss_output()} or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS model object.}
+\code{r4ss::ss_output()} or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS output
+(the tibble returned by \code{FIMS::get_estimates()}, a \code{FIMSFit} object, or an .rds file with either).}
 
 \item{fleet_names}{Character vector. Names of fleets in the assessment model as
 shortened in the output file or names of fleets to replace the native fleet
diff --git a/man/conout_rceattle.Rd b/man/conout_rceattle.Rd
index 85612b51..e053b2a8 100644
--- a/man/conout_rceattle.Rd
+++ b/man/conout_rceattle.Rd
@@ -8,7 +8,8 @@ conout_rceattle(file, out_new, fleet_names = NULL)
 }
 \arguments{
 \item{file}{Path. Assessment model output file path. Currently accepts SS3
-\code{r4ss::ss_output()} or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS model object.}
+\code{r4ss::ss_output()} or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS output
+(the tibble returned by \code{FIMS::get_estimates()}, a \code{FIMSFit} object, or an .rds file with either).}
 
 \item{fleet_names}{Character vector. Names of fleets in the assessment model as
 shortened in the output file or names of fleets to replace the native fleet
diff --git a/man/conout_ss3.Rd b/man/conout_ss3.Rd
index fe5de1e8..303b7fbf 100644
--- a/man/conout_ss3.Rd
+++ b/man/conout_ss3.Rd
@@ -8,7 +8,8 @@ conout_ss3(file, out_new, fleet_names = NULL)
 }
 \arguments{
 \item{file}{Path. Assessment model output file path. Currently accepts SS3
-\code{r4ss::ss_output()} or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS model object.}
+\code{r4ss::ss_output()} or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS output
+(the tibble returned by \code{FIMS::get_estimates()}, a \code{FIMSFit} object, or an .rds file with either).}
 
 \item{fleet_names}{Character vector. Names of fleets in the assessment model as
 shortened in the output file or names of fleets to replace the native fleet
diff --git a/man/convert_output.Rd b/man/convert_output.Rd
index 782b50e7..b409306e 100644
--- a/man/convert_output.Rd
+++ b/man/convert_output.Rd
@@ -4,11 +4,18 @@
 \alias{convert_output}
 \title{Convert Output}
 \usage{
-convert_output(file, model = NULL, fleet_names = NULL, save_dir = NULL)
+convert_output(
+  file,
+  model = NULL,
+  fleet_names = NULL,
+  save_dir = NULL,
+  start_year = NULL
+)
 }
 \arguments{
 \item{file}{Path. Assessment model output file path. Currently accepts SS3
-\code{r4ss::ss_output()} or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS model object.}
+\code{r4ss::ss_output()} or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS output
+(the tibble returned by \code{FIMS::get_estimates()}, a \code{FIMSFit} object, or an .rds file with either).}
 
 \item{model}{String. Assessment model used in evaluation
 
@@ -28,6 +35,12 @@ saved file name is not indicated, the file will be saved as "std_output.rda"
 in the file path provided.
 
 Default: NULL}
+
+\item{start_year}{Number. First year of the model, for FIMS output only.
+FIMS reports years as indices (1 = first model year), which are converted
+to years when \code{start_year} is given.
+
+Default: NULL (years are left as indices)}
 }
 \value{
 A reformatted and standardized tibble of assessment model results
@@ -41,6 +54,14 @@ Format stock assessment output file to a standardized tibble.
 The resulting object is a long-format, transformed and machine readable
 version of a model output file. Converted data frame is always returned.
 It will also be saved if save_dir is not NULL.
+
+For FIMS, each data point (index, landings, and age and length
+compositions) is split into an observed row (e.g., "index_observed") and a
+predicted row (e.g., "index_predicted"), as for the other models. Observed
+rows hold the observation error (a log-scale standard error for lognormal
+data) in \code{uncertainty} and the negative log-likelihood in \code{likelihood}.
+Other estimates are relabeled with \code{fims_var_names.csv} (e.g.,
+"expected_recruitment" becomes "recruitment").
 }
 \examples{
 \dontrun{
diff --git a/man/plot_runs_test.Rd b/man/plot_runs_test.Rd
new file mode 100644
index 00000000..a3687cac
--- /dev/null
+++ b/man/plot_runs_test.Rd
@@ -0,0 +1,83 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_runs_test.R
+\name{plot_runs_test}
+\alias{plot_runs_test}
+\title{Plot residual runs tests for indices of abundance}
+\usage{
+plot_runs_test(
+  dat,
+  fleet = NULL,
+  module = NULL,
+  interactive = TRUE,
+  ncol = NULL,
+  make_rda = FALSE,
+  figures_dir = getwd()
+)
+}
+\arguments{
+\item{dat}{A tibble from \code{\link[=convert_output]{convert_output()}} with labels "index_observed"
+and "index_predicted".}
+
+\item{fleet}{Character vector. Fleet(s) to plot. Default: NULL (all fleets)}
+
+\item{module}{Character vector. (Optional) Module name found in \code{dat$module_name}.
+If selecting >1 module, place them in a vector like c("module1", "module2").
+
+Default: NULL
+
+If the interactive and >1 module_name is found, user will select the
+module_name in the console. @seealso \code{\link[=filter_data]{filter_data()}}}
+
+\item{interactive}{Logical. TRUE/FALSE; indicate whether the environment is interactive.
+
+Default: \code{FALSE}}
+
+\item{ncol}{Number. Number of columns of panels. Default: NULL}
+
+\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,
+the rda will be exported to the folder indicated in the argument "figures_dir".
+
+Default: \code{FALSE}.}
+
+\item{figures_dir}{Path. Path to the "figures" folder.
+
+Default: \code{getwd()}
+
+The folder is created within the path if it does not exist.}
+}
+\value{
+A \code{ggplot2} object.
+}
+\description{
+Plot the residuals of the fit to each index over time with three-sigma
+control limits and a runs test for non-random residuals (Carvalho et al.
+2017), as in \code{JABBA} and \code{ss3diags} (\code{SSplotRunstest()}).
+}
+\details{
+Residuals are \eqn{\log(o) - \log(p)} for observed \eqn{o} and
+predicted \eqn{p}, assuming lognormal errors. The one-sided runs test
+checks whether residual signs are random over time (normal
+approximation, as in \code{ss3diags}). The shaded band spans the three-sigma
+limits, with sigma from the average moving range; it is green when the
+test passes (p >= 0.05) and red when it fails. Residuals outside the
+limits are filled red.
+
+No observations are dropped silently: a fleet observed in more than one
+area, season, subseason, month, or sex becomes separate series (e.g.,
+"Survey (season 2)"), observations without a positive observed and
+predicted value are reported in a message, and output from more than one
+model or with values that cannot be told apart is an error. Years with no
+observed index (BAM reports every year) are not observations.
+}
+\examples{
+plot_runs_test(stockplotr::example_data, interactive = FALSE)
+}
+\references{
+Carvalho, F., et al. 2017. Can diagnostic tests help identify
+model misspecification in integrated stock assessments? Fisheries Research
+192: 28-40.
+}
+\seealso{
+\code{\link[=plot_index]{plot_index()}}, \code{\link[=convert_output]{convert_output()}}
+}
diff --git a/tests/testthat/fixtures/fims_output/README.md b/tests/testthat/fixtures/fims_output/README.md
new file mode 100644
index 00000000..a16f69d9
--- /dev/null
+++ b/tests/testthat/fixtures/fims_output/README.md
@@ -0,0 +1,7 @@
+# FIMS fixture
+
+`fims_estimates.rds` is a subset of `FIMS::estimates_big_optimized` (the output of
+`FIMS::get_estimates()` for the `data_big` example model shipped with FIMS,
+NOAA-FIMS/FIMS commit 506b290). Rows kept: every data row (non-missing
+`observed` or `likelihood`), all fixed and random effects, and the spawning
+biomass, recruitment, biomass, index, and catch derived quantities.
diff --git a/tests/testthat/fixtures/fims_output/fims_estimates.rds b/tests/testthat/fixtures/fims_output/fims_estimates.rds
new file mode 100644
index 00000000..e2f03011
Binary files /dev/null and b/tests/testthat/fixtures/fims_output/fims_estimates.rds differ
diff --git a/tests/testthat/helper-fit_captions.R b/tests/testthat/helper-fit_captions.R
new file mode 100644
index 00000000..8d6d8306
--- /dev/null
+++ b/tests/testthat/helper-fit_captions.R
@@ -0,0 +1,38 @@
+# Evaluate plotting code with `make_rda = TRUE` in a temporary directory and
+# return the saved rda (a list with the figure, caption, and alt text), with
+# the key quantities recorded in key_quantities.csv as `kqs` (a named
+# vector of the filled-in values).
+# `captions` is an optional captions_alt_text.csv to start from, and `dir` the
+# figures_dir the plot was saved to (default: the temporary directory).
+fit_rda <- function(code, topic_label, captions = NULL, dir = NULL) {
+  wd <- tempfile("fit_rda")
+  dir.create(wd)
+  old <- setwd(wd)
+  on.exit(
+    {
+      setwd(old)
+      unlink(wd, recursive = TRUE)
+    },
+    add = TRUE
+  )
+  if (!is.null(captions)) {
+    utils::write.csv(captions, "captions_alt_text.csv", row.names = FALSE)
+  }
+  suppressMessages(invisible(utils::capture.output(force(code))))
+  rda <- NULL
+  out <- if (is.null(dir)) wd else dir
+  load(file.path(out, "figures", paste0(topic_label, "_figure.rda")))
+  kqs <- utils::read.csv(file.path(out, "key_quantities.csv"))
+  kqs <- kqs[!is.na(kqs$value) & kqs$value != "", , drop = FALSE]
+  rda$kqs <- stats::setNames(as.character(kqs$value), kqs$key_quantity)
+  rda
+}
+
+# Expect that every placeholder starting with `prefix` (e.g., "runs.test.")
+# was replaced in the caption and alt text
+expect_tags_filled <- function(rda, prefix) {
+  expect_s3_class(rda$figure, "ggplot")
+  expect_no_match(c(rda$caption, rda$alt_text), prefix, fixed = TRUE)
+  expect_gt(nchar(rda$caption), 50)
+  expect_gt(nchar(rda$alt_text), 100)
+}
diff --git a/tests/testthat/test-convert_output_fims.R b/tests/testthat/test-convert_output_fims.R
new file mode 100644
index 00000000..1b25ecc7
--- /dev/null
+++ b/tests/testthat/test-convert_output_fims.R
@@ -0,0 +1,143 @@
+# FIMS::get_estimates() output for the FIMS data_big example
+# (see fixtures/fims_output/README.md)
+fims_estimates <- function() {
+  readRDS(test_path("fixtures", "fims_output", "fims_estimates.rds"))
+}
+convert_fims <- function(dat = fims_estimates(), ...) {
+  suppressMessages(convert_output(dat, model = "fims", ...))
+}
+
+test_that("FIMS data points are split into observed and predicted rows", {
+  est <- fims_estimates()
+  out <- convert_fims(est, start_year = 2001)
+  counts <- table(out$label[grepl("_(observed|predicted)$", out$label)])
+  expect_equal(
+    as.vector(counts[c("index_observed", "index_predicted", "landings_observed", "landings_predicted")]),
+    c(30, 30, 30, 30)
+  )
+  expect_equal(
+    as.vector(counts[c("age_composition_observed", "length_composition_observed")]),
+    c(sum(est$label == "agecomp_expected" & !is.na(est$observed)), sum(est$label == "lengthcomp_expected" & !is.na(est$observed)))
+  )
+
+  rows <- est[est$label == "index_expected" & !is.na(est$observed), ]
+  obs <- out[out$label == "index_observed", ]
+  pred <- out[out$label == "index_predicted", ]
+  expect_equal(obs$estimate, rows$observed)
+  expect_equal(pred$estimate, rows$estimated)
+  # Observation error and negative log-likelihood of each observation
+  # Lognormal observation error is reported as a CV
+  expect_equal(obs$uncertainty, sqrt(exp(exp(rows$log_sd)^2) - 1))
+  expect_equal(unique(obs$uncertainty_label), "cv")
+  expect_equal(obs$likelihood, -rows$likelihood)
+  expect_equal(unique(obs$fleet), "survey1")
+  expect_equal(obs$year, rows$year_i + 2000)
+
+  # Age compositions keep their bins and counts
+  comp <- out[out$label == "age_composition_observed", ]
+  expect_true(all(!is.na(comp$age)))
+  expect_equal(sum(comp$estimate), sum(est$observed[est$label == "agecomp_expected"], na.rm = TRUE))
+  expect_true(all(!is.na(out$length_bins[out$label == "length_composition_observed"])))
+})
+
+test_that("FIMS expected values for fleets without data keep their FIMS label", {
+  out <- convert_fims()
+  # fleet1 is the fishery, which has no index
+  expect_false(any(out$label == "index_predicted" & out$fleet == "fleet1"))
+  expect_true(any(out$label == "index_expected" & out$fleet == "fleet1"))
+})
+
+test_that("FIMS observations coded as missing are not observations", {
+  est <- fims_estimates()
+  idx <- which(est$label == "index_expected" & !is.na(est$observed))
+  est$observed[idx[1:2]] <- -999
+  out <- convert_fims(est)
+  expect_equal(sum(out$label == "index_observed"), length(idx) - 2)
+  expect_equal(sum(out$label == "index_predicted"), length(idx))
+})
+
+test_that("FIMS lognormal intervals are on the original scale", {
+  est <- fims_estimates()
+  est <- est[est$label == "landings_expected" | est$label == "catch_expected", ]
+  row <- which(!is.na(est$observed))[1]
+  est$observed[row] <- 100
+  est$log_sd[row] <- log(0.2)
+  out <- convert_fims(est)
+  obs <- out[out$label == "landings_observed", ][1, ]
+  expect_equal(obs$uncertainty_label, "cv")
+  # The CV gives back the lognormal sigma, and so the lognormal interval
+  sigma <- sqrt(log(1 + obs$uncertainty^2))
+  expect_equal(sigma, 0.2)
+  expect_equal(100 * exp(c(-1.96, 1.96) * sigma), c(67.57, 147.99), tolerance = 1e-4)
+})
+
+test_that("FIMS estimates that could not be calculated are NA", {
+  est <- fims_estimates()
+  ssb <- which(est$label == "spawning_biomass")
+  est$estimated[ssb[1]] <- -999
+  est$uncertainty[ssb[1]] <- -999
+  out <- convert_fims(est)
+  out_ssb <- out[out$label == "spawning_biomass", ]
+  # The row is kept, with a missing estimate
+  expect_equal(nrow(out_ssb), length(ssb))
+  expect_true(is.na(out_ssb$estimate[1]) && is.na(out_ssb$uncertainty[1]))
+  expect_false(any(out$estimate %in% -999))
+})
+
+test_that("FIMS estimates are relabeled and get years and eras", {
+  est <- fims_estimates()
+  out <- convert_fims(est, start_year = 2001)
+  ssb <- out[out$label == "spawning_biomass", ]
+  expect_equal(ssb$estimate, est$estimated[est$label == "spawning_biomass"])
+  expect_equal(ssb$uncertainty, est$uncertainty[est$label == "spawning_biomass"])
+  expect_equal(range(ssb$year), c(2001, 2031))
+  # Spawning biomass at the start of the year after the model ends
+  expect_equal(ssb$era[ssb$year == 2031], "fore")
+  expect_true(all(ssb$era[ssb$year <= 2030] == "time"))
+
+  expect_equal(sum(out$label == "recruitment"), sum(est$label == "expected_recruitment"))
+  expect_equal(sum(out$label == "log_recruitment_deviations"), sum(est$label == "log_devs"))
+  expect_false(any(c("expected_recruitment", "log_devs") %in% out$label))
+
+  # Without start_year, years are the FIMS year indices
+  expect_message(convert_output(est, model = "fims"), "start_year")
+  expect_equal(range(convert_fims(est)$year[convert_fims(est)$label == "spawning_biomass"]), c(1, 31))
+})
+
+test_that("FIMS output is read from an rds file", {
+  path <- tempfile(fileext = ".rds")
+  saveRDS(fims_estimates(), path)
+  on.exit(unlink(path), add = TRUE)
+  expect_equal(nrow(convert_fims(path)), nrow(convert_fims()))
+})
+
+test_that("a FIMSFit object needs FIMS", {
+  skip_if(requireNamespace("FIMS", quietly = TRUE))
+  methods::setClass("FIMSFit", representation(x = "numeric"), where = environment())
+  expect_error(
+    suppressMessages(convert_output(methods::new("FIMSFit", x = 1))),
+    "FIMS"
+  )
+})
+
+test_that("plots work with converted FIMS output", {
+  out <- convert_fims(start_year = 2001)
+  p <- suppressMessages(plot_index(out, interactive = FALSE))
+  b <- ggplot2::ggplot_build(p)
+  expect_equal(nrow(b$data[[1]]), 30)
+  expect_equal(range(b$data[[1]]$x), c(2001, 2030))
+  for (f in list(plot_landings, plot_spawning_biomass, plot_biomass, plot_recruitment, plot_recruitment_deviations)) {
+    p <- suppressMessages(suppressWarnings(f(out, interactive = FALSE)))
+    expect_s3_class(p, "ggplot")
+    expect_gt(sum(sapply(ggplot2::ggplot_build(p)$data, nrow)), 0)
+  }
+})
+
+test_that("plot_index plots only observed and predicted indices", {
+  # pollock (Rceattle) has other labels that contain "index", such as
+  # index_log_q, and ss3_simple has a single index fleet
+  for (dat in list(stockplotr::pollock, stockplotr::ss3_simple, stockplotr::example_data)) {
+    p <- suppressMessages(plot_index(dat, interactive = FALSE))
+    expect_no_error(ggplot2::ggplot_build(p))
+  }
+})
diff --git a/tests/testthat/test-plot_runs_test.R b/tests/testthat/test-plot_runs_test.R
new file mode 100644
index 00000000..47e52745
--- /dev/null
+++ b/tests/testthat/test-plot_runs_test.R
@@ -0,0 +1,127 @@
+fims_converted <- function() {
+  est <- readRDS(test_path("fixtures", "fims_output", "fims_estimates.rds"))
+  suppressMessages(convert_output(est, model = "fims", start_year = 2001))
+}
+
+test_that("residual_control estimates limits and the runs test", {
+  ctl <- residual_control(c(1, -1, 1, -1, 1, -1))
+  # All moving ranges are 2, so sigma = 2 / 1.128
+  expect_equal(ctl$ucl, 3 * 2 / 1.128)
+  expect_equal(ctl$lcl, -ctl$ucl)
+  # Alternating signs have the most runs; long runs of one sign the fewest
+  expect_gt(ctl$p, 0.9)
+  expect_lt(residual_control(c(-3:-1, 1:3))$p, 0.05)
+  expect_true(is.na(residual_control(1)$lcl))
+  expect_equal(residual_runs_test(rep(1, 5)), 2 * 0.5^5)
+  # Normal approximation: 6 positive, 4 negative, 5 runs
+  x <- c(1, 1, -1, 1, -1, -1, -1, 1, 1, 1)
+  mu <- 2 * 6 * 4 / 10 + 1
+  v <- 2 * 6 * 4 * (2 * 6 * 4 - 10) / (10^2 * 9)
+  expect_equal(residual_runs_test(x), stats::pnorm((5 - mu) / sqrt(v)))
+})
+
+test_that("index residuals match each example model output", {
+  datasets <- list(
+    example_data = stockplotr::example_data,
+    pollock = stockplotr::pollock,
+    ss3_simple = stockplotr::ss3_simple,
+    vermilion_snapper = stockplotr::vermilion_snapper,
+    fims = fims_converted()
+  )
+  for (name in names(datasets)) {
+    dat <- datasets[[name]]
+    # ss3_simple also has catch labeled as an index in the "catch" module
+    module <- if (name == "ss3_simple") "cpue" else NULL
+    res <- suppressMessages(index_residuals(dat, module = module, interactive = FALSE))
+    rows <- dat[dat$label %in% c("index_observed", "index_predicted"), ]
+    if (!is.null(module)) rows <- rows[rows$module_name == module, ]
+    obs <- rows[rows$label == "index_observed" & !is.na(rows$estimate), ]
+    pred <- rows[rows$label == "index_predicted", ]
+    obs <- obs[order(obs$fleet, obs$year), ]
+    pred <- pred[match(paste(obs$fleet, obs$year), paste(pred$fleet, pred$year)), ]
+    expect_equal(nrow(res), nrow(obs), label = name)
+    expect_equal(res$residual, log(obs$estimate) - log(pred$estimate), label = name)
+    # Log-scale standard errors, with CVs (BAM and FIMS) converted
+    se <- ifelse(obs$uncertainty_label == "cv", sqrt(log(1 + obs$uncertainty^2)), obs$uncertainty)
+    expect_equal(res$se, se, label = name)
+
+    p <- suppressMessages(plot_runs_test(dat, module = module, interactive = FALSE))
+    expect_no_error(ggplot2::ggplot_build(p))
+  }
+})
+
+test_that("index residuals keep every observation", {
+  two_seasons <- data.frame(
+    label = rep(c("index_observed", "index_predicted"), each = 6),
+    fleet = "Survey",
+    year = rep(rep(2001:2003, 2), 2),
+    season = rep(rep(1:2, each = 3), 2),
+    estimate = c(1, 2, 3, 10, 20, 30, 1.1, 2.1, 2.9, 11, 19, 31),
+    uncertainty = 0.2,
+    uncertainty_label = "se"
+  )
+  res <- index_residuals(two_seasons, interactive = FALSE)
+  expect_setequal(unique(res$series), c("Survey (season 1)", "Survey (season 2)"))
+  # Exact copies are combined; values that cannot be told apart are an error
+  expect_equal(nrow(index_residuals(rbind(two_seasons, two_seasons), interactive = FALSE)), 6)
+  two_seasons$season <- NULL
+  expect_error(index_residuals(two_seasons, interactive = FALSE), "cannot be told apart")
+
+  # Observations without a predicted value are reported
+  no_pred <- data.frame(
+    label = c(rep("index_observed", 3), rep("index_predicted", 2)),
+    fleet = "a", year = c(1:3, 1:2), estimate = c(1, 2, 3, 1, 2),
+    uncertainty = NA, uncertainty_label = NA
+  )
+  expect_message(res <- index_residuals(no_pred, interactive = FALSE), "1 observation")
+  expect_equal(res$year, 1:2)
+
+  expect_error(
+    index_residuals(list(stockplotr::example_data, stockplotr::example_data), interactive = FALSE),
+    "2 models"
+  )
+  expect_error(index_residuals(stockplotr::example_data, fleet = "none", interactive = FALSE), "No index found")
+})
+
+test_that("plot_runs_test flags outliers and failed runs tests", {
+  p <- plot_runs_test(stockplotr::example_data, interactive = FALSE)
+  # The Triennial index in 2004 is outside the control limits
+  expect_equal(p$data$year[p$data$outside], 2004)
+  expect_equal(p$data$fleet[p$data$outside], "Triennial")
+
+  p <- suppressMessages(plot_runs_test(stockplotr::pollock, interactive = FALSE))
+  ctl <- p$layers[[1]]$data
+  expect_true(all(c("Passed", "Failed (p < 0.05)") %in% ctl$result))
+  # Long fleet names are wrapped, not cut
+  expect_true(all(grepl("\n", ctl$panel)))
+})
+
+test_that("plot_runs_test saves a caption and alt text that describe the tests", {
+  rda <- fit_rda(plot_runs_test(stockplotr::example_data, interactive = FALSE, make_rda = TRUE), "runs_test")
+  expect_tags_filled(rda, "runs.test.")
+  expect_match(rda$alt_text, "passed for 2 of 2 indices")
+  expect_match(rda$alt_text, "One residual is outside the control limits: Triennial in 2004.")
+  # The results are recorded in key_quantities.csv
+  resid <- index_residuals(stockplotr::example_data, interactive = FALSE)
+  kqs <- rda$kqs[grepl("^runs[.]test[.]", names(rda$kqs))]
+  expect_setequal(
+    names(kqs),
+    paste0("runs.test.", c("start.year", "end.year", "n.indices", "n.passed", "failed", "p.values", "n.outside"))
+  )
+  expect_equal(kqs[["runs.test.start.year"]], as.character(min(resid$year)))
+  expect_equal(kqs[["runs.test.n.indices"]], "2")
+  expect_equal(kqs[["runs.test.n.passed"]], "2")
+  expect_equal(kqs[["runs.test.failed"]], "none")
+  expect_equal(kqs[["runs.test.n.outside"]], "1")
+  expect_match(kqs[["runs.test.p.values"]], "^[^;]+: [0-9.]+; [^;]+: [0-9.]+$")
+
+  # save_all_plots() saves to a figures folder other than the working directory
+  figs <- tempfile("figs")
+  on.exit(unlink(figs, recursive = TRUE), add = TRUE)
+  rda <- fit_rda(
+    plot_runs_test(stockplotr::example_data, interactive = FALSE, make_rda = TRUE, figures_dir = figs),
+    "runs_test",
+    dir = figs
+  )
+  expect_tags_filled(rda, "runs.test.")
+})
diff --git a/tests/testthat/test-save_all_plots.R b/tests/testthat/test-save_all_plots.R
index 19eb6112..0e4ebe73 100644
--- a/tests/testthat/test-save_all_plots.R
+++ b/tests/testthat/test-save_all_plots.R
@@ -10,6 +10,7 @@ fig_base_temp_files <- c(
   "catch_comp_figure.rda",
   "recruitment_deviations_figure.rda",
   "recruitment_figure.rda",
+  "runs_test_figure.rda",
   "selectivity_figure.rda",
   "spawning_biomass_figure.rda",
   "spr_figure.rda",
diff --git a/tests/testthat/test-utils_fit_captions.R b/tests/testthat/test-utils_fit_captions.R
new file mode 100644
index 00000000..680882e0
--- /dev/null
+++ b/tests/testthat/test-utils_fit_captions.R
@@ -0,0 +1,50 @@
+test_that("caption text helpers format lists and panels", {
+  expect_equal(fit_list_text("A"), "A")
+  expect_equal(fit_list_text(c("A", "B")), "A and B")
+  expect_equal(fit_list_text(c("A", "B", "C", "A")), "A, B, and C")
+  expect_equal(fit_panel_text("a"), "A single panel for a.")
+  expect_equal(fit_panel_text(c("a", "b")), "A set of 2 panels, one for each fleet (a and b).")
+  expect_equal(fit_tag_value(c("costs $5", "a\\b")), "costs 5 ab")
+})
+
+test_that("fit_add_caption_rows adds missing rows to an existing captions file", {
+  dir <- tempfile("captions")
+  dir.create(dir)
+  old <- setwd(dir)
+  on.exit(
+    {
+      setwd(old)
+      unlink(dir, recursive = TRUE)
+    },
+    add = TRUE
+  )
+  template <- utils::read.csv(
+    system.file("resources", "captions_alt_text_template.csv", package = "stockplotr")
+  )
+  utils::write.csv(template[template$label != "runs_test", ], "captions_alt_text.csv", row.names = FALSE)
+  fit_add_caption_rows("runs_test")
+  caps <- utils::read.csv("captions_alt_text.csv")
+  expect_equal(sum(caps$label == "runs_test"), 1)
+  expect_equal(nrow(caps), nrow(template))
+  # Rows already present are not duplicated
+  fit_add_caption_rows("runs_test")
+  expect_equal(nrow(utils::read.csv("captions_alt_text.csv")), nrow(template))
+})
+
+test_that("fit_add_kq_rows adds missing rows to an existing key quantities file", {
+  dir <- tempfile("kqs")
+  dir.create(dir)
+  on.exit(unlink(dir, recursive = TRUE), add = TRUE)
+  template <- utils::read.csv(
+    system.file("resources", "key_quantity_template.csv", package = "stockplotr")
+  )
+  old <- template[!grepl("^runs[.]test[.]", template$key_quantity), ]
+  utils::write.csv(old, file.path(dir, "key_quantities.csv"), row.names = FALSE)
+  fit_add_kq_rows(c("runs.test.n.passed", "runs.test.failed"), dir = dir)
+  kqs <- utils::read.csv(file.path(dir, "key_quantities.csv"))
+  expect_equal(nrow(kqs), nrow(old) + 2)
+  expect_true(all(c("runs.test.n.passed", "runs.test.failed") %in% kqs$key_quantity))
+  # Rows already present are not duplicated
+  fit_add_kq_rows("runs.test.n.passed", dir = dir)
+  expect_equal(nrow(utils::read.csv(file.path(dir, "key_quantities.csv"))), nrow(old) + 2)
+})