diff --git a/NAMESPACE b/NAMESPACE
index a89a4c9c..9ffdce67 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -12,6 +12,15 @@ export(plot_abundance_at_age)
 export(plot_biomass)
 export(plot_biomass_at_age)
 export(plot_catch_comp)
+export(plot_comp)
+export(plot_comp_andre)
+export(plot_comp_bubble)
+export(plot_comp_francis)
+export(plot_comp_residual_dist)
+export(plot_comp_residual_pattern)
+export(plot_comp_residuals)
+export(plot_comp_sample_size)
+export(plot_comp_sex_ratio)
 export(plot_discard)
 export(plot_fishing_mortality)
 export(plot_index)
@@ -24,6 +33,7 @@ export(plot_selectivity)
 export(plot_spawning_biomass)
 export(plot_stock_recruitment)
 export(plot_timeseries)
+export(prepare_comp_data)
 export(process_data)
 export(process_table)
 export(reference_line)
diff --git a/R/data.R b/R/data.R
index aa90d1e6..f3c4bef4 100644
--- a/R/data.R
+++ b/R/data.R
@@ -41,3 +41,65 @@
 #' }
 #'
 "example_data"
+
+#' Big Skate length composition data
+#'
+#' Length compositions and fits from the 2019 Big Skate stock assessment for
+#' the U.S. West Coast (`BigSkate_2019` in the Stock Synthesis test models),
+#' the model used for the r4ss example plots at
+#' \url{https://r4ss.github.io/r4ss/BigSkate/}. The data are the `lendbase`
+#' element of `r4ss::SS_output()` with columns renamed to follow
+#' stockplotr conventions. See `data-raw/bigskate_comps.R`.
+#'
+#' Compositions are sex-specific, with female and male proportions summing to
+#' one within each year, fleet, and partition. `Fishery_current` has separate
+#' retained and discarded compositions.
+#'
+#' @format A tibble with 3626 rows and 11 variables:
+#' \describe{
+#'   \item{comp_type}{Type of composition ("length")}
+#'   \item{fleet}{Fleet name: "Fishery_current", "WCGBTS", or "Triennial"}
+#'   \item{year}{Year of the observation (`Yr` in r4ss)}
+#'   \item{month}{Month of the observation}
+#'   \item{sex}{"female" or "male"}
+#'   \item{part}{Partition: 0 = whole catch, 1 = discard, 2 = retained}
+#'   \item{bin}{Lower edge of the length bin (cm)}
+#'   \item{observed}{Observed proportion (`Obs` in r4ss)}
+#'   \item{expected}{Expected proportion (`Exp` in r4ss)}
+#'   \item{pearson}{Pearson residual reported by Stock Synthesis}
+#'   \item{n}{Adjusted input sample size (`Nsamp_adj` in r4ss)}
+#' }
+#' @source Taylor, I.G., Gertseva, V., Stephens, A., and Bizzarro, J. (2019). Status of
+#' Big Skate (Beringraja binoculata) off the U.S. Pacific Coast in 2019.
+#' Pacific Fishery Management Council, Portland, OR.
+#' @seealso [bigskate_caal], [prepare_comp_data()], [plot_comp()]
+"bigskate_lencomp"
+
+#' Big Skate conditional age-at-length data
+#'
+#' Conditional age-at-length data and fits from the 2019 Big Skate stock
+#' assessment (see [bigskate_lencomp]). The data are the `condbase` element
+#' of `r4ss::SS_output()` with columns renamed to follow stockplotr conventions.
+#' Each age composition is for one sex and one length bin, and sums to one.
+#' See `data-raw/bigskate_comps.R`.
+#'
+#' @format A tibble with 6560 rows and 12 variables:
+#' \describe{
+#'   \item{comp_type}{Type of composition ("age")}
+#'   \item{fleet}{Fleet name: "Fishery_current" or "WCGBTS"}
+#'   \item{year}{Year of the observation (`Yr` in r4ss)}
+#'   \item{month}{Month of the observation}
+#'   \item{sex}{"female" or "male"}
+#'   \item{part}{Partition: 0 = whole catch, 1 = discard, 2 = retained}
+#'   \item{length_bin}{Length bin the ages are conditioned on (`Lbin_lo` in r4ss)}
+#'   \item{bin}{Age}
+#'   \item{observed}{Observed proportion at age within the length bin}
+#'   \item{expected}{Expected proportion at age within the length bin}
+#'   \item{pearson}{Pearson residual reported by Stock Synthesis}
+#'   \item{n}{Adjusted input sample size (`Nsamp_adj` in r4ss)}
+#' }
+#' @source Taylor, I.G., Gertseva, V., Stephens, A., and Bizzarro, J. (2019). Status of
+#' Big Skate (Beringraja binoculata) off the U.S. Pacific Coast in 2019.
+#' Pacific Fishery Management Council, Portland, OR.
+#' @seealso [bigskate_lencomp], [prepare_comp_data()], [plot_comp_andre()]
+"bigskate_caal"
diff --git a/R/plot_comp.R b/R/plot_comp.R
new file mode 100644
index 00000000..851320c2
--- /dev/null
+++ b/R/plot_comp.R
@@ -0,0 +1,742 @@
+#' Plot composition distributions
+#'
+#' Plot age or length compositions as distributions (proportion by bin), one
+#' panel per year, optionally overlaying the model-expected values. Setting
+#' `aggregate = TRUE` instead sums the compositions across years so that each
+#' fleet is summarized in a single panel.
+#'
+#' @param dat A data frame of composition data in any format accepted by
+#' [prepare_comp_data()], such as FIMS output, an r4ss `lendbase`/`agedbase`
+#' table, or a tidy data frame like [bigskate_lencomp].
+#' @param show_expected Logical. TRUE/FALSE; overlay the expected values as
+#' lines. Ignored if `dat` contains no expected values.
+#'
+#' Default: TRUE
+#'
+#' @param aggregate Logical. TRUE/FALSE; aggregate compositions across years
+#' within each fleet. Years are weighted by their sample size when sample
+#' sizes are available.
+#'
+#' Default: FALSE
+#'
+#' @param show_n Logical. TRUE/FALSE; label each panel with its adjusted
+#' input sample size ("N adj.") and, when expected values are shown, its
+#' effective sample size ("N eff."), as in r4ss. Not used when
+#' `aggregate = TRUE`.
+#'
+#' Default: TRUE
+#'
+#' @param ncol Number. Number of facet columns. NULL lets `ggplot2` choose.
+#'
+#' Default: NULL
+#'
+#' @inheritParams prepare_comp_data
+#'
+#' @returns A `ggplot2` object.
+#'
+#' @details When compositions are sex-specific (both "female" and "male"
+#' values in the `sex` column), females are plotted above zero and males
+#' below zero. Bars show observed proportions and lines show expected
+#' proportions.
+#'
+#' If more than one fleet is plotted without aggregating, panels are
+#' created for every fleet and year combination; use `fleet` to focus on one
+#' fleet at a time.
+#'
+#' @section r4ss equivalent: `r4ss::SSplotComps()` subplot 1 (by year) and
+#' subplot 21 (aggregated across years).
+#'
+#' @seealso [prepare_comp_data()], [plot_comp_bubble()],
+#' [plot_comp_residuals()], [plot_comp_sex_ratio()], [plot_comp_francis()]
+#' @export
+#'
+#' @examples
+#' # Observations only
+#' plot_comp(bigskate_lencomp, fleet = "WCGBTS", show_expected = FALSE)
+#'
+#' # Observed and expected
+#' plot_comp(bigskate_lencomp, fleet = "Triennial")
+#'
+#' # All fleets aggregated across years
+#' plot_comp(bigskate_lencomp, aggregate = TRUE)
+plot_comp <- function(
+  dat,
+  show_expected = TRUE,
+  aggregate = FALSE,
+  comp_type = NULL,
+  fleet = NULL,
+  show_n = TRUE,
+  ncol = NULL
+) {
+  comps <- prepare_comp_data(dat, comp_type = comp_type, fleet = fleet)
+  bin_label <- comp_bin_label(comps)
+
+  if (aggregate) {
+    comps <- aggregate_comp_years(comps)
+  }
+  show_expected <- show_expected && any(!is.na(comps$expected))
+  two_sex <- all(c("female", "male") %in% comps$sex)
+
+  # Mirror males below zero for sex-specific comps
+  comps <- comps |>
+    dplyr::mutate(
+      direction = dplyr::if_else(two_sex & .data$sex == "male", -1, 1),
+      obs_plot = .data$observed * .data$direction,
+      exp_plot = .data$expected * .data$direction
+    )
+
+  bins <- sort(unique(comps$bin))
+  width <- if (length(bins) > 1) min(diff(bins)) else 1
+
+  plt <- ggplot2::ggplot(comps, ggplot2::aes(x = .data$bin))
+  if (length(unique(comps$sex)) > 1) {
+    plt <- plt +
+      ggplot2::geom_col(
+        ggplot2::aes(y = .data$obs_plot, fill = .data$sex),
+        width = width * 0.9,
+        alpha = 0.6,
+        position = if (two_sex) "identity" else "dodge"
+      )
+  } else {
+    plt <- plt +
+      ggplot2::geom_col(
+        ggplot2::aes(y = .data$obs_plot),
+        width = width * 0.9,
+        alpha = 0.6,
+        fill = "gray50"
+      )
+  }
+  if (show_expected) {
+    exp_dat <- dplyr::filter(comps, !is.na(.data$exp_plot))
+    if (length(unique(comps$sex)) > 1) {
+      plt <- plt +
+        ggplot2::geom_line(
+          data = exp_dat,
+          ggplot2::aes(y = .data$exp_plot, group = .data$sex, color = .data$sex),
+          linewidth = 0.7
+        )
+    } else {
+      plt <- plt +
+        ggplot2::geom_line(
+          data = exp_dat,
+          ggplot2::aes(y = .data$exp_plot, group = .data$sex),
+          color = "#0085CA",
+          linewidth = 0.7
+        )
+    }
+  }
+  if (two_sex) {
+    plt <- plt + ggplot2::geom_hline(yintercept = 0, linewidth = 0.3)
+  }
+
+  # Sample size labels, as in r4ss: adjusted input (N adj.) and effective
+  # (N eff.) sample sizes
+  if (show_n && !aggregate && any(!is.na(comps$n))) {
+    n_lab <- comps |>
+      dplyr::filter(!is.na(.data$n)) |>
+      dplyr::distinct(.data$comp_type, .data$fleet, .data$year, .data$n) |>
+      dplyr::mutate(label = paste0("N adj.=", signif(.data$n, 3)))
+    if (show_expected) {
+      n_lab <- n_lab |>
+        dplyr::left_join(
+          dplyr::select(comp_effn_table(comps), "comp_type", "fleet", "year", "effn"),
+          by = c("comp_type", "fleet", "year")
+        ) |>
+        dplyr::mutate(label = paste0(.data$label, "\nN eff.=", signif(.data$effn, 3)))
+    }
+    plt <- plt +
+      ggplot2::geom_text(
+        data = n_lab,
+        ggplot2::aes(x = Inf, y = Inf, label = .data$label),
+        hjust = 1.1,
+        vjust = 1.5,
+        size = 2.8,
+        inherit.aes = FALSE
+      )
+  }
+
+  ylab <- if (two_sex) "Proportion (females above, males below)" else "Proportion"
+  plt <- plt +
+    ggplot2::scale_y_continuous(labels = function(x) abs(x)) +
+    ggplot2::labs(
+      x = bin_label,
+      y = ylab,
+      fill = if (length(unique(comps$sex)) > 1) "Sex",
+      color = if (length(unique(comps$sex)) > 1 && show_expected) "Sex"
+    ) +
+    theme_noaa()
+
+  # Facets
+  facet_vars <- if (aggregate) "fleet" else c("fleet", "year")
+  if (!aggregate && length(unique(comps$fleet)) == 1) facet_vars <- "year"
+  if (length(unique(comps$comp_type)) > 1) facet_vars <- c("comp_type", facet_vars)
+  if (!aggregate && length(unique(comps$fleet)) > 1) {
+    cli::cli_alert_info(
+      "Plotting {length(unique(comps$fleet))} fleets by year; use {.arg fleet} to plot one fleet at a time."
+    )
+  }
+  plt <- plt + comp_facet(facet_vars, ncol = ncol)
+  if (!aggregate) {
+    plt <- plt + ggplot2::theme(
+      strip.text = ggplot2::element_text(size = 9),
+      axis.text = ggplot2::element_text(size = 8)
+    )
+  }
+  plt
+}
+
+#------------------------------------------------------------------------------
+
+#' Plot composition bubble plots
+#'
+#' Bubble plot of compositions by year (x-axis) and bin (y-axis). Bubbles can
+#' show the observed (or expected) proportions, which makes strong and weak
+#' cohorts easy to follow through time, or the Pearson residuals, which show
+#' blocks or stripes of misfit among years or bins.
+#'
+#' @inheritParams plot_comp
+#' @param value String. Value to plot. Options: "observed", "expected", or
+#' "pearson".
+#'
+#' Default: "observed"
+#'
+#' @param max_size Number. Size of the largest bubble.
+#'
+#' Default: 6
+#'
+#' @param conditional Logical. TRUE/FALSE; plot conditional age-at-length
+#' data (see [prepare_comp_data()]), with age on the y-axis, length bin on
+#' the x-axis, and one panel per year.
+#'
+#' Default: FALSE
+#'
+#' @returns A `ggplot2` object.
+#'
+#' @details Panels are created for each fleet, and for each sex when
+#' compositions are sex-specific. For Pearson residuals, filled bubbles are
+#' positive residuals (observed > expected) and open bubbles are negative
+#' residuals (observed < expected), with bubble area proportional to the
+#' absolute value of the residual. The maximum absolute residual for each
+#' panel is included in the panel label.
+#'
+#' @section r4ss equivalent: `r4ss::SSplotComps()` subplot 2 (bubble plot of
+#' observations or residuals) and subplot 3 (conditional age-at-length).
+#'
+#' @seealso [prepare_comp_data()], [plot_comp_residuals()], [plot_comp()]
+#' @export
+#'
+#' @examples
+#' # Observed proportions
+#' plot_comp_bubble(bigskate_lencomp)
+#'
+#' # Pearson residuals
+#' plot_comp_bubble(bigskate_lencomp, value = "pearson")
+#'
+#' # Residuals for conditional age-at-length data
+#' plot_comp_bubble(
+#'   dplyr::filter(bigskate_caal, year %in% 2016:2018),
+#'   value = "pearson",
+#'   fleet = "WCGBTS",
+#'   conditional = TRUE
+#' )
+plot_comp_bubble <- function(
+  dat,
+  value = c("observed", "pearson", "expected"),
+  comp_type = NULL,
+  fleet = NULL,
+  max_size = 6,
+  ncol = NULL,
+  conditional = FALSE
+) {
+  value <- match.arg(value)
+  comps <- prepare_comp_data(
+    dat,
+    comp_type = if (conditional) "age" else comp_type,
+    fleet = fleet,
+    conditional = conditional
+  )
+  bin_label <- comp_bin_label(comps)
+  # Conditional age-at-length: age by length bin, one panel per year
+  x_var <- if (conditional) "length_bin" else "year"
+  comps$value <- comps[[value]]
+  comps <- dplyr::filter(comps, !is.na(.data$value))
+  if (nrow(comps) == 0) {
+    if (value == "pearson") {
+      cli::cli_abort(c(
+        "No Pearson residuals available.",
+        "i" = "Supply a {.field pearson} column, or {.field expected} values and sample sizes ({.field n}) so residuals can be calculated."
+      ))
+    }
+    cli::cli_abort("No {value} values found in {.arg dat}.")
+  }
+
+  facet_vars <- c(
+    if (length(unique(comps$comp_type)) > 1) "comp_type",
+    if (!conditional || length(unique(comps$fleet)) > 1) "fleet",
+    if (length(unique(comps$sex)) > 1) "sex",
+    if (conditional) "year"
+  )
+  if (value == "pearson") {
+    # Add the max absolute residual to a single panel label, ordered by
+    # comp type, fleet, and sex
+    comps <- comps |>
+      dplyr::arrange(.data$comp_type, .data$fleet, .data$sex, .data$year) |>
+      dplyr::group_by(dplyr::across(dplyr::all_of(facet_vars))) |>
+      dplyr::mutate(
+        panel = paste0(
+          paste(unlist(dplyr::cur_group()), collapse = ", "),
+          "\n(max = ", format(round(max(abs(.data$value)), 2), nsmall = 2), ")"
+        )
+      ) |>
+      dplyr::ungroup() |>
+      dplyr::mutate(panel = factor(.data$panel, levels = unique(.data$panel))) |>
+      dplyr::mutate(
+        sign = factor(
+          dplyr::if_else(.data$value >= 0, "Observed > expected", "Observed < expected"),
+          levels = c("Observed > expected", "Observed < expected")
+        )
+      )
+    plt <- ggplot2::ggplot(
+      comps,
+      ggplot2::aes(x = .data[[x_var]], y = .data$bin, size = abs(.data$value), fill = .data$sign)
+    ) +
+      ggplot2::geom_point(shape = 21, color = "black", alpha = 0.8) +
+      ggplot2::scale_fill_manual(
+        values = c("Observed > expected" = "gray20", "Observed < expected" = "white"),
+        drop = FALSE
+      ) +
+      ggplot2::labs(size = "|Pearson residual|", fill = NULL) +
+      ggplot2::guides(fill = ggplot2::guide_legend(override.aes = list(size = 4)))
+    facet_vars <- "panel"
+  } else {
+    plt <- ggplot2::ggplot(
+      comps,
+      ggplot2::aes(x = .data[[x_var]], y = .data$bin, size = .data$value)
+    ) +
+      ggplot2::geom_point(shape = 21, color = "black", fill = "gray40", alpha = 0.5) +
+      ggplot2::labs(size = cap_first_letter(paste(value, "proportion")))
+  }
+
+  plt +
+    ggplot2::scale_size_area(max_size = max_size) +
+    ggplot2::scale_x_continuous(
+      breaks = if (conditional) ggplot2::waiver() else axis_breaks(comps$year),
+      guide = ggplot2::guide_axis(minor.ticks = TRUE)
+    ) +
+    ggplot2::labs(x = if (conditional) "Length bin" else "Year", y = bin_label) +
+    theme_noaa() +
+    comp_facet(facet_vars, ncol = ncol)
+}
+
+#------------------------------------------------------------------------------
+
+#' Plot composition Pearson residuals
+#'
+#' A convenience wrapper around [plot_comp_bubble()] with
+#' `value = "pearson"`. Residuals are taken from `dat` when supplied or
+#' calculated by [prepare_comp_data()] assuming a multinomial distribution.
+#'
+#' @inheritParams plot_comp_bubble
+#' @param ... Additional arguments passed to [plot_comp_bubble()].
+#'
+#' @returns A `ggplot2` object.
+#' @section r4ss equivalent: `r4ss::SSplotComps()` subplot 2 for residuals by
+#' fleet and subplot 24 (residuals compared across fleets).
+#'
+#' @seealso [plot_comp_bubble()], [prepare_comp_data()]
+#' @export
+#'
+#' @examples
+#' plot_comp_residuals(bigskate_lencomp)
+#'
+#' # Compare residuals across fleets for females
+#' plot_comp_residuals(
+#'   dplyr::filter(bigskate_lencomp, sex == "female"),
+#'   ncol = 1
+#' )
+plot_comp_residuals <- function(dat, ...) {
+  plot_comp_bubble(dat, value = "pearson", ...)
+}
+
+#------------------------------------------------------------------------------
+
+#' Plot sex ratios implied by composition data
+#'
+#' Plot the proportion female, female / (female + male), by bin for the
+#' observed and expected compositions. Trends in the lack of fit to sex ratios
+#' can be easier to see here than in residuals, especially for bins with
+#' small proportions.
+#'
+#' @inheritParams plot_comp
+#' @param interval Number. Coverage of the interval around the observed sex
+#' ratio, calculated as a Jeffreys interval using the sample size in each bin
+#' (the sample size multiplied by the proportion of the composition in the
+#' bin). Set to NULL to omit intervals. Intervals are only shown when sample
+#' sizes are available.
+#'
+#' Default: 0.75
+#'
+#' @returns A `ggplot2` object.
+#'
+#' @details Requires sex-specific compositions with both "female" and
+#' "male" values in the `sex` column of `dat`. Points show observed ratios
+#' and lines show expected ratios. As in r4ss, the smallest observed
+#' proportion in each composition (such as the small constant Stock Synthesis
+#' adds to every bin) is subtracted from the observed and expected values
+#' first, so bins with no fish observed are omitted. When `aggregate = TRUE`, female and male
+#' compositions are summed across years (weighted by sample size when
+#' available) before calculating ratios.
+#'
+#' @section r4ss equivalent: `r4ss::SSplotSexRatio()` with
+#' `sexratio.option = 2` (female:total).
+#'
+#' @seealso [prepare_comp_data()], [plot_comp()]
+#' @export
+#'
+#' @examples
+#' plot_comp_sex_ratio(bigskate_lencomp, fleet = "WCGBTS")
+#' plot_comp_sex_ratio(bigskate_lencomp, aggregate = TRUE)
+plot_comp_sex_ratio <- function(
+  dat,
+  aggregate = FALSE,
+  interval = 0.75,
+  comp_type = NULL,
+  fleet = NULL,
+  ncol = NULL
+) {
+  comps <- prepare_comp_data(dat, comp_type = comp_type, fleet = fleet)
+  bin_label <- comp_bin_label(comps)
+  if (!all(c("female", "male") %in% comps$sex)) {
+    cli::cli_abort("Sex ratios require compositions with both {.val female} and {.val male} values in the {.field sex} column.")
+  }
+
+  # Females and males side by side, weighted by sample size when available.
+  # As in r4ss, the smallest observation in each composition (e.g., the
+  # constant Stock Synthesis adds to every bin) is subtracted first, so that
+  # bins without fish do not appear as 50:50 sex ratios.
+  wide <- comps |>
+    dplyr::filter(.data$sex %in% c("female", "male")) |>
+    dplyr::group_by(.data$comp_type, .data$fleet, .data$year) |>
+    dplyr::mutate(
+      min_obs = min(.data$observed, na.rm = TRUE),
+      observed = pmax(.data$observed - .data$min_obs, 0),
+      expected = pmax(.data$expected - .data$min_obs, 0)
+    ) |>
+    dplyr::ungroup() |>
+    dplyr::mutate(w = dplyr::if_else(is.na(.data$n), 1, .data$n)) |>
+    dplyr::select("comp_type", "fleet", "year", "sex", "bin", "observed", "expected", "n", "w") |>
+    tidyr::pivot_wider(names_from = "sex", values_from = c("observed", "expected"))
+  group_vars <- if (aggregate) c("comp_type", "fleet", "bin") else c("comp_type", "fleet", "year", "bin")
+
+  ratios <- wide |>
+    dplyr::group_by(dplyr::across(dplyr::all_of(group_vars))) |>
+    dplyr::summarise(
+      f_obs = sum(.data$observed_female * .data$w, na.rm = TRUE),
+      m_obs = sum(.data$observed_male * .data$w, na.rm = TRUE),
+      f_exp = sum(.data$expected_female * .data$w, na.rm = TRUE),
+      m_exp = sum(.data$expected_male * .data$w, na.rm = TRUE),
+      has_n = all(!is.na(.data$n)),
+      has_exp = any(!is.na(.data$expected_female)),
+      .groups = "drop"
+    ) |>
+    dplyr::mutate(
+      ratio_obs = .data$f_obs / (.data$f_obs + .data$m_obs),
+      ratio_exp = dplyr::if_else(.data$has_exp, .data$f_exp / (.data$f_exp + .data$m_exp), NA_real_),
+      # weighted sums equal the number of fish in the bin when w is the sample size
+      # (at least 1, as in r4ss)
+      n_bin = dplyr::if_else(
+        .data$has_n & .data$f_obs + .data$m_obs > 0,
+        pmax(.data$f_obs + .data$m_obs, 1),
+        NA_real_
+      )
+    )
+  ratios$ratio_obs[!is.finite(ratios$ratio_obs)] <- NA_real_
+  ratios$ratio_exp[!is.finite(ratios$ratio_exp)] <- NA_real_
+
+  show_interval <- !is.null(interval) && any(!is.na(ratios$n_bin))
+  if (show_interval) {
+    ci <- jeffreys_interval(ratios$ratio_obs, ratios$n_bin, level = interval)
+    ratios$lower <- ci$lower
+    ratios$upper <- ci$upper
+  }
+
+  plt <- ggplot2::ggplot(ratios, ggplot2::aes(x = .data$bin)) +
+    ggplot2::geom_hline(yintercept = 0.5, linetype = "dashed", color = "gray60")
+
+  if (show_interval) {
+    plt <- plt +
+      ggplot2::geom_linerange(
+        data = dplyr::filter(ratios, !is.na(.data$ratio_obs), !is.na(.data$lower)),
+        ggplot2::aes(ymin = .data$lower, ymax = .data$upper),
+        color = "gray50"
+      )
+  }
+  plt <- plt +
+    ggplot2::geom_point(
+      data = dplyr::filter(ratios, !is.na(.data$ratio_obs)),
+      ggplot2::aes(y = .data$ratio_obs, shape = "Observed"),
+      size = 1.5
+    )
+  if (any(!is.na(ratios$ratio_exp))) {
+    plt <- plt +
+      ggplot2::geom_line(
+        data = dplyr::filter(ratios, !is.na(.data$ratio_exp)),
+        ggplot2::aes(y = .data$ratio_exp, linetype = "Expected"),
+        color = "#0085CA",
+        linewidth = 0.8
+      )
+  }
+
+  facet_vars <- if (aggregate) "fleet" else c("fleet", "year")
+  if (!aggregate && length(unique(ratios$fleet)) == 1) facet_vars <- "year"
+  if (length(unique(ratios$comp_type)) > 1) facet_vars <- c("comp_type", facet_vars)
+
+  plt <- plt +
+    ggplot2::scale_y_continuous(limits = c(0, 1), breaks = c(0, 0.5, 1)) +
+    ggplot2::labs(
+      x = bin_label,
+      y = "Fraction female",
+      shape = NULL,
+      linetype = NULL,
+      caption = if (!is.null(interval) && any(!is.na(ratios$n_bin))) {
+        paste0("Vertical lines: ", round(interval * 100), "% Jeffreys intervals for observed ratios")
+      }
+    ) +
+    theme_noaa() +
+    comp_facet(facet_vars, ncol = ncol)
+  if (!aggregate) {
+    plt <- plt + ggplot2::theme(
+      strip.text = ggplot2::element_text(size = 9),
+      axis.text = ggplot2::element_text(size = 8)
+    )
+  }
+  plt
+}
+
+#------------------------------------------------------------------------------
+
+#' Plot mean age or length with Francis (TA1.8) data weights
+#'
+#' Time series of the mean age or length in the observed and expected
+#' compositions. Developed as a diagnostic for the Francis (2011) data
+#' weighting method (TA1.8), the plot is also a useful summary of trends in
+#' compositions and lack of fit among years.
+#'
+#' @inheritParams plot_comp
+#' @param interval Number. Coverage of the interval around the observed mean.
+#'
+#' Default: 0.95
+#'
+#' @param show_francis Logical. TRUE/FALSE; add the Francis TA1.8 weight for
+#' each fleet to the panel label and show thinner, capped intervals
+#' representing sample sizes after applying that weight.
+#'
+#' Default: TRUE
+#'
+#' @returns A `ggplot2` object.
+#'
+#' @details For each year and fleet (with sexes combined), the observed mean
+#' is \eqn{\bar{x}_o = \sum_b x_b o_b} and the expected mean is
+#' \eqn{\bar{x}_e = \sum_b x_b e_b}, where \eqn{x_b} is the bin and \eqn{o_b}
+#' and \eqn{e_b} are the observed and expected proportions. Assuming a
+#' multinomial distribution, the standard error of the mean is
+#' \eqn{\sqrt{v_e / N}}, where \eqn{v_e = \sum_b x_b^2 e_b - \bar{x}_e^2} is
+#' the variance of the expected composition and \eqn{N} is the sample size.
+#' Intervals are drawn around the observed means.
+#'
+#' The Francis weight is \eqn{1 / \mathrm{var}(r)} where
+#' \eqn{r = (\bar{x}_o - \bar{x}_e) / SE} are the standardized residuals of
+#' the mean across years. As in Stock Synthesis and
+#' `r4ss::SSMethod.TA1.8()`, all partitions of a fleet (e.g., retained and
+#' discarded) share one weight, and its 95% interval is calculated by
+#' bootstrapping the standardized residuals (1000 resamples), so use
+#' [set.seed()] for reproducible intervals. The weight assumes a multinomial likelihood and has not been
+#' adapted for Dirichlet-multinomial or other likelihoods.
+#'
+#' Without expected values, only the observed means are plotted. Without
+#' sample sizes, no intervals or weights are shown.
+#'
+#' @references Francis, R.I.C.C. (2011). Data weighting in statistical
+#' fisheries stock assessment models. Canadian Journal of Fisheries and
+#' Aquatic Sciences 68: 1124-1138.
+#'
+#' @section r4ss equivalent: `r4ss::SSplotComps()` subplots 8 and 9, which
+#' use `r4ss::SSMethod.TA1.8()`.
+#'
+#' @seealso [prepare_comp_data()], [plot_comp()]
+#' @export
+#'
+#' @examples
+#' plot_comp_francis(bigskate_lencomp)
+plot_comp_francis <- function(
+  dat,
+  interval = 0.95,
+  show_francis = TRUE,
+  comp_type = NULL,
+  fleet = NULL,
+  ncol = NULL
+) {
+  comps <- prepare_comp_data(dat, comp_type = comp_type, fleet = fleet)
+  bin_label <- comp_bin_label(comps)
+  means <- comp_mean_table(comps)
+  z <- stats::qnorm(1 - (1 - interval) / 2)
+  has_exp <- any(!is.na(means$mean_exp))
+  has_se <- any(!is.na(means$se))
+
+  # Francis weights by fleet. As in Stock Synthesis and
+  # r4ss::SSMethod.TA1.8(), partitions of a fleet share one weight.
+  weights <- NULL
+  if (show_francis && has_exp && has_se) {
+    means$base_fleet <- base_fleet(means$fleet)
+    weights <- means |>
+      dplyr::filter(!is.na(.data$std_resid)) |>
+      dplyr::group_by(.data$comp_type, .data$base_fleet) |>
+      dplyr::filter(dplyr::n() > 1) |>
+      dplyr::summarise(
+        n_parts = length(unique(.data$fleet)),
+        weight = 1 / stats::var(.data$std_resid),
+        weight_ci = list(francis_boot_ci(.data$std_resid)),
+        .groups = "drop"
+      ) |>
+      dplyr::mutate(
+        weight_label = paste0(
+          "TA1.8 weight", dplyr::if_else(.data$n_parts > 1, " (all partitions)", ""),
+          " = ", signif(.data$weight, 3), " (",
+          signif(purrr::map_dbl(.data$weight_ci, 1), 3), "-",
+          signif(purrr::map_dbl(.data$weight_ci, 2), 3), ")"
+        )
+      )
+    means <- means |>
+      dplyr::left_join(
+        dplyr::select(weights, "comp_type", "base_fleet", "weight", "weight_label"),
+        by = c("comp_type", "base_fleet")
+      ) |>
+      dplyr::mutate(
+        fleet = dplyr::if_else(
+          is.na(.data$weight_label),
+          .data$fleet,
+          paste0(.data$fleet, "\n", .data$weight_label)
+        ),
+        se_adj = .data$se / sqrt(.data$weight)
+      )
+  }
+
+  plt <- ggplot2::ggplot(means, ggplot2::aes(x = .data$year))
+  if (has_se) {
+    plt <- plt +
+      ggplot2::geom_linerange(
+        ggplot2::aes(
+          ymin = .data$mean_obs - z * .data$se,
+          ymax = .data$mean_obs + z * .data$se
+        ),
+        color = "gray40",
+        linewidth = 0.8
+      )
+    if (!is.null(weights) && nrow(weights) > 0) {
+      plt <- plt +
+        ggplot2::geom_errorbar(
+          ggplot2::aes(
+            ymin = .data$mean_obs - z * .data$se_adj,
+            ymax = .data$mean_obs + z * .data$se_adj
+          ),
+          color = "gray20",
+          linewidth = 0.3,
+          width = 0.4
+        )
+    }
+  }
+  plt <- plt +
+    ggplot2::geom_point(ggplot2::aes(y = .data$mean_obs, shape = "Observed"), size = 2)
+  if (has_exp) {
+    plt <- plt +
+      ggplot2::geom_line(
+        data = dplyr::filter(means, !is.na(.data$mean_exp)),
+        ggplot2::aes(y = .data$mean_exp, linetype = "Expected"),
+        color = "#0085CA",
+        linewidth = 0.8
+      )
+  }
+
+  caption <- NULL
+  if (has_se) {
+    caption <- paste0(round(interval * 100), "% intervals based on current sample sizes")
+    if (!is.null(weights) && nrow(weights) > 0) {
+      caption <- paste0(
+        caption,
+        "; thinner intervals (with capped ends) show result of further adjusting sample sizes (Francis TA1.8)"
+      )
+    }
+  }
+
+  facet_vars <- c(if (length(unique(means$comp_type)) > 1) "comp_type", "fleet")
+  plt +
+    ggplot2::scale_x_continuous(
+      breaks = axis_breaks(means$year),
+      guide = ggplot2::guide_axis(minor.ticks = TRUE)
+    ) +
+    ggplot2::labs(
+      x = "Year",
+      y = paste("Mean", tolower(bin_label)),
+      shape = NULL,
+      linetype = NULL,
+      caption = caption
+    ) +
+    theme_noaa() +
+    comp_facet(facet_vars, ncol = ncol, scales = "free_y")
+}
+
+# Bootstrap 95% interval for the Francis weight, following
+# r4ss::SSMethod.TA1.8() (1000 resamples of the standardized residuals)
+francis_boot_ci <- function(std_resid, n_boot = 1000) {
+  boot <- matrix(
+    sample(std_resid, n_boot * length(std_resid), replace = TRUE),
+    nrow = length(std_resid)
+  )
+  w <- apply(boot, 2, function(x) 1 / stats::var(x))
+  as.vector(stats::quantile(w, c(0.025, 0.975), na.rm = TRUE))
+}
+
+# Mean bin of observed and expected compositions by year and fleet. As in
+# r4ss::SSMethod.TA1.8(), sexes and partitions of a fleet in the same year
+# are combined into one composition, with the mean sample size.
+comp_mean_table <- function(comps) {
+  comps <- comps |>
+    dplyr::mutate(
+      part_name = fleet_partition(.data$fleet),
+      fleet = base_fleet(.data$fleet)
+    )
+  comps |>
+    dplyr::group_by(.data$comp_type, .data$fleet) |>
+    dplyr::mutate(
+      fleet = if (length(unique(stats::na.omit(.data$part_name))) > 1) {
+        paste0(.data$fleet, " (", paste(sort(unique(stats::na.omit(.data$part_name))), collapse = " & "), ")")
+      } else {
+        .data$fleet
+      }
+    ) |>
+    dplyr::group_by(.data$comp_type, .data$fleet, .data$year, .data$bin) |>
+    dplyr::summarise(
+      observed = sum(.data$observed, na.rm = TRUE),
+      expected = if (all(is.na(.data$expected))) NA_real_ else sum(.data$expected, na.rm = TRUE),
+      n = mean(.data$n),
+      .groups = "drop"
+    ) |>
+    dplyr::group_by(.data$comp_type, .data$fleet, .data$year) |>
+    dplyr::summarise(
+      mean_obs = sum(.data$bin * .data$observed) / sum(.data$observed),
+      mean_exp = sum(.data$bin * .data$expected) / sum(.data$expected),
+      # Use the expected distribution for the variance when available
+      var_bin = if (all(is.na(.data$expected))) {
+        sum(.data$bin^2 * .data$observed) / sum(.data$observed) - .data$mean_obs[1]^2
+      } else {
+        sum(.data$bin^2 * .data$expected) / sum(.data$expected) - .data$mean_exp[1]^2
+      },
+      n = mean(.data$n),
+      .groups = "drop"
+    ) |>
+    dplyr::filter(is.finite(.data$mean_obs)) |>
+    dplyr::mutate(
+      se = sqrt(.data$var_bin / .data$n),
+      std_resid = (.data$mean_obs - .data$mean_exp) / .data$se
+    )
+}
diff --git a/R/plot_comp_diagnostics.R b/R/plot_comp_diagnostics.R
new file mode 100644
index 00000000..6d96e008
--- /dev/null
+++ b/R/plot_comp_diagnostics.R
@@ -0,0 +1,599 @@
+#' Plot input versus effective sample size of compositions
+#'
+#' Compare the input (or adjusted input) sample size of each composition with
+#' its effective sample size, the sample size implied by how well the model
+#' fits the composition. This is the diagnostic behind the McAllister and
+#' Ianelli (1997) data weighting method.
+#'
+#' @inheritParams plot_comp
+#' @param dat Composition data in any format accepted by
+#' [prepare_comp_data()], or [convert_output()] output from Stock Synthesis
+#' that contains the `FIT_LEN_COMPS` or `FIT_AGE_COMPS` tables (as in
+#' `stockplotr::example_data`), from which the adjusted input (`nsamp_adj`)
+#' and effective (`effn`) sample sizes are taken.
+#' @param show_years Logical. TRUE/FALSE; label each point with its year.
+#'
+#' Default: FALSE
+#'
+#' @returns A `ggplot2` object.
+#'
+#' @details For each composition (year and fleet), the effective sample size
+#' is
+#' \deqn{N_{eff} = \sum_b e_b (1 - e_b) / \sum_b (o_b - e_b)^2}
+#' where \eqn{o_b} and \eqn{e_b} are the observed and expected proportions in
+#' bin \eqn{b} (McAllister and Ianelli 1997), the same calculation used by
+#' Stock Synthesis.
+#'
+#' Each panel shows a 1:1 line, a loess smoother (dashed red, when there are
+#' more than six distinct input sample sizes), the arithmetic mean of the
+#' input sample sizes (vertical dashed line), and the harmonic mean of the
+#' effective sample sizes (horizontal dashed line). The ratio of the harmonic
+#' mean effective sample size to the mean input sample size, the
+#' McAllister-Ianelli weight, is included in each panel label. Ratios below
+#' one suggest that input sample sizes are too large relative to the fit. As
+#' in Stock Synthesis, the weight is calculated across all partitions of a
+#' fleet; multiplying it by the current variance adjustment gives the
+#' `Recommend_var_adj` reported by Stock Synthesis.
+#'
+#' Effective sample sizes from this formula are not appropriate for
+#' compositions fit with a Dirichlet-multinomial likelihood.
+#'
+#' @references McAllister, M.K., and Ianelli, J.N. (1997). Bayesian stock
+#' assessment using catch-age data and the sampling-importance resampling
+#' algorithm. Canadian Journal of Fisheries and Aquatic Sciences 54: 284-300.
+#'
+#' @section r4ss equivalent: `r4ss::SSplotComps()` subplot 7 ("N-EffN
+#' comparison").
+#'
+#' @seealso [prepare_comp_data()], [plot_comp_francis()]
+#' @export
+#'
+#' @examples
+#' plot_comp_sample_size(bigskate_lencomp)
+#'
+#' # Stock Synthesis output converted with convert_output()
+#' plot_comp_sample_size(stockplotr::example_data, comp_type = "length")
+plot_comp_sample_size <- function(
+  dat,
+  comp_type = NULL,
+  fleet = NULL,
+  show_years = FALSE,
+  ncol = NULL
+) {
+  if (is_ss3_fit_summary(dat)) {
+    if (!is.null(comp_type)) comp_type <- match.arg(comp_type, c("age", "length"))
+    sizes <- comp_sample_size_from_ss3(dat, comp_type = comp_type)
+    if (!is.null(fleet)) {
+      sizes <- sizes[sizes$fleet %in% as.character(fleet), , drop = FALSE]
+    }
+    if (nrow(sizes) == 0) {
+      cli::cli_abort("No sample sizes remain after filtering by {.arg comp_type} and {.arg fleet}.")
+    }
+  } else {
+    comps <- prepare_comp_data(dat, comp_type = comp_type, fleet = fleet)
+    sizes <- comp_effn_table(comps)
+  }
+  sizes <- dplyr::filter(sizes, is.finite(.data$n), is.finite(.data$effn))
+  if (nrow(sizes) == 0) {
+    cli::cli_abort(c(
+      "Effective sample sizes could not be calculated.",
+      "i" = "Supply {.field expected} values and sample sizes ({.field n})."
+    ))
+  }
+
+  facet_vars <- c(if (length(unique(sizes$comp_type)) > 1) "comp_type", "fleet")
+  # Harmonic mean effective N and mean input N by panel, and the
+  # McAllister-Ianelli weight by fleet. As in Stock Synthesis, partitions of
+  # a fleet share one weight.
+  sizes$base_fleet <- base_fleet(sizes$fleet)
+  mi <- sizes |>
+    dplyr::group_by(.data$comp_type, .data$base_fleet) |>
+    dplyr::summarise(
+      n_parts = length(unique(.data$fleet)),
+      mi_weight = (1 / mean(1 / .data$effn)) / mean(.data$n),
+      .groups = "drop"
+    )
+  summ <- sizes |>
+    dplyr::group_by(dplyr::across(dplyr::all_of(unique(c("comp_type", facet_vars, "base_fleet"))))) |>
+    dplyr::summarise(
+      mean_n = mean(.data$n),
+      hmean_effn = 1 / mean(1 / .data$effn),
+      .groups = "drop"
+    ) |>
+    dplyr::left_join(mi, by = c("comp_type", "base_fleet")) |>
+    dplyr::mutate(
+      panel = paste0(
+        .data$fleet, "\nMcAllister-Ianelli weight",
+        dplyr::if_else(.data$n_parts > 1, " (all partitions)", ""),
+        " = ", signif(.data$mi_weight, 3)
+      )
+    )
+  sizes <- dplyr::left_join(sizes, summ, by = unique(c("comp_type", facet_vars, "base_fleet")))
+
+  # Loess smoother by panel when there are enough distinct sample sizes
+  smooth <- lapply(split(sizes, sizes$panel), function(d) {
+    if (length(unique(d$n)) <= 6 || diff(range(d$n)) <= 2) {
+      return(NULL)
+    }
+    fit <- suppressWarnings(stats::loess(effn ~ n, data = d, degree = 1))
+    d <- d[order(d$n), ]
+    d$fit <- stats::predict(fit, newdata = data.frame(n = d$n))
+    d
+  })
+  smooth <- dplyr::bind_rows(smooth)
+
+  plt <- ggplot2::ggplot(sizes, ggplot2::aes(x = .data$n, y = .data$effn)) +
+    ggplot2::geom_abline(slope = 1, intercept = 0, color = "black") +
+    ggplot2::geom_vline(
+      data = summ,
+      ggplot2::aes(xintercept = .data$mean_n),
+      linetype = "22",
+      color = "#0085CA"
+    ) +
+    ggplot2::geom_hline(
+      data = summ,
+      ggplot2::aes(yintercept = .data$hmean_effn),
+      linetype = "22",
+      color = "#0085CA"
+    ) +
+    ggplot2::geom_point(color = "gray30", size = 2)
+  if (nrow(smooth) > 0) {
+    plt <- plt +
+      ggplot2::geom_line(
+        data = smooth,
+        ggplot2::aes(x = .data$n, y = .data$fit),
+        color = "#D55E00",
+        linetype = "dashed"
+      )
+  }
+  if (show_years) {
+    plt <- plt +
+      ggplot2::geom_text(
+        ggplot2::aes(label = .data$year),
+        hjust = -0.2,
+        size = 2.8
+      )
+  }
+  plt +
+    ggplot2::expand_limits(x = 0, y = 0) +
+    ggplot2::labs(
+      x = "Adjusted input sample size (N adj.)",
+      y = "Effective sample size (N eff.)",
+      caption = "Solid line: 1:1. Dashed blue lines: arithmetic mean N adj. and harmonic mean N eff."
+    ) +
+    theme_noaa() +
+    comp_facet(c(if ("comp_type" %in% facet_vars) "comp_type", "panel"), ncol = ncol, scales = "free")
+}
+
+#------------------------------------------------------------------------------
+
+#' Plot the distribution of composition residuals
+#'
+#' Quantile-quantile plot or histogram of Pearson residuals for each fleet,
+#' to check for heavy tails, skewness, and outliers.
+#'
+#' @inheritParams plot_comp
+#' @param type String. Type of plot. Options: "qq" (normal quantile-quantile
+#' plot) or "histogram" (with a standard normal density overlaid).
+#'
+#' Default: "qq"
+#'
+#' @returns A `ggplot2` object.
+#'
+#' @details The standard deviation of the normalized residuals (SDNR) for
+#' each fleet is included in the panel label. Values well above one suggest
+#' that the compositions are given more weight (larger sample sizes) than the
+#' fit supports. Pearson residuals for compositions are not expected to be
+#' exactly normally distributed, especially for bins with small expected
+#' proportions or small sample sizes, so treat departures in the tails with
+#' some caution.
+#'
+#' @seealso [plot_comp_residuals()], [plot_comp_residual_pattern()],
+#' [prepare_comp_data()]
+#' @export
+#'
+#' @examples
+#' plot_comp_residual_dist(bigskate_lencomp)
+#' plot_comp_residual_dist(bigskate_lencomp, type = "histogram")
+plot_comp_residual_dist <- function(
+  dat,
+  type = c("qq", "histogram"),
+  comp_type = NULL,
+  fleet = NULL,
+  ncol = NULL
+) {
+  type <- match.arg(type)
+  comps <- comp_residuals_or_abort(dat, comp_type = comp_type, fleet = fleet)
+  facet_vars <- c(if (length(unique(comps$comp_type)) > 1) "comp_type", "fleet")
+  comps <- comps |>
+    dplyr::group_by(dplyr::across(dplyr::all_of(facet_vars))) |>
+    dplyr::mutate(
+      panel = paste0(.data$fleet, " (SDNR = ", format(round(stats::sd(.data$pearson), 2), nsmall = 2), ")")
+    ) |>
+    dplyr::ungroup()
+
+  if (type == "qq") {
+    plt <- ggplot2::ggplot(comps, ggplot2::aes(sample = .data$pearson)) +
+      ggplot2::geom_qq_line(color = "#0085CA", linewidth = 0.8) +
+      ggplot2::geom_qq(size = 1, alpha = 0.5) +
+      ggplot2::labs(x = "Theoretical quantiles (standard normal)", y = "Pearson residual")
+  } else {
+    plt <- ggplot2::ggplot(comps, ggplot2::aes(x = .data$pearson)) +
+      ggplot2::geom_histogram(
+        ggplot2::aes(y = ggplot2::after_stat(.data$density)),
+        bins = 30,
+        fill = "gray60",
+        color = "white"
+      ) +
+      ggplot2::stat_function(fun = stats::dnorm, color = "#0085CA", linewidth = 0.8) +
+      ggplot2::labs(x = "Pearson residual", y = "Density")
+  }
+  plt +
+    theme_noaa() +
+    comp_facet(c(if ("comp_type" %in% facet_vars) "comp_type", "panel"), ncol = ncol, scales = "free")
+}
+
+#------------------------------------------------------------------------------
+
+#' Plot patterns in composition residuals by bin or year
+#'
+#' Plot the mean Pearson residual for each bin (across years) or each year
+#' (across bins), with a runs test for non-random patterns in the signs of
+#' the residuals. These plots summarize the stripes that can appear in
+#' residual bubble plots, for example, consistently negative residuals at
+#' older ages caused by mis-specified selectivity.
+#'
+#' @inheritParams plot_comp
+#' @param by String. Summarize residuals by "bin" (age or length) or by
+#' "year".
+#'
+#' Default: "bin"
+#'
+#' @returns A `ggplot2` object.
+#'
+#' @details Vertical lines are plus and minus two standard errors of the mean
+#' residual. The runs test (Wald-Wolfowitz) is applied to the signs of the
+#' residuals ordered by year (when `by = "bin"`) or by bin (when
+#' `by = "year"`). Filled points have fewer runs than expected by chance
+#' (one-sided p < 0.05), which indicates a non-random pattern of residuals.
+#' When all residuals have the same sign, the p-value is the probability of
+#' that happening by chance, \eqn{2 \times 0.5^n}.
+#' Residuals are summarized separately for each fleet and sex.
+#'
+#' @seealso [plot_comp_residuals()], [plot_comp_residual_dist()],
+#' [prepare_comp_data()]
+#' @export
+#'
+#' @examples
+#' plot_comp_residual_pattern(bigskate_lencomp)
+#' plot_comp_residual_pattern(bigskate_lencomp, by = "year")
+plot_comp_residual_pattern <- function(
+  dat,
+  by = c("bin", "year"),
+  comp_type = NULL,
+  fleet = NULL,
+  ncol = NULL
+) {
+  by <- match.arg(by)
+  comps <- comp_residuals_or_abort(dat, comp_type = comp_type, fleet = fleet)
+  bin_label <- comp_bin_label(comps)
+  order_var <- if (by == "bin") "year" else "bin"
+
+  summ <- comps |>
+    dplyr::arrange(.data[[order_var]]) |>
+    dplyr::group_by(.data$comp_type, .data$fleet, .data$sex, .data[[by]]) |>
+    dplyr::summarise(
+      mean_resid = mean(.data$pearson),
+      se = stats::sd(.data$pearson) / sqrt(dplyr::n()),
+      runs_p = runs_test(.data$pearson),
+      .groups = "drop"
+    ) |>
+    dplyr::mutate(
+      runs = factor(
+        dplyr::if_else(!is.na(.data$runs_p) & .data$runs_p < 0.05, "Non-random (p < 0.05)", "Random"),
+        levels = c("Non-random (p < 0.05)", "Random")
+      )
+    )
+
+  facet_vars <- c(
+    if (length(unique(summ$comp_type)) > 1) "comp_type",
+    "fleet",
+    if (length(unique(summ$sex)) > 1) "sex"
+  )
+
+  ggplot2::ggplot(summ, ggplot2::aes(x = .data[[by]], y = .data$mean_resid)) +
+    ggplot2::geom_hline(yintercept = 0, color = "gray50") +
+    ggplot2::geom_linerange(
+      ggplot2::aes(
+        ymin = .data$mean_resid - 2 * .data$se,
+        ymax = .data$mean_resid + 2 * .data$se
+      ),
+      color = "gray40",
+      na.rm = TRUE
+    ) +
+    ggplot2::geom_point(ggplot2::aes(fill = .data$runs), shape = 21, size = 2.5) +
+    ggplot2::scale_fill_manual(
+      values = c("Non-random (p < 0.05)" = "gray10", "Random" = "white"),
+      drop = FALSE
+    ) +
+    ggplot2::labs(
+      x = if (by == "bin") bin_label else "Year",
+      y = "Mean Pearson residual",
+      fill = paste("Runs test across", if (by == "bin") "years" else "bins")
+    ) +
+    theme_noaa() +
+    comp_facet(facet_vars, ncol = ncol)
+}
+
+#------------------------------------------------------------------------------
+
+#' Plot Andre's mean age at length for conditional age-at-length data
+#'
+#' For conditional age-at-length data, plot the observed and expected mean
+#' age in each length bin, and the standard error of the mean age, for each
+#' year. These are "Andre's plots" from `r4ss::SSplotComps()` (subplot 12).
+#'
+#' @inheritParams plot_comp
+#' @param dat Conditional age-at-length data in any format accepted by
+#' [prepare_comp_data()] with `conditional = TRUE`: a tidy data frame with a
+#' `length_bin` column (see [bigskate_caal]), or an
+#' `r4ss::SS_output()` list (its `condbase` element is used).
+#' @param interval Number. Coverage of the intervals around the observed
+#' values.
+#'
+#' Default: 0.90
+#'
+#' @returns A `ggplot2` object with one row of panels per year: mean age at
+#' length on the left and standard error of the mean age on the right.
+#'
+#' @details Within each year and length bin, the observed and expected mean
+#' ages are \eqn{\sum_a a\, o_a} and \eqn{\sum_a a\, e_a}, and the standard
+#' errors are the corresponding standard deviations of age divided by
+#' \eqn{\sqrt{N}}, where \eqn{N} is the number of fish aged in the length
+#' bin. Following r4ss, the interval for the mean is the observed mean plus or
+#' minus the normal quantile times the observed standard error, and the
+#' interval for the standard error is based on the chi-square distribution.
+#' Observed values outside the intervals, or expected values that
+#' consistently differ from the observed values, suggest mis-specified growth
+#' or variability in length at age.
+#'
+#' Length bins with no observations or no sample size are omitted. Plot one
+#' fleet at a time with `fleet`, and subset years by filtering `dat`.
+#'
+#' @section r4ss equivalent: `r4ss::SSplotComps()` subplot 12 (Andre's plots of
+#' mean age and standard deviation for conditional age-at-length data).
+#'
+#' @seealso [prepare_comp_data()], [plot_comp_bubble()]
+#' @export
+#'
+#' @examples
+#' plot_comp_andre(
+#'   dplyr::filter(bigskate_caal, year %in% 2016:2018),
+#'   fleet = "WCGBTS"
+#' )
+plot_comp_andre <- function(
+  dat,
+  interval = 0.90,
+  fleet = NULL,
+  ncol = NULL
+) {
+  comps <- prepare_comp_data(dat, comp_type = "age", fleet = fleet, conditional = TRUE)
+  if (length(unique(comps$fleet)) > 1) {
+    cli::cli_alert_info(
+      "Plotting {length(unique(comps$fleet))} fleets; use {.arg fleet} to plot one fleet at a time."
+    )
+  }
+  z <- stats::qnorm(1 - (1 - interval) / 2)
+  stats_df <- comps |>
+    dplyr::group_by(.data$fleet, .data$year, .data$sex, .data$length_bin) |>
+    dplyr::filter(
+      max(.data$observed, na.rm = TRUE) > 1e-4,
+      !is.na(dplyr::first(.data$n)),
+      dplyr::first(.data$n) > 0
+    ) |>
+    dplyr::summarise(
+      N = dplyr::first(.data$n),
+      obs_mean = sum(.data$bin * .data$observed) / sum(.data$observed),
+      obs_sd = sqrt(max(0, sum(.data$bin^2 * .data$observed) / sum(.data$observed) - .data$obs_mean[1]^2)),
+      exp_mean = sum(.data$bin * .data$expected) / sum(.data$expected),
+      exp_sd = sqrt(max(0, sum(.data$bin^2 * .data$expected) / sum(.data$expected) - .data$exp_mean[1]^2)),
+      .groups = "drop"
+    ) |>
+    dplyr::mutate(
+      obs_se = .data$obs_sd / sqrt(.data$N),
+      exp_se = .data$exp_sd / sqrt(.data$N)
+    )
+  if (nrow(stats_df) == 0) {
+    cli::cli_abort("No conditional age-at-length data with observations and sample sizes found.")
+  }
+
+  mean_df <- stats_df |>
+    dplyr::transmute(
+      .data$fleet, .data$year, .data$sex, .data$length_bin,
+      statistic = "Mean age",
+      observed = .data$obs_mean,
+      expected = .data$exp_mean,
+      lower = .data$obs_mean - z * .data$obs_se,
+      upper = .data$obs_mean + z * .data$obs_se
+    )
+  se_df <- stats_df |>
+    dplyr::transmute(
+      .data$fleet, .data$year, .data$sex, .data$length_bin,
+      statistic = "SE of mean age",
+      observed = .data$obs_se,
+      expected = .data$exp_se,
+      lower = dplyr::if_else(
+        .data$N > 1,
+        .data$obs_se * sqrt(pmax(.data$N - 1, 0) / stats::qchisq(1 - (1 - interval) / 2, .data$N)),
+        NA_real_
+      ),
+      upper = dplyr::if_else(
+        .data$N > 1,
+        .data$obs_se * sqrt(pmax(.data$N - 1, 0) / stats::qchisq((1 - interval) / 2, .data$N)),
+        NA_real_
+      )
+    )
+  plot_df <- dplyr::bind_rows(mean_df, se_df) |>
+    dplyr::mutate(statistic = factor(.data$statistic, levels = c("Mean age", "SE of mean age")))
+
+  row_vars <- c(
+    if (length(unique(plot_df$fleet)) > 1) "fleet",
+    if (length(unique(plot_df$sex)) > 1) "sex",
+    "year"
+  )
+
+  ggplot2::ggplot(plot_df, ggplot2::aes(x = .data$length_bin)) +
+    ggplot2::geom_ribbon(
+      ggplot2::aes(ymin = .data$lower, ymax = .data$upper),
+      fill = "gray90",
+      na.rm = TRUE
+    ) +
+    ggplot2::geom_line(ggplot2::aes(y = .data$lower), linetype = "dotted", na.rm = TRUE) +
+    ggplot2::geom_line(ggplot2::aes(y = .data$upper), linetype = "dotted", na.rm = TRUE) +
+    ggplot2::geom_line(
+      ggplot2::aes(y = .data$expected, linetype = "Expected"),
+      color = "#0085CA",
+      linewidth = 0.8,
+      na.rm = TRUE
+    ) +
+    ggplot2::geom_point(ggplot2::aes(y = .data$observed, shape = "Observed"), size = 1.5) +
+    ggplot2::facet_grid(
+      rows = ggplot2::vars(!!!rlang::syms(row_vars)),
+      cols = ggplot2::vars(.data$statistic),
+      scales = "free_y",
+      labeller = ggplot2::labeller(.multi_line = FALSE)
+    ) +
+    ggplot2::expand_limits(y = 0) +
+    ggplot2::labs(
+      x = "Length bin",
+      y = NULL,
+      shape = NULL,
+      linetype = NULL,
+      caption = paste0(
+        "Left: mean age at length (obs. and exp.) with ", round(interval * 100),
+        "% intervals based on adding ", round(z, 2), " SE of mean to the data.\n",
+        "Right: SE of mean age at length (obs. and exp.) with ", round(interval * 100),
+        "% intervals based on the chi-square distribution."
+      )
+    ) +
+    theme_noaa() +
+    ggplot2::theme(
+      strip.text = ggplot2::element_text(size = 9),
+      axis.text = ggplot2::element_text(size = 8)
+    )
+}
+
+#------------------------------------------------------------------------------
+# Internal helpers
+#------------------------------------------------------------------------------
+
+# Prepare comps and keep rows with residuals, or abort with a hint
+comp_residuals_or_abort <- function(dat, comp_type = NULL, fleet = NULL) {
+  comps <- prepare_comp_data(dat, comp_type = comp_type, fleet = fleet)
+  comps <- dplyr::filter(comps, is.finite(.data$pearson))
+  if (nrow(comps) == 0) {
+    cli::cli_abort(c(
+      "No Pearson residuals available.",
+      "i" = "Supply a {.field pearson} column, or {.field expected} values and sample sizes ({.field n}) so residuals can be calculated."
+    ))
+  }
+  comps
+}
+
+# One-sided Wald-Wolfowitz runs test on the signs of x (in order). Returns
+# the p-value for fewer runs than expected (positive autocorrelation).
+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))
+}
+
+# Input and effective (McAllister-Ianelli) sample size by composition
+comp_effn_table <- function(comps) {
+  comps |>
+    dplyr::filter(!is.na(.data$expected)) |>
+    dplyr::group_by(.data$comp_type, .data$fleet, .data$year) |>
+    dplyr::summarise(
+      n = dplyr::first(.data$n),
+      effn = sum(.data$expected * (1 - .data$expected)) /
+        sum((.data$observed - .data$expected)^2),
+      .groups = "drop"
+    )
+}
+
+# Does dat contain Stock Synthesis FIT_*_COMPS tables from convert_output()?
+is_ss3_fit_summary <- function(dat) {
+  is.data.frame(dat) &&
+    all(c("module_name", "label", "estimate") %in% names(dat)) &&
+    any(dat$module_name %in% c("FIT_LEN_COMPS", "FIT_AGE_COMPS", "FIT_SIZE_COMPS") &
+      dat$label %in% "effn")
+}
+
+# Sample sizes from the Stock Synthesis FIT_*_COMPS tables in
+# convert_output() output. Each composition is identified by its module,
+# fleet, year, and other indexing columns; compositions that cannot be
+# separated (e.g., conditional age-at-length, which lack the length bin in
+# the converted output) are dropped.
+comp_sample_size_from_ss3 <- function(dat, comp_type = NULL) {
+  keys <- intersect(
+    c("module_name", "fleet", "year", "sex", "season", "subseason", "month", "time"),
+    names(dat)
+  )
+  fit <- dat[dat$module_name %in% c("FIT_LEN_COMPS", "FIT_AGE_COMPS", "FIT_SIZE_COMPS") &
+    dat$label %in% c("nsamp_adj", "effn", "part", "ageerr"), c(keys, "label", "estimate")]
+  # convert_output() keeps the values from each row of the SS3 table
+  # together, in order, so a new composition starts at each "part" (or
+  # "nsamp_adj" in output without partitions)
+  first_label <- if ("part" %in% fit$label) "part" else "nsamp_adj"
+  fit <- fit |>
+    dplyr::group_by(.data$module_name) |>
+    dplyr::mutate(comp_id = cumsum(.data$label == first_label)) |>
+    dplyr::ungroup()
+  wide <- fit |>
+    dplyr::group_by(.data$module_name, .data$comp_id) |>
+    dplyr::filter(all(c("nsamp_adj", "effn") %in% .data$label), !any(duplicated(.data$label))) |>
+    dplyr::ungroup() |>
+    tidyr::pivot_wider(names_from = "label", values_from = "estimate") |>
+    dplyr::filter(!is.na(.data$year), !is.na(.data$fleet))
+  if (!"part" %in% names(wide)) wide$part <- NA_real_
+  if (!"ageerr" %in% names(wide)) wide$ageerr <- NA_real_
+  wide$comp_type <- dplyr::if_else(wide$module_name == "FIT_AGE_COMPS", "age", "length")
+  if (!is.null(comp_type)) wide <- wide[wide$comp_type == comp_type, , drop = FALSE]
+  # Several compositions for the same fleet, year, sex, and partition are
+  # conditional age-at-length data, which lack the length bin in the
+  # converted output
+  wide <- wide |>
+    dplyr::group_by(dplyr::across(dplyr::all_of(c(keys, "part", "ageerr")))) |>
+    dplyr::mutate(n_comps = dplyr::n()) |>
+    dplyr::ungroup()
+  if (any(wide$n_comps > 1)) {
+    cli::cli_alert_info(
+      "Dropped {sum(wide$n_comps > 1)} composition{?s} that appear to be conditional age-at-length data."
+    )
+  }
+  wide |>
+    dplyr::filter(.data$n_comps == 1) |>
+    dplyr::mutate(
+      year = suppressWarnings(as.numeric(.data$year)),
+      n = .data$nsamp_adj
+    ) |>
+    dplyr::group_by(.data$comp_type) |>
+    dplyr::mutate(fleet = label_partition(.data$fleet, .data$part)) |>
+    dplyr::ungroup() |>
+    dplyr::select("comp_type", "fleet", "year", "n", "effn")
+}
diff --git a/R/utils_comp.R b/R/utils_comp.R
new file mode 100644
index 00000000..ebf02a20
--- /dev/null
+++ b/R/utils_comp.R
@@ -0,0 +1,548 @@
+#' Prepare composition data for plotting
+#'
+#' Reshape age or length composition data (observed and, optionally,
+#' expected proportions) into the tidy format used by the composition
+#' plotting functions ([plot_comp()], [plot_comp_bubble()],
+#' [plot_comp_residuals()], [plot_comp_sex_ratio()], and [plot_comp_francis()]).
+#'
+#' @param dat A data frame of composition data. The following formats are
+#' recognized automatically:
+#' \describe{
+#'   \item{Tidy composition data}{Columns `year`, `bin`, and `observed`, with
+#'   optional columns `fleet`, `sex`, `part` (partition: 0 = whole catch,
+#'   1 = discard, 2 = retained), `expected`, `n` (input or adjusted sample
+#'   size), `pearson`, and `comp_type`. See [bigskate_lencomp].}
+#'   \item{FIMS output}{Output from FIMS (e.g., `FIMS::get_estimates()` or
+#'   [convert_output()] with `model = "fims"`) that contains `observed` and
+#'   `expected` columns. Rows whose `label` contains "comp" are kept and the
+#'   bin is taken from the `age` or `length` column.}
+#'   \item{stockplotr standard output}{Long-format output from
+#'   [convert_output()] where observed and expected compositions are stored in
+#'   the `estimate` column with labels such as
+#'   "age_composition_proportion"/"age_composition_observed" and
+#'   "age_composition_expected"/"age_composition_predicted" (and the
+#'   corresponding "length_" labels).}
+#'   \item{r4ss output}{A list returned by `r4ss::SS_output()` (length and
+#'   age compositions are combined and fleet names are added) or its
+#'   `lendbase` or `agedbase` elements (columns `Yr`, `Fleet`, `Sex`, `Bin`,
+#'   `Obs`, `Exp`, `Pearson`, `Nsamp_adj`, ...). Conditional
+#'   age-at-length rows are dropped from `agedbase` so that only marginal age
+#'   compositions remain, unless `conditional = TRUE`.}
+#' }
+#' For conditional age-at-length data (`conditional = TRUE`), tidy data also
+#' need a `length_bin` column giving the length bin that each age
+#' composition is conditioned on (see [bigskate_caal]).
+#' @param comp_type String. Type of composition to keep when `dat` contains
+#' more than one. Options: "age", "length", or NULL (keep all).
+#'
+#' Default: NULL
+#'
+#' @param fleet Character or numeric vector. Fleet(s) to keep, matching
+#' either the fleet name (e.g., "Fishery_current") or the fleet name with its
+#' partition (e.g., "Fishery_current (retained)"). NULL keeps all fleets.
+#'
+#' Default: NULL
+#'
+#' @param normalize Logical. TRUE/FALSE; rescale observed and expected values
+#' so that they sum to one within each year, fleet, and composition type
+#' (and length bin for conditional data). Sex-specific compositions are
+#' normalized jointly across sexes, unless each sex already sums to one on
+#' its own (single-sex compositions). Use TRUE when `dat` contains counts
+#' rather than proportions.
+#'
+#' Default: TRUE
+#'
+#' @param conditional Logical. TRUE/FALSE; keep conditional age-at-length
+#' data instead of marginal compositions. Conditional compositions are
+#' normalized within each length bin. Used by [plot_comp_andre()] and
+#' [plot_comp_bubble()].
+#'
+#' Default: FALSE
+#'
+#' @returns A tibble with one row per year, fleet, sex, and bin and columns
+#' `comp_type`, `fleet`, `year`, `sex`, `bin`, `observed`, `expected`, `n`,
+#' and `pearson` (plus `length_bin` when `conditional = TRUE`). Columns that
+#' cannot be determined from `dat` are `NA`.
+#'
+#' @details When Pearson residuals are not supplied, they are calculated
+#' assuming a multinomial distribution as
+#' \deqn{r = (o - e) / \sqrt{e (1 - e) / N}}
+#' where \eqn{o} and \eqn{e} are the observed and expected proportions and
+#' \eqn{N} is the sample size in the `n` column. If `n` is not available,
+#' residuals are left as `NA`.
+#'
+#' @seealso [plot_comp()], [plot_comp_bubble()], [plot_comp_residuals()],
+#' [plot_comp_sex_ratio()], [plot_comp_francis()]
+#' @export
+#'
+#' @examples
+#' comps <- prepare_comp_data(bigskate_lencomp, fleet = "WCGBTS")
+#' head(comps)
+#'
+#' # Data in the stockplotr long format are reshaped automatically
+#' long <- data.frame(
+#'   label = rep(c("age_composition_proportion", "age_composition_expected"), each = 3),
+#'   year = 2020,
+#'   fleet = "fleet1",
+#'   age = rep(1:3, 2),
+#'   estimate = c(0.2, 0.5, 0.3, 0.25, 0.45, 0.3)
+#' )
+#' prepare_comp_data(long)
+prepare_comp_data <- function(
+  dat,
+  comp_type = NULL,
+  fleet = NULL,
+  normalize = TRUE,
+  conditional = FALSE
+) {
+  # Full r4ss::SS_output() list: combine length and marginal age comps, or
+  # use the conditional age-at-length data
+  if (!is.data.frame(dat) && is.list(dat) &&
+    any(c("lendbase", "agedbase", "condbase") %in% names(dat))) {
+    dbase <- if (conditional) {
+      dat[["condbase"]]
+    } else {
+      dplyr::bind_rows(dat[["lendbase"]], dat[["agedbase"]])
+    }
+    if (is.null(dbase) || nrow(dbase) == 0) {
+      cli::cli_abort("No {if (conditional) 'conditional age-at-length' else 'composition'} data found in the r4ss output.")
+    }
+    if (!is.null(dat[["FleetNames"]]) && nrow(dbase) > 0) {
+      dbase$Fleet_Name <- dat[["FleetNames"]][dbase$Fleet]
+    }
+    dat <- dbase
+  }
+  if (!is.data.frame(dat)) {
+    cli::cli_abort("{.arg dat} must be a data frame or an {.code r4ss::SS_output()} list.")
+  }
+  # Return early if the data are already prepared
+  prepared <- inherits(dat, "stockplotr_comp")
+  if (prepared) {
+    out <- dat
+    if (conditional != ("length_bin" %in% names(out))) {
+      cli::cli_abort("{.arg dat} was prepared with {.code conditional = {!conditional}}.")
+    }
+  } else {
+    out <- switch(detect_comp_format(dat),
+      "r4ss" = comp_from_r4ss(dat, conditional = conditional),
+      "obs_exp_cols" = comp_from_obs_exp_cols(dat),
+      "long" = comp_from_long(dat),
+      "tidy" = comp_from_tidy(dat),
+      cli::cli_abort(c(
+        "Could not find composition data in {.arg dat}.",
+        "i" = "Supply a data frame with columns {.field year}, {.field bin}, and {.field observed} (see {.code ?prepare_comp_data})."
+      ))
+    )
+    has_lbin <- any(!is.na(out$length_bin))
+    if (conditional && !has_lbin) {
+      cli::cli_abort(c(
+        "No conditional age-at-length data found in {.arg dat}.",
+        "i" = "Tidy data need a {.field length_bin} column."
+      ))
+    }
+    if (!conditional && has_lbin && length(unique(stats::na.omit(out$length_bin))) > 1) {
+      cli::cli_abort(c(
+        "{.arg dat} contains conditional age-at-length data ({.field length_bin}).",
+        "i" = "Use {.code conditional = TRUE}, e.g., with {.fn plot_comp_andre}."
+      ))
+    }
+    if (!conditional) out$length_bin <- NULL
+  }
+
+  # Filter to comp type and fleet
+  if (!is.null(comp_type)) {
+    keep_type <- match.arg(comp_type, c("age", "length"))
+    out$comp_type[is.na(out$comp_type)] <- keep_type
+    out <- out[out$comp_type == keep_type, , drop = FALSE]
+  } else if (length(unique(stats::na.omit(out$comp_type))) > 1) {
+    cli::cli_alert_info(
+      "Both age and length compositions found; use {.arg comp_type} to plot one type at a time."
+    )
+  }
+  if (!is.null(fleet)) {
+    # Match either the fleet name or the fleet name with its partition,
+    # e.g., "Fishery_current" or "Fishery_current (retained)"
+    fleet_base <- base_fleet(out$fleet)
+    keep <- out$fleet %in% as.character(fleet) | fleet_base %in% as.character(fleet)
+    out <- out[keep, , drop = FALSE]
+  }
+  if (nrow(out) == 0) {
+    cli::cli_abort("No composition data remain after filtering by {.arg comp_type} and {.arg fleet}.")
+  }
+
+  if (normalize && !prepared) {
+    # Compositions are normalized jointly across sexes unless each sex
+    # already sums to one on its own (single-sex compositions, as is common
+    # for conditional age-at-length data)
+    out <- out |>
+      dplyr::group_by(dplyr::across(dplyr::any_of(
+        c("comp_type", "fleet", "year", "length_bin")
+      ))) |>
+      dplyr::mutate(
+        sex_group = if (length(unique(.data$sex)) > 1 &&
+          all(abs(tapply(.data$observed, .data$sex, sum, na.rm = TRUE) - 1) < 1e-3)) {
+          .data$sex
+        } else {
+          "joint"
+        }
+      ) |>
+      dplyr::group_by(.data$sex_group, .add = TRUE) |>
+      dplyr::mutate(
+        observed = .data$observed / sum(.data$observed, na.rm = TRUE),
+        expected = .data$expected / sum(.data$expected, na.rm = TRUE)
+      ) |>
+      dplyr::ungroup() |>
+      dplyr::select(-"sex_group") |>
+      # 0/0 when a year has no observations
+      dplyr::mutate(
+        observed = dplyr::if_else(is.nan(.data$observed), NA_real_, .data$observed),
+        expected = dplyr::if_else(is.nan(.data$expected), NA_real_, .data$expected)
+      )
+  }
+
+  # Calculate multinomial Pearson residuals if not supplied
+  need_pearson <- is.na(out$pearson) & !is.na(out$expected) & !is.na(out$n)
+  if (any(need_pearson)) {
+    out$pearson[need_pearson] <- with(
+      out[need_pearson, ],
+      (observed - expected) / sqrt(expected * (1 - expected) / n)
+    )
+    out$pearson[!is.finite(out$pearson)] <- NA_real_
+  }
+
+  out <- out |>
+    dplyr::arrange(dplyr::across(dplyr::any_of(
+      c("comp_type", "fleet", "year", "length_bin", "sex", "bin")
+    )))
+  class(out) <- unique(c("stockplotr_comp", class(tibble::as_tibble(out))))
+  out
+}
+
+#------------------------------------------------------------------------------
+# Internal helpers for composition data
+#------------------------------------------------------------------------------
+
+# Columns returned by prepare_comp_data()
+comp_cols <- c(
+  "comp_type", "fleet", "year", "sex", "bin",
+  "observed", "expected", "n", "pearson"
+)
+
+# Find a column in dat ignoring case; returns NULL if absent
+find_col <- function(dat, options) {
+  nm <- names(dat)
+  for (o in options) {
+    hit <- nm[tolower(nm) == tolower(o)]
+    if (length(hit) > 0) {
+      return(hit[1])
+    }
+  }
+  NULL
+}
+
+# Pull a column (by any of several names) or return a default vector
+pull_col <- function(dat, options, default = NA) {
+  col <- find_col(dat, options)
+  if (is.null(col)) {
+    return(rep(default, nrow(dat)))
+  }
+  dat[[col]]
+}
+
+detect_comp_format <- function(dat) {
+  if (!is.null(find_col(dat, "Yr")) && !is.null(find_col(dat, "Bin")) &&
+    !is.null(find_col(dat, "Obs"))) {
+    return("r4ss")
+  }
+  has_obs <- !is.null(find_col(dat, "observed"))
+  has_bin <- !is.null(find_col(dat, c("bin", "age", "length", "length_bins")))
+  if (has_obs && !is.null(find_col(dat, "label")) &&
+    any(grepl("comp", dat[[find_col(dat, "label")]], ignore.case = TRUE)) && has_bin) {
+    return("obs_exp_cols")
+  }
+  if (!is.null(find_col(dat, "label")) && !is.null(find_col(dat, "estimate")) &&
+    any(grepl("comp", dat[[find_col(dat, "label")]], ignore.case = TRUE) &
+      grepl("observed|proportion|expected|predicted", dat[[find_col(dat, "label")]], ignore.case = TRUE))) {
+    return("long")
+  }
+  if (has_obs && has_bin && !is.null(find_col(dat, "year"))) {
+    return("tidy")
+  }
+  "unknown"
+}
+
+# Standardize sex codes to "female", "male", or "combined"
+standardize_sex <- function(x) {
+  if (all(is.na(x))) {
+    return(rep("combined", length(x)))
+  }
+  x <- tolower(as.character(x))
+  dplyr::case_when(
+    is.na(x) ~ "combined",
+    x %in% c("1", "f", "female", "females") ~ "female",
+    x %in% c("2", "m", "male", "males") ~ "male",
+    x %in% c("0", "3", "both", "combined", "unsexed", "u", "all") ~ "combined",
+    TRUE ~ x
+  )
+}
+
+comp_type_from_label <- function(label) {
+  dplyr::case_when(
+    grepl("len|size", label, ignore.case = TRUE) ~ "length",
+    grepl("age", label, ignore.case = TRUE) ~ "age",
+    TRUE ~ NA_character_
+  )
+}
+
+finish_comp <- function(out) {
+  out <- tibble::as_tibble(out)
+  out$fleet <- as.character(out$fleet)
+  out$sex <- standardize_sex(out$sex)
+  if (!"length_bin" %in% names(out)) out$length_bin <- NA_real_
+  for (col in c("year", "bin", "observed", "expected", "n", "pearson", "length_bin")) {
+    out[[col]] <- suppressWarnings(as.numeric(out[[col]]))
+  }
+  out$fleet[is.na(out$fleet)] <- "All"
+  out <- dplyr::filter(out, !is.na(.data$bin), !is.na(.data$year))
+  out[, c(comp_cols, "length_bin")]
+}
+
+comp_from_tidy <- function(dat) {
+  finish_comp(data.frame(
+    comp_type = pull_col(dat, "comp_type", NA_character_),
+    fleet = label_partition(
+      as.character(pull_col(dat, c("fleet", "fleet_name"), "All")),
+      pull_col(dat, "part")
+    ),
+    year = pull_col(dat, "year"),
+    sex = pull_col(dat, "sex"),
+    bin = pull_col(dat, c("bin", "age", "length", "length_bins")),
+    observed = pull_col(dat, "observed"),
+    expected = pull_col(dat, c("expected", "predicted")),
+    n = pull_col(dat, c("n", "sample_size", "nsamp", "nsamp_adj")),
+    pearson = pull_col(dat, c("pearson", "pearson_residual", "residual")),
+    length_bin = pull_col(dat, c("length_bin", "lbin", "lbin_lo"))
+  ))
+}
+
+comp_from_obs_exp_cols <- function(dat) {
+  label_col <- find_col(dat, "label")
+  dat <- dat[grepl("comp", dat[[label_col]], ignore.case = TRUE), , drop = FALSE]
+  type <- comp_type_from_label(dat[[label_col]])
+  age <- pull_col(dat, "age")
+  len <- pull_col(dat, c("length", "length_bins"))
+  bin <- pull_col(dat, "bin")
+  bin <- dplyr::case_when(
+    !is.na(bin) ~ suppressWarnings(as.numeric(bin)),
+    type == "length" ~ suppressWarnings(as.numeric(len)),
+    TRUE ~ suppressWarnings(as.numeric(age))
+  )
+  finish_comp(data.frame(
+    comp_type = type,
+    fleet = pull_col(dat, c("fleet", "fleet_name", "module_id"), "All"),
+    year = pull_col(dat, c("year", "year_i")),
+    sex = pull_col(dat, "sex"),
+    bin = bin,
+    observed = pull_col(dat, "observed"),
+    expected = pull_col(dat, c("expected", "predicted")),
+    n = pull_col(dat, c("n", "sample_size", "nsamp", "nsamp_adj")),
+    pearson = pull_col(dat, c("pearson", "pearson_residual"))
+  ))
+}
+
+comp_from_long <- function(dat) {
+  label_col <- find_col(dat, "label")
+  dat <- dat[grepl("comp", dat[[label_col]], ignore.case = TRUE), , drop = FALSE]
+  lab <- dat[[label_col]]
+  value_type <- dplyr::case_when(
+    grepl("expected|predicted", lab, ignore.case = TRUE) ~ "expected",
+    grepl("observed|proportion", lab, ignore.case = TRUE) ~ "observed",
+    TRUE ~ NA_character_
+  )
+  type <- comp_type_from_label(lab)
+  bin <- dplyr::if_else(
+    type == "length",
+    suppressWarnings(as.numeric(pull_col(dat, c("length_bins", "length", "bin")))),
+    suppressWarnings(as.numeric(pull_col(dat, c("age", "bin"))))
+  )
+  long <- data.frame(
+    comp_type = type,
+    fleet = as.character(pull_col(dat, c("fleet", "fleet_name"), "All")),
+    year = suppressWarnings(as.numeric(pull_col(dat, "year"))),
+    sex = as.character(pull_col(dat, "sex")),
+    bin = bin,
+    value_type = value_type,
+    value = suppressWarnings(as.numeric(dat[[find_col(dat, "estimate")]]))
+  ) |>
+    dplyr::filter(!is.na(.data$value_type))
+  wide <- long |>
+    tidyr::pivot_wider(
+      names_from = "value_type",
+      values_from = "value",
+      values_fn = function(x) sum(x, na.rm = TRUE)
+    )
+  if (!"observed" %in% names(wide)) {
+    cli::cli_abort("No observed compositions found in {.arg dat}.")
+  }
+  if (!"expected" %in% names(wide)) wide$expected <- NA_real_
+  wide$n <- NA_real_
+  wide$pearson <- NA_real_
+  finish_comp(wide)
+}
+
+comp_from_r4ss <- function(dat, conditional = FALSE) {
+  dat <- as.data.frame(dat)
+  # Drop conditional age-at-length data, keeping marginal age comps
+  lo <- find_col(dat, "Lbin_lo")
+  hi <- find_col(dat, "Lbin_hi")
+  kind <- pull_col(dat, "Kind")
+  if (!is.null(lo) && !is.null(hi)) {
+    is_age <- if (all(is.na(kind))) rep(TRUE, nrow(dat)) else toupper(kind) %in% "AGE"
+    range <- dat[[hi]] - dat[[lo]]
+    multiple_ranges <- any(is_age) && length(unique(range[is_age])) > 1
+    if (conditional) {
+      # Keep age rows conditioned on a length bin narrower than the full range
+      keep <- is_age & (if (multiple_ranges) range < max(range[is_age], na.rm = TRUE) else TRUE)
+      dat <- dat[keep, , drop = FALSE]
+      kind <- kind[keep]
+    } else if (multiple_ranges) {
+      cli::cli_alert_info("Conditional age-at-length data removed; only marginal age compositions are kept.")
+      keep <- !is_age | range == max(range[is_age], na.rm = TRUE)
+      dat <- dat[keep, , drop = FALSE]
+      kind <- kind[keep]
+    }
+  }
+  type <- dplyr::case_when(
+    toupper(kind) %in% c("LEN", "SIZE") ~ "length",
+    toupper(kind) == "AGE" ~ "age",
+    !is.null(lo) ~ "age",
+    TRUE ~ "length"
+  )
+  fleet <- as.character(pull_col(dat, c("Fleet_Name", "Fleet_name", "Fleet"), "All"))
+  # Add partition to fleet name when a fleet has more than one partition
+  fleet <- label_partition(fleet, pull_col(dat, "Part"))
+  out <- finish_comp(data.frame(
+    comp_type = type,
+    fleet = fleet,
+    year = pull_col(dat, "Yr"),
+    sex = r4ss_sex(dat),
+    bin = pull_col(dat, "Bin"),
+    observed = pull_col(dat, "Obs"),
+    expected = pull_col(dat, "Exp"),
+    # Dirichlet-multinomial sample size when available, as in r4ss
+    n = dplyr::coalesce(
+      suppressWarnings(as.numeric(pull_col(dat, "Nsamp_DM"))),
+      suppressWarnings(as.numeric(pull_col(dat, c("Nsamp_adj", "N"))))
+    ),
+    pearson = pull_col(dat, "Pearson"),
+    length_bin = if (conditional && !is.null(lo)) dat[[lo]] else NA_real_
+  ))
+  # Single-sex male compositions are reported twice in the Stock Synthesis
+  # composition database (once for each internal sex code), with identical
+  # values, so keep one copy
+  dplyr::distinct(out)
+}
+
+# Add the partition to fleet names, as in r4ss, when a fleet has data from
+# more than one partition (0 = whole catch, 1 = discard, 2 = retained)
+label_partition <- function(fleet, part) {
+  fleet <- as.character(fleet)
+  if (all(is.na(part))) {
+    return(fleet)
+  }
+  part_name <- c("0" = "whole catch", "1" = "discard", "2" = "retained")[as.character(part)]
+  multi_part <- tapply(part, fleet, function(x) length(unique(stats::na.omit(x))) > 1)
+  ifelse(
+    !is.na(part_name) & multi_part[fleet] %in% TRUE,
+    paste0(fleet, " (", part_name, ")"),
+    fleet
+  )
+}
+
+# Fleet name without the partition added by label_partition()
+base_fleet <- function(fleet) {
+  sub(" \\((whole catch|discard|retained)\\)$", "", fleet)
+}
+
+# Partition name added by label_partition(), or NA
+fleet_partition <- function(fleet) {
+  pattern <- "^.* \\((whole catch|discard|retained)\\)$"
+  ifelse(grepl(pattern, fleet), sub(pattern, "\\1", fleet), NA_character_)
+}
+
+# Sex from an r4ss composition database. r4ss adds a lowercase "sex"
+# column (0 = unknown, 1 = female, 2 = male); the uppercase "Sex" column is
+# an internal Stock Synthesis code where 1 means female or unknown, so
+# single-sex male compositions (Sexes = 2) can have Sex = 1.
+r4ss_sex <- function(dat) {
+  if ("sex" %in% names(dat)) {
+    return(dat[["sex"]])
+  }
+  sexes <- pull_col(dat, "Sexes")
+  sex <- pull_col(dat, "Sex")
+  if (all(is.na(sexes))) {
+    return(sex)
+  }
+  dplyr::if_else(sexes == 3, sex, sexes)
+}
+
+# Axis label for the bin variable
+comp_bin_label <- function(comps) {
+  types <- unique(stats::na.omit(comps$comp_type))
+  if (length(types) == 1 && types == "age") {
+    "Age"
+  } else if (length(types) == 1 && types == "length") {
+    "Length"
+  } else {
+    "Bin"
+  }
+}
+
+# Jeffreys interval for a binomial proportion
+jeffreys_interval <- function(p, n, level = 0.75) {
+  alpha <- (1 - level) / 2
+  x <- p * n
+  lower <- ifelse(x <= 0, 0, stats::qbeta(alpha, x + 0.5, n - x + 0.5))
+  upper <- ifelse(x >= n, 1, stats::qbeta(1 - alpha, x + 0.5, n - x + 0.5))
+  list(lower = lower, upper = upper)
+}
+
+# Aggregate compositions across years within fleet (and sex). Years are
+# weighted by sample size when available.
+aggregate_comp_years <- function(comps) {
+  comps |>
+    dplyr::group_by(.data$comp_type, .data$fleet, .data$year) |>
+    dplyr::mutate(w = dplyr::if_else(is.na(.data$n), 1, .data$n)) |>
+    dplyr::group_by(.data$comp_type, .data$fleet, .data$sex, .data$bin) |>
+    dplyr::summarise(
+      observed = sum(.data$observed * .data$w, na.rm = TRUE),
+      expected = sum(.data$expected * .data$w, na.rm = TRUE),
+      has_exp = any(!is.na(.data$expected)),
+      .groups = "drop"
+    ) |>
+    dplyr::mutate(n = NA_real_) |>
+    dplyr::group_by(.data$comp_type, .data$fleet) |>
+    dplyr::mutate(
+      observed = .data$observed / sum(.data$observed),
+      expected = dplyr::if_else(.data$has_exp, .data$expected / sum(.data$expected), NA_real_)
+    ) |>
+    dplyr::ungroup() |>
+    dplyr::select(-"has_exp")
+}
+
+# Facet a composition plot by the variables that vary in the data
+comp_facet <- function(vars, ncol = NULL, scales = "fixed") {
+  if (length(vars) == 0) {
+    return(NULL)
+  }
+  list(
+    ggplot2::facet_wrap(
+      ggplot2::vars(!!!rlang::syms(vars)),
+      ncol = ncol,
+      scales = scales,
+      labeller = ggplot2::labeller(.multi_line = FALSE)
+    ),
+    # Fleet names with partitions and weights make long panel labels
+    ggplot2::theme(strip.text = ggplot2::element_text(size = 10))
+  )
+}
diff --git a/data-raw/bigskate_comps.R b/data-raw/bigskate_comps.R
new file mode 100644
index 00000000..f2a96d5d
--- /dev/null
+++ b/data-raw/bigskate_comps.R
@@ -0,0 +1,57 @@
+## Code to prepare `bigskate_lencomp` and `bigskate_caal`
+# Composition data and fits from the 2019 Big Skate stock assessment
+# (BigSkate_2019 in the Stock Synthesis test models), the model used for the
+# r4ss example plots at https://r4ss.github.io/r4ss/BigSkate/.
+#
+# The test model includes Report.sso but not CompReport.sso, so the model is
+# rerun at its maximum likelihood estimates (from ss3.par, without
+# estimation) to write the composition database, then read with
+# r4ss::SS_output().
+
+model_dir <- file.path(tempdir(), "bigskate")
+dir.create(model_dir, showWarnings = FALSE)
+
+# Get the model and the Stock Synthesis executable
+r4ss::download_models(dir = model_dir)
+bigskate_dir <- file.path(model_dir, "models", "BigSkate_2019")
+r4ss::get_ss3_exe(dir = bigskate_dir)
+
+# Start from the estimated parameters in ss3.par and do not estimate
+starter <- r4ss::SS_readstarter(file.path(bigskate_dir, "starter.ss"), verbose = FALSE)
+starter[["init_values_src"]] <- 1
+r4ss::SS_writestarter(starter, dir = bigskate_dir, overwrite = TRUE, verbose = FALSE)
+r4ss::run(dir = bigskate_dir, extras = "-maxfn 0 -phase 50 -nohess", exe = "ss3", skipfinished = FALSE)
+
+replist <- r4ss::SS_output(
+  bigskate_dir,
+  verbose = FALSE, printstats = FALSE, covar = FALSE, forecast = FALSE
+)
+
+# Reshape r4ss composition databases to stockplotr's snake_case columns
+tidy_dbase <- function(dbase, comp_type, fleet_names) {
+  out <- tibble::tibble(
+    comp_type = comp_type,
+    fleet = fleet_names[dbase[["Fleet"]]],
+    year = dbase[["Yr"]],
+    month = dbase[["Month"]],
+    # r4ss's lowercase "sex" (0 = unknown, 1 = female, 2 = male), not the
+    # Stock Synthesis "Sex" code, which is 1 for male-only compositions
+    sex = c("1" = "female", "2" = "male", "0" = "combined")[as.character(dbase[["sex"]])],
+    part = dbase[["Part"]],
+    length_bin = dbase[["Lbin_lo"]],
+    bin = dbase[["Bin"]],
+    observed = dbase[["Obs"]],
+    expected = dbase[["Exp"]],
+    pearson = dbase[["Pearson"]],
+    n = dbase[["Nsamp_adj"]]
+  )
+  if (comp_type == "length") out[["length_bin"]] <- NULL
+  # Single-sex male compositions are reported twice in the composition
+  # database (once for each internal sex code) with identical values
+  dplyr::distinct(out)
+}
+
+bigskate_lencomp <- tidy_dbase(replist[["lendbase"]], "length", replist[["FleetNames"]])
+bigskate_caal <- tidy_dbase(replist[["condbase"]], "age", replist[["FleetNames"]])
+
+usethis::use_data(bigskate_lencomp, bigskate_caal, overwrite = TRUE)
diff --git a/data/bigskate_caal.rda b/data/bigskate_caal.rda
new file mode 100644
index 00000000..b9606829
Binary files /dev/null and b/data/bigskate_caal.rda differ
diff --git a/data/bigskate_lencomp.rda b/data/bigskate_lencomp.rda
new file mode 100644
index 00000000..6eb32a50
Binary files /dev/null and b/data/bigskate_lencomp.rda differ
diff --git a/inst/WORDLIST b/inst/WORDLIST
index 2d22da2f..62494c59 100644
--- a/inst/WORDLIST
+++ b/inst/WORDLIST
@@ -1,37 +1,25 @@
-BAM
-BMSY
-Birthseason
-CPUE
-FIMS
-FLIMIT
-FMSY
-FTARGET
-Fmort
-Hennen
-Lifecycle
-MSY
-Modelling
-NOAA
-ORCID
-Petrale
-Pre
-Preformatted
-Rceattle
-Rda
-Rmarkdown
-Rmd
-Spawners
-StockSmart
-Subseason
 aa
+agedbase
 alttext
+andre
 asar
 assmt
+BAM
 bam
+Beringraja
+BigSkate
+bigskate
+binoculata
 birthseas
+Birthseason
+Bizzarro
+BMSY
+caal
 cha
 cmd
+condbase
 conjuction
+CPUE
 csv
 customizations
 cv
@@ -39,41 +27,83 @@ dat
 dataframe's
 diagonostic
 dplyr
+EffN
+effn
 facetting
+FIMS
 fims
 flextable
+FLIMIT
+Fmort
+FMSY
+francis
+FTARGET
+Gertseva
 getwd
 ggplot
 gp
+Hennen
 hline
 href
+Ianelli
 is’
+Jeffreys
 kableExtra
 kobe
 kqs
+Lbin
 len
+lencomp
+lendbase
+Lifecycle
+loess
+McAllister
+mis
+Modelling
+MSY
 msy
 natively
+NOAA
 noaa
+Nsamp
+nsamp
 nsim
 obsvpred
+ORCID
+pearson
+Petrale
 pkgdown
+Pre
+Preformatted
 preformatted
+qq
+Rceattle
 rceattle
+Rda
 rda
 rda's
 rdat
+Rmarkdown
+Rmd
 rowid
 satf
 sd
+SDNR
 se
 seealso
+sexratio
 socio
+Spawners
 ss
+SSMethod
 sso
+SSplotComps
+SSplotSexRatio
 stddev
+StockSmart
 str
 stringr
+Subseason
 subseason
 theming
 tibble
@@ -82,5 +112,7 @@ tibbles
 tidyr
 tidyselect
 unfished
+WCGBTS
+Wolfowitz
 yaml
 yml
diff --git a/man/bigskate_caal.Rd b/man/bigskate_caal.Rd
new file mode 100644
index 00000000..a7549b5a
--- /dev/null
+++ b/man/bigskate_caal.Rd
@@ -0,0 +1,42 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data.R
+\docType{data}
+\name{bigskate_caal}
+\alias{bigskate_caal}
+\title{Big Skate conditional age-at-length data}
+\format{
+A tibble with 6560 rows and 12 variables:
+\describe{
+  \item{comp_type}{Type of composition ("age")}
+  \item{fleet}{Fleet name: "Fishery_current" or "WCGBTS"}
+  \item{year}{Year of the observation (`Yr` in r4ss)}
+  \item{month}{Month of the observation}
+  \item{sex}{"female" or "male"}
+  \item{part}{Partition: 0 = whole catch, 1 = discard, 2 = retained}
+  \item{length_bin}{Length bin the ages are conditioned on (`Lbin_lo` in r4ss)}
+  \item{bin}{Age}
+  \item{observed}{Observed proportion at age within the length bin}
+  \item{expected}{Expected proportion at age within the length bin}
+  \item{pearson}{Pearson residual reported by Stock Synthesis}
+  \item{n}{Adjusted input sample size (`Nsamp_adj` in r4ss)}
+}
+}
+\source{
+Taylor, I.G., Gertseva, V., Stephens, A., and Bizzarro, J. (2019). Status of
+Big Skate (Beringraja binoculata) off the U.S. Pacific Coast in 2019.
+Pacific Fishery Management Council, Portland, OR.
+}
+\usage{
+bigskate_caal
+}
+\description{
+Conditional age-at-length data and fits from the 2019 Big Skate stock
+assessment (see [bigskate_lencomp]). The data are the `condbase` element
+of `r4ss::SS_output()` with columns renamed to follow stockplotr conventions.
+Each age composition is for one sex and one length bin, and sums to one.
+See `data-raw/bigskate_comps.R`.
+}
+\seealso{
+[bigskate_lencomp], [prepare_comp_data()], [plot_comp_andre()]
+}
+\keyword{datasets}
diff --git a/man/bigskate_lencomp.Rd b/man/bigskate_lencomp.Rd
new file mode 100644
index 00000000..9450782a
--- /dev/null
+++ b/man/bigskate_lencomp.Rd
@@ -0,0 +1,47 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data.R
+\docType{data}
+\name{bigskate_lencomp}
+\alias{bigskate_lencomp}
+\title{Big Skate length composition data}
+\format{
+A tibble with 3626 rows and 11 variables:
+\describe{
+  \item{comp_type}{Type of composition ("length")}
+  \item{fleet}{Fleet name: "Fishery_current", "WCGBTS", or "Triennial"}
+  \item{year}{Year of the observation (`Yr` in r4ss)}
+  \item{month}{Month of the observation}
+  \item{sex}{"female" or "male"}
+  \item{part}{Partition: 0 = whole catch, 1 = discard, 2 = retained}
+  \item{bin}{Lower edge of the length bin (cm)}
+  \item{observed}{Observed proportion (`Obs` in r4ss)}
+  \item{expected}{Expected proportion (`Exp` in r4ss)}
+  \item{pearson}{Pearson residual reported by Stock Synthesis}
+  \item{n}{Adjusted input sample size (`Nsamp_adj` in r4ss)}
+}
+}
+\source{
+Taylor, I.G., Gertseva, V., Stephens, A., and Bizzarro, J. (2019). Status of
+Big Skate (Beringraja binoculata) off the U.S. Pacific Coast in 2019.
+Pacific Fishery Management Council, Portland, OR.
+}
+\usage{
+bigskate_lencomp
+}
+\description{
+Length compositions and fits from the 2019 Big Skate stock assessment for
+the U.S. West Coast (`BigSkate_2019` in the Stock Synthesis test models),
+the model used for the r4ss example plots at
+\url{https://r4ss.github.io/r4ss/BigSkate/}. The data are the `lendbase`
+element of `r4ss::SS_output()` with columns renamed to follow
+stockplotr conventions. See `data-raw/bigskate_comps.R`.
+}
+\details{
+Compositions are sex-specific, with female and male proportions summing to
+one within each year, fleet, and partition. `Fishery_current` has separate
+retained and discarded compositions.
+}
+\seealso{
+[bigskate_caal], [prepare_comp_data()], [plot_comp()]
+}
+\keyword{datasets}
diff --git a/man/plot_comp.Rd b/man/plot_comp.Rd
new file mode 100644
index 00000000..88e6e410
--- /dev/null
+++ b/man/plot_comp.Rd
@@ -0,0 +1,92 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_comp.R
+\name{plot_comp}
+\alias{plot_comp}
+\title{Plot composition distributions}
+\usage{
+plot_comp(
+  dat,
+  show_expected = TRUE,
+  aggregate = FALSE,
+  comp_type = NULL,
+  fleet = NULL,
+  show_n = TRUE,
+  ncol = NULL
+)
+}
+\arguments{
+\item{dat}{A data frame of composition data in any format accepted by
+[prepare_comp_data()], such as FIMS output, an r4ss `lendbase`/`agedbase`
+table, or a tidy data frame like [bigskate_lencomp].}
+
+\item{show_expected}{Logical. TRUE/FALSE; overlay the expected values as
+lines. Ignored if `dat` contains no expected values.
+
+Default: TRUE}
+
+\item{aggregate}{Logical. TRUE/FALSE; aggregate compositions across years
+within each fleet. Years are weighted by their sample size when sample
+sizes are available.
+
+Default: FALSE}
+
+\item{comp_type}{String. Type of composition to keep when `dat` contains
+more than one. Options: "age", "length", or NULL (keep all).
+
+Default: NULL}
+
+\item{fleet}{Character or numeric vector. Fleet(s) to keep, matching
+either the fleet name (e.g., "Fishery_current") or the fleet name with its
+partition (e.g., "Fishery_current (retained)"). NULL keeps all fleets.
+
+Default: NULL}
+
+\item{show_n}{Logical. TRUE/FALSE; label each panel with its adjusted
+input sample size ("N adj.") and, when expected values are shown, its
+effective sample size ("N eff."), as in r4ss. Not used when
+`aggregate = TRUE`.
+
+Default: TRUE}
+
+\item{ncol}{Number. Number of facet columns. NULL lets `ggplot2` choose.
+
+Default: NULL}
+}
+\value{
+A `ggplot2` object.
+}
+\description{
+Plot age or length compositions as distributions (proportion by bin), one
+panel per year, optionally overlaying the model-expected values. Setting
+`aggregate = TRUE` instead sums the compositions across years so that each
+fleet is summarized in a single panel.
+}
+\details{
+When compositions are sex-specific (both "female" and "male"
+values in the `sex` column), females are plotted above zero and males
+below zero. Bars show observed proportions and lines show expected
+proportions.
+
+If more than one fleet is plotted without aggregating, panels are
+created for every fleet and year combination; use `fleet` to focus on one
+fleet at a time.
+}
+\section{r4ss equivalent}{
+ `r4ss::SSplotComps()` subplot 1 (by year) and
+subplot 21 (aggregated across years).
+}
+
+\examples{
+# Observations only
+plot_comp(bigskate_lencomp, fleet = "WCGBTS", show_expected = FALSE)
+
+# Observed and expected
+plot_comp(bigskate_lencomp, fleet = "Triennial")
+
+# All fleets aggregated across years
+plot_comp(bigskate_lencomp, aggregate = TRUE)
+}
+\seealso{
+[prepare_comp_data()], [plot_comp_bubble()],
+[plot_comp_residuals()], [plot_comp_sex_ratio()], [plot_comp_francis()]
+}
diff --git a/man/plot_comp_andre.Rd b/man/plot_comp_andre.Rd
new file mode 100644
index 00000000..87daf96b
--- /dev/null
+++ b/man/plot_comp_andre.Rd
@@ -0,0 +1,67 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_comp_diagnostics.R
+\name{plot_comp_andre}
+\alias{plot_comp_andre}
+\title{Plot Andre's mean age at length for conditional age-at-length data}
+\usage{
+plot_comp_andre(dat, interval = 0.9, fleet = NULL, ncol = NULL)
+}
+\arguments{
+\item{dat}{Conditional age-at-length data in any format accepted by
+[prepare_comp_data()] with `conditional = TRUE`: a tidy data frame with a
+`length_bin` column (see [bigskate_caal]), or an
+`r4ss::SS_output()` list (its `condbase` element is used).}
+
+\item{interval}{Number. Coverage of the intervals around the observed
+values.
+
+Default: 0.90}
+
+\item{fleet}{Character or numeric vector. Fleet(s) to keep, matching
+either the fleet name (e.g., "Fishery_current") or the fleet name with its
+partition (e.g., "Fishery_current (retained)"). NULL keeps all fleets.
+
+Default: NULL}
+
+\item{ncol}{Number. Number of facet columns. NULL lets `ggplot2` choose.
+
+Default: NULL}
+}
+\value{
+A `ggplot2` object with one row of panels per year: mean age at
+length on the left and standard error of the mean age on the right.
+}
+\description{
+For conditional age-at-length data, plot the observed and expected mean
+age in each length bin, and the standard error of the mean age, for each
+year. These are "Andre's plots" from `r4ss::SSplotComps()` (subplot 12).
+}
+\details{
+Within each year and length bin, the observed and expected mean
+ages are \eqn{\sum_a a\, o_a} and \eqn{\sum_a a\, e_a}, and the standard
+errors are the corresponding standard deviations of age divided by
+\eqn{\sqrt{N}}, where \eqn{N} is the number of fish aged in the length
+bin. Following r4ss, the interval for the mean is the observed mean plus or
+minus the normal quantile times the observed standard error, and the
+interval for the standard error is based on the chi-square distribution.
+Observed values outside the intervals, or expected values that
+consistently differ from the observed values, suggest mis-specified growth
+or variability in length at age.
+
+Length bins with no observations or no sample size are omitted. Plot one
+fleet at a time with `fleet`, and subset years by filtering `dat`.
+}
+\section{r4ss equivalent}{
+ `r4ss::SSplotComps()` subplot 12 (Andre's plots of
+mean age and standard deviation for conditional age-at-length data).
+}
+
+\examples{
+plot_comp_andre(
+  dplyr::filter(bigskate_caal, year \%in\% 2016:2018),
+  fleet = "WCGBTS"
+)
+}
+\seealso{
+[prepare_comp_data()], [plot_comp_bubble()]
+}
diff --git a/man/plot_comp_bubble.Rd b/man/plot_comp_bubble.Rd
new file mode 100644
index 00000000..3e87ce6b
--- /dev/null
+++ b/man/plot_comp_bubble.Rd
@@ -0,0 +1,91 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_comp.R
+\name{plot_comp_bubble}
+\alias{plot_comp_bubble}
+\title{Plot composition bubble plots}
+\usage{
+plot_comp_bubble(
+  dat,
+  value = c("observed", "pearson", "expected"),
+  comp_type = NULL,
+  fleet = NULL,
+  max_size = 6,
+  ncol = NULL,
+  conditional = FALSE
+)
+}
+\arguments{
+\item{dat}{A data frame of composition data in any format accepted by
+[prepare_comp_data()], such as FIMS output, an r4ss `lendbase`/`agedbase`
+table, or a tidy data frame like [bigskate_lencomp].}
+
+\item{value}{String. Value to plot. Options: "observed", "expected", or
+"pearson".
+
+Default: "observed"}
+
+\item{comp_type}{String. Type of composition to keep when `dat` contains
+more than one. Options: "age", "length", or NULL (keep all).
+
+Default: NULL}
+
+\item{fleet}{Character or numeric vector. Fleet(s) to keep, matching
+either the fleet name (e.g., "Fishery_current") or the fleet name with its
+partition (e.g., "Fishery_current (retained)"). NULL keeps all fleets.
+
+Default: NULL}
+
+\item{max_size}{Number. Size of the largest bubble.
+
+Default: 6}
+
+\item{ncol}{Number. Number of facet columns. NULL lets `ggplot2` choose.
+
+Default: NULL}
+
+\item{conditional}{Logical. TRUE/FALSE; plot conditional age-at-length
+data (see [prepare_comp_data()]), with age on the y-axis, length bin on
+the x-axis, and one panel per year.
+
+Default: FALSE}
+}
+\value{
+A `ggplot2` object.
+}
+\description{
+Bubble plot of compositions by year (x-axis) and bin (y-axis). Bubbles can
+show the observed (or expected) proportions, which makes strong and weak
+cohorts easy to follow through time, or the Pearson residuals, which show
+blocks or stripes of misfit among years or bins.
+}
+\details{
+Panels are created for each fleet, and for each sex when
+compositions are sex-specific. For Pearson residuals, filled bubbles are
+positive residuals (observed > expected) and open bubbles are negative
+residuals (observed < expected), with bubble area proportional to the
+absolute value of the residual. The maximum absolute residual for each
+panel is included in the panel label.
+}
+\section{r4ss equivalent}{
+ `r4ss::SSplotComps()` subplot 2 (bubble plot of
+observations or residuals) and subplot 3 (conditional age-at-length).
+}
+
+\examples{
+# Observed proportions
+plot_comp_bubble(bigskate_lencomp)
+
+# Pearson residuals
+plot_comp_bubble(bigskate_lencomp, value = "pearson")
+
+# Residuals for conditional age-at-length data
+plot_comp_bubble(
+  dplyr::filter(bigskate_caal, year \%in\% 2016:2018),
+  value = "pearson",
+  fleet = "WCGBTS",
+  conditional = TRUE
+)
+}
+\seealso{
+[prepare_comp_data()], [plot_comp_residuals()], [plot_comp()]
+}
diff --git a/man/plot_comp_francis.Rd b/man/plot_comp_francis.Rd
new file mode 100644
index 00000000..97585fb5
--- /dev/null
+++ b/man/plot_comp_francis.Rd
@@ -0,0 +1,92 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_comp.R
+\name{plot_comp_francis}
+\alias{plot_comp_francis}
+\title{Plot mean age or length with Francis (TA1.8) data weights}
+\usage{
+plot_comp_francis(
+  dat,
+  interval = 0.95,
+  show_francis = TRUE,
+  comp_type = NULL,
+  fleet = NULL,
+  ncol = NULL
+)
+}
+\arguments{
+\item{dat}{A data frame of composition data in any format accepted by
+[prepare_comp_data()], such as FIMS output, an r4ss `lendbase`/`agedbase`
+table, or a tidy data frame like [bigskate_lencomp].}
+
+\item{interval}{Number. Coverage of the interval around the observed mean.
+
+Default: 0.95}
+
+\item{show_francis}{Logical. TRUE/FALSE; add the Francis TA1.8 weight for
+each fleet to the panel label and show thinner, capped intervals
+representing sample sizes after applying that weight.
+
+Default: TRUE}
+
+\item{comp_type}{String. Type of composition to keep when `dat` contains
+more than one. Options: "age", "length", or NULL (keep all).
+
+Default: NULL}
+
+\item{fleet}{Character or numeric vector. Fleet(s) to keep, matching
+either the fleet name (e.g., "Fishery_current") or the fleet name with its
+partition (e.g., "Fishery_current (retained)"). NULL keeps all fleets.
+
+Default: NULL}
+
+\item{ncol}{Number. Number of facet columns. NULL lets `ggplot2` choose.
+
+Default: NULL}
+}
+\value{
+A `ggplot2` object.
+}
+\description{
+Time series of the mean age or length in the observed and expected
+compositions. Developed as a diagnostic for the Francis (2011) data
+weighting method (TA1.8), the plot is also a useful summary of trends in
+compositions and lack of fit among years.
+}
+\details{
+For each year and fleet (with sexes combined), the observed mean
+is \eqn{\bar{x}_o = \sum_b x_b o_b} and the expected mean is
+\eqn{\bar{x}_e = \sum_b x_b e_b}, where \eqn{x_b} is the bin and \eqn{o_b}
+and \eqn{e_b} are the observed and expected proportions. Assuming a
+multinomial distribution, the standard error of the mean is
+\eqn{\sqrt{v_e / N}}, where \eqn{v_e = \sum_b x_b^2 e_b - \bar{x}_e^2} is
+the variance of the expected composition and \eqn{N} is the sample size.
+Intervals are drawn around the observed means.
+
+The Francis weight is \eqn{1 / \mathrm{var}(r)} where
+\eqn{r = (\bar{x}_o - \bar{x}_e) / SE} are the standardized residuals of
+the mean across years. As in Stock Synthesis and
+`r4ss::SSMethod.TA1.8()`, all partitions of a fleet (e.g., retained and
+discarded) share one weight, and its 95% interval is calculated by
+bootstrapping the standardized residuals (1000 resamples), so use
+[set.seed()] for reproducible intervals. The weight assumes a multinomial likelihood and has not been
+adapted for Dirichlet-multinomial or other likelihoods.
+
+Without expected values, only the observed means are plotted. Without
+sample sizes, no intervals or weights are shown.
+}
+\section{r4ss equivalent}{
+ `r4ss::SSplotComps()` subplots 8 and 9, which
+use `r4ss::SSMethod.TA1.8()`.
+}
+
+\examples{
+plot_comp_francis(bigskate_lencomp)
+}
+\references{
+Francis, R.I.C.C. (2011). Data weighting in statistical
+fisheries stock assessment models. Canadian Journal of Fisheries and
+Aquatic Sciences 68: 1124-1138.
+}
+\seealso{
+[prepare_comp_data()], [plot_comp()]
+}
diff --git a/man/plot_comp_residual_dist.Rd b/man/plot_comp_residual_dist.Rd
new file mode 100644
index 00000000..b057510a
--- /dev/null
+++ b/man/plot_comp_residual_dist.Rd
@@ -0,0 +1,63 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_comp_diagnostics.R
+\name{plot_comp_residual_dist}
+\alias{plot_comp_residual_dist}
+\title{Plot the distribution of composition residuals}
+\usage{
+plot_comp_residual_dist(
+  dat,
+  type = c("qq", "histogram"),
+  comp_type = NULL,
+  fleet = NULL,
+  ncol = NULL
+)
+}
+\arguments{
+\item{dat}{A data frame of composition data in any format accepted by
+[prepare_comp_data()], such as FIMS output, an r4ss `lendbase`/`agedbase`
+table, or a tidy data frame like [bigskate_lencomp].}
+
+\item{type}{String. Type of plot. Options: "qq" (normal quantile-quantile
+plot) or "histogram" (with a standard normal density overlaid).
+
+Default: "qq"}
+
+\item{comp_type}{String. Type of composition to keep when `dat` contains
+more than one. Options: "age", "length", or NULL (keep all).
+
+Default: NULL}
+
+\item{fleet}{Character or numeric vector. Fleet(s) to keep, matching
+either the fleet name (e.g., "Fishery_current") or the fleet name with its
+partition (e.g., "Fishery_current (retained)"). NULL keeps all fleets.
+
+Default: NULL}
+
+\item{ncol}{Number. Number of facet columns. NULL lets `ggplot2` choose.
+
+Default: NULL}
+}
+\value{
+A `ggplot2` object.
+}
+\description{
+Quantile-quantile plot or histogram of Pearson residuals for each fleet,
+to check for heavy tails, skewness, and outliers.
+}
+\details{
+The standard deviation of the normalized residuals (SDNR) for
+each fleet is included in the panel label. Values well above one suggest
+that the compositions are given more weight (larger sample sizes) than the
+fit supports. Pearson residuals for compositions are not expected to be
+exactly normally distributed, especially for bins with small expected
+proportions or small sample sizes, so treat departures in the tails with
+some caution.
+}
+\examples{
+plot_comp_residual_dist(bigskate_lencomp)
+plot_comp_residual_dist(bigskate_lencomp, type = "histogram")
+}
+\seealso{
+[plot_comp_residuals()], [plot_comp_residual_pattern()],
+[prepare_comp_data()]
+}
diff --git a/man/plot_comp_residual_pattern.Rd b/man/plot_comp_residual_pattern.Rd
new file mode 100644
index 00000000..f14db601
--- /dev/null
+++ b/man/plot_comp_residual_pattern.Rd
@@ -0,0 +1,67 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_comp_diagnostics.R
+\name{plot_comp_residual_pattern}
+\alias{plot_comp_residual_pattern}
+\title{Plot patterns in composition residuals by bin or year}
+\usage{
+plot_comp_residual_pattern(
+  dat,
+  by = c("bin", "year"),
+  comp_type = NULL,
+  fleet = NULL,
+  ncol = NULL
+)
+}
+\arguments{
+\item{dat}{A data frame of composition data in any format accepted by
+[prepare_comp_data()], such as FIMS output, an r4ss `lendbase`/`agedbase`
+table, or a tidy data frame like [bigskate_lencomp].}
+
+\item{by}{String. Summarize residuals by "bin" (age or length) or by
+"year".
+
+Default: "bin"}
+
+\item{comp_type}{String. Type of composition to keep when `dat` contains
+more than one. Options: "age", "length", or NULL (keep all).
+
+Default: NULL}
+
+\item{fleet}{Character or numeric vector. Fleet(s) to keep, matching
+either the fleet name (e.g., "Fishery_current") or the fleet name with its
+partition (e.g., "Fishery_current (retained)"). NULL keeps all fleets.
+
+Default: NULL}
+
+\item{ncol}{Number. Number of facet columns. NULL lets `ggplot2` choose.
+
+Default: NULL}
+}
+\value{
+A `ggplot2` object.
+}
+\description{
+Plot the mean Pearson residual for each bin (across years) or each year
+(across bins), with a runs test for non-random patterns in the signs of
+the residuals. These plots summarize the stripes that can appear in
+residual bubble plots, for example, consistently negative residuals at
+older ages caused by mis-specified selectivity.
+}
+\details{
+Vertical lines are plus and minus two standard errors of the mean
+residual. The runs test (Wald-Wolfowitz) is applied to the signs of the
+residuals ordered by year (when `by = "bin"`) or by bin (when
+`by = "year"`). Filled points have fewer runs than expected by chance
+(one-sided p < 0.05), which indicates a non-random pattern of residuals.
+When all residuals have the same sign, the p-value is the probability of
+that happening by chance, \eqn{2 \times 0.5^n}.
+Residuals are summarized separately for each fleet and sex.
+}
+\examples{
+plot_comp_residual_pattern(bigskate_lencomp)
+plot_comp_residual_pattern(bigskate_lencomp, by = "year")
+}
+\seealso{
+[plot_comp_residuals()], [plot_comp_residual_dist()],
+[prepare_comp_data()]
+}
diff --git a/man/plot_comp_residuals.Rd b/man/plot_comp_residuals.Rd
new file mode 100644
index 00000000..84ad13c5
--- /dev/null
+++ b/man/plot_comp_residuals.Rd
@@ -0,0 +1,40 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_comp.R
+\name{plot_comp_residuals}
+\alias{plot_comp_residuals}
+\title{Plot composition Pearson residuals}
+\usage{
+plot_comp_residuals(dat, ...)
+}
+\arguments{
+\item{dat}{A data frame of composition data in any format accepted by
+[prepare_comp_data()], such as FIMS output, an r4ss `lendbase`/`agedbase`
+table, or a tidy data frame like [bigskate_lencomp].}
+
+\item{...}{Additional arguments passed to [plot_comp_bubble()].}
+}
+\value{
+A `ggplot2` object.
+}
+\description{
+A convenience wrapper around [plot_comp_bubble()] with
+`value = "pearson"`. Residuals are taken from `dat` when supplied or
+calculated by [prepare_comp_data()] assuming a multinomial distribution.
+}
+\section{r4ss equivalent}{
+ `r4ss::SSplotComps()` subplot 2 for residuals by
+fleet and subplot 24 (residuals compared across fleets).
+}
+
+\examples{
+plot_comp_residuals(bigskate_lencomp)
+
+# Compare residuals across fleets for females
+plot_comp_residuals(
+  dplyr::filter(bigskate_lencomp, sex == "female"),
+  ncol = 1
+)
+}
+\seealso{
+[plot_comp_bubble()], [prepare_comp_data()]
+}
diff --git a/man/plot_comp_sample_size.Rd b/man/plot_comp_sample_size.Rd
new file mode 100644
index 00000000..3fa3d042
--- /dev/null
+++ b/man/plot_comp_sample_size.Rd
@@ -0,0 +1,90 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_comp_diagnostics.R
+\name{plot_comp_sample_size}
+\alias{plot_comp_sample_size}
+\title{Plot input versus effective sample size of compositions}
+\usage{
+plot_comp_sample_size(
+  dat,
+  comp_type = NULL,
+  fleet = NULL,
+  show_years = FALSE,
+  ncol = NULL
+)
+}
+\arguments{
+\item{dat}{Composition data in any format accepted by
+[prepare_comp_data()], or [convert_output()] output from Stock Synthesis
+that contains the `FIT_LEN_COMPS` or `FIT_AGE_COMPS` tables (as in
+`stockplotr::example_data`), from which the adjusted input (`nsamp_adj`)
+and effective (`effn`) sample sizes are taken.}
+
+\item{comp_type}{String. Type of composition to keep when `dat` contains
+more than one. Options: "age", "length", or NULL (keep all).
+
+Default: NULL}
+
+\item{fleet}{Character or numeric vector. Fleet(s) to keep, matching
+either the fleet name (e.g., "Fishery_current") or the fleet name with its
+partition (e.g., "Fishery_current (retained)"). NULL keeps all fleets.
+
+Default: NULL}
+
+\item{show_years}{Logical. TRUE/FALSE; label each point with its year.
+
+Default: FALSE}
+
+\item{ncol}{Number. Number of facet columns. NULL lets `ggplot2` choose.
+
+Default: NULL}
+}
+\value{
+A `ggplot2` object.
+}
+\description{
+Compare the input (or adjusted input) sample size of each composition with
+its effective sample size, the sample size implied by how well the model
+fits the composition. This is the diagnostic behind the McAllister and
+Ianelli (1997) data weighting method.
+}
+\details{
+For each composition (year and fleet), the effective sample size
+is
+\deqn{N_{eff} = \sum_b e_b (1 - e_b) / \sum_b (o_b - e_b)^2}
+where \eqn{o_b} and \eqn{e_b} are the observed and expected proportions in
+bin \eqn{b} (McAllister and Ianelli 1997), the same calculation used by
+Stock Synthesis.
+
+Each panel shows a 1:1 line, a loess smoother (dashed red, when there are
+more than six distinct input sample sizes), the arithmetic mean of the
+input sample sizes (vertical dashed line), and the harmonic mean of the
+effective sample sizes (horizontal dashed line). The ratio of the harmonic
+mean effective sample size to the mean input sample size, the
+McAllister-Ianelli weight, is included in each panel label. Ratios below
+one suggest that input sample sizes are too large relative to the fit. As
+in Stock Synthesis, the weight is calculated across all partitions of a
+fleet; multiplying it by the current variance adjustment gives the
+`Recommend_var_adj` reported by Stock Synthesis.
+
+Effective sample sizes from this formula are not appropriate for
+compositions fit with a Dirichlet-multinomial likelihood.
+}
+\section{r4ss equivalent}{
+ `r4ss::SSplotComps()` subplot 7 ("N-EffN
+comparison").
+}
+
+\examples{
+plot_comp_sample_size(bigskate_lencomp)
+
+# Stock Synthesis output converted with convert_output()
+plot_comp_sample_size(stockplotr::example_data, comp_type = "length")
+}
+\references{
+McAllister, M.K., and Ianelli, J.N. (1997). Bayesian stock
+assessment using catch-age data and the sampling-importance resampling
+algorithm. Canadian Journal of Fisheries and Aquatic Sciences 54: 284-300.
+}
+\seealso{
+[prepare_comp_data()], [plot_comp_francis()]
+}
diff --git a/man/plot_comp_sex_ratio.Rd b/man/plot_comp_sex_ratio.Rd
new file mode 100644
index 00000000..65ca2579
--- /dev/null
+++ b/man/plot_comp_sex_ratio.Rd
@@ -0,0 +1,80 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_comp.R
+\name{plot_comp_sex_ratio}
+\alias{plot_comp_sex_ratio}
+\title{Plot sex ratios implied by composition data}
+\usage{
+plot_comp_sex_ratio(
+  dat,
+  aggregate = FALSE,
+  interval = 0.75,
+  comp_type = NULL,
+  fleet = NULL,
+  ncol = NULL
+)
+}
+\arguments{
+\item{dat}{A data frame of composition data in any format accepted by
+[prepare_comp_data()], such as FIMS output, an r4ss `lendbase`/`agedbase`
+table, or a tidy data frame like [bigskate_lencomp].}
+
+\item{aggregate}{Logical. TRUE/FALSE; aggregate compositions across years
+within each fleet. Years are weighted by their sample size when sample
+sizes are available.
+
+Default: FALSE}
+
+\item{interval}{Number. Coverage of the interval around the observed sex
+ratio, calculated as a Jeffreys interval using the sample size in each bin
+(the sample size multiplied by the proportion of the composition in the
+bin). Set to NULL to omit intervals. Intervals are only shown when sample
+sizes are available.
+
+Default: 0.75}
+
+\item{comp_type}{String. Type of composition to keep when `dat` contains
+more than one. Options: "age", "length", or NULL (keep all).
+
+Default: NULL}
+
+\item{fleet}{Character or numeric vector. Fleet(s) to keep, matching
+either the fleet name (e.g., "Fishery_current") or the fleet name with its
+partition (e.g., "Fishery_current (retained)"). NULL keeps all fleets.
+
+Default: NULL}
+
+\item{ncol}{Number. Number of facet columns. NULL lets `ggplot2` choose.
+
+Default: NULL}
+}
+\value{
+A `ggplot2` object.
+}
+\description{
+Plot the proportion female, female / (female + male), by bin for the
+observed and expected compositions. Trends in the lack of fit to sex ratios
+can be easier to see here than in residuals, especially for bins with
+small proportions.
+}
+\details{
+Requires sex-specific compositions with both "female" and
+"male" values in the `sex` column of `dat`. Points show observed ratios
+and lines show expected ratios. As in r4ss, the smallest observed
+proportion in each composition (such as the small constant Stock Synthesis
+adds to every bin) is subtracted from the observed and expected values
+first, so bins with no fish observed are omitted. When `aggregate = TRUE`, female and male
+compositions are summed across years (weighted by sample size when
+available) before calculating ratios.
+}
+\section{r4ss equivalent}{
+ `r4ss::SSplotSexRatio()` with
+`sexratio.option = 2` (female:total).
+}
+
+\examples{
+plot_comp_sex_ratio(bigskate_lencomp, fleet = "WCGBTS")
+plot_comp_sex_ratio(bigskate_lencomp, aggregate = TRUE)
+}
+\seealso{
+[prepare_comp_data()], [plot_comp()]
+}
diff --git a/man/prepare_comp_data.Rd b/man/prepare_comp_data.Rd
new file mode 100644
index 00000000..0941e908
--- /dev/null
+++ b/man/prepare_comp_data.Rd
@@ -0,0 +1,108 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils_comp.R
+\name{prepare_comp_data}
+\alias{prepare_comp_data}
+\title{Prepare composition data for plotting}
+\usage{
+prepare_comp_data(
+  dat,
+  comp_type = NULL,
+  fleet = NULL,
+  normalize = TRUE,
+  conditional = FALSE
+)
+}
+\arguments{
+\item{dat}{A data frame of composition data. The following formats are
+recognized automatically:
+\describe{
+  \item{Tidy composition data}{Columns `year`, `bin`, and `observed`, with
+  optional columns `fleet`, `sex`, `part` (partition: 0 = whole catch,
+  1 = discard, 2 = retained), `expected`, `n` (input or adjusted sample
+  size), `pearson`, and `comp_type`. See [bigskate_lencomp].}
+  \item{FIMS output}{Output from FIMS (e.g., `FIMS::get_estimates()` or
+  [convert_output()] with `model = "fims"`) that contains `observed` and
+  `expected` columns. Rows whose `label` contains "comp" are kept and the
+  bin is taken from the `age` or `length` column.}
+  \item{stockplotr standard output}{Long-format output from
+  [convert_output()] where observed and expected compositions are stored in
+  the `estimate` column with labels such as
+  "age_composition_proportion"/"age_composition_observed" and
+  "age_composition_expected"/"age_composition_predicted" (and the
+  corresponding "length_" labels).}
+  \item{r4ss output}{A list returned by `r4ss::SS_output()` (length and
+  age compositions are combined and fleet names are added) or its
+  `lendbase` or `agedbase` elements (columns `Yr`, `Fleet`, `Sex`, `Bin`,
+  `Obs`, `Exp`, `Pearson`, `Nsamp_adj`, ...). Conditional
+  age-at-length rows are dropped from `agedbase` so that only marginal age
+  compositions remain, unless `conditional = TRUE`.}
+}
+For conditional age-at-length data (`conditional = TRUE`), tidy data also
+need a `length_bin` column giving the length bin that each age
+composition is conditioned on (see [bigskate_caal]).}
+
+\item{comp_type}{String. Type of composition to keep when `dat` contains
+more than one. Options: "age", "length", or NULL (keep all).
+
+Default: NULL}
+
+\item{fleet}{Character or numeric vector. Fleet(s) to keep, matching
+either the fleet name (e.g., "Fishery_current") or the fleet name with its
+partition (e.g., "Fishery_current (retained)"). NULL keeps all fleets.
+
+Default: NULL}
+
+\item{normalize}{Logical. TRUE/FALSE; rescale observed and expected values
+so that they sum to one within each year, fleet, and composition type
+(and length bin for conditional data). Sex-specific compositions are
+normalized jointly across sexes, unless each sex already sums to one on
+its own (single-sex compositions). Use TRUE when `dat` contains counts
+rather than proportions.
+
+Default: TRUE}
+
+\item{conditional}{Logical. TRUE/FALSE; keep conditional age-at-length
+data instead of marginal compositions. Conditional compositions are
+normalized within each length bin. Used by [plot_comp_andre()] and
+[plot_comp_bubble()].
+
+Default: FALSE}
+}
+\value{
+A tibble with one row per year, fleet, sex, and bin and columns
+`comp_type`, `fleet`, `year`, `sex`, `bin`, `observed`, `expected`, `n`,
+and `pearson` (plus `length_bin` when `conditional = TRUE`). Columns that
+cannot be determined from `dat` are `NA`.
+}
+\description{
+Reshape age or length composition data (observed and, optionally,
+expected proportions) into the tidy format used by the composition
+plotting functions ([plot_comp()], [plot_comp_bubble()],
+[plot_comp_residuals()], [plot_comp_sex_ratio()], and [plot_comp_francis()]).
+}
+\details{
+When Pearson residuals are not supplied, they are calculated
+assuming a multinomial distribution as
+\deqn{r = (o - e) / \sqrt{e (1 - e) / N}}
+where \eqn{o} and \eqn{e} are the observed and expected proportions and
+\eqn{N} is the sample size in the `n` column. If `n` is not available,
+residuals are left as `NA`.
+}
+\examples{
+comps <- prepare_comp_data(bigskate_lencomp, fleet = "WCGBTS")
+head(comps)
+
+# Data in the stockplotr long format are reshaped automatically
+long <- data.frame(
+  label = rep(c("age_composition_proportion", "age_composition_expected"), each = 3),
+  year = 2020,
+  fleet = "fleet1",
+  age = rep(1:3, 2),
+  estimate = c(0.2, 0.5, 0.3, 0.25, 0.45, 0.3)
+)
+prepare_comp_data(long)
+}
+\seealso{
+[plot_comp()], [plot_comp_bubble()], [plot_comp_residuals()],
+[plot_comp_sex_ratio()], [plot_comp_francis()]
+}
diff --git a/pkgdown/_pkgdown.yml b/pkgdown/_pkgdown.yml
index fd62ae0e..51f9f217 100644
--- a/pkgdown/_pkgdown.yml
+++ b/pkgdown/_pkgdown.yml
@@ -27,6 +27,8 @@ navbar:
         href: articles/how-caps-alttext-are-made.html
       - text: Standardizing Assessment Model Output
         href: articles/convert_output_details.html
+      - text: Composition data diagnostics
+        href: articles/composition_diagnostics.html
       - text: Automating data extraction for SIS
         href: articles/extract-sis-data.html
         
diff --git a/tests/testthat/test-plot_comp.R b/tests/testthat/test-plot_comp.R
new file mode 100644
index 00000000..d8fed15f
--- /dev/null
+++ b/tests/testthat/test-plot_comp.R
@@ -0,0 +1,233 @@
+r4ss_output <- readRDS(
+  test_path("fixtures", "r4ss_output", "simple_small", "simple_r4ss.rds")
+)
+
+has_layer <- function(p, geom_class) {
+  any(vapply(p$layers, function(l) inherits(l$geom, geom_class), logical(1)))
+}
+
+# prepare_comp_data ----------------------------------------------------------
+test_that("prepare_comp_data returns standard columns from tidy data", {
+  comps <- prepare_comp_data(stockplotr::bigskate_lencomp)
+  expect_s3_class(comps, "stockplotr_comp")
+  expect_named(
+    comps,
+    c("comp_type", "fleet", "year", "sex", "bin", "observed", "expected", "n", "pearson")
+  )
+  # proportions sum to one across sexes within year, fleet, and partition
+  sums <- comps |>
+    dplyr::summarise(total = sum(observed), .by = c(fleet, year))
+  expect_true(all(abs(sums$total - 1) < 1e-6))
+  # Stock Synthesis Pearson residuals are kept
+  expect_equal(sort(comps$pearson), sort(stockplotr::bigskate_lencomp$pearson))
+})
+
+test_that("prepare_comp_data calculates multinomial Pearson residuals", {
+  comps <- prepare_comp_data(
+    dplyr::select(stockplotr::bigskate_lencomp, -pearson)
+  )
+  i <- 10
+  expect_equal(
+    comps$pearson[i],
+    (comps$observed[i] - comps$expected[i]) /
+      sqrt(comps$expected[i] * (1 - comps$expected[i]) / comps$n[i])
+  )
+})
+
+test_that("prepare_comp_data filters by fleet and comp_type", {
+  comps <- prepare_comp_data(stockplotr::bigskate_lencomp, fleet = "WCGBTS")
+  expect_equal(unique(comps$fleet), "WCGBTS")
+  expect_error(
+    prepare_comp_data(stockplotr::bigskate_lencomp, fleet = "not a fleet"),
+    "No composition data"
+  )
+  expect_error(
+    prepare_comp_data(stockplotr::bigskate_lencomp, comp_type = "age"),
+    "No composition data"
+  )
+})
+
+test_that("partitions are added to fleet names as in r4ss", {
+  comps <- prepare_comp_data(stockplotr::bigskate_lencomp)
+  expect_setequal(
+    unique(comps$fleet),
+    c("Fishery_current (discard)", "Fishery_current (retained)", "WCGBTS", "Triennial")
+  )
+  # a fleet can be selected with or without its partition
+  both <- prepare_comp_data(stockplotr::bigskate_lencomp, fleet = "Fishery_current")
+  expect_equal(length(unique(both$fleet)), 2)
+  retained <- prepare_comp_data(stockplotr::bigskate_lencomp, fleet = "Fishery_current (retained)")
+  expect_equal(unique(retained$fleet), "Fishery_current (retained)")
+})
+
+test_that("single-sex compositions are normalized separately", {
+  one_sex <- data.frame(
+    year = 2000,
+    sex = rep(c("female", "male"), each = 3),
+    bin = rep(1:3, 2),
+    observed = c(0.2, 0.3, 0.5, 0.1, 0.1, 0.8)
+  )
+  comps <- prepare_comp_data(one_sex)
+  expect_equal(comps$observed, one_sex$observed)
+  # joint (two-sex) compositions are normalized together
+  joint <- one_sex
+  joint$observed <- joint$observed / 2
+  expect_equal(prepare_comp_data(joint)$observed, joint$observed)
+})
+
+test_that("prepare_comp_data normalizes counts", {
+  counts <- data.frame(year = 2000, bin = 1:4, observed = c(10, 20, 30, 40))
+  comps <- prepare_comp_data(counts)
+  expect_equal(comps$observed, c(0.1, 0.2, 0.3, 0.4))
+  expect_equal(unique(comps$sex), "combined")
+  expect_true(all(is.na(comps$expected)))
+})
+
+test_that("prepare_comp_data reshapes stockplotr long format", {
+  long <- data.frame(
+    label = rep(c("age_composition_proportion", "age_composition_expected"), each = 3),
+    year = 2020,
+    fleet = "fleet1",
+    age = rep(1:3, 2),
+    estimate = c(0.2, 0.5, 0.3, 0.25, 0.45, 0.3)
+  )
+  comps <- prepare_comp_data(long)
+  expect_equal(nrow(comps), 3)
+  expect_equal(comps$observed, c(0.2, 0.5, 0.3))
+  expect_equal(comps$expected, c(0.25, 0.45, 0.3))
+  expect_equal(unique(comps$comp_type), "age")
+})
+
+test_that("prepare_comp_data reads FIMS-style observed/expected columns", {
+  fims <- data.frame(
+    module_name = "Fleet",
+    label = rep(c("agecomp_expected", "lengthcomp_expected"), each = 4),
+    fleet = "survey1",
+    year = 1,
+    age = c(1:4, rep(NA, 4)),
+    length = c(rep(NA, 4), seq(10, 40, by = 10)),
+    observed = c(1, 2, 3, 4, 4, 3, 2, 1),
+    expected = c(1.5, 2, 3, 3.5, 4, 3, 2, 1)
+  )
+  comps <- prepare_comp_data(fims)
+  expect_setequal(comps$comp_type, c("age", "length"))
+  len <- prepare_comp_data(fims, comp_type = "length")
+  expect_equal(len$bin, c(10, 20, 30, 40))
+  expect_equal(sum(len$observed), 1)
+})
+
+test_that("prepare_comp_data reads r4ss output", {
+  comps <- prepare_comp_data(r4ss_output)
+  expect_setequal(comps$comp_type, c("age", "length"))
+  expect_true("FISHERY" %in% comps$fleet)
+  expect_setequal(unique(comps$sex), c("female", "male"))
+  # r4ss Pearson residuals are kept rather than recalculated
+  len <- prepare_comp_data(r4ss_output$lendbase)
+  expect_equal(
+    sort(len$pearson),
+    sort(as.numeric(r4ss_output$lendbase$Pearson))
+  )
+})
+
+# plot_comp ------------------------------------------------------------------
+test_that("plot_comp makes observation-only and fit plots", {
+  p_obs <- plot_comp(stockplotr::bigskate_lencomp, fleet = "WCGBTS", show_expected = FALSE)
+  expect_s3_class(p_obs, "ggplot")
+  expect_true(has_layer(p_obs, "GeomCol"))
+  expect_false(has_layer(p_obs, "GeomLine"))
+
+  p_fit <- plot_comp(stockplotr::bigskate_lencomp, fleet = "Triennial")
+  expect_true(has_layer(p_fit, "GeomLine"))
+  # males are mirrored below zero
+  expect_true(any(p_fit$data$obs_plot < 0))
+  # r4ss-style sample size labels
+  text_layer <- p_fit$layers[vapply(p_fit$layers, function(l) inherits(l$geom, "GeomText"), logical(1))][[1]]
+  expect_true(all(grepl("^N adj\\.=.*N eff\\.=", text_layer$data$label)))
+})
+
+test_that("plot_comp aggregates across years", {
+  p <- plot_comp(stockplotr::bigskate_lencomp, aggregate = TRUE)
+  expect_s3_class(p, "ggplot")
+  expect_false("year" %in% names(p$data))
+  sums <- p$data |>
+    dplyr::summarise(total = sum(observed), .by = fleet)
+  expect_true(all(abs(sums$total - 1) < 1e-8))
+})
+
+test_that("plot_comp works with r4ss output", {
+  expect_no_error(plot_comp(r4ss_output, comp_type = "length", fleet = "FISHERY"))
+})
+
+# plot_comp_bubble / plot_comp_residuals --------------------------------------
+test_that("bubble plots show observations and residuals", {
+  p <- plot_comp_bubble(stockplotr::bigskate_lencomp)
+  expect_s3_class(p, "ggplot")
+  expect_true(has_layer(p, "GeomPoint"))
+
+  r <- plot_comp_residuals(stockplotr::bigskate_lencomp)
+  expect_s3_class(r, "ggplot")
+  expect_true(all(c("Observed > expected", "Observed < expected") %in% levels(r$data$sign)))
+  expect_true(any(grepl("max =", r$data$panel)))
+})
+
+test_that("residual plot errors without expected values or sample sizes", {
+  obs_only <- dplyr::select(stockplotr::bigskate_lencomp, -expected, -n, -pearson)
+  expect_error(plot_comp_residuals(obs_only), "No Pearson residuals")
+})
+
+# plot_comp_sex_ratio ---------------------------------------------------------
+test_that("plot_comp_sex_ratio calculates ratios and intervals", {
+  p <- plot_comp_sex_ratio(stockplotr::bigskate_lencomp, fleet = "WCGBTS")
+  expect_s3_class(p, "ggplot")
+  expect_true(has_layer(p, "GeomLinerange"))
+  ratios <- p$data
+  ok <- !is.na(ratios$ratio_obs)
+  expect_true(all(ratios$ratio_obs[ok] >= 0 & ratios$ratio_obs[ok] <= 1))
+  expect_true(all(ratios$lower[ok] <= ratios$ratio_obs[ok] + 1e-12))
+  expect_true(all(ratios$upper[ok] >= ratios$ratio_obs[ok] - 1e-12))
+
+  # bins without fish (only the constant Stock Synthesis adds) are omitted
+  # rather than shown as 50:50
+  expect_false(any(ratios$ratio_obs[ok] == 0.5 & ratios$bin[ok] > 180))
+
+  expect_no_error(plot_comp_sex_ratio(stockplotr::bigskate_lencomp, aggregate = TRUE))
+})
+
+test_that("plot_comp_sex_ratio requires two sexes", {
+  combined <- dplyr::filter(stockplotr::bigskate_lencomp, sex == "female")
+  expect_error(plot_comp_sex_ratio(combined), "both")
+})
+
+# plot_comp_francis --------------------------------------------------------------
+test_that("plot_comp_francis matches r4ss Francis weights", {
+  # Values from r4ss::SSMethod.TA1.8() for the simple_small model
+  means <- comp_mean_table(
+    prepare_comp_data(r4ss_output, comp_type = "length", fleet = "FISHERY")
+  )
+  expect_equal(1 / stats::var(means$std_resid), 1.3108142, tolerance = 1e-6)
+  means <- comp_mean_table(
+    prepare_comp_data(r4ss_output, comp_type = "age", fleet = "SURVEY1")
+  )
+  expect_equal(1 / stats::var(means$std_resid), 2.667424, tolerance = 1e-6)
+  # Big Skate: partitions of a fleet are combined, as in r4ss
+  means <- comp_mean_table(
+    prepare_comp_data(stockplotr::bigskate_lencomp, fleet = "Fishery_current")
+  )
+  expect_equal(1 / stats::var(means$std_resid), 1.0019480, tolerance = 1e-6)
+  expect_equal(unique(means$fleet), "Fishery_current (discard & retained)")
+  means <- comp_mean_table(
+    prepare_comp_data(stockplotr::bigskate_lencomp, fleet = "WCGBTS")
+  )
+  expect_equal(1 / stats::var(means$std_resid), 1.0000112, tolerance = 1e-6)
+})
+
+test_that("plot_comp_francis plots observed and expected means", {
+  p <- plot_comp_francis(stockplotr::bigskate_lencomp)
+  expect_s3_class(p, "ggplot")
+  expect_true(has_layer(p, "GeomLine"))
+  expect_true(has_layer(p, "GeomErrorbar"))
+  expect_true(all(grepl("TA1.8 weight", p$data$fleet)))
+
+  no_francis <- plot_comp_francis(stockplotr::bigskate_lencomp, show_francis = FALSE)
+  expect_false(has_layer(no_francis, "GeomErrorbar"))
+})
diff --git a/tests/testthat/test-plot_comp_diagnostics.R b/tests/testthat/test-plot_comp_diagnostics.R
new file mode 100644
index 00000000..f0429789
--- /dev/null
+++ b/tests/testthat/test-plot_comp_diagnostics.R
@@ -0,0 +1,182 @@
+r4ss_output <- readRDS(
+  test_path("fixtures", "r4ss_output", "simple_small", "simple_r4ss.rds")
+)
+
+has_layer <- function(p, geom_class) {
+  any(vapply(p$layers, function(l) inherits(l$geom, geom_class), logical(1)))
+}
+
+# plot_comp_sample_size -------------------------------------------------------
+test_that("effective sample sizes match Stock Synthesis", {
+  comps <- prepare_comp_data(r4ss_output, comp_type = "length")
+  effn <- comp_effn_table(comps)
+  ss3 <- r4ss_output$lendbase |>
+    dplyr::distinct(Yr, Fleet, effN) |>
+    dplyr::mutate(fleet = r4ss_output$FleetNames[Fleet])
+  joined <- dplyr::inner_join(effn, ss3, by = c("year" = "Yr", "fleet"))
+  expect_gt(nrow(joined), 0)
+  expect_equal(joined$effn, joined$effN, tolerance = 1e-3)
+})
+
+test_that("plot_comp_sample_size works with composition data", {
+  p <- plot_comp_sample_size(stockplotr::bigskate_lencomp)
+  expect_s3_class(p, "ggplot")
+  expect_true(has_layer(p, "GeomAbline"))
+  expect_true(all(grepl("McAllister-Ianelli weight", p$data$panel)))
+  # Matches Stock Synthesis: Recommend_var_adj = Curr_Var_Adj * weight, with
+  # partitions of a fleet combined (Length_Comp_Fit_Summary for Big Skate)
+  weights <- dplyr::distinct(p$data, base_fleet, mi_weight)
+  ss3 <- data.frame(
+    base_fleet = c("Fishery_current", "WCGBTS", "Triennial"),
+    curr_var_adj = c(0.239885, 0.067420, 1),
+    recommend_var_adj = c(0.107271, 0.637793, 0.870899)
+  )
+  joined <- dplyr::inner_join(weights, ss3, by = "base_fleet")
+  expect_equal(nrow(joined), 3)
+  # (tolerance allows for the precision of values printed by Stock Synthesis)
+  expect_equal(joined$mi_weight * joined$curr_var_adj, joined$recommend_var_adj, tolerance = 1e-3)
+  expect_no_error(plot_comp_sample_size(r4ss_output, comp_type = "age", show_years = TRUE))
+})
+
+test_that("plot_comp_sample_size works with converted SS3 output", {
+  p <- plot_comp_sample_size(stockplotr::example_data, comp_type = "length")
+  expect_s3_class(p, "ggplot")
+  expect_true(all(p$data$comp_type == "length"))
+  # retained and discard compositions are separate panels
+  expect_true(any(grepl("retained", p$data$fleet)))
+  expect_true(any(grepl("discard", p$data$fleet)))
+  # conditional age-at-length comps are dropped from age comps
+  expect_message(
+    plot_comp_sample_size(stockplotr::example_data, comp_type = "age"),
+    "conditional age-at-length"
+  )
+})
+
+test_that("plot_comp_sample_size errors without expected values", {
+  obs_only <- dplyr::select(stockplotr::bigskate_lencomp, -expected)
+  expect_error(plot_comp_sample_size(obs_only), "could not be calculated")
+})
+
+# residual distribution and patterns -----------------------------------------
+test_that("plot_comp_residual_dist makes qq plots and histograms", {
+  qq <- plot_comp_residual_dist(stockplotr::bigskate_lencomp)
+  expect_s3_class(qq, "ggplot")
+  expect_true(has_layer(qq, "GeomPoint"))
+  expect_true(all(grepl("SDNR", qq$data$panel)))
+
+  hist <- plot_comp_residual_dist(stockplotr::bigskate_lencomp, type = "histogram")
+  expect_true(has_layer(hist, "GeomBar"))
+})
+
+test_that("plot_comp_residual_pattern summarizes by bin and year", {
+  p_bin <- plot_comp_residual_pattern(stockplotr::bigskate_lencomp)
+  expect_s3_class(p_bin, "ggplot")
+  expect_true("bin" %in% names(p_bin$data))
+  # mean residual by bin matches a direct calculation
+  direct <- stockplotr::bigskate_lencomp |>
+    dplyr::filter(fleet == "WCGBTS", sex == "female", bin == 50)
+  check <- dplyr::filter(p_bin$data, fleet == "WCGBTS", sex == "female", bin == 50)
+  expect_equal(check$mean_resid, mean(direct$pearson))
+
+  p_year <- plot_comp_residual_pattern(stockplotr::bigskate_lencomp, by = "year")
+  expect_true("year" %in% names(p_year$data))
+})
+
+test_that("residual plots error without residuals", {
+  obs_only <- dplyr::select(stockplotr::bigskate_lencomp, -expected, -n, -pearson)
+  expect_error(plot_comp_residual_dist(obs_only), "No Pearson residuals")
+  expect_error(plot_comp_residual_pattern(obs_only), "No Pearson residuals")
+})
+
+test_that("runs_test detects patterns", {
+  expect_lt(runs_test(c(rep(1, 10), rep(-1, 10))), 0.05)
+  expect_gt(runs_test(rep(c(1, -1), 10)), 0.5)
+  expect_equal(runs_test(rep(-1, 8)), 2 * 0.5^8)
+  expect_true(is.na(runs_test(c(1, -1))))
+})
+
+# conditional age-at-length --------------------------------------------------
+test_that("prepare_comp_data handles conditional age-at-length data", {
+  caal <- prepare_comp_data(stockplotr::bigskate_caal, conditional = TRUE)
+  expect_true("length_bin" %in% names(caal))
+  # single-sex compositions stay normalized within each sex and length bin
+  sums <- caal |>
+    dplyr::summarise(total = sum(observed), .by = c(fleet, year, sex, length_bin))
+  expect_true(all(abs(sums$total - 1) < 1e-6))
+  expect_equal(caal$observed, stockplotr::bigskate_caal$observed[
+    order(stockplotr::bigskate_caal$fleet, stockplotr::bigskate_caal$year,
+      stockplotr::bigskate_caal$length_bin, stockplotr::bigskate_caal$sex,
+      stockplotr::bigskate_caal$bin)
+  ], tolerance = 1e-6)
+  # marginal plots refuse conditional data
+  expect_error(prepare_comp_data(stockplotr::bigskate_caal), "conditional = TRUE")
+  expect_error(plot_comp(stockplotr::bigskate_caal), "conditional = TRUE")
+  # and conditional plots need a length bin
+  expect_error(
+    prepare_comp_data(stockplotr::bigskate_lencomp, conditional = TRUE),
+    "length_bin"
+  )
+})
+
+test_that("prepare_comp_data reads conditional data from r4ss dbase", {
+  dbase <- r4ss_output$agedbase
+  # make half of the age comps conditional on a narrow length bin
+  dbase$Lbin_hi[seq_len(nrow(dbase) / 2)] <- dbase$Lbin_lo[seq_len(nrow(dbase) / 2)]
+  caal <- prepare_comp_data(dbase, conditional = TRUE)
+  expect_equal(nrow(caal), nrow(dbase) / 2)
+  marginal <- suppressMessages(prepare_comp_data(dbase))
+  expect_equal(nrow(marginal), nrow(dbase) / 2)
+  # r4ss output without conditional data
+  expect_error(prepare_comp_data(r4ss_output, conditional = TRUE), "No conditional")
+})
+
+test_that("plot_comp_andre shows mean and SE of age at length", {
+  dat <- dplyr::filter(stockplotr::bigskate_caal, year %in% 2017:2018)
+  p <- plot_comp_andre(dat, fleet = "WCGBTS")
+  expect_s3_class(p, "ggplot")
+  expect_setequal(as.character(unique(p$data$statistic)), c("Mean age", "SE of mean age"))
+  expect_true(has_layer(p, "GeomRibbon"))
+  # observed mean ages increase with length
+  means <- dplyr::filter(p$data, statistic == "Mean age", year == 2018, sex == "female")
+  expect_gt(stats::cor(means$length_bin, means$observed), 0.8)
+  # one value per year, sex, and length bin, calculated as in r4ss
+  z <- dplyr::filter(
+    stockplotr::bigskate_caal,
+    fleet == "WCGBTS", year == 2018, sex == "female", length_bin == 60
+  )
+  check <- dplyr::filter(p$data, year == 2018, sex == "female", length_bin == 60)
+  obs_mean <- sum(z$bin * z$observed) / sum(z$observed)
+  expect_equal(check$observed[check$statistic == "Mean age"], obs_mean)
+  expect_equal(
+    check$expected[check$statistic == "SE of mean age"],
+    sqrt(sum(z$bin^2 * z$expected) / sum(z$expected) -
+      (sum(z$bin * z$expected) / sum(z$expected))^2) / sqrt(z$n[1])
+  )
+})
+
+test_that("conditional bubble plots use length bins on the x-axis", {
+  dat <- dplyr::filter(stockplotr::bigskate_caal, year %in% 2017:2018)
+  p <- plot_comp_residuals(dat, fleet = "WCGBTS", conditional = TRUE)
+  expect_s3_class(p, "ggplot")
+  expect_true(any(grepl("^female, 2017\n\\(max", levels(p$data$panel))))
+  expect_no_error(plot_comp_bubble(dat, fleet = "WCGBTS", conditional = TRUE))
+})
+
+test_that("r4ss sex codes are read correctly", {
+  # single-sex male compositions have the Stock Synthesis code Sex = 1 but
+  # r4ss sex = 2, and appear twice in the composition database
+  dbase <- data.frame(
+    Yr = 2000, Fleet = 1, Kind = "AGE", Part = 0,
+    Sexes = 2, Sex = c(1, 1, 2, 2), sex = 2,
+    Lbin_lo = 50, Lbin_hi = 50, Bin = c(1, 2, 1, 2),
+    Obs = c(0.4, 0.6, 0.4, 0.6), Exp = c(0.5, 0.5, 0.5, 0.5),
+    Pearson = c(-0.2, 0.2, -0.2, 0.2), Nsamp_adj = 10
+  )
+  comps <- prepare_comp_data(dbase, conditional = TRUE)
+  expect_equal(nrow(comps), 2)
+  expect_equal(unique(comps$sex), "male")
+  expect_equal(comps$observed, c(0.4, 0.6))
+  # without the r4ss sex column, sex is derived from Sexes and Sex
+  comps <- prepare_comp_data(dplyr::select(dbase, -sex), conditional = TRUE)
+  expect_equal(unique(comps$sex), "male")
+})
diff --git a/vignettes/composition_diagnostics.Rmd b/vignettes/composition_diagnostics.Rmd
new file mode 100644
index 00000000..bed02976
--- /dev/null
+++ b/vignettes/composition_diagnostics.Rmd
@@ -0,0 +1,268 @@
+---
+title: "Composition data diagnostics"
+output:
+  html_document:
+    toc: true
+vignette: >
+  %\VignetteIndexEntry{Composition data diagnostics}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>",
+  fig.width = 9,
+  fig.height = 7,
+  out.width = "100%"
+)
+```
+
+```{r setup}
+library(stockplotr)
+```
+
+Age and length compositions are often the most informative, and the hardest
+to fit, data in a stock assessment. {stockplotr} has a set of plots for
+composition data based on the figures that
+[r4ss](https://github.com/r4ss/r4ss) makes for Stock Synthesis models. Each
+plot shows the data in a different way, and each has an r4ss equivalent:
+
+| Plot | Function | What it shows | r4ss equivalent |
+|------|----------|---------------|-----------------|
+| Distributions of observations | `plot_comp(show_expected = FALSE)` | The most intuitive way to look at compositions | `SSplotComps()` subplot 1, `datonly = TRUE` |
+| Bubble plot of observations | `plot_comp_bubble()` | Cohort strength over time | `SSplotComps()` subplot 2, `datonly = TRUE` |
+| Distributions of observed and expected values | `plot_comp()` | An intuitive way to look at fit | `SSplotComps()` subplots 1 and 21 |
+| Pearson residuals | `plot_comp_residuals()` | Blocks or stripes of misfit among years or bins | `SSplotComps()` subplots 2 and 24 |
+| Sex ratios | `plot_comp_sex_ratio()` | Trends in misfit for bins with small proportions, which have small residuals | `SSplotSexRatio()` |
+| Mean age or length (Francis) | `plot_comp_francis()` | Trends among years; Francis (2011) TA1.8 data weights | `SSplotComps()` subplots 8 and 9, `SSMethod.TA1.8()` |
+| N-EffN comparison | `plot_comp_sample_size()` | McAllister-Ianelli data weights | `SSplotComps()` subplot 7 |
+| Residual distribution | `plot_comp_residual_dist()` | Heavy tails, skewness, and outliers; SDNR | |
+| Residual patterns | `plot_comp_residual_pattern()` | Bias by bin or year, with a runs test | |
+| Andre's plots | `plot_comp_andre()` | Mean age at length for conditional age-at-length data | `SSplotComps()` subplot 12 |
+| Conditional age-at-length residuals | `plot_comp_residuals(conditional = TRUE)` | Misfit by length and age within years | `SSplotComps()` subplot 3 |
+
+Each function can plot fleets separately or summarize all fleets. Use `fleet`
+to choose fleets, and `aggregate = TRUE` in `plot_comp()` and
+`plot_comp_sex_ratio()` to combine all years.
+
+## Composition data
+
+The plotting functions take composition data in one of several formats and
+standardize it with `prepare_comp_data()`:
+
+* A tidy data frame with columns `year`, `bin`, and `observed`, and,
+  optionally, `fleet`, `sex`, `part` (partition), `expected`, `n` (sample
+  size), `pearson`, and `comp_type` ("age" or "length").
+* FIMS output with `observed` and `expected` columns, such as the output of
+  `FIMS::get_estimates()`.
+* {stockplotr} standard output from `convert_output()` with labels such as
+  `age_composition_proportion` and `age_composition_expected`.
+* Output from `r4ss::SS_output()`, either the whole list or its `lendbase`,
+  `agedbase`, or `condbase` data frames.
+
+This vignette uses data from the 2019 Big Skate stock assessment for the U.S.
+West Coast, the same model used for the
+[r4ss example plots](https://r4ss.github.io/r4ss/BigSkate/), so the figures
+below can be compared directly with the r4ss versions. `bigskate_lencomp`
+holds the length compositions (`lendbase` from `r4ss::SS_output()`) and
+`bigskate_caal` the conditional age-at-length data (`condbase`), with columns
+renamed to follow stockplotr conventions.
+
+```{r}
+head(bigskate_lencomp)
+```
+
+The length compositions are sex-specific. `Fishery_current` has separate
+retained (`part` = 2) and discarded (`part` = 1) compositions, which are
+labeled "Fishery_current (retained)" and "Fishery_current (discard)", as in
+r4ss. Use `fleet = "Fishery_current"` to select both.
+
+If Pearson residuals are not supplied, `prepare_comp_data()` calculates them
+assuming a multinomial distribution,
+$r = (o - e) / \sqrt{e (1 - e) / N}$, using the sample size `n`. The Big
+Skate data include the residuals reported by Stock Synthesis.
+
+## Observations
+
+Distributions are the most intuitive way to look at compositions. For
+sex-specific data, females are shown above zero and males below.
+
+```{r obs-dist, fig.height = 9}
+plot_comp(bigskate_lencomp, fleet = "WCGBTS", show_expected = FALSE)
+```
+
+Bubble plots make it easier to see patterns across years, such as strong
+cohorts moving through the population.
+
+```{r obs-bubble}
+plot_comp_bubble(bigskate_lencomp)
+```
+
+## Fits to the data
+
+Adding the expected values to the distribution plot shows the fit to each
+year. As in r4ss, each panel shows the adjusted input sample size (N adj.)
+and the effective sample size (N eff.).
+
+```{r fit-dist, fig.height = 4}
+plot_comp(bigskate_lencomp, fleet = "Triennial")
+```
+
+Aggregating across years gives an overview of each fleet. Years are weighted
+by their adjusted input sample sizes.
+
+```{r fit-agg, fig.height = 6}
+plot_comp(bigskate_lencomp, aggregate = TRUE)
+```
+
+### Pearson residuals
+
+Pearson residuals are a more statistically rigorous way to look at fit.
+Filled bubbles show where more fish were observed than expected, and open
+bubbles where fewer were observed. The largest residual in each panel is
+included in the label.
+
+```{r resid, fig.height = 8}
+plot_comp_residuals(bigskate_lencomp)
+```
+
+### Sex ratios
+
+The fraction female in each bin, female / (female + male), can show trends
+in misfit for bins with small proportions, since those bins have small
+residuals. Vertical lines are Jeffreys intervals based on the sample size in
+each bin. For Big Skate, the largest fish are almost all female.
+
+```{r sex-ratio, fig.height = 9}
+plot_comp_sex_ratio(bigskate_lencomp, fleet = "WCGBTS")
+```
+
+```{r sex-ratio-agg, fig.height = 6}
+plot_comp_sex_ratio(bigskate_lencomp, aggregate = TRUE)
+```
+
+### Mean length and Francis weights
+
+The time series of mean age or length was developed as a diagnostic for the
+Francis (2011) data weighting method (TA1.8), but it is also a useful summary
+of trends in the data and the fit. Thick intervals are based on the current
+sample sizes, and thin, capped intervals show the result of further adjusting
+sample sizes by the TA1.8 weight in each panel label. As in
+`r4ss::SSMethod.TA1.8()`, partitions of a fleet are combined into one
+composition per year. The weights near one show that this model was already
+tuned with the Francis method.
+
+```{r francis, fig.height = 5}
+set.seed(1) # the interval for the Francis weight is bootstrapped
+plot_comp_francis(bigskate_lencomp)
+```
+
+### N-EffN comparison
+
+The effective sample size measures how well the model fits each composition
+(McAllister and Ianelli 1997). If the input sample sizes suit the fit, the
+points should scatter around the 1:1 line. The ratio of the harmonic mean
+effective sample size to the mean input sample size, shown in each panel
+label, is the McAllister-Ianelli weight. As in Stock Synthesis, it is
+calculated across all partitions of a fleet, and multiplying it by the
+current variance adjustment gives the `Recommend_var_adj` that Stock
+Synthesis reports.
+
+```{r sample-size, fig.height = 6}
+plot_comp_sample_size(bigskate_lencomp)
+```
+
+For Stock Synthesis models, the same plot can be made directly from
+`convert_output()` output, using the adjusted input and effective sample
+sizes that Stock Synthesis reports:
+
+```{r sample-size-ss3}
+plot_comp_sample_size(stockplotr::example_data, comp_type = "length")
+```
+
+### Residual distribution and patterns
+
+A quantile-quantile plot shows whether residuals have heavier tails than
+expected. The standard deviation of the normalized residuals (SDNR) is shown
+in each panel label; values well above one suggest the data are given more
+weight than the fit supports.
+
+```{r resid-qq, fig.height = 6}
+plot_comp_residual_dist(bigskate_lencomp)
+```
+
+Averaging residuals by bin or year turns stripes in the bubble plot into a
+single line. Filled points fail a runs test, meaning the residual signs in
+that bin are not randomly ordered across years.
+
+```{r resid-pattern, fig.height = 8}
+plot_comp_residual_pattern(bigskate_lencomp)
+```
+
+## Conditional age-at-length data
+
+Conditional age-at-length data are age compositions for each length bin.
+Use `conditional = TRUE` in `prepare_comp_data()` and `plot_comp_bubble()`,
+and supply a `length_bin` column in tidy data, as in `bigskate_caal`. The Big
+Skate conditional data are single-sex compositions, which are normalized
+separately for each sex.
+
+`plot_comp_andre()` makes "Andre's plots" from r4ss: the observed and
+expected mean age in each length bin (left), and the standard error of the
+mean age (right), one row per year and sex.
+
+```{r caal-andre, fig.height = 10}
+caal_recent <- dplyr::filter(bigskate_caal, year %in% 2016:2018)
+plot_comp_andre(caal_recent, fleet = "WCGBTS")
+```
+
+Pearson residuals for conditional data are plotted by length bin and age for
+each year:
+
+```{r caal-resid, fig.height = 6}
+plot_comp_residuals(caal_recent, fleet = "WCGBTS", conditional = TRUE)
+```
+
+## Using model output
+
+With FIMS, pass the estimates directly:
+
+```{r fims, eval = FALSE}
+estimates <- FIMS::get_estimates(fit)
+plot_comp_residuals(estimates, comp_type = "age")
+plot_comp_francis(estimates, comp_type = "length")
+```
+
+FIMS reports years as indices (`year_i`); add a `year` column to plot
+calendar years. FIMS estimates do not currently include composition sample
+sizes. Add an `n` column, or supply Pearson residuals in a `pearson` column,
+to get residuals and intervals.
+
+With Stock Synthesis, use the output of `r4ss::SS_output()`. Pearson residuals
+and adjusted sample sizes are taken from the model output, the Francis
+weights match `r4ss::SSMethod.TA1.8()`, and the effective sample sizes and
+McAllister-Ianelli weights match those reported by Stock Synthesis.
+Conditional age-at-length data come from the `condbase` element:
+
+```{r ss3, eval = FALSE}
+replist <- r4ss::SS_output("path/to/BigSkate_2019")
+plot_comp(replist, comp_type = "length", fleet = "WCGBTS")
+plot_comp_francis(replist, comp_type = "length")
+plot_comp_andre(replist, fleet = "WCGBTS")
+```
+
+## References
+
+Francis, R.I.C.C. (2011). Data weighting in statistical fisheries stock
+assessment models. Canadian Journal of Fisheries and Aquatic Sciences 68:
+1124-1138.
+
+McAllister, M.K., and Ianelli, J.N. (1997). Bayesian stock assessment using
+catch-age data and the sampling-importance resampling algorithm. Canadian
+Journal of Fisheries and Aquatic Sciences 54: 284-300.
+
+Taylor, I.G., Gertseva, V., Stephens, A., and Bizzarro, J. (2019). Status of
+Big Skate (*Beringraja binoculata*) off the U.S. Pacific Coast in 2019.
+Pacific Fishery Management Council, Portland, OR.