From 040775e7e545df12c16c234149709498f084da4c Mon Sep 17 00:00:00 2001
From: Dominique Makowski <dom.mak19@gmail.com>
Date: Mon, 25 Mar 2024 17:47:20 +0000
Subject: [PATCH] Add report.compare.loo

---
 DESCRIPTION               |  1 +
 NAMESPACE                 |  1 +
 R/report.compare.loo.R    | 59 +++++++++++++++++++++++++++++++++++++++
 man/report.compare.loo.Rd | 39 ++++++++++++++++++++++++++
 4 files changed, 100 insertions(+)
 create mode 100644 R/report.compare.loo.R
 create mode 100644 man/report.compare.loo.Rd

diff --git a/DESCRIPTION b/DESCRIPTION
index 35d6c831..a38edffe 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -107,6 +107,7 @@ Collate:
     'report.stanreg.R'
     'report.brmsfit.R'
     'report.character.R'
+    'report.compare.loo.R'
     'report.compare_performance.R'
     'report.data.frame.R'
     'report.default.R'
diff --git a/NAMESPACE b/NAMESPACE
index a3cb2ff5..16c56e96 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -36,6 +36,7 @@ S3method(report,bayesfactor_inclusion)
 S3method(report,bayesfactor_models)
 S3method(report,brmsfit)
 S3method(report,character)
+S3method(report,compare.loo)
 S3method(report,compare_performance)
 S3method(report,data.frame)
 S3method(report,default)
diff --git a/R/report.compare.loo.R b/R/report.compare.loo.R
new file mode 100644
index 00000000..fc9d8f70
--- /dev/null
+++ b/R/report.compare.loo.R
@@ -0,0 +1,59 @@
+#' Reporting Bayesian Model Comparison
+#'
+#' Automatically report the results of Bayesian model comparison using the `loo` package.
+#'
+#' @param x An object of class [loo_compare()].
+#'
+#' @examplesIf require("brms", quietly = TRUE)
+#' \donttest{
+#' library(brms)
+#'
+#' m1 <- brms::brm(mpg ~ qsec, data=mtcars)
+#' m2 <- brms::brm(mpg ~ qsec + drat, data=mtcars)
+#'
+#' x <- brms::loo_compare(brms::add_criterion(m1, "loo"),
+#'                        brms::add_criterion(m2, "loo"),
+#'                        model_names=c("m1", "m2"))
+#' report(x)
+#' }
+#'
+#' @details
+#' The rule of thumb is that the models are "very similar" if |elpd_diff| (the
+#' absolute value of elpd_diff) is less than 4 (Sivula, Magnusson and Vehtari, 2020).
+#' If superior to 4, then one can use the SE to obtain a standardized difference
+#' (Z-diff) and interpret it as such, assuming that the difference is normally
+#' distributed.
+#'
+#' @return Objects of class [report_text()].
+#' @export
+report.compare.loo <- function(x, ...) {
+  # https://stats.stackexchange.com/questions/608881/how-to-interpret-elpd-diff-of-bayesian-loo-estimate-in-bayesian-logistic-regress
+  # https://users.aalto.fi/%7Eave/CV-FAQ.html#12_What_is_the_interpretation_of_ELPD__elpd_loo__elpd_diff
+  # https://users.aalto.fi/%7Eave/CV-FAQ.html#se_diff
+
+  # The difference in expected log predictive density (elpd) between each model
+  # and the best model as well as the standard error of this difference (assuming
+  # the difference is approximately normal).
+  x <- as.data.frame(x)
+  # The values in the first row are 0s because the models are ordered from best to worst according to their elpd.
+  modnames <- rownames(x)
+
+  #
+  z_elpd_diff <- x$elpd_diff / x$se_diff
+
+  text <- "The difference in predictive accuracy, as index by Expected Log Predictive Density (ELPD), suggests that '"
+  text <- paste0(text, modnames[1], "' is the best model, followed by '")
+
+  for (m in 2:(nrow(x))) {
+    text <- paste0(text, modnames[m ],
+                   "' (diff = ",
+                   insight::format_value(x$elpd_diff[m]),
+                   # " ± ",
+                   # insight::format_value(x$se_diff[m]),
+                   ", Z-diff = ",
+                   insight::format_value(z_elpd_diff[m]),
+                   ")")
+  }
+  class(text) <- c("report_text", class(text))
+  text
+}
\ No newline at end of file
diff --git a/man/report.compare.loo.Rd b/man/report.compare.loo.Rd
new file mode 100644
index 00000000..e03ef453
--- /dev/null
+++ b/man/report.compare.loo.Rd
@@ -0,0 +1,39 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/report.compare.loo.R
+\name{report.compare.loo}
+\alias{report.compare.loo}
+\title{Reporting Bayesian Model Comparison}
+\usage{
+\method{report}{compare.loo}(x)
+}
+\arguments{
+\item{x}{An object of class \code{\link[=loo_compare]{loo_compare()}}.}
+}
+\value{
+Objects of class \code{\link[=report_text]{report_text()}}.
+}
+\description{
+Automatically report the results of Bayesian model comparison using the \code{loo} package.
+}
+\details{
+The rule of thumb is that the models are "very similar" if |elpd_diff| (the
+absolute value of elpd_diff) is less than 4 (Sivula, Magnusson and Vehtari, 2020).
+If superior to 4, then one can use the SE to obtain a standardized difference
+(Z-diff) and interpret it as such, assuming that the difference is normally
+distributed.
+}
+\examples{
+\dontshow{if (require("brms", quietly = TRUE)) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf}
+\donttest{
+library(brms)
+
+m1 <- brms::brm(mpg ~ qsec, data=mtcars)
+m2 <- brms::brm(mpg ~ qsec + wt, data=mtcars)
+
+x <- brms::loo_compare(brms::add_criterion(m1, "loo"),
+                       brms::add_criterion(m2, "loo"),
+                       model_names=c("m1", "m2"))
+report(x)
+}
+\dontshow{\}) # examplesIf}
+}