Update loo_compare() to support new pred_measure API

.gitignore

@@ -30,4 +30,6 @@ release-prep.R
 
 # personal maintainer scratch (not shared)
 internal-notes/
+notes/loo_se.pdf
+notes/loo-compare-se-diff.md
 CRAN-SUBMISSION

DESCRIPTION

@@ -64,4 +64,3 @@ LazyData: TRUE
 Roxygen: list(markdown = TRUE)
 SystemRequirements: pandoc (>= 1.12.3), pandoc-citeproc
 Config/roxygen2/version: 8.0.0
-RoxygenNote: 7.3.3

NAMESPACE

@@ -48,6 +48,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)
@@ -158,6 +160,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,8 +7,111 @@
 * 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()].
-* Update user messages in `print()` by @ishaan-arora-1, @florence-bockting in 
-#328.
+* `loo_compare()` is now called `model_compare()`, and `loo_compare()` is
+  deprecated. It still works and is still a generic, so `loo_compare` methods
+  registered by other packages (e.g. `loo_compare.brmsfit` in **brms**,
+  `loo_compare.stanreg` in **rstanarm**) keep dispatching, but it now warns
+  (once per session) and is frozen at its previous behaviour: `"loo"`, `"waic"`, and `"kfold"` objects
+  compared on ELPD only, returning exactly what `model_compare()` returns for
+  those inputs. Comparing `pred_measure` results, or using `rank_by` or
+  `custom_se_fn`, requires `model_compare()`.
+* `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.
+* `print()` on a `pred_measure` comparison now always names the reference model
+  each set of differences was computed against, including the default case where
+  every measure keeps its own best model. The line is written to standard output
+  rather than emitted as a message, so it survives `suppressMessages()` and
+  knitr chunks with `message = FALSE`. For `loo_pred_measure()` comparisons the
+  PSIS-LOO Pareto k diagnostic moved out of the per-measure difference tables
+  into a block above them: it describes a model's LOO approximation, not any one
+  measure or pair of models, so it is now reported once per model. The
+  `diag_elpd` column on the returned object is unchanged, and `print()` for
+  plain `"loo"` comparisons is unchanged. Printed sentences are wrapped so that
+  no line of output exceeds 80 characters.
+* Bug fix: the `diag_diff` column tested `elpd_diff > -4` rather than
+  `|elpd_diff| < 4`, so every positive difference was flagged as a small one.
+  This was only reachable with `rank_by` set to a model that is not the best
+  model, where differences against the reference can be positive.
+* 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 `rank_by` 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_info`, 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_info`
+  attribute, and models that disagree on it cannot be compared.
+* Measure values are always stored on the measure's own scale; `loss` is the
+  single setting describing a measure's orientation. Measures no longer take a
+  `higher_is_better` argument, `control` no longer accepts one, and results no
+  longer carry a `measure_higher_is_better` attribute. Comparisons were already
+  invariant to the setting --- `model_compare()` reports every measure on the
+  utility scale regardless --- so no comparison output changes. Negate the
+  values yourself if you want a loss stored as a utility.
+* Bug fix: `rps` was registered as a utility and `srps` as a loss, but
+  `measure_rps()` returns the Gneiting & Raftery (2007) unscaled score, where
+  lower is better, and `measure_srps()` the Bolin & Wallin (2023) scaled score,
+  where higher is better. `model_compare()` therefore flipped `srps_diff` and
+  left `rps_diff` unflipped, and `rank_by = "rps"` ordered models backwards.
+  The registry now records `rps` as a loss and `srps` as a utility.
+* The per-measure metadata attribute `measure_compare_meta` is now called
+  `measure_info`.
+* `control` entries now work for custom measures instead of failing with
+  `object 'measure_<name>' of mode 'function' was not found`. 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.
+* `print()` on a `pred_measure` comparison now accepts `simplify`. With 
+  `simplify = FALSE` each measure table also shows that measure's per-model estimate and standard error, and an ELPD table also shows `p` and `se_p`.
+* `model_compare(A = m1, B = m2)` now works.
+* Update user messages in `print()` by @ishaan-arora-1, @florence-bockting in #328.
 
 # loo 2.10.1
 

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/helpers.R

@@ -182,18 +182,34 @@ loo_cores <- function(cores) {
   return(cores)
 }
 
-
-# nocov start
-# release reminders (for devtools)
-release_questions <- function() {
-  c(
-    "Have you updated references?",
-    "Have you updated inst/CITATION?",
-    "Have you updated the vignettes?"
-  )
-}
-# nocov end
-
 is_constant <- function(x, tol = .Machine$double.eps) {
   abs(max(x) - min(x)) < tol
 }
+
+#' Issue a deprecation warning the first time it is triggered in a session
+#'
+#' Repeated calls with the same `id` are silent, so a script calling a
+#' deprecated function in a loop is not flooded with warnings. `old` is passed
+#' on explicitly so the message does not depend on which method called this.
+#' Which `id`s have already warned is kept in `state`, an environment created
+#' once when the package is built and private to this function.
+#'
+#' @noRd
+#' @param id Identifier for the deprecation; one warning per `id` per session.
+#' @param new,old Name of the replacement and deprecated function.
+#' @return `TRUE` if a warning was issued, `FALSE` otherwise, invisibly.
+#'
+.deprecate_once <- local({
+  state <- new.env(parent = emptyenv())
+  function(id, new, old = id) {
+    if (isTRUE(state[[id]])) {
+      return(invisible(FALSE))
+    }
+    state[[id]] <- TRUE
+    warning(
+      "\n'", old, "' is deprecated. Use '", new, "' instead.\n",
+      call. = TRUE, immediate. = TRUE
+    )
+    invisible(TRUE)
+  }
+})

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