Indicator Comparison: Doctor Visits: Smoothed Adj CLI vs JHU COVID-19 Cases

This notebook evaluates a candidate indicator against a guiding indicator (presumed ground truth) to determine whether the candidate has value for nowcasting or forecasting.

---

1 Data Loading

# Load guiding indicator
if (!is.null(params$guiding_csv)) {
    df_guiding <- read.csv(params$guiding_csv) %>%
        mutate(time_value = as.Date(time_value)) %>%
        filter(time_value >= as.Date(params$start_day) & time_value <= as.Date(params$end_day)) %>%
        as_epi_df()
} else {
    # Determine if new API can be used for guiding source
    new_api_sources <- c("nssp", "nhsn", "pophive", "nwss", "claims_data_outpatient")
    use_new_api <- FALSE
    if (params$guiding_source %in% new_api_sources) {
        meta_available <- tryCatch(
            {
                meta_res <- epidata_meta(params$guiding_source)
                params$guiding_indicator %in% meta_res[[params$guiding_source]]$signals
            },
            error = function(e) FALSE
        )
        if (meta_available) use_new_api <- TRUE
    }
    if (use_new_api) {
        df_guiding <- epidata_snapshot(
            source = params$guiding_source,
            signals = params$guiding_indicator,
            geo_type = params$geo_type
        ) %>%
            rename(any_of(c(time_value = "reference_time", version = "report_time")))
        if ("fill_method" %in% colnames(df_guiding)) {
            df_guiding <- df_guiding %>%
                filter(fill_method == "source" | !any(fill_method == "source"))
        }
        if ("version" %in% colnames(df_guiding)) {
            df_guiding <- df_guiding %>%
                group_by(geo_value, time_value) %>%
                slice_max(version, n = 1, with_ties = FALSE) %>%
                ungroup()
        }
        df_guiding <- df_guiding %>%
            filter(
                time_value >= as.Date(params$start_day),
                time_value <= as.Date(params$end_day)
            ) %>%
            rename(any_of(c(time_value = "reference_time"))) %>%
            select(geo_value, time_value, value) %>%
            as_epi_df()
    } else {
        df_guiding <- pub_covidcast(
            params$guiding_source, params$guiding_indicator,
            params$geo_type, params$time_type,
            time_values = time_range(params$start_day, params$end_day),
            as_of = Sys.Date()
        ) %>% as_epi_df()
    }
}

cat(sprintf(
    "- **Summary**: %s rows across %s locations\n",
    label_comma()(nrow(df_guiding)), label_comma()(n_distinct(df_guiding$geo_value))
))

• Summary: 104,896 rows across 3,278 locations

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

    # 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))
        candidate_archive <- df_candidate_raw %>%
            select(geo_value, time_value, version, value) %>%
            as_epi_archive()
        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()
    }
} 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$candidate_source %in% new_api_sources) {
        meta_available <- tryCatch(
            {
                meta_res <- epidata_meta(params$candidate_source)
                params$candidate_indicator %in% meta_res[[params$candidate_source]]$signals
            },
            error = function(e) FALSE
        )
        if (meta_available) use_new_api <- TRUE
    }
    if (use_new_api) {
        df_candidate <- epidata_snapshot(
            source = params$candidate_source,
            signals = params$candidate_indicator,
            geo_type = params$geo_type
        ) %>%
            rename(any_of(c(time_value = "reference_time", version = "report_time")))
        if ("fill_method" %in% colnames(df_candidate)) {
            df_candidate <- df_candidate %>%
                filter(fill_method == "source" | !any(fill_method == "source"))
        }
        if ("version" %in% colnames(df_candidate)) {
            df_candidate <- df_candidate %>%
                group_by(geo_value, time_value) %>%
                slice_max(version, n = 1, with_ties = FALSE) %>%
                ungroup()
        }
        df_candidate <- df_candidate %>%
            filter(
                time_value >= as.Date(params$start_day),
                time_value <= as.Date(params$end_day)
            ) %>%
            rename(any_of(c(time_value = "reference_time"))) %>%
            select(geo_value, time_value, value) %>%
            as_epi_df()
    } else {
        df_candidate <- pub_covidcast(
            params$candidate_source, params$candidate_indicator,
            params$geo_type, params$time_type,
            time_values = time_range(params$start_day, params$end_day),
            as_of = Sys.Date()
        ) %>% as_epi_df()
    }
}

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

• Summary: 67,798 rows across 2,330 locations

# Infer time type and geo type from guiding (defaults/overrides resolved here before archive query)
env_time_type <- attributes(df_guiding)$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_guiding)$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
}

# Fetch revision history from API if not already built from CSV
if (is.null(candidate_archive) && is.null(params$candidate_csv)) {
    all_geos <- unique(df_candidate$geo_value)
    archive_geos <- random_sample(all_geos, params$max_archive_locs)
    # Adjust lookback window
    lookback <- switch(params$time_type,
        "day" = lubridate::weeks(31),
        "week" = lubridate::weeks(150),
        "month" = lubridate::months(36),
        lubridate::weeks(31)
    )

    archive_start <- as.character(max(
        as.Date(params$start_day),
        as.Date(params$end_day) - lookback
    ))

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

    # Request only for the sampled locations and the relevant time range
    # to avoid server-side timeouts/errors on large indicators.
    # Determine if new API can be used for candidate source archive
    new_api_sources <- c("nssp", "nhsn", "pophive", "nwss", "claims_data_outpatient")
    use_new_api <- FALSE
    if (params$candidate_source %in% new_api_sources) {
        meta_available <- tryCatch(
            {
                meta_res <- epidata_meta(params$candidate_source)
                params$candidate_indicator %in% meta_res[[params$candidate_source]]$signals
            },
            error = function(e) FALSE
        )
        if (meta_available) use_new_api <- TRUE
    }
    if (use_new_api) {
        df_archive_raw <- epidata_archive(
            source = params$candidate_source,
            signals = params$candidate_indicator,
            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)
            )
    } else {
        df_archive_raw <- pub_covidcast(
            source = params$candidate_source,
            signal = params$candidate_indicator,
            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, value)

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

# Retrieve time_type (already resolved in load-candidate chunk)
env_time_type <- params$time_type

# When time_type is weekly, align all dates to Sundays and aggregate
if (env_time_type == "week") {
    df_guiding <- df_guiding %>%
        mutate(time_value = round_to_sunday(time_value)) %>%
        group_by(geo_value, time_value) %>%
        summarize(value = mean(value, na.rm = TRUE), .groups = "drop") %>%
        as_epi_df()

    df_candidate <- df_candidate %>%
        mutate(time_value = round_to_sunday(time_value)) %>%
        group_by(geo_value, time_value) %>%
        summarize(value = mean(value, na.rm = TRUE), .groups = "drop") %>%
        as_epi_df()

    cat("- Aligned both datasets to Sunday-start weeks  \n")
}

# Combine into a single epi_df
grid_start <- as.Date(params$start_day)
grid_end <- as.Date(params$end_day)
if (env_time_type == "week") {
    grid_start <- round_to_sunday(grid_start)
    grid_end <- round_to_sunday(grid_end)
}

df <- df_guiding %>%
    rename(guiding = value) %>%
    full_join(
        df_candidate %>% select(geo_value, time_value, value) %>%
            rename(candidate = value),
        by = c("geo_value", "time_value")
    ) %>%
    group_by(geo_value) %>%
    complete(
        time_value = seq.Date(
            from = grid_start,
            to = grid_end,
            by = env_time_type
        )
    ) %>%
    ungroup() %>%
    as_epi_df()

n_locations <- n_distinct(df$geo_value)
cat(sprintf(
    "- **Combined (full join)**: %s rows, %s locations (%s to %s)  \n",
    label_comma()(nrow(df)), label_comma()(n_locations),
    min(df$time_value), max(df$time_value)
))

• Combined (full join): 104,960 rows, 3,280 locations (2023-01-01 to 2023-02-01)

# Intersected where both indicators must be present
df_inter <- df %>%
    filter(!is.na(candidate), !is.na(guiding))

cat(sprintf(
    "- **Intersection (both present)**: %s rows (%.1f%% of full grid)\n\n",
    label_comma()(nrow(df_inter)), 100 * nrow(df_inter) / nrow(df)
))

• Intersection (both present): 67,737 rows (64.5% of full grid)

start_day <- min(df$time_value)
end_day <- max(df$time_value)
n_inter_locations <- n_distinct(df_inter$geo_value)

# Define adaptive units and lag ranges for downstream analysis
time_unit <- switch(env_time_type %||% "unknown",
    week = "weeks",
    month = "months",
    year = "years",
    "days" # default
)
lag_unit_label <- sprintf("Lag (%s)", time_unit)

# Standard discrete lags to check
lag_values <- switch(env_time_type %||% "unknown",
    week = -2:2,
    month = -2:2,
    year = -1:1,
    c(-14, -7, 0, 7, 14) # default
)

# Systematic sweep range (includes negative lags/leads)
lags_sweep <- switch(env_time_type %||% "unknown",
    week = -5:5,
    month = -5:5,
    year = -2:2,
    -21:21 # default
)

---

2 Exploratory Data Analysis

The data is aggregated at the county level. The intersection of candidate and guiding datasets contains 2328 valid locations measured at a day-level frequency.

2.1 Time Series Overview

This section provides a direct visual comparison of the guiding and candidate indicators over time for a subset of observation locations. We start by looking at the raw values for both indicators across a sample of locations. Then, a third plot overlays both indicators after normalizing them.

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

df_long <- df_inter %>%
    filter(geo_value %in% sample_geos) %>%
    pivot_longer(
        cols = c(candidate, guiding),
        names_to = "indicator", values_to = "value"
    ) %>%
    group_by(geo_value, indicator) %>%
    mutate(value_norm = (value - mean(value, na.rm = TRUE)) / sd(value, na.rm = TRUE)) %>%
    ungroup() %>%
    mutate(indicator = recode(indicator,
        candidate = params$candidate_name,
        guiding = params$guiding_name
    ))

Note

Showing a random sample of 60 locations (out of 2328 total) to keep plots readable. You can adjust this with the max_locations_plot parameter.

In this normalized view, overlapping patterns, lead-lag relationships, or divergent behavior become easier to spot. For instance, if the candidate curve tends to rise before the guiding curve, we might have a useful leading relationship.

autoplot(df_long, value_norm,
.interactive = T, .facet_by = "geo_value", .facet_to_dropdown = T) %>%
plotly::layout(yaxis = list(title = "Normalized Value")) %>%
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")

df_quantiles_wide <- df_inter %>%
    pivot_longer(
        cols = c(candidate, guiding),
        names_to = "indicator", values_to = "value"
    ) %>%
    group_by(time_value, indicator) %>%
    summarize(vals = list(value), .groups = "drop") %>%
    mutate(geo_value = indicator) %>%
    as_epi_df()

# 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) %>%
    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 = recode(geo_value,
            candidate = params$candidate_name,
            guiding = params$guiding_name
        )
    )

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", labeller = label_wrap_html(width = 25)) +
    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_inter_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) for each indicator.

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 step chart below tracks the number of geographic locations with valid data on any given date for the candidate indicator, the guiding indicator, and instances where both are available together. The dashed line indicates when both indicators are available simultaneously.

On average, both indicators are simultaneously available in 2116.8 locations per day.

ggplot(coverage, aes(x = time_value)) +
    geom_area(
        aes(
            y = candidate_available, color = params$candidate_name,
            fill = params$candidate_name
        ),
        linewidth = 0.7, alpha = 0.05
    ) +
    geom_area(
        aes(
            y = guiding_available, color = params$guiding_name,
            fill = params$guiding_name
        ),
        linewidth = 0.7, alpha = 0.05
    ) +
    geom_area(
        aes(y = both_available, color = "Both", fill = "Both"),
        linewidth = 0.7, alpha = 0.05, linetype = "dashed"
    ) +
    scale_x_date(breaks = breaks_pretty(), label = label_date_short()) +
    scale_y_continuous(labels = label_comma(), limits = c(0, NA)) +
    scale_colour_manual(
        name = "Indicator",
        values = dis_pal_list,
        breaks = c(
            params$guiding_name,
            params$candidate_name, "Both"
        )
    ) +
    scale_fill_manual(
        name = "Indicator",
        values = dis_pal_list,
        breaks = c(
            params$guiding_name,
            params$candidate_name, "Both"
        )
    ) +
    labs(
        title = "Data coverage over time",
        subtitle = sprintf("Number of locations with non-missing values per %s", env_time_type),
        x = "Date", y = "Number of locations"
    ) +
    guides(
        color = guide_legend(override.aes = list(alpha = 1, fill = NA)),
        fill = guide_legend(override.aes = list(alpha = 0.2))
    )

2.3.1 Summary Statistics by Location

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

Excluded locations (no candidate data): 02220 (Sitka City and Borough), 06049 (Modoc County), 08111 (San Juan County), 15005 (Kalawao County), 16067 (Minidoka County), 17003 (Alexander County), 17139 (Moultrie County), 20039 (Decatur County), 21209 (Scott County), 26013 (Baraga County), 27075 (Lake County), 38055 (McLean County), 40059 (Harper County), 45065 (McCormick County), 54065 (Morgan County), etc.

Excluded locations (no guiding data): 78010 (Saint Croix Island), 78030 (Saint Thomas Island)

Zero variance locations (candidate flatlined): 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 (Candidate) vs 3278 locations (Guiding)
• 5-25% missing: 213 locations (Candidate) vs 0 locations (Guiding)
• 25-50% missing: 140 locations (Candidate) vs 0 locations (Guiding)
• 50-99.9% missing: 171 locations (Candidate) vs 0 locations (Guiding)
• 100% missing: 950 locations (Candidate) vs 2 locations (Guiding)

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 location using the intersected data
noise_stats <- df_inter %>%
    arrange(geo_value, time_value) %>%
    group_by(geo_value) %>%
    summarize(
        cand_cv = sd(candidate, na.rm = TRUE) /
            abs(mean(candidate, na.rm = TRUE)),
        guid_cv = sd(guiding, na.rm = TRUE) /
            abs(mean(guiding, na.rm = TRUE)),
        .groups = "drop"
    ) %>%
    # Remove locations where metrics couldn't be computed
    filter(if_all(c(cand_cv, guid_cv), ~ is.finite(.x)))

# Clip extreme outliers for robust visualization
cv_upper <- quantile(c(noise_stats$cand_cv, noise_stats$guid_cv), 0.99, na.rm = TRUE)
noise_stats_plot <- noise_stats %>%
    filter(cand_cv <= cv_upper, guid_cv <= cv_upper)

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

# Summary Table
noise_summary <- tibble(
    Metric = "Coefficient of Variation (CV)",
    `Median (candidate)` = median(noise_stats$cand_cv),
    `Median (guiding)` = median(noise_stats$guid_cv),
    `Median Difference` = median(noise_stats$cand_cv - noise_stats$guid_cv),
    `% Smoother (candidate < guiding)` = mean(noise_stats$cand_cv < noise_stats$guid_cv) * 100
)

DT::datatable(noise_summary,
    caption = sprintf(
        "Comparative CV summary across locations.%s\nA negative Median Difference and high %% Smoother indicate the candidate is typically less volatile.",
        if (n_excluded > 0) sprintf(" (%d extreme outlier(s) excluded from plots.)", n_excluded) else ""
    ),
    rownames = FALSE, options = list(scrollX = TRUE, paging = FALSE, dom = "t")
) %>% DT::formatRound(columns = 2:5, digits = 3)

# Calculate CV per location using the intersected data
# Reshape for plotting
noise_long <- noise_stats_plot %>%
    pivot_longer(-geo_value, names_to = "indicator", values_to = "value") %>%
    mutate(
        indicator = recode(indicator,
            cand_cv = params$candidate_name,
            guid_cv = params$guiding_name
        )
    )

# Histogram comparison
p_hist <- ggplot(noise_long, aes(x = value, fill = indicator)) +
    geom_histogram(bins = 30, color = "white", alpha = 0.7, position = "identity") +
    facet_wrap(~indicator, ncol = 1, labeller = label_wrap_html(width = 25)) +
    scale_x_continuous(labels = label_percent()) +
    labs(
        title = "CV Distribution by Location",
        subtitle = "Lower values indicate a more stable indicator",
        x = "Coefficient of Variation (CV)", y = "Count", fill = "Indicator"
    ) +
    theme(legend.position = "top")

# Paired Scatter Plot
p_scatter <- ggplot(noise_stats_plot, aes(x = guid_cv, y = cand_cv)) +
    geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "grey50") +
    geom_point(aes(color = cand_cv < guid_cv), alpha = 0.7, size = 2) +
    scale_x_continuous(labels = label_percent()) +
    scale_y_continuous(labels = label_percent()) +
    scale_color_manual(
        name = "Result",
        values = c("TRUE" = dis_pal_list[2], "FALSE" = dis_pal_list[4]),
        labels = c("TRUE" = "Candidate More Stable", "FALSE" = "Guiding More Stable")
    ) +
    labs(
        title = "Paired CV Comparison",
        subtitle = "Points below the dashed line represent locations where the candidate is more stable.",
        x = sprintf("CV: %s (Guiding)", params$guiding_name),
        y = sprintf("CV: %s (Candidate)", params$candidate_name)
    )

p_hist / p_scatter + plot_layout(heights = c(1.2, 1))

3 Correlation Analysis

All correlation metrics are computed on the intersection of the candidate and guiding indicators (rows where both values are non-missing) using the latest available revisions unless otherwise indicated.

3.1 Correlation over Time

Here we calculate the Spearman rank correlation between the candidate and guiding indicators across all available locations for each individual date. We also overlay a LOESS smoothed trend line to help capture broader patterns over time.

Values near one tell us there is a strong positive correlation across areas on that specific day. When you see dips or sustained drops in the line, it suggests that the candidate decoupled from the guiding indicator during that period.

Two additional panels are included below the correlation plot to support interpretation. The middle panel shows smoothed quantile trends for both indicators. The bottom panel shows the number of locations where both indicators were available on each date.

mean_rho <- NA
cor_by_time <- epi_cor(df_inter, candidate, guiding,
    cor_by = time_value, method = "spearman"
)

mean_rho <- mean(cor_by_time$cor, na.rm = TRUE)

p_cor <- ggplot(cor_by_time, aes(x = time_value)) +
    geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
    geom_hline(
        yintercept = c(0.3, 0.4, 0.6, 0.7, 0.8, 0.9),
        linetype = "dotted", color = "gray85", linewidth = 0.3
    ) +
    geom_line(aes(y = cor, color = "Correlation"), linewidth = 0.8) +
    geom_smooth(aes(y = cor, color = "LOESS smoothed trend"),
        method = "loess", span = 0.2, se = TRUE,
        fill = dis_pal_list[2], alpha = 0.15
    ) +
    scale_x_date(
        limits = range(cor_by_time$time_value, na.rm = TRUE),
        breaks = breaks_pretty(n = 6), label = label_date_short()
    ) +
    scale_y_continuous(
        name = expression(paste("Spearman ", rho)),
        limits = c(floor(min(cor_by_time$cor, na.rm = TRUE) * 10) / 10, 1)
    ) +
    scale_color_manual(
        name = "Metric",
        values = c(
            "Correlation" = dis_pal_list[1],
            "LOESS smoothed trend" = dis_pal_list[2]
        )
    ) +
    labs(
        title = "Correlation over time",
        subtitle = bquote("Spearman" ~ rho ~ "across" ~ .(n_inter_locations) ~ "locations | Mean =" ~ .(sprintf("%.3f", mean_rho))),
        x = "Date", y = expression(paste("Spearman ", rho, " across locations"))
    ) +
    theme(
        legend.position = "none",
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank()
    )

p_cor <- ggplot(cor_by_time, aes(x = time_value)) +
    geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
    geom_hline(
        yintercept = c(0.3, 0.4, 0.6, 0.7, 0.8, 0.9),
        linetype = "dotted", color = "gray85", linewidth = 0.3
    ) +
    geom_line(aes(y = cor, color = "Correlation"), linewidth = 0.8) +
    geom_smooth(aes(y = cor, color = "LOESS smoothed trend"),
        method = "loess", span = 0.2, se = TRUE,
        fill = dis_pal_list[2], alpha = 0.15
    ) +
    scale_x_date(
        limits = range(cor_by_time$time_value, na.rm = TRUE),
        breaks = breaks_pretty(n = 6), label = label_date_short()
    ) +
    scale_y_continuous(
        name = expression(paste("Spearman ", rho)),
        limits = c(floor(min(cor_by_time$cor, na.rm = TRUE) * 10) / 10, 1)
    ) +
    scale_color_manual(
        name = "Metric",
        values = c(
            "Correlation" = dis_pal_list[1],
            "LOESS smoothed trend" = dis_pal_list[2]
        )
    ) +
    labs(
        title = "Diagnostic Comparison: Correlation, Quantiles, and Coverage",
        subtitle = bquote("Correlation over time | Mean" ~ rho == .(sprintf("%.3f", mean_rho)) ~ "(latest revisions)"),
        x = NULL
    ) +
    theme(
        legend.position = "none",
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank()
    )

p_cor / p_quant / p_cov + plot_layout(heights = c(1.2, 0.4, 0.4))

3.2 Distribution of Temporal Correlations

This histogram shows the distribution of the Spearman correlation computed within each geographic location over the full time period. The dashed line marks the median.

cor_by_geo <- epi_cor(df_inter, candidate, guiding,
    cor_by = geo_value, method = "spearman"
)

med_rho <- median(cor_by_geo$cor, na.rm = TRUE)

if (n_inter_locations > 1) {
    p_rho <- ggplot(cor_by_geo, aes(x = cor)) +
        geom_histogram(bins = 20, alpha = 0.6, color = "white") +
        geom_vline(
            xintercept = med_rho, linetype = "dashed",
            color = dis_pal_list[2], linewidth = 1
        ) +
        scale_x_continuous(breaks = breaks_pretty()) +
        labs(
            subtitle = bquote("Median Spearman" ~ rho == .(sprintf("%.3f", med_rho)) ~ "across" ~ .(nrow(cor_by_geo)) ~ "locations"),
            x = expression(Spearman ~ rho), y = "Count"
        )

    p_rho + labs(title = "Distribution of intra-location temporal correlations")
} 

if (!(n_inter_locations > 1)) {
    cat(sprintf("For the single location %s, the temporal Spearman correlation is **%.3f**\n", unique(df_inter$geo_value)[1], med_rho))
}

3.3 Lagged Correlations

To explore lead and lag relationships, we shift the candidate indicator by various intervals (-14, -7, 0, 7, 14 days). Specifically, a positive lag (e.g., +7 days) compares the candidate indicator from 7 days ago against the guiding indicator today. This allows us to see if the candidate acts as a leading indicator.

3.3.1 Correlation over time, at different days

This plot expands on the basic correlation by showing multiple smoothed trend lines, each representing a different lag applied to the candidate. If the curve for a positive lag is consistently higher than the zero-lag curve, it is a strong indication that the candidate indicator provides a lead over the guiding indicator.

cor_lagged_time <- map(lag_values, ~ {
    epi_cor(df_inter, candidate, guiding,
        cor_by = time_value, dt1 = -.x, method = "spearman"
    ) %>%
        mutate(lag = .x)
}) %>%
    list_rbind() %>%
    mutate(lag = as.factor(lag))

if (n_inter_locations > 1) {
    ggplot(cor_lagged_time, aes(x = time_value, y = cor, color = lag)) +
        geom_line(alpha = 0.4) +
        geom_smooth(method = "loess", span = 0.2, se = FALSE, linewidth = 1.2) +
        geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
        scale_x_date(breaks = breaks_pretty(), label = label_date_short()) +
        scale_y_continuous(limits = c(-1, 1), breaks = breaks_pretty()) +
        labs(
            title = sprintf("Correlation over time at different %s", time_unit),
            subtitle = "Transparent line: Raw Correlation, Solid line: LOESS smoothed trend",
            x = "Date", y = expression(paste("Spearman ", rho)),
            color = lag_unit_label
        ) +
        guides(color = guide_legend(nrow = 1, override.aes = list(linewidth = 1.5, alpha = 1)))
}

3.3.2 Real-time lagged correlations

Same as above, but the candidate (Doctor Visits: Smoothed Adj CLI) values are drawn from epix_as_of(candidate_archive, v) at each version date v, so the correlation reflects what was observable in real time rather than using the final revised values.

versions <- sort(unique(candidate_archive$DT$version))

cor_lagged_realtime <- map_dfr(versions, function(vv) {
    snap <- epix_as_of(candidate_archive, vv)
    if (env_time_type == "week") {
        snap <- snap %>%
            mutate(time_value = round_to_sunday(time_value)) %>%
            group_by(geo_value, time_value) %>%
            summarize(value = mean(value, na.rm = TRUE), .groups = "drop")
    }
    valid_tvs <- snap$time_value[snap$time_value %in% df_guiding$time_value]
    if (length(valid_tvs) == 0) {
        return(tibble(version = vv, lag = lag_values, cor = NA_real_))
    }
    latest_tv <- max(valid_tvs, na.rm = TRUE)
    map_dfr(lag_values, function(lg) {
        target_tv <- switch(env_time_type,
            week  = latest_tv - lubridate::weeks(lg),
            month = latest_tv - lubridate::months(lg),
            year  = latest_tv - lubridate::years(lg),
            latest_tv - lg
        )
        df_cand <- snap %>%
            filter(time_value == target_tv) %>%
            rename(candidate = value)
        df_guid <- df_guiding %>%
            filter(time_value == latest_tv) %>%
            select(geo_value, guiding = value)
        df_joined <- inner_join(df_cand, df_guid, by = "geo_value")
        if (n_distinct(df_joined$geo_value) < 3) {
            return(tibble(version = vv, lag = lg, cor = NA_real_))
        }
        tibble(version = vv, lag = lg, cor = cor(df_joined$candidate, df_joined$guiding,
            method = "spearman", use = "complete.obs"
        ))
    })
}) %>%
    mutate(lag = as.factor(lag))

ggplot(cor_lagged_realtime, aes(x = version, y = cor, color = lag)) +
    geom_line(alpha = 0.4) +
    geom_smooth(method = "loess", span = 0.2, se = FALSE, linewidth = 1.2) +
    geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
    scale_x_date(breaks = breaks_pretty(), label = label_date_short()) +
    scale_y_continuous(limits = c(-1, 1), breaks = breaks_pretty()) +
    labs(
        title = sprintf("Real-time correlation over time at different %s", time_unit),
        subtitle = "Candidate values as-published at each version date",
        x = "Version date", y = expression(paste("Spearman ", rho)),
        color = lag_unit_label
    ) +
    guides(color = guide_legend(nrow = 1, override.aes = list(linewidth = 1.5, alpha = 1)))

3.3.3 Per-location temporal correlation distributions, at different lags

Here we show density distributions of the temporal Spearman correlations for each location, calculated across several different lead and lag shifts.

The goal is to see which density curve has the highest correlation values. The lag associated with that curve is the one that best aligns the temporal dynamics of the two indicators for the majority of areas. This complements the previous plot by confirming whether the optimal lag remains consistent over time within individual locations.

cor_lagged_geo <- map(lag_values, ~ {
    epi_cor(df_inter, candidate, guiding,
        cor_by = geo_value, dt1 = -.x, method = "spearman"
    ) %>%
        mutate(lag = .x)
}) %>%
    list_rbind() %>%
    mutate(lag = as.factor(lag))

ggplot(cor_lagged_geo, aes(x = cor)) +
    geom_density(aes(fill = lag, color = lag), alpha = 0.15) +
    xlim(-1, 1) +
    labs(
        title = sprintf("Temporal correlations at different %s", time_unit),
        subtitle = "Positive lag = candidate leads the guiding indicator",
        x = expression(paste("Spearman ", rho)), y = "Density",
        fill = lag_unit_label, color = lag_unit_label
    ) +
    guides(
        fill = guide_legend(override.aes = list(alpha = 0.35))
    )

3.4 Systematic Lag Analysis

In this section, we sweep through a wider range of lags (from -21 to 21 days) to pinpoint the optimal lead time.

if (n_inter_locations > 1) {
    cat("The plot displays summary statistics like the mean, median, and interquartile range of the temporal correlations across all locations, calculated for every possible lag in the window. A vertical line marks the optimal lag, which has the highest mean correlation.\n\n")
} else {
    cat("The plot displays the temporal correlation for the single location, calculated for every possible lag in the window. A vertical line marks the optimal lag, which has the highest correlation.\n\n")
}

The plot displays summary statistics like the mean, median, and interquartile range of the temporal correlations across all locations, calculated for every possible lag in the window. A vertical line marks the optimal lag, which has the highest mean correlation.

cat("If the optimal lag is positive, the candidate's past values are most predictive of the guiding indicator's current values (which makes them valuable for forecasting).\n\n")

If the optimal lag is positive, the candidate’s past values are most predictive of the guiding indicator’s current values (which makes them valuable for forecasting).

lags <- lags_sweep

cor_sweep <- map(lags, ~ {
    epi_cor(df_inter, candidate, guiding,
        cor_by = geo_value, dt1 = -.x, method = "spearman"
    ) %>%
        mutate(lag = .x)
}) %>% list_rbind()

lag_summary <- cor_sweep %>%
    group_by(lag) %>%
    summarize(
        mean_cor = mean(cor, na.rm = TRUE),
        median_cor = median(cor, na.rm = TRUE),
        q25 = quantile(cor, 0.25, na.rm = TRUE),
        q75 = quantile(cor, 0.75, na.rm = TRUE),
        .groups = "drop"
    )

optimal_lag <- lag_summary$lag[which.max(lag_summary$mean_cor)]

if (n_inter_locations > 1) {
    ggplot(lag_summary, aes(x = lag)) +
        geom_ribbon(aes(ymin = q25, ymax = q75, fill = "IQR"), alpha = 0.2) +
        geom_line(aes(y = mean_cor, color = "Mean"), linewidth = 1.2) +
        geom_point(aes(y = mean_cor, color = "Mean"), size = 2, fill = "transparent") +
        geom_line(aes(y = median_cor, color = "Median"),
            linetype = "dashed", linewidth = 0.8
        ) +
        geom_vline(xintercept = optimal_lag, linetype = "dotted", color = "black") +
        annotate("text",
            x = optimal_lag, y = -1, label = sprintf("Optimal: %d", optimal_lag),
            hjust = -0.1, vjust = -0.5, size = 4
        ) +
        scale_x_continuous(breaks = breaks_pretty()) +
        scale_y_continuous(limits = c(-1, 1), breaks = breaks_pretty()) +
        labs(
            title = sprintf("Summary metrics vs. %s", lag_unit_label),
            x = lag_unit_label, y = expression(paste("Spearman ", rho)),
            color = "", fill = ""
        )
} else {
    ggplot(lag_summary, aes(x = lag)) +
        geom_line(aes(y = mean_cor, color = "Correlation"), linewidth = 1.2) +
        geom_point(aes(y = mean_cor, color = "Correlation"), size = 2, fill = "transparent") +
        geom_vline(xintercept = optimal_lag, linetype = "dotted", color = "black") +
        annotate("text",
            x = optimal_lag, y = -1, label = sprintf("Optimal: %d", optimal_lag),
            hjust = -0.1, vjust = -0.5, size = 4
        ) +
        scale_x_continuous(breaks = breaks_pretty()) +
        scale_y_continuous(limits = c(-1, 1), breaks = breaks_pretty()) +
        labs(
            title = sprintf("Correlation vs. %s", lag_unit_label),
            x = lag_unit_label, y = expression(paste("Spearman ", rho)),
            color = ""
        )
}

3.5 Choropleth Map at Optimal Lag

This geographic map colors each county according to its temporal Spearman correlation with the guiding indicator, specifically evaluated at the optimal lag we found in the previous section.

cor_optimal <- epi_cor(df_inter, candidate, guiding,
    cor_by = geo_value, dt1 = -optimal_lag,
    method = "spearman"
)

if (n_inter_locations > 1 && params$geo_type %in% c("state", "county")) {
    # Basic check to see if we have valid-looking state/county codes before trying geographic match
    looks_like_state <- params$geo_type == "state" && all(nchar(cor_optimal$geo_value) == 2)
    looks_like_county <- params$geo_type == "county" && all(nchar(cor_optimal$geo_value) == 5)

    if (looks_like_state || looks_like_county) {
        cor_choro <- cor_optimal %>%
            rename(value = cor) %>%
            mutate(
                time_value = as.Date(params$start_day),
                issue = as.Date(params$start_day),
                data_source = "evaluation",
                signal = "correlation"
            )

        attributes(cor_choro)$metadata$geo_type <- params$geo_type
        class(cor_choro) <- c("covidcast_signal", "data.frame")

        # Plot object evaluating as the last statement in the conditional
        print(plot(cor_choro,
            range = c(-1, 1), choro_col = c(dis_pal_list[3], "white", dis_pal_list[2]),
            title = sprintf("Per-location correlation at optimal lag (%d %s)", optimal_lag, time_unit)
        )) # For some reason this plot is only showing up when using print()
    }
}

# Summary table of correlations at optimal lag
cor_top <- cor_optimal %>%
    mutate(geo_name = geo_to_name(geo_value)) %>%
    arrange(desc(cor)) %>%
    select(geo_value, geo_name, cor)

DT::datatable(cor_top,
    colnames = c("Geo-level", "Name", "Spearman Correlation"),
    caption = sprintf("Strongest local correlations at %d %s lag", optimal_lag, time_unit),
    rownames = FALSE, options = list(scrollX = TRUE, scrollY = "400px", scrollCollapse = TRUE, paging = TRUE, dom = "tp")
) %>% DT::formatRound(columns = 3, digits = 3)

---

4 Summary

if (n_inter_locations > 1) {
    findings <- tibble(
        Metric = c(
            "Mean daily spatial correlation (across locations)",
            sprintf("Optimal lead time (%s)", time_unit),
            "Mean temporal correlation (at optimal lag)",
            "Median temporal correlation (lag 0)",
            "Median temporal correlation (optimal lag)",
            "Locations with |temporal ρ| > 0.5 (optimal lag)"
        ),
        Value = c(
            sprintf("%.3f", mean_rho),
            as.character(optimal_lag),
            sprintf("%.3f", max(lag_summary$mean_cor)),
            sprintf("%.3f", med_rho),
            sprintf("%.3f", median(cor_optimal$cor, na.rm = TRUE)),
            sprintf(
                "%d / %d (%.0f%%)",
                sum(abs(cor_optimal$cor) > 0.5, na.rm = TRUE),
                nrow(cor_optimal),
                100 * mean(abs(cor_optimal$cor) > 0.5, na.rm = TRUE)
            )
        )
    )
} else {
    findings <- tibble(
        Metric = c(
            "Temporal rank correlation (ρ)",
            sprintf("Optimal lead time (%s)", time_unit),
            "Temporal ρ at optimal lag"
        ),
        Value = c(
            sprintf("%.3f", med_rho),
            as.character(optimal_lag),
            sprintf("%.3f", max(lag_summary$mean_cor))
        )
    )
}

DT::datatable(findings, caption = "Key evaluation metrics", rownames = FALSE, options = list(scrollX = TRUE, paging = FALSE, dom = "t"))