-
-
Notifications
You must be signed in to change notification settings - Fork 70
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge branch 'main' into joss-paper-draft
- Loading branch information
Showing
15 changed files
with
316 additions
and
220 deletions.
There are no files selected for viewing
This file was deleted.
Oops, something went wrong.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,87 @@ | ||
| #' Reporting `BFBayesFactor` objects from the `BayesFactor` package | ||
| #' | ||
| #' Interpretation of the Bayes factor output from the `BayesFactor` package. | ||
| #' | ||
| #' @param x An object of class `BFBayesFactor`. | ||
| #' @param h0,h1 Names of the null and alternative hypotheses. | ||
| #' @param table A `parameters` table (this argument is meant for internal use). | ||
| #' @param ... Other arguments to be passed to [effectsize::interpret_bf] and [insight::format_bf]. | ||
| #' | ||
| #' @examplesIf requireNamespace("BayesFactor", quietly = TRUE) | ||
| #' \donttest{ | ||
| #' library(BayesFactor) | ||
| #' | ||
| #' rez <- BayesFactor::ttestBF(iris$Sepal.Width, iris$Sepal.Length) | ||
| #' report_statistics(rez, exact = TRUE) # Print exact BF | ||
| #' report(rez, h0 = "the null hypothesis", h1 = "the alternative") | ||
| #' | ||
| #' rez <- BayesFactor::correlationBF(iris$Sepal.Width, iris$Sepal.Length) | ||
| #' report(rez) | ||
| #' } | ||
| #' | ||
| #' @export | ||
| report.BFBayesFactor <- function(x, h0 = "H0", h1 = "H1", ...) { | ||
| if (inherits("BFlinearModel", class(x@numerator[[1]]))) { | ||
| return(report(bayestestR::bayesfactor_models(x), ...)) | ||
| } | ||
|
|
||
| if (length(x@numerator) > 1) { | ||
| insight::format_alert( | ||
| "Multiple `BFBayesFactor` models detected - reporting for the first numerator model.", | ||
| "See help(\"get_parameters\", package = \"insight\")." | ||
| ) | ||
| x <- x[1] | ||
| } | ||
|
|
||
| param <- parameters::parameters(x[1], ...) | ||
| bf <- param$BF | ||
| other_dir <- ifelse(bf < 1, "h0", "h1") | ||
|
|
||
|
|
||
| if (other_dir == "h1") { | ||
| other_text <- paste0( | ||
| "There is ", | ||
| effectsize::interpret_bf(bf, ...), | ||
| " ", | ||
| h1, | ||
| " over ", | ||
| h0, | ||
| " (", report_statistics(x, ...), ")." | ||
| ) | ||
| } else { | ||
| other_text <- paste0( | ||
| "There is ", | ||
| effectsize::interpret_bf(1 / bf, ...), | ||
| " ", | ||
| h0, | ||
| " over ", | ||
| h1, | ||
| " (", report_statistics(x, ...), ")." | ||
| ) | ||
| } | ||
| other_text | ||
| } | ||
|
|
||
|
|
||
|
|
||
| #' @rdname report.BFBayesFactor | ||
| #' @export | ||
| report_statistics.BFBayesFactor <- function(x, table = NULL, ...) { | ||
| if (is.null(table)) { | ||
| if (length(x@numerator) > 1) { | ||
| insight::format_alert( | ||
| "Multiple `BFBayesFactor` models detected - reporting for the first numerator model.", | ||
| "See help(\"get_parameters\", package = \"insight\")." | ||
| ) | ||
| x <- x[1] | ||
| } | ||
| table <- parameters::parameters(x, ...) | ||
| } | ||
|
|
||
| bf <- table$BF | ||
| other_text <- ifelse(bf < 1, | ||
| insight::format_bf(1 / bf, name = "BF01", ...), | ||
| insight::format_bf(bf, name = "BF10", ...) | ||
| ) | ||
| other_text | ||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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 | ||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.
Oops, something went wrong.
Oops, something went wrong.