|
| 1 | +#' Bayesian Model Comparison |
| 2 | +#' |
| 3 | +#' Make a table of Bayesian model comparisons using the `loo` package. |
| 4 | +#' |
| 5 | +#' @param model An object of class [brms::loo_compare]. |
| 6 | +#' @param include_IC Whether to include the information criteria (IC). |
| 7 | +#' @param include_ENP Whether to include the effective number of parameters (ENP). |
| 8 | +#' @param ... Additional arguments (not used for now). |
| 9 | +#' |
| 10 | +# nolint start |
| 11 | +#' @examplesIf all(insight::check_if_installed(c("brms", "RcppEigen", "BH"), quietly = TRUE)) |
| 12 | +# nolint end |
| 13 | +#' \donttest{ |
| 14 | +#' library(brms) |
| 15 | +#' |
| 16 | +#' m1 <- brms::brm(mpg ~ qsec, data = mtcars) |
| 17 | +#' m2 <- brms::brm(mpg ~ qsec + drat, data = mtcars) |
| 18 | +#' m3 <- brms::brm(mpg ~ qsec + drat + wt, data = mtcars) |
| 19 | +#' |
| 20 | +#' x <- suppressWarnings(brms::loo_compare( |
| 21 | +#' brms::add_criterion(m1, "loo"), |
| 22 | +#' brms::add_criterion(m2, "loo"), |
| 23 | +#' brms::add_criterion(m3, "loo"), |
| 24 | +#' model_names = c("m1", "m2", "m3") |
| 25 | +#' )) |
| 26 | +#' model_parameters(x) |
| 27 | +#' model_parameters(x, include_IC = FALSE, include_ENP = TRUE) |
| 28 | +#' } |
| 29 | +#' |
| 30 | +#' @details |
| 31 | +#' The rule of thumb is that the models are "very similar" if |elpd_diff| (the |
| 32 | +#' absolute value of elpd_diff) is less than 4 (Sivula, Magnusson and Vehtari, 2020). |
| 33 | +#' If superior to 4, then one can use the SE to obtain a standardized difference |
| 34 | +#' (Z-diff) and interpret it as such, assuming that the difference is normally |
| 35 | +#' distributed. The corresponding p-value is then calculated as `2 * pnorm(-abs(Z-diff))`. |
| 36 | +#' However, note that if the raw ELPD difference is small (less than 4), it doesn't |
| 37 | +#' make much sense to rely on its standardized value: it is not very useful to |
| 38 | +#' conclude that a model is much better than another if both models make very |
| 39 | +#' similar predictions. |
| 40 | +#' |
| 41 | +#' @return Objects of `parameters_model`. |
| 42 | +#' @export |
| 43 | +model_parameters.compare.loo <- function(model, include_IC = TRUE, include_ENP = FALSE, ...) { |
| 44 | + # nolint start |
| 45 | + # https://stats.stackexchange.com/questions/608881/how-to-interpret-elpd-diff-of-bayesian-loo-estimate-in-bayesian-logistic-regress |
| 46 | + # nolint end |
| 47 | + # https://users.aalto.fi/%7Eave/CV-FAQ.html#12_What_is_the_interpretation_of_ELPD__elpd_loo__elpd_diff |
| 48 | + # https://users.aalto.fi/%7Eave/CV-FAQ.html#se_diff |
| 49 | + |
| 50 | + # The difference in expected log predictive density (elpd) between each model |
| 51 | + # and the best model as well as the standard error of this difference (assuming |
| 52 | + # the difference is approximately normal). |
| 53 | + |
| 54 | + # The values in the first row are 0s because the models are ordered from best to worst according to their elpd. |
| 55 | + x <- as.data.frame(model) |
| 56 | + |
| 57 | + out <- data.frame(Name = rownames(x), stringsAsFactors = FALSE) |
| 58 | + if ("looic" %in% colnames(x)) { |
| 59 | + if (include_IC) out$LOOIC <- x[["looic"]] |
| 60 | + if (include_ENP) out$ENP <- x[["p_loo"]] |
| 61 | + out$ELPD <- x[["elpd_loo"]] |
| 62 | + } else { |
| 63 | + if (include_IC) out$WAIC <- x[["waic"]] |
| 64 | + if (include_ENP) out$ENP <- x[["p_waic"]] |
| 65 | + out$ELPD <- x[["elpd_waic"]] |
| 66 | + } |
| 67 | + |
| 68 | + out$Difference <- x[["elpd_diff"]] |
| 69 | + out$Difference_SE <- x[["se_diff"]] |
| 70 | + |
| 71 | + z_elpd_diff <- x[["elpd_diff"]] / x[["se_diff"]] |
| 72 | + out$p <- 2 * stats::pnorm(-abs(z_elpd_diff)) |
| 73 | + |
| 74 | + class(out) <- c("parameters_model", "data.frame") |
| 75 | + out |
| 76 | +} |
0 commit comments