add diff diagnostics

R/loo_compare.R

@@ -2,33 +2,30 @@
 #'
 #' @description Compare fitted models based on [ELPD][loo-glossary].
 #'
-#'   By default the print method shows only the most important information. Use
-#'   `print(..., simplify=FALSE)` to print a more detailed summary.
-#'
 #' @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.
 #'
-#' @return A matrix with class `"compare.loo"` that has its own
-#'   print method. See the **Details** section.
+#' @return A data frame with class `"compare.loo"` that has its own
+#'   print method. See the **Details** and **Examples** sections.
 #'
 #' @details
 #'   When comparing two fitted models, we can estimate the difference in their
 #'   expected predictive accuracy by the difference in
 #'   [`elpd_loo`][loo-glossary] or `elpd_waic` (or multiplied by \eqn{-2}, if
 #'   desired, to be on the deviance scale).
 #'
-#'   When using `loo_compare()`, the returned matrix will have one row per model
-#'   and several columns of estimates. The values in the
-#'   [`elpd_diff`][loo-glossary] and [`se_diff`][loo-glossary] columns of the
-#'   returned matrix are computed by making pairwise comparisons between each
-#'   model and the model with the largest ELPD (the model in the first row). For
-#'   this reason the `elpd_diff` column will always have the value `0` in the
-#'   first row (i.e., the difference between the preferred model and itself) and
-#'   negative values in subsequent rows for the remaining models.
+#' ## `elpd_diff` and `se_diff`
+#'   When using `loo_compare()`, the returned data frame will have one row per
+#'   model and several columns of estimates. The values of
+#'   [`elpd_diff`][loo-glossary] and [`se_diff`][loo-glossary] are computed by
+#'   making pairwise comparisons between each model and the model with the
+#'   largest ELPD (the model listed first). Therefore, the first `elpd_diff`
+#'   value will always be `0` (i.e., the difference between the preferred model
+#'   and itself) and the rest of the values will be negative.
 #'
 #'   To compute the standard error of the difference in [ELPD][loo-glossary] ---
 #'   which should not be expected to equal the difference of the standard errors
@@ -41,9 +38,41 @@
 #'   standard approach of comparing differences of deviances to a Chi-squared
 #'   distribution, a practice derived for Gaussian linear models or
 #'   asymptotically, and which only applies to nested models in any case.
-#'   Sivula et al. (2022) discuss the conditions when the normal
-#'   approximation used for SE and `se_diff` is good.
 #'
+#' ## `p_worse`, `diag_diff`, and `diag_elpd`
+#'   The values in the `p_worse` column show the probability of each model
+#'   having worse ELPD than the best model. These probabilities are computed
+#'   with a normal approximation using the values from `elpd_diff` and
+#'   `se_diff`. Sivula et al. (2025) present the conditions when the normal
+#'   approximation used for SE and `se_diff` is good, and the column
+#'   `diag_diff` contains possible diagnostic messages:
+#'
+#'   * `N < 100` (small data)
+#'   * `|elpd_diff| < 4` (models make similar predictions)
+#'   * `khat > 0.5` (possible outliers)
+#'
+#'   If any of these diagnostic messages is shown, the error distribution is
+#'   skewed or thick tailed and the normal approximation based on `elpd_diff`
+#'   and `se_diff` is not well calibrated. The probabilities `p_worse` are
+#'   likely to be too large (small data or similar predictions) or too small
+#'   (outliers). `elpd_diff` and `se_diff` are still indicative of the
+#'   differences and uncertainties, and for example, if `|elpd_diff|` is
+#'   many times larger than `se_diff` the difference is quite certain.
+#'   While `khat > 0.5` indicates possibility of outliers, it is also
+#'   possible that both models compared seem to be well specified based
+#'   on model checking, but the pointwise ELPD differences have such thick
+#'   tails that the normal approximation for the sum is not good (Vehtari et
+#'   al., 2024). Note that `khat`s in `diag_diff` column are different from
+#'   the `khat`s in PSIS-LOO  diagnostic. 
+#'
+#'   The column `diag_elpd` shows diagnostic for the pointwise ELPD
+#'   computations for each model. If `k khat_psis > 0.7` is shown,
+#'   where `k` is the number of high high Pareto k values in Pareto
+#'   smoothed importance sampling computation, then there may be
+#'   significant bias in `elpd_diff` favoring models with a large
+#'   number of high Pareto k values.
+#'
+#' ## 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
 #'   estimate whether the differences in predictive performance are potentially
@@ -52,7 +81,7 @@
 #'   selection process. In that case users are recommended to avoid model
 #'   selection based on LOO-CV, and instead to favor model averaging/stacking or
 #'   projection predictive inference.
-#' 
+#'
 #' @seealso
 #' * The [FAQ page](https://mc-stan.org/loo/articles/online-only/faq.html) on
 #'   the __loo__ website for answers to frequently asked questions.
@@ -68,12 +97,8 @@
 #' comp <- loo_compare(loo1, loo2, loo3)
 #' print(comp, digits = 2)
 #'
-#' # show more details with simplify=FALSE
-#' # (will be the same for all models in this artificial example)
-#' print(comp, simplify = FALSE, digits = 3)
-#'
 #' # can use a list of objects with custom names
-#' # will use apple, banana, and cherry, as the names in the output
+#' # the names will be used in the output
 #' loo_compare(list("apple" = loo1, "banana" = loo2, "cherry" = loo3))
 #'
 #' \dontrun{
@@ -101,54 +126,122 @@ loo_compare.default <- function(x, ...) {
     loos <- x
   }
 
-  # If subsampling is used
+  # if subsampling is used
   if (any(sapply(loos, inherits, "psis_loo_ss"))) {
     return(loo_compare.psis_loo_ss_list(loos))
   }
 
+  # 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)
-
-  # compute elpd_diff and se_elpd_diff relative to best model
   rnms <- rownames(comp)
   diffs <- mapply(FUN = elpd_diffs, loos[ord[1]], loos[ord])
+  colnames(diffs) <- rnms
   elpd_diff <- apply(diffs, 2, sum)
   se_diff <- apply(diffs, 2, se_elpd_diff)
-  comp <- cbind(elpd_diff = elpd_diff, se_diff = se_diff, comp)
-  rownames(comp) <- rnms
 
-  # run order statistics-based checks on models
+  # compute probabilities that a model has worse elpd than the best model
+  # using a normal approximation (Sivula et al., 2025)
+  p_worse <- stats::pnorm(0, elpd_diff, se_diff)
+  p_worse[elpd_diff == 0] <- NA
+
+  # diagnostics to assess whether the normal approximation can be trusted
+  # * N < 100: small data (Sivula et al., 2025)
+  # * |elpd_diff| < 4: similar predictions (Sivula et al., 2025)
+  # * khat_diff > 0.5: possible outliers in differences (Sivula et al., 2025; Vehtari et al., 2024)
+  N <- nrow(diffs)
+  if (N < 100) {
+    diag_diff <- rep("N < 100", length(elpd_diff))
+    diag_diff[elpd_diff == 0] <- ""
+  } else {
+    diag_diff <- rep("", length(elpd_diff))
+    diag_diff[elpd_diff > -4 & elpd_diff != 0] <- "|elpd_diff| < 4"
+    khat_diff <- rep(NA, length(elpd_diff))
+    khat_diff[elpd_diff != 0] <- apply(
+      diffs[, elpd_diff != 0, drop = FALSE], 2,
+      function(x) ifelse(length(unique(x)) <= 20, NA, posterior::pareto_khat(x, tail = "both")
+    ))
+    diag_diff[khat_diff > 0.5] <- "khat_diff > 0.5"
+  }
+
+  # get khats for PSIS
+  khat_psis <- sapply(loos[ord],
+                      \(loo) {
+                        k <- loo$diagnostics[["pareto_k"]]
+                        if (is.null(k)) {
+                          out = ""
+                        } else {
+                          S <- dim(loo)[1]
+                          khat_threshold <- ps_khat_threshold(S)
+                          K <- sum(k > khat_threshold)
+                          if (K==0) {
+                            out <- ""
+                          } else {
+                            out <- paste0(K, " khat_psis > ", round(khat_threshold, 2))
+                          }
+                        }
+                        out
+                      }
+                      )
+
+  comp <- cbind(
+    data.frame(
+      model = rnms,
+      elpd_diff = elpd_diff,
+      se_diff = se_diff,
+      p_worse = p_worse,
+      diag_diff = diag_diff,
+      diag_elpd = khat_psis
+    ),
+    as.data.frame(comp)
+  )
+  rownames(comp) <- NULL
+
+  # run order statistics-based checks for many model comparisons
   loo_order_stat_check(loos, ord)
 
   class(comp) <- c("compare.loo", class(comp))
-  return(comp)
+  comp
 }
 
 #' @rdname loo_compare
 #' @export
 #' @param digits For the print method only, the number of digits to use when
 #'   printing.
-#' @param simplify For the print method only, should only the essential columns
-#'   of the summary matrix be printed? The entire matrix is always returned, but
-#'   by default only the most important columns are printed.
-print.compare.loo <- function(x, ..., digits = 1, simplify = TRUE) {
-  xcopy <- x
-  if (inherits(xcopy, "old_compare.loo")) {
-    if (NCOL(xcopy) >= 2 && simplify) {
-      patts <- "^elpd_|^se_diff|^p_|^waic$|^looic$"
-      xcopy <- xcopy[, grepl(patts, colnames(xcopy))]
-    }
-  } else if (NCOL(xcopy) >= 2 && simplify) {
-     xcopy <- xcopy[, c("elpd_diff", "se_diff")]
+#' @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) {
+  if (inherits(x, "old_compare.loo")) {
+    return(unclass(x))
   }
-  print(.fr(xcopy, digits), quote = FALSE)
+  if (!inherits(x, "data.frame")) {
+    class(x) <- c(class(x), "data.frame")
+  }
+  if (!all(c("model", "elpd_diff", "se_diff") %in% colnames(x))) {
+    print(as.data.frame(x))
+    return(x)
+  }
+  x2 <- cbind(
+    model = x$model,
+    .fr(x[, c("elpd_diff", "se_diff")], digits)
+  )
+  if (p_worse && "p_worse" %in% colnames(x)) {
+    x2 <- cbind(
+      x2,
+      p_worse = .fr(x[, "p_worse"], digits = 2),
+      diag_diff = x[, "diag_diff"],
+      diag_elpd = x[, "diag_elpd"]
+    )
+  }
+  print(x2, quote = FALSE, row.names = FALSE)
   invisible(x)
 }
 
 
-
 # internal ----------------------------------------------------------------
 
 #' Compute pointwise elpd differences
@@ -172,7 +265,6 @@ se_elpd_diff <- function(diffs) {
   sqrt(N) * sd(diffs)
 }
 
-
 #' Perform checks on `"loo"` objects before comparison
 #' @noRd
 #' @param loos List of `"loo"` objects.

R/loo_compare.psis_loo_ss_list.R

@@ -173,14 +173,9 @@ loo_compare_checks.psis_loo_ss_list <- function(loos) {
 
 #' @rdname loo_compare
 #' @export
-print.compare.loo_ss <- function(x, ..., digits = 1, simplify = TRUE) {
+print.compare.loo_ss <- function(x, ..., digits = 1) {
   xcopy <- x
-  if (inherits(xcopy, "old_compare.loo")) {
-    if (NCOL(xcopy) >= 2 && simplify) {
-      patts <- "^elpd_|^se_diff|^p_|^waic$|^looic$"
-      xcopy <- xcopy[, grepl(patts, colnames(xcopy))]
-    }
-  } else if (NCOL(xcopy) >= 2 && simplify) {
+  if (NCOL(xcopy) >= 2) {
     xcopy <- xcopy[, c("elpd_diff", "se_diff", "subsampling_se_diff")]
   }
   print(.fr(xcopy, digits), quote = FALSE)