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..1913b517 100644 --- a/R/loo-glossary.R +++ b/R/loo-glossary.R @@ -254,4 +254,80 @@ #' #' 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`. When the overall estimate is a sum or mean of pointwise +#' contributions, these are computed from paired pointwise differences 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) using the same +#' approach as `elpd_diff` and `se_diff` (Eq 24 in VGG2017 for sums; the mean +#' analogue for means). 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`. For sum- and +#' mean-based measures, the reference model also has `m_se_diff = 0`; for +#' `estimates_only` measures (e.g. `r2`, `mse`, `rmse`), `m_se_diff` is `NA`. +#' 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()]). +#' +#' For measures where pointwise values do not define the overall estimate (e.g. +#' `r2`, `mse`, `rmse`), `m_diff` is the difference between overall estimates +#' (on a utility scale) and `m_se_diff` is `NA`. +#' +#' 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 paired differences are aggregated: `"sum"`, +#' `"mean"`, `"estimates_only"`, or `"auto"` (inferred at compare time for +#' custom measures) +#' +#' Built-in measures take `loss` and `diff_method` from the package measure +#' registry; custom measures default to `loss = FALSE` and `diff_method = "auto"`. +#' [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, `sign_converted_measures` lists loss measures whose sign was +#' flipped onto the utility scale, and `measures_no_pointwise_se` lists measures +#' for which `{measure}_se_diff` is unavailable (overall estimate not defined +#' from pointwise values). 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..9f2ee882 100644 --- a/R/loo_compare.R +++ b/R/loo_compare.R @@ -1,17 +1,40 @@ #' 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"`). +#' Attribute `measures_no_pointwise_se` lists measures without pointwise-based +#' `{measure}_se_diff` values. +#' #' @details #' When comparing two fitted models, we can estimate the difference in their #' expected predictive accuracy by the difference in @@ -67,9 +90,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 +154,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 +270,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 +313,157 @@ 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\")." + ) + } + } + + .warn_measures_no_pointwise_se(attr(x, "measures_no_pointwise_se")) + + 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 +471,445 @@ 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") +} -# internal ---------------------------------------------------------------- +#' 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 +} + +#' 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() + measures_no_pointwise_se <- list() + for (col in compare_cols) { + bare <- .display_name(col) + method <- .measure_pointwise_diff_method(loos_ord, col) + if (method == "estimates_only") { + measures_no_pointwise_se[[bare]] <- bare + } + 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 + ) + + attr(comp, "measures_no_pointwise_se") <- unique(unlist(measures_no_pointwise_se)) + 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 + ) + } + non_null <- metas[has_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 + ) +} + +#' Warn when `se_diff` is unavailable for compared measures +#' @noRd +.warn_measures_no_pointwise_se <- function(measures) { + if (!length(measures)) { + return(invisible(NULL)) + } + warning( + paste0( + "se_diff unavailable for: ", + paste(measures, collapse = ", "), + "." + ), + call. = FALSE + ) + invisible(NULL) +} + +#' 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 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 + loss <- NULL + + if (!is.null(loos)) { + meta <- .get_measure_compare_meta(loos, bare) + if (!is.null(meta)) { + higher_is_better <- meta$higher_is_better + loss <- meta$loss + } 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)) + } + + if (is.null(loss)) { + spec <- .measure_spec[[bare]] + loss <- if (!is.null(spec)) { + isTRUE(spec$loss) + } else { + bare %in% c("ic", "mae", "mse", "rmse", "brier", "srps") + } + } + + isTRUE(loss) +} + +#' 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", + " reported 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, and `"estimates_only"` +#' when pointwise values do not define the overall estimate. +#' @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) + } + + spec <- .measure_spec[[bare]] + if (!is.null(spec) && identical(spec$diff_method, "estimates_only")) { + return("estimates_only") + } + 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" +} + +#' 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) + + if (method == "estimates_only") { + est_utility <- function(estimates) { + val <- estimates[col, "Estimate"] + if (flip) -val else val + } + return(c( + diff = est_utility(cmp$estimates) - est_utility(ref$estimates), + se = NA_real_ + )) + } + + 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 +935,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 +967,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 (all(sapply(loos, is.kfold))) { + if (!kfold_checks) { + return(invisible(NULL)) + } + + 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 +1032,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 +1112,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..e7537e52 100644 --- a/R/pred_measure-builtin.R +++ b/R/pred_measure-builtin.R @@ -940,20 +940,52 @@ 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 = "estimates_only"), + rmse = list(fun = measure_rmse, loss = TRUE, diff_method = "estimates_only"), + mse = list(fun = measure_mse, loss = TRUE, diff_method = "estimates_only"), + acc = list(fun = measure_acc, loss = FALSE, diff_method = "mean"), + bacc = list(fun = measure_bacc, loss = FALSE, diff_method = "estimates_only"), + 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") ) +#' 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 + )) + } + + list( + higher_is_better = higher_is_better, + loss = FALSE, + diff_method = "auto" + ) +} + #' Supported predictive measure names #' #' A character vector of measure names that can be passed to the `measure` @@ -991,6 +1023,7 @@ 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 ) } diff --git a/R/pred_measure-compute.R b/R/pred_measure-compute.R index 9fbbcd54..14585023 100644 --- a/R/pred_measure-compute.R +++ b/R/pred_measure-compute.R @@ -165,7 +165,9 @@ 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") ) pointwise <- .merge_matrix( source = source, @@ -184,7 +186,7 @@ do_pred_measure <- function( save_psis = save_psis ) - .add_attributes( + predperf_res <- .add_attributes( save_psis, predperf_res, y, @@ -197,6 +199,7 @@ do_pred_measure <- function( predperf, source ) + predperf_res } # internal helper functions --------------------------------------------------- @@ -323,10 +326,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 +540,11 @@ 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`). #' #' @return Updated matrix with `name` as a row or column name. #' @@ -551,7 +560,15 @@ 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 +) { is_row <- margin == 1 bind_fn <- if (is_row) rbind else cbind name_updated <- .measure_result_name(source, name) @@ -562,8 +579,43 @@ 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) + } + + 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 +638,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 +651,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 +676,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 +721,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 +757,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 +823,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.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/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..7b867ba1 100644 --- a/man/insample_pred_measure.Rd +++ b/man/insample_pred_measure.Rd @@ -67,7 +67,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 +106,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..f31b3936 100644 --- a/man/kfold_pred_measure.Rd +++ b/man/kfold_pred_measure.Rd @@ -78,6 +78,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..e7b0b0e7 100644 --- a/man/loo-glossary.Rd +++ b/man/loo-glossary.Rd @@ -270,6 +270,88 @@ 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}. When the overall estimate is a sum or mean of pointwise +contributions, these are computed from paired pointwise differences 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) using the same +approach as \code{elpd_diff} and \code{se_diff} (Eq 24 in VGG2017 for sums; the mean +analogue for means). 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}. For sum- and +mean-based measures, the reference model also has \code{m_se_diff = 0}; for +\code{estimates_only} measures (e.g. \code{r2}, \code{mse}, \code{rmse}), \code{m_se_diff} is \code{NA}. +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()}}). + +For measures where pointwise values do not define the overall estimate (e.g. +\code{r2}, \code{mse}, \code{rmse}), \code{m_diff} is the difference between overall estimates +(on a utility scale) and \code{m_se_diff} is \code{NA}. + +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 paired differences are aggregated: \code{"sum"}, +\code{"mean"}, \code{"estimates_only"}, or \code{"auto"} (inferred at compare time for +custom measures) +} + +Built-in measures take \code{loss} and \code{diff_method} from the package measure +registry; custom measures default to \code{loss = FALSE} and \code{diff_method = "auto"}. +\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, \code{sign_converted_measures} lists loss measures whose sign was +flipped onto the utility scale, and \code{measures_no_pointwise_se} lists measures +for which \verb{\{measure\}_se_diff} is unavailable (overall estimate not defined +from pointwise values). 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..b7d9c8f0 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,38 @@ 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"}). +Attribute \code{measures_no_pointwise_se} lists measures without pointwise-based +\verb{\{measure\}_se_diff} values. } \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 +129,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 +192,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..389a48ae 100644 --- a/man/loo_pred_measure.Rd +++ b/man/loo_pred_measure.Rd @@ -113,6 +113,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..fe8880da 100644 --- a/man/pred_measure.Rd +++ b/man/pred_measure.Rd @@ -69,7 +69,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 +101,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/test_pred_measure.Rd b/man/test_pred_measure.Rd index 0b3710e9..6b93f97b 100644 --- a/man/test_pred_measure.Rd +++ b/man/test_pred_measure.Rd @@ -79,6 +79,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 *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..ec2506f5 100644 --- a/tests/testthat/_snaps/compare.md +++ b/tests/testthat/_snaps/compare.md @@ -1,3 +1,131 @@ +# 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"). + Condition + Warning: + se_diff unavailable for: r2. + +# 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"). + Condition + Warning: + se_diff unavailable for: r2. + +--- + + 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 NA + m2 0.01 NA + m1 -0.20 NA + + -- 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. + Condition + Warning: + se_diff unavailable for: r2. + +--- + + Code + print(comp, measures = c("r2", "mae")) + Output + + -- r2 (vs m3) -- + model r2_diff se_diff + m3 0.0 NA + m2 0.0 NA + m1 -0.2 NA + + -- 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. + Condition + Warning: + se_diff unavailable for: r2. + +--- + + 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"). + Condition + Warning: + se_diff unavailable for: r2. + # 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..411e7999 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,458 @@ 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_true(all(is.na(comp$r2_se_diff))) + expect_true(all(is.na(comp$mse_se_diff))) + expect_false("r2_loo_diff" %in% colnames(comp)) + expect_false("mse_p_worse" %in% colnames(comp)) + expect_warning( + print(comp), + "se_diff unavailable for: r2, mse" + ) + + 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, "estimates_only") + expect_equal(attr(pm1, "measure_compare_meta")$r2$diff_method, "estimates_only") + 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", "estimates_only", 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", "estimates_only", 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"), "estimates_only") + expect_equal(loo:::.measure_pointwise_diff_method(loos, "r2_loo"), "estimates_only") + expect_equal(loo:::.measure_pointwise_diff_method(loos, "elpd_loo"), "sum") + expect_true(is.na(pair_mse["se"])) +}) + +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