Revision Analysis: Hospital Admissions: Smoothed COVID-19

This notebook analyses the revision behaviour of a single time series stored as an epi_archive. The central question is:

How long, on average, does a series continue to receive revisions after its first publication?

---

1 Data Loading

archive <- NULL

if (!is.null(params$input_dir)) {
    path <- here::here(params$input_dir)

    if (fs::is_dir(path)) {
        cli::cli_inform("Reading data from directory: {.path {params$input_dir}}")
        raw <- arrow::open_dataset(path) %>%
            select(any_of(c("geo_value", "time_value", "version", "issue", "value"))) %>%
            filter(
                time_value >= as.Date(params$start_day),
                time_value <= as.Date(params$end_day)
            ) %>%
            collect()
    } else if (grepl("\\.parquet$", path)) {
        cli::cli_inform("Reading Parquet file: {.path {params$input_dir}}")
        raw <- arrow::read_parquet(path)
    } else {
        cli::cli_inform("Reading CSV file: {.path {params$input_dir}}")
        raw <- read_csv(path)
    }

    raw <- raw %>% mutate(time_value = as.Date(time_value))

    # Accept either 'version' or 'issue' as the version column.
    if ("issue" %in% names(raw) && !"version" %in% names(raw)) {
        raw <- rename(raw, version = issue)
    }

    if (!"version" %in% names(raw)) {
        cli::cli_abort(
            "The input data must contain a 'version' or 'issue' column. ",
            "See the README for the required schema."
        )
    }

    # If it was a single file (CSV/Parquet), we still need to filter/mutate metadata
    if (!fs::is_dir(path)) {
        raw <- raw %>%
            mutate(version = as.Date(version)) %>%
            filter(
                time_value >= as.Date(params$start_day),
                time_value <= as.Date(params$end_day)
            )
    }

    archive <- raw %>%
        select(geo_value, time_value, version, value) %>%
        as_epi_archive()

    msg_verb <- "Loaded"
} else {
    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)
                params$signal %in% meta_res[[params$source]]$signals
            },
            error = function(e) {
                FALSE
            }
        )
        if (meta_available) {
            use_new_api <- TRUE
        }
    }

    if (use_new_api) {
        df_raw <- epidata_archive(
            source = params$source,
            signals = params$signal,
            geo_type = params$geo_type
        ) %>%
            rename(any_of(c(time_value = "reference_time", version = "report_time"))) %>%
            filter(
                time_value >= as.Date(params$start_day),
                time_value <= as.Date(params$end_day)
            )
    } else {
        df_raw <- pub_covidcast(
            source = params$source,
            signal = params$signal,
            geo_type = params$geo_type,
            time_type = params$time_type,
            time_values = epirange(params$start_day, params$end_day),
            issues = "*"
        )
    }

    if ("issue" %in% colnames(df_raw) && !"version" %in% colnames(df_raw)) {
        df_raw <- rename(df_raw, version = issue)
    }

    if (nrow(df_raw) == 0 || !"version" %in% colnames(df_raw)) {
        cli::cli_abort(
            "No versioned data returned from the API. ",
            "The signal may not support revision history."
        )
    }

    archive <- df_raw %>%
        select(geo_value, time_value, version, value) %>%
        as_epi_archive()

    msg_verb <- "Fetched"
}

#> Reading data from directory: 'revision_analysis/data/revisions'

# Metadata calculation
n_rows <- nrow(archive$DT)
n_locations <- n_distinct(archive$DT$geo_value)
n_versions <- n_distinct(archive$DT$version)
version_start <- min(archive$DT$version)
version_end <- max(archive$DT$version)

# Data loading summary message
cli::cli_inform(c("i" = "{msg_verb} {label_comma()(n_rows)} rows | {label_comma()(n_locations)} locations | {label_comma()(n_versions)} versions ({version_start} – {version_end})"))

#> ℹ Loaded 11,483,645 rows | 51 locations | 1,815 versions (2020-06-21 –
#>   2026-04-02)

# Reclaim memory
invisible(gc())

---

2 Revision count overview

The plot below shows how many revisions were published for each time value, with the distribution across locations shown as a boxplot.

---

3 Overall revision summary

revision_analysis() provides a concise summary of how a signal matures. The results describe the spread of revisions, which measures the difference between the minimum and maximum values ever published for an observation. Scaling this by the signal’s value gives the relative spread. This shows the magnitude of revisions as a percentage. The report also presents convergence through the time it takes for data to reach within 20% of its final value. This convergence metric identifies the median duration before an observation enters and stays in its final neighborhood.

rev_result <- epiprocess::revision_analysis(
    archive,
    value,
    min_waiting_period = 0
)
rev_beh <- rev_result$revision_behavior

# Track which (geo_value, time_value) pairs were included vs. excluded.
# revision_analysis() drops pairs that have no version past the waiting
# period — typically the most recent time values with only an initial report.
all_pairs <- archive$DT %>% distinct(geo_value, time_value)
analyzed_pairs <- rev_beh %>% distinct(geo_value, time_value)
n_excluded <- nrow(anti_join(all_pairs, analyzed_pairs, by = c("geo_value", "time_value")))
n_analyzed <- nrow(analyzed_pairs)

# Global metrics for follow-on plots/narrative
max_lag_days <- as_days(rev_beh$max_lag)
no_rev_pct <- mean(rev_beh$n_revisions == 0, na.rm = TRUE)
conv_days_all <- as_days(rev_beh$lag_near_latest)
invisible(gc())

# Pre-calculate relative differences to final for all observations.
# We also find the first release version for age calculation.
convergence_df <- archive$DT %>%
    semi_join(analyzed_pairs, by = c("geo_value", "time_value")) %>%
    group_by(geo_value, time_value) %>%
    mutate(
        version_first = min(version, na.rm = TRUE),
        value_latest = value[which.max(version)]
    ) %>%
    ungroup() %>%
    filter(value_latest != 0, !is.na(value_latest), !is.na(value)) %>%
    mutate(
        days_since_first = as.numeric(version - version_first, units = "days"),
        abs_rel_diff = abs(value - value_latest) / abs(value_latest)
    )
invisible(gc())

Important

revision_analysis() requires at least one version past the waiting period to assess revision behavior. 16060 time×location pair(s) were excluded because they only have an initial report with no subsequent versions yet. The 54543 remaining pairs are used for all analyses below.

• Time range: 2020-02-01 to 2023-12-31
• Min lag (days to first version): Min.: 3, 1st Qu.: 3, Median: 4, Mean: 3.7, 3rd Qu.: 4, Max.: 5
• No revisions: 0%
• Quick revisions (≤ 7 days): 0%
• Few revisions (≤ 3): 0%
• Minor relative revisions (< 0.1): 1%
• Major absolute revisions (> 1.961): 19%
• Days until within 20% of latest: Min.: 3, 1st Qu.: 6, Median: 18, Mean: 55.1, 3rd Qu.: 35, Max.: 1182

Here we collect a number of statistics for each location. This includes:

• Avg # revisions. The average number of times a value was updated after first release across all dates for that location.
• Min lag. The minimum number of days between the event date and the first publication.
• Max lag. The maximum number of days between the event date and the final revision.
• Avg spread. The average absolute difference between the minimum and maximum published values.
• Avg rel. spread. The average relative difference between the minimum and maximum published values.
• Med. rel. spread. The median relative spread across all observation dates.
• Med. window (d). The median number of days from first publication to the final revision.
• Avg conv. (d). The average days until the value stays within a certain percentage (default 20%) of the final value and never leaves it.

---

4 Reporting Latency

4.1 Issue date cadence

issue_dates <- sort(unique(archive$DT$version))
gaps <- as.numeric(diff(issue_dates), units = "days")

tibble(gap = gaps) %>%
    ggplot(aes(x = gap)) +
    geom_histogram(bins = 30, color = "white") +
    scale_x_continuous(breaks = breaks_pretty()) +
    labs(
        title = "Cadence of archive updates",
        subtitle = sprintf(
            "%d unique issue dates | median gap: %.0f d",
            n_versions, median(gaps)
        ),
        x = "Gap (days)", y = "Count"
    )

4.2 Reporting latency by geography

Each coloured line is a different issue date. Darker shades are older. Vertical spread at a single time_value indicates back-revision.

backfill_df <- archive$DT %>%
    filter(!is.na(value)) %>%
    group_by(geo_value, time_value) %>%
    mutate(
        backfill_gap = as.numeric(version - time_value, units = "days"),
        is_initial_release = version == min(version, na.rm = TRUE)
    ) %>%
    ungroup()

# Identify worst locations for deeper drill-down
worst_geos <- location_stats %>%
    arrange(desc(rel_spread)) %>%
    slice_head(n = params$n_worst %||% 12) %>%
    pull(geo_value)
n_worst_actual <- length(worst_geos)

backfill_geom <- backfill_df %>%
    filter(geo_value %in% worst_geos) %>%
    mutate(geo_value = factor(geo_value, levels = worst_geos))

ggplot(backfill_geom, aes(x = backfill_gap, fill = is_initial_release)) +
    geom_histogram(bins = 40, position = "identity", alpha = 0.7) +
    facet_wrap(~geo_value, ncol = 3, scales = "free_y") +
    scale_fill_manual(
        values = c("TRUE" = dis_pal[2], "FALSE" = dis_pal[1]),
        labels = c("TRUE" = "Initial Release", "FALSE" = "Subsequent Revision"),
        name = "Release type"
    ) +
    scale_x_continuous(breaks = breaks_pretty()) +
    labs(
        title = "Backfill gap distribution (worst locations)",
        subtitle = "Initial release latency (blue) vs. subsequent revision depth (red)",
        x = "Backfill gap (days)",
        y = "Count"
    ) +
    theme(legend.position = "top", strip.text.x = element_text(size = 8))

4.3 Visualizing version history

Each coloured line represents a different issue date. Vertical spread at a single time_value indicates the magnitude of the backfill corrections over time.

min_v <- min(issue_dates)
max_v <- max(issue_dates)

# selected_versions <- data.frame(version = issue_dates) %>%
#     mutate(age = as.numeric(version - min_v)) %>%
#     filter(
#         age <= 7 |
#             (age <= 21 & age %% 2 == 0) |
#             (age <= 365 & age %% 30 == 0) |
#             (age %% 60 == 0) |
#             version == max_v
#     ) %>%
#     pull(version)

# archive %>%
#     autoplot(value,
#         .versions = selected_versions,
#         .interactive = FALSE
#     ) +
#     scale_colour_viridis_c(
#         option = "plasma",
#         begin = 0, end = 0.8,
#         labels = scales::label_date(),
#         trans = "date"
#     ) +
#     theme(legend.position = "bottom")

---

5 Revision duration & maturation

This section assesses how long it takes for a signal to stabilize. We consider both the time until the final revision and the time until the value stays persistently within a given percentage of its eventually finalized value.

# Define thresholds for assessment (relative differences)
thresholds <- c(0.5, 0.2, 0.1, 0.05, 0.02)

# Calculate stabilization lags for each threshold
stability_multi <- map_dfr(thresholds, function(tau) {
    convergence_df %>%
        group_by(geo_value, time_value) %>%
        arrange(desc(version)) %>%
        summarize(
            # Find the lag of the last version that was ABOVE the threshold.
            # The next version is when it became "stable".
            lag_stable = {
                over <- which(abs_rel_diff > tau)
                if (length(over) == 0) {
                    0
                } else {
                    max(days_since_first[over])
                }
            },
            threshold = label_percent()(tau),
            threshold_num = tau,
            .groups = "drop"
        )
})

# Combine with the "final revision" window (max_lag from original behavior)
rev_window_df <- rev_beh %>%
    filter(n_revisions > 0) %>%
    transmute(
        geo_value, time_value,
        lag_stable = as_days(max_lag),
        threshold = "Final revision",
        threshold_num = 0
    )

stability_combined <- bind_rows(stability_multi, rev_window_df) %>%
    mutate(threshold = factor(threshold, levels = c("50%", "20%", "10%", "5%", "2%", "Final revision")))
invisible(gc())

5.1 Distribution of stabilization times

How long does it typically take to reach each stability threshold?

ggplot(stability_combined, aes(x = lag_stable, fill = threshold)) +
    geom_histogram(bins = 40, color = "white") +
    facet_wrap(~threshold, ncol = 1, scales = "free_y") +
    scale_x_continuous(breaks = breaks_pretty()) +
    labs(
        title = "Stabilization time distribution",
        subtitle = "Days from first publication to persistent stability",
        x = "Days since first release",
        y = "Number of observations",
        fill = "Threshold"
    )

5.2 Comparative ECDF of stabilization

By day X, what fraction of observations have settled within the specified threshold and never leave it?

ggplot(stability_combined, aes(x = lag_stable, color = threshold)) +
    stat_ecdf(linewidth = 0.8) +
    scale_x_continuous(breaks = breaks_pretty()) +
    scale_y_continuous(labels = label_percent()) +
    labs(
        title = "ECDF of stabilization time by threshold",
        subtitle = "Showing when observations reach and stay within X% of the final value",
        x = "Days from first publication to stabilization",
        y = "Cumulative fraction of observations",
        color = "Threshold"
    )

5.3 Stabilization summary

---

6 Bias

This section explores if revisions follow predictable trends based on timing, geography, or signal magnitude.

6.1 Censoring bias

The plot below shows how the speed of data stabilization changes over time. By grouping observations into time buckets, we can identify trends in reporting maturity for each threshold. The dashed line represents the censoring boundary, which marks the maximum possible lag available given the current date. As recent observations approach this boundary, their stabilization times are naturally truncated. This suggests that these data points may still be subject to further revisions.

max_v <- max(archive$DT$version)

# Determine best aggregation unit based on date range
date_range_days <- as.numeric(max(stability_combined$time_value) - min(stability_combined$time_value))
agg_unit <- case_when(
    date_range_days > 365 * 2 ~ "quarter",
    date_range_days > 365 ~ "month",
    date_range_days > 100 ~ "2-week",
    date_range_days > 40 ~ "week",
    date_range_days > 20 ~ "3-day",
    TRUE ~ "day"
)

# Plot stabilization lag over time across all thresholds
stability_combined %>%
    mutate(time_bucket = floor_date(time_value, unit = agg_unit)) %>%
    ggplot(aes(x = time_bucket, y = lag_stable, color = threshold, fill = threshold)) +
    stat_summary(
        fun.data = median_hilow, fun.args = list(conf.int = 0.5),
        geom = "ribbon", alpha = 0.2, color = NA
    ) +
    stat_summary(
        fun = median, geom = "line"
    ) +
    geom_function(
        fun = function(x) as.numeric(max_v - as.Date(x, origin = "1970-01-01")),
        linetype = "dashed", color = "gray50", linewidth = 0.5
    ) +
    annotate(
        "text",
        x = min(stability_combined$time_value),
        y = as_days(max_v - min(stability_combined$time_value)),
        label = "Censoring Boundary", color = "gray50",
        vjust = 1, hjust = -1, size = 3
    ) +
    scale_x_date(breaks = breaks_pretty()) +
    scale_y_continuous(breaks = breaks_pretty()) +
    labs(
        title = "Stabilization distribution over time",
        subtitle = sprintf(
            "Distribution of days to persistent stability grouped by %s
            Lines indicate the median, while ribbons show the IQR.",
            agg_unit
        ),
        x = "Observation date (aggregated)",
        y = "Days since first publication",
        color = "Threshold",
        fill = "Threshold"
    )

6.2 Signal magnitude bias

This section examines if the relative scale of revisions is influenced by the magnitude of the signal itself. This is visualized by comparing the maximum observed value against the relative spread (total revision range divided by the signal’s final value) on a log-log scale. This analysis helps determine if the signal becomes more volatile at specific scales, for instance, if low-value periods are subject to proportionately larger revisions than high-value periods.

rev_beh %>%
    filter(rel_spread > 0) %>%
    ggplot(aes(x = max_value, y = rel_spread)) +
    geom_point(alpha = 0.2, size = 1) +
    geom_smooth() +
    scale_x_log10(labels = label_comma()) +
    scale_y_log10(labels = label_percent()) +
    labs(
        title = "Signal magnitude vs. relative spread",
        subtitle = "Log-log scale",
        x = "Maximum observed value (log scale)", y = "Relative spread (log scale)"
    )

7 By Location Bias

7.1 Initial Bias

Here we characterise whether each location tends to under- or over-report in its first release.

For each observation we compute:

\[\text{initial bias} = \frac{v_\text{first} - v_\text{latest}}{v_\text{latest}}\]

A negative value means the first release was below the latest (under-reporting); a positive value means it was above (over-reporting).

# For each (geo_value, time_value) we need the *first* and *latest* value.
# revision_behavior already carries `spread` (max - min) and `rel_spread`.
# We reconstruct initial vs. latest from the raw archive DT.
first_vals <- archive$DT %>%
    group_by(geo_value, time_value) %>%
    slice_min(version, n = 1, with_ties = FALSE) %>%
    ungroup() %>%
    select(geo_value, time_value, value_first = value)

latest_vals <- archive$DT %>%
    group_by(geo_value, time_value) %>%
    slice_max(version, n = 1, with_ties = FALSE) %>%
    ungroup() %>%
    select(geo_value, time_value, value_latest = value)

bias_df <- first_vals %>%
    inner_join(latest_vals, by = c("geo_value", "time_value")) %>%
    # Standard ordering
    select(geo_value, time_value, value_first, value_latest) %>%
    filter(value_latest != 0, !is.na(value_latest), !is.na(value_first)) %>%
    mutate(initial_bias = (value_first - value_latest) / abs(value_latest))

bias_by_loc <- bias_df %>%
    group_by(geo_value) %>%
    summarize(
        median_bias = median(initial_bias, na.rm = TRUE),
        pct_under   = mean(initial_bias < 0, na.rm = TRUE),
        pct_over    = mean(initial_bias > 0, na.rm = TRUE),
        .groups     = "drop"
    ) %>%
    arrange(median_bias)

bias_df %>%
    left_join(bias_by_loc %>% select(geo_value, median_bias), by = "geo_value") %>%
    mutate(
        geo_value = factor(geo_value, levels = bias_by_loc$geo_value),
        direction = case_when(
            median_bias < 0 ~ "Under-reporting",
            median_bias > 0 ~ "Over-reporting",
            TRUE ~ "No bias"
        )
    ) %>%
    ggplot(aes(x = geo_value, y = initial_bias, fill = direction)) +
    geom_boxplot(outlier.shape = NA, alpha = 0.7) +
    geom_hline(yintercept = 0, linewidth = 0.4, color = "gray40") +
    scale_y_continuous(labels = label_percent(), breaks = breaks_pretty()) +
    coord_cartesian(ylim = c(
        min(0, quantile(bias_df$initial_bias, 0.025, na.rm = TRUE)),
        max(0, quantile(bias_df$initial_bias, 0.975, na.rm = TRUE))
    )) +
    labs(
        title = "Initial-release bias distribution by location",
        subtitle = "Distribution of (first_value − latest_value) / latest_value",
        x = NULL, y = "Relative bias", fill = NULL
    ) +
    theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 7))

The table below provides a detailed breakdown of these maturation patterns for each location. It includes the median relative bias along with the percentage of observations that were under-reported (initial release was lower than the final value) or over-reported (initial release was higher).

7.2 Total revision

While the previous section focused on the first release, this analysis pools every version published for an observation. This provides a comprehensive view of the total volatility and systematic bias across the entire maturation history for each location.

total_bias_df <- archive$DT %>%
    inner_join(latest_vals, by = c("geo_value", "time_value")) %>%
    filter(value_latest != 0, !is.na(value_latest)) %>%
    mutate(rel_bias = (value - value_latest) / abs(value_latest))

total_bias_by_loc <- total_bias_df %>%
    group_by(geo_value) %>%
    summarise(
        median_bias = median(rel_bias, na.rm = TRUE),
        sd_bias     = sd(rel_bias, na.rm = TRUE),
        .groups     = "drop"
    ) %>%
    arrange(median_bias)

total_bias_df %>%
    left_join(total_bias_by_loc %>% select(geo_value, median_bias), by = "geo_value") %>%
    mutate(
        geo_value = factor(geo_value, levels = total_bias_by_loc$geo_value),
        direction = case_when(
            median_bias < 0 ~ "Under-reporting",
            median_bias > 0 ~ "Over-reporting",
            TRUE ~ "No bias"
        )
    ) %>%
    ggplot(aes(x = geo_value, y = rel_bias, fill = direction)) +
    geom_boxplot(outlier.shape = NA, alpha = 0.7) +
    geom_hline(yintercept = 0, linewidth = 0.4, color = "gray40") +
    scale_y_continuous(labels = label_percent(), breaks = breaks_pretty()) +
    coord_cartesian(ylim = c(
        min(0, quantile(total_bias_df$rel_bias, 0.025, na.rm = TRUE)),
        max(0, quantile(total_bias_df$rel_bias, 0.975, na.rm = TRUE))
    )) +
    scale_fill_manual(
        values = c(
            "Under-reporting" = dis_pal[1],
            "Over-reporting"  = dis_pal[2],
            "No bias"         = "gray70"
        )
    ) +
    labs(
        title = "Total revision bias distribution by location",
        subtitle = "Distribution of (all_versions − latest_value) / latest_value",
        x = NULL, y = "Relative bias", fill = NULL
    ) +
    theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 7))

8 Geographic variation in revision window

While the overall signal may stabilize within a consistent timeframe, reporting maturity often varies by geography. The plot below highlights how long it takes for each location to reach various stages of stability (50%, 20%, 10%, and 5% of their final value).

# Re-sort locations by their median stabilization lag at the tightest threshold (5%)
geo_sorting <- stability_combined %>%
    filter(threshold == "5%") %>%
    group_by(geo_value) %>%
    summarize(med_lag = median(lag_stable, na.rm = TRUE), .groups = "drop") %>%
    arrange(desc(med_lag)) %>%
    pull(geo_value)

# Fallback to ensure no NAs if sorting fails to match
if (length(geo_sorting) == 0) geo_sorting <- sort(unique(stability_combined$geo_value))

rev_by_loc <- rev_beh %>%
    filter(n_revisions > 0) %>%
    group_by(geo_value) %>%
    summarize(
        n_obs = n(),
        avg_revisions = mean(n_revisions, na.rm = TRUE),
        median_window_days = median(as_days(max_lag), na.rm = TRUE),
        p90_window_days = quantile(as_days(max_lag), 0.9, na.rm = TRUE),
        avg_rel_spread = mean(rel_spread, na.rm = TRUE),
        median_conv_days = median(as_days(lag_near_latest), na.rm = TRUE),
        .groups = "drop"
    )

# Plot maturation stages across all thresholds defined in Section 5.1
stability_combined %>%
    mutate(
        geo_value = factor(geo_value, levels = geo_sorting)
    ) %>%
    ggplot(aes(x = geo_value, y = lag_stable, color = threshold)) +
    stat_summary(
        fun.data = median_hilow, fun.args = list(conf.int = 0.5),
        geom = "pointrange", position = position_dodge(width = 0.7)
    ) +
    scale_y_continuous(breaks = breaks_pretty()) +
    labs(
        title = "Maturation stages by location",
        subtitle = "Days from first publication to reach 50%, 20%, 10%, and 5% stability (Dots = Median, Lines = IQR)",
        x = "Location", y = "Days since first publication",
        color = "Threshold"
    ) +
    theme(
        axis.text.x = element_text(angle = 45, hjust = 1, size = 8)
    )

---

9 Convergence trajectory

Here, it is measured how closely revisions approach the final reported value over time, using the absolute relative difference to the latest known value.

9.1 Overall distance to final value

To understand how quickly the data stabilizes, we measure its distance to the final reported value over time. For each observation, we identify its initial release and the version we currently consider the ground truth. We then calculate the absolute relative difference (|value - latest_value| / |latest_value|) across every version of that data point. By averaging these errors based on the number of days since the first release, we can see how the gap closes over time. A value of 0% indicates full convergence, while a 10% difference suggests that, on average, the data point is still moving toward its eventual value even several days after its initial publication.

overall_conv <- convergence_df %>%
    group_by(days_since_first) %>%
    summarize(
        median_rel_diff = median(abs_rel_diff, na.rm = TRUE),
        n_obs = n(),
        .groups = "drop"
    ) %>%
    # Cap days to a reasonable plotting limit if it gets too sparse at the tail
    filter(n_obs > 5)

p1 <- ggplot(overall_conv, aes(x = days_since_first, y = median_rel_diff)) +
    geom_line(color = dis_pal[1], linewidth = 1) +
    scale_y_continuous(labels = label_percent(), breaks = breaks_pretty()) +
    labs(
        title = "Distance to final value across all locations",
        subtitle = "Median absolute relative error vs. days since first publication",
        x = "Days since first publication",
        y = "Median absolute rel. difference"
    )

p2 <- ggplot(overall_conv, aes(x = days_since_first, y = n_obs)) +
    geom_col(fill = "gray70") +
    scale_y_continuous(labels = label_comma(), breaks = breaks_pretty()) +
    labs(
        x = "Days since first release",
        y = "Number of observations"
    )

p1 / p2 + plot_layout(heights = c(3, 1))

9.2 Convergence by worst locations

# Identify worst locations by relative spread
n_worst <- params$n_worst %||% 12
worst_geos <- location_stats %>%
    arrange(desc(rel_spread)) %>%
    slice_head(n = n_worst) %>%
    pull(geo_value)
n_worst_actual <- length(worst_geos)

worst_geom_conv <- convergence_df %>%
    filter(geo_value %in% worst_geos) %>%
    group_by(geo_value, days_since_first) %>%
    summarize(
        median_rel_diff = median(abs_rel_diff, na.rm = TRUE),
        n_obs = n(),
        .groups = "drop"
    ) %>%
    filter(n_obs > 1) %>% # Filter extremely sparse points
    mutate(geo_value = factor(geo_value, levels = worst_geos))

ggplot(worst_geom_conv, aes(x = days_since_first)) +
    geom_col(
        aes(y = n_obs / max(n_obs) * max(median_rel_diff, na.rm = TRUE)),
        fill = "gray80", alpha = 0.5
    ) +
    geom_line(aes(y = median_rel_diff), color = dis_pal[1], linewidth = 0.8) +
    facet_wrap(~geo_value, ncol = 3, scales = "free_y") +
    scale_y_continuous(
        labels = label_percent(),
        breaks = breaks_pretty(),
        name = "Median absolute rel. diff."
    ) +
    labs(
        title = sprintf("Distance to final value in the worst %d locations", n_worst_actual),
        subtitle = "Line: median relative error; Bars: observation count vs. days since first publication",
        x = "Days since first publication"
    ) +
    theme(strip.text.x = element_text(size = 8))

---

10 Summary

For Hospital Admissions: Smoothed COVID-19 (hospital-admissions / smoothed_covid19_from_claims), the version history covering 2020-01-01 to 2023-12-31 shows the following revision and backfill profile:

• 0% of observations were never revised after their initial report.
• The median time to persistent stability is 1031 days: half of all revised observations reached their final value within 1031 days of their first publication.
• The 90th percentile for maturity is 1496 days (~214 weeks).
• 45% of locations exhibit systematic under-reporting in their initial release.
• Median days to reach stability: 50%: 0 days | 20%: 12 days | 10%: 37 days | 5%: 68 days | 2%: 70 days | Final revision: 1031 days.