From 44c92ad16bb756d6bf50d4cd6bad4a8bdc536ad7 Mon Sep 17 00:00:00 2001 From: Dominique Makowski Date: Mon, 1 Apr 2024 11:33:08 +0100 Subject: [PATCH] [Feature] Add report.compare.loo (#419) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Add report.compare.loo * more option for loo-comapre * fixes * Update WORDLIST * fix build_reference_index, styler * fix lints * styler * bump version * update readme, and #no lint start for long URL lint * styler # no lint start with space * # nolint start remove space..... --------- Co-authored-by: Mattan S. Ben-Shachar Co-authored-by: Rémi Thériault <13123390+rempsyc@users.noreply.github.com> --- DESCRIPTION | 3 +- NAMESPACE | 1 + NEWS.md | 6 +++ R/report.compare.loo.R | 111 ++++++++++++++++++++++++++++++++++++++ _pkgdown.yml | 1 + inst/WORDLIST | 8 ++- man/report.compare.loo.Rd | 45 ++++++++++++++++ 7 files changed, 173 insertions(+), 2 deletions(-) create mode 100644 R/report.compare.loo.R create mode 100644 man/report.compare.loo.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 35d6c831..96dfeec9 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: report Type: Package Title: Automated Reporting of Results and Statistical Models -Version: 0.5.8.1 +Version: 0.5.8.2 Authors@R: c(person(given = "Dominique", family = "Makowski", @@ -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/NEWS.md b/NEWS.md index 4a52b5b5..00b229cb 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,9 @@ +# report 0.5.9 + +Minor changes + +* `report` now supports reporting of Bayesian model comparison with variables of class `brms::loo_compare`. + # report 0.5.8 New features diff --git a/R/report.compare.loo.R b/R/report.compare.loo.R new file mode 100644 index 00000000..7225dae4 --- /dev/null +++ b/R/report.compare.loo.R @@ -0,0 +1,111 @@ +#' Reporting Bayesian Model Comparison +#' +#' Automatically report the results of Bayesian model comparison using the `loo` package. +#' +#' @param x An object of class [brms::loo_compare]. +#' @param index type if index to report - expected log pointwise predictive +#' density (ELPD) or information criteria (IC). +#' @param ... Additional arguments (not used for now). +#' +#' @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, index = c("ELPD", "IC"), ...) { + # nolint start + # https://stats.stackexchange.com/questions/608881/how-to-interpret-elpd-diff-of-bayesian-loo-estimate-in-bayesian-logistic-regress + # nolint end + # 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). + index <- match.arg(index) + 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) + + elpd_diff <- x[["elpd_diff"]] + ic_diff <- -2 * elpd_diff + + z_elpd_diff <- elpd_diff / x[["se_diff"]] + z_ic_diff <- -z_elpd_diff + + if ("looic" %in% colnames(x)) { + type <- "LOO" + ENP <- x[["p_loo"]] + } else { + type <- "WAIC" + ENP <- x[["p_waic"]] + } + + if (index == "ELPD") { + index_label <- sprintf("Expected Log Predictive Density (ELPD-%s)", type) + } else if (type == "LOO") { + index_label <- "Leave-One-Out CV Information Criterion (LOOIC)" + } else { + index_label <- "Widely Applicable Information Criterion (WAIC)" + } + + out_text <- sprintf( + paste( + "The difference in predictive accuracy, as index by %s, suggests that '%s' ", + "is the best model (effective number of parameters (ENP) = %.2f), followed by" + ), + index_label, modnames[1], ENP[1] + ) + + if (index == "ELPD") { + other_texts <- sprintf( + "'%s' (diff = %.2f, ENP = %.2f, z-diff = %.2f)", + modnames[-1], + elpd_diff[-1], + ENP[-1], + z_elpd_diff[-1] + ) + } else { + other_texts <- sprintf( + "'%s' (diff = %.2f, ENP = %.2f, z-diff = %.2f)", + modnames[-1], + ic_diff[-1], + ENP[-1], + z_ic_diff[-1] + ) + } + + sep <- "." + nothermods <- length(other_texts) + if (nothermods > 1L) { + if (nothermods == 2L) { + sep <- c(" and ", sep) + } else { + sep <- c(rep(", ", length = nothermods - 2), ", and ", sep) + } + } + + other_texts <- paste0(other_texts, sep, collapse = "") + + out_text <- paste(out_text, other_texts, collapse = "") + class(text) <- c("report_text", class(text)) + out_text +} diff --git a/_pkgdown.yml b/_pkgdown.yml index e8dc6725..a2d84eec 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -91,6 +91,7 @@ reference: - report.stanreg - report.test_performance - report.estimate_contrasts + - report.compare.loo - title: Report Non-Statistical Objects desc: | diff --git a/inst/WORDLIST b/inst/WORDLIST index 055ef665..af7640ca 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -4,9 +4,12 @@ BibLaTeX CMD CSL Dom +ELPD ESS Gabry Hotfix +IC +Magnusson Mattan Newcombe ORCID @@ -16,16 +19,20 @@ Rhat SEM SEXIT Shachar +Sivula +Vehtari amongst anova bmwiernik brms eXistence easystats +elpd github htest https ivreg +lifecycle mattansb pacakges participants’ @@ -42,4 +49,3 @@ unarchive versicolor virginica ’s -lifecycle diff --git a/man/report.compare.loo.Rd b/man/report.compare.loo.Rd new file mode 100644 index 00000000..47c2ff3b --- /dev/null +++ b/man/report.compare.loo.Rd @@ -0,0 +1,45 @@ +% 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, index = c("ELPD", "IC"), ...) +} +\arguments{ +\item{x}{An object of class \link[brms:loo_compare.brmsfit]{brms::loo_compare}.} + +\item{index}{type if index to report - expected log pointwise predictive +density (ELPD) or information criteria (IC).} + +\item{...}{Additional arguments (not used for now).} +} +\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 + drat, 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} +}