Update loo_compare() to support new pred_measure API

.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

NAMESPACE

@@ -47,6 +47,8 @@ S3method(loo_moment_match,default)
 S3method(loo_predictive_metric,matrix)
 S3method(loo_scrps,matrix)
 S3method(loo_subsample,"function")
+S3method(model_compare,default)
+S3method(model_compare,psis_loo_ss_list)
 S3method(nobs,psis_loo_ss)
 S3method(plot,loo)
 S3method(plot,psis)
@@ -159,6 +161,7 @@ export(measure_r2)
 export(measure_rmse)
 export(measure_rps)
 export(measure_srps)
+export(model_compare)
 export(nlist)
 export(obs_idx)
 export(pareto_k_ids)

NEWS.md

@@ -7,6 +7,74 @@
 * 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()` is now called `model_compare()`. `loo_compare()` remains a
+  working alias and is still a generic, so existing code and methods registered
+  by other packages keep working.
+* `model_compare()` supports every `pred_measure` result --- from
+  `loo_pred_measure()`, `kfold_pred_measure()`, `test_pred_measure()`, and
+  `insample_pred_measure()` --- with 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.
+  Measures are matched on bare names, so the source suffix (`_loo`, `_kfold`,
+  `_test`, or none for in-sample) is handled transparently. All models in one
+  call must share an evaluation source; mixing them is an error. `diag_elpd` is
+  reported only for `loo_pred_measure()` comparisons, k-fold comparisons warn on
+  differing `K`, and in-sample comparisons warn that they are optimistically
+  biased.
+* In `model_compare()`, when `rank_by` is not supplied each measure is now
+  compared against *its own* best model, so e.g. `mse_diff` can be relative to a
+  different model than `elpd_diff`. Rows are still ordered by `"elpd"`, and each
+  `{measure}_diff` column has exactly one `0` entry, at that measure's best
+  model. Supplying `rank_by` keeps the previous behaviour of pinning the
+  top-ranked model as the single reference for every measure. The reference used
+  per measure is recorded in the `compare_reference` attribute and shown by
+  `print(x, measures = "all")`. Each printed measure table is now sorted by its
+  own difference, so the best model on that measure is always the first row and
+  the differences run in decreasing order.
+* In `model_compare()`, `rank_by` now also accepts a **model name** (as shown in
+  the `model` column). The named model then becomes the single reference model
+  for every measure, whether or not it is the best one, while rows stay ordered
+  by `"elpd"`. This form also works for plain `"loo"` comparisons, where
+  `elpd_diff` is then relative to the named model. The pinned model is recorded
+  in the new `compare_ref_model` attribute; a name matching both a measure and a
+  model is treated as the measure, with a warning.
+* The `diff_method` value `"pairwise"` is now called `"measure_specific"`. Every
+  model comparison is pairwise, and the `"sum"`/`"mean"` methods are the ones
+  built from paired pointwise differences, so the old name described the one
+  category that is *not* a paired pointwise standard deviation. `diff_method` is
+  read-only metadata in `measure_compare_meta`, so no user code that passes
+  arguments is affected; `custom_se_fn` still accepts only `"sum"` and `"mean"`.
+* Bug fix: in `model_compare()`, `mlpd` was registered with
+  `diff_method = "sum"` although its estimate is the *mean* of the pointwise log
+  predictive densities. `mlpd_diff` was therefore reported as the sum of the
+  pointwise differences (a factor of `N` too large, and inconsistent with the
+  difference of the reported `mlpd` estimates), with a matching `mlpd_se_diff`.
+  `mlpd` now uses `diff_method = "mean"`.
+* `model_compare()` gains a `custom_se_fn` argument controlling how the standard
+  error of a difference is computed for a **custom** measure. It accepts a
+  function called as `custom_se_fn(ref, cmp)`, the shorthands `"sum"` and
+  `"mean"` for the paired pointwise formulas, or `NULL` to report the difference
+  with an `NA` standard error; with two or more custom measures, pass a list
+  named by measure. It is required whenever a custom measure is compared ---
+  nothing is inferred from a measure's values any more, and custom measures now
+  carry `diff_method = "custom"` instead of `"auto"`. The previous
+  `attr(my_fun, "se_diff_fun")` route has been removed.
+* A custom measure can declare that it is a loss with
+  `attr(my_fun, "measure_loss") <- TRUE`, alongside `attr(my_fun,
+  "measure_name")`. `model_compare()` then flips its differences onto the
+  utility scale and orders `rank_by` by the lowest loss, as it does for built-in
+  loss measures; without the declaration a custom measure is still treated as a
+  utility. The declaration is recorded as `loss` in the `measure_compare_meta`
+  attribute, and models that disagree on it cannot be compared.
+* `control` entries now work for custom measures instead of failing with
+  `object 'measure_<name>' of mode 'function' was not found`. `higher_is_better`
+  is accepted for any measure, built-in or custom, and selects the scale its
+  values are stored on. A `control` entry naming no measure being computed now
+  warns instead of aborting.
+* `print()` on a multi-measure comparison now labels each measure's standard
+  error column `{measure}_se_diff`, matching the column name on the returned
+  data frame, instead of relabelling it `se_diff`. ELPD tables are unchanged.
 
 # loo 2.10.0
 

R/compare.R

@@ -1,6 +1,6 @@
 #' Model comparison (deprecated, old version)
 #'
-#' **This function is deprecated**. Please use the new [loo_compare()] function
+#' **This function is deprecated**. Please use the new [model_compare()] function
 #' instead. See `vignette("migration-guide", package = "loo")` for details.
 #' `compare()` and `R/compare.R` are scheduled for removal in a future release.
 #'
@@ -60,7 +60,7 @@
 #' }
 #'
 compare <- function(..., x = list()) {
-  .Deprecated("loo_compare")
+  .Deprecated("model_compare")
   dots <- list(...)
   if (length(dots)) {
     if (length(x)) {

R/kfold-generic.R

@@ -7,7 +7,7 @@
 #'
 #'   The **Value** section below describes the objects that `kfold()`
 #'   methods should return in order to be compatible with
-#'   [loo_compare()] and the **loo** package print methods.
+#'   [model_compare()] and the **loo** package print methods.
 #'
 #'
 #' @name kfold-generic
@@ -25,7 +25,7 @@
 #'
 #'   It is important for the object to have at least these classes and
 #'   components so that it is compatible with other functions like
-#'   [loo_compare()] and `print()` methods.
+#'   [model_compare()] and `print()` methods.
 #'
 NULL
 

R/loo-glossary.R

@@ -161,7 +161,7 @@
 #'
 #'     p_worse = pnorm(0, elpd_diff, se_diff).
 #'
-#' The best-ranked model (the first row in the `loo_compare()` output, where
+#' The best-ranked model (the first row in the `model_compare()` output, where
 #' `elpd_diff = 0`) always receives `NA`, since the comparison is defined
 #' relative to that model.
 #'
@@ -176,15 +176,15 @@
 #' appear more clearly worse than the data actually support. Conversely, when
 #' `elpd_diff` is biased due to an unreliable LOO approximation, `p_worse` can
 #' point in the wrong direction entirely. When any of these conditions are
-#' present, `diag_diff` or `diag_elpd` will be flagged in the `loo_compare()`
+#' present, `diag_diff` or `diag_elpd` will be flagged in the `model_compare()`
 #' output. 
 #' For further guidance, see the sections below and the case study on
 #' [Uncertainty in Bayesian LOO-CV Model Comparison](
 #' https://users.aalto.fi/~ave/casestudies/LOO_uncertainty/loo_uncertainty.html).
 #'
 #' @section `diag_diff` (pairwise comparison diagnostics):
 #'
-#' `diag_diff` is a diagnostic column in the `loo_compare()` output for each
+#' `diag_diff` is a diagnostic column in the `model_compare()` output for each
 #' model comparison against the current reference model. It flags conditions
 #' under which the normal approximation behind `se_diff` and `p_worse` is likely
 #' to be poorly calibrated. The column contains a short label when a condition
@@ -210,7 +210,7 @@
 #'
 #' The conditions flagged by `diag_diff` are not independent: they tend to
 #' co-occur, and when they do, some flags carry more information than others.
-#' `loo_compare()` therefore follows a priority hierarchy and shows only the
+#' `model_compare()` therefore follows a priority hierarchy and shows only the
 #' most critical flag in the table output.
 #'
 #' The hierarchy is as follows:
@@ -233,7 +233,7 @@
 #' 
 #' @section `diag_elpd`:
 #'
-#' `diag_elpd` is a diagnostic column in the `loo_compare()` output that flags
+#' `diag_elpd` is a diagnostic column in the `model_compare()` output that flags
 #' when the PSIS-LOO approximation for an individual model is unreliable. Unlike
 #' `diag_diff`, which concerns the *comparison* between models, `diag_elpd`
 #' concerns the quality of the `elpd_loo` estimate for each model individually.
@@ -254,4 +254,125 @@
 #'
 #' 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
+#' `model_compare()`, paired differences are computed for every predictive
+#' measure common to all models. Rows are ordered by the `rank_by` argument
+#' (default `"elpd"`). By default each measure is compared against the model
+#' that is best on that measure, so different difference columns may use
+#' different reference models. Supplying `rank_by` pins the top-ranked model as
+#' the single reference for all difference columns.
+#'
+#' ### `{measure}_diff` and `{measure}_se_diff`
+#'
+#' For each non-ELPD measure `m`, `model_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 a built-in measure is a transformation of such quantities, it 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.
+#' * For custom measures it comes from the `custom_se_fn` argument of
+#'   [model_compare()], which must be supplied whenever a custom measure is
+#'   compared. It is `NA` when `custom_se_fn` is `NULL` for that measure.
+#'
+#' 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, `model_compare()` emits a
+#' short message naming those measures (see [model_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 (`TRUE`
+#' or `FALSE`). Measures left at their natural orientation have no entry, which
+#' reads as `NULL`; `elpd` is always `NULL`. Used by [model_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 [model_compare()]. Each entry
+#' is a list with:
+#'
+#' * `higher_is_better` — the orientation setting used when the measure was
+#'   computed (`NULL`, `TRUE`, or `FALSE`), that is, which scale the values were
+#'   stored on
+#' * `loss` — whether the measure itself is a loss (lower is better), regardless
+#'   of the scale its values are stored on
+#' * `diff_method` — how the standard error of the difference is obtained:
+#'   `"sum"` or `"mean"` (paired pointwise differences),
+#'   `"measure_specific"` (the built-in measure's own `se_diff_fun`), or
+#'   `"custom"`. Nothing is inferred from a measure's values. Under `"custom"` the standard error is resolved at
+#'   comparison time from the `custom_se_fn` argument of [model_compare()],
+#'   which supplies either a function, the `"sum"`/`"mean"` pointwise formulas,
+#'   or `NULL` for an `NA` standard error. A missing standard error is not an
+#'   error state — the difference itself is still reported.
+#'   `"estimates_only"` appears only for legacy objects that carry no
+#'   comparison metadata at all.
+#' * `se_diff_fun` — for built-in measures with
+#'   `diff_method = "measure_specific"`, the name of the built-in implementation
+#'   used. Custom measures never store a function here.
+#' * `extra` — optional list of auxiliary data the measure stored for the
+#'   standard error of its difference, 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 [model_compare()] runs; `bacc` stores the class
+#'   index of each observation, which its pointwise values do not determine).
+#'   Custom measures return it as an `extra` element, and it is passed on to
+#'   `custom_se_fn`. 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 always get `diff_method = "custom"`
+#' and take `loss` from `attr(my_fun, "measure_loss") <- TRUE`, which declares
+#' that lower values are better; without it they are treated as utilities (see
+#' [insample_pred_measure()]).
+#' [model_compare()] requires all models to provide matching metadata for each
+#' shared measure; mismatched `higher_is_better` or `measure_loss` settings, or
+#' missing metadata on some models, produce an error.
+#'
+#' ### `rank_by`, `compare_measures`, and related attributes`
+#'
+#' The `rank_by` argument takes either a measure name or a model name. A
+#' measure name selects which measure determines model ordering, and pins the
+#' top-ranked model as the single reference model for all pairwise differences.
+#' A model name keeps the `"elpd"` ordering but pins that model as the single
+#' reference model, whether or not it is the best one. When `rank_by` is
+#' omitted, models are ordered by `"elpd"` and
+#' each measure is compared against its own best model; attribute `rank_by` is
+#' set only when `rank_by` names a measure, and attribute `compare_ref_model`
+#' only when it names a model. Attribute `compare_reference`
+#' is a named character vector recording the reference model used for each
+#' measure. 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.
+#' Each printed table is sorted by its own measure, best model first, so the
+#' same model need not lead every table.
 NULL

R/loo.R

@@ -95,7 +95,7 @@
 #'  * [psis()] for the underlying Pareto Smoothed Importance Sampling (PSIS)
 #'    procedure used in the LOO-CV approximation.
 #'  * [pareto-k-diagnostic] for convenience functions for looking at diagnostics.
-#'  * [loo_compare()] for model comparison.
+#'  * [model_compare()] for model comparison.
 #'
 #' @template loo-and-psis-references
 #'

R/loo_approximate_posterior.R

@@ -27,7 +27,7 @@
 #'  }
 #' }
 #'
-#' @seealso [loo()], [psis()], [loo_compare()]
+#' @seealso [loo()], [psis()], [model_compare()]
 #' @template loo-large-data-references
 #'
 #' @export loo_approximate_posterior