Dataset Comparison

Doctor Visits (Outpatient Claims): Legacy vs. New ETL Pipeline

1 Doctor Visits (hrr)

This document compares the Doctor Visits (Outpatient Claims) indicator between the legacy covidcast pipeline and the new etl pipeline, both derived from claims data at the hrr geo level.

1.1 Methodological Differences

• The most visible difference is in how smoothing is applied. The legacy pipeline uses a one-sided Gaussian linear filter (h = 100) applied to the rate after division (smooth.py:1, sensor.py:24), which produces a wide, lagged smoother that blunts peaks and troughs by several days. The new pipeline instead computes a 7-day trailing sum on raw counts before dividing (common/calculate_parallel_seven_day_sum.sql:83), using a fixed RANGE BETWEEN INTERVAL '6 day' PRECEDING AND CURRENT ROW window (common/calculate_parallel_seven_day_sum.sql:92).
• The minimum observation requirements also differ: the legacy pipeline requires at least 3 days with counts (MIN_RECENT_OBS, config.py:40), while the new requires at least 5 of 7 (min_allowed_week_size = 5, common/calculate_parallel_seven_day_sum.sql:5).
• The old pipeline extends the denominator window by up to 7 days until it accumulates 500 visits (MIN_CUM_VISITS, sensor.py:81, config.py:37), while the new uses a fixed window.
• The rate formula also differs. The legacy pipeline applies a Jeffreys prior ((num + 0.5) / (denom + 1) in sensor.py:163) to inflate values slightly in sparse geographies. The new pipeline uses a plain ratio (claims/censor_ratio_insert.sql:15), so values will tend to be lower where counts are small.
• Denominator censoring also differs. The old pipeline suppresses when the smoothed denominator falls below 100 or fewer than 3 observation days are available (sensor.py:83, config.py:39), while the new deletes any row where the 7-day summed denominator is below CENSOR_THRESHOLD = 100 (claims/censor_ratio_insert.sql:4, constants.py:48).
• The COVID-19 case definition changes over time in the new pipeline but not the old. The legacy pipeline always combines Covid_like + Flu_like + Mixed (config.py:20, sensor.py:188). The new pipeline distinguishes two eras (claims/censor_ratio_insert.sql:15, constants.py:49). Before COVID_INCLUDES_UNKNOWN_END_DATE = 2021-06-01, flu_like codes are treated as pandemic unknowns and included in pct_covid as (covid_like + flu_like) / denominator. After that date, only covid_like / denominator is used.
• The legacy pipeline caps how far back each run generates estimates via n_backfill_days (run.py:85, default n_backfill_days = 70 in params.json.template:9), so each run only overwrites the most recent 70 days of estimates. This is visible in the coverage plot (Section 1.3.3) as a hard upper bound on the span of versions retained per (geo_value, time_value) observation. The new pipeline does not apply such a cap.
• The old pipeline seems to map to HRRs, rather than extracting them, via GeoMapper (geo_maps.py:107). In the new pipeline, both county and HRR exist in the raw source, according to constants.py:45, and they’re used directly at those levels here claims_task_groups.py:81.

---

1.2 Data Loading

old_dir <- here::here("dataset_comparison", "data", geo, "old")
new_dir <- here::here("dataset_comparison", "data", geo, "new")

df_old_smoothed <- load_archive(fs::path(old_dir, "monthly", glue("old_{indicator_meta$slug}")))
df_new <- load_archive(fs::path(new_dir, glue("new_{indicator_meta$slug}_{geo}_smoothed.parquet")), recurse = T)

v_old <- df_old_smoothed |> summarise(v = max(version, na.rm = TRUE)) |> collect() |> pull(v) |> as.Date()
v_new <- df_new         |> summarise(v = max(version, na.rm = TRUE)) |> collect() |> pull(v) |> as.Date()
v_shared <- min(v_old, v_new)

# Snapshot datasets (as of latest shared version)
old_snap <- summarize_revisions_explicit(df_old_smoothed, version_limit = v_shared) |>
  mutate(geo_value = as.character(geo_value))
new_snap <- summarize_revisions_explicit(df_new, version_limit = v_shared) |>
  mutate(geo_value = as.character(geo_value))

# Smoothed new pipeline snapshot 
new_snap_smoothed <- summarize_revisions_explicit(
  df_new |> mutate(value = coalesce(value_smoothed, value)),
  version_limit = v_shared
) |> mutate(geo_value = as.character(geo_value))

---

1.3 Dataset Overview

1.3.1 Time Range and Observation Counts

arch_old <- df_old_smoothed |>
  summarise(
    n = n(),
    n_geos = n_distinct(geo_value),
    date_min = min(time_value, na.rm = TRUE),
    date_max = max(time_value, na.rm = TRUE)
  ) |>
  collect() |>
  mutate(across(c(date_min, date_max), as.Date))

arch_new <- df_new |>
  summarise(
    n = n(),
    n_geos = n_distinct(geo_value),
    date_min = min(time_value, na.rm = TRUE),
    date_max = max(time_value, na.rm = TRUE)
  ) |>
  collect() |>
  mutate(across(c(date_min, date_max), as.Date))

tibble(
  View = c("Archive", glue("Snapshot (≤ {v_shared})")),
  `Old obs` = c(arch_old$n,                          nrow(old_snap)),
  `New obs` = c(arch_new$n,                          nrow(new_snap)),
  `Old geos` = c(arch_old$n_geos,                     n_distinct(old_snap$geo_value)),
  `New geos` = c(arch_new$n_geos,                     n_distinct(new_snap$geo_value)),
  `Old date range` = c(
    glue("{arch_old$date_min} – {arch_old$date_max}"),
    glue("{min(old_snap$time_value)} – {max(old_snap$time_value)}")
  ),
  `New date range` = c(
    glue("{arch_new$date_min} – {arch_new$date_max}"),
    glue("{min(new_snap$time_value)} – {max(new_snap$time_value)}")
  )
) |>
  knitr::kable(format.args = list(big.mark = ","))

View	Old obs	New obs	Old geos	New geos	Old date range	New date range
Archive	38,732,454	192,233,523	306	304	2020-02-01 – 2026-04-23	2019-10-08 – 2026-04-12
Snapshot (≤ 2026-04-12)	690,843	722,627	306	304	2020-02-01 – 2026-04-07	2019-10-08 – 2026-04-11

1.3.2 Geographic Coverage Over Time

This section shows how many geo units are reporting in each pipeline on each time_value date, using the snapshot as of 2026-04-12. Gaps or divergences indicate dates where one pipeline censors more aggressively.

bind_rows(
  old_snap |> count(time_value) |> mutate(pipeline = "Old"),
  new_snap |> count(time_value) |> mutate(pipeline = "New")
) |>
  ggplot(aes(x = time_value, y = n, color = pipeline)) +
  geom_line(alpha = 0.8) +
  scale_color_manual(values = dis_pal_list[c(2, 1)]) +
  labs(
    title    = glue("{indicator_meta$display_name}: {geo} units reporting per date"),
    subtitle = "Gaps = dates where that pipeline has no data",
    x = NULL, y = glue("# {geo}s with a value"), color = NULL
  )

geos_old <- old_snap |> pull(geo_value) |> unique() |> as.character()
geos_new <- new_snap |> pull(geo_value) |> unique() |> as.character()

only_old <- vctrs::vec_set_difference(geos_old, geos_new)
only_new <- vctrs::vec_set_difference(geos_new, geos_old)

cat(glue::glue(
  "- {geo}s in old: {length(geos_old)} | in new: {length(geos_new)}\n",
  "- Only in old ({length(only_old)}): {if (length(only_old)) paste(only_old, collapse = ', ') else '—'}\n",
  "- Only in new ({length(only_new)}): {if (length(only_new)) paste(only_new, collapse = ', ') else '—'}\n"
))

• hrrs in old: 306 | in new: 304
• Only in old (2): 351, 417
• Only in new (0): —

1.3.3 Coverage: Locations × Versions Over Time

This section shows the total number of (geo, date) pairs published in each version. A growing trend reflects normal backfill accumulation; sudden drops flag versions where a pipeline lost coverage.

cov_old <- df_old_smoothed |>
  group_by(version) |>
  summarise(n_obs = n(), .groups = "drop") |>
  collect() |>
  mutate(version = as.Date(version), pipeline = "Old")

cov_new <- df_new |>
  group_by(version) |>
  summarise(n_obs = n(), .groups = "drop") |>
  collect() |>
  mutate(version = as.Date(version), pipeline = "New")

cov <- bind_rows(cov_old, cov_new)

cov |>
  group_by(pipeline) |>
  arrange(version) |>
  mutate(n_obs_ma = slider::slide_dbl(n_obs, mean, .before = 30, .complete = FALSE)) |>
  ungroup() |>
  ggplot(aes(x = version, color = pipeline)) +
  geom_line(aes(y = n_obs), alpha = 0.25) +
  geom_line(aes(y = n_obs_ma)) +
  scale_color_manual(values = dis_pal_list[c(2, 1)]) +
  scale_y_continuous(labels = scales::comma) +
  labs(
    title    = glue("{indicator_meta$display_name} ({geo}): coverage per version"),
    subtitle = "Total (geo × date) published in each version; bold line (30-version rolling mean)",
    x = NULL, y = "# geo × date observations", color = NULL
  )

---

1.4 Within-Pipeline Revision Behavior

Before comparing the two pipelines, we look at how much each pipeline revises its own estimates.

bind_rows(
  old_snap |> mutate(pipeline = "Old"),
  new_snap |> mutate(pipeline = "New")
) |>
  mutate(
    revision_pct = (val_latest - val_first) / (abs(val_first) + 1e-6)
  ) |>
  filter(is.finite(revision_pct), abs(revision_pct) < 10) |>
  ggplot(aes(x = revision_pct, fill = pipeline)) +
  geom_histogram(binwidth = 0.05, alpha = 0.7, position = "identity") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_x_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = dis_pal_list[c(2, 1)]) +
  labs(
    title    = "Within-pipeline revision: (latest − first) / |first|",
    x = "Revision magnitude", y = "Count", fill = NULL
  )

bind_rows(
  old_snap |> mutate(pipeline = "Old"),
  new_snap |> mutate(pipeline = "New")
) |>
  mutate(rev_pct = (val_latest - val_first) / (abs(val_first) + 1e-6)) |>
  filter(is.finite(rev_pct)) |>
  group_by(pipeline) |>
  summarise(
    median_abs_rev  = median(abs(rev_pct), na.rm = TRUE),
    pct_unrevised   = mean(rev_pct == 0, na.rm = TRUE),
    pct_revised_gt5 = mean(abs(rev_pct) > 0.05, na.rm = TRUE),
    n               = n(),
    .groups = "drop"
  ) |>
  knitr::kable(digits = 3,
               col.names = c("Pipeline", "Median |revision|",
                              "% unrevised", "% revised >5%", "N"))

Pipeline	Median |revision|	% unrevised	% revised >5%	N
New	0.428	0.044	0.885	722627
Old	0.330	0.019	0.896	690843

---

1.5 Initial Comparison

This section provides an initial direct comparison between the two datasets to guide the next steps needed to identify the sources of differences.

Before reviewing the results, we define the error metrics used to compare the new pipeline values (\(y_i\)) against the old pipeline values (\(x_i\)) across \(n\) observations.

The Mean Absolute Error (MAE), sometimes referred to as Mean Absolute Difference (MAD), measures the average magnitude of the differences between the two pipelines. \[ \text{MAE} = \frac{1}{n} \sum_{i=1}^n |y_i - x_i| \]

The Root Mean Square Error (RMSE) squares these differences before averaging them. This approach places a higher penalty on large discrepancies compared to the MAE. \[ \text{RMSE} = \sqrt{\frac{1}{n} \sum_{i=1}^n (y_i - x_i)^2} \]

The Mean Absolute Scaled Error (MASE) scales the mean absolute difference by the average day-to-day variation of the old pipeline (\(x_i\)). Using the old pipeline as the scale makes the metric interpretable as how large is the inter-pipeline difference relative to how much the old pipeline itself typically moves between consecutive days. \[ \text{MASE} = \frac{\frac{1}{n} \sum_{i=1}^n |y_i - x_i|}{\frac{1}{n-1} \sum_{i=2}^n |x_i - x_{i-1}|} \] where \(y_i\) is the new pipeline value and \(x_i\) is the old pipeline value at time \(i\).

1.5.1 Shared Snapshot (as of 2026-04-12)

This comparison uses the state of both signals as they were known at the latest common issue date (2026-04-12) to ensure a fair comparison.

joined_snap <- calculate_differences(old_snap, new_snap)
metrics_snap <- get_summary_metrics(joined_snap)

1.5.1.1 Summary Metrics (Snapshot)

metrics_snap |>
  knitr::kable(digits = 3,
               col.names = c("Pearson r", "Spearman ρ", "MAD", "RMSE", "MASE", "N"))

Pearson r	Spearman ρ	MAD	RMSE	MASE	N
0.804	0.911	0.998	2.35	4.682	686189

1.5.1.2 Snapshot Scatter Plot w/ IQR-Scaled Difference

IQR_wrapper <- function(x) IQR(x, na.rm = TRUE)

joined_snap_iqr <- joined_snap |>
  group_by(geo_value) |>
  arrange(time_value) |>
  mutate(
    IQR_val = slider::slide_dbl(val_latest_old, IQR_wrapper, .before = 59, .complete = FALSE)
  ) |>
  tidyr::fill(IQR_val, .direction = "downup") |>
  mutate(
    IQR_backup = suppressWarnings(min(IQR_val[IQR_val > 0], na.rm = TRUE)),
    IQR_val    = dplyr::if_else(IQR_val == 0 | is.na(IQR_val), IQR_backup, IQR_val),
    abs_diff_scaled = abs_diff / IQR_val
  ) |>
  ungroup()

The scatter shows where the two pipelines agree or diverge. To gives sense of whether a given discrepancy is meaningful relative to the typical variability of that indicator, the plot addresses that by colouring each point by |New − Old| / IQR, where the IQR is computed as a 60-day rolling interquartile range of the old pipeline’s values for that geography.

breaks_iqr <- round(
  expm1(seq(
    log1p(min(joined_snap_iqr$abs_diff_scaled[joined_snap_iqr$abs_diff_scaled > 0], na.rm = TRUE)),
    log1p(max(joined_snap_iqr$abs_diff_scaled, na.rm = TRUE)),
    length.out = 4
  )), 1)

joined_snap_iqr |>
  filter(abs_diff_scaled > 0) |>
  ggplot(aes(x = val_latest_old, y = val_latest_new, color = abs_diff_scaled)) +
  geom_point(size = 0.8, alpha = 0.1) +
  geom_point(data = ~ filter(.x, abs_diff_scaled > 50), size = 0.8, alpha = 0.5) +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "black") +
  scale_color_gradient(
    low    = "gray90",
    high   = dis_pal_list[1],
    name   = "|Diff| / IQR",
    trans  = "log1p",
    breaks = breaks_iqr
  ) +
  scale_x_log1p() +
  scale_y_log1p() +
  labs(
    title    = glue("{indicator_meta$display_name} ({geo}): Snapshot View as of {v_shared}"),
    subtitle = "|New − Old| / rolling 60-day IQR of Old pipeline",
    x        = "Old Pipeline (log1p scale)",
    y        = "New Pipeline (log1p scale)"
  )

1.5.1.3 Top 1,000 Observations by IQR-Scaled Difference

joined_snap_iqr |>
  select(
geo_value, time_value,
    val_old = val_latest_old,
    val_new = val_latest_new,
    abs_diff,
    IQR_val,
    abs_diff_scaled
  ) |>
    mutate(across(where(is.double), \(x) round(x, 4))) |>
  slice_max(abs_diff_scaled, n = 1000) |>
  DT::datatable(
    caption  = glue("Geo-level: {geo} · Snapshot as of {v_shared}"),
    rownames = FALSE,
    colnames = c("Geo", "Date", "Old", "New", "Diff", "IQR", "|Diff|/IQR"),
    options  = list(scrollX = TRUE, pageLength = 20)
  )

NoteRelevant Observations

• From the previous plots, it can be observed that most values do not align with the one-to-one diagonal line and tend to be larger for the old pipeline. As noted by David, the old pipeline data “should be ~1/7 as large as the new pipeline”. Consequently, the next steps will be performed using a smoothed version of the new pipeline data.
• An initial comparison using state-level data showed smaller differences. Here, MAD and RMSE are 0.998 and 2.35, respectively, compared to 1.203 and 2.736 at the state level. This improvement likely results from the new pipeline’s direct use of the identified HRR regions. The next step is to perform the analysis at the county level to eliminate differences caused by the Geomapper used in the previous pipeline (this is not true for hospital admissions).
• Examining the IQR-scaled differences closely reveals that a handful of observations have very high values. These are primarily from 2020, a period in which the rolling 60-day IQR was nearly zero. To mitigate this, future analysis will begin in May, 2020 as discussed with Nat.
• As noted in Section 1.3.3, the legacy pipeline caps revision history at n_backfill_days = 70. Subsequent sections that analyze revision behavior (archive view, pre/post comparison) are therefore limited to this window.

---

1.6 Shared Snapshot (may 2020 onward w/ Smoothed New Pipeline)

snap_start <- as.Date("2020-05-01")
joined_snap_sm <- calculate_differences(
  old_snap,
  new_snap_smoothed
)
metrics_snap_may  <- get_summary_metrics(joined_snap_sm |> 
  filter(time_value >= snap_start))

joined_snap_sm_iqr <- joined_snap_sm |>
  group_by(geo_value) |>
  arrange(time_value) |>
  mutate(
    IQR_val = slider::slide_dbl(val_latest_old, IQR_wrapper, .before = 59, .complete = FALSE)
  ) |>
  tidyr::fill(IQR_val, .direction = "downup") |>
  mutate(
    IQR_backup = suppressWarnings(min(IQR_val[IQR_val > 0], na.rm = TRUE)),
    IQR_val = dplyr::if_else(IQR_val == 0 | is.na(IQR_val), IQR_backup, IQR_val),
    abs_diff_scaled = abs(abs_diff) / IQR_val
  ) |>
    filter(time_value >= snap_start) |> 
  ungroup()

1.6.1 Gaussian Smoother Effect

Before comparing against the old pipeline, this section shows how the Gaussian smoother changes the new pipeline’s own values. Each point is a (geo, date) pair. The x-axis is the raw value, the y-axis the smoothed value.

smoother_effect <- new_snap |>
  filter(time_value >= snap_start) |>
  inner_join(
    new_snap_smoothed |>
      filter(time_value >= snap_start) |>
      select(geo_value, time_value, val_latest_smoothed = val_latest),
    by = c("geo_value", "time_value")
  ) |>
  mutate(smooth_diff = val_latest_smoothed - val_latest)

smoother_effect |>
  ggplot(aes(x = val_latest, y = val_latest_smoothed)) +
  geom_point(alpha = 0.05, size = 0.7, color = dis_pal_list[2]) +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "red") +
  scale_x_log1p() +
  scale_y_log1p() +
  labs(
    title = glue("{indicator_meta$display_name} ({geo}): Smoother Effect on New Pipeline"),
    x = "New Pipeline (raw value, log1p scale)",
    y = "New Pipeline (smoothed value, log1p scale)"
  )

# Pick geos with the largest median absolute smoother adjustment
sample_geos <- smoother_effect |>
  group_by(geo_value) |>
  summarise(med_adj = median(abs(smooth_diff), na.rm = TRUE), .groups = "drop") |>
  slice_max(med_adj, n = 6) |>
  pull(geo_value)

smoother_effect |>
  filter(geo_value %in% sample_geos) |>
  left_join(
old_snap |> filter(time_value >= snap_start) |> select(geo_value, time_value, val_old = val_latest),
    by = c("geo_value", "time_value")
  ) |>
    select(geo_value, time_value, `New (raw)` = val_latest, `New (smoothed)` = val_latest_smoothed, Old = val_old) |>
  tidyr::pivot_longer(c(`New (raw)`, `New (smoothed)`, Old), names_to = "series", values_to = "value") |>
  mutate(series = factor(series, levels = c("Old", "New (raw)", "New (smoothed)"))) |>
  ggplot(aes(x = time_value, y = value, color = series)) +
  geom_line(alpha = 0.8) +
  scale_color_manual(values = c(
    "Old"            = dis_pal_list[1],
    "New (raw)"      = dis_pal_list[2],
    "New (smoothed)" = dis_pal_list[3]
  )) +
  facet_wrap(~ geo_value, scales = "free_y", ncol = 2) +
  labs(
    title    = glue("{indicator_meta$display_name} ({geo}): Top 6 most affected geos"),
    subtitle = "Geos selected by largest median",
    x        = NULL, y = "Value", color = NULL
  )

1.6.2 Summary Metrics

metrics_snap_may |>
  knitr::kable(digits = 3,
               col.names = c("Pearson r", "Spearman ρ", "MAD", "RMSE", "MASE", "N"))

Pearson r	Spearman ρ	MAD	RMSE	MASE	N
0.798	0.921	1.018	2.385	4.687	658829

1.6.3 Snapshot Scatter Plot w/ IQR-Scaled Difference

breaks <- round(
expm1(seq(log1p(min(joined_snap_sm_iqr$abs_diff_scaled[joined_snap_sm_iqr$abs_diff_scaled >0])), 
log1p(max(joined_snap_sm_iqr$abs_diff_scaled)), length.out = 4)), 
1)

joined_snap_sm_iqr |>
  filter(abs_diff_scaled > 0) |> 
  ggplot(aes(x = val_latest_old, y = val_latest_new, color = abs_diff_scaled)) +
  geom_point(size = 0.8, alpha = 0.1) +
  geom_point(data = ~ filter(.x, abs_diff_scaled > 50), size = 0.8, alpha = 0.5) +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "black") +
  scale_color_gradient(
    low = "gray90",
    high = dis_pal_list[1],
    name = "|Diff| / IQR",
    trans = "log1p",
    breaks = breaks
  ) +
  scale_x_log1p() +
  scale_y_log1p() +
  labs(
    title = glue("{indicator_meta$display_name} ({geo}): Snapshot as of {v_shared}"),
    subtitle = "|New (smoothed) − Old| / rolling 60-day IQR of Old",
    x = "Old Pipeline (log1p scale)",
    y = "New Pipeline (smoothed, log1p scale)"
  )

1.6.4 Top Observations by IQR-Scaled Difference

joined_snap_sm_iqr |>
  select(
    geo_value, time_value,
    val_old = val_latest_old,
    val_new = val_latest_new,
    abs_diff, IQR_val, abs_diff_scaled
  ) |>
  mutate(across(where(is.double), \(x) round(x, 4))) |>
  arrange(desc(abs_diff_scaled)) |>
  slice_head(n = 1000) |>
  DT::datatable(
    caption  = glue("Geo-level: {geo} · Snapshot as of {v_shared} · may 2020 onward · smoothed new pipeline"),
    rownames = FALSE,
    colnames = c("Geo", "Date", "Old", "New (smoothed)", "Diff", "IQR", "|Diff|/IQR"),
    options  = list(scrollX = TRUE, pageLength = 20)
  )

---

1.6.5 Full Archive Overlap

This section compares the two pipelines across their shared historical overlap from May 2020 onward. We join the datasets on geo_value, time_value, and version, using the Gaussian-smoothed values from the new pipeline. This shows how closely the pipelines agreed at every point in time they both published.

Note that this version is not scaled by the IQR due to computational expense.

archive_start <- as.Date("2020-05-01")

# Prepare each Arrow dataset with renamed columns before joining
ds_old_prep <- df_old_smoothed |>
  filter(time_value >= archive_start) |>
  mutate(geo_value = cast(geo_value, utf8()),
         version   = cast(version, date32())) |>
  select(geo_value, time_value, version, value) |>
  rename(value_old = value)

# cast(version, date32()) collapses multiple dttm values to the same date,
# so deduplicate after collecting to avoid many-to-many join expansion.
ds_new_prep <- df_new |>
  filter(time_value >= archive_start) |>
  mutate(
    geo_value = cast(geo_value, utf8()),
    version   = cast(version, date32()),
    value_new = coalesce(value_smoothed, value)
  ) |>
  select(geo_value, time_value, version, value_new) |>
  collect() |>
  mutate(across(c(time_value, version), as.Date)) |>
  distinct(geo_value, time_value, version, .keep_all = TRUE)

ds_old_prep_collected <- ds_old_prep |>
  collect() |>
  mutate(across(c(time_value, version), as.Date))

# Join in R after deduplication
archive_joined <- ds_old_prep_collected |>
  inner_join(ds_new_prep, by = c("geo_value", "time_value", "version")) |>
  mutate(
    val_latest_old = value_old,
    val_latest_new = value_new,
    abs_diff = val_latest_new - val_latest_old,
    rel_diff = abs_diff / (abs(val_latest_old) + 1e-6),
    diff     = abs_diff
  ) |>
  select(geo_value, time_value, version,
         val_latest_old, val_latest_new, abs_diff, rel_diff, diff)

# Free the Arrow dataset handles
rm(ds_old_prep, ds_old_prep_collected, ds_new_prep, df_old_smoothed, df_new)
invisible(gc())

archive_metrics <- get_summary_metrics(archive_joined)

# Compute rolling IQR on a stratified sample to avoid memory problems
sampled_versions <- archive_joined |>
  distinct(geo_value, version) |>
  group_by(geo_value) |>
  mutate(row_id = sample(n())) |>
  filter(row_id <= 20L) |>
  ungroup() |>
  select(-row_id)

archive_iqr <- archive_joined |>
  semi_join(sampled_versions, by = c("geo_value", "version")) |>
  mutate(abs_diff = abs(diff)) |>
  group_by(geo_value, version) |>
  arrange(time_value) |>
  mutate(
    IQR_val = slider::slide_dbl(val_latest_old, IQR_wrapper, .before = 59, .complete = FALSE)
  ) |>
  tidyr::fill(IQR_val, .direction = "downup") |>
  mutate(
    IQR_backup      = suppressWarnings(min(IQR_val[IQR_val > 0], na.rm = TRUE)),
    IQR_val         = dplyr::if_else(IQR_val == 0 | is.na(IQR_val), IQR_backup, IQR_val),
    abs_diff_scaled = abs_diff / IQR_val
  ) |>
  ungroup()

rm(sampled_versions)
invisible(gc())

1.6.5.1 Summary Metrics (Full Archive)

archive_metrics |>
  knitr::kable(digits = 3,
               col.names = c("Pearson r", "Spearman ρ", "MAD", "RMSE", "MASE", "N"))

Pearson r	Spearman ρ	MAD	RMSE	MASE	N
0.796	0.914	1.098	2.509	4.746	32766160

1.6.5.2 Archive Scatter Plot

archive_joined |>
  # Since we are only looking at the shape of the difference here, 
  # a sample should be enough.
  slice_sample(n = min(5e5, nrow(archive_joined))) |> 
  ggplot(aes(x = val_latest_old, y = val_latest_new)) +
  geom_point(alpha = 0.05, color = dis_pal_list[2]) +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "red") +
  scale_x_log1p() +
  scale_y_log1p() +
  labs(
    title    = glue("{indicator_meta$display_name} ({geo}): Full Archive Historical Agreement"),
    subtitle = glue("Pearson r = {round(archive_metrics$pearson, 3)}, N = {scales::comma(archive_metrics$n)} | New pipeline smoothed, from {archive_start}"),
    x = "Old Pipeline (log1p scale)",
    y = "New Pipeline (smoothed, log1p scale)"
  )

NoteRelevant Observations

• We can see that smoothing the new pipeline did not reduced the difference between the pipelines. There is still a difference from the way that values are procesed since the new pipeline applies a 7-day triling sum before division, but it can be observed from the plots with the locations with bigger median absolute difference that values are in different levels. So even applying a mooving average to the new pipeline would not explain the difference.
• From the plots for locations with larger median absolute differences, it is clear that the differences are concentrated at the beginning, where the old pipeline produces clearly larger values. This can be explained by the change in methodology in the old pipeline, in which unknowns were not incorporated into the new pipeline after June 2021, and by the removal of mixed cases. To test this, the next section splits the data before and after June 2021.
• Another relevant observation is that Jeffreys prior can inflate the values for the old pipeline, but the effect is probably small given the \(N \geq 100\) denominator floor.

1.7 Pre- vs Post-2021-06 Comparison

1.7.1 Snapshot

if (!is.null(joined_snap_sm_iqr)) {
  add_era_naive(joined_snap_sm_iqr |> filter(time_value >= snap_start), "geo_value") |>
    era_summarise() |>
    knitr::kable(digits = 3, col.names = era_col_names)
}

Era	Pearson r	Spearman ρ	MAD	Bias (new−old)	MASE	N
post-2021-06	0.877	0.902	0.454	-0.390	3.264	538445
pre-2021-06	0.706	0.807	3.540	-3.504	6.242	120384

1.7.2 Archive

add_era_naive(archive_iqr, c("geo_value", "version")) |>
  era_summarise() |>
  knitr::kable(digits = 3, col.names = era_col_names)

Era	Pearson r	Spearman ρ	MAD	Bias (new−old)	MASE	N
post-2021-06	0.890	0.901	0.458	-0.403	3.127	297119
pre-2021-06	0.729	0.775	3.733	-3.725	6.322	74251

1.7.3 Mean IQR-Scaled Absolute Difference Over Time

ts_data <- bind_rows(
  if (!is.null(joined_snap_sm_iqr))
    joined_snap_sm_iqr |> filter(!is.na(abs_diff_scaled), is.finite(abs_diff_scaled)) |>
      mutate(view = "Snapshot"),
  if (!is.null(archive_iqr))
    archive_iqr |> filter(!is.na(abs_diff_scaled), is.finite(abs_diff_scaled)) |>
      mutate(view = "Archive")
)

if (nrow(ts_data) > 0) {
  ts_data |>
    group_by(view, time_value) |>
    summarise(mean_abs_diff_scaled = mean(abs_diff_scaled, na.rm = TRUE), .groups = "drop") |>
    ggplot(aes(x = time_value, y = mean_abs_diff_scaled)) +
    geom_line(alpha = 0.5, color = dis_pal_list[1]) +
    geom_smooth(method = "loess", span = 0.1, se = FALSE,
                color = dis_pal_list[2], linewidth = 1.2) +
    geom_hline(yintercept = 0, linetype = "dashed") +
    geom_vline(xintercept = as.Date("2021-06-01"),
               linetype = "dotted", color = "gray40") +
    annotate("text", x = as.Date("2021-06-01"), y = Inf,
             label = "COVID def. switch", vjust = 2, hjust = -0.05,
             size = 3, color = "gray40") +
    facet_wrap(~view, ncol = 1, scales = "free_y") +
    labs(
      title    = glue("{indicator_meta$display_name} ({geo}): Mean IQR-Scaled Absolute Difference Over Time"),
      subtitle = "Cross-geography mean",
      x = NULL, y = "Mean |new − old| / IQR(old)"
    )
}

1.8 Geographic Variation

1.8.1 Snapshot

if (!is.null(joined_snap_sm_iqr)) {
  joined_snap_sm_iqr |>
    group_by(geo_value) |>
    arrange(time_value) |>
    mutate(naive_diff = val_latest_old - lag(val_latest_old)) |>
    summarise(
      mase = mean(abs_diff, na.rm = TRUE) / mean(abs(naive_diff), na.rm = TRUE),
      bias = mean(rel_diff, na.rm = TRUE),
      .groups = "drop"
    ) |>
    geo_mase_plot("Snapshot view")
}

1.8.2 Archive

archive_iqr |>
  group_by(geo_value, version) |>
  arrange(time_value) |>
  mutate(naive_diff = val_latest_old - lag(val_latest_old)) |>
  group_by(geo_value) |>
  summarise(
    mase = mean(abs_diff, na.rm = TRUE) / mean(abs(naive_diff), na.rm = TRUE),
    bias = mean(rel_diff, na.rm = TRUE),
    .groups = "drop"
  ) |>
  geo_mase_plot("Archive view")

---

1.9 Final Notes

NoteRelevant Observations

• Apparently, the most consequential methodological difference is not smoothing but the COVID definition. The legacy pipeline always includes Mixed + Covid_like + Flu_like, whereas the new pipeline uses covid_like + flu_like before June 2021 and covid_like alone afterward. The pre/post split shows larger errors in the pre-2021-06 period for both snapshot and archive views.
• The old pipeline applies a one-sided Gaussian filter to the computed rate. The new pipeline uses a 7-day trailing sum on raw counts. Applying a Gaussian filter to the new pipeline’s output does not eliminate the level difference.
• At the HRR level, MAD and RMSE are lower than at state level. This is consistent with the new pipeline using HRR geographies directly rather than GeoMapper-aggregated HRRs.
• At the minimum denominator floor (N = 100), the Jeffreys prior adds at most ~0.5 pp upward bias. This effect is visible mainly in sparse geographies near the censoring threshold and does not explain the bulk of the difference.
• The archive-view MASE and IQR-scaled metrics align with the snapshot view, indicating that disagreements are stable over time and not driven by revisions.