Summary
pairwise_comparison_one_group() realigns the same scores from scratch for every model pair. Pivoting the scores into a model-by-unit matrix once, then reading each pair's overlapping forecasts off two columns, gives identical output and a large speedup that grows with the number of models. Sharing a validated prototype and benchmark evidence here rather than a PR, so you can decide whether and how you'd like to adopt it.
Where it hurts
On the FluSight forecast hub, generating the evaluation dashboard data is dominated by relative-skill pairwise comparisons. In a hubverse benchmark of that build (full-season, roughly 50 models):
| stage |
time |
% |
relative skill (get_pairwise_comparisons) |
427.5 s |
62.2% |
├ of which wilcox.test |
16.4 s |
2.4% |
score() |
167.5 s |
24.4% |
| load |
32.5 s |
4.7% |
The cost is the merging, not the statistics: forderv (data.table ordering inside merge.data.table inside compare_forecasts()) is 43.8% of total self-time on its own. For context, the downstream CI build this feeds currently runs around 3.5 hours, and this is the single biggest lever.
Root cause
pairwise_comparison_one_group() forms all choose(n, 2) model pairs and calls compare_forecasts() once per pair. Each call does a fresh join to find the overlapping forecasts:
overlap <- merge(a, b, by = merge_by, allow.cartesian = TRUE)
overlap <- unique(overlap)
For 95 models that is 4465 joins per metric per group, each re-keying and re-sorting the same underlying data. The alignment is redone per pair rather than once.
Proposed approach
Pivot once into a wide matrix (rows = forecast unit minus compare, columns = comparators, values = the metric). Then each pair is two columns; the overlap is the rows where both are non-NA:
forecast_unit <- get_forecast_unit(scores)
merge_by <- setdiff(forecast_unit, compare)
# one dedup + one uniqueness check, instead of unique() inside every merge
scores <- unique(as.data.table(scores))
if (anyDuplicated(scores, by = forecast_unit)) {
cli_abort("Non-unique forecast unit per comparator; cannot pivot safely.")
}
wide <- dcast(scores, as.formula(paste(paste(merge_by, collapse = " + "), "~", compare)),
value.var = metric)
mat <- as.matrix(wide[, comparators, with = FALSE])
colnames(mat) <- comparators
# per pair: overlap = rows where both models forecast
ok <- !is.na(mat[, i]) & !is.na(mat[, j])
ratio <- sum(mat[ok, i]) / sum(mat[ok, j])Everything downstream (ordering, p.adjust, mirroring, geometric-mean theta, baseline scaling, final merge) is untouched, so the only thing that changes is how ratio and p-value are obtained.
Full drop-in pairwise_comparison_one_group() (only the marked block differs from current)
pairwise_comparison_one_group <- function(scores,
metric,
baseline,
compare = "model",
by,
test_type = c("non_parametric", "permutation", NULL),
one_sided = FALSE,
n_permutations = 999) {
if (!(compare %in% names(scores))) {
cli_abort("pairwise comparisons require a column as given by `compare`")
}
comparators <- unique(scores[[compare]])
if (length(comparators) < 2) {
cli_abort(c("!" = "There are not enough comparators to do any comparison"))
}
combinations <- as.data.table(t(combn(comparators, m = 2)))
colnames(combinations) <- c("..compare", "compare_against")
# ==== CHANGED: pivot once, then read each pair off two columns =============
forecast_unit <- get_forecast_unit(scores)
merge_by <- setdiff(forecast_unit, compare)
# done once here, versus unique() inside every per-pair merge previously
scores <- unique(as.data.table(scores))
if (anyDuplicated(scores, by = forecast_unit)) {
cli_abort("Non-unique forecast unit per comparator; cannot pivot safely.")
}
wide <- dcast(
scores,
as.formula(paste(paste(merge_by, collapse = " + "), "~", compare)),
value.var = metric
)
mat <- as.matrix(wide[, comparators, with = FALSE])
colnames(mat) <- comparators
pair_pval <- function(x, y) {
if (is.null(test_type)) return(NA_real_)
tt <- match.arg(test_type)
if (tt == "permutation") {
permutation_test(x, y, n_permutation = n_permutations,
one_sided = one_sided, comparison_mode = "difference")
} else {
wilcox.test(x, y, paired = TRUE)$p.value
}
}
stat <- vapply(seq_len(nrow(combinations)), function(k) {
x <- mat[, as.character(combinations$..compare[k])]
y <- mat[, as.character(combinations$compare_against[k])]
ok <- !is.na(x) & !is.na(y)
if (!any(ok)) return(c(ratio = NA_real_, pval = NA_real_))
x <- x[ok]; y <- y[ok]
c(ratio = sum(x) / sum(y), pval = pair_pval(x, y))
}, numeric(2))
combinations[, ratio := stat["ratio", ]]
combinations[, pval := stat["pval", ]]
# ==== end changed block; everything below is unchanged ====================
combinations <- combinations[order(ratio)]
combinations[, adj_pval := p.adjust(pval)]
combinations_mirrored <- copy(combinations)
setnames(combinations_mirrored,
old = c("..compare", "compare_against"),
new = c("compare_against", "..compare"))
combinations_mirrored[, ratio := 1 / ratio]
combinations_equal <- data.table(
..compare = comparators, compare_against = comparators,
ratio = 1, pval = 1, adj_pval = 1
)
result <- rbindlist(
list(combinations, combinations_mirrored, combinations_equal),
use.names = TRUE
)
result[, `:=`(
..compare = as.character(..compare),
compare_against = as.character(compare_against)
)]
setnames(result, old = "..compare", new = compare)
result[, theta := geometric_mean(ratio), by = compare]
if (!is.null(baseline)) {
baseline_theta <- unique(result[get(compare) == baseline, ]$theta)
if (length(baseline_theta) == 0) {
cli_abort("Baseline comparator {.var {baseline}} missing.")
}
result[, rel_to_baseline := theta / baseline_theta]
}
cols_to_keep <- unique(c(by, compare))
cols_to_remove <- colnames(scores)[!(colnames(scores) %in% cols_to_keep)]
scores[, eval(cols_to_remove) := NULL]
scores <- unique(scores)
out <- merge(scores, result, by = compare, all = TRUE)
setnames(out, old = c("ratio", "theta"),
new = c("mean_scores_ratio", paste(metric, "relative_skill", sep = "_")))
if (!is.null(baseline)) {
setnames(out, old = "rel_to_baseline",
new = paste(metric, "scaled_relative_skill", sep = "_"))
}
out[]
}The test options are shown as explicit arguments for clarity; in-tree they arrive via ... exactly as today. Equivalence was validated on the default non-parametric path.
Results
Validated against the current implementation (this prototype only swaps the per-pair merge; all later steps are the existing code):
Identical output. Byte-for-byte equal ratios, p-values, adjusted p-values, relative skill and scaled relative skill across the real FluSight scores objects (aggregate plus four disaggregations across three evaluation sets, up to ~123k output rows) and synthetic cases, with and without a baseline.
Speedup, real FluSight scores (aggregate call, ~210k rows, 48 models, 1128 pairs):
|
median |
memory allocated |
| current |
13.8 s |
19.0 GB |
| pivot |
3.6 s |
1.9 GB |
Roughly 3.8x faster and 10x less allocation, and the gap widens with model count (synthetic sweep: 2.7x at 20 models, 3.8x at 40, 4.1x at 60, 4.6x at 80), since the pivot stays a single operation while the pair count grows as n squared.
One correctness note
The current per-pair unique(overlap) guards against exact duplicate rows. The prototype does the equivalent once up front (unique() plus an anyDuplicated() assert on the forecast unit); on real hubverse data there are no duplicates, so it is a no-op safety net, but the assert makes the pivot's one-score-per-unit assumption explicit rather than letting dcast() silently aggregate.
Aside: a smaller, separate saving
Unrelated to the merge, but noticed alongside it: add_relative_skill() drops the pval and adj_pval columns before returning (it keeps only relative skill), yet it calls through with the default test_type = "non_parametric", so the Wilcoxon tests are computed and then discarded. In this benchmark that is the full 16.4 s (2.4%) of wilcox.test time, all of it thrown away. Defaulting add_relative_skill(test_type = NULL), or otherwise letting relative-skill-only callers opt out, would remove it. Much smaller than the merge and a default/behaviour choice that is yours, so flagging only as a heads-up.
Summary
pairwise_comparison_one_group()realigns the same scores from scratch for every model pair. Pivoting the scores into a model-by-unit matrix once, then reading each pair's overlapping forecasts off two columns, gives identical output and a large speedup that grows with the number of models. Sharing a validated prototype and benchmark evidence here rather than a PR, so you can decide whether and how you'd like to adopt it.Where it hurts
On the FluSight forecast hub, generating the evaluation dashboard data is dominated by relative-skill pairwise comparisons. In a hubverse benchmark of that build (full-season, roughly 50 models):
get_pairwise_comparisons)wilcox.testscore()The cost is the merging, not the statistics:
forderv(data.table ordering insidemerge.data.tableinsidecompare_forecasts()) is 43.8% of total self-time on its own. For context, the downstream CI build this feeds currently runs around 3.5 hours, and this is the single biggest lever.Root cause
pairwise_comparison_one_group()forms allchoose(n, 2)model pairs and callscompare_forecasts()once per pair. Each call does a fresh join to find the overlapping forecasts:For 95 models that is 4465 joins per metric per group, each re-keying and re-sorting the same underlying data. The alignment is redone per pair rather than once.
Proposed approach
Pivot once into a wide matrix (rows = forecast unit minus
compare, columns = comparators, values = the metric). Then each pair is two columns; the overlap is the rows where both are non-NA:Everything downstream (ordering,
p.adjust, mirroring, geometric-meantheta, baseline scaling, final merge) is untouched, so the only thing that changes is how ratio and p-value are obtained.Full drop-in
pairwise_comparison_one_group()(only the marked block differs from current)The test options are shown as explicit arguments for clarity; in-tree they arrive via
...exactly as today. Equivalence was validated on the default non-parametric path.Results
Validated against the current implementation (this prototype only swaps the per-pair merge; all later steps are the existing code):
Identical output. Byte-for-byte equal ratios, p-values, adjusted p-values, relative skill and scaled relative skill across the real FluSight scores objects (aggregate plus four disaggregations across three evaluation sets, up to ~123k output rows) and synthetic cases, with and without a baseline.
Speedup, real FluSight scores (aggregate call, ~210k rows, 48 models, 1128 pairs):
Roughly 3.8x faster and 10x less allocation, and the gap widens with model count (synthetic sweep: 2.7x at 20 models, 3.8x at 40, 4.1x at 60, 4.6x at 80), since the pivot stays a single operation while the pair count grows as n squared.
One correctness note
The current per-pair
unique(overlap)guards against exact duplicate rows. The prototype does the equivalent once up front (unique()plus ananyDuplicated()assert on the forecast unit); on real hubverse data there are no duplicates, so it is a no-op safety net, but the assert makes the pivot's one-score-per-unit assumption explicit rather than lettingdcast()silently aggregate.Aside: a smaller, separate saving
Unrelated to the merge, but noticed alongside it:
add_relative_skill()drops thepvalandadj_pvalcolumns before returning (it keeps only relative skill), yet it calls through with the defaulttest_type = "non_parametric", so the Wilcoxon tests are computed and then discarded. In this benchmark that is the full 16.4 s (2.4%) ofwilcox.testtime, all of it thrown away. Defaultingadd_relative_skill(test_type = NULL), or otherwise letting relative-skill-only callers opt out, would remove it. Much smaller than the merge and a default/behaviour choice that is yours, so flagging only as a heads-up.