diff --git a/NAMESPACE b/NAMESPACE
index 1d639a5e2..9dbe658ee 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -242,6 +242,7 @@ importFrom(scoringRules,es_sample)
 importFrom(scoringRules,logs_sample)
 importFrom(scoringRules,rps_probs)
 importFrom(scoringRules,vs_sample)
+importFrom(stats,as.formula)
 importFrom(stats,cor)
 importFrom(stats,mad)
 importFrom(stats,median)
diff --git a/NEWS.md b/NEWS.md
index d720d70c7..016f79d1b 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,7 @@
 # scoringutils (development version)
 
+- `get_pairwise_comparisons()`, and therefore `add_relative_skill()`, is now substantially faster and uses much less memory. Scores are pivoted once into a forecast unit by comparator matrix instead of being merged separately for every pair of comparators. Results are unchanged. Scores with more than one row per forecast unit and comparator now produce an informative error instead of silently comparing duplicated rows (#1221, thanks to @annakrystalli for the analysis and prototype).
+- `add_relative_skill()` no longer runs a statistical test for each pair of comparators by default (`test_type = NULL`), as it does not return the resulting p-values. This removes unnecessary computation and spurious warnings from `wilcox.test()` when scores are tied. Relative skill scores are unchanged. A test can still be requested via the new `test_type` argument (#1222, thanks to @annakrystalli).
 - Fixed several validation and messaging issues (#1211): `as_forecast_sample()`, `as_forecast_quantile()` and `as_forecast_multivariate_sample()` now error at validation time when `observed` or `predicted` are not numeric, instead of failing later inside `score()`; the rounding warning in `as_forecast_quantile()` now correctly states that quantile levels are rounded to 9 digits; and `get_pit_histogram()` for quantile-based forecasts now displays its full warning message and actually falls back to the quantiles present in the forecast when requested quantiles are missing, instead of returning an empty or incorrect result.
 - Added `filter_scores()` and `impute_missing_scores()` for handling missing forecasts before summarisation. `filter_scores()` removes target combinations with insufficient model coverage, while `impute_missing_scores()` fills in missing scores using configurable strategies (worst, mean, NA, or reference model). Both use a strategy function pattern for extensibility. See `vignette("handling-missing-forecasts")` for details (#1122).
 - Added `plot_discrimination()` to visualise the discrimination ability of binary forecasts by plotting the distribution of predicted probabilities, stratified by the observed outcome. The function requires a `forecast_binary` object (created with `as_forecast_binary()`) (#942).
diff --git a/R/pairwise-comparisons.R b/R/pairwise-comparisons.R
index 741507396..059a189ac 100644
--- a/R/pairwise-comparisons.R
+++ b/R/pairwise-comparisons.R
@@ -246,9 +246,15 @@ get_pairwise_comparisons <- function(
 #' [get_pairwise_comparisons()] splits the data into arbitrary subgroups
 #' specified by the user (e.g. if pairwise comparison should be done separately
 #' for different forecast targets) and then the actual pairwise comparison for
-#' that subgroup is managed from [pairwise_comparison_one_group()]. In order to
-#' actually do the comparison between two models over a subset of common
-#' forecasts it calls [compare_forecasts()].
+#' that subgroup is managed from [pairwise_comparison_one_group()].
+#'
+#' Internally, the scores are pivoted once into a matrix with one row per
+#' forecast unit (excluding the `compare` column) and one column per
+#' comparator (see [.pivot_scores()]). The set of overlapping forecasts for
+#' any pair of comparators is then simply the set of rows for which both
+#' columns are non-missing. This avoids re-joining the scores for every pair
+#' of comparators and is considerably faster than calling
+#' [compare_forecasts()] for every pair, while giving identical results.
 #' @inherit get_pairwise_comparisons params return
 #' @importFrom cli cli_abort
 #' @importFrom data.table setnames
@@ -285,16 +291,29 @@ pairwise_comparison_one_group <- function(scores,
   combinations <- as.data.table(t(combn(comparators, m = 2)))
   colnames(combinations) <- c("..compare", "compare_against")
 
-  combinations[, c("ratio", "pval") := compare_forecasts(
-    compare = compare,
-    scores = scores,
-    name_comparator1 = ..compare,
-    name_comparator2 = compare_against,
-    metric = metric,
-    ...
-  ),
-  by = seq_len(NROW(combinations))
-  ]
+  # pivot the scores once into a forecast unit x comparator matrix. For every
+  # pair of comparators the overlapping forecasts are the rows where both
+  # columns are non-missing.
+  score_matrix <- .pivot_scores(scores, compare = compare, metric = metric)
+  idx1 <- match(as.character(combinations$..compare), colnames(score_matrix))
+  idx2 <- match(
+    as.character(combinations$compare_against), colnames(score_matrix)
+  )
+
+  ratios <- rep(NA_real_, nrow(combinations))
+  pvals <- rep(NA_real_, nrow(combinations))
+  for (i in seq_len(nrow(combinations))) {
+    values_x <- score_matrix[, idx1[i]]
+    values_y <- score_matrix[, idx2[i]]
+    overlap <- !is.na(values_x) & !is.na(values_y)
+    if (!any(overlap)) {
+      next
+    }
+    comparison <- .compare_scores(values_x[overlap], values_y[overlap], ...)
+    ratios[i] <- comparison$mean_scores_ratio
+    pvals[i] <- comparison$pval
+  }
+  combinations[, `:=`(ratio = ratios, pval = pvals)]
 
   combinations <- combinations[order(ratio)]
   combinations[, adj_pval := p.adjust(pval)]
@@ -375,18 +394,134 @@ pairwise_comparison_one_group <- function(scores,
   return(out[])
 }
 
+#' @title Pivot scores into a forecast unit by comparator matrix
+#'
+#' @description
+#' Pivots a set of scores into a matrix with one row per forecast unit
+#' (excluding the `compare` column) and one column per comparator. Entries
+#' are the values of `metric` and are `NA` where a comparator did not
+#' provide a forecast for a given forecast unit.
+#'
+#' The function is used by [pairwise_comparison_one_group()] to align the
+#' scores of all comparators once, rather than once per pair of comparators.
+#' It errors if there is more than one score for the same forecast unit and
+#' comparator, as the scores could then not be pivoted unambiguously.
+#' @inheritParams get_pairwise_comparisons
+#' @returns A numeric matrix with one row per forecast unit and one column
+#'   per comparator. Column names are the comparators (as character).
+#' @importFrom data.table as.data.table dcast
+#' @importFrom stats as.formula
+#' @importFrom cli cli_abort
+#' @keywords internal
+.pivot_scores <- function(scores, compare = "model", metric) {
+  forecast_unit <- get_forecast_unit(scores)
+  merge_by <- setdiff(forecast_unit, compare)
+
+  # remove exact duplicates once here, instead of inside every pairwise merge
+  scores <- unique(as.data.table(scores))
+  if (anyDuplicated(scores, by = forecast_unit) > 0) {
+    #nolint start: object_usage_linter
+    cli_abort(
+      c(
+        `!` = "Found more than one score for the same forecast unit and
+        element of {.var {compare}}.",
+        i = "Pairwise comparisons require exactly one score per forecast unit
+        and comparator. Consider summarising the scores first using
+        {.fn summarise_scores}."
+      )
+    )
+    #nolint end
+  }
+
+  # column names are wrapped in backticks so that non-syntactic names work
+  if (length(merge_by) == 0) {
+    lhs <- "."
+  } else {
+    lhs <- paste0("`", merge_by, "`", collapse = " + ")
+  }
+  pivot_formula <- as.formula(paste0(lhs, " ~ `", compare, "`"))
+  wide <- dcast(
+    scores[, c(merge_by, compare, metric), with = FALSE],
+    pivot_formula,
+    value.var = metric
+  )
+
+  # the first columns of the wide table are the columns in `merge_by` (or a
+  # single placeholder column if `merge_by` is empty)
+  value_cols <- setdiff(names(wide), c(merge_by, "."))
+  score_matrix <- as.matrix(wide[, value_cols, with = FALSE])
+  colnames(score_matrix) <- value_cols
+  return(score_matrix)
+}
+
+#' @title Compare two aligned vectors of scores
+#'
+#' @description
+#' Computes the mean score ratio and, optionally, a p-value for two vectors
+#' of scores that have already been aligned, i.e. where `values_x[i]` and
+#' `values_y[i]` are the scores of two comparators for the same forecast
+#' unit. This is the shared computational core of
+#' [pairwise_comparison_one_group()] and [compare_forecasts()].
+#' @param values_x Numeric vector of scores of the first comparator.
+#' @param values_y Numeric vector of scores of the second comparator, aligned
+#'   with `values_x`.
+#' @inheritParams compare_forecasts
+#' @inherit compare_forecasts return
+#' @importFrom stats wilcox.test
+#' @keywords internal
+.compare_scores <- function(
+  values_x,
+  values_y,
+  one_sided = FALSE,
+  test_type = c("non_parametric", "permutation", NULL),
+  n_permutations = 999
+) {
+  # calculate ratio to of average scores achieved by both comparator.
+  # this should be equivalent to theta_ij in Johannes Bracher's document.
+  # ratio < 1 --> comparator 1 is better.
+  # note we could also take mean(values_x) / mean(values_y), as it cancels out
+  ratio <- sum(values_x) / sum(values_y)
+
+  # If test_type is NULL, return NA for p-value
+  if (is.null(test_type)) {
+    pval <- NA_real_
+  } else {
+    # test whether the ratio is significantly different from one
+    # equivalently, one can test whether the difference between the two values
+    # is significantly different from zero.
+    test_type <- match.arg(test_type)
+    if (test_type == "permutation") {
+      # adapted from the surveillance package
+      pval <- permutation_test(values_x, values_y,
+        n_permutation = n_permutations,
+        one_sided = one_sided,
+        comparison_mode = "difference"
+      )
+    } else {
+      # this probably needs some more thought
+      # alternative: do a paired t-test on ranks?
+      pval <- wilcox.test(values_x, values_y, paired = TRUE)$p.value
+    }
+  }
+
+  return(list(
+    mean_scores_ratio = ratio,
+    pval = pval
+  ))
+}
+
 #' @title Compare a subset of common forecasts
 #'
 #' @description
-#' This function compares two comparators based on the subset of forecasts for which
-#' both comparators have made a prediction. It gets called
-#' from [pairwise_comparison_one_group()], which handles the
-#' comparison of multiple comparators on a single set of forecasts (there are no
-#' subsets of forecasts to be distinguished). [pairwise_comparison_one_group()]
-#' in turn gets called from from [get_pairwise_comparisons()] which can handle
-#' pairwise comparisons for a set of forecasts with multiple subsets, e.g.
-#' pairwise comparisons for one set of forecasts, but done separately for two
-#' different forecast targets.
+#' This function compares two comparators based on the subset of forecasts for
+#' which both comparators have made a prediction. The overlapping forecasts
+#' are found by merging the scores of the two comparators on the forecast
+#' unit. The actual comparison is then done by [.compare_scores()].
+#'
+#' [pairwise_comparison_one_group()] no longer calls this function for every
+#' pair of comparators (it aligns the scores of all comparators at once using
+#' [.pivot_scores()] instead), but it is kept as a simple reference
+#' implementation of the comparison between two comparators.
 #' @inheritParams get_pairwise_comparisons
 #' @param name_comparator1 Character, name of the first comparator
 #' @param name_comparator2 Character, name of the comparator to compare against
@@ -440,37 +575,12 @@ compare_forecasts <- function(scores,
   values_x <- overlap[[paste0(metric, ".x")]]
   values_y <- overlap[[paste0(metric, ".y")]]
 
-  # calculate ratio to of average scores achieved by both comparator.
-  # this should be equivalent to theta_ij in Johannes Bracher's document.
-  # ratio < 1 --> comparator 1 is better.
-  # note we could also take mean(values_x) / mean(values_y), as it cancels out
-  ratio <- sum(values_x) / sum(values_y)
-
-  # If test_type is NULL, return NA for p-value
-  if (is.null(test_type)) {
-    pval <- NA_real_
-  } else {
-    # test whether the ratio is significantly different from one
-    # equivalently, one can test whether the difference between the two values
-    # is significantly different from zero.
-    test_type <- match.arg(test_type)
-    if (test_type == "permutation") {
-      # adapted from the surveillance package
-      pval <- permutation_test(values_x, values_y,
-        n_permutation = n_permutations,
-        one_sided = one_sided,
-        comparison_mode = "difference"
-      )
-    } else {
-      # this probably needs some more thought
-      # alternative: do a paired t-test on ranks?
-      pval <- wilcox.test(values_x, values_y, paired = TRUE)$p.value
-    }
-  }
-
-  return(list(
-    mean_scores_ratio = ratio,
-    pval = pval
+  return(.compare_scores(
+    values_x = values_x,
+    values_y = values_y,
+    one_sided = one_sided,
+    test_type = test_type,
+    n_permutations = n_permutations
   ))
 }
 
@@ -563,7 +673,17 @@ permutation_test <- function(scores1,
 #' Relative skill will be calculated for the aggregation level specified in
 #' `by`.
 #'
+#' Unlike [get_pairwise_comparisons()], this function does not return
+#' p-values. By default no statistical test is therefore run for the
+#' pairwise comparisons (`test_type = NULL`), which avoids unnecessary
+#' computation. Relative skill scores do not depend on `test_type`.
+#'
 #' @inheritParams get_pairwise_comparisons
+#' @param test_type Character, either "non_parametric", "permutation", or
+#'   `NULL` (the default). Determines which kind of test is run for the
+#'   pairwise comparisons. As p-values are not returned by
+#'   `add_relative_skill()`, no test is run by default. See
+#'   [compare_forecasts()] for more information.
 #' @export
 #' @keywords scoring
 add_relative_skill <- function(
@@ -572,6 +692,7 @@ add_relative_skill <- function(
   by = NULL,
   metric = intersect(c("wis", "crps", "brier_score"), names(scores)),
   baseline = NULL,
+  test_type = NULL,
   ...
 ) {
 
@@ -583,6 +704,7 @@ add_relative_skill <- function(
     baseline = baseline,
     compare = compare,
     by = by,
+    test_type = test_type,
     ...
   )
 
diff --git a/man/add_relative_skill.Rd b/man/add_relative_skill.Rd
index c57ae7852..031be8b4b 100644
--- a/man/add_relative_skill.Rd
+++ b/man/add_relative_skill.Rd
@@ -10,6 +10,7 @@ add_relative_skill(
   by = NULL,
   metric = intersect(c("wis", "crps", "brier_score"), names(scores)),
   baseline = NULL,
+  test_type = NULL,
   ...
 )
 }
@@ -37,6 +38,12 @@ given, then a scaled relative skill with respect to the baseline will be
 returned. By default (\code{NULL}), relative skill will not be scaled with
 respect to a baseline model.}
 
+\item{test_type}{Character, either "non_parametric", "permutation", or
+\code{NULL} (the default). Determines which kind of test is run for the
+pairwise comparisons. As p-values are not returned by
+\code{add_relative_skill()}, no test is run by default. See
+\code{\link[=compare_forecasts]{compare_forecasts()}} for more information.}
+
 \item{...}{Additional arguments for the comparison between two models. See
 \code{\link[=compare_forecasts]{compare_forecasts()}} for more information.}
 }
@@ -48,4 +55,10 @@ the computation of relative skill, see \code{\link[=get_pairwise_comparisons]{ge
 Relative skill will be calculated for the aggregation level specified in
 \code{by}.
 }
+\details{
+Unlike \code{\link[=get_pairwise_comparisons]{get_pairwise_comparisons()}}, this function does not return
+p-values. By default no statistical test is therefore run for the
+pairwise comparisons (\code{test_type = NULL}), which avoids unnecessary
+computation. Relative skill scores do not depend on \code{test_type}.
+}
 \keyword{scoring}
diff --git a/man/compare_forecasts.Rd b/man/compare_forecasts.Rd
index f7c658beb..ad7ce9b88 100644
--- a/man/compare_forecasts.Rd
+++ b/man/compare_forecasts.Rd
@@ -48,15 +48,15 @@ A list with mean score ratios and p-values for the comparison
 between two comparators
 }
 \description{
-This function compares two comparators based on the subset of forecasts for which
-both comparators have made a prediction. It gets called
-from \code{\link[=pairwise_comparison_one_group]{pairwise_comparison_one_group()}}, which handles the
-comparison of multiple comparators on a single set of forecasts (there are no
-subsets of forecasts to be distinguished). \code{\link[=pairwise_comparison_one_group]{pairwise_comparison_one_group()}}
-in turn gets called from from \code{\link[=get_pairwise_comparisons]{get_pairwise_comparisons()}} which can handle
-pairwise comparisons for a set of forecasts with multiple subsets, e.g.
-pairwise comparisons for one set of forecasts, but done separately for two
-different forecast targets.
+This function compares two comparators based on the subset of forecasts for
+which both comparators have made a prediction. The overlapping forecasts
+are found by merging the scores of the two comparators on the forecast
+unit. The actual comparison is then done by \code{\link[=.compare_scores]{.compare_scores()}}.
+
+\code{\link[=pairwise_comparison_one_group]{pairwise_comparison_one_group()}} no longer calls this function for every
+pair of comparators (it aligns the scores of all comparators at once using
+\code{\link[=.pivot_scores]{.pivot_scores()}} instead), but it is kept as a simple reference
+implementation of the comparison between two comparators.
 }
 \author{
 Johannes Bracher, \email{johannes.bracher@kit.edu}
diff --git a/man/dot-compare_scores.Rd b/man/dot-compare_scores.Rd
new file mode 100644
index 000000000..e0106a778
--- /dev/null
+++ b/man/dot-compare_scores.Rd
@@ -0,0 +1,43 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/pairwise-comparisons.R
+\name{.compare_scores}
+\alias{.compare_scores}
+\title{Compare two aligned vectors of scores}
+\usage{
+.compare_scores(
+  values_x,
+  values_y,
+  one_sided = FALSE,
+  test_type = c("non_parametric", "permutation", NULL),
+  n_permutations = 999
+)
+}
+\arguments{
+\item{values_x}{Numeric vector of scores of the first comparator.}
+
+\item{values_y}{Numeric vector of scores of the second comparator, aligned
+with \code{values_x}.}
+
+\item{one_sided}{Boolean, default is \code{FALSE}, whether two conduct a one-sided
+instead of a two-sided test to determine significance in a pairwise
+comparison.}
+
+\item{test_type}{Character, either "non_parametric" (the default), "permutation",
+or NULL. This determines which kind of test shall be conducted to determine
+p-values. If NULL, no test will be conducted and p-values will be NA.}
+
+\item{n_permutations}{Numeric, the number of permutations for a
+permutation test. Default is 999.}
+}
+\value{
+A list with mean score ratios and p-values for the comparison
+between two comparators
+}
+\description{
+Computes the mean score ratio and, optionally, a p-value for two vectors
+of scores that have already been aligned, i.e. where \code{values_x[i]} and
+\code{values_y[i]} are the scores of two comparators for the same forecast
+unit. This is the shared computational core of
+\code{\link[=pairwise_comparison_one_group]{pairwise_comparison_one_group()}} and \code{\link[=compare_forecasts]{compare_forecasts()}}.
+}
+\keyword{internal}
diff --git a/man/dot-pivot_scores.Rd b/man/dot-pivot_scores.Rd
new file mode 100644
index 000000000..fe144d5fc
--- /dev/null
+++ b/man/dot-pivot_scores.Rd
@@ -0,0 +1,37 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/pairwise-comparisons.R
+\name{.pivot_scores}
+\alias{.pivot_scores}
+\title{Pivot scores into a forecast unit by comparator matrix}
+\usage{
+.pivot_scores(scores, compare = "model", metric)
+}
+\arguments{
+\item{scores}{An object of class \code{scores} (a data.table with
+scores and an additional attribute \code{metrics} as produced by \code{\link[=score]{score()}}).}
+
+\item{compare}{Character vector with a single colum name that defines the
+elements for the pairwise comparison. For example, if this is set to
+"model" (the default), then elements of the "model" column will be
+compared.}
+
+\item{metric}{A string with the name of the metric for which
+a relative skill shall be computed. By default this is either "crps",
+"wis" or "brier_score" if any of these are available.}
+}
+\value{
+A numeric matrix with one row per forecast unit and one column
+per comparator. Column names are the comparators (as character).
+}
+\description{
+Pivots a set of scores into a matrix with one row per forecast unit
+(excluding the \code{compare} column) and one column per comparator. Entries
+are the values of \code{metric} and are \code{NA} where a comparator did not
+provide a forecast for a given forecast unit.
+
+The function is used by \code{\link[=pairwise_comparison_one_group]{pairwise_comparison_one_group()}} to align the
+scores of all comparators once, rather than once per pair of comparators.
+It errors if there is more than one score for the same forecast unit and
+comparator, as the scores could then not be pivoted unambiguously.
+}
+\keyword{internal}
diff --git a/man/pairwise_comparison_one_group.Rd b/man/pairwise_comparison_one_group.Rd
index 2982a753f..4c4bd3bfb 100644
--- a/man/pairwise_comparison_one_group.Rd
+++ b/man/pairwise_comparison_one_group.Rd
@@ -54,8 +54,14 @@ multiple models involved. It gets called from \code{\link[=get_pairwise_comparis
 \code{\link[=get_pairwise_comparisons]{get_pairwise_comparisons()}} splits the data into arbitrary subgroups
 specified by the user (e.g. if pairwise comparison should be done separately
 for different forecast targets) and then the actual pairwise comparison for
-that subgroup is managed from \code{\link[=pairwise_comparison_one_group]{pairwise_comparison_one_group()}}. In order to
-actually do the comparison between two models over a subset of common
-forecasts it calls \code{\link[=compare_forecasts]{compare_forecasts()}}.
+that subgroup is managed from \code{\link[=pairwise_comparison_one_group]{pairwise_comparison_one_group()}}.
+
+Internally, the scores are pivoted once into a matrix with one row per
+forecast unit (excluding the \code{compare} column) and one column per
+comparator (see \code{\link[=.pivot_scores]{.pivot_scores()}}). The set of overlapping forecasts for
+any pair of comparators is then simply the set of rows for which both
+columns are non-missing. This avoids re-joining the scores for every pair
+of comparators and is considerably faster than calling
+\code{\link[=compare_forecasts]{compare_forecasts()}} for every pair, while giving identical results.
 }
 \keyword{internal}
diff --git a/tests/testthat/test-pairwise_comparison.R b/tests/testthat/test-pairwise_comparison.R
index 543a496e9..ed3ccd0a2 100644
--- a/tests/testthat/test-pairwise_comparison.R
+++ b/tests/testthat/test-pairwise_comparison.R
@@ -605,3 +605,156 @@ test_that("add_relative_skill() works without warnings when not computing p-valu
   expect_type(scores_w_rel_skill$ae_median_relative_skill, "double")
   expect_false(anyNA(scores_w_rel_skill$ae_median_relative_skill))
 })
+
+# tests for the pivot-based implementation ------------------------------------
+
+test_that(".pivot_scores() pivots scores into a forecast unit x model matrix", {
+  scores <- data.table::copy(scores_quantile)
+  # drop some forecasts for one model to create missing overlap
+  scores <- scores[!(model == "EuroCOVIDhub-ensemble" & location == "DE")]
+
+  score_matrix <- .pivot_scores(scores, compare = "model", metric = "wis")
+  forecast_unit <- setdiff(get_forecast_unit(scores), "model")
+  n_units <- nrow(unique(scores[, forecast_unit, with = FALSE]))
+
+  expect_true(is.matrix(score_matrix))
+  expect_type(score_matrix, "double")
+  expect_identical(nrow(score_matrix), n_units)
+  expect_setequal(colnames(score_matrix), unique(scores$model))
+  expect_identical(
+    sum(!is.na(score_matrix[, "EuroCOVIDhub-ensemble"])),
+    nrow(scores[model == "EuroCOVIDhub-ensemble"])
+  )
+  expect_identical(
+    sum(!is.na(score_matrix)),
+    nrow(scores)
+  )
+})
+
+test_that(".pivot_scores() errors with duplicated scores per forecast unit", {
+  scores <- data.table::copy(scores_quantile)
+  metrics <- get_metrics(scores)
+  duplicated_scores <- rbind(scores, scores[1:5][, wis := wis + 1])
+  duplicated_scores <- new_scores(duplicated_scores, metrics = metrics)
+  expect_error(
+    .pivot_scores(duplicated_scores, compare = "model", metric = "wis"),
+    "more than one score for the same forecast unit"
+  )
+  # exact duplicates are removed rather than causing an error
+  exact_duplicates <- new_scores(rbind(scores, scores[1:5]), metrics = metrics)
+  expect_no_error(
+    .pivot_scores(exact_duplicates, compare = "model", metric = "wis")
+  )
+})
+
+test_that("get_pairwise_comparisons() matches per-pair compare_forecasts()", {
+  scores <- data.table::copy(scores_quantile)
+  # drop some forecasts so that models do not overlap perfectly and so that
+  # one pair of models has no overlapping forecasts at all
+  scores <- scores[!(model == "EuroCOVIDhub-ensemble" & location == "DE")]
+  scores <- scores[!(model == "EuroCOVIDhub-ensemble" & target_type == "Cases")]
+  scores <- scores[!(model == "epiforecasts-EpiNow2" & target_type == "Deaths")]
+
+  reference <- function(scores, ...) {
+    models <- unique(scores$model)
+    pairs <- data.table::as.data.table(t(combn(models, m = 2)))
+    data.table::setnames(pairs, c("model", "compare_against"))
+    pairs[, c("mean_scores_ratio", "pval") := compare_forecasts(
+      scores = scores,
+      compare = "model",
+      name_comparator1 = model,
+      name_comparator2 = compare_against,
+      metric = "wis",
+      ...
+    ), by = seq_len(nrow(pairs))]
+    pairs[]
+  }
+
+  check_against_reference <- function(pairwise, reference) {
+    merged <- merge(
+      pairwise, reference,
+      by = c("model", "compare_against"), suffixes = c("", "_ref")
+    )
+    expect_identical(nrow(merged), nrow(reference))
+    expect_identical(merged$mean_scores_ratio, merged$mean_scores_ratio_ref)
+    expect_identical(merged$pval, merged$pval_ref)
+  }
+
+  # default non-parametric test, no grouping
+  pairwise <- get_pairwise_comparisons(scores, metric = "wis")
+  check_against_reference(pairwise, reference(scores))
+
+  # pairs without any overlap get NA
+  expect_true(anyNA(pairwise$mean_scores_ratio))
+
+  # grouping via `by`
+  pairwise_by <- get_pairwise_comparisons(
+    scores, metric = "wis", by = "target_type"
+  )
+  reference_by <- scores[, reference(.SD), by = "target_type"]
+  merged <- merge(
+    pairwise_by, reference_by,
+    by = c("target_type", "model", "compare_against"),
+    suffixes = c("", "_ref")
+  )
+  expect_identical(nrow(merged), nrow(reference_by))
+  expect_identical(merged$mean_scores_ratio, merged$mean_scores_ratio_ref)
+  expect_identical(merged$pval, merged$pval_ref)
+
+  # no test
+  pairwise_notest <- get_pairwise_comparisons(
+    scores, metric = "wis", test_type = NULL
+  )
+  check_against_reference(pairwise_notest, reference(scores, test_type = NULL))
+  expect_true(all(is.na(pairwise_notest[model != compare_against]$pval)))
+
+  # permutation test with the same seed gives the same p-values
+  set.seed(42)
+  pairwise_perm <- get_pairwise_comparisons(
+    scores, metric = "wis", test_type = "permutation", n_permutations = 50
+  )
+  set.seed(42)
+  reference_perm <- reference(
+    scores, test_type = "permutation", n_permutations = 50
+  )
+  check_against_reference(pairwise_perm, reference_perm)
+})
+
+test_that("get_pairwise_comparisons() works when `compare` is a factor", {
+  scores <- data.table::copy(scores_quantile)
+  scores[, model := factor(model)]
+  pairwise_factor <- get_pairwise_comparisons(scores, metric = "wis")
+  pairwise_character <- get_pairwise_comparisons(
+    data.table::copy(scores_quantile), metric = "wis"
+  )
+  expect_identical(
+    pairwise_factor[order(model, compare_against)],
+    pairwise_character[order(model, compare_against)]
+  )
+})
+
+test_that("add_relative_skill() skips the test by default", {
+  scores <- data.table::copy(scores_quantile)
+
+  # count calls to wilcox.test() rather than relying on its warnings, which
+  # differ across R versions
+  calls <- new.env()
+  calls$n <- 0L
+  testthat::local_mocked_bindings(
+    wilcox.test = function(...) {
+      calls$n <- calls$n + 1L
+      list(p.value = NA_real_)
+    },
+    .package = "scoringutils"
+  )
+
+  without_test <- expect_no_warning(add_relative_skill(scores, metric = "wis"))
+  expect_identical(calls$n, 0L)
+  expect_false("pval" %in% colnames(without_test))
+
+  with_test <- add_relative_skill(
+    scores, metric = "wis", test_type = "non_parametric"
+  )
+  expect_gt(calls$n, 0L)
+  expect_identical(without_test, with_test)
+})