Indicator Exploratory Analysis: Doctor Visits: Smoothed Adj CLI

This notebook evaluates a candidate indicator to understand its characteristics, coverage, missingness, and versioning/revision behavior.

---

1 Data Loading

# Load candidate indicator
candidate_archive <- NULL
if (!is.null(params$input_csv)) {
    df_candidate_raw <- read.csv(params$input_csv) %>%
        mutate(time_value = as.Date(time_value)) %>%
        filter(time_value >= as.Date(params$start_day) & time_value <= as.Date(params$end_day))

    # Add signal column if missing (required for long signal format)
    if (!("signal" %in% colnames(df_candidate_raw))) {
        df_candidate_raw$signal <- if (length(signals_vec) > 0) signals_vec[1] else "value"
    }

    # CSV may include version/issue columns for revision analysis
    has_versions <- any(c("version", "issue") %in% colnames(df_candidate_raw))
    if (has_versions) {
        df_candidate_raw <- df_candidate_raw %>%
            {
                if ("issue" %in% names(.)) rename(., version = issue) else .
            } %>%
            mutate(version = as.Date(version))
        csv_archive_df  <- df_candidate_raw %>%
            select(geo_value, time_value, version, any_of("signal"), value) %>%
            group_by(geo_value, time_value, version, pick(any_of("signal"))) %>%
            slice_max(value, n = 1, with_ties = FALSE) %>%
            ungroup()

        candidate_archive <- csv_archive_df %>%
            as_epi_archive(signal_format = "long")
        df_candidate <- epix_as_of(candidate_archive, candidate_archive$versions_end)
        cat(sprintf(
            "- CSV contains version history: %s versions\n",
            label_comma()(n_distinct(df_candidate_raw$version))
        ))
    } else {
        df_candidate <- df_candidate_raw %>%
            as_epi_df(signal_format = "long")
    }
} else {
    # Determine if new API can be used for candidate source
    new_api_sources <- c("nssp", "nhsn", "pophive", "nwss", "claims_data_outpatient")
    use_new_api <- FALSE
    if (params$source %in% new_api_sources) {
        meta_available <- tryCatch(
            {
                meta_res <- epidata_meta(params$source)
                all(signals_vec %in% meta_res[[params$source]]$signals)
            },
            error = function(e) FALSE
        )
        if (meta_available) use_new_api <- TRUE
    }
    if (use_new_api) {
        # epidata_snapshot does not accept multiple signals
        df_candidate <- purrr::map(signals_vec, ~ {
            epidata_snapshot(
                source = params$source,
                signals = .x,
                geo_type = params$geo_type
            ) %>%
                rename(any_of(c(time_value = "reference_time", version = "report_time"))) %>%
                {
                    if ("fill_method" %in% colnames(.)) {
                        filter(., fill_method == "source" | !any(fill_method == "source"))
                    } else .
                } %>%
                # {
                #     if ("version" %in% colnames(.)) {
                #         group_by(., geo_value, time_value) %>%
                #             slice_max(version, n = 1, with_ties = FALSE) %>%
                #             ungroup()
                #     } else .
                # } %>%
                filter(
                    time_value >= as.Date(params$start_day),
                    time_value <= as.Date(params$end_day)
                ) %>%
                select(geo_value, any_of("signal"), time_value, value)
        }) %>%
            bind_rows() %>%
            as_epi_df(signal_format = "long")
    } else {
        # pub_covidcast accepts a signal vector and returns a signal column
        df_candidate <- pub_covidcast(
            params$source, signals_vec,
            params$geo_type, params$time_type,
            time_values = time_range(params$start_day, params$end_day),
            as_of = Sys.Date()
        ) %>% as_epi_df(signal_format = "long")
    }
}

cat(sprintf(
    "- **Summary**: %s rows across %s locations and %d signal(s)\n",
    label_comma()(nrow(df_candidate)),
    label_comma()(n_distinct(df_candidate$geo_value)),
    length(signals_vec)
))

• Summary: 67,798 rows across 2,330 locations and 1 signal(s)

# Fetch revision history from API if not already built from CSV
if (is.null(candidate_archive) && is.null(params$input_csv)) {
    # Only fetch version history if the source supports revisions
    new_api_sources <- c("nssp", "nhsn", "pophive", "nwss", "claims_data_outpatient")
    use_new_api <- FALSE
    if (params$source %in% new_api_sources) {
        meta_available <- tryCatch(
            {
                meta_res <- epidata_meta(params$source)
                all(signals_vec %in% meta_res[[params$source]]$signals)
            },
            error = function(e) FALSE
        )
        if (meta_available) use_new_api <- TRUE
    }

    # Fetch archive for revision behavior
    archive_geos <- unique(df_candidate$geo_value)
    # Limit number of locations for version history queries to avoid API limits
    if (length(archive_geos) > params$max_locations_plot) {
        archive_geos <- random_sample(archive_geos, params$max_locations_plot)
    }

    # Determine window of issues available
    archive_start <- as.Date(params$start_day) - 60

    archive_info <- list(
        n_sampled = length(archive_geos),
        n_total = length(unique(df_candidate$geo_value)),
        time_value_start = archive_start,
        time_value_end = params$end_day
    )

    if (use_new_api) {
        # epidata_archive does not accept multiple signals
        df_archive_raw <- purrr::map(signals_vec, ~ {
            epidata_archive(
                source = params$source,
                signals = .x,
                geo_type = params$geo_type
            ) %>%
                rename(any_of(c(time_value = "reference_time", version = "report_time"))) %>%
                filter(
                    geo_value %in% archive_geos,
                    time_value >= as.Date(archive_start),
                    time_value <= as.Date(params$end_day)
                )
        }) %>% bind_rows()
    } else {
        # pub_covidcast accepts a signal vector and returns a signal column
        df_archive_raw <- pub_covidcast(
            source = params$source,
            signal = signals_vec,
            geo_type = params$geo_type,
            time_type = params$time_type,
            geo_values = archive_geos,
            time_values = time_range(archive_start, params$end_day),
            issues = "*",
        )
    }

    if (nrow(df_archive_raw) > 0) {
        df_archive_raw <- df_archive_raw %>%
            rename(any_of(c(time_value = "reference_time", version = "issue"))) %>%
            rename(any_of(c(version = "report_time"))) %>%
            select(geo_value, time_value, version, any_of("signal"), value) #%>%
            # Deduplicate: sources like NWSS can have multiple rows per
            # (geo, time, version) from different measurement sources.
            # group_by(geo_value, time_value, version, across(any_of("signal"))) %>%
            # slice_max(value, n = 1, with_ties = FALSE) %>%
            # ungroup()

        if ("version" %in% colnames(df_archive_raw)) {
            archive_info$version_start <- min(df_archive_raw$version)
            archive_info$version_end <- max(df_archive_raw$version)
            candidate_archive <- df_archive_raw %>%
                as_epi_archive(signal_format = "long")
        }
    }
}

# Determine time type and number of locations
env_time_type <- attributes(df_candidate)$metadata$time_type
if (is.null(env_time_type) || length(env_time_type) != 1 || is.na(env_time_type) || !is.character(env_time_type)) {
    env_time_type <- params$time_type
}
if (is.null(env_time_type) || !nzchar(env_time_type)) {
    env_time_type <- "day"
}
params$time_type <- env_time_type

env_geo_type <- attributes(df_candidate)$metadata$geo_type
if (is.null(params$geo_type) || length(params$geo_type) != 1 || is.na(params$geo_type) || !is.character(params$geo_type)) {
    params$geo_type <- env_geo_type
}

n_locations <- n_distinct(df_candidate$geo_value)

2 Exploratory Data Analysis

The data is aggregated at the county level. The dataset contains 2330 valid locations measured at a day-level frequency.

2.1 Time Series Overview

This section provides a visual overview of the candidate indicator over time for a subset of observation locations. We look at the raw values for the indicator across a sample of locations.

# Pick a diverse sample of locations for visualization
all_geos <- sort(unique(df_candidate$geo_value))
n_plot <- min(length(all_geos), params$max_locations_plot)
sample_geos <- random_sample(all_geos, n_plot)

df_sample <- df_candidate %>%
    filter(geo_value %in% sample_geos)

autoplot(df_sample, value,
        .interactive = TRUE, .facet_by = "geo_value", .facet_to_dropdown = TRUE) %>%
        plotly::toWebGL() %>%
        plotly::partial_bundle()

2.2 Quantile Trends

This plot displays the distribution of both indicators over time across all included locations. The solid line represents the median, the inner ribbon shows the 25th to 75th percentile range, and the outer ribbon covers the 5th to 95th percentile range. We smooth the data with a moving average to handle day-of-week effects, and we use different scales for the top and bottom panels so each indicator’s distribution is clearly visible.

By observing the width of these ribbons, you can get a sense of geographic variability. Narrow ribbons mean the values are tightly clustered across locations, while wide ribbons show more spread. This helps us assess the overall stability and systemic trends of the indicators.

cat("\n")

# Aggregate location values per (signal, time_value); use signal name as geo_value
# so epi_slide partitions by signal automatically.
df_quantiles_wide <- df_candidate %>%
    group_by(signal, time_value) %>%
    summarize(
        vals = list(value), 
        geo_value = "all",
        .groups = "drop"
    ) %>%
    as_epi_df(signal_format = "long")

# Handle time frequency differences for .window_size.
window_size <- switch(env_time_type %||% "unknown",
    day   = 7L, # 7-day smoothing
    week  = as.difftime(3, units = "weeks"), # 3-week smoothing
    month = as.difftime(3, units = "months"), # 3-month smoothing
    year  = as.difftime(3, units = "years"), # 3-year smoothing
    Inf # Default to cumulative for unknown types
    # (too sparse for sliding averages)
)

df_quantiles_wide <- df_quantiles_wide %>%
    group_by(geo_value, signal) %>%
    epi_slide(
        ~ tibble(
            `0.05` = quantile(unlist(.x$vals), 0.05, na.rm = TRUE),
            `0.25` = quantile(unlist(.x$vals), 0.25, na.rm = TRUE),
            `0.5`  = quantile(unlist(.x$vals), 0.5, na.rm = TRUE),
            `0.75` = quantile(unlist(.x$vals), 0.75, na.rm = TRUE),
            `0.95` = quantile(unlist(.x$vals), 0.95, na.rm = TRUE)
        ),
        .window_size = window_size
    ) %>%
    ungroup() %>%
    mutate(indicator = signal_labels[signal])

ggplot(df_quantiles_wide, aes(x = time_value)) +
    geom_ribbon(aes(ymin = `0.05`, ymax = `0.95`, fill = "5%-95%"), alpha = 0.2) +
    geom_ribbon(aes(ymin = `0.25`, ymax = `0.75`, fill = "25%-75%"), alpha = 0.4) +
    geom_line(aes(y = `0.5`, color = "Median"), linewidth = 1) +
    facet_wrap(~indicator, ncol = 1, scales = "free_y") +
    scale_x_date(breaks = breaks_pretty(), label = label_date_short()) +
    scale_y_continuous(labels = scales::label_number(scale_cut = scales::cut_si(""))) +
    scale_fill_manual(values = c("5%-95%" = dis_pal_list[2], "25%-75%" = dis_pal_list[2])) +
    scale_color_manual(values = c("Median" = dis_pal_list[2])) +
    labs(
        title = "Smoothed Quantile Trends",
        subtitle = sprintf("Across all %d locations, smoothed window: %s", n_locations, format(window_size)),
        x = "Date", y = "Value (smoothed)",
        fill = "Range", color = "Trend"
    ) +
    theme(legend.position = "bottom")

2.2.1 Cross-Sectional Geographic Outliers

While the quantile plots summarize the main body of the distribution, they may hide the behavior of individual highly erratic locations. Here, we track how often individual locations jump into the extreme upper tail (above the cross-sectional 95th percentile recorded on that specific date).

Note

This approach compares locations against one another on each date. If the indicator relies on raw counts rather than population-standardized rates or percentages, this table will frequently just identify locations with the largest populations. In those cases, this metric reflects size rather than outliers.

2.3 Coverage and Missingness

This section explores the availability of data across time and space. The area chart below tracks the number of geographic locations with valid data on any given date for the candidate indicator.

ggplot(coverage, aes(x = time_value, y = n_available,
                     color = indicator, linetype = indicator)) +
    geom_line(linewidth = 0.8) +
    scale_x_date(breaks = breaks_pretty(), label = label_date_short()) +
    scale_y_continuous(labels = label_comma(), limits = c(0, NA)) +
    labs(
        title = "Data coverage over time",
        subtitle = sprintf("Number of locations with non-missing values per %s%s", env_time_type,
                           if (env_time_type != params$time_type && slide_window > 1) sprintf(" (%d-%s rolling average)", slide_window, env_time_type) else ""),
        x = "Date", y = "Number of locations",
        color = "Indicator", linetype = "Indicator"
    ) +
    {if (length(signals_vec) > 6) facet_wrap(~indicator, ncol = 2, labeller = label_wrap_html(width = 25)) else NULL}

2.3.1 Summary Statistics by Location

This table summarizes the available data across all included locations. Locations with no candidate information are not considered.

Zero variance locations (flatlined) — Doctor Visits: Smoothed Adj CLI: 01047 (Dallas County), 05055 (Greene County), 13189 (McDuffie County), 16007 (Bear Lake County), 19125 (Marion County), 20159 (Rice County), 21039 (Carlisle County), 27039 (Dodge County), 29105 (Laclede County), 29171 (Putnam County), 31101 (Keith County), 39141 (Ross County), 42127 (Wayne County), 46033 (Custer County)

Missingness Overview (All Locations):

• <5% missing: 1806 locations
• 5-25% missing: 213 locations
• 25-50% missing: 140 locations
• 50-99.9% missing: 171 locations
• 100% missing: 0 locations

2.4 Noise Assessment

We calculate the Coefficient of Variation (CV) for each location to compare indicator stability across geographic areas.

\[\text{CV}_i = \frac{\text{SD}(X_i)}{|\bar{X}_i|}\]

# Calculate CV per (signal, location)
noise_stats <- df_candidate %>%
    arrange(signal, geo_value, time_value) %>%
    group_by(signal, geo_value) %>%
    summarize(
        cand_cv = sd(value, na.rm = TRUE) / abs(mean(value, na.rm = TRUE)),
        .groups = "drop"
    ) %>%
    filter(is.finite(cand_cv))

# Clip extreme outliers for robust visualization
cv_upper <- quantile(noise_stats$cand_cv, 0.99, na.rm = TRUE)
noise_stats_plot <- noise_stats %>%
    filter(cand_cv <= cv_upper) %>%
    mutate(signal_label = signal_labels[signal])

n_excluded <- nrow(noise_stats) - nrow(noise_stats_plot)

# Summary Table — one row per signal
noise_summary <- noise_stats %>%
    group_by(signal) %>%
    summarize(`Median CV` = median(cand_cv), .groups = "drop") %>%
    mutate(Signal = signal_labels[signal]) %>%
    select(Signal, `Median CV`)

DT::datatable(noise_summary,
    caption = sprintf(
        "CV summary across locations.%s",
        if (n_excluded > 0) sprintf(" (%d extreme outlier(s) excluded from plot.)", n_excluded) else ""
    ),
    rownames = FALSE, options = list(scrollX = TRUE, paging = FALSE, dom = "t")
) %>% DT::formatRound(columns = 2, digits = 3)

# Histogram (faceted by signal when multi_signal)
ggplot(noise_stats_plot, aes(x = cand_cv)) +
    geom_histogram(bins = 30, color = "white", fill = dis_pal_list[2], alpha = 0.7) +
    scale_x_continuous(labels = label_percent()) +
    {if (multi_signal) facet_wrap(~signal_label, scales = "free", labeller = label_wrap_html(width = 25)) else NULL} +
    labs(
        title = "CV Distribution by Location",
        subtitle = "Lower values indicate a more stable indicator",
        x = "Coefficient of Variation (CV)", y = "Count"
    ) +
    theme(legend.position = "none")

2.5 Candidate Versioning

This section examines how published values change over time as new data arrives. Each observation may be revised multiple times after its initial release, and the statistics below summarize both how often and how much those revisions occur. Relative revisions measure the magnitude of a revision relative to the signal’s own scale. A relative revision of 0.1 means the observation shifted by 10% of its own magnitude across revisions, while values near 1 indicate revisions as large as the signal itself.

The revision_summary is under revision and may produce misleading results in certain edge cases. Use the summary statistics below with caution while development is ongoing.

Revision analysis is performed on a sample of 60 out of 2330 locations, with versions from 2022-11-06 to 2023-04-15, and time from 2022-11-02 to 2023-02-01.

• Time range: 2022-11-02 to 2023-02-01
• Min lag (days to first version): Min.: 4, 1st Qu.: 4, Median: 4, Mean: 4.1, 3rd Qu.: 4, Max.: 5
• Total number of versions: 151
• No revisions: 1%
• Quick revisions (<= 3 days): 0%
• Few revisions (<= 3): 2%
• Minor relative revisions (< 0.1): 0%
• Major absolute revisions (> 0.395): 60%
• Days until within 20% of latest: Min.: 4, 1st Qu.: 11, Median: 20, Mean: 25.3, 3rd Qu.: 35, Max.: 73

The time series below visualize the candidate indicator(s) (Doctor Visits: Smoothed Adj CLI) with line colors corresponding to the issue date, meaning the date the data was initially published or subsequently revised. Darker colors denote older versions, while lighter colors represent more recent updates.

If the lines perfectly overlap, it means the data for that location is not subject to revisions. Where you see a vertical spread along a single date on the x-axis, the data has been revised over time.

rev_by_loc <- rev_details$revision_behavior %>%
    filter(n_revisions != 0) %>%
    group_by(signal, geo_value) %>%
    summarize(
        n_versions = n(),
        n_rev = mean(n_revisions, na.rm = TRUE),
        min_lag = min(min_lag, na.rm = TRUE),
        max_lag = max(max_lag, na.rm = TRUE),
        spread = mean(spread, na.rm = TRUE),
        rel_spread = mean(rel_spread, na.rm = TRUE),
        lag_near_latest = mean(lag_near_latest, na.rm = TRUE),
        .groups = "drop"
    ) %>%
    mutate(signal_label = signal_labels[signal]) %>%
    arrange(desc(rel_spread))

if (nrow(rev_by_loc) > 0) {
    cat("The table below summarizes revision behavior per location. Column definitions:\n\n")
    cat("- Min Lag. Minimum days between the event date and its first publication.\n")
    cat("- Max Lag. Maximum days between the event date and its final revision.\n")
    cat("- Avg. Spread. Average absolute difference between the minimum and maximum published values.\n")
    cat("- Avg. Rel. Spread. Average relative difference between the minimum and maximum published values (spread divided by the maximum value).\n")
    cat(sprintf("- Lag Near Latest. Average days until the value stays within %d%% of the latest published value.\n\n", rev_details$within_latest * 100))

    rev_by_loc_display <- if (multi_signal) {
        rev_by_loc %>% select(signal_label, geo_value, n_versions:lag_near_latest)
    } else {
        rev_by_loc %>% select(geo_value, n_versions:lag_near_latest)
    }
    DT::datatable(rev_by_loc_display,
        colnames = c(
            if (multi_signal) "Signal" else NULL,
            "Location", "Versions", "Avg. Revisions", "Min Lag", "Max Lag",
            "Avg. Spread", "Avg. Rel. Spread", "Lag Near Latest"
        ),
        caption = "Candidate revision behavior by location",
        rownames = FALSE, options = list(scrollX = TRUE, scrollY = "400px", scrollCollapse = TRUE, paging = TRUE, dom = "tp")
    ) %>% DT::formatRound(columns = if (multi_signal) 3:9 else 2:8, digits = 2)
}

The table below summarizes revision behavior per location. Column definitions:

• Min Lag. Minimum days between the event date and its first publication.
• Max Lag. Maximum days between the event date and its final revision.
• Avg. Spread. Average absolute difference between the minimum and maximum published values.
• Avg. Rel. Spread. Average relative difference between the minimum and maximum published values (spread divided by the maximum value).
• Lag Near Latest. Average days until the value stays within 20% of the latest published value.

if (nrow(rev_by_loc) > 0 && n_locations > 1) {
    # Pick top locations by max rel_spread across signals
    rev_geos <- rev_by_loc %>%
        group_by(geo_value) %>%
        summarize(rel_spread = max(rel_spread, na.rm = TRUE), .groups = "drop") %>%
        arrange(desc(rel_spread)) %>%
        #sample_geos
        head(params$max_locations_plot) %>%
        pull(geo_value)

    # Limit to last year for clarity
    # For multi-signal, loop per signal to avoid autoplot dropdown issues with
    # the signal other-key (same pattern as eda-overview).
    archive_filtered <- archive %>%
        filter(geo_value %in% rev_geos) %>%
        filter(time_value >= max(time_value) - years(1))

    if (!multi_signal) {
        n_geos_rev <- n_distinct(archive_filtered$DT$geo_value)
        autoplot(archive_filtered, .versions = "month",
            .interactive = TRUE, .facet_to_dropdown = n_geos_rev > 1) %>%
            plotly::toWebGL() %>%
            plotly::partial_bundle()
    } else {
        rev_plots <- purrr::map(signals_vec, function(sig) {
            arch_sig_dt <- archive_filtered$DT %>%
                as_tibble() %>%
                filter(signal == sig) %>%
                select(-signal)
            if (nrow(arch_sig_dt) == 0) return(NULL)
            arch_sig <- as_epi_archive(arch_sig_dt)
            n_geos_sig <- n_distinct(arch_sig$DT$geo_value)
            p <- autoplot(arch_sig, .versions = "month",
                .interactive = TRUE, .facet_to_dropdown = n_geos_sig > 1) %>%
                plotly::toWebGL() %>%
                plotly::partial_bundle()
            htmltools::tagList(
                htmltools::tags$h5(signal_labels[sig]),
                p
            )
        }) %>%
            purrr::compact()
        htmltools::tagList(rev_plots)
    }
}

2.6 Signal sanity checks

This section is by its nature more difficult than others to do systematically. Things to look for:

• If it’s a percentile, is it between 0 and 100?
• If it’s a rate, is it roughly proportional to population? (requires external population data — inspect manually)
• Are the actual values reasonable? E.g. if it’s a percentage, does it ever approach 100 outside of literal pandemics?
• If it’s for a seasonal disease, are rates highest in the winter and lowest in the summer?
• Are geo-aggregated coarser levels plausible aggregations? (checked below for state-level indicators vs HHS/nation)

The histogram and table below show the distribution of all observed values and the value range per location for the candidate indicator, to assess whether the range is plausible given the indicator’s expected scale.

df_candidate_plot <- df_candidate %>%
    mutate(indicator = signal_labels[signal])

ggplot(df_candidate_plot, aes(x = value)) +
    geom_histogram(bins = 40, alpha = 0.8, color = "white") +
    scale_x_continuous(labels = label_number(scale_cut = cut_si(""))) +
    scale_y_continuous(labels = label_number(scale_cut = cut_si(""))) +
    {if (multi_signal) facet_wrap(~indicator, scales = "free", labeller = label_wrap_html(width = 25)) else NULL} +
    labs(
        title = if (multi_signal) "Value distribution by signal"
                else sprintf("Value distribution: %s", signal_labels[signals_vec]),
        subtitle = sprintf(
            "%s observations across %d locations",
            label_comma()(nrow(df_candidate)), n_locations
        ),
        x = "Value", y = "Count"
    )

The plot below shows the seasonal pattern of the indicator, with each panel representing one epidemiological year (June–May). The line traces the monthly median across all locations, the inner band covers the interquartile range, and the outer band spans the full observed range.

season_month_labels <- function(x) {
c("Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec", "Jan", "Feb", "Mar", "Apr", "May")[x]
}

df_seasonal <- df_candidate %>%
    mutate(
        yr = as.integer(format(time_value, "%Y")),
        month = as.integer(format(time_value, "%m")),
        season_start = if_else(month >= 6L, yr, yr - 1L),
        season = paste0(season_start, "-", substr(as.character(season_start + 1L), 3L, 4L)),
        season_month = if_else(month >= 6L, month - 5L, month + 7L)
    ) %>%
    group_by(signal, season, season_month) %>%
    summarize(
        median_val = median(value, na.rm = TRUE),
        min_val = min(value, na.rm = TRUE),
        max_val = max(value, na.rm = TRUE),
        p25 = quantile(value, 0.25, na.rm = TRUE),
        p75 = quantile(value, 0.75, na.rm = TRUE),
        .groups = "drop"
    ) %>%
    # Keep only signal-season combos with >= 3 months of data
    semi_join(
        count(., signal, season) %>% filter(n >= 3),
        by = c("signal", "season")
    ) %>%
    mutate(signal_label = signal_labels[signal])

ggplot(df_seasonal, aes(x = season_month, color = season, fill = season, group = season)) +
    geom_ribbon(aes(ymin = min_val, ymax = max_val), alpha = 0.10, linewidth = 0.2) +
    geom_ribbon(aes(ymin = p25, ymax = p75), alpha = 0.30, linewidth = 0.3) +
    geom_line(aes(y = median_val), linewidth = 0.8) +
    geom_point(aes(y = median_val), size = 1.5) +
    {
        if (multi_signal) facet_grid(signal_label ~ season, scales = "free_y")
        else facet_wrap(~season, ncol = 2, scales = "free_y")
    } +
    scale_x_continuous(labels = season_month_labels) +
    scale_y_continuous(labels = label_number(scale_cut = cut_si(""))) +
    labs(
        title = if (multi_signal) "Seasonal patterns by signal"
                else sprintf("Seasonal pattern: %s", signal_labels[signals_vec]),
        subtitle = "Median by month across all locations — inner band: IQR, outer band: full range",
        x = "Month", y = "Median value"
    ) +
    theme(legend.position = "none")

The table below flags locations whose seasonal peak is at least twice the typical (median) peak across all locations in the same season.