diff --git a/.gitignore b/.gitignore
index f0f29c32..dbd669d2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -30,3 +30,5 @@ release-prep.R
 
 # personal maintainer scratch (not shared)
 internal-notes/
+notes/loo_se.pdf
+notes/loo-compare-se-diff.md
\ No newline at end of file
diff --git a/NEWS.md b/NEWS.md
index 78be699d..9a02a982 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -7,6 +7,10 @@
 * New predictive performance API: `insample_pred_measure()`, `loo_pred_measure()`,
   `kfold_pred_measure()`, `test_pred_measure()`, and `pred_measure()` with
   built-in measures via `measure_*()` and [supported_measures_list()].
+* `loo_compare()` now supports `loo_pred_measure` objects: paired differences
+  for all measures common to the compared models, optional `rank_by` ranking,
+  utility-scale sign conversion for loss measures, and
+  `print(compare, measures = ...)` for multi-measure tables by @florence-bockting in #380.
 
 # loo 2.10.0
 
diff --git a/R/loo-glossary.R b/R/loo-glossary.R
index bb1cb809..17c756bb 100644
--- a/R/loo-glossary.R
+++ b/R/loo-glossary.R
@@ -254,4 +254,108 @@
 #'
 #' See for further information on Pareto-k values the "Pareto k estimates"
 #' section.
+#'
+#' @section Multi-measure model comparisons:
+#'
+#' When comparing [`loo_pred_measure()`][loo_pred_measure] objects with
+#' `loo_compare()`, paired differences are computed for every predictive
+#' measure common to all models. Models are ranked by the `rank_by` argument
+#' (default `"elpd"`); the top-ranked model is the reference for all difference
+#' columns.
+#'
+#' ### `{measure}_diff` and `{measure}_se_diff`
+#'
+#' For each non-ELPD measure `m`, `loo_compare()` adds columns `m_diff` and
+#' `m_se_diff`. In all cases `m_diff` is the difference between the two overall
+#' estimates on a utility scale (higher is better; loss measures such as MSE,
+#' Brier score, and SRPS have their sign flipped from the raw loss orientation).
+#' Measures already returned on a utility scale (e.g. ELPD, CRPS/RPS) are not
+#' sign-flipped. Negative `m_diff` values then indicate worse performance than
+#' the reference model, which has `m_diff = 0`.
+#'
+#' How `m_se_diff` is obtained depends on the measure:
+#'
+#' * When the overall estimate is a sum or mean of pointwise contributions, it
+#'   is computed from paired pointwise differences using the same approach as
+#'   `elpd_diff` and `se_diff` (Eq 24 in VGG2017 for sums; the mean analogue for
+#'   means). This covers ELPD, `mlpd`, `ic`, `mae`, `mse`, `acc`, `brier`, and
+#'   the ranked probability scores.
+#' * When it is a transformation of such quantities, the measure supplies its
+#'   own delta-method standard error (`se_diff_fun`). For `rmse` this is the
+#'   first-order bivariate Taylor approximation propagated from the MSE scale,
+#'   which requires the covariance between the two models' pointwise squared
+#'   errors and is therefore not a paired pointwise standard deviation. For
+#'   `r2` it is the trivariate analogue, which additionally propagates the
+#'   uncertainty in the baseline `MSE(y)` shared by both models.
+#' * Otherwise `m_se_diff` is `NA`. This applies only to custom measures whose
+#'   estimate is neither a sum nor a mean of their pointwise values and which do
+#'   not attach an `se_diff_fun`.
+#'
+#' The reference model has `m_se_diff = 0` whenever an `m_se_diff` is available.
+#' Attribute `measure_higher_is_better` on each `*_pred_measure()`
+#' result records the `higher_is_better` setting used when each measure was
+#' computed; when stored values are on a loss scale, `loo_compare()` emits a
+#' short message naming those measures (see [loo_compare()]).
+#'
+#' ELPD-family measures use the column names `elpd_diff` and `se_diff` rather
+#' than a prefixed form. Only ELPD comparisons include `p_worse` and `diag_diff`;
+#' these diagnostics do not apply to other predictive measures.
+#'
+#' ### `measure_higher_is_better`
+#'
+#' Attribute on all `*_pred_measure()` and [pred_measure()] results: a named
+#' list recording the `higher_is_better` setting used for each measure (`NULL`,
+#' `TRUE`, or `FALSE`; `elpd` is always `NULL`). Used by [loo_compare()] with
+#' `measure_compare_meta` to decide whether paired differences need a sign flip
+#' when converting to a utility scale.
+#'
+#' ### `measure_compare_meta`
+#'
+#' Attribute on all `*_pred_measure()` and [pred_measure()] results: a named
+#' list of per-measure comparison metadata used by [loo_compare()]. Each entry
+#' is a list with:
+#'
+#' * `higher_is_better` — the orientation setting used when the measure was
+#'   computed (`NULL`, `TRUE`, or `FALSE`)
+#' * `loss` — whether stored values are on a loss scale (lower is better)
+#' * `diff_method` — how the standard error of the difference is obtained:
+#'   `"sum"` or `"mean"` (paired pointwise differences), `"pairwise"` (the
+#'   measure's own `se_diff_fun`), `"estimates_only"` (unavailable, `NA`), or
+#'   `"auto"` (inferred at compare time for custom measures). No built-in
+#'   measure declares `"estimates_only"`; it is reached only when
+#'   autodetection under `"auto"` cannot establish that the estimate is a sum
+#'   or mean of its pointwise values. It is not an error state — the
+#'   difference itself is still reported, and only its standard error is
+#'   marked unavailable.
+#' * `se_diff_fun` — for `diff_method = "pairwise"`, either the name of a
+#'   built-in implementation or, for custom measures, the function itself
+#' * `extra` — optional list of auxiliary data the measure stored for its
+#'   `se_diff_fun`, present only for measures that need it (`r2` stores the
+#'   pointwise baseline `(y_i - mean(y))^2`, which `y` no longer supplies by
+#'   the time [loo_compare()] runs; `bacc` stores the class index of each
+#'   observation, which its pointwise values do not determine). It is excluded
+#'   from the metadata consistency check below, since it varies with the data
+#'   rather than with how the measure was configured.
+#'
+#' Built-in measures take `loss`, `diff_method`, and `se_diff_fun` from the
+#' package measure registry; custom measures default to `loss = FALSE` and
+#' `diff_method = "auto"` unless they attach an `se_diff_fun` (see
+#' [insample_pred_measure()]).
+#' [loo_compare()] requires all models to provide matching metadata for each
+#' shared measure; mismatched `higher_is_better` settings or missing metadata on
+#' some models produce an error.
+#'
+#' ### `rank_by`, `compare_measures`, and related attributes`
+#'
+#' The `rank_by` argument selects which measure determines model ordering and
+#' the reference model for all pairwise differences. When `rank_by` is omitted,
+#' models are ranked by `"elpd"`; attribute `rank_by` is set only when `rank_by`
+#' is passed explicitly. Attribute `compare_measures` lists all measures that
+#' were compared, and `sign_converted_measures` lists loss measures whose sign
+#' was flipped onto the utility scale. The print method shows the ranking
+#' measure by default
+#' (`"elpd"` when `rank_by` was not set); use `print(x, measures = "all")` or
+#' `print(x, measures = c("rmse", "r2"))` to display additional measure tables.
+#' Printed tables label the standard-error column `se_diff` even for non-ELPD
+#' measures; the data frame columns remain `{measure}_se_diff`.
 NULL
diff --git a/R/loo_compare.R b/R/loo_compare.R
index ce5bea4a..49dc09bc 100644
--- a/R/loo_compare.R
+++ b/R/loo_compare.R
@@ -1,17 +1,38 @@
 #' Model comparison
 #'
-#' @description Compare fitted models based on [ELPD][loo-glossary].
+#' @description Compare fitted models based on [ELPD][loo-glossary] or, for
+#'   [`loo_pred_measure()`][loo_pred_measure] results, multiple predictive
+#'   performance measures at once.
 #'
 #' @export
-#' @param x An object of class `"loo"` or a list of such objects. If a list is
-#'   used then the list names will be used as the model names in the output. See
-#'   **Examples**.
-#' @param ... Additional objects of class `"loo"`, if not passed in as a single
-#'   list.
+#' @param x An object of class `"loo"` or `"loo_pred_measure"`, or a list of
+#'   such objects. If a list is used then the list names will be used as the
+#'   model names in the output. See **Examples**.
+#' @param ... Additional objects of class `"loo"` or `"loo_pred_measure"`, if not
+#'   passed in as a single list.
+#' @param rank_by For [`loo_pred_measure()`][loo_pred_measure] comparisons only,
+#'   the bare measure name used to rank models and define the reference model
+#'   for all pairwise differences (default `"elpd"`). For example,
+#'   `rank_by = "mse"` ranks models by predictive MSE (best/lowest MSE first) and
+#'   computes all measure differences relative to that model on a utility scale
+#'   (higher is better; loss measures such as MSE have their sign flipped).
 #'
 #' @return A data frame with class `"compare.loo"` that has its own
 #'   print method. See the **Details** and **Examples** sections.
 #'
+#'   For classic `"loo"` / `"waic"` / `"kfold"` comparisons, the returned
+#'   columns are unchanged from previous versions.
+#'
+#'   For [`loo_pred_measure()`][loo_pred_measure] comparisons, the data frame
+#'   additionally contains `{measure}_diff` and `{measure}_se_diff` columns for
+#'   every predictive measure common to all models (e.g. `rmse_diff`,
+#'   `rmse_se_diff`). ELPD-family measures use `elpd_diff` and `se_diff`.
+#'   `p_worse` and `diag_diff` are computed for ELPD only. Per-model PSIS
+#'   diagnostics appear in `diag_elpd`. Attributes `compare_measures` and
+#'   `sign_converted_measures` record which measures were compared and which
+#'   loss measures had their sign flipped for comparison. Attribute `rank_by` is
+#'   set when `rank_by` was passed explicitly (default ranking is by `"elpd"`).
+#'
 #' @details
 #'   When comparing two fitted models, we can estimate the difference in their
 #'   expected predictive accuracy by the difference in
@@ -67,9 +88,41 @@
 #'   computation, then there may be significant bias in `elpd_diff` favoring
 #'   models with a large number of high Pareto k values.
 #'
+#' ## Comparing `loo_pred_measure` objects
+#'   When all inputs are [`loo_pred_measure()`][loo_pred_measure] objects,
+#'   `loo_compare()` computes paired differences for every predictive measure
+#'   present in all models. Models are ranked by `rank_by` (default `"elpd"`);
+#'   the top-ranked model is the reference for all `{measure}_diff` columns.
+#'   Measures may use different orientations in their raw form (e.g. ELPD and
+#'   CRPS/RPS are returned on a utility scale where higher is better, while MSE
+#'   and Brier score are loss measures where lower is better). For comparison,
+#'   all `{measure}_diff` values are reported on a common utility scale (higher
+#'   is better). Loss measures have their sign flipped from the raw loss
+#'   orientation so that negative `{measure}_diff` values indicate worse
+#'   performance than the reference model. Each `*_pred_measure()` result stores
+#'   attribute `measure_higher_is_better`, a named list recording the
+#'   `higher_is_better` setting used when each measure was computed. When loss measures are compared
+#'   on a utility scale, `loo_compare()` emits a short message naming the affected
+#'   measures, for example:
+#'   "For model comparison, differences for mse are reported on a utility scale
+#'   (higher is better)."
+#'
+#'   `p_worse` and `diag_diff` are computed for ELPD-family measures only. Other
+#'   measures receive `{measure}_diff` and `{measure}_se_diff` from paired
+#'   pointwise contributions when the overall estimate is a sum or mean of those
+#'   contributions (using the same standard error formula as `se_diff`). For
+#'   measures where pointwise values do not define the overall estimate (e.g.
+#'   `r2`, `mse`, `rmse`), `{measure}_diff` is the difference between overall
+#'   estimates and `{measure}_se_diff` is `NA`. When models were fit with
+#'   different `measure` sets, only measures common to all models are compared; a
+#'   warning lists omitted measures. Use `print(x, measures = "all")` to display
+#'   diff tables for every compared measure; see [loo-glossary] for column
+#'   definitions.
+#'
 #' ## Warnings for many model comparisons
 #'   If more than \eqn{11} models are compared, we internally recompute the model
-#'   differences using the median model by ELPD as the baseline model. We then
+#'   differences using the median model (by ELPD, or by `rank_by` for
+#'   `loo_pred_measure` comparisons) as the baseline model. We then
 #'   estimate whether the differences in predictive performance are potentially
 #'   due to chance as described by McLatchie and Vehtari (2023). This will flag
 #'   a warning if it is deemed that there is a risk of over-fitting due to the
@@ -99,39 +152,84 @@
 #' \dontrun{
 #' # works for waic (and kfold) too
 #' loo_compare(waic(LL), waic(LL - 10))
+#'
+#' # compare multiple predictive measures from loo_pred_measure()
+#' if (requireNamespace("brms", quietly = TRUE)) {
+#'   fit1 <- brms::brm(
+#'     Reaction ~ Days, data = lme4::sleepstudy,
+#'     refresh = 0, chains = 2, iter = 1000
+#'   )
+#'   fit2 <- brms::brm(
+#'     Reaction ~ poly(Days, 2), data = lme4::sleepstudy,
+#'     refresh = 0, chains = 2, iter = 1000
+#'   )
+#'   pm1 <- loo_pred_measure(
+#'     loo = loo(fit1, save_psis = TRUE),
+#'     y = fit1$data$Reaction,
+#'     mupred = brms::posterior_epred(fit1),
+#'     measure = c("rmse", "r2")
+#'   )
+#'   pm2 <- loo_pred_measure(
+#'     loo = loo(fit2, save_psis = TRUE),
+#'     y = fit2$data$Reaction,
+#'     mupred = brms::posterior_epred(fit2),
+#'     measure = c("rmse", "r2")
+#'   )
+#'   comp <- loo_compare(pm1, pm2)
+#'   print(comp)                      # ranked by elpd (default)
+#'   print(comp, measures = "all")    # all measure diff tables
+#'   loo_compare(pm1, pm2, rank_by = "rmse")
+#' }
 #' }
 #'
-loo_compare <- function(x, ...) {
+loo_compare <- function(x, ..., rank_by = NULL) {
   UseMethod("loo_compare")
 }
 
 #' @rdname loo_compare
 #' @export
-loo_compare.default <- function(x, ...) {
-  if (is.loo(x)) {
-    dots <- list(...)
-    loos <- c(list(x), dots)
-  } else {
-    if (!is.list(x) || !length(x)) {
-      stop("'x' must be a list if not a 'loo' object.")
-    }
-    if (length(list(...))) {
-      stop("If 'x' is a list then '...' should not be specified.")
-    }
-    loos <- x
-  }
+loo_compare.default <- function(x, ..., rank_by = NULL) {
+  loos <- .loo_compare_inputs(x, ...)
 
   # if subsampling is used
   if (any(sapply(loos, inherits, "psis_loo_ss"))) {
     return(loo_compare.psis_loo_ss_list(loos))
   }
 
+  if (all(vapply(loos, is.loo_pred_measure, logical(1)))) {
+    return(compare_loo_pred_measure(loos, rank_by = rank_by))
+  }
+
+  if (any(vapply(loos, is.loo_pred_measure, logical(1)))) {
+    stop(
+      "Cannot mix 'loo_pred_measure' objects with other 'loo' objects. ",
+      "Compare models using 'loo_pred_measure()' for each model.",
+      call. = FALSE
+    )
+  }
+
+  if (any(vapply(loos, inherits, what = "pred_measure", logical(1)))) {
+    stop(
+      "'loo_compare' for predictive measures requires 'loo_pred_measure' objects. ",
+      "Use loo_pred_measure() instead of insample_pred_measure(), ",
+      "kfold_pred_measure(), or test_pred_measure().",
+      call. = FALSE
+    )
+  }
+
+  if (!is.null(rank_by)) {
+    warning(
+      "`rank_by` is only used for `loo_pred_measure` comparisons and will be ignored.",
+      call. = FALSE
+    )
+  }
+
   # run pre-comparison checks
   loo_compare_checks(loos)
 
   # compute elpd_diff and se_elpd_diff relative to best model
-  comp <- loo_compare_matrix(loos)
   ord <- loo_compare_order(loos)
+  comp <- loo_compare_matrix(loos, ord = ord)
   rnms <- rownames(comp)
   diffs <- mapply(FUN = elpd_diffs, loos[ord[1]], loos[ord])
   colnames(diffs) <- rnms
@@ -170,13 +268,31 @@ loo_compare.default <- function(x, ...) {
 #' @param p_worse For the print method only, should we include the normal
 #'   approximation based probability of each model having worse performance than
 #'   the best model? The default is `TRUE`.
-print.compare.loo <- function(x, ..., digits = 1, p_worse = TRUE) {
+#' @param measures For `loo_pred_measure` comparisons only, which measures to
+#'   print diff tables for. `NULL` (default) prints only the ranking measure
+#'   (`"elpd"` when `rank_by` was not set, otherwise `rank_by`);
+#'   `"all"` prints all compared measures; or a character vector of measure
+#'   names (e.g. `c("elpd", "mse")`). Printed tables use the column name
+#'   `se_diff` for the standard error of the difference even when the data
+#'   frame column is `{measure}_se_diff`.
+print.compare.loo <- function(x, ..., digits = 1, p_worse = TRUE, measures = NULL) {
   if (inherits(x, "old_compare.loo")) {
     return(unclass(x))
   }
   if (!inherits(x, "data.frame")) {
     class(x) <- c(class(x), "data.frame")
   }
+
+  compare_measures <- attr(x, "compare_measures")
+  if (!is.null(compare_measures)) {
+    return(.print_compare_loo_pred_measure(
+      x,
+      digits = digits,
+      p_worse = p_worse,
+      measures = measures
+    ))
+  }
+
   if (!all(c("model", "elpd_diff", "se_diff") %in% colnames(x))) {
     print(as.data.frame(x))
     return(x)
@@ -195,9 +311,155 @@ print.compare.loo <- function(x, ..., digits = 1, p_worse = TRUE) {
   }
   print(x2, quote = FALSE, row.names = FALSE)
 
-  # show glossary for diagnostic flags
-  has_diag <- any(nzchar(x[["diag_diff"]], keepNA = FALSE), na.rm = TRUE) ||
-              any(nzchar(x[["diag_elpd"]], keepNA = FALSE), na.rm = TRUE)
+  .print_compare_diag_message(x, p_worse = p_worse)
+  invisible(x)
+}
+
+
+# internal ----------------------------------------------------------------
+
+#' Print `compare.loo` results from `loo_pred_measure` comparisons
+#' @noRd
+.print_compare_loo_pred_measure <- function(x, digits, p_worse, measures) {
+  rank_by <- attr(x, "rank_by")
+  compare_measures <- attr(x, "compare_measures")
+  ref_model <- x$model[[1L]]
+  primary_measure <- if (is.null(rank_by)) "elpd" else rank_by
+
+  measures_to_print <- if (is.null(measures)) {
+    primary_measure
+  } else if (identical(measures, "all")) {
+    compare_measures
+  } else {
+    measures
+  }
+
+  unknown <- setdiff(measures_to_print, compare_measures)
+  if (length(unknown)) {
+    stop(
+      paste0(
+        "Unknown measure(s) in `measures`: ",
+        paste(unknown, collapse = ", "),
+        ". Available measures: ",
+        paste(compare_measures, collapse = ", ")
+      ),
+      call. = FALSE
+    )
+  }
+
+  if (identical(measures, "all") && length(compare_measures) > 4L) {
+    message(
+      "Printing ", length(compare_measures), " measure comparisons; ",
+      "consider `measures = c(...)`."
+    )
+  }
+
+  if (is.null(measures) && !is.null(rank_by)) {
+    message(
+      "Models ranked by ", rank_by, " (reference: ", ref_model, ")."
+    )
+  }
+
+  show_diag_elpd_primary <- is.null(measures) || identical(measures, "all")
+  for (i in seq_along(measures_to_print)) {
+    measure <- measures_to_print[[i]]
+    if (!is.null(measures)) {
+      cat("\n-- ", measure, " (vs ", ref_model, ") --\n", sep = "")
+    }
+    .print_compare_measure_table(
+      x,
+      measure = measure,
+      digits = digits,
+      p_worse = p_worse,
+      show_diag_elpd = show_diag_elpd_primary &&
+        identical(measure, primary_measure) &&
+        i == match(primary_measure, measures_to_print)
+    )
+  }
+
+  has_diag_msg <- .print_compare_diag_message(
+    x,
+    p_worse = p_worse,
+    measures = measures_to_print
+  )
+
+  if (is.null(measures)) {
+    other <- setdiff(compare_measures, primary_measure)
+    if (length(other)) {
+      message(
+        if (has_diag_msg) "\n",
+        "Other measures compared: ",
+        paste(other, collapse = ", "),
+        ". Use print(x, measures = \"all\")."
+      )
+    }
+  }
+
+  invisible(x)
+}
+
+#' Print one measure's comparison table
+#' @noRd
+.print_compare_measure_table <- function(x, measure, digits, p_worse, show_diag_elpd) {
+  if (.is_elpd_measure(measure)) {
+    diff_col <- "elpd_diff"
+    se_col <- "se_diff"
+    diff_name <- "elpd_diff"
+    se_name <- "se_diff"
+  } else {
+    diff_col <- paste0(measure, "_diff")
+    se_col <- paste0(measure, "_se_diff")
+    diff_name <- diff_col
+    se_name <- "se_diff"
+  }
+
+  if (!all(c(diff_col, se_col) %in% colnames(x))) {
+    stop(
+      "Comparison columns for measure '", measure, "' are missing.",
+      call. = FALSE
+    )
+  }
+
+  x2 <- data.frame(
+    model = x$model,
+    diff = unname(.fr(x[[diff_col]], digits)),
+    se_diff = unname(.fr(x[[se_col]], digits)),
+    check.names = FALSE,
+    stringsAsFactors = FALSE
+  )
+  names(x2)[2:3] <- c(diff_name, se_name)
+
+  if (.is_elpd_measure(measure) && p_worse && "p_worse" %in% colnames(x)) {
+    x2$p_worse <- unname(.fr(x[["p_worse"]], digits = 2))
+    x2$diag_diff <- x[["diag_diff"]]
+  }
+  if (show_diag_elpd && "diag_elpd" %in% colnames(x)) {
+    x2$diag_elpd <- x[["diag_elpd"]]
+  }
+
+  print(x2, quote = FALSE, row.names = FALSE)
+}
+
+#' Print diagnostic glossary message for compare output
+#' @noRd
+.print_compare_diag_message <- function(x, p_worse, measures = NULL) {
+  diag_cols <- c("diag_elpd")
+  if (is.null(measures) || "elpd" %in% measures) {
+    diag_cols <- c("diag_diff", diag_cols)
+  } else if (!is.null(measures)) {
+    elpd_measures <- measures[vapply(measures, .is_elpd_measure, logical(1))]
+    if (length(elpd_measures)) {
+      diag_cols <- c("diag_diff", diag_cols)
+    }
+  }
+
+  has_diag <- any(
+    vapply(
+      intersect(diag_cols, colnames(x)),
+      function(col) any(nzchar(x[[col]], keepNA = FALSE), na.rm = TRUE),
+      logical(1)
+    )
+  )
   if (has_diag && p_worse) {
     message(
       "\nDiagnostic flags present.\n",
@@ -205,11 +467,522 @@ print.compare.loo <- function(x, ..., digits = 1, p_worse = TRUE) {
       "or https://mc-stan.org/loo/reference/loo-glossary.html."
     )
   }
-  invisible(x)
+  invisible(has_diag && p_worse)
+}
+
+#' Is an object a PSIS-LOO predictive measure result?
+#' @noRd
+is.loo_pred_measure <- function(x) {
+  inherits(x, "loo_pred_measure")
 }
 
+#' Normalize `loo_compare()` inputs to a list of model results
+#' @noRd
+.loo_compare_inputs <- function(x, ...) {
+  if (is.loo(x) || inherits(x, "pred_measure")) {
+    dots <- list(...)
+    return(c(list(x), dots))
+  }
+  if (!is.list(x) || !length(x)) {
+    stop(
+      "'x' must be a list if not a 'loo' or 'pred_measure' object.",
+      call. = FALSE
+    )
+  }
+  if (length(list(...))) {
+    stop("If 'x' is a list then '...' should not be specified.", call. = FALSE)
+  }
+  x
+}
 
-# internal ----------------------------------------------------------------
+#' Compare `loo_pred_measure` objects (multi-measure path)
+#' @noRd
+#' @param loos List of `loo_pred_measure` objects.
+#' @param rank_by Bare measure name used to order models.
+compare_loo_pred_measure <- function(loos, rank_by = NULL) {
+  loo_compare_checks(
+    loos,
+    class_check = is.loo_pred_measure,
+    class_msg = "All inputs must have class 'loo_pred_measure'.",
+    kfold_checks = FALSE
+  )
+  .compare_metadata_check(loos)
+  .warn_omitted_compare_measures(loos)
+
+  rank_measure <- .resolve_rank_measure(loos, rank_by)
+  compare_cols <- .compare_pointwise_cols(loos)
+  .inform_compare_sign_conversion(compare_cols, loos)
+  ord <- loo_compare_order(loos, rank_measure$internal)
+  loos_ord <- loos[ord]
+  ref_loo <- loos_ord[[1L]]
+
+  comp <- loo_compare_matrix(
+    loos_ord,
+    bare_names = TRUE,
+    ord = seq_along(loos_ord)
+  )
+  rnms <- rownames(comp)
+  n_obs <- nrow(loos_ord[[1L]]$pointwise)
+
+  diff_cols <- list()
+  for (col in compare_cols) {
+    bare <- .display_name(col)
+    method <- .measure_pointwise_diff_method(loos_ord, col)
+    pair_stats <- vapply(
+      loos_ord,
+      .pair_measure_stats,
+      FUN.VALUE = c(diff = 0, se = 0),
+      ref = ref_loo,
+      col = col,
+      method = method,
+      loos = loos_ord
+    )
+    measure_diff <- pair_stats["diff", ]
+    measure_se <- pair_stats["se", ]
+
+    if (.is_elpd_measure(col)) {
+      diff_cols$elpd_diff <- measure_diff
+      diff_cols$se_diff <- measure_se
+      p_worse <- stats::pnorm(0, measure_diff, measure_se)
+      p_worse[measure_diff == 0] <- NA_real_
+      diff_cols$p_worse <- p_worse
+      diff_cols$diag_diff <- diag_diff(n_obs, measure_diff)
+    } else {
+      diff_cols[[paste0(bare, "_diff")]] <- measure_diff
+      diff_cols[[paste0(bare, "_se_diff")]] <- measure_se
+    }
+  }
+
+  comp <- cbind(
+    data.frame(
+      model = rnms,
+      diff_cols,
+      diag_elpd = diag_elpd(loos_ord),
+      stringsAsFactors = FALSE
+    ),
+    as.data.frame(comp)
+  )
+  rownames(comp) <- NULL
+
+  loo_order_stat_check(
+    loos_ord,
+    seq_along(loos_ord),
+    rank_col = rank_measure$internal
+  )
+
+  if (!is.null(rank_by)) {
+    attr(comp, "rank_by") <- rank_measure$bare
+  }
+  attr(comp, "compare_measures") <- .compare_measures(loos)
+  attr(comp, "sign_converted_measures") <- .compare_sign_converted_measures(
+    compare_cols,
+    loos
+  )
+  class(comp) <- c("compare.loo", class(comp))
+  comp
+}
+
+#' Strip `_loo` suffix for `loo_compare` display names
+#' @noRd
+.display_name <- function(col) {
+  sub("_loo$", "", col)
+}
+
+#' Map bare measure name to `pointwise` column name
+#' @noRd
+.pointwise_col <- function(name, cols) {
+  if (name %in% cols) {
+    return(name)
+  }
+  internal <- paste0(name, "_loo")
+  if (internal %in% cols) {
+    return(internal)
+  }
+  stop(
+    paste0(
+      "Measure '", name, "' not found in all models. ",
+      "Available measures: ",
+      paste(vapply(cols, .display_name, character(1)), collapse = ", ")
+    ),
+    call. = FALSE
+  )
+}
+
+#' Common `pointwise` columns across models, excluding complexity terms
+#' @noRd
+.compare_pointwise_cols <- function(loos) {
+  cols <- Reduce(
+    intersect,
+    lapply(loos, function(x) colnames(x$pointwise))
+  )
+  cols[!grepl("^p_", cols)]
+}
+
+#' Check that comparison metadata is consistent across models
+#' @noRd
+.compare_metadata_check <- function(loos) {
+  bare_measures <- .compare_measures(loos)
+  if (!length(bare_measures)) {
+    return(invisible(NULL))
+  }
+
+  for (bare in bare_measures) {
+    metas <- lapply(loos, function(x) {
+      compare_meta <- attr(x, "measure_compare_meta")
+      if (is.null(compare_meta)) {
+        return(NULL)
+      }
+      compare_meta[[bare]]
+    })
+    has_meta <- !vapply(metas, is.null, logical(1))
+    if (any(has_meta) && !all(has_meta)) {
+      stop(
+        "Not all models provide comparison metadata for measure '",
+        bare,
+        "'. Recompute all inputs with the current version of `loo_pred_measure()`.",
+        call. = FALSE
+      )
+    }
+    # `extra` holds per-measure auxiliary data (for `r2`, the pointwise
+    # baseline derived from `y`), which legitimately differs when models are
+    # fitted to different data. That case is already reported by the `yhash`
+    # warning, so comparing `extra` here would only mislabel it as a
+    # `higher_is_better` disagreement.
+    non_null <- lapply(metas[has_meta], function(meta) {
+      meta$extra <- NULL
+      meta
+    })
+    if (length(non_null) > 1L) {
+      ref <- non_null[[1L]]
+      inconsistent <- vapply(
+        non_null[-1L],
+        function(meta) !identical(meta, ref),
+        logical(1)
+      )
+      if (any(inconsistent)) {
+        stop(
+          "Models disagree on comparison metadata for measure '",
+          bare,
+          "'. Ensure all models use the same `higher_is_better` settings for each measure.",
+          call. = FALSE
+        )
+      }
+    }
+  }
+
+  invisible(NULL)
+}
+
+#' Warn when models do not share the same predictive measures
+#' @noRd
+.warn_omitted_compare_measures <- function(loos) {
+  model_names <- find_model_names(loos)
+  if (anyDuplicated(model_names)) {
+    model_names <- make.unique(model_names, sep = "_")
+  }
+  by_model <- stats::setNames(
+    lapply(loos, function(x) {
+      cols <- colnames(x$pointwise)
+      cols <- cols[!grepl("^p_", cols)]
+      unname(vapply(cols, .display_name, character(1)))
+    }),
+    model_names
+  )
+  common <- Reduce(intersect, by_model)
+  omitted <- setdiff(unique(unlist(by_model)), common)
+  if (!length(omitted)) {
+    return(invisible(NULL))
+  }
+  omitted <- sort(omitted)
+
+  omitted_detail <- vapply(
+    omitted,
+    function(measure) {
+      present <- names(by_model)[vapply(
+        by_model,
+        function(measures) measure %in% measures,
+        logical(1)
+      )]
+      paste0(measure, " (", paste(present, collapse = ", "), ")")
+    },
+    character(1)
+  )
+
+  warning(
+    paste0(
+      "Omitted measures: ",
+      paste(omitted_detail, collapse = ", "),
+      ". Compared: ",
+      paste(common, collapse = ", "),
+      "."
+    ),
+    call. = FALSE
+  )
+}
+
+#' Bare measure names available for comparison across models
+#' @noRd
+.compare_measures <- function(loos) {
+  cols <- .compare_pointwise_cols(loos)
+  unname(vapply(cols, .display_name, character(1)))
+}
+
+#' Resolve `rank_by` to bare and internal `pointwise` column names
+#' @noRd
+.resolve_rank_measure <- function(loos, rank_by = NULL) {
+  cols <- .compare_pointwise_cols(loos)
+  bare <- if (is.null(rank_by)) "elpd" else rank_by
+  internal <- .pointwise_col(bare, cols)
+  list(
+    bare = .display_name(internal),
+    internal = internal
+  )
+}
+
+#' Is a measure an ELPD-family measure (for `p_worse` / `diag_diff`)?
+#' @noRd
+.is_elpd_measure <- function(name) {
+  grepl("^elpd", .display_name(name))
+}
+
+#' Look up per-measure comparison metadata on a result object
+#' @noRd
+.get_measure_compare_meta <- function(loos, bare) {
+  compare_meta <- attr(loos[[1L]], "measure_compare_meta")
+  if (is.null(compare_meta)) {
+    return(NULL)
+  }
+  compare_meta[[bare]]
+}
+
+#' Whether a measure is intrinsically a loss (natural scale: lower is better)
+#'
+#' Unlike `.measure_lower_is_better()` this ignores `higher_is_better`, so it
+#' describes the measure itself rather than the scale its values are stored on.
+#' @noRd
+.measure_is_loss <- function(name, loos = NULL) {
+  bare <- .display_name(name)
+
+  if (!is.null(loos)) {
+    meta <- .get_measure_compare_meta(loos, bare)
+    if (!is.null(meta) && !is.null(meta$loss)) {
+      return(isTRUE(meta$loss))
+    }
+  }
+
+  spec <- .measure_spec[[bare]]
+  if (!is.null(spec)) {
+    return(isTRUE(spec$loss))
+  }
+  bare %in% c("ic", "mae", "mse", "rmse", "brier", "srps")
+}
+
+#' Whether stored values are on a loss scale (lower is better)
+#' @noRd
+.measure_lower_is_better <- function(name, loos = NULL) {
+  bare <- .display_name(name)
+  higher_is_better <- NULL
+
+  if (!is.null(loos)) {
+    meta <- .get_measure_compare_meta(loos, bare)
+    if (!is.null(meta)) {
+      higher_is_better <- meta$higher_is_better
+    } else {
+      hib_attr <- attr(loos[[1L]], "measure_higher_is_better")
+      if (!is.null(hib_attr) && bare %in% names(hib_attr)) {
+        higher_is_better <- hib_attr[[bare]]
+      }
+    }
+  }
+
+  if (!is.null(higher_is_better)) {
+    return(!isTRUE(higher_is_better))
+  }
+
+  .measure_is_loss(name, loos)
+}
+
+#' Sign converting stored measure values to the measure's natural scale
+#'
+#' `higher_is_better` may have negated the stored values (see
+#' `.create_measure_structure()`). Delta-method standard errors are derived on
+#' the natural scale (e.g. RMSE positive), so they must be undone first.
+#' @noRd
+.measure_natural_sign <- function(name, loos = NULL) {
+  stored_lower_is_better <- .measure_lower_is_better(name, loos)
+  if (identical(stored_lower_is_better, .measure_is_loss(name, loos))) 1 else -1
+}
+
+#' Bare names of measures whose sign is flipped for `loo_compare()`
+#' @noRd
+.compare_sign_converted_measures <- function(cols, loos) {
+  bare <- vapply(cols, .display_name, character(1))
+  unique(bare[vapply(
+    cols,
+    function(col) .measure_lower_is_better(col, loos),
+    logical(1)
+  )])
+}
+
+#' Inform when measure signs are flipped for comparison
+#' @noRd
+.inform_compare_sign_conversion <- function(cols, loos) {
+  converted <- .compare_sign_converted_measures(cols, loos)
+  if (!length(converted)) {
+    return(invisible(NULL))
+  }
+  message(
+    "For model comparison, differences for ",
+    paste(converted, collapse = ", "),
+    " ",
+    if (length(converted) == 1L) "is" else "are",
+    "\nreported on a utility scale (higher is better)."
+  )
+  invisible(NULL)
+}
+
+#' How to aggregate paired pointwise differences for a measure
+#'
+#' Returns `"sum"` when the overall estimate equals the sum of pointwise
+#' contributions, `"mean"` when it equals the mean, `"pairwise"` when the
+#' measure supplies its own `se_diff_fun`, and `"estimates_only"` when pointwise
+#' values do not define the overall estimate and no `se_diff_fun` is available.
+#' @noRd
+.measure_pointwise_diff_method <- function(loos, col) {
+  bare <- .display_name(col)
+  meta <- .get_measure_compare_meta(loos, bare)
+  if (!is.null(meta) && !identical(meta$diff_method, "auto")) {
+    return(meta$diff_method)
+  }
+
+  if (.is_elpd_measure(col) || bare == "ic") {
+    return("sum")
+  }
+
+  ref <- loos[[1L]]
+  est <- ref$estimates[col, "Estimate"]
+  pw <- ref$pointwise[, col, drop = TRUE]
+  if (!length(pw) || !is.finite(est)) {
+    return("estimates_only")
+  }
+
+  tol <- sqrt(.Machine$double.eps) * max(abs(c(est, pw)), na.rm = TRUE)
+  if (isTRUE(all.equal(est, sum(pw), tolerance = tol, check.attributes = FALSE))) {
+    return("sum")
+  }
+  if (isTRUE(all.equal(est, mean(pw), tolerance = tol, check.attributes = FALSE))) {
+    return("mean")
+  }
+  "estimates_only"
+}
+
+#' Resolve a measure's `se_diff_fun`
+#'
+#' Built-in measures name an entry of `.se_diff_funs`; custom measures store the
+#' function itself (see `.measure_compare_meta()`).
+#' @noRd
+.measure_se_diff_fun <- function(loos, col) {
+  bare <- .display_name(col)
+  meta <- .get_measure_compare_meta(loos, bare)
+
+  fun <- meta$se_diff_fun
+  if (is.null(fun)) {
+    fun <- .measure_spec[[bare]]$se_diff_fun
+  }
+  if (is.character(fun)) {
+    fun <- .se_diff_funs[[fun]]
+  }
+  if (!is.function(fun)) {
+    stop(
+      paste0(
+        "No 'se_diff_fun' available for measure '", bare, "'."
+      ),
+      call. = FALSE
+    )
+  }
+  fun
+}
+
+#' Assemble one model's inputs for an `se_diff_fun`
+#'
+#' Every element describes the single model `x`, including `extra`, which is
+#' read from that model's own comparison metadata rather than the reference
+#' model's.
+#' @noRd
+#' @param sgn Sign restoring the measure's natural scale, see
+#'   `.measure_natural_sign()`.
+.se_diff_input <- function(x, col, sgn) {
+  list(
+    estimate = sgn * x$estimates[col, "Estimate"],
+    se = x$estimates[col, "SE"],
+    pointwise = sgn * x$pointwise[, col, drop = TRUE],
+    extra = .get_measure_compare_meta(list(x), .display_name(col))$extra
+  )
+}
+
+#' Validate the value returned by an `se_diff_fun`
+#' @noRd
+.validate_se_diff <- function(se, col) {
+  if (!is.numeric(se) || length(se) != 1L) {
+    stop(
+      paste0(
+        "The 'se_diff_fun' for measure '", .display_name(col),
+        "' must return a numeric scalar."
+      ),
+      call. = FALSE
+    )
+  }
+  unname(se)
+}
+
+#' Paired measure difference and SE for one model vs a reference
+#' @noRd
+.pair_measure_stats <- function(cmp, ref, col, method = NULL, loos = list(ref)) {
+  if (is.null(method)) {
+    method <- .measure_pointwise_diff_method(c(list(ref, cmp)), col)
+  }
+
+  flip <- .measure_lower_is_better(col, loos)
+  est_utility <- function(estimates) {
+    val <- estimates[col, "Estimate"]
+    if (flip) -val else val
+  }
+
+  if (method == "estimates_only") {
+    return(c(
+      diff = est_utility(cmp$estimates) - est_utility(ref$estimates),
+      se = NA_real_
+    ))
+  }
+
+  if (method == "pairwise") {
+    se_diff_fun <- .measure_se_diff_fun(loos, col)
+    sgn <- .measure_natural_sign(col, loos)
+    se <- se_diff_fun(
+      ref = .se_diff_input(ref, col, sgn),
+      cmp = .se_diff_input(cmp, col, sgn)
+    )
+    return(c(
+      diff = est_utility(cmp$estimates) - est_utility(ref$estimates),
+      se = .validate_se_diff(se, col)
+    ))
+  }
+
+  to_utility <- function(pointwise) {
+    x <- pointwise[, col, drop = TRUE]
+    if (flip) -x else x
+  }
+  diffs <- to_utility(cmp$pointwise) - to_utility(ref$pointwise)
+
+  diff <- if (method == "sum") sum(diffs) else mean(diffs)
+  se <- if (method == "sum") {
+    se_elpd_diff(diffs)
+  } else {
+    N <- length(diffs)
+    if (N <= 1L) 0 else stats::sd(diffs) / sqrt(N)
+  }
+  c(diff = diff, se = se)
+}
 
 #' Compute pointwise elpd differences
 #' @noRd
@@ -235,14 +1008,22 @@ se_elpd_diff <- function(diffs) {
 #' Perform checks on `"loo"` objects before comparison
 #' @noRd
 #' @param loos List of `"loo"` objects.
+#' @param class_check Function returning `TRUE` for valid input objects.
+#' @param class_msg Error message when `class_check` fails.
+#' @param kfold_checks If `TRUE`, run k-fold comparison warnings.
 #' @return Nothing, just possibly throws errors/warnings.
-loo_compare_checks <- function(loos) {
+loo_compare_checks <- function(
+  loos,
+  class_check = is.loo,
+  class_msg = "All inputs should have class 'loo'.",
+  kfold_checks = TRUE
+) {
   ## errors
   if (length(loos) <= 1L) {
-    stop("'loo_compare' requires at least two models.", call.=FALSE)
+    stop("'loo_compare' requires at least two models.", call. = FALSE)
   }
-  if (!all(sapply(loos, is.loo))) {
-    stop("All inputs should have class 'loo'.", call.=FALSE)
+  if (!all(vapply(loos, class_check, logical(1)))) {
+    stop(class_msg, call. = FALSE)
   }
 
   Ns <- vapply(loos, function(x) nrow(x$pointwise), integer(1))
@@ -259,26 +1040,37 @@ loo_compare_checks <- function(loos) {
   ## warnings
 
   yhash <- lapply(loos, attr, which = "yhash")
-  yhash_ok <- sapply(yhash, function(x) { # ok only if all yhash are same (all NULL is ok)
+  yhash_ok <- vapply(yhash, function(x) {
     isTRUE(all.equal(x, yhash[[1]]))
-  })
+  }, logical(1))
   if (!all(yhash_ok)) {
-    warning("Not all models have the same y variable. ('yhash' attributes do not match)",
-            call. = FALSE)
+    warning(
+      "Not all models have the same y variable. ('yhash' attributes do not match)",
+      call. = FALSE
+    )
+  }
+
+  if (!kfold_checks) {
+    return(invisible(NULL))
   }
 
-  if (all(sapply(loos, is.kfold))) {
+  if (all(vapply(loos, is.kfold, logical(1)))) {
     Ks <- unlist(lapply(loos, attr, which = "K"))
     if (!all(Ks == Ks[1])) {
-      warning("Not all kfold objects have the same K value. ",
-              "For a more accurate comparison use the same number of folds. ",
-              call. = FALSE)
+      warning(
+        "Not all kfold objects have the same K value. ",
+        "For a more accurate comparison use the same number of folds. ",
+        call. = FALSE
+      )
     }
-  } else if (any(sapply(loos, is.kfold)) && any(sapply(loos, is.psis_loo))) {
-    warning("Comparing LOO-CV to K-fold-CV. ",
-            "For a more accurate comparison use the same number of folds ",
-            "or loo for all models compared.",
-            call. = FALSE)
+  } else if (any(vapply(loos, is.kfold, logical(1))) &&
+      any(vapply(loos, is.psis_loo, logical(1)))) {
+    warning(
+      "Comparing LOO-CV to K-fold-CV. ",
+      "For a more accurate comparison use the same number of folds ",
+      "or loo for all models compared.",
+      call. = FALSE
+    )
   }
 }
 
@@ -313,46 +1105,79 @@ find_model_names <- function(x) {
 }
 
 
-#' Compute the loo_compare matrix
+#' Build estimates table for `loo_compare()` ordering and matrix output
 #' @noRd
-#' @param loos List of `"loo"` objects.
-loo_compare_matrix <- function(loos){
-  tmp <- sapply(loos, function(x) {
+.loo_compare_estimates_table <- function(loos, bare_names = FALSE) {
+  sapply(loos, function(x) {
     est <- x$estimates
-    setNames(c(est), nm = c(rownames(est), paste0("se_", rownames(est))))
+    rows <- if (bare_names) .display_name(rownames(est)) else rownames(est)
+    setNames(c(est), nm = c(rows, paste0("se_", rows)))
   })
+}
+
+#' Compute the loo_compare matrix
+#' @noRd
+#' @param loos List of `"loo"` objects.
+#' @param bare_names If `TRUE`, strip `_loo` suffixes from estimate row names.
+#' @param ord Optional model ordering indices; computed from ELPD when `NULL`.
+loo_compare_matrix <- function(loos, bare_names = FALSE, ord = NULL) {
+  tmp <- .loo_compare_estimates_table(loos, bare_names = bare_names)
   colnames(tmp) <- find_model_names(loos)
-  rnms <- rownames(tmp)
-  comp <- tmp
-  ord <- loo_compare_order(loos)
-  comp <- t(comp)[ord, ]
-  patts <- c("elpd", "p_", "^waic$|^looic$", "^se_waic$|^se_looic$")
-  col_ord <- unlist(sapply(patts, function(p) grep(p, colnames(comp))),
-                    use.names = FALSE)
-  comp <- comp[, col_ord]
+  comp <- t(tmp)
+
+  if (is.null(ord)) {
+    ord <- loo_compare_order(loos)
+  }
+  comp <- comp[ord, , drop = FALSE]
+
+  patts <- if (bare_names) {
+    c("^elpd$", "^p$", "^se_elpd$", "^se_p$")
+  } else {
+    c("elpd", "p_", "^waic$|^looic$", "^se_waic$|^se_looic$")
+  }
+  col_ord <- unique(unlist(
+    lapply(patts, function(p) grep(p, colnames(comp))),
+    use.names = FALSE
+  ))
+  if (bare_names) {
+    other <- setdiff(seq_len(ncol(comp)), col_ord)
+    comp <- comp[, c(col_ord, other), drop = FALSE]
+  } else {
+    comp <- comp[, col_ord, drop = FALSE]
+  }
   comp
 }
 
 #' Computes the order of loos for comparison
 #' @noRd
 #' @param loos List of `"loo"` objects.
-loo_compare_order <- function(loos){
-  tmp <- sapply(loos, function(x) {
-    est <- x$estimates
-    setNames(c(est), nm = c(rownames(est), paste0("se_", rownames(est))))
-  })
-  colnames(tmp) <- find_model_names(loos)
-  rnms <- rownames(tmp)
-  ord <- order(tmp[grep("^elpd", rnms), ], decreasing = TRUE)
-  ord
+#' @param rank_col Optional internal `pointwise` column name used for ranking.
+loo_compare_order <- function(loos, rank_col = NULL) {
+  if (is.null(rank_col)) {
+    tmp <- .loo_compare_estimates_table(loos, bare_names = FALSE)
+    colnames(tmp) <- find_model_names(loos)
+    rnms <- rownames(tmp)
+    return(order(tmp[grep("^elpd", rnms), ], decreasing = TRUE))
+  }
+
+  est_row <- vapply(loos, function(x) {
+    val <- x$estimates[rank_col, "Estimate"]
+    if (.measure_lower_is_better(rank_col, loos)) -val else val
+  }, numeric(1))
+  order(est_row, decreasing = TRUE)
 }
 
 #' Perform checks on `"loo"` objects __after__ comparison
 #' @noRd
 #' @param loos List of `"loo"` objects.
 #' @param ord List of `"loo"` object orderings.
+#' @param measure_diff Optional precomputed model differences for the rank
+#'   measure; computed from the median model when `NULL`.
+#' @param rank_col Optional internal `pointwise` column name used for the
+#'   median-baseline differences when `measure_diff` is `NULL` and inputs are not
+#'   classic `"loo"` objects.
 #' @return Nothing, just possibly throws errors/warnings.
-loo_order_stat_check <- function(loos, ord) {
+loo_order_stat_check <- function(loos, ord, measure_diff = NULL, rank_col = NULL) {
 
   ## breaks
 
@@ -360,32 +1185,49 @@ loo_order_stat_check <- function(loos, ord) {
     # procedure cannot be diagnosed for fewer than ten candidate models
     # (total models = worst model + ten candidates)
     # break from function
-    return(NULL)
+    return(invisible(NULL))
   }
 
   ## warnings
 
-  # compute the elpd differences from the median model
-  baseline_idx <- middle_idx(ord)
-  diffs <- mapply(FUN = elpd_diffs, loos[ord[baseline_idx]], loos[ord])
-  elpd_diff <- apply(diffs, 2, sum)
+  if (is.null(measure_diff)) {
+    # compute differences from the median model
+    baseline_idx <- middle_idx(ord)
+    ref_loo <- loos[[ord[baseline_idx]]]
+    if (is.null(rank_col)) {
+      diffs <- mapply(FUN = elpd_diffs, loos[ord[baseline_idx]], loos[ord])
+      measure_diff <- apply(diffs, 2, sum)
+    } else {
+      method <- .measure_pointwise_diff_method(loos, rank_col)
+      measure_diff <- vapply(
+        loos[ord],
+        .pair_measure_stats,
+        FUN.VALUE = c(diff = 0, se = 0),
+        ref = ref_loo,
+        col = rank_col,
+        method = method,
+        loos = loos
+      )["diff", ]
+    }
+  }
 
   # estimate the standard deviation of the upper-half-normal
-  diff_median <- stats::median(elpd_diff)
-  elpd_diff_trunc <- elpd_diff[elpd_diff >= diff_median]
-  n_models <- sum(!is.na(elpd_diff_trunc))
-  candidate_sd <- sqrt(1 / n_models * sum(elpd_diff_trunc^2, na.rm = TRUE))
+  diff_median <- stats::median(measure_diff)
+  measure_diff_trunc <- measure_diff[measure_diff >= diff_median]
+  n_models <- sum(!is.na(measure_diff_trunc))
+  candidate_sd <- sqrt(1 / n_models * sum(measure_diff_trunc^2, na.rm = TRUE))
 
   # estimate expected best diff under null hypothesis
   K <- length(loos) - 1
   order_stat <- order_stat_heuristic(K, candidate_sd)
 
-  if (max(elpd_diff) <= order_stat) {
+  if (max(measure_diff) <= order_stat) {
     # flag warning if we suspect no model is theoretically better than the baseline
     warning("Difference in performance potentially due to chance. ",
             "See McLatchie and Vehtari (2023) for details.",
             call. = FALSE)
   }
+  invisible(NULL)
 }
 
 #' Returns the middle index of a vector
diff --git a/R/pred_measure-builtin.R b/R/pred_measure-builtin.R
index fda28a32..760bf6a1 100644
--- a/R/pred_measure-builtin.R
+++ b/R/pred_measure-builtin.R
@@ -443,14 +443,25 @@ measure_bacc <- function(
     acc_i <- (mupred_hat == y) * 1L
   }
   
+  # recomputed rather than reused from the branch above, which the
+  # precomputed-`pointwise` branch never enters
+  classes <- sort(unique(y))
+  K <- length(classes)
+  class_id <- match(y, classes)
+  n_c <- tabulate(class_id, nbins = K)
+
   acc_c <- vapply(classes, function(c) mean(acc_i[y == c]), numeric(1))
-  n_c <- tabulate(match(y, classes))
-  bacc_i <- acc_i / (K * n_c[match(y, classes)])
-  
+  bacc_i <- acc_i / (K * n_c[class_id])
+
   res <- list(
     estimate = mean(acc_c),
     se = sqrt(sum(acc_c * (1 - acc_c) / n_c)) / K,
-    pointwise = bacc_i
+    pointwise = bacc_i,
+    # the class strata are what makes a difference of balanced accuracies more
+    # than a mean of pointwise differences, and `.se_diff_bacc()` cannot
+    # recover them from `bacc_i`: every misclassified observation stores a
+    # zero whatever its class
+    extra = list(class_id = class_id)
   )
   .create_measure_structure(
     res, higher_is_better, "bacc", n_draws = n_draws, n_obs = n_obs
@@ -685,6 +696,33 @@ measure_rmse <- function(
   )
 }
 
+#' Delta-method standard error of an R-squared quantity
+#'
+#' The R-squared of one model and the difference in R-squared between two
+#' models have the same form: a mean of squared-error contributions divided by
+#' the model-independent baseline `MSE(y)`. Their standard errors are therefore
+#' the same first-order Taylor approximation, evaluated either at one model's
+#' pointwise squared errors or at the pointwise differences between two models.
+#'
+#' Writing `c = mean(sqe) / MSE(y)`, the three-term expansion
+#' `Var[MSE] - 2 * c * Cov[MSE, MSE(y)] + c^2 * Var[MSE(y)]`, scaled by
+#' `MSE(y)^-1`, is exactly the standard error of the mean of
+#' `sqe_i - c * mse_y_i`. That is the form used here: it needs one variance
+#' rather than three moments, it cannot go negative under the square root, and
+#' it is exactly `0` when `sqe` is identically zero, as it is when a model is
+#' compared against itself.
+#'
+#' @noRd
+#' @param sqe Pointwise squared errors of one model, or pointwise differences
+#'   in squared error between two models.
+#' @param mse_y_i Pointwise baseline `(y_i - mean(y))^2`.
+#' @return Numeric scalar standard error.
+.se_r2_delta <- function(sqe, mse_y_i) {
+  mse_y_hat <- mean(mse_y_i)
+  scaled <- sqe - (mean(sqe) / mse_y_hat) * mse_y_i
+  sqrt(var(scaled) / length(sqe)) / mse_y_hat
+}
+
 #' Predictive R-squared (`r2`)
 #'
 #' Computes predictive R-squared as one minus the ratio of prediction MSE to
@@ -721,22 +759,18 @@ measure_r2 <- function(
   
   mse_y_i <- (y - mean(y))^2
   mse_y_hat <- mean(mse_y_i)
-   
-  var_mse_hat <- mse_res$estimate[2]^2     
-  cov_mse_msey <- stats::cov(sqe_i, mse_y_i) / n_obs              
-  var_mse_y_hat <- var(mse_y_i) / n_obs 
-  
-  t1 <- var_mse_hat
-  t2 <- -2 * (mse_hat / mse_y_hat) * cov_mse_msey
-  t3 <- (mse_hat^2 / mse_y_hat^2) * var_mse_y_hat
-  se_r2 <- sqrt(t1 + t2 + t3) * (1 / mse_y_hat)
-  
+
   est_r2 <- 1 - mse_hat / mse_y_hat
-  
+  se_r2 <- .se_r2_delta(sqe_i, mse_y_i)
+
   res <- list(
     estimate = est_r2,
     se = se_r2,
-    pointwise = sqe_i
+    pointwise = sqe_i,
+    # `loo_compare()` needs the baseline to propagate uncertainty into the
+    # standard error of an r2 difference; `y` is gone by then. See
+    # `.se_diff_r2()`.
+    extra = list(mse_y_i = mse_y_i)
   )
   .create_measure_structure(
     res, higher_is_better, "r2", n_draws = n_draws, n_obs = n_obs
@@ -940,20 +974,252 @@ measure_srps <- function(y, ypred, log_weights = NULL, pointwise = NULL,
 # @param measure The measure used.
 # @return The measure specification.
 .measure_spec <- list(
-  elpd = list(fun = measure_elpd, loss = FALSE),
-  ic = list(fun = measure_ic, loss = TRUE),
-  mlpd = list(fun = measure_mlpd, loss = FALSE),
-  mae = list(fun = measure_mae, loss = TRUE),
-  r2 = list(fun = measure_r2, loss = FALSE),
-  rmse = list(fun = measure_rmse, loss = TRUE),
-  mse = list(fun = measure_mse, loss = TRUE),
-  acc = list(fun = measure_acc, loss = FALSE),
-  bacc = list(fun = measure_bacc, loss = FALSE),
-  rps = list(fun = measure_rps, loss = FALSE),
-  srps = list(fun = measure_srps, loss = TRUE),
-  brier = list(fun = measure_brier, loss = TRUE)
+  elpd = list(fun = measure_elpd, loss = FALSE, diff_method = "sum"),
+  ic = list(fun = measure_ic, loss = TRUE, diff_method = "sum"),
+  mlpd = list(fun = measure_mlpd, loss = FALSE, diff_method = "sum"),
+  mae = list(fun = measure_mae, loss = TRUE, diff_method = "mean"),
+  r2 = list(
+    fun = measure_r2,
+    loss = FALSE,
+    diff_method = "pairwise",
+    se_diff_fun = "r2"
+  ),
+  rmse = list(
+    fun = measure_rmse,
+    loss = TRUE,
+    diff_method = "pairwise",
+    se_diff_fun = "rmse"
+  ),
+  mse = list(fun = measure_mse, loss = TRUE, diff_method = "mean"),
+  acc = list(fun = measure_acc, loss = FALSE, diff_method = "mean"),
+  bacc = list(
+    fun = measure_bacc,
+    loss = FALSE,
+    diff_method = "pairwise",
+    se_diff_fun = "bacc"
+  ),
+  rps = list(fun = measure_rps, loss = FALSE, diff_method = "mean"),
+  srps = list(fun = measure_srps, loss = TRUE, diff_method = "mean"),
+  brier = list(fun = measure_brier, loss = TRUE, diff_method = "mean")
 )
 
+# pairwise standard errors -----------------------------
+#
+# Measures whose overall estimate is not a sum or mean of pointwise
+# contributions cannot use the paired pointwise standard error. They register a
+# `se_diff_fun` in `.measure_spec`, naming an entry of `.se_diff_funs` below.
+# Custom measures may instead attach a function directly via
+# `attr(my_fun, "se_diff_fun")`.
+#
+# Such a function receives `ref` and `cmp`, each a list with the elements
+# `estimate`, `se`, `pointwise`, and `extra` for one model, always on the
+# measure's natural scale (`higher_is_better` sign flips are undone before the
+# call and reapplied to the difference afterwards), and returns the standard
+# error of the difference as a numeric scalar. The difference itself is always
+# `estimate_cmp - estimate_ref` and is computed by `loo_compare()`.
+
+#' Standard error of an RMSE difference
+#'
+#' First-order bivariate Taylor (delta method) approximation of the standard
+#' error of \eqn{RMSE(M_cmp) - RMSE(M_ref)}, propagated from the MSE scale on
+#' which the pointwise squared errors live.
+#'
+#' @noRd
+#' @param ref,cmp Per-model inputs; `pointwise` holds squared errors.
+#' @return Numeric scalar standard error.
+.se_diff_rmse <- function(ref, cmp) {
+  sqe_ref <- ref$pointwise
+  sqe_cmp <- cmp$pointwise
+  n <- length(sqe_ref)
+  mse_ref <- mean(sqe_ref)
+  mse_cmp <- mean(sqe_cmp)
+
+  # a perfect predictor leaves the ratios below undefined; `measure_rmse()`
+  # reports a zero standard error in that case, so do the same here
+  if (n <= 1L || mse_ref <= 0 || mse_cmp <= 0) {
+    return(0)
+  }
+
+  se_mse_ref <- sqrt(stats::var(sqe_ref) / n)
+  se_mse_cmp <- sqrt(stats::var(sqe_cmp) / n)
+  cov_mse <- sum((sqe_cmp - mse_cmp) * (sqe_ref - mse_ref)) / (n * (n - 1))
+
+  # standard errors of the two MSEs relative to their own scale, and their
+  # correlation
+  rel_ref <- se_mse_ref / sqrt(mse_ref)
+  rel_cmp <- se_mse_cmp / sqrt(mse_cmp)
+  rho <- if (se_mse_ref <= 0 || se_mse_cmp <= 0) {
+    0
+  } else {
+    cov_mse / (se_mse_cmp * se_mse_ref)
+  }
+
+  # algebraically `rel_cmp^2 + rel_ref^2 - 2 * rho * rel_cmp * rel_ref`, but
+  # written so that two models with the same squared errors cancel exactly
+  # rather than leaving rounding noise behind. Both terms are non-negative,
+  # because the correlation cannot exceed one.
+  v <- (rel_cmp - rel_ref)^2 + 2 * rel_cmp * rel_ref * max(1 - rho, 0)
+
+  0.5 * sqrt(v)
+}
+
+#' Standard error of an R-squared difference
+#'
+#' First-order trivariate Taylor (delta method) approximation of the standard
+#' error of \eqn{R^2(M_cmp) - R^2(M_ref)}. The difference equals
+#' \eqn{-MSE(M_cmp, M_ref) / MSE(y)}, so it is the same expansion as the
+#' single-model standard error in `measure_r2()` with the pointwise squared
+#' errors replaced by their pointwise differences; both go through
+#' `.se_r2_delta()`.
+#'
+#' @noRd
+#' @param ref,cmp Per-model inputs; `pointwise` holds squared errors and
+#'   `extra$mse_y_i` the baseline `(y_i - mean(y))^2` stored by `measure_r2()`.
+#' @return Numeric scalar standard error, or `NA_real_` when the baseline is
+#'   unavailable.
+.se_diff_r2 <- function(ref, cmp) {
+  # the baseline is a property of `y`, so either model's copy will do; models
+  # fitted to different `y` are already reported by the `yhash` warning
+  mse_y_i <- ref$extra$mse_y_i
+  if (is.null(mse_y_i)) {
+    mse_y_i <- cmp$extra$mse_y_i
+  }
+
+  # objects computed before the baseline was stored cannot support the
+  # covariance terms; report the difference without a standard error rather
+  # than refusing the whole comparison
+  if (!is.numeric(mse_y_i) || length(mse_y_i) != length(ref$pointwise)) {
+    return(NA_real_)
+  }
+
+  .se_r2_delta(cmp$pointwise - ref$pointwise, mse_y_i)
+}
+
+#' Standard error of a balanced-accuracy difference
+#'
+#' Balanced accuracy averages class-wise accuracies, so a difference of two
+#' balanced accuracies is a difference of two stratified means, not a mean of
+#' pointwise differences. The two sources of dependence separate: the class
+#' strata are disjoint sets of observations and so contribute independent
+#' variances, while within a stratum both models score the *same* `n_c`
+#' observations and are therefore paired. Writing
+#' \eqn{d_i = acc_i(M_cmp) - acc_i(M_ref)},
+#'
+#' \deqn{SE = \frac{1}{K} \sqrt{\sum_c Var(d_i : i \in c) / n_c}}
+#'
+#' which is the difference-analogue of the single-model
+#' \eqn{\sqrt{\sum_c acc_c (1 - acc_c) / n_c} / K} in `measure_bacc()`: the
+#' per-stratum binomial variance replaced by the paired-difference variance.
+#' Within a stratum this is the McNemar variance of a paired difference of
+#' proportions, \eqn{(b + c)/n_c^2 - (b - c)^2/n_c^3} in discordant-pair form,
+#' up to the \eqn{n_c/(n_c - 1)} of the sample variance.
+#'
+#' For binary outcomes balanced accuracy is \eqn{(sens + spec)/2}, so this is
+#' the estimand of Newcombe (2001) at a mixing parameter of one half, and half
+#' the difference of two Youden indices in a paired design (Chen et al., 2015).
+#'
+#' @references
+#' Newcombe, R. G. (2001). Simultaneous comparison of sensitivity and
+#' specificity of two tests in the paired design: a straightforward graphical
+#' approach. *Statistics in Medicine*, 20(6):907--915.
+#'
+#' Chen, F., Xue, Y., Tan, M. T., and Chen, P. (2015). Efficient statistical
+#' tests to compare Youden index: accounting for contingency correlation.
+#' *Statistics in Medicine*, 34(9):1560--1576.
+#'
+#' @noRd
+#' @param ref,cmp Per-model inputs; `pointwise` holds `acc_i / (K * n_c)` and
+#'   `extra$class_id` the class index stored by `measure_bacc()`.
+#' @return Numeric scalar standard error, or `NA_real_` when the class strata
+#'   are unavailable.
+.se_diff_bacc <- function(ref, cmp) {
+  # the strata are a property of `y`, so either model's copy will do; models
+  # fitted to different `y` are already reported by the `yhash` warning
+  class_id <- ref$extra$class_id
+  if (is.null(class_id)) {
+    class_id <- cmp$extra$class_id
+  }
+
+  n <- length(ref$pointwise)
+  # objects computed before the strata were stored cannot be stratified;
+  # report the difference without a standard error rather than refusing the
+  # whole comparison
+  if (!is.numeric(class_id) || length(class_id) != n) {
+    return(NA_real_)
+  }
+
+  n_c <- tabulate(class_id)
+  K <- length(n_c)
+  # undo the `acc_i / (K * n_c)` scaling to recover the 0/1 accuracies, so the
+  # variances below are on the natural per-observation scale
+  d <- (cmp$pointwise - ref$pointwise) * (K * n_c[class_id])
+
+  # a stratum holding a single observation supports no variance estimate and
+  # contributes nothing, matching `measure_bacc()`, where `acc_c` is then 0 or
+  # 1 and its binomial variance vanishes
+  var_c <- vapply(
+    seq_len(K),
+    function(k) {
+      if (n_c[k] < 2L) 0 else stats::var(d[class_id == k]) / n_c[k]
+    },
+    numeric(1)
+  )
+
+  sqrt(sum(var_c)) / K
+}
+
+# registry of built-in `se_diff_fun` implementations, referenced by name from
+# `.measure_spec` so that stored objects carry a string rather than a closure
+.se_diff_funs <- list(
+  rmse = .se_diff_rmse,
+  r2 = .se_diff_r2,
+  bacc = .se_diff_bacc
+)
+
+#' Return comparison metadata for a measure
+#' @noRd
+#' @param measure_entry Normalized measure entry, or a built-in measure name.
+#' @param higher_is_better Value of `higher_is_better` used for this measure.
+.measure_compare_meta <- function(measure_entry, higher_is_better = NULL) {
+  if (is.character(measure_entry)) {
+    measure_entry <- list(
+      name = measure_entry,
+      type = "builtin",
+      key = measure_entry
+    )
+  }
+
+  if (measure_entry$type == "builtin") {
+    entry <- .measure_spec[[measure_entry$key]]
+    if (is.null(entry)) {
+      return(NULL)
+    }
+    return(list(
+      higher_is_better = higher_is_better,
+      loss = isTRUE(entry$loss),
+      diff_method = entry$diff_method,
+      se_diff_fun = entry$se_diff_fun
+    ))
+  }
+
+  se_diff_fun <- attr(measure_entry$key, "se_diff_fun", exact = TRUE)
+  if (!is.null(se_diff_fun) && !is.function(se_diff_fun)) {
+    cli::cli_abort(c(
+      "Attribute {.code se_diff_fun} of custom measure",
+      "{.val {measure_entry$name}} must be a function.",
+      "i" = "It is called as {.code se_diff_fun(ref, cmp)} and must return a",
+      " " = "numeric scalar."
+    ))
+  }
+
+  list(
+    higher_is_better = higher_is_better,
+    loss = FALSE,
+    diff_method = if (is.null(se_diff_fun)) "auto" else "pairwise",
+    se_diff_fun = se_diff_fun
+  )
+}
+
 #' Supported predictive measure names
 #'
 #' A character vector of measure names that can be passed to the `measure`
@@ -991,6 +1257,26 @@ supported_measures_list <- names(.measure_spec)
     out,
     class = c("measure", "loo"),
     measure = measure_name,
-    dims = c(n_draws, n_obs)
+    dims = c(n_draws, n_obs),
+    higher_is_better = higher_is_better,
+    # not affected by the sign flip above: `extra` carries auxiliary data for
+    # `se_diff_fun()`, which always works on the measure's natural scale
+    compare_extra = res$extra
   )
 }
+
+#' Auxiliary data a measure stores for its `se_diff_fun`
+#'
+#' Built-in measures carry it as the `compare_extra` attribute added by
+#' `.create_measure_structure()`; custom measures return it as an `extra`
+#' element of their result list.
+#' @noRd
+#' @param res A measure result.
+#' @return A list, or `NULL` when the measure stores nothing.
+.measure_compare_extra <- function(res) {
+  extra <- attr(res, "compare_extra", exact = TRUE)
+  if (is.null(extra)) {
+    extra <- res$extra
+  }
+  extra
+}
diff --git a/R/pred_measure-compute.R b/R/pred_measure-compute.R
index 9fbbcd54..0b555f33 100644
--- a/R/pred_measure-compute.R
+++ b/R/pred_measure-compute.R
@@ -34,6 +34,24 @@
 #'   Custom functions are called with any of `y`, `ypred`, `mupred`, `ylp`, and
 #'   `log_weights` that appear in their formals, plus arguments from `control`.
 #'   They must return a list with  `estimates` and `pointwise`.
+#'
+#'   For [loo_compare()], the standard error of a difference is derived from
+#'   paired pointwise differences when the overall estimate is the sum or the
+#'   mean of `pointwise`; otherwise it is `NA`. A custom measure that is a
+#'   transformation of pointwise quantities (as `rmse` is of squared errors) can
+#'   supply its own by setting
+#'   `attr(my_fun, "se_diff_fun") <- function(ref, cmp) ...`. That function
+#'   receives one list per model with elements `estimate`, `se`, `pointwise`,
+#'   and `extra`, always on the measure's natural scale, and must return the
+#'   standard error of the difference as a numeric scalar. Note that the
+#'   function is stored on the result object, so it carries its enclosing
+#'   environment into any saved copy of that object.
+#'
+#'   `extra` is for anything the standard error needs that the pointwise values
+#'   do not carry. Return it as an additional list element `extra` from the
+#'   measure function and it is stored alongside the estimates; the built-in
+#'   `r2` uses it for the baseline `(y_i - mean(y))^2`, which cannot be
+#'   recovered once `y` is out of scope.
 #' @param measure_name For a single custom function, set
 #'   `attr(my_fun, "measure_name") <- "my_metric"` before passing `my_fun` to
 #'   `measure`.
@@ -165,7 +183,10 @@ do_pred_measure <- function(
       mat = estimates,
       name = entry$name,
       values = .measure_estimate_se(sel_measure),
-      margin = 1
+      margin = 1,
+      measure_entry = entry,
+      higher_is_better = attr(sel_measure, "higher_is_better"),
+      extra = .measure_compare_extra(sel_measure)
     )
     pointwise <- .merge_matrix(
       source = source,
@@ -184,7 +205,7 @@ do_pred_measure <- function(
     save_psis = save_psis
   )
   
-  .add_attributes(
+  predperf_res <- .add_attributes(
     save_psis,
     predperf_res, 
     y, 
@@ -197,6 +218,7 @@ do_pred_measure <- function(
     predperf, 
     source
   )
+  predperf_res
   }
 
 # internal helper functions ---------------------------------------------------
@@ -323,10 +345,11 @@ do_pred_measure <- function(
     base_measure
 ) {
   if (measure_entry$type == "builtin") {
-    measure_fun <- .measure_spec[[measure_entry$key]]$fun
-    if (is.null(measure_fun)) {
+    spec <- .measure_spec[[measure_entry$key]]
+    if (is.null(spec)) {
       cli::cli_abort("Unknown built-in measure {.val {measure_entry$key}}.")
     }
+    measure_fun <- spec$fun
   } else {
     measure_fun <- measure_entry$key
   }
@@ -536,6 +559,14 @@ do_pred_measure <- function(
 #'   `(estimate, se)`; for `margin = 2`, length-`n` pointwise vector.
 #' @param margin `1` to merge along rows (estimates table), `2` along columns
 #'   (pointwise table).
+#' @param measure_entry Optional normalized measure entry; when merging an
+#'   estimates row (`margin = 1`), comparison metadata for [loo_compare()] is
+#'   recorded from this entry and `higher_is_better`.
+#' @param higher_is_better Optional logical or `NULL`; records the orientation
+#'   used for `name` when merging an estimates row (`margin = 1`).
+#' @param extra Optional list of auxiliary data the measure stores for its
+#'   `se_diff_fun` (see `.measure_compare_extra()`); recorded in the comparison
+#'   metadata when merging an estimates row (`margin = 1`).
 #'
 #' @return Updated matrix with `name` as a row or column name.
 #'
@@ -551,7 +582,16 @@ do_pred_measure <- function(
 }
 
 #' @noRd
-.merge_matrix <- function(source, mat, name, values, margin) {
+.merge_matrix <- function(
+  source,
+  mat,
+  name,
+  values,
+  margin,
+  measure_entry = NULL,
+  higher_is_better = NULL,
+  extra = NULL
+) {
   is_row <- margin == 1
   bind_fn <- if (is_row) rbind else cbind
   name_updated <- .measure_result_name(source, name)
@@ -562,8 +602,46 @@ do_pred_measure <- function(
     matrix(values, ncol = 1, dimnames = list(NULL, name_updated))
   }
 
-  if (is.null(mat)) return(new_slice)
-  bind_fn(mat, new_slice)
+  compare_meta <- if (is_row && !is.null(measure_entry)) {
+    .measure_compare_meta(measure_entry, higher_is_better)
+  }
+  if (!is.null(compare_meta) && !is.null(extra)) {
+    compare_meta$extra <- extra
+  }
+
+  old_higher_is_better <- if (is_row && !is.null(mat)) {
+    attr(mat, "measure_higher_is_better")
+  }
+  old_compare_meta <- if (is_row && !is.null(mat)) {
+    attr(mat, "measure_compare_meta")
+  }
+
+  mat <- if (is.null(mat)) new_slice else bind_fn(mat, new_slice)
+
+  if (is_row) {
+    if (!is.null(measure_entry) || !is.null(old_higher_is_better)) {
+      measure_higher_is_better <- old_higher_is_better
+      if (is.null(measure_higher_is_better)) {
+        measure_higher_is_better <- list()
+      }
+      if (!is.null(measure_entry)) {
+        measure_higher_is_better[[name]] <- higher_is_better
+      }
+      attr(mat, "measure_higher_is_better") <- measure_higher_is_better
+    }
+    if (!is.null(compare_meta) || !is.null(old_compare_meta)) {
+      measure_compare_meta <- old_compare_meta
+      if (is.null(measure_compare_meta)) {
+        measure_compare_meta <- list()
+      }
+      if (!is.null(compare_meta)) {
+        measure_compare_meta[[name]] <- compare_meta
+      }
+      attr(mat, "measure_compare_meta") <- measure_compare_meta
+    }
+  }
+
+  mat
 }
 
 #' Construct the S3 predictive measure result object
@@ -586,7 +664,9 @@ do_pred_measure <- function(
 #' @param save_psis Logical; if `TRUE`, include `psis_object` in the result.
 #'
 #' @return A list with elements `estimates`, `pointwise`, and optionally
-#'   `diagnostics`, `psis_object`, and `log_weights`. Class attributes are added
+#'   `diagnostics`, `psis_object`, and `log_weights`. Attributes
+#'   `measure_higher_is_better` and `measure_compare_meta` record per-measure
+#'   metadata for measures added in the current call. Class attributes are added
 #'   by \code{.add_attributes()}.
 #'
 #' @noRd
@@ -597,6 +677,18 @@ do_pred_measure <- function(
   psis_object,
   save_psis
 ) {
+  measure_higher_is_better <- attr(estimates, "measure_higher_is_better")
+  if (is.null(measure_higher_is_better)) {
+    measure_higher_is_better <- list()
+  }
+  attr(estimates, "measure_higher_is_better") <- NULL
+
+  measure_compare_meta <- attr(estimates, "measure_compare_meta")
+  if (is.null(measure_compare_meta)) {
+    measure_compare_meta <- list()
+  }
+  attr(estimates, "measure_compare_meta") <- NULL
+
   output_list <- list(
     estimates = estimates,
     pointwise = pointwise
@@ -610,23 +702,32 @@ do_pred_measure <- function(
   if (!is.null(psis_object)) {
     output_list$log_weights <- psis_object$log_weights
   }
-  
-  structure(output_list)
+
+  structure(
+    output_list,
+    measure_higher_is_better = measure_higher_is_better,
+    measure_compare_meta = measure_compare_meta
+  )
 }
 
 #' Attach S3 classes and metadata attributes to a result
 #'
 #' @description
-#' Sets `class`, `source`, and `dims` attributes on a predictive measure object.
+#' Sets `class`, `source`, `dims`, `measure_higher_is_better`, and
+#' `measure_compare_meta` attributes on a predictive measure object.
 #'
 #' When updating an existing result (`predperf` is not `NULL`), copies attributes
 #' from `predperf` and refreshes `dims` from newly supplied input matrices.
+#' Merges `measure_higher_is_better` and `measure_compare_meta` from the prior result
+#' with any new entries supplied on `predperf_res` (from
+#' \code{.build_pred_measure()}).
 #' When `save_psis = FALSE`, clears any stored `psis_object` from the prior
 #' result.
 #'
 #' For new objects, copies relevant attributes from `loo` or `kfold` inputs
 #' (e.g. `yhash`, `model_name`, fold structure) and assigns a source-specific
-#' subclass (`"insample_pred_measure"`, `"loo_pred_measure"`, etc.).
+#' subclass (`"insample_pred_measure"`, `"loo_pred_measure"`, etc.). Sets
+#' `measure_higher_is_better` and `measure_compare_meta`, seeding `elpd` defaults.
 #'
 #' @param save_psis Logical; when `FALSE` and accumulating, clears stored
 #'   `psis_object` from the prior result.
@@ -646,7 +747,28 @@ do_pred_measure <- function(
 #' @return The updated `predperf_res` with class and attributes set.
 #'
 #' @noRd
-.add_attributes <- function(save_psis, predperf_res, y, ypred, mupred, ylp, ylp_test, kfold, loo, predperf, source) {
+.add_attributes <- function(
+  save_psis,
+  predperf_res,
+  y,
+  ypred,
+  mupred,
+  ylp,
+  ylp_test,
+  kfold,
+  loo,
+  predperf,
+  source
+) {
+  new_higher_is_better <- attr(predperf_res, "measure_higher_is_better")
+  if (is.null(new_higher_is_better)) {
+    new_higher_is_better <- list()
+  }
+  new_compare_meta <- attr(predperf_res, "measure_compare_meta")
+  if (is.null(new_compare_meta)) {
+    new_compare_meta <- list()
+  }
+
   if (!is.null(predperf)) {
     if (isFALSE(save_psis)) {
       predperf$psis_object <- NULL
@@ -661,6 +783,25 @@ do_pred_measure <- function(
       dim(ylp)
     }
     attr(predperf_res, "dims") <- dims
+    measure_higher_is_better <- attr(predperf, "measure_higher_is_better")
+    if (is.null(measure_higher_is_better)) {
+      measure_higher_is_better <- list()
+    }
+    if (length(new_higher_is_better)) {
+      measure_higher_is_better[names(new_higher_is_better)] <- new_higher_is_better
+    }
+    attr(predperf_res, "measure_higher_is_better") <- measure_higher_is_better
+    compare_meta <- attr(predperf, "measure_compare_meta")
+    if (is.null(compare_meta)) {
+      compare_meta <- list()
+    }
+    if (is.null(compare_meta$elpd)) {
+      compare_meta$elpd <- .measure_compare_meta("elpd")
+    }
+    if (length(new_compare_meta)) {
+      compare_meta[names(new_compare_meta)] <- new_compare_meta
+    }
+    attr(predperf_res, "measure_compare_meta") <- compare_meta
     
     return(predperf_res)
   }
@@ -708,6 +849,19 @@ do_pred_measure <- function(
   }
   attr(predperf_res, "class") <- classes
   attr(predperf_res, "source") <- source
+  measure_higher_is_better <- list(elpd = NULL)
+  if (length(new_higher_is_better)) {
+    measure_higher_is_better[names(new_higher_is_better)] <- new_higher_is_better
+  }
+  attr(predperf_res, "measure_higher_is_better") <- measure_higher_is_better
+  compare_meta <- list(
+    elpd = .measure_compare_meta("elpd")
+  )
+  if (length(new_compare_meta)) {
+    compare_meta[names(new_compare_meta)] <- new_compare_meta
+  }
+  attr(predperf_res, "measure_compare_meta") <- compare_meta
   
   return(predperf_res)
-}
\ No newline at end of file
+}
+
diff --git a/R/pred_measure-helpers.R b/R/pred_measure-helpers.R
index 761f2a6e..ce6e47d1 100644
--- a/R/pred_measure-helpers.R
+++ b/R/pred_measure-helpers.R
@@ -236,6 +236,13 @@
       "length {.val {n_obs}}, not {.val {length(res$pointwise)}}."
     ))
   }
+  if (!is.null(res$extra) && !is.list(res$extra)) {
+    cli::cli_abort(c(
+      "{.field extra} from custom measure {.val {measure_name}} must be a list.",
+      "i" = "It is handed to {.code se_diff_fun(ref, cmp)} as
+             {.code ref$extra} and {.code cmp$extra}."
+    ))
+  }
 
   invisible(res)
 }
diff --git a/R/pred_measure.R b/R/pred_measure.R
index 829dfb59..2a5da6d2 100644
--- a/R/pred_measure.R
+++ b/R/pred_measure.R
@@ -28,7 +28,11 @@
 #' }
 #'
 #' The attribute `source` is `"insample"`. Attribute `dims` gives posterior
-#' draws × observations. Use [print()] for a readable summary table.
+#' draws × observations. Attribute `measure_higher_is_better` records the
+#' `higher_is_better` setting used for each measure; see section below. Use [print()]
+#' for a readable summary table.
+#'
+#' @template measure-higher-is-better-attribute
 #'
 #' @details
 #' **Input requirements by measure.** Supply only the inputs each measure
@@ -52,6 +56,11 @@
 #' declared in the function signature among `y`, `ypred`, `mupred`, `ylp`, and
 #' `log_weights` are supplied automatically.
 #'
+#' Custom measures are assumed to be on a utility scale (higher is better) in
+#' [loo_compare()]. For a custom loss measure, pass
+#' `control = list(my_measure = list(higher_is_better = TRUE))` or negate values
+#' in the custom function so that [loo_compare()] ranks models correctly.
+#'
 #' @examples
 #' \donttest{
 #' if (requireNamespace("brms", quietly = TRUE)) {
@@ -151,6 +160,8 @@ insample_pred_measure <- function(
 #'
 #' Measure names carry a `_loo` suffix (e.g. `elpd_loo`, `crps_loo`).
 #'
+#' @template measure-higher-is-better-attribute
+#'
 #' @details
 #' **Three equivalent input patterns:**
 #'
@@ -241,6 +252,8 @@ loo_pred_measure <- function(
 #' list contains `estimates` and `pointwise`; measure names carry a `_kfold`
 #' suffix (e.g. `elpd_kfold`, `crps_kfold`).
 #'
+#' @template measure-higher-is-better-attribute
+#'
 #' @details
 #' For distributional measures on held-out folds, obtain posterior predictions
 #' with `brms::kfold_predict()` and pass the resulting `yrep` matrices as
@@ -323,6 +336,8 @@ kfold_pred_measure <- function(
 #' `elpd_test`, `crps_test`). Attribute `dims` reflects the test-set size
 #' (from `ylp_test`), not the training data.
 #'
+#' @template measure-higher-is-better-attribute
+#'
 #' @details
 #' The base summary `elpd_test` is computed from `ylp_test` on the holdout
 #' observations only.
@@ -401,7 +416,10 @@ test_pred_measure <- function(
 #' An updated object of the same class as `predperf`, with new rows in
 #' `estimates` and columns in `pointwise` for each requested measure. Base
 #' summaries (`elpd` and LOO/k-fold complexity terms such as `p_loo`) are not
-#' recomputed.
+#' recomputed. Attribute `measure_higher_is_better` is updated for any newly added
+#' measures.
+#'
+#' @template measure-higher-is-better-attribute
 #'
 #' @details
 #' **Typical workflow:**
diff --git a/_pkgdown.yml b/_pkgdown.yml
index b3d56da9..f1adbeae 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -66,6 +66,7 @@ articles:
     contents:
       - articles-online-only/pred-measure-workflow
       - articles-online-only/overview-measures
+      - articles-online-only/comparison-standard-errors
   - title: Frequently asked questions
     contents:
       - faq
diff --git a/man-roxygen/measure-higher-is-better-attribute.R b/man-roxygen/measure-higher-is-better-attribute.R
new file mode 100644
index 00000000..5500c9e7
--- /dev/null
+++ b/man-roxygen/measure-higher-is-better-attribute.R
@@ -0,0 +1,16 @@
+#' @section `measure_higher_is_better` attribute:
+#' All `*_pred_measure()` and [pred_measure()] results include attribute
+#' `measure_higher_is_better`: a named list recording the `higher_is_better`
+#' setting used for each measure (`NULL`, `TRUE`, or `FALSE` per bare measure
+#' name; `elpd` is always `NULL`). Built-in loss measures such as MSE are stored
+#' on a loss scale by default (`NULL`); pass
+#' `control = list(mse = list(higher_is_better = TRUE))` to store values on a
+#' utility scale. When measures are added incrementally with [pred_measure()],
+#' the attribute is updated for newly computed measures.
+#'
+#' Attribute `measure_compare_meta` records per-measure comparison metadata
+#' (`higher_is_better`, `loss`, and `diff_method`) used by [loo_compare()].
+#' Built-in measures take `loss` and `diff_method` from the package measure
+#' registry; custom measures default to `diff_method = "auto"` and `loss = FALSE`
+#' (utility scale). For custom loss measures, use `higher_is_better = TRUE` in
+#' `control` so [loo_compare()] interprets the sign correctly. See [loo-glossary].
diff --git a/man/insample_pred_measure.Rd b/man/insample_pred_measure.Rd
index 98b3598e..48b93c1e 100644
--- a/man/insample_pred_measure.Rd
+++ b/man/insample_pred_measure.Rd
@@ -42,7 +42,25 @@ functions (custom measures), e.g. \code{list("rps", my_metric = my_fun)}.
 }
 Custom functions are called with any of \code{y}, \code{ypred}, \code{mupred}, \code{ylp}, and
 \code{log_weights} that appear in their formals, plus arguments from \code{control}.
-They must return a list with  \code{estimates} and \code{pointwise}.}
+They must return a list with  \code{estimates} and \code{pointwise}.
+
+For \code{\link[=loo_compare]{loo_compare()}}, the standard error of a difference is derived from
+paired pointwise differences when the overall estimate is the sum or the
+mean of \code{pointwise}; otherwise it is \code{NA}. A custom measure that is a
+transformation of pointwise quantities (as \code{rmse} is of squared errors) can
+supply its own by setting
+\code{attr(my_fun, "se_diff_fun") <- function(ref, cmp) ...}. That function
+receives one list per model with elements \code{estimate}, \code{se}, \code{pointwise},
+and \code{extra}, always on the measure's natural scale, and must return the
+standard error of the difference as a numeric scalar. Note that the
+function is stored on the result object, so it carries its enclosing
+environment into any saved copy of that object.
+
+\code{extra} is for anything the standard error needs that the pointwise values
+do not carry. Return it as an additional list element \code{extra} from the
+measure function and it is stored alongside the estimates; the built-in
+\code{r2} uses it for the baseline \code{(y_i - mean(y))^2}, which cannot be
+recovered once \code{y} is out of scope.}
 
 \item{group_ids}{Optional vector of group identifiers for grouped summaries
 (reserved; not yet implemented).}
@@ -67,7 +85,9 @@ per measure).}
 }
 
 The attribute \code{source} is \code{"insample"}. Attribute \code{dims} gives posterior
-draws × observations. Use \code{\link[=print]{print()}} for a readable summary table.
+draws × observations. Attribute \code{measure_higher_is_better} records the
+\code{higher_is_better} setting used for each measure; see section below. Use \code{\link[=print]{print()}}
+for a readable summary table.
 }
 \description{
 Compute predictive performance measures on the same data used to fit the
@@ -104,7 +124,31 @@ article for definitions and orientation (higher vs lower is better).
 \code{measure_name} and return \code{estimate}, \code{se}, and \code{pointwise}. Only arguments
 declared in the function signature among \code{y}, \code{ypred}, \code{mupred}, \code{ylp}, and
 \code{log_weights} are supplied automatically.
+
+Custom measures are assumed to be on a utility scale (higher is better) in
+\code{\link[=loo_compare]{loo_compare()}}. For a custom loss measure, pass
+\code{control = list(my_measure = list(higher_is_better = TRUE))} or negate values
+in the custom function so that \code{\link[=loo_compare]{loo_compare()}} ranks models correctly.
 }
+\section{\code{measure_higher_is_better} attribute}{
+
+All \verb{*_pred_measure()} and \code{\link[=pred_measure]{pred_measure()}} results include attribute
+\code{measure_higher_is_better}: a named list recording the \code{higher_is_better}
+setting used for each measure (\code{NULL}, \code{TRUE}, or \code{FALSE} per bare measure
+name; \code{elpd} is always \code{NULL}). Built-in loss measures such as MSE are stored
+on a loss scale by default (\code{NULL}); pass
+\code{control = list(mse = list(higher_is_better = TRUE))} to store values on a
+utility scale. When measures are added incrementally with \code{\link[=pred_measure]{pred_measure()}},
+the attribute is updated for newly computed measures.
+
+Attribute \code{measure_compare_meta} records per-measure comparison metadata
+(\code{higher_is_better}, \code{loss}, and \code{diff_method}) used by \code{\link[=loo_compare]{loo_compare()}}.
+Built-in measures take \code{loss} and \code{diff_method} from the package measure
+registry; custom measures default to \code{diff_method = "auto"} and \code{loss = FALSE}
+(utility scale). For custom loss measures, use \code{higher_is_better = TRUE} in
+\code{control} so \code{\link[=loo_compare]{loo_compare()}} interprets the sign correctly. See \link{loo-glossary}.
+}
+
 \examples{
 \donttest{
 if (requireNamespace("brms", quietly = TRUE)) {
diff --git a/man/kfold_pred_measure.Rd b/man/kfold_pred_measure.Rd
index aa0588e4..c5018ec4 100644
--- a/man/kfold_pred_measure.Rd
+++ b/man/kfold_pred_measure.Rd
@@ -42,7 +42,25 @@ functions (custom measures), e.g. \code{list("rps", my_metric = my_fun)}.
 }
 Custom functions are called with any of \code{y}, \code{ypred}, \code{mupred}, \code{ylp}, and
 \code{log_weights} that appear in their formals, plus arguments from \code{control}.
-They must return a list with  \code{estimates} and \code{pointwise}.}
+They must return a list with  \code{estimates} and \code{pointwise}.
+
+For \code{\link[=loo_compare]{loo_compare()}}, the standard error of a difference is derived from
+paired pointwise differences when the overall estimate is the sum or the
+mean of \code{pointwise}; otherwise it is \code{NA}. A custom measure that is a
+transformation of pointwise quantities (as \code{rmse} is of squared errors) can
+supply its own by setting
+\code{attr(my_fun, "se_diff_fun") <- function(ref, cmp) ...}. That function
+receives one list per model with elements \code{estimate}, \code{se}, \code{pointwise},
+and \code{extra}, always on the measure's natural scale, and must return the
+standard error of the difference as a numeric scalar. Note that the
+function is stored on the result object, so it carries its enclosing
+environment into any saved copy of that object.
+
+\code{extra} is for anything the standard error needs that the pointwise values
+do not carry. Return it as an additional list element \code{extra} from the
+measure function and it is stored alongside the estimates; the built-in
+\code{r2} uses it for the baseline \code{(y_i - mean(y))^2}, which cannot be
+recovered once \code{y} is out of scope.}
 
 \item{kfold}{A \code{kfold} object from \code{\link[brms:kfold]{brms::kfold()}}. Supplies ELPD summaries
 and fold structure for \code{\link[=kfold_pred_measure]{kfold_pred_measure()}}; pass \code{y}, \code{ypred}, and/or
@@ -78,6 +96,25 @@ with \code{brms::kfold_predict()} and pass the resulting \code{yrep} matrices as
 \code{ypred} and/or \code{mupred}. See the sleep-study workflow in
 \href{https://mc-stan.org/loo/articles/articles-online-only/pred-measure-workflow.html}{pred-measure workflow article}.
 }
+\section{\code{measure_higher_is_better} attribute}{
+
+All \verb{*_pred_measure()} and \code{\link[=pred_measure]{pred_measure()}} results include attribute
+\code{measure_higher_is_better}: a named list recording the \code{higher_is_better}
+setting used for each measure (\code{NULL}, \code{TRUE}, or \code{FALSE} per bare measure
+name; \code{elpd} is always \code{NULL}). Built-in loss measures such as MSE are stored
+on a loss scale by default (\code{NULL}); pass
+\code{control = list(mse = list(higher_is_better = TRUE))} to store values on a
+utility scale. When measures are added incrementally with \code{\link[=pred_measure]{pred_measure()}},
+the attribute is updated for newly computed measures.
+
+Attribute \code{measure_compare_meta} records per-measure comparison metadata
+(\code{higher_is_better}, \code{loss}, and \code{diff_method}) used by \code{\link[=loo_compare]{loo_compare()}}.
+Built-in measures take \code{loss} and \code{diff_method} from the package measure
+registry; custom measures default to \code{diff_method = "auto"} and \code{loss = FALSE}
+(utility scale). For custom loss measures, use \code{higher_is_better = TRUE} in
+\code{control} so \code{\link[=loo_compare]{loo_compare()}} interprets the sign correctly. See \link{loo-glossary}.
+}
+
 \examples{
 \donttest{
 if (requireNamespace("brms", quietly = TRUE)) {
diff --git a/man/loo-glossary.Rd b/man/loo-glossary.Rd
index 58529429..b4089796 100644
--- a/man/loo-glossary.Rd
+++ b/man/loo-glossary.Rd
@@ -270,6 +270,117 @@ section.
 }
 }
 
+\section{Multi-measure model comparisons}{
+
+
+When comparing \code{\link[=loo_pred_measure]{loo_pred_measure()}} objects with
+\code{loo_compare()}, paired differences are computed for every predictive
+measure common to all models. Models are ranked by the \code{rank_by} argument
+(default \code{"elpd"}); the top-ranked model is the reference for all difference
+columns.
+\subsection{\verb{\{measure\}_diff} and \verb{\{measure\}_se_diff}}{
+
+For each non-ELPD measure \code{m}, \code{loo_compare()} adds columns \code{m_diff} and
+\code{m_se_diff}. In all cases \code{m_diff} is the difference between the two overall
+estimates on a utility scale (higher is better; loss measures such as MSE,
+Brier score, and SRPS have their sign flipped from the raw loss orientation).
+Measures already returned on a utility scale (e.g. ELPD, CRPS/RPS) are not
+sign-flipped. Negative \code{m_diff} values then indicate worse performance than
+the reference model, which has \code{m_diff = 0}.
+
+How \code{m_se_diff} is obtained depends on the measure:
+\itemize{
+\item When the overall estimate is a sum or mean of pointwise contributions, it
+is computed from paired pointwise differences using the same approach as
+\code{elpd_diff} and \code{se_diff} (Eq 24 in VGG2017 for sums; the mean analogue for
+means). This covers ELPD, \code{mlpd}, \code{ic}, \code{mae}, \code{mse}, \code{acc}, \code{brier}, and
+the ranked probability scores.
+\item When it is a transformation of such quantities, the measure supplies its
+own delta-method standard error (\code{se_diff_fun}). For \code{rmse} this is the
+first-order bivariate Taylor approximation propagated from the MSE scale,
+which requires the covariance between the two models' pointwise squared
+errors and is therefore not a paired pointwise standard deviation. For
+\code{r2} it is the trivariate analogue, which additionally propagates the
+uncertainty in the baseline \code{MSE(y)} shared by both models.
+\item Otherwise \code{m_se_diff} is \code{NA}. This applies only to custom measures whose
+estimate is neither a sum nor a mean of their pointwise values and which do
+not attach an \code{se_diff_fun}.
+}
+
+The reference model has \code{m_se_diff = 0} whenever an \code{m_se_diff} is available.
+Attribute \code{measure_higher_is_better} on each \verb{*_pred_measure()}
+result records the \code{higher_is_better} setting used when each measure was
+computed; when stored values are on a loss scale, \code{loo_compare()} emits a
+short message naming those measures (see \code{\link[=loo_compare]{loo_compare()}}).
+
+ELPD-family measures use the column names \code{elpd_diff} and \code{se_diff} rather
+than a prefixed form. Only ELPD comparisons include \code{p_worse} and \code{diag_diff};
+these diagnostics do not apply to other predictive measures.
+}
+
+\subsection{\code{measure_higher_is_better}}{
+
+Attribute on all \verb{*_pred_measure()} and \code{\link[=pred_measure]{pred_measure()}} results: a named
+list recording the \code{higher_is_better} setting used for each measure (\code{NULL},
+\code{TRUE}, or \code{FALSE}; \code{elpd} is always \code{NULL}). Used by \code{\link[=loo_compare]{loo_compare()}} with
+\code{measure_compare_meta} to decide whether paired differences need a sign flip
+when converting to a utility scale.
+}
+
+\subsection{\code{measure_compare_meta}}{
+
+Attribute on all \verb{*_pred_measure()} and \code{\link[=pred_measure]{pred_measure()}} results: a named
+list of per-measure comparison metadata used by \code{\link[=loo_compare]{loo_compare()}}. Each entry
+is a list with:
+\itemize{
+\item \code{higher_is_better} — the orientation setting used when the measure was
+computed (\code{NULL}, \code{TRUE}, or \code{FALSE})
+\item \code{loss} — whether stored values are on a loss scale (lower is better)
+\item \code{diff_method} — how the standard error of the difference is obtained:
+\code{"sum"} or \code{"mean"} (paired pointwise differences), \code{"pairwise"} (the
+measure's own \code{se_diff_fun}), \code{"estimates_only"} (unavailable, \code{NA}), or
+\code{"auto"} (inferred at compare time for custom measures). No built-in
+measure declares \code{"estimates_only"}; it is reached only when
+autodetection under \code{"auto"} cannot establish that the estimate is a sum
+or mean of its pointwise values. It is not an error state — the
+difference itself is still reported, and only its standard error is
+marked unavailable.
+\item \code{se_diff_fun} — for \code{diff_method = "pairwise"}, either the name of a
+built-in implementation or, for custom measures, the function itself
+\item \code{extra} — optional list of auxiliary data the measure stored for its
+\code{se_diff_fun}, present only for measures that need it (\code{r2} stores the
+pointwise baseline \code{(y_i - mean(y))^2}, which \code{y} no longer supplies by
+the time \code{\link[=loo_compare]{loo_compare()}} runs; \code{bacc} stores the class index of each
+observation, which its pointwise values do not determine). It is excluded
+from the metadata consistency check below, since it varies with the data
+rather than with how the measure was configured.
+}
+
+Built-in measures take \code{loss}, \code{diff_method}, and \code{se_diff_fun} from the
+package measure registry; custom measures default to \code{loss = FALSE} and
+\code{diff_method = "auto"} unless they attach an \code{se_diff_fun} (see
+\code{\link[=insample_pred_measure]{insample_pred_measure()}}).
+\code{\link[=loo_compare]{loo_compare()}} requires all models to provide matching metadata for each
+shared measure; mismatched \code{higher_is_better} settings or missing metadata on
+some models produce an error.
+}
+
+\subsection{\code{rank_by}, \code{compare_measures}, and related attributes`}{
+
+The \code{rank_by} argument selects which measure determines model ordering and
+the reference model for all pairwise differences. When \code{rank_by} is omitted,
+models are ranked by \code{"elpd"}; attribute \code{rank_by} is set only when \code{rank_by}
+is passed explicitly. Attribute \code{compare_measures} lists all measures that
+were compared, and \code{sign_converted_measures} lists loss measures whose sign
+was flipped onto the utility scale. The print method shows the ranking
+measure by default
+(\code{"elpd"} when \code{rank_by} was not set); use \code{print(x, measures = "all")} or
+\code{print(x, measures = c("rmse", "r2"))} to display additional measure tables.
+Printed tables label the standard-error column \code{se_diff} even for non-ELPD
+measures; the data frame columns remain \verb{\{measure\}_se_diff}.
+}
+}
+
 \references{
 Vehtari, A., Gelman, A., and Gabry, J. (2017). Practical Bayesian model
 evaluation using leave-one-out cross-validation and WAIC.
diff --git a/man/loo_compare.Rd b/man/loo_compare.Rd
index c54b9f86..e07f83ec 100644
--- a/man/loo_compare.Rd
+++ b/man/loo_compare.Rd
@@ -8,21 +8,28 @@
 \alias{print.compare.loo_ss}
 \title{Model comparison}
 \usage{
-loo_compare(x, ...)
+loo_compare(x, ..., rank_by = NULL)
 
-\method{loo_compare}{default}(x, ...)
+\method{loo_compare}{default}(x, ..., rank_by = NULL)
 
-\method{print}{compare.loo}(x, ..., digits = 1, p_worse = TRUE)
+\method{print}{compare.loo}(x, ..., digits = 1, p_worse = TRUE, measures = NULL)
 
 \method{print}{compare.loo_ss}(x, ..., digits = 1)
 }
 \arguments{
-\item{x}{An object of class \code{"loo"} or a list of such objects. If a list is
-used then the list names will be used as the model names in the output. See
-\strong{Examples}.}
+\item{x}{An object of class \code{"loo"} or \code{"loo_pred_measure"}, or a list of
+such objects. If a list is used then the list names will be used as the
+model names in the output. See \strong{Examples}.}
 
-\item{...}{Additional objects of class \code{"loo"}, if not passed in as a single
-list.}
+\item{...}{Additional objects of class \code{"loo"} or \code{"loo_pred_measure"}, if not
+passed in as a single list.}
+
+\item{rank_by}{For \code{\link[=loo_pred_measure]{loo_pred_measure()}} comparisons only,
+the bare measure name used to rank models and define the reference model
+for all pairwise differences (default \code{"elpd"}). For example,
+\code{rank_by = "mse"} ranks models by predictive MSE (best/lowest MSE first) and
+computes all measure differences relative to that model on a utility scale
+(higher is better; loss measures such as MSE have their sign flipped).}
 
 \item{digits}{For the print method only, the number of digits to use when
 printing.}
@@ -30,13 +37,36 @@ printing.}
 \item{p_worse}{For the print method only, should we include the normal
 approximation based probability of each model having worse performance than
 the best model? The default is \code{TRUE}.}
+
+\item{measures}{For \code{loo_pred_measure} comparisons only, which measures to
+print diff tables for. \code{NULL} (default) prints only the ranking measure
+(\code{"elpd"} when \code{rank_by} was not set, otherwise \code{rank_by});
+\code{"all"} prints all compared measures; or a character vector of measure
+names (e.g. \code{c("elpd", "mse")}). Printed tables use the column name
+\code{se_diff} for the standard error of the difference even when the data
+frame column is \verb{\{measure\}_se_diff}.}
 }
 \value{
 A data frame with class \code{"compare.loo"} that has its own
 print method. See the \strong{Details} and \strong{Examples} sections.
+
+For classic \code{"loo"} / \code{"waic"} / \code{"kfold"} comparisons, the returned
+columns are unchanged from previous versions.
+
+For \code{\link[=loo_pred_measure]{loo_pred_measure()}} comparisons, the data frame
+additionally contains \verb{\{measure\}_diff} and \verb{\{measure\}_se_diff} columns for
+every predictive measure common to all models (e.g. \code{rmse_diff},
+\code{rmse_se_diff}). ELPD-family measures use \code{elpd_diff} and \code{se_diff}.
+\code{p_worse} and \code{diag_diff} are computed for ELPD only. Per-model PSIS
+diagnostics appear in \code{diag_elpd}. Attributes \code{compare_measures} and
+\code{sign_converted_measures} record which measures were compared and which
+loss measures had their sign flipped for comparison. Attribute \code{rank_by} is
+set when \code{rank_by} was passed explicitly (default ranking is by \code{"elpd"}).
 }
 \description{
-Compare fitted models based on \link[=loo-glossary]{ELPD}.
+Compare fitted models based on \link[=loo-glossary]{ELPD} or, for
+\code{\link[=loo_pred_measure]{loo_pred_measure()}} results, multiple predictive
+performance measures at once.
 }
 \details{
 When comparing two fitted models, we can estimate the difference in their
@@ -97,10 +127,44 @@ computation, then there may be significant bias in \code{elpd_diff} favoring
 models with a large number of high Pareto k values.
 }
 
+\subsection{Comparing \code{loo_pred_measure} objects}{
+
+When all inputs are \code{\link[=loo_pred_measure]{loo_pred_measure()}} objects,
+\code{loo_compare()} computes paired differences for every predictive measure
+present in all models. Models are ranked by \code{rank_by} (default \code{"elpd"});
+the top-ranked model is the reference for all \verb{\{measure\}_diff} columns.
+Measures may use different orientations in their raw form (e.g. ELPD and
+CRPS/RPS are returned on a utility scale where higher is better, while MSE
+and Brier score are loss measures where lower is better). For comparison,
+all \verb{\{measure\}_diff} values are reported on a common utility scale (higher
+is better). Loss measures have their sign flipped from the raw loss
+orientation so that negative \verb{\{measure\}_diff} values indicate worse
+performance than the reference model. Each \verb{*_pred_measure()} result stores
+attribute \code{measure_higher_is_better}, a named list recording the
+\code{higher_is_better} setting used when each measure was computed. When loss measures are compared
+on a utility scale, \code{loo_compare()} emits a short message naming the affected
+measures, for example:
+"For model comparison, differences for mse are reported on a utility scale
+(higher is better)."
+
+\code{p_worse} and \code{diag_diff} are computed for ELPD-family measures only. Other
+measures receive \verb{\{measure\}_diff} and \verb{\{measure\}_se_diff} from paired
+pointwise contributions when the overall estimate is a sum or mean of those
+contributions (using the same standard error formula as \code{se_diff}). For
+measures where pointwise values do not define the overall estimate (e.g.
+\code{r2}, \code{mse}, \code{rmse}), \verb{\{measure\}_diff} is the difference between overall
+estimates and \verb{\{measure\}_se_diff} is \code{NA}. When models were fit with
+different \code{measure} sets, only measures common to all models are compared; a
+warning lists omitted measures. Use \code{print(x, measures = "all")} to display
+diff tables for every compared measure; see \link{loo-glossary} for column
+definitions.
+}
+
 \subsection{Warnings for many model comparisons}{
 
 If more than \eqn{11} models are compared, we internally recompute the model
-differences using the median model by ELPD as the baseline model. We then
+differences using the median model (by ELPD, or by \code{rank_by} for
+\code{loo_pred_measure} comparisons) as the baseline model. We then
 estimate whether the differences in predictive performance are potentially
 due to chance as described by McLatchie and Vehtari (2023). This will flag
 a warning if it is deemed that there is a risk of over-fitting due to the
@@ -126,6 +190,34 @@ loo_compare(list("apple" = loo1, "banana" = loo2, "cherry" = loo3))
 \dontrun{
 # works for waic (and kfold) too
 loo_compare(waic(LL), waic(LL - 10))
+
+# compare multiple predictive measures from loo_pred_measure()
+if (requireNamespace("brms", quietly = TRUE)) {
+  fit1 <- brms::brm(
+    Reaction ~ Days, data = lme4::sleepstudy,
+    refresh = 0, chains = 2, iter = 1000
+  )
+  fit2 <- brms::brm(
+    Reaction ~ poly(Days, 2), data = lme4::sleepstudy,
+    refresh = 0, chains = 2, iter = 1000
+  )
+  pm1 <- loo_pred_measure(
+    loo = loo(fit1, save_psis = TRUE),
+    y = fit1$data$Reaction,
+    mupred = brms::posterior_epred(fit1),
+    measure = c("rmse", "r2")
+  )
+  pm2 <- loo_pred_measure(
+    loo = loo(fit2, save_psis = TRUE),
+    y = fit2$data$Reaction,
+    mupred = brms::posterior_epred(fit2),
+    measure = c("rmse", "r2")
+  )
+  comp <- loo_compare(pm1, pm2)
+  print(comp)                      # ranked by elpd (default)
+  print(comp, measures = "all")    # all measure diff tables
+  loo_compare(pm1, pm2, rank_by = "rmse")
+}
 }
 
 }
diff --git a/man/loo_pred_measure.Rd b/man/loo_pred_measure.Rd
index eeae51b3..89b0d215 100644
--- a/man/loo_pred_measure.Rd
+++ b/man/loo_pred_measure.Rd
@@ -44,7 +44,25 @@ functions (custom measures), e.g. \code{list("rps", my_metric = my_fun)}.
 }
 Custom functions are called with any of \code{y}, \code{ypred}, \code{mupred}, \code{ylp}, and
 \code{log_weights} that appear in their formals, plus arguments from \code{control}.
-They must return a list with  \code{estimates} and \code{pointwise}.}
+They must return a list with  \code{estimates} and \code{pointwise}.
+
+For \code{\link[=loo_compare]{loo_compare()}}, the standard error of a difference is derived from
+paired pointwise differences when the overall estimate is the sum or the
+mean of \code{pointwise}; otherwise it is \code{NA}. A custom measure that is a
+transformation of pointwise quantities (as \code{rmse} is of squared errors) can
+supply its own by setting
+\code{attr(my_fun, "se_diff_fun") <- function(ref, cmp) ...}. That function
+receives one list per model with elements \code{estimate}, \code{se}, \code{pointwise},
+and \code{extra}, always on the measure's natural scale, and must return the
+standard error of the difference as a numeric scalar. Note that the
+function is stored on the result object, so it carries its enclosing
+environment into any saved copy of that object.
+
+\code{extra} is for anything the standard error needs that the pointwise values
+do not carry. Return it as an additional list element \code{extra} from the
+measure function and it is stored alongside the estimates; the built-in
+\code{r2} uses it for the baseline \code{(y_i - mean(y))^2}, which cannot be
+recovered once \code{y} is out of scope.}
 
 \item{loo}{A \code{\link[=loo]{loo()}} result, computed with
 \code{save_psis = TRUE} so that PSIS weights are available for additional
@@ -113,6 +131,25 @@ adding measures incrementally, call \code{\link[=pred_measure]{pred_measure()}}
 an existing result; use \code{save_psis = TRUE} on the initial call so weights
 are stored.
 }
+\section{\code{measure_higher_is_better} attribute}{
+
+All \verb{*_pred_measure()} and \code{\link[=pred_measure]{pred_measure()}} results include attribute
+\code{measure_higher_is_better}: a named list recording the \code{higher_is_better}
+setting used for each measure (\code{NULL}, \code{TRUE}, or \code{FALSE} per bare measure
+name; \code{elpd} is always \code{NULL}). Built-in loss measures such as MSE are stored
+on a loss scale by default (\code{NULL}); pass
+\code{control = list(mse = list(higher_is_better = TRUE))} to store values on a
+utility scale. When measures are added incrementally with \code{\link[=pred_measure]{pred_measure()}},
+the attribute is updated for newly computed measures.
+
+Attribute \code{measure_compare_meta} records per-measure comparison metadata
+(\code{higher_is_better}, \code{loss}, and \code{diff_method}) used by \code{\link[=loo_compare]{loo_compare()}}.
+Built-in measures take \code{loss} and \code{diff_method} from the package measure
+registry; custom measures default to \code{diff_method = "auto"} and \code{loss = FALSE}
+(utility scale). For custom loss measures, use \code{higher_is_better = TRUE} in
+\code{control} so \code{\link[=loo_compare]{loo_compare()}} interprets the sign correctly. See \link{loo-glossary}.
+}
+
 \examples{
 \donttest{
 if (requireNamespace("brms", quietly = TRUE)) {
diff --git a/man/pred_measure.Rd b/man/pred_measure.Rd
index 380ca194..4096d7d1 100644
--- a/man/pred_measure.Rd
+++ b/man/pred_measure.Rd
@@ -44,7 +44,25 @@ functions (custom measures), e.g. \code{list("rps", my_metric = my_fun)}.
 }
 Custom functions are called with any of \code{y}, \code{ypred}, \code{mupred}, \code{ylp}, and
 \code{log_weights} that appear in their formals, plus arguments from \code{control}.
-They must return a list with  \code{estimates} and \code{pointwise}.}
+They must return a list with  \code{estimates} and \code{pointwise}.
+
+For \code{\link[=loo_compare]{loo_compare()}}, the standard error of a difference is derived from
+paired pointwise differences when the overall estimate is the sum or the
+mean of \code{pointwise}; otherwise it is \code{NA}. A custom measure that is a
+transformation of pointwise quantities (as \code{rmse} is of squared errors) can
+supply its own by setting
+\code{attr(my_fun, "se_diff_fun") <- function(ref, cmp) ...}. That function
+receives one list per model with elements \code{estimate}, \code{se}, \code{pointwise},
+and \code{extra}, always on the measure's natural scale, and must return the
+standard error of the difference as a numeric scalar. Note that the
+function is stored on the result object, so it carries its enclosing
+environment into any saved copy of that object.
+
+\code{extra} is for anything the standard error needs that the pointwise values
+do not carry. Return it as an additional list element \code{extra} from the
+measure function and it is stored alongside the estimates; the built-in
+\code{r2} uses it for the baseline \code{(y_i - mean(y))^2}, which cannot be
+recovered once \code{y} is out of scope.}
 
 \item{predperf}{An existing predictive measure object (class
 \code{"pred_measure"}) to update. When supplied, base density summaries and
@@ -69,7 +87,8 @@ measure's summary function (e.g. \code{list(new_measure = list(add_arg = 10))}).
 An updated object of the same class as \code{predperf}, with new rows in
 \code{estimates} and columns in \code{pointwise} for each requested measure. Base
 summaries (\code{elpd} and LOO/k-fold complexity terms such as \code{p_loo}) are not
-recomputed.
+recomputed. Attribute \code{measure_higher_is_better} is updated for any newly added
+measures.
 }
 \description{
 Extend a \code{"pred_measure"} object with additional measures \strong{without
@@ -100,6 +119,25 @@ When extending a LOO result, ensure the initial call used \code{save_psis = TRUE
 (or that \code{predperf} already contains a \code{psis_object}) so LOO weights are
 available for additional measures.
 }
+\section{\code{measure_higher_is_better} attribute}{
+
+All \verb{*_pred_measure()} and \code{\link[=pred_measure]{pred_measure()}} results include attribute
+\code{measure_higher_is_better}: a named list recording the \code{higher_is_better}
+setting used for each measure (\code{NULL}, \code{TRUE}, or \code{FALSE} per bare measure
+name; \code{elpd} is always \code{NULL}). Built-in loss measures such as MSE are stored
+on a loss scale by default (\code{NULL}); pass
+\code{control = list(mse = list(higher_is_better = TRUE))} to store values on a
+utility scale. When measures are added incrementally with \code{\link[=pred_measure]{pred_measure()}},
+the attribute is updated for newly computed measures.
+
+Attribute \code{measure_compare_meta} records per-measure comparison metadata
+(\code{higher_is_better}, \code{loss}, and \code{diff_method}) used by \code{\link[=loo_compare]{loo_compare()}}.
+Built-in measures take \code{loss} and \code{diff_method} from the package measure
+registry; custom measures default to \code{diff_method = "auto"} and \code{loss = FALSE}
+(utility scale). For custom loss measures, use \code{higher_is_better = TRUE} in
+\code{control} so \code{\link[=loo_compare]{loo_compare()}} interprets the sign correctly. See \link{loo-glossary}.
+}
+
 \examples{
 \donttest{
 if (requireNamespace("brms", quietly = TRUE)) {
diff --git a/man/pred_measure_params.Rd b/man/pred_measure_params.Rd
index 60ba9954..26a05399 100644
--- a/man/pred_measure_params.Rd
+++ b/man/pred_measure_params.Rd
@@ -54,7 +54,25 @@ functions (custom measures), e.g. \code{list("rps", my_metric = my_fun)}.
 }
 Custom functions are called with any of \code{y}, \code{ypred}, \code{mupred}, \code{ylp}, and
 \code{log_weights} that appear in their formals, plus arguments from \code{control}.
-They must return a list with  \code{estimates} and \code{pointwise}.}
+They must return a list with  \code{estimates} and \code{pointwise}.
+
+For \code{\link[=loo_compare]{loo_compare()}}, the standard error of a difference is derived from
+paired pointwise differences when the overall estimate is the sum or the
+mean of \code{pointwise}; otherwise it is \code{NA}. A custom measure that is a
+transformation of pointwise quantities (as \code{rmse} is of squared errors) can
+supply its own by setting
+\code{attr(my_fun, "se_diff_fun") <- function(ref, cmp) ...}. That function
+receives one list per model with elements \code{estimate}, \code{se}, \code{pointwise},
+and \code{extra}, always on the measure's natural scale, and must return the
+standard error of the difference as a numeric scalar. Note that the
+function is stored on the result object, so it carries its enclosing
+environment into any saved copy of that object.
+
+\code{extra} is for anything the standard error needs that the pointwise values
+do not carry. Return it as an additional list element \code{extra} from the
+measure function and it is stored alongside the estimates; the built-in
+\code{r2} uses it for the baseline \code{(y_i - mean(y))^2}, which cannot be
+recovered once \code{y} is out of scope.}
 
 \item{predperf}{An existing predictive measure object (class
 \code{"pred_measure"}) to update. When supplied, base density summaries and
diff --git a/man/test_pred_measure.Rd b/man/test_pred_measure.Rd
index 0b3710e9..d499fa3c 100644
--- a/man/test_pred_measure.Rd
+++ b/man/test_pred_measure.Rd
@@ -47,7 +47,25 @@ functions (custom measures), e.g. \code{list("rps", my_metric = my_fun)}.
 }
 Custom functions are called with any of \code{y}, \code{ypred}, \code{mupred}, \code{ylp}, and
 \code{log_weights} that appear in their formals, plus arguments from \code{control}.
-They must return a list with  \code{estimates} and \code{pointwise}.}
+They must return a list with  \code{estimates} and \code{pointwise}.
+
+For \code{\link[=loo_compare]{loo_compare()}}, the standard error of a difference is derived from
+paired pointwise differences when the overall estimate is the sum or the
+mean of \code{pointwise}; otherwise it is \code{NA}. A custom measure that is a
+transformation of pointwise quantities (as \code{rmse} is of squared errors) can
+supply its own by setting
+\code{attr(my_fun, "se_diff_fun") <- function(ref, cmp) ...}. That function
+receives one list per model with elements \code{estimate}, \code{se}, \code{pointwise},
+and \code{extra}, always on the measure's natural scale, and must return the
+standard error of the difference as a numeric scalar. Note that the
+function is stored on the result object, so it carries its enclosing
+environment into any saved copy of that object.
+
+\code{extra} is for anything the standard error needs that the pointwise values
+do not carry. Return it as an additional list element \code{extra} from the
+measure function and it is stored alongside the estimates; the built-in
+\code{r2} uses it for the baseline \code{(y_i - mean(y))^2}, which cannot be
+recovered once \code{y} is out of scope.}
 
 \item{group_ids}{Optional vector of group identifiers for grouped summaries
 (reserved; not yet implemented).}
@@ -79,6 +97,25 @@ training fit.
 The base summary \code{elpd_test} is computed from \code{ylp_test} on the holdout
 observations only.
 }
+\section{\code{measure_higher_is_better} attribute}{
+
+All \verb{*_pred_measure()} and \code{\link[=pred_measure]{pred_measure()}} results include attribute
+\code{measure_higher_is_better}: a named list recording the \code{higher_is_better}
+setting used for each measure (\code{NULL}, \code{TRUE}, or \code{FALSE} per bare measure
+name; \code{elpd} is always \code{NULL}). Built-in loss measures such as MSE are stored
+on a loss scale by default (\code{NULL}); pass
+\code{control = list(mse = list(higher_is_better = TRUE))} to store values on a
+utility scale. When measures are added incrementally with \code{\link[=pred_measure]{pred_measure()}},
+the attribute is updated for newly computed measures.
+
+Attribute \code{measure_compare_meta} records per-measure comparison metadata
+(\code{higher_is_better}, \code{loss}, and \code{diff_method}) used by \code{\link[=loo_compare]{loo_compare()}}.
+Built-in measures take \code{loss} and \code{diff_method} from the package measure
+registry; custom measures default to \code{diff_method = "auto"} and \code{loss = FALSE}
+(utility scale). For custom loss measures, use \code{higher_is_better = TRUE} in
+\code{control} so \code{\link[=loo_compare]{loo_compare()}} interprets the sign correctly. See \link{loo-glossary}.
+}
+
 \examples{
 \donttest{
 if (requireNamespace("brms", quietly = TRUE)) {
diff --git a/notes/developer-notes.Rmd b/notes/developer-notes.Rmd
index ae190ae4..2298e5c7 100644
--- a/notes/developer-notes.Rmd
+++ b/notes/developer-notes.Rmd
@@ -20,10 +20,10 @@ knitr::opts_chunk$set(
 
 > **Status:** In Progress  
 > **Base branch:** `loo-v3.0.0`  
-> **Compare branch:** `pred_measure`  
+> **Compare branch:** `pred_measure` (+ `integrate-loo_compare`)  
 > **Related PR:** [#363](https://github.com/stan-dev/loo/pull/363)  
 > **Contributors:** @florence-bockting, @avehtari, @VisruthSK, @jgabry  
-> **Last updated:** 2026-07-03
+> **Last updated:** 2026-07-07
 
 These notes document internal design decisions and ongoing work for the
 `pred_measure` feature. This PR **adds** the new API.
@@ -37,7 +37,7 @@ These notes document internal design decisions and ongoing work for the
 | Scoring rules (`measure_rps`) | Done |
 | Documentation | In progress |
 | `group_ids` grouping | Not started |
-| `loo_compare` integration | Not started |
+| `loo_compare` integration | Done (`integrate-loo_compare`) |
 
 ## Scope of this PR (`pred_measure` vs `loo-v3.0.0`)
 
@@ -52,6 +52,7 @@ These notes document internal design decisions and ongoing work for the
 - `vignettes/migration-guide.Rmd`
 - Website-only articles: `overview-measures.Rmd`, `pred-measure-workflow.Rmd`
 - Test suite + pre-fitted fixtures + `test_data_generation.R`
+- `loo_compare()` multi-measure path for `loo_pred_measure` objects
 
 ### Changed on existing code (implementations retained)
 
@@ -61,6 +62,9 @@ These notes document internal design decisions and ongoing work for the
 - `elpd()` — refactored to `.elpd_matrix_impl()` to avoid double deprecation
   warnings
 - Minor doc cross-references in `compare.R`, `psislw.R`
+- `loo_compare()` — extended for `loo_pred_measure` objects: `rank_by`, multi-measure
+  paired diffs, updated `print.compare.loo(measures = ...)`; classic `loo` path unchanged
+- `R/loo-glossary.R` — multi-measure comparison columns (`{measure}_diff`, `rank_by`, etc.)
 - `NEWS.md`, `NAMESPACE`, `_pkgdown.yml`, pkgdown CI workflow
 
 ## Design decisions (resolved)
@@ -111,6 +115,29 @@ Design choices **internal to `pred_measure`**:
 - [ ] Verify `kfold_pred_measure()` with categorical/multinomial examples
   end-to-end (penguins fixture exists; confirm test/doc coverage)
 
+### D4: `loo_compare()` for `loo_pred_measure` objects
+
+**Decision:** Extend existing `loo_compare()`, not a separate function.
+
+- When all inputs are `loo_pred_measure` objects, compute paired differences for
+  every measure common to all models
+- Rank models by `rank_by` (default `"elpd"`); top-ranked model is the reference
+  for all `{measure}_diff` columns
+- ELPD-family measures keep `elpd_diff` / `se_diff`; other measures use
+  `{measure}_diff` / `{measure}_se_diff`
+- `p_worse` and `diag_diff` apply to ELPD only; `diag_elpd` per model as before
+- Loss measures (MSE, RMSE, MAE, IC) compared on a utility scale (negated) so
+  worse models have negative diffs, consistent with ELPD
+- Pointwise SEs use the same paired formula as ELPD when the overall estimate is
+  a sum or mean of pointwise contributions; otherwise `{measure}_se_diff` is `NA`
+  (e.g. `r2`, `mse`, `rmse`)
+- Reuse `elpd_diffs`, `se_elpd_diff`, `diag_diff`, `diag_elpd`, and many-model
+  order-statistic check (with `rank_by` when applicable)
+- `print.compare.loo(measures = ...)` shows one or all measure diff tables
+
+Implemented on branch `integrate-loo_compare`; tests in `test_compare.R` with
+fixture `test_data_roaches_compare.Rds`.
+
 ## Open decisions
 
 ### D1: Sign convention for pointwise estimates
@@ -138,7 +165,7 @@ Design choices **internal to `pred_measure`**:
 - [x] In `*_pred_measure()`, compute `elpd` as base when `ylp` supplied;
   require explicit `measure = "ic"` for information criterion
 - [x] Document and test deprecated vs new API comparisons *(see appendix)*
-- [ ] Provide an interface to `loo_compare` and verify consistency
+- [x] Provide an interface to `loo_compare` and verify consistency
 - [ ] Resolve `r_eff` handling *(see D3)*
 
 ### Implementation
@@ -156,7 +183,8 @@ Design choices **internal to `pred_measure`**:
 - [x] Online-only articles published via `_pkgdown.yml`
 - [ ] Formula derivations article (`pred_measure-formulas.Rmd`)
 - [ ] Detailed per-measure descriptions (derivations where appropriate)
-- [ ] Extend glossary (`R/loo-glossary.R`) — measure, metric, score, utility, loss
+- [x] Extend glossary (`R/loo-glossary.R`) — multi-measure `loo_compare` columns
+- [ ] Extend glossary further — measure, metric, score, utility, loss (general terms)
 
 ### Grouping via `group_ids`
 
@@ -170,8 +198,8 @@ Design choices **internal to `pred_measure`**:
 
 - Rename `ic` → `information_criteria` for clarity?
 - Should `measure_elpd()` also return `ic`, or keep them separate?
-- What defines class `"loo"` on measure objects? (e.g. deprecated `elpd_generic`
-  inherits `"loo"`)
+- What defines class `"loo"` on measure objects? `loo_pred_measure` inherits `"loo"`
+  (see `integrate-loo_compare`); deprecated `elpd_generic` also inherits `"loo"`.
 - Should `elpd` always be computed when `ylp` is supplied, or allow
   `loo_pred_measure()` for non-ELPD measures only?
 
diff --git a/notes/developer-notes.md b/notes/developer-notes.md
index 0dd96f5c..28cb6b5d 100644
--- a/notes/developer-notes.md
+++ b/notes/developer-notes.md
@@ -14,10 +14,10 @@ Developer Notes: `pred_measure` Feature
 
 > **Status:** In Progress  
 > **Base branch:** `loo-v3.0.0`  
-> **Compare branch:** `add-pred-measure`  
+> **Compare branch:** `pred_measure` (+ `integrate-loo_compare`)  
 > **Related PR:** [\#363](https://github.com/stan-dev/loo/pull/363)  
 > **Contributors:** @florence-bockting, @avehtari, @VisruthSK, @jgabry  
-> **Last updated:** 2026-07-03
+> **Last updated:** 2026-07-07
 
 These notes document internal design decisions and ongoing work for the
 `pred_measure` feature. This PR **adds** the new API on top of
@@ -39,7 +39,7 @@ For the merge summary, see the PR description
 | Scoring rules (`measure_rps`) | Done                   |
 | Documentation                 | In progress            |
 | `group_ids` grouping          | Not started            |
-| `loo_compare` integration     | Not started            |
+| `loo_compare` integration     | Done (`integrate-loo_compare`) |
 
 ------------------------------------------------------------------------
 
@@ -59,6 +59,8 @@ For the merge summary, see the PR description
 - Website-only articles: `overview-measures.Rmd`,
   `pred-measure-workflow.Rmd`
 - Test suite + pre-fitted fixtures + `test_data_generation.R`
+- `loo_compare()` multi-measure path for `loo_pred_measure` objects
+  (`integrate-loo_compare`)
 
 ### Changed on existing code (implementations retained)
 
@@ -69,6 +71,11 @@ For the merge summary, see the PR description
 - `elpd()` — refactored to `.elpd_matrix_impl()` to avoid double
   deprecation warnings
 - Minor doc cross-references in `compare.R`, `psislw.R`
+- `loo_compare()` — extended for `loo_pred_measure` objects: `rank_by`,
+  multi-measure paired diffs, updated `print.compare.loo(measures = ...)`;
+  classic `loo` path unchanged
+- `R/loo-glossary.R` — multi-measure comparison columns (`{measure}_diff`,
+  `rank_by`, etc.)
 - `NEWS.md`, `NAMESPACE`, `_pkgdown.yml`, pkgdown CI workflow
 
 ------------------------------------------------------------------------
@@ -127,17 +134,53 @@ Design choices **internal to `pred_measure`** (not a migration from
   examples end-to-end (penguins fixture exists; confirm test/doc
   coverage)
 
+### D4: `loo_compare()` for `loo_pred_measure` objects
+
+**Decision:** Extend existing `loo_compare()`
+
+- When all inputs are `loo_pred_measure` objects, compute paired
+  differences for every measure common to all models
+- Rank models by `rank_by` (default `"elpd"`); top-ranked model is the
+  reference for all `{measure}_diff` columns
+- ELPD-family measures keep `elpd_diff` / `se_diff`; other measures use
+  `{measure}_diff` / `{measure}_se_diff`
+- `p_worse` and `diag_diff` apply to ELPD only; `diag_elpd` per model as
+  before
+- Loss measures (MSE, RMSE, MAE, IC, Brier score, SRPS) compared on a common
+  utility scale (higher is better): sign flipped from the raw loss orientation
+  so worse models have negative diffs, consistent with ELPD. Orientation is
+  read from `measure_compare_meta` and `measure_higher_is_better` on each
+  `*_pred_measure()` result; attribute `sign_converted_measures` records
+  affected measures. A short message is emitted at compare time; full
+  interpretation is in `?loo_compare` / `?loo-glossary`.
+- Pointwise SEs use the same paired formula as ELPD when the overall
+  estimate is a sum or mean of pointwise contributions; otherwise
+  `{measure}_se_diff` is `NA` (e.g. `r2`, `mse`, `rmse`)
+- Reuse `elpd_diffs`, `se_elpd_diff`, `diag_diff`, `diag_elpd`, and
+  many-model order-statistic check (with `rank_by` when applicable)
+- `print.compare.loo(measures = ...)` shows one or all measure diff tables
+
+Implemented on branch `integrate-loo_compare`; tests in `test_compare.R`
+with fixture `test_data_roaches_compare.Rds`.
+
 ------------------------------------------------------------------------
 
 ## Open decisions
 
 ### D1: Sign convention for pointwise estimates
 
-- **Context:** Measures differ in orientation (`rps`: lower is better;
-  `srps`: higher is better). Aligning orientations may help comparisons.
-- **Options:** `lower_is_better`, `orientation = "utility" / "loss"`,
-  `revert_sign` (currently internal on some `measure_*()` functions)
-- **Decision:** *pending*
+- **Context:** Measures differ in orientation (e.g. ELPD/CRPS on a utility
+  scale; MSE and Brier score as losses). `loo_compare()` aligns them for
+  paired differences.
+- **Decision (for `loo_compare`):** Each `*_pred_measure()` result stores
+  `higher_is_better` per measure in `measure_higher_is_better` and comparison
+  metadata in `measure_compare_meta` (`higher_is_better`, `loss`,
+  `diff_method`). Built-in loss measures are sign-flipped for utility-scale
+  `{measure}_diff` when stored on a loss scale (`higher_is_better = NULL` with
+  `loss = TRUE`).
+- **Still open:** Whether to expose orientation metadata on `*_pred_measure()`
+  results themselves beyond the attributes above (e.g. when
+  `higher_is_better = TRUE` in `control`).
 
 ### D3: Handling of `r_eff`
 
@@ -163,7 +206,7 @@ Design choices **internal to `pred_measure`** (not a migration from
   require explicit `measure = "ic"` for information criterion
 - [x] Document and test deprecated vs new API comparisons *(see
   appendix)*
-- [ ] Provide an interface to `loo_compare` and verify consistency
+- [x] Provide an interface to `loo_compare` and verify consistency
 - [ ] Resolve `r_eff` handling *(see D3)*
 
 ### Implementation
@@ -181,8 +224,9 @@ Design choices **internal to `pred_measure`** (not a migration from
 - [x] Online-only articles published via `_pkgdown.yml`
 - [ ] Formula derivations article (`pred_measure-formulas.Rmd`)
 - [ ] Detailed per-measure descriptions (derivations where appropriate)
-- [ ] Extend glossary (`R/loo-glossary.R`) — measure, metric, score,
-  utility, loss
+- [x] Extend glossary (`R/loo-glossary.R`) — multi-measure `loo_compare` columns
+- [ ] Extend glossary further — measure, metric, score, utility, loss
+  (general terms)
 
 ### Grouping via `group_ids`
 
@@ -200,8 +244,9 @@ Design choices **internal to `pred_measure`** (not a migration from
 
 - Rename `ic` → `information_criteria` for clarity?
 - Should `measure_elpd()` also return `ic`, or keep them separate?
-- What defines class `"loo"` on measure objects? (e.g. deprecated
-  `elpd_generic` inherits `"loo"`)
+- What defines class `"loo"` on measure objects? `loo_pred_measure` inherits
+  `"loo"` (see `integrate-loo_compare`); deprecated `elpd_generic` also
+  inherits `"loo"`.
 - Should `elpd` always be computed when `ylp` is supplied, or allow
   `loo_pred_measure()` for non-ELPD measures only?
 
@@ -236,7 +281,7 @@ estimators.
 
 | Deprecated        | New workflow                                    | Notes                               |
 |-------------------|-------------------------------------------------|-------------------------------------|
-| `crps(x, x2, y)`  | `measure_rps(y, ypred = x, revert_sign = TRUE)` | Sign flip on unscaled score         |
+| `crps(x, x2, y)`  | `measure_rps(y, ypred = x, higher_is_better = FALSE)` | Sign flip on unscaled score         |
 | `scrps(x, x2, y)` | `measure_srps(y, ypred = x)`                    | Same sign convention                |
 | `loo_crps(...)`   | `loo_pred_measure(..., measure = "rps")`        | Additional LOO weighting difference |
 | `loo_scrps(...)`  | `loo_pred_measure(..., measure = "srps")`       | Additional LOO weighting difference |
@@ -245,7 +290,7 @@ estimators.
 
 1.  **Sign convention (unscaled only).** `crps()` returns
     `0.5·EXX − EXy` (utility: higher is better). Default `measure_rps()`
-    negates this; use `revert_sign = TRUE` to match `crps()`. Scaled
+    negates this; use `higher_is_better = FALSE` to match `crps()`. Scaled
     scores (`scrps` / `measure_srps`) already share the formula
     `−EXy/EXX − 0.5·log(EXX)`.
 
@@ -295,7 +340,7 @@ replications</figcaption>
 </figure>
 
 *Figure: 200 simulations (S = 100, n = 30). Left — `crps()` vs
-`measure_rps(revert_sign = TRUE)`; right — `scrps()` vs
+`measure_rps(higher_is_better = FALSE)`; right — `scrps()` vs
 `measure_srps()`.*
 
 #### LOO outcome comparison
diff --git a/tests/testthat/_snaps/compare.md b/tests/testthat/_snaps/compare.md
index 8abb5cf3..c0735f69 100644
--- a/tests/testthat/_snaps/compare.md
+++ b/tests/testthat/_snaps/compare.md
@@ -1,3 +1,116 @@
+# loo_compare works with three loo_pred_measure models
+
+    Code
+      print(comp)
+    Message
+      Models ranked by mae (reference: C).
+    Output
+       model mae_diff se_diff        diag_elpd
+           C      0.0     0.0 13 k_psis > 0.62
+           B     -0.1     0.7  9 k_psis > 0.62
+           A     -7.6     1.7 25 k_psis > 0.62
+    Message
+      
+      Diagnostic flags present.
+      See ?`loo-glossary` (sections `diag_diff` and `diag_elpd`)
+      or https://mc-stan.org/loo/reference/loo-glossary.html.
+      
+      Other measures compared: elpd, r2. Use print(x, measures = "all").
+
+# loo_compare informs when measure signs are converted
+
+    Code
+      comp <- loo_compare(pm1, pm2)
+    Message
+      For model comparison, differences for mse is reported on a utility scale (higher is better).
+
+# print.compare.loo works for loo_pred_measure comparisons
+
+    Code
+      print(comp)
+    Output
+       model elpd_diff se_diff p_worse diag_diff        diag_elpd
+          m3       0.0     0.0      NA           13 k_psis > 0.62
+          m2    -219.0   221.6    0.84            9 k_psis > 0.62
+          m1   -2764.4   634.3    1.00           25 k_psis > 0.62
+    Message
+      
+      Diagnostic flags present.
+      See ?`loo-glossary` (sections `diag_diff` and `diag_elpd`)
+      or https://mc-stan.org/loo/reference/loo-glossary.html.
+      
+      Other measures compared: r2, mae. Use print(x, measures = "all").
+
+---
+
+    Code
+      print(comp, measures = "all", digits = 2)
+    Output
+      
+      -- elpd (vs m3) --
+       model elpd_diff se_diff p_worse diag_diff        diag_elpd
+          m3      0.00    0.00      NA           13 k_psis > 0.62
+          m2   -218.98  221.65    0.84            9 k_psis > 0.62
+          m1  -2764.42  634.26    1.00           25 k_psis > 0.62
+      
+      -- r2 (vs m3) --
+       model r2_diff se_diff
+          m3    0.00    0.00
+          m2    0.01    0.07
+          m1   -0.20    0.13
+      
+      -- mae (vs m3) --
+       model mae_diff se_diff
+          m3     0.00    0.00
+          m2    -0.08    0.71
+          m1    -7.64    1.68
+    Message
+      
+      Diagnostic flags present.
+      See ?`loo-glossary` (sections `diag_diff` and `diag_elpd`)
+      or https://mc-stan.org/loo/reference/loo-glossary.html.
+
+---
+
+    Code
+      print(comp, measures = c("r2", "mae"))
+    Output
+      
+      -- r2 (vs m3) --
+       model r2_diff se_diff
+          m3     0.0     0.0
+          m2     0.0     0.1
+          m1    -0.2     0.1
+      
+      -- mae (vs m3) --
+       model mae_diff se_diff
+          m3      0.0     0.0
+          m2     -0.1     0.7
+          m1     -7.6     1.7
+    Message
+      
+      Diagnostic flags present.
+      See ?`loo-glossary` (sections `diag_diff` and `diag_elpd`)
+      or https://mc-stan.org/loo/reference/loo-glossary.html.
+
+---
+
+    Code
+      print(comp_mae)
+    Message
+      Models ranked by mae (reference: m2).
+    Output
+       model mae_diff se_diff        diag_elpd
+          m2      0.0     0.0  9 k_psis > 0.62
+          m1     -7.6     1.5 25 k_psis > 0.62
+    Message
+      
+      Diagnostic flags present.
+      See ?`loo-glossary` (sections `diag_diff` and `diag_elpd`)
+      or https://mc-stan.org/loo/reference/loo-glossary.html.
+      
+      Other measures compared: elpd, r2. Use print(x, measures = "all").
+
 # loo_compare returns expected results (2 models)
 
     WAoAAAACAAQEAgACAwAAAAMTAAAADAAAABAAAAACAAQACQAAAAZtb2RlbDEABAAJAAAABm1v
diff --git a/tests/testthat/data-for-tests/test_data_generation.R b/tests/testthat/data-for-tests/test_data_generation.R
index 3b05b39a..a7fa7a76 100644
--- a/tests/testthat/data-for-tests/test_data_generation.R
+++ b/tests/testthat/data-for-tests/test_data_generation.R
@@ -125,6 +125,52 @@ get_roaches_res <- function() {
   )
 }
 
+get_roaches_compare_res <- function() {
+  data(roaches, package = "rstanarm")
+  roaches$sqrt_roach1 <- sqrt(roaches$roach1)
+  
+  fit_p <- brm(
+    y ~ sqrt_roach1 + treatment + senior + offset(log(exposure2)),
+    data = roaches,
+    family = poisson,
+    prior = prior(normal(0, 1), class = b),
+    chains = 2,
+    iter = 400,
+    refresh = 0,
+    seed = SEED
+  )
+
+  fit_p <- add_criterion(
+    fit_p,
+    criterion = "loo",
+    moment_match = TRUE,
+    overwrite = TRUE
+  )
+
+  fit_p_m1 <- update(fit_p, formula = y ~ treatment + senior) |>
+    add_criterion(criterion = "loo", moment_match = TRUE)
+  fit_p_m2 <- update(fit_p, formula = y ~ sqrt_roach1 + senior)  |>
+    add_criterion(criterion = "loo", moment_match = TRUE)
+  fit_p_m3 <- update(fit_p, formula = y ~ sqrt_roach1 + treatment) |>
+    add_criterion(criterion = "loo", moment_match = TRUE)
+
+  return(list(
+    y = fit_p$data$y,
+    loo_p = fit_p$loo,
+    mupred = brms::posterior_epred(fit_p),
+    ylp = brms::log_lik(fit_p),
+    loo_p_m1 = fit_p_m1$loo,
+    mupred_m1 = brms::posterior_epred(fit_p_m1),
+    ylp_m1 = brms::log_lik(fit_p_m1),
+    loo_p_m2 = fit_p_m2$loo,
+    mupred_m2 = brms::posterior_epred(fit_p_m2),
+    ylp_m2 = brms::log_lik(fit_p_m2),
+    loo_p_m3 = fit_p_m3$loo,
+    mupred_m3 = brms::posterior_epred(fit_p_m3),
+    ylp_m3 = brms::log_lik(fit_p_m3)
+  ))
+}
+
 get_sleep_test_train_res <- function() {
   # specifically for testing test_pred_measure
   data("sleepstudy", package = "lme4")
@@ -256,10 +302,12 @@ generate_test_data <- function(vignettes_only = FALSE) {
   full_binomial <- get_binomial_res()
   full_sleep <- get_sleep_res()
   full_sleep_test <- get_sleep_test_train_res()
+  full_roaches_compare <- get_roaches_compare_res()
 
   if (!vignettes_only) {
     test_path <- "tests/testthat/data-for-tests/"
     saveRDS(full_roaches$res, paste0(test_path, "test_data_roaches.Rds"))
+    saveRDS(full_roaches_compare, paste0(test_path, "test_data_roaches_compare.Rds"))
     saveRDS(full_binary$res, paste0(test_path, "test_data_binary.Rds"))
     saveRDS(full_penguins$res, paste0(test_path, "test_data_penguins.Rds"))
     saveRDS(full_binomial$res, paste0(test_path, "test_data_binomial.Rds"))
diff --git a/tests/testthat/data-for-tests/test_data_roaches_compare.Rds b/tests/testthat/data-for-tests/test_data_roaches_compare.Rds
new file mode 100644
index 00000000..efcf2dac
Binary files /dev/null and b/tests/testthat/data-for-tests/test_data_roaches_compare.Rds differ
diff --git a/tests/testthat/test_compare.R b/tests/testthat/test_compare.R
index 4a00a10a..8d522ea5 100644
--- a/tests/testthat/test_compare.R
+++ b/tests/testthat/test_compare.R
@@ -10,7 +10,7 @@ test_that("loo_compare throws appropriate errors", {
   w3 <- suppressWarnings(waic(LLarr[,, -1]))
   w4 <- suppressWarnings(waic(LLarr[,, -(1:2)]))
 
-  expect_error(loo_compare(2, 3), "must be a list if not a 'loo' object")
+  expect_error(loo_compare(2, 3), "must be a list if not a 'loo' or 'pred_measure' object")
   expect_error(
     loo_compare(w1, w2, x = list(w1, w2)),
     "If 'x' is a list then '...' should not be specified"
@@ -35,6 +35,987 @@ test_that("loo_compare throws appropriate errors", {
   )
 })
 
+test_that("loo_compare dispatches loo_pred_measure inputs", {
+  res <- readRDS("data-for-tests/test_data_roaches_compare.Rds")
+  pm1 <- loo_pred_measure(
+    loo = res$loo_p_m1,
+    y = res$y,
+    mupred = res$mupred_m1,
+    ylp = res$ylp_m1,
+    measure = c("r2", "mse")
+  )
+  pm2 <- loo_pred_measure(
+    loo = res$loo_p_m1,
+    y = res$y,
+    mupred = res$mupred_m1,
+    ylp = res$ylp_m1,
+    measure = c("r2", "mse")
+  )
+
+  comp <- suppressMessages(loo_compare(pm1, pm2))
+  expect_s3_class(comp, "compare.loo")
+  expect_null(attr(comp, "rank_by"))
+  expect_true(all(c("elpd_diff", "se_diff", "p_worse", "diag_diff") %in% colnames(comp)))
+  expect_true(all(c("r2_diff", "r2_se_diff", "mse_diff", "mse_se_diff") %in% colnames(comp)))
+  expect_false(anyNA(comp$r2_se_diff))
+  expect_false(anyNA(comp$mse_se_diff))
+  expect_false("r2_loo_diff" %in% colnames(comp))
+  expect_false("mse_p_worse" %in% colnames(comp))
+
+  expect_error(
+    loo_compare(w1, pm1),
+    "Cannot mix 'loo_pred_measure' objects with other 'loo' objects",
+    fixed = TRUE
+  )
+  expect_error(
+    loo_compare(
+      insample_pred_measure(ylp = res$ylp_m1),
+      insample_pred_measure(ylp = res$ylp_m1)
+    ),
+    "requires 'loo_pred_measure' objects",
+    fixed = TRUE
+  )
+  expect_error(loo_compare(pm1), "requires at least two models", fixed = TRUE)
+  expect_null(attr(loo_compare(w1, w2), "rank_by"))
+})
+
+test_that("loo_compare warns when predictive measures differ across models", {
+  res <- readRDS("data-for-tests/test_data_roaches_compare.Rds")
+  pm1 <- loo_pred_measure(
+    loo = res$loo_p_m1,
+    y = res$y,
+    mupred = res$mupred_m1,
+    ylp = res$ylp_m1,
+    measure = c("r2", "mse")
+  )
+  pm2 <- loo_pred_measure(
+    loo = res$loo_p_m2,
+    y = res$y,
+    mupred = res$mupred_m2,
+    ylp = res$ylp_m2,
+    measure = c("r2", "mae")
+  )
+
+  expect_warning(
+    comp <- suppressMessages(loo_compare(list(m1 = pm1, m2 = pm2))),
+    "Omitted measures: mae \\(m2\\), mse \\(m1\\)"
+  )
+  expect_equal(attr(comp, "compare_measures"), c("elpd", "r2"))
+  expect_false("mse_diff" %in% colnames(comp))
+  expect_false("mae_diff" %in% colnames(comp))
+})
+
+test_that("loo_compare works with three loo_pred_measure models", {
+  res <- readRDS("data-for-tests/test_data_roaches_compare.Rds")
+  pm1 <- loo_pred_measure(
+    loo = res$loo_p_m1,
+    y = res$y,
+    mupred = res$mupred_m1,
+    ylp = res$ylp_m1,
+    measure = c("r2", "mae")
+  )
+  pm2 <- loo_pred_measure(
+    loo = res$loo_p_m2,
+    y = res$y,
+    mupred = res$mupred_m2,
+    ylp = res$ylp_m2,
+    measure = c("r2", "mae")
+  )
+  pm3 <- loo_pred_measure(
+    loo = res$loo_p_m3,
+    y = res$y,
+    mupred = res$mupred_m3,
+    ylp = res$ylp_m3,
+    measure = c("r2", "mae")
+  )
+
+  comp <- loo_compare(
+    list("A" = pm1, "B" = pm2, "C" = pm3),
+    rank_by = "mae"
+  )
+  expect_snapshot(print(comp))
+  expect_equal(nrow(comp), 3L)
+  expect_equal(comp$model, c("C", "B", "A"))
+  expect_equal(attr(comp, "rank_by"), "mae")
+  expect_equal(attr(comp, "compare_measures"), c("elpd", "r2", "mae"))
+  expect_equal(comp$mae_diff[1L], 0)
+  expect_true(all(comp$mae_diff[-1L] < 0))
+  expect_true(all(comp$elpd_diff[-1L] <= 0))
+  expect_equal(attr(comp, "sign_converted_measures"), c("mae"))
+})
+
+test_that("loo_compare informs when measure signs are converted", {
+  res <- readRDS("data-for-tests/test_data_roaches_compare.Rds")
+  pm1 <- loo_pred_measure(
+    loo = res$loo_p_m1,
+    y = res$y,
+    mupred = res$mupred_m1,
+    ylp = res$ylp_m1,
+    measure = c("r2", "mse")
+  )
+  pm2 <- loo_pred_measure(
+    loo = res$loo_p_m2,
+    y = res$y,
+    mupred = res$mupred_m2,
+    ylp = res$ylp_m2,
+    measure = c("r2", "mse")
+  )
+
+  expect_snapshot(comp <- loo_compare(pm1, pm2))
+  expect_equal(attr(comp, "sign_converted_measures"), "mse")
+
+  pm_elpd <- loo_pred_measure(
+    loo = res$loo_p_m1,
+    y = res$y,
+    mupred = res$mupred_m1,
+    ylp = res$ylp_m1
+  )
+  expect_no_message(loo_compare(pm_elpd, pm_elpd))
+})
+
+test_that("loo_compare rank_by changes order for loo_pred_measure", {
+  res <- readRDS("data-for-tests/test_data_roaches_compare.Rds")
+  pm1 <- loo_pred_measure(
+    loo = res$loo_p_m1,
+    y = res$y,
+    mupred = res$mupred_m1,
+    ylp = res$ylp_m1,
+    measure = c("r2", "mae")
+  )
+  pm2 <- loo_pred_measure(
+    loo = res$loo_p_m2,
+    y = res$y,
+    mupred = res$mupred_m2,
+    ylp = res$ylp_m2,
+    measure = c("r2", "mae")
+  )
+
+  comp_elpd <- loo_compare(pm1, pm2, rank_by = "elpd")
+  comp_mse <- loo_compare(pm1, pm2, rank_by = "mae")
+  expect_equal(attr(comp_elpd, "rank_by"), "elpd")
+  expect_equal(attr(comp_mse, "rank_by"), "mae")
+  expect_equal(comp_elpd$elpd_diff[1L], 0)
+  expect_equal(comp_mse$mae_diff[1L], 0)
+})
+
+test_that("print.compare.loo works for loo_pred_measure comparisons", {
+  res <- readRDS("data-for-tests/test_data_roaches_compare.Rds")
+  pm1 <- loo_pred_measure(
+    loo = res$loo_p_m1,
+    y = res$y,
+    mupred = res$mupred_m1,
+    ylp = res$ylp_m1,
+    measure = c("r2", "mae")
+  )
+  pm2 <- loo_pred_measure(
+    loo = res$loo_p_m2,
+    y = res$y,
+    mupred = res$mupred_m2,
+    ylp = res$ylp_m2,
+    measure = c("r2", "mae")
+  )
+  pm3 <- loo_pred_measure(
+    loo = res$loo_p_m3,
+    y = res$y,
+    mupred = res$mupred_m3,
+    ylp = res$ylp_m3,
+    measure = c("r2", "mae")
+  )
+
+  comp <- suppressMessages(loo_compare(list(m1 = pm1, m2 = pm2, m3 = pm3)))
+  expect_snapshot(print(comp))
+  expect_snapshot(print(comp, measures = "all", digits = 2))
+  expect_snapshot(print(comp, measures = c("r2", "mae")))
+
+  comp_mae <- suppressMessages(loo_compare(list(m1 = pm1, m2 = pm2), rank_by = "mae"))
+  expect_snapshot(print(comp_mae))
+
+  expect_error(
+    print(comp, measures = "foo"),
+    "Unknown measure\\(s\\) in `measures`"
+  )
+})
+
+test_that("loo_compare measure helpers work as expected", {
+  res <- readRDS("data-for-tests/test_data_roaches_compare.Rds")
+  pm1 <- loo_pred_measure(
+    loo = res$loo_p_m1,
+    y = res$y,
+    mupred = res$mupred_m1,
+    ylp = res$ylp_m1,
+    measure = c("r2", "mse")
+  )
+  pm2 <- loo_pred_measure(
+    loo = res$loo_p_m2,
+    y = res$y,
+    mupred = res$mupred_m2,
+    ylp = res$ylp_m2,
+    measure = c("r2", "mse")
+  )
+  loos <- list(pm1, pm2)
+  cols <- loo:::.compare_pointwise_cols(loos)
+
+  expect_equal(cols, c("elpd_loo", "r2_loo", "mse_loo"))
+  expect_equal(loo:::.compare_measures(loos), c("elpd", "r2", "mse"))
+  expect_equal(loo:::.pointwise_col("mse", cols), "mse_loo")
+  expect_equal(loo:::.pointwise_col("elpd", cols), "elpd_loo")
+  expect_equal(loo:::.display_name("rmse_loo"), "rmse")
+  expect_equal(loo:::.resolve_rank_measure(loos, NULL)$bare, "elpd")
+  expect_equal(loo:::.resolve_rank_measure(loos, "mse")$internal, "mse_loo")
+  expect_true(loo:::.is_elpd_measure("elpd_loo"))
+  expect_false(loo:::.is_elpd_measure("mse_loo"))
+  expect_equal(attr(pm1, "measure_higher_is_better")$mse, NULL)
+  expect_equal(attr(pm1, "measure_higher_is_better")$r2, NULL)
+  expect_equal(attr(pm1, "measure_higher_is_better")$elpd, NULL)
+  expect_equal(attr(pm1, "measure_compare_meta")$elpd$diff_method, "sum")
+  expect_equal(attr(pm1, "measure_compare_meta")$mse$loss, TRUE)
+  expect_equal(attr(pm1, "measure_compare_meta")$mse$diff_method, "mean")
+  expect_equal(attr(pm1, "measure_compare_meta")$r2$diff_method, "pairwise")
+  expect_equal(attr(pm1, "measure_compare_meta")$r2$se_diff_fun, "r2")
+  expect_equal(
+    attr(pm1, "measure_compare_meta")$r2$extra$mse_y_i,
+    (res$y - mean(res$y))^2
+  )
+  # only measures that need it carry `extra`
+  expect_null(attr(pm1, "measure_compare_meta")$mse$extra)
+  expect_true(loo:::.measure_lower_is_better("mse_loo", loos))
+  expect_false(loo:::.measure_lower_is_better("r2_loo", loos))
+  expect_true(loo:::.measure_lower_is_better("mse_loo"))
+  expect_false(loo:::.measure_lower_is_better("r2_loo"))
+  expect_equal(
+    loo:::.compare_sign_converted_measures(c("elpd_loo", "mse_loo", "r2_loo"), loos),
+    c("mse")
+  )
+
+  pair_stats_elpd <- loo:::.pair_measure_stats(
+    pm2, pm1, "elpd_loo", "sum", loos = loos
+  )
+  expect_equal(unname(pair_stats_elpd["se"]), loo:::se_elpd_diff(
+    pm2$pointwise[, "elpd_loo"] - pm1$pointwise[, "elpd_loo"]
+  ))
+  expect_equal(
+    unname(loo:::.pair_measure_stats(pm1, pm1, "elpd_loo", "sum", loos = loos)["diff"]),
+    0
+  )
+
+  pair_mse <- loo:::.pair_measure_stats(
+    pm2, pm1, "mse_loo", "mean", loos = loos
+  )
+  expect_equal(
+    unname(pair_mse["diff"]),
+    pm1$estimates["mse_loo", "Estimate"] - pm2$estimates["mse_loo", "Estimate"]
+  )
+  pair_r2 <- loo:::.pair_measure_stats(
+    pm2, pm1, "r2_loo", "pairwise", loos = loos
+  )
+  expect_equal(
+    unname(pair_r2["diff"]),
+    pm2$estimates["r2_loo", "Estimate"] - pm1$estimates["r2_loo", "Estimate"]
+  )
+  expect_equal(loo:::.measure_pointwise_diff_method(loos, "mse_loo"), "mean")
+  expect_equal(loo:::.measure_pointwise_diff_method(loos, "r2_loo"), "pairwise")
+  expect_equal(loo:::.measure_pointwise_diff_method(loos, "elpd_loo"), "sum")
+  expect_equal(
+    unname(pair_mse["se"]),
+    stats::sd(
+      pm2$pointwise[, "mse_loo"] - pm1$pointwise[, "mse_loo"]
+    ) / sqrt(nrow(pm1$pointwise))
+  )
+})
+
+test_that("rmse differences use the delta-method standard error", {
+  res <- readRDS("data-for-tests/test_data_roaches_compare.Rds")
+  pm1 <- loo_pred_measure(
+    loo = res$loo_p_m1,
+    y = res$y,
+    mupred = res$mupred_m1,
+    ylp = res$ylp_m1,
+    measure = "rmse"
+  )
+  pm2 <- loo_pred_measure(
+    loo = res$loo_p_m2,
+    y = res$y,
+    mupred = res$mupred_m2,
+    ylp = res$ylp_m2,
+    measure = "rmse"
+  )
+  loos <- list(pm1, pm2)
+
+  expect_equal(loo:::.measure_pointwise_diff_method(loos, "rmse_loo"), "pairwise")
+
+  pair <- loo:::.pair_measure_stats(pm2, pm1, "rmse_loo", loos = loos)
+
+  # rmse is a loss, so the reported difference is on the utility scale
+  expect_equal(
+    unname(pair["diff"]),
+    pm1$estimates["rmse_loo", "Estimate"] - pm2$estimates["rmse_loo", "Estimate"]
+  )
+
+  # first-order bivariate Taylor approximation propagated from the MSE scale,
+  # using the covariance between the two models' pointwise squared errors
+  sqe1 <- pm1$pointwise[, "rmse_loo"]
+  sqe2 <- pm2$pointwise[, "rmse_loo"]
+  n <- length(sqe1)
+  mse1 <- mean(sqe1)
+  mse2 <- mean(sqe2)
+  cov_mse <- sum((sqe2 - mse2) * (sqe1 - mse1)) / (n * (n - 1))
+  expected_se <- 0.5 * sqrt(
+    (var(sqe2) / n) / mse2 +
+      (var(sqe1) / n) / mse1 -
+      2 * cov_mse / sqrt(mse2 * mse1)
+  )
+  expect_equal(unname(pair["se"]), expected_se)
+
+  # the standard error is a proper paired quantity, not a sum of the two
+  # per-model standard errors
+  expect_lt(
+    unname(pair["se"]),
+    pm1$estimates["rmse_loo", "SE"] + pm2$estimates["rmse_loo", "SE"]
+  )
+
+  # a model compared against itself has zero difference and zero uncertainty
+  self <- loo:::.pair_measure_stats(pm1, pm1, "rmse_loo", loos = loos)
+  expect_equal(unname(self["diff"]), 0)
+  expect_equal(unname(self["se"]), 0)
+
+  comp <- suppressMessages(loo_compare(pm1, pm2))
+  expect_false(anyNA(comp$rmse_se_diff))
+  expect_equal(comp$rmse_se_diff[1], 0)
+})
+
+test_that("r2 differences use the delta-method standard error", {
+  res <- readRDS("data-for-tests/test_data_roaches_compare.Rds")
+  make <- function(loo, mupred, ylp) {
+    loo_pred_measure(
+      loo = loo, y = res$y, mupred = mupred, ylp = ylp, measure = "r2"
+    )
+  }
+  pm1 <- make(res$loo_p_m1, res$mupred_m1, res$ylp_m1)
+  pm2 <- make(res$loo_p_m2, res$mupred_m2, res$ylp_m2)
+  loos <- list(pm1, pm2)
+
+  expect_equal(loo:::.measure_pointwise_diff_method(loos, "r2_loo"), "pairwise")
+
+  pair <- loo:::.pair_measure_stats(pm2, pm1, "r2_loo", loos = loos)
+
+  # r2 is already a utility, so the difference is reported as stored
+  expect_equal(
+    unname(pair["diff"]),
+    pm2$estimates["r2_loo", "Estimate"] - pm1$estimates["r2_loo", "Estimate"]
+  )
+
+  # first-order trivariate Taylor approximation, written out term by term as in
+  # the derivation rather than in the collapsed single-variance form the
+  # implementation uses
+  sqe1 <- pm1$pointwise[, "r2_loo"]
+  sqe2 <- pm2$pointwise[, "r2_loo"]
+  d <- sqe2 - sqe1
+  n <- length(d)
+  mse_diff <- mean(d)
+  msey_i <- (res$y - mean(res$y))^2
+  mse_y <- mean(msey_i)
+  t1 <- var(d) / n
+  t2 <- -2 * (mse_diff / mse_y) *
+    (sum((d - mse_diff) * (msey_i - mse_y)) / (n * (n - 1)))
+  t3 <- (mse_diff^2 / mse_y^2) * (var(msey_i) / n)
+  expect_equal(unname(pair["se"]), sqrt(t1 + t2 + t3) / mse_y)
+
+  # the difference is also the negative MSE difference over the baseline
+  expect_equal(unname(pair["diff"]), -mse_diff / mse_y)
+
+  # the uncertainty in a difference does not depend on which model is the
+  # reference, even though the difference itself changes sign
+  swapped <- loo:::.pair_measure_stats(pm1, pm2, "r2_loo", loos = loos)
+  expect_equal(unname(swapped["se"]), unname(pair["se"]))
+  expect_equal(unname(swapped["diff"]), -unname(pair["diff"]))
+
+  # a model compared against itself has zero difference and zero uncertainty
+  self <- loo:::.pair_measure_stats(pm1, pm1, "r2_loo", loos = loos)
+  expect_equal(unname(self["diff"]), 0)
+  expect_equal(unname(self["se"]), 0)
+
+  comp <- suppressMessages(loo_compare(pm1, pm2))
+  expect_false(anyNA(comp$r2_se_diff))
+  expect_equal(comp$r2_se_diff[1], 0)
+
+  # the single-model standard error is the same expansion evaluated at one
+  # model's squared errors
+  t1 <- var(sqe1) / n
+  t2 <- -2 * (mean(sqe1) / mse_y) *
+    (sum((sqe1 - mean(sqe1)) * (msey_i - mse_y)) / (n * (n - 1)))
+  t3 <- (mean(sqe1)^2 / mse_y^2) * (var(msey_i) / n)
+  expect_equal(
+    unname(pm1$estimates["r2_loo", "SE"]),
+    sqrt(t1 + t2 + t3) / mse_y
+  )
+})
+
+test_that("r2 se_diff is invariant to `higher_is_better`", {
+  res <- readRDS("data-for-tests/test_data_roaches_compare.Rds")
+  make <- function(loo, mupred, ylp, higher_is_better) {
+    loo_pred_measure(
+      loo = loo,
+      y = res$y,
+      mupred = mupred,
+      ylp = ylp,
+      measure = "r2",
+      control = list(r2 = list(higher_is_better = higher_is_better))
+    )
+  }
+
+  utility <- list(
+    make(res$loo_p_m1, res$mupred_m1, res$ylp_m1, TRUE),
+    make(res$loo_p_m2, res$mupred_m2, res$ylp_m2, TRUE)
+  )
+  loss <- list(
+    make(res$loo_p_m1, res$mupred_m1, res$ylp_m1, FALSE),
+    make(res$loo_p_m2, res$mupred_m2, res$ylp_m2, FALSE)
+  )
+
+  # r2 is naturally a utility, so `higher_is_better = FALSE` is the flipped one
+  expect_equal(loo:::.measure_natural_sign("r2_loo", utility), 1)
+  expect_equal(loo:::.measure_natural_sign("r2_loo", loss), -1)
+
+  # the baseline is a property of `y`, so it is never sign-flipped
+  expect_equal(
+    attr(loss[[1L]], "measure_compare_meta")$r2$extra$mse_y_i,
+    attr(utility[[1L]], "measure_compare_meta")$r2$extra$mse_y_i
+  )
+
+  pair_utility <- loo:::.pair_measure_stats(
+    utility[[2L]], utility[[1L]], "r2_loo", loos = utility
+  )
+  pair_loss <- loo:::.pair_measure_stats(
+    loss[[2L]], loss[[1L]], "r2_loo", loos = loss
+  )
+  expect_equal(unname(pair_loss["se"]), unname(pair_utility["se"]))
+  expect_equal(unname(pair_loss["diff"]), unname(pair_utility["diff"]))
+})
+
+test_that("r2 reports the difference without an se when the baseline is gone", {
+  res <- readRDS("data-for-tests/test_data_roaches_compare.Rds")
+  make <- function(loo, mupred, ylp) {
+    loo_pred_measure(
+      loo = loo, y = res$y, mupred = mupred, ylp = ylp, measure = "r2"
+    )
+  }
+  pm1 <- make(res$loo_p_m1, res$mupred_m1, res$ylp_m1)
+  pm2 <- make(res$loo_p_m2, res$mupred_m2, res$ylp_m2)
+
+  # objects computed before the baseline was stored
+  drop_baseline <- function(x) {
+    meta <- attr(x, "measure_compare_meta")
+    meta$r2$extra <- NULL
+    attr(x, "measure_compare_meta") <- meta
+    x
+  }
+  stale1 <- drop_baseline(pm1)
+  stale2 <- drop_baseline(pm2)
+
+  pair <- loo:::.pair_measure_stats(
+    stale2, stale1, "r2_loo", loos = list(stale1, stale2)
+  )
+  expect_equal(
+    unname(pair["diff"]),
+    pm2$estimates["r2_loo", "Estimate"] - pm1$estimates["r2_loo", "Estimate"]
+  )
+  expect_true(is.na(pair["se"]))
+
+  comp <- suppressMessages(loo_compare(stale1, stale2))
+  expect_false(anyNA(comp$r2_diff))
+  expect_true(all(is.na(comp$r2_se_diff)))
+
+  # one stale model does not cost the others their standard error: the
+  # baseline is shared, so the other model's copy is used, and the metadata
+  # check ignores `extra` rather than reporting it as a disagreement
+  mixed <- suppressMessages(loo_compare(stale1, pm2))
+  expect_false(anyNA(mixed$r2_se_diff))
+})
+
+test_that("rmse se_diff is invariant to `higher_is_better`", {
+  res <- readRDS("data-for-tests/test_data_roaches_compare.Rds")
+  make <- function(loo, mupred, ylp, higher_is_better) {
+    loo_pred_measure(
+      loo = loo,
+      y = res$y,
+      mupred = mupred,
+      ylp = ylp,
+      measure = "rmse",
+      control = list(rmse = list(higher_is_better = higher_is_better))
+    )
+  }
+
+  loss <- list(
+    make(res$loo_p_m1, res$mupred_m1, res$ylp_m1, FALSE),
+    make(res$loo_p_m2, res$mupred_m2, res$ylp_m2, FALSE)
+  )
+  utility <- list(
+    make(res$loo_p_m1, res$mupred_m1, res$ylp_m1, TRUE),
+    make(res$loo_p_m2, res$mupred_m2, res$ylp_m2, TRUE)
+  )
+
+  # stored on opposite scales, so the natural scale must be restored before the
+  # square roots and ratios of the delta method are applied
+  expect_equal(loo:::.measure_natural_sign("rmse_loo", loss), 1)
+  expect_equal(loo:::.measure_natural_sign("rmse_loo", utility), -1)
+
+  pair_loss <- loo:::.pair_measure_stats(
+    loss[[2L]], loss[[1L]], "rmse_loo", loos = loss
+  )
+  pair_utility <- loo:::.pair_measure_stats(
+    utility[[2L]], utility[[1L]], "rmse_loo", loos = utility
+  )
+  expect_equal(unname(pair_loss["se"]), unname(pair_utility["se"]))
+  expect_equal(unname(pair_loss["diff"]), unname(pair_utility["diff"]))
+})
+
+# two balanced-accuracy measures over the same three-class outcome: the second
+# model has probability mass shifted towards the first (and largest) class, so
+# the two disagree on a subset of observations and the class strata are
+# unbalanced, which is where balanced accuracy differs from plain accuracy
+.make_bacc_pms <- function(bias = 0.6, higher_is_better = NULL) {
+  res <- readRDS("data-for-tests/test_data_penguins.Rds")
+  y <- as.integer(res$y)
+  set.seed(4321)
+  ylp <- matrix(
+    rnorm(nrow(res$mupred) * ncol(res$mupred)),
+    nrow = nrow(res$mupred)
+  )
+  biased <- res$mupred
+  biased[, , 1L] <- biased[, , 1L] + bias
+  biased <- sweep(biased, c(1, 2), apply(biased, c(1, 2), sum), "/")
+
+  make <- function(mupred) {
+    suppressWarnings(loo_pred_measure(
+      ylp = ylp,
+      y = y,
+      mupred = mupred,
+      measure = "bacc",
+      control = list(bacc = list(higher_is_better = higher_is_better))
+    ))
+  }
+  list(pm1 = make(res$mupred), pm2 = make(biased), y = y)
+}
+
+test_that("bacc differences use the stratified paired standard error", {
+  fx <- .make_bacc_pms()
+  pm1 <- fx$pm1
+  pm2 <- fx$pm2
+  loos <- list(pm1, pm2)
+
+  expect_equal(loo:::.measure_pointwise_diff_method(loos, "bacc_loo"), "pairwise")
+  expect_equal(
+    attr(pm1, "measure_compare_meta")$bacc$se_diff_fun,
+    "bacc"
+  )
+
+  pair <- loo:::.pair_measure_stats(pm2, pm1, "bacc_loo", loos = loos)
+
+  # bacc is already a utility, so the difference is reported as stored
+  expect_equal(
+    unname(pair["diff"]),
+    pm2$estimates["bacc_loo", "Estimate"] - pm1$estimates["bacc_loo", "Estimate"]
+  )
+
+  # recover the 0/1 accuracies and check the standard error against the
+  # McNemar discordant-count form of the paired difference of proportions,
+  # written out per stratum rather than in the pointwise-variance form the
+  # implementation uses
+  class_id <- attr(pm1, "measure_compare_meta")$bacc$extra$class_id
+  n_c <- tabulate(class_id)
+  K <- length(n_c)
+  acc1 <- round(pm1$pointwise[, "bacc_loo"] * K * n_c[class_id])
+  acc2 <- round(pm2$pointwise[, "bacc_loo"] * K * n_c[class_id])
+  expect_true(all(acc1 %in% c(0, 1)) && all(acc2 %in% c(0, 1)))
+
+  v <- 0
+  for (k in seq_len(K)) {
+    in_k <- class_id == k
+    b <- sum(acc2[in_k] == 1 & acc1[in_k] == 0)
+    cc <- sum(acc2[in_k] == 0 & acc1[in_k] == 1)
+    nk <- n_c[k]
+    # the Wald paired-proportion variance, scaled by nk / (nk - 1) to match the
+    # sample variance the implementation takes
+    v <- v + ((b + cc) / nk^2 - (b - cc)^2 / nk^3) * (nk / (nk - 1))
+  }
+  expect_equal(unname(pair["se"]), sqrt(v) / K)
+
+  # the strata carry information the pointwise vector alone does not: pooling
+  # them would give a different answer, so this is not the `"mean"` path
+  d <- (pm2$pointwise[, "bacc_loo"] - pm1$pointwise[, "bacc_loo"])
+  expect_false(isTRUE(all.equal(
+    unname(pair["se"]),
+    sd(d) / sqrt(length(d))
+  )))
+
+  # the uncertainty in a difference does not depend on which model is the
+  # reference, even though the difference itself changes sign
+  swapped <- loo:::.pair_measure_stats(pm1, pm2, "bacc_loo", loos = loos)
+  expect_equal(unname(swapped["se"]), unname(pair["se"]))
+  expect_equal(unname(swapped["diff"]), -unname(pair["diff"]))
+
+  # a model compared against itself has zero difference and zero uncertainty
+  self <- loo:::.pair_measure_stats(pm1, pm1, "bacc_loo", loos = loos)
+  expect_equal(unname(self["diff"]), 0)
+  expect_equal(unname(self["se"]), 0)
+
+  comp <- suppressMessages(loo_compare(pm1, pm2))
+  expect_false(anyNA(comp$bacc_se_diff))
+  expect_equal(comp$bacc_se_diff[1], 0)
+})
+
+test_that("bacc se_diff is invariant to `higher_is_better`", {
+  utility <- .make_bacc_pms(higher_is_better = TRUE)
+  loss <- .make_bacc_pms(higher_is_better = FALSE)
+
+  # bacc is naturally a utility, so `higher_is_better = FALSE` is the flipped one
+  expect_equal(
+    loo:::.measure_natural_sign("bacc_loo", list(utility$pm1, utility$pm2)), 1
+  )
+  expect_equal(
+    loo:::.measure_natural_sign("bacc_loo", list(loss$pm1, loss$pm2)), -1
+  )
+
+  # the strata are a property of `y`, so they are never sign-flipped
+  expect_equal(
+    attr(loss$pm1, "measure_compare_meta")$bacc$extra$class_id,
+    attr(utility$pm1, "measure_compare_meta")$bacc$extra$class_id
+  )
+
+  pair_utility <- loo:::.pair_measure_stats(
+    utility$pm2, utility$pm1, "bacc_loo", loos = list(utility$pm1, utility$pm2)
+  )
+  pair_loss <- loo:::.pair_measure_stats(
+    loss$pm2, loss$pm1, "bacc_loo", loos = list(loss$pm1, loss$pm2)
+  )
+  expect_equal(unname(pair_loss["se"]), unname(pair_utility["se"]))
+  expect_equal(unname(pair_loss["diff"]), unname(pair_utility["diff"]))
+})
+
+test_that("bacc reports the difference without an se when the strata are gone", {
+  fx <- .make_bacc_pms()
+
+  # objects computed before the class strata were stored
+  drop_strata <- function(x) {
+    meta <- attr(x, "measure_compare_meta")
+    meta$bacc$extra <- NULL
+    attr(x, "measure_compare_meta") <- meta
+    x
+  }
+  stale1 <- drop_strata(fx$pm1)
+  stale2 <- drop_strata(fx$pm2)
+
+  pair <- loo:::.pair_measure_stats(
+    stale2, stale1, "bacc_loo", loos = list(stale1, stale2)
+  )
+  expect_equal(
+    unname(pair["diff"]),
+    fx$pm2$estimates["bacc_loo", "Estimate"] -
+      fx$pm1$estimates["bacc_loo", "Estimate"]
+  )
+  expect_true(is.na(pair["se"]))
+
+  # the strata are shared, so a stale model paired with a current one still
+  # gets a standard error, from whichever copy survives. `pm1` is the better
+  # model and so heads the table; staleness in `pm2` costs nothing at all
+  mixed <- suppressMessages(loo_compare(fx$pm1, stale2))
+  expect_false(anyNA(mixed$bacc_se_diff))
+
+  # but a stale model at the head of the table has no second copy to fall back
+  # on for its own row, which is a comparison against itself
+  mixed_stale_first <- suppressMessages(loo_compare(stale1, fx$pm2))
+  expect_true(is.na(mixed_stale_first$bacc_se_diff[1L]))
+  expect_false(is.na(mixed_stale_first$bacc_se_diff[2L]))
+})
+
+test_that("custom measures can supply their own se_diff_fun", {
+  res <- readRDS("data-for-tests/test_data_roaches_compare.Rds")
+
+  # a custom rmse, whose estimate is neither a sum nor a mean of `pointwise`
+  my_rmse <- function(y, mupred) {
+    sqe_i <- (y - colMeans(mupred))^2
+    list(
+      estimate = sqrt(mean(sqe_i)),
+      se = sqrt(var(sqe_i) / length(sqe_i)) / (2 * sqrt(mean(sqe_i))),
+      pointwise = sqe_i
+    )
+  }
+  attr(my_rmse, "measure_name") <- "my_rmse"
+
+  make <- function(loo, mupred, ylp, fun) {
+    loo_pred_measure(
+      loo = loo, y = res$y, mupred = mupred, ylp = ylp, measure = fun
+    )
+  }
+
+  # without an se_diff_fun the difference has no standard error
+  plain <- list(
+    make(res$loo_p_m1, res$mupred_m1, res$ylp_m1, my_rmse),
+    make(res$loo_p_m2, res$mupred_m2, res$ylp_m2, my_rmse)
+  )
+  expect_equal(
+    loo:::.measure_pointwise_diff_method(plain, "my_rmse_loo"),
+    "estimates_only"
+  )
+  expect_true(is.na(
+    loo:::.pair_measure_stats(plain[[2L]], plain[[1L]], "my_rmse_loo", loos = plain)["se"]
+  ))
+
+  attr(my_rmse, "se_diff_fun") <- loo:::.se_diff_rmse
+  with_se <- list(
+    make(res$loo_p_m1, res$mupred_m1, res$ylp_m1, my_rmse),
+    make(res$loo_p_m2, res$mupred_m2, res$ylp_m2, my_rmse)
+  )
+  expect_equal(
+    attr(with_se[[1L]], "measure_compare_meta")$my_rmse$diff_method,
+    "pairwise"
+  )
+  expect_equal(
+    loo:::.measure_pointwise_diff_method(with_se, "my_rmse_loo"),
+    "pairwise"
+  )
+
+  pair <- loo:::.pair_measure_stats(
+    with_se[[2L]], with_se[[1L]], "my_rmse_loo", loos = with_se
+  )
+  expect_false(is.na(pair["se"]))
+  expect_gt(unname(pair["se"]), 0)
+
+  # a custom measure can carry its own auxiliary data through to its
+  # se_diff_fun, and each model receives its own copy
+  my_scaled <- function(y, mupred) {
+    sqe_i <- (y - colMeans(mupred))^2
+    list(
+      estimate = sqrt(mean(sqe_i)),
+      se = sqrt(var(sqe_i) / length(sqe_i)) / (2 * sqrt(mean(sqe_i))),
+      pointwise = sqe_i,
+      extra = list(scale = stats::sd(y), n_used = length(y))
+    )
+  }
+  attr(my_scaled, "measure_name") <- "my_scaled"
+  attr(my_scaled, "se_diff_fun") <- function(ref, cmp) {
+    stopifnot(
+      identical(ref$extra$n_used, length(ref$pointwise)),
+      identical(cmp$extra$scale, ref$extra$scale)
+    )
+    loo:::.se_diff_rmse(ref, cmp) / ref$extra$scale
+  }
+  scaled <- list(
+    make(res$loo_p_m1, res$mupred_m1, res$ylp_m1, my_scaled),
+    make(res$loo_p_m2, res$mupred_m2, res$ylp_m2, my_scaled)
+  )
+  expect_equal(
+    attr(scaled[[1L]], "measure_compare_meta")$my_scaled$extra,
+    list(scale = stats::sd(res$y), n_used = length(res$y))
+  )
+  pair_scaled <- loo:::.pair_measure_stats(
+    scaled[[2L]], scaled[[1L]], "my_scaled_loo", loos = scaled
+  )
+  expect_equal(
+    unname(pair_scaled["se"]),
+    unname(pair["se"]) / stats::sd(res$y)
+  )
+
+  # `extra` that is not a list is rejected at compute time
+  my_bad_extra <- function(y, mupred) {
+    sqe_i <- (y - colMeans(mupred))^2
+    list(
+      estimate = mean(sqe_i),
+      se = sqrt(var(sqe_i) / length(sqe_i)),
+      pointwise = sqe_i,
+      extra = 1
+    )
+  }
+  attr(my_bad_extra, "measure_name") <- "my_bad_extra"
+  expect_error(
+    make(res$loo_p_m1, res$mupred_m1, res$ylp_m1, my_bad_extra),
+    "must be a list"
+  )
+
+  # a custom measure declaring an se_diff_fun that returns nonsense is caught
+  attr(my_rmse, "se_diff_fun") <- function(ref, cmp) c(1, 2)
+  bad <- list(
+    make(res$loo_p_m1, res$mupred_m1, res$ylp_m1, my_rmse),
+    make(res$loo_p_m2, res$mupred_m2, res$ylp_m2, my_rmse)
+  )
+  expect_error(
+    loo:::.pair_measure_stats(bad[[2L]], bad[[1L]], "my_rmse_loo", loos = bad),
+    "must return a numeric scalar"
+  )
+})
+
+test_that("loo_compare errors on inconsistent measure metadata", {
+  res <- readRDS("data-for-tests/test_data_roaches_compare.Rds")
+  pm1 <- loo_pred_measure(
+    loo = res$loo_p_m1,
+    y = res$y,
+    mupred = res$mupred_m1,
+    ylp = res$ylp_m1,
+    measure = "mse",
+    control = list(mse = list(higher_is_better = NULL))
+  )
+  pm2 <- loo_pred_measure(
+    loo = res$loo_p_m2,
+    y = res$y,
+    mupred = res$mupred_m2,
+    ylp = res$ylp_m2,
+    measure = "mse",
+    control = list(mse = list(higher_is_better = TRUE))
+  )
+
+  expect_error(
+    suppressMessages(loo_compare(pm1, pm2)),
+    "disagree on comparison metadata for measure 'mse'"
+  )
+})
+
+test_that("loo_compare errors when compare metadata is missing on some models", {
+  res <- readRDS("data-for-tests/test_data_roaches_compare.Rds")
+  pm1 <- loo_pred_measure(
+    loo = res$loo_p_m1,
+    y = res$y,
+    mupred = res$mupred_m1,
+    ylp = res$ylp_m1,
+    measure = "mse"
+  )
+  pm2 <- loo_pred_measure(
+    loo = res$loo_p_m2,
+    y = res$y,
+    mupred = res$mupred_m2,
+    ylp = res$ylp_m2,
+    measure = "mse"
+  )
+  compare_meta <- attr(pm2, "measure_compare_meta")
+  compare_meta$mse <- NULL
+  attr(pm2, "measure_compare_meta") <- compare_meta
+
+  expect_error(
+    suppressMessages(loo_compare(pm1, pm2)),
+    "Not all models provide comparison metadata for measure 'mse'"
+  )
+})
+
+test_that("loo_compare warns when rank_by is ignored for classic loo objects", {
+  expect_warning(
+    loo_compare(w1, w2, rank_by = "mse"),
+    "`rank_by` is only used for `loo_pred_measure` comparisons"
+  )
+})
+
+.make_compare_pm <- function(res, model = 1L, measure, extra_args = list()) {
+  suffix <- model
+  args <- c(
+    list(
+      loo = res[[paste0("loo_p_m", suffix)]],
+      y = res$y,
+      mupred = res[[paste0("mupred_m", suffix)]],
+      ylp = res[[paste0("ylp_m", suffix)]],
+      measure = measure
+    ),
+    extra_args
+  )
+  do.call(loo_pred_measure, args)
+}
+
+.make_compare_pm_synthetic <- function(measure) {
+  if (measure == "brier") {
+    res_binary <- readRDS("data-for-tests/test_data_binary.Rds")
+    ylp <- matrix(
+      rnorm(nrow(res_binary$ypred) * ncol(res_binary$ypred)),
+      nrow = nrow(res_binary$ypred)
+    )
+    return(loo_pred_measure(
+      ylp = ylp,
+      y = res_binary$y,
+      ypred = res_binary$ypred,
+      measure = measure
+    ))
+  }
+  if (measure %in% c("acc", "bacc")) {
+    res_cat <- readRDS("data-for-tests/test_data_penguins.Rds")
+    ylp <- matrix(
+      rnorm(nrow(res_cat$mupred) * ncol(res_cat$mupred)),
+      nrow = nrow(res_cat$mupred)
+    )
+    return(loo_pred_measure(
+      ylp = ylp,
+      y = as.integer(res_cat$y),
+      mupred = res_cat$mupred,
+      measure = measure
+    ))
+  }
+  stop("Unsupported synthetic measure: ", measure)
+}
+
+test_that("loo_compare works for all built-in measures", {
+  res <- readRDS("data-for-tests/test_data_roaches_compare.Rds")
+  res_roaches <- readRDS("data-for-tests/test_data_roaches.Rds")
+  roaches_measures <- c(
+    "ic", "mlpd", "mae", "r2", "rmse", "mse"
+  )
+  for (measure in roaches_measures) {
+    pm1 <- .make_compare_pm(res, 1L, measure)
+    pm2 <- .make_compare_pm(res, 2L, measure)
+    comp <- suppressMessages(loo_compare(pm1, pm2))
+    expect_true(paste0(measure, "_diff") %in% colnames(comp), info = measure)
+    expect_equal(attr(comp, "compare_measures"), c("elpd", measure), info = measure)
+  }
+
+  for (measure in c("rps", "srps")) {
+    pm1 <- loo_pred_measure(
+      loo = res$loo_p_m1,
+      y = res$y,
+      ypred = res_roaches$ypred,
+      ylp = res$ylp_m1,
+      measure = measure
+    )
+    pm2 <- loo_pred_measure(
+      loo = res$loo_p_m2,
+      y = res$y,
+      ypred = res_roaches$ypred,
+      ylp = res$ylp_m2,
+      measure = measure
+    )
+    comp <- suppressMessages(loo_compare(pm1, pm2))
+    expect_true(paste0(measure, "_diff") %in% colnames(comp), info = measure)
+    expect_equal(attr(comp, "compare_measures"), c("elpd", measure), info = measure)
+  }
+
+  for (measure in c("brier", "acc", "bacc")) {
+    pm1 <- .make_compare_pm_synthetic(measure)
+    pm2 <- .make_compare_pm_synthetic(measure)
+    comp <- suppressMessages(loo_compare(pm1, pm2))
+    expect_true(paste0(measure, "_diff") %in% colnames(comp), info = measure)
+    expect_equal(attr(comp, "compare_measures"), c("elpd", measure), info = measure)
+  }
+})
+
+.make_many_compare_pms <- function(res, n, noise_scale = 0.01) {
+  lapply(seq_len(n), function(i) {
+    loo_pred_measure(
+      loo = res$loo_p_m1,
+      y = res$y,
+      mupred = res$mupred_m1 + rnorm(length(res$y), 0, noise_scale * i),
+      ylp = res$ylp_m1,
+      measure = "mae"
+    )
+  })
+}
+
+test_that("loo_compare warns for many loo_pred_measure models", {
+  res <- readRDS("data-for-tests/test_data_roaches_compare.Rds")
+  set.seed(123)
+  pm_list <- .make_many_compare_pms(res, 25L)
+  expect_warning(
+    suppressMessages(loo_compare(pm_list)),
+    "Difference in performance potentially due to chance. See McLatchie and Vehtari (2023) for details.",
+    fixed = TRUE
+  )
+
+  pm_list_short <- .make_many_compare_pms(res, 4L)
+  expect_no_warning(suppressMessages(loo_compare(pm_list_short)))
+})
+
 test_that("loo_compare throws appropriate warnings", {
   w3 <- w1
   w4 <- w2
diff --git a/vignettes/articles-online-only/comparison-standard-errors.Rmd b/vignettes/articles-online-only/comparison-standard-errors.Rmd
new file mode 100644
index 00000000..3a84d8ae
--- /dev/null
+++ b/vignettes/articles-online-only/comparison-standard-errors.Rmd
@@ -0,0 +1,591 @@
+---
+title: "Differences and their standard errors in model comparison"
+author: "Florence Bockting"
+date: "`r Sys.Date()`"
+output:
+  rmarkdown::html_vignette:
+    toc: true
+    toc_depth: 3
+params:
+  EVAL: TRUE #!r identical(Sys.getenv("NOT_CRAN"), "true")
+vignette: >
+  %\VignetteIndexEntry{comparison-standard-errors}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r, child="../children/SEE-ONLINE.txt", eval = if (isTRUE(exists("params"))) !params$EVAL else TRUE}
+```
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>",
+  message = FALSE,
+  warning = FALSE,
+  eval = params$EVAL,
+  fig.width = 7,
+  fig.height = 4,
+  fig.align = "center"
+)
+devtools::load_all(".")
+```
+
+## Introduction
+
+When you call `loo_compare()` on results from `loo_pred_measure()`, every
+measure gets two columns: a difference and its standard error.
+For example, `loo_compare` based on `loo_pred_measure` with measures `elpd`,
+`mse`, and `rmse`, will produce:
+
+```
+model  elpd_diff  se_diff   mse_diff  mse_se_diff   rmse_diff  rmse_se_diff
+```
+
+The difference (`elpd_diff`, `mse_diff`, `rmse_diff`) tells you *which* model
+did better and the standard error provides an estimate in how far you can trust
+this estimate. In this vignette we explain:
+
+1. The approaches for computing the point estimate and standard error of the
+   difference for the different measures.
+2. How to provide a formula for computing the standard error of a difference for
+   a custom measure.
+
+## A simple example
+
+In the following, we want to compare three Poisson regression models on the same data set.
+For this we use the `roaches` data from **rstanarm**. It records a pest-management trial in
+262 apartment buildings. The outcome `y` is the number of roaches caught in
+traps after treatment and ranges from 0 to 357. The predictors are:
+
+| Predictor | Description |
+| :-- | :---- |
+| `roach1` | pre-treatment roach count (used as `sqrt_roach1 = sqrt(roach1)`) |
+| `treatment` | whether the building received pest management |
+| `senior` | whether the building is restricted to elderly residents |
+| `exposure2` | number of days the traps were exposed |
+
+The three models that we will use in the following comparison tasks are Poisson
+regressions with different predictors:
+
++ `m1`: `y ~ treatment + senior`
++ `m2`: `y ~ sqrt_roach1 + senior`
++ `m3`: `y ~ sqrt_roach1 + treatment`
+
+For each model we compute its predictive performance via `loo_pred_measure()` with
+measures `elpd`, `mae`, `mse`, `rmse`, and `r2`.
+
+```{r load-data}
+res <- readRDS("../../tests/testthat/data-for-tests/test_data_roaches_compare.Rds")
+
+fit_measure <- function(i, measure) {
+  loo_pred_measure(
+    y       = res$y,
+    mupred  = res[[paste0("mupred_m", i)]],
+    ylp     = res[[paste0("ylp_m", i)]],
+    measure = measure
+  )
+}
+
+m1 <- fit_measure(1, c("mae", "mse", "rmse", "r2"))
+m2 <- fit_measure(2, c("mae", "mse", "rmse", "r2"))
+m3 <- fit_measure(3, c("mae", "mse", "rmse", "r2"))
+
+comp <- loo_compare(list(m1 = m1, m2 = m2, m3 = m3))
+print(comp, measures = "all")
+```
+
+## A closer look at the model comparison output
+
+### Estimated difference between two models
+The estimates `elpd_diff`, `mae_diff`, `rmse_diff`, and `r2_diff` 
+represent the differences of the estimates of the two models which enter the 
+comparison (e.g., `m1$estimates["elpd_loo", "Estimate"] - m3$estimates["elpd_loo", "Estimate"]`). 
+Whereby the first model in this difference is the comparison model and the
+second model is the reference model.
+
+```{r point-estimate}
+est <- function(x, measure) x$estimates[measure, "Estimate"]
+
+data.frame(
+  measure  = c(
+    "elpd_diff", "mse_diff", "mae_diff", "rmse_diff", "r2_diff"
+  ),
+  by_hand  = c(
+    est(m1, "elpd_loo") - est(m3, "elpd_loo"),
+    est(m1, "mse_loo")  - est(m3, "mse_loo"),
+    est(m1, "mae_loo")  - est(m3, "mae_loo"),
+    est(m1, "rmse_loo") - est(m3, "rmse_loo"),
+    est(m1, "r2_loo")   - est(m3, "r2_loo")
+  ),
+  reported = c(
+    comp$elpd_diff[comp$model == "m1"],
+    comp$mse_diff[comp$model == "m1"],
+    comp$mae_diff[comp$model == "m1"],
+    comp$rmse_diff[comp$model == "m1"],
+    comp$r2_diff[comp$model == "m1"]
+  )
+)
+```
+
+Note that the reported `mae_diff`, `mse_diff`, and `rmse_diff` rows differ in 
+sign from the hand-computed value, and this is deliberate. `mae`, `mse`, and `rmse`
+are **losses**, while `elpd` and `r2` are **utilities**. For losses, the stimated value
+is interpreted as lower is better, while for utilities, higher is better.
+
+However, `loo_compare()` reports every difference on a **utility scale** 
+(higher is better). Measures that are on a loss scale are converted to a utility scale.
+In consequnence, a negative number always means "this model is 
+worse than the reference model", whatever the measure. 
+
+`loo_compare()` informs the user when a measure is converted to a utility scale
+for the model comparison task:
+
+```{r sign-message, warning=TRUE, message=TRUE}
+comp <- loo_compare(list(m1 = m1, m2 = m2, m3 = m3))
+```
+
+### Standard error of the difference between two models
+
+Depending on the measure there are different approaches to compute 
+the standard error (SE) of the difference. Which formula a measure uses 
+is recorded on the result object, in the `measure_compare_meta` 
+attribute (see `?"loo-glossary"`):
+
+```{r diff-method}
+meta <- attr(m1, "measure_compare_meta")
+sapply(meta, function(x) x$diff_method)
+```
+
+There are four possible values:
+
+| `diff_method` | SE of the difference | Measures |
+| :--- | :--- | :--- |
+| `"sum"` | `sqrt(N) * sd(d_i)` | `elpd`, `ic`, `mlpd` |
+| `"mean"` | `sd(d_i) / sqrt(N)` | `mae`, `mse`, `acc`, `brier`, `rps`, `srps` |
+| `"pairwise"` | the measure's own `se_diff_fun` | `rmse`, `r2`, `bacc` |
+| `"auto"` | inferred at comparison time | custom measures |
+
+whereby `N` is the number of observations and `d_i = pointwise_cmp - pointwise_ref`
+is the paired difference of the pointwise estimates of the two models 
+(e.g., `m1$pointwise[, "elpd_loo"]` and `m3$pointwise[, "elpd_loo"]`). 
+
+A fifth value, `"estimates_only"`, means no standard error is available. In 
+this case the estimate for the difference is still reported but the standard 
+error gets `NA`. This can only happen when `"auto"` cannot work out how a
+custom measure's pointwise values relate to its estimate (see here[TODO] for more
+details on custom measures).
+
+#### Method: "sum"
+
+For `elpd`, `ic`, and `mlpd` the estimate is a sum of pointwise log predictive
+densities, so the standard error of a difference scales up with `sqrt(N)`:
+
+```{r sum-method}
+n <- attr(m1, "dims")[2]
+d_elpd <- m1$pointwise[, "elpd_loo"] - m3$pointwise[, "elpd_loo"]
+
+c(by_hand     = sqrt(n) * sd(d_elpd),
+  reported = comp$se_diff[comp$model == "m1"])
+```
+
+#### Method: "mean"
+
+For `mae`, `mse`, and most scoring rules, the estimate is a mean of pointwise
+squared errors, so the same standard deviation is divided by `sqrt(N)` instead:
+
+```{r mean-method}
+d_mse <- m1$pointwise[, "mse_loo"] - m3$pointwise[, "mse_loo"]
+
+c(hand     = sd(d_mse) / sqrt(n),
+  reported = comp$mse_se_diff[comp$model == "m1"])
+```
+
+#### Method: "pairwise"
+
+**RMSE**
+
+For RMSE, the delta method [REF-TODO] on the MSE scale is used 
+and propagated through the square root:
+
+$$
+SE_{\Delta RMSE} = \frac{1}{2}\sqrt{
+  \frac{SE_{MSE(a)}^2}{MSE(a)} + \frac{SE_{MSE(b)}^2}{MSE(b)}
+  - \frac{2\,\mathrm{Cov}[MSE(a), MSE(b)]}{RMSE(a)\,RMSE(b)}
+}
+$$
+
+Everything on the right hand side is computable from the stored pointwise
+squared errors:
+
+```{r pairwise-rmse}
+sqe_ref <- m3$pointwise[, "rmse_loo"]
+sqe_cmp <- m1$pointwise[, "rmse_loo"]
+
+mse_ref <- mean(sqe_ref)
+mse_cmp <- mean(sqe_cmp)
+se_ref  <- sqrt(var(sqe_ref) / n)
+se_cmp  <- sqrt(var(sqe_cmp) / n)
+cov_rc  <- sum((sqe_cmp - mse_cmp) * (sqe_ref - mse_ref)) / (n * (n - 1))
+
+se_rmse <- 0.5 * sqrt(
+  se_cmp^2 / mse_cmp + se_ref^2 / mse_ref -
+    2 * cov_rc / (sqrt(mse_cmp) * sqrt(mse_ref))
+)
+
+c(by_hand  = se_rmse,
+  reported = comp$rmse_se_diff[comp$model == "m1"])
+```
+
+**R2**
+
+For R2, the delta method with one more term is applied, because its 
+denominator `MSE(y)` is estimated from the same data and carries its own
+uncertainty. It needs the pointwise baseline `(y_i - mean(y))^2`, which
+`loo_compare()` no longer has access to by the time it runs, so 
+`measure_r2()` stores it on the object:
+
+```{r r2-extra}
+str(attr(m1, "measure_compare_meta")$r2$extra)
+```
+
+```{r pairwise-r2}
+mse_y_i <- attr(m1, "measure_compare_meta")$r2$extra$mse_y_i
+mse_y   <- mean(mse_y_i)
+
+d_sqe  <- m1$pointwise[, "r2_loo"] - m3$pointwise[, "r2_loo"]
+scaled <- d_sqe - (mean(d_sqe) / mse_y) * mse_y_i
+
+c(by_hand  = sd(scaled) / sqrt(n) / mse_y,
+  reported = comp$r2_se_diff[comp$model == "m1"])
+```
+
+**bacc** 
+The balanced accuracy (bacc) is the average of the per-class accuracies, which 
+gives each class equal weight no matter how rare it is. That makes it the 
+standard choice for imbalanced outcomes.
+
+To explain the standard error of the difference for bacc, we need to introduce
+another example where the outcome variable is binary. Therfore, we simulate 200 
+observations from a logistic model with a single continuous predictor `x`. 
+
+```{r bacc-data}
+set.seed(2024)
+n_obs   <- 200
+n_draws <- 400
+
+x <- rnorm(n_obs)
+y <- rbinom(n_obs, 1, plogis(-1.4 + 1.6 * x))
+table(y)
+```
+
+TODO from here
+
+## Custom measures
+
+A custom measure is any function returning a list with `estimate`, `se`, and
+`pointwise`. By default it gets `diff_method = "auto"`, and `loo_compare()`
+works out what to do by checking how the estimate relates to the pointwise
+values.
+
+### When autodetection works
+
+If your estimate really is the mean of your pointwise values, autodetection
+finds that and you get a standard error for free. Here is a Huber-type loss:
+
+```{r custom-auto}
+huber <- function(y, mupred) {
+  delta <- 10
+  r <- y - colMeans(mupred)
+  l <- ifelse(abs(r) <= delta, 0.5 * r^2, delta * (abs(r) - 0.5 * delta))
+  list(estimate = mean(l), se = sd(l) / sqrt(length(l)), pointwise = l)
+}
+attr(huber, "measure_name") <- "huber"
+
+h1 <- fit_measure(1, huber)
+h3 <- fit_measure(3, huber)
+
+attr(h1, "measure_compare_meta")$huber$diff_method
+
+comp_h <- loo_compare(list(m3 = h3, m1 = h1))
+comp_h[, c("model", "huber_diff", "huber_se_diff")]
+```
+
+```{r custom-auto-check}
+d_huber <- h1$pointwise[, "huber_loo"] - h3$pointwise[, "huber_loo"]
+
+c(hand     = sd(d_huber) / sqrt(n),
+  reported = comp_h$huber_se_diff[comp_h$model == "m1"])
+```
+
+Two things to note. First, autodetection confirmed that `mean(pointwise)`
+reproduces the estimate, and used the `"mean"` formula. Second, the reported
+difference is **positive** even though a larger Huber loss is worse. Custom
+measures are always treated as utilities — higher is better — because the
+package has no way to know otherwise. If your measure is a loss, either return
+its negative, or read the difference column with that in mind.
+
+### When autodetection cannot work
+
+Now a measure whose estimate is *not* a sum or mean of its pointwise values —
+RMSE normalised by the standard deviation of the outcome:
+
+```{r custom-estimates-only}
+nrmse <- function(y, mupred) {
+  sqe <- (y - colMeans(mupred))^2
+  list(
+    estimate  = sqrt(mean(sqe)) / sd(y),
+    se        = sqrt(var(sqe) / length(sqe)) / (2 * sqrt(mean(sqe))) / sd(y),
+    pointwise = sqe,
+    extra     = list(sd_y = sd(y))
+  )
+}
+attr(nrmse, "measure_name") <- "nrmse"
+
+n1 <- fit_measure(1, nrmse)
+n3 <- fit_measure(3, nrmse)
+
+attr(n1, "measure_compare_meta")$nrmse$diff_method
+
+comp_n <- loo_compare(list(m3 = n3, m1 = n1))
+comp_n[, c("model", "nrmse_diff", "nrmse_se_diff")]
+```
+
+The stored `diff_method` is still `"auto"`, because `"auto"` is not a formula —
+it is an instruction to work one out at comparison time. Here that attempt
+failed: neither `sum(pointwise)` nor `mean(pointwise)` reproduces the estimate,
+so `loo_compare()` resolved it to `"estimates_only"` and reported the difference
+with an `NA` standard error.
+
+This is not an error. It is the package declining to print a number it cannot
+justify, which is the correct behaviour when the alternative is a confident
+wrong one.
+
+It is worth being aware of what autodetection can and cannot check. It verifies
+that the estimate equals the sum or the mean of the pointwise values, which is
+an *algebraic* fact about the numbers in front of it. The standard error it then
+applies needs something stronger: that those pointwise values behave as `N`
+exchangeable contributions whose sample standard deviation describes how much
+the estimate would vary. Balanced accuracy is a measure where the first holds
+and the second does not — its pointwise values do sum to its estimate, but
+using the `"sum"` formula would throw away the class stratification, which is
+why it registers its own function instead. If your measure has structure of that
+kind, supply an `se_diff_fun` rather than relying on autodetection.
+
+### Supplying your own `se_diff_fun`
+
+To get a real standard error, attach a function to the measure. It is called as
+`se_diff_fun(ref, cmp)`, where each argument is a list describing **one** model
+with elements `estimate`, `se`, `pointwise`, and `extra`, and it returns a
+single number.
+
+Normalised RMSE is RMSE divided by a quantity that does not depend on the
+model, so its delta-method standard error is the RMSE one divided by `sd(y)` —
+which is exactly what the measure stashed in `extra` for this purpose:
+
+```{r custom-se-diff-fun}
+attr(nrmse, "se_diff_fun") <- function(ref, cmp) {
+  n <- length(ref$pointwise)
+  mse_ref <- mean(ref$pointwise)
+  mse_cmp <- mean(cmp$pointwise)
+  se_ref  <- sqrt(var(ref$pointwise) / n)
+  se_cmp  <- sqrt(var(cmp$pointwise) / n)
+  cov_rc  <- sum((cmp$pointwise - mse_cmp) * (ref$pointwise - mse_ref)) /
+    (n * (n - 1))
+
+  rel_ref <- se_ref / sqrt(mse_ref)
+  rel_cmp <- se_cmp / sqrt(mse_cmp)
+  rho     <- cov_rc / (se_ref * se_cmp)
+
+  # algebraically the same as the RMSE formula shown earlier, but written so
+  # that a model compared against itself cancels to exactly zero
+  se_rmse <- 0.5 * sqrt(
+    (rel_cmp - rel_ref)^2 + 2 * rel_cmp * rel_ref * max(1 - rho, 0)
+  )
+
+  se_rmse / ref$extra$sd_y
+}
+
+n1 <- fit_measure(1, nrmse)
+n3 <- fit_measure(3, nrmse)
+
+attr(n1, "measure_compare_meta")$nrmse$diff_method
+
+comp_n <- loo_compare(list(m3 = n3, m1 = n1))
+comp_n[, c("model", "nrmse_diff", "nrmse_se_diff")]
+```
+
+Attaching the function switched `diff_method` from `"auto"` to `"pairwise"`,
+and the standard error column is populated.
+
+Three practical points:
+
+- **Write the formula so that identical inputs cancel exactly.** The literal
+  three-term version of the RMSE delta method can produce a tiny negative
+  number under the square root when both models have the same squared errors,
+  which is the first row of every comparison table. The rearranged form above
+  returns exactly `0` there.
+- **Your function always sees the measure's natural scale.** If you set an
+  orientation that flips the stored values, the package undoes it before
+  calling you and reapplies it to the difference afterwards. You never have to
+  think about sign conventions.
+- **`extra` is how you get data to your function.** Anything your measure
+  returns as `extra` is stored on the object and handed back to
+  `se_diff_fun()`. This is the only way to reach data, like `y`, that is no
+  longer available when `loo_compare()` runs.
+
+One caveat about storage: a custom `se_diff_fun` is stored on the result object
+as a closure, so `saveRDS()` will serialise its enclosing environment. If you
+define it inline in a script where a large draws matrix is in scope, that matrix
+goes into the file. Define such functions in a clean environment if object size
+matters. Built-in measures avoid this by storing a name instead of a function.
+
+## When measures disagree
+
+Everything so far has been about computing the two columns. This section is
+about what they are for.
+
+Different measures do not have to agree on which model is best, and in the
+roaches comparison they do not. Ranking the three models by each measure
+separately:
+
+```{r disagree-ranking}
+estimates <- sapply(
+  list(m1 = m1, m2 = m2, m3 = m3),
+  function(x) x$estimates[, "Estimate"]
+)
+round(estimates, 3)
+```
+
+```{r disagree-best}
+# elpd and r2 are utilities (higher is better); mae, mse and rmse are losses
+higher_is_better <- c(elpd_loo = TRUE, p_loo = NA, mae_loo = FALSE,
+                      mse_loo = FALSE, rmse_loo = FALSE, r2_loo = TRUE)
+
+ranking <- function(measure) {
+  v <- sort(estimates[measure, ], decreasing = higher_is_better[measure])
+  paste(names(v), collapse = " > ")
+}
+
+measures <- c("elpd_loo", "mae_loo", "mse_loo", "rmse_loo", "r2_loo")
+data.frame(measure = measures, ranking = sapply(measures, ranking),
+           row.names = NULL)
+```
+
+Two different answers. `elpd` and `mae` put **m3** first; `mse`, `rmse`, and
+`r2` put **m2** first. A reader who ran only one of these would come away with a
+different "best model".
+
+This is not a malfunction, and the pattern is a common one:
+
+- `mse`, `rmse`, and `r2` agree with each other because they are all monotone
+  functions of the same pointwise squared errors — they cannot disagree, and
+  reporting all three is reporting one thing three times.
+- Squared error is dominated by the largest residuals. `mae` weights every
+  observation's error linearly and so is far less sensitive to a handful of
+  badly predicted buildings, and `elpd` scores the whole predictive
+  distribution rather than a point summary. When those two part company from
+  the squared-error family, it usually means a few extreme observations are
+  driving the ranking.
+
+So which model is best? The standard error column answers that, and the answer
+is that the question is not settled by this data:
+
+```{r disagree-se}
+row_m2 <- comp[comp$model == "m2", ]
+
+data.frame(
+  measure = c("elpd", "mae", "mse", "rmse", "r2"),
+  diff    = c(row_m2$elpd_diff, row_m2$mae_diff, row_m2$mse_diff,
+              row_m2$rmse_diff, row_m2$r2_diff),
+  se_diff = c(row_m2$se_diff, row_m2$mae_se_diff, row_m2$mse_se_diff,
+              row_m2$rmse_se_diff, row_m2$r2_se_diff)
+) |>
+  transform(ratio = round(abs(diff) / se_diff, 2))
+```
+
+Every one of these differences is **smaller than its own standard error**. The
+`elpd` difference is about one standard error; the rest are a fifth of one. Both
+orderings are noise. The measures are not really contradicting each other — they
+are all saying that `m2` and `m3` are indistinguishable on 262 observations, and
+the apparent ranking flip is which side of zero the noise happened to land on.
+
+Note also that all three squared-error measures produce nearly the same ratio.
+For `mse` and `rmse` that is exact up to the delta-method approximation: a
+first-order expansion scales a difference and its standard error by the same
+derivative, so the ratio between them survives the square root. `r2` is very
+slightly different (`0.199` against `0.202` at full precision) because its
+standard error also carries the uncertainty in the estimated baseline `MSE(y)`,
+which is not simply a transformation of `MSE`.
+
+`m1` is the one genuine finding here — it drops `sqrt_roach1`, and every measure
+puts it last by a margin large relative to its standard error. That is a
+conclusion the data supports; "m2 beats m3 on RMSE" is not.
+
+The practical advice:
+
+1. **Choose your measure before you look**, based on what you want the model to
+   do well. Picking the measure that favours your preferred model after the fact
+   is a way of fooling yourself, and with several measures on the table it is an
+   easy mistake to make.
+2. **Read the difference against its own standard error**, not against zero. A
+   ranking is only meaningful if the differences that produce it are.
+3. **Do not read `mse`, `rmse`, and `r2` as three independent votes.** They are
+   one measure in three costumes.
+
+## Summary
+
+```{r summary-table, echo=FALSE}
+knitr::kable(
+  data.frame(
+    Measure = c("elpd, ic, mlpd", "mae, mse, acc, brier, rps, srps",
+                "rmse", "r2", "bacc", "custom (default)",
+                "custom (with se_diff_fun)"),
+    `diff_method` = c("sum", "mean", "pairwise", "pairwise", "pairwise",
+                      "auto", "pairwise"),
+    `Standard error` = c(
+      "sqrt(N) * sd(d_i)",
+      "sd(d_i) / sqrt(N)",
+      "delta method on the MSE scale",
+      "delta method, including baseline uncertainty",
+      "stratified by class, paired within class",
+      "detected as sum or mean, otherwise NA",
+      "whatever your function returns"
+    ),
+    check.names = FALSE
+  )
+)
+```
+
+The things worth remembering:
+
+1. The difference is always `estimate_cmp - estimate_ref`, reported so that
+   negative means worse.
+2. The standard error is always paired, because the models are evaluated on the
+   same observations.
+3. `NA` in a standard error column means "not available", not "zero" and not
+   "an error occurred".
+4. A custom measure gets a standard error for free if it is a mean or a sum of
+   its pointwise values, and needs an `se_diff_fun` otherwise.
+
+## References
+
+Newcombe, R. G. (2001). Simultaneous comparison of sensitivity and specificity
+of two tests in the paired design: a straightforward graphical approach.
+*Statistics in Medicine*, 20(6):907–915.
+
+Chen, F., Xue, Y., Tan, M. T., and Chen, P. (2015). Efficient statistical tests
+to compare Youden index: accounting for contingency correlation. *Statistics in
+Medicine*, 34(9):1560–1576.
+
+Vehtari, A., Gelman, A., and Gabry, J. (2017). Practical Bayesian model
+evaluation using leave-one-out cross-validation and WAIC. *Statistics and
+Computing*, 27(5):1413–1432.
+
+## See also
+
+- [Computing predictive performance measures](pred-measure-workflow.html) — how
+  to produce the objects compared here.
+- [Overview of scores and metrics](overview-measures.html) — definitions and
+  formulas for each measure.
+- `?"loo-glossary"` — the `measure_compare_meta` attribute in full.