From da69f9e1cff1aea161a685fb1780698f4bbf13c6 Mon Sep 17 00:00:00 2001 From: "Mattan S. Ben-Shachar" Date: Thu, 26 Oct 2023 23:30:15 +0300 Subject: [PATCH] use new dlogspline fix #624 --- DESCRIPTION | 2 +- R/bayesfactor_models.R | 44 ++++----------- R/bayesfactor_parameters.R | 57 +++++++++----------- R/bayesfactor_restricted.R | 4 +- tests/testthat/test-bayesfactor_parameters.R | 15 +++--- 5 files changed, 44 insertions(+), 78 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index ae192a25e..4fe4c524c 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -92,7 +92,7 @@ Suggests: knitr, lavaan, lme4, - logspline, + logspline (>= 2.1.21), MASS, mclust, mediation, diff --git a/R/bayesfactor_models.R b/R/bayesfactor_models.R index 1991ee0ca..1f6ad7079 100644 --- a/R/bayesfactor_models.R +++ b/R/bayesfactor_models.R @@ -333,29 +333,16 @@ bayesfactor_models.default <- function(..., denominator = 1, verbose = TRUE) { #' @export bayesfactor_models.stanreg <- function(..., denominator = 1, verbose = TRUE) { - insight::check_if_installed("rstanarm") - - # Organize the models and their names mods <- list(...) - denominator <- list(denominator) - - cl <- match.call(expand.dots = FALSE) - names(mods) <- sapply(cl$`...`, insight::safe_deparse) - names(denominator) <- insight::safe_deparse(cl$denominator) - - mods <- .cleanup_BF_models(mods, denominator, cl) - denominator <- attr(mods, "denominator", exact = TRUE) - - .bayesfactor_models_stan(mods, denominator = denominator, verbose = verbose) -} - - -#' @export -bayesfactor_models.brmsfit <- function(..., denominator = 1, verbose = TRUE) { - insight::check_if_installed("brms") + if (inherits(mods[[1]], "stanreg")) { + insight::check_if_installed("rstanarm") + } else if (inherits(mods[[1]], "brmsfit")) { + insight::check_if_installed("brms") + } else if (inherits(mods[[1]], "blavaan")) { + insight::check_if_installed("blavaan") + } # Organize the models and their names - mods <- list(...) denominator <- list(denominator) cl <- match.call(expand.dots = FALSE) @@ -370,22 +357,11 @@ bayesfactor_models.brmsfit <- function(..., denominator = 1, verbose = TRUE) { #' @export -bayesfactor_models.blavaan <- function(..., denominator = 1, verbose = TRUE) { - insight::check_if_installed("blavaan") +bayesfactor_models.brmsfit <- bayesfactor_models.stanreg - # Organize the models and their names - mods <- list(...) - denominator <- list(denominator) - cl <- match.call(expand.dots = FALSE) - names(mods) <- sapply(cl$`...`, insight::safe_deparse) - names(denominator) <- insight::safe_deparse(cl$denominator) - - mods <- .cleanup_BF_models(mods, denominator, cl) - denominator <- attr(mods, "denominator", exact = TRUE) - - .bayesfactor_models_stan(mods, denominator = denominator, verbose = verbose) -} +#' @export +bayesfactor_models.blavaan <- bayesfactor_models.stanreg #' @export diff --git a/R/bayesfactor_parameters.R b/R/bayesfactor_parameters.R index 77d98db70..0aa82978a 100644 --- a/R/bayesfactor_parameters.R +++ b/R/bayesfactor_parameters.R @@ -418,9 +418,9 @@ bayesfactor_parameters.data.frame <- function(posterior, } - sdbf <- numeric(ncol(posterior)) + sdlogbf <- numeric(ncol(posterior)) for (par in seq_along(posterior)) { - sdbf[par] <- .bayesfactor_parameters( + sdlogbf[par] <- .logbayesfactor_parameters( posterior[[par]], prior[[par]], direction = direction, @@ -431,7 +431,7 @@ bayesfactor_parameters.data.frame <- function(posterior, bf_val <- data.frame( Parameter = colnames(posterior), - log_BF = log(sdbf), + log_BF = sdlogbf, stringsAsFactors = FALSE ) @@ -471,23 +471,23 @@ bayesfactor_parameters.rvar <- bayesfactor_parameters.draws #' @keywords internal -.bayesfactor_parameters <- function(posterior, - prior, - direction = 0, - null = 0, - ...) { +.logbayesfactor_parameters <- function(posterior, + prior, + direction = 0, + null = 0, + ...) { stopifnot(length(null) %in% c(1, 2)) if (isTRUE(all.equal(posterior, prior))) { - return(1) + return(0) } insight::check_if_installed("logspline") if (length(null) == 1) { - relative_density <- function(samples) { + relative_loglikelihood <- function(samples) { f_samples <- .logspline(samples, ...) - d_samples <- logspline::dlogspline(null, f_samples) + d_samples <- logspline::dlogspline(null, f_samples, log = TRUE) if (direction < 0) { norm_samples <- logspline::plogspline(null, f_samples) @@ -497,36 +497,29 @@ bayesfactor_parameters.rvar <- bayesfactor_parameters.draws norm_samples <- 1 } - d_samples / norm_samples + d_samples - log(norm_samples) } - - return(relative_density(prior) / relative_density(posterior)) } else if (length(null) == 2) { null <- sort(null) null[is.infinite(null)] <- 1.797693e+308 * sign(null[is.infinite(null)]) - f_prior <- .logspline(prior, ...) - f_posterior <- .logspline(posterior, ...) - - h0_prior <- diff(logspline::plogspline(null, f_prior)) - h0_post <- diff(logspline::plogspline(null, f_posterior)) + relative_loglikelihood <- function(samples) { + f_samples <- .logspline(samples, ...) + p_samples <- diff(logspline::plogspline(null, f_samples)) - BF_null_full <- h0_post / h0_prior + if (direction < 0) { + norm_samples <- logspline::plogspline(min(null), f_samples) + } else if (direction > 0) { + norm_samples <- 1 - logspline::plogspline(max(null), f_samples) + } else { + norm_samples <- 1 - p_samples + } - if (direction < 0) { - h1_prior <- logspline::plogspline(min(null), f_prior) - h1_post <- logspline::plogspline(min(null), f_posterior) - } else if (direction > 0) { - h1_prior <- 1 - logspline::plogspline(max(null), f_prior) - h1_post <- 1 - logspline::plogspline(max(null), f_posterior) - } else { - h1_prior <- 1 - h0_prior - h1_post <- 1 - h0_post + log(p_samples) - log(norm_samples) } - BF_alt_full <- h1_post / h1_prior - - return(BF_alt_full / BF_null_full) } + + relative_loglikelihood(prior) - relative_loglikelihood(posterior) } # Bad Methods ------------------------------------------------------------- diff --git a/R/bayesfactor_restricted.R b/R/bayesfactor_restricted.R index 454d9adff..699a8d522 100644 --- a/R/bayesfactor_restricted.R +++ b/R/bayesfactor_restricted.R @@ -221,13 +221,13 @@ bayesfactor_restricted.data.frame <- function(posterior, hypothesis, prior = NUL posterior_p <- sapply(posterior_l, mean) prior_p <- sapply(prior_l, mean) - BF <- posterior_p / prior_p + log_BF <- log(posterior_p) - log(prior_p) res <- data.frame( Hypothesis = hypothesis, p_prior = prior_p, p_posterior = posterior_p, - log_BF = log(BF) + log_BF = log_BF ) attr(res, "bool_results") <- list(posterior = posterior_l, prior = prior_l) diff --git a/tests/testthat/test-bayesfactor_parameters.R b/tests/testthat/test-bayesfactor_parameters.R index 9556483b8..9a3538b79 100644 --- a/tests/testthat/test-bayesfactor_parameters.R +++ b/tests/testthat/test-bayesfactor_parameters.R @@ -1,9 +1,5 @@ test_that("bayesfactor_parameters data frame", { - skip_if_offline() - skip_if_not_or_load_if_installed("rstanarm") - skip_if_not_or_load_if_installed("BayesFactor") - skip_if_not_or_load_if_installed("httr") - skip_if_not_or_load_if_installed("brms") + skip_if_not_or_load_if_installed("logspline", "2.1.21") Xprior <- data.frame( x = distribution_normal(1e4), @@ -59,10 +55,9 @@ test_that("bayesfactor_parameters data frame", { test_that("bayesfactor_parameters RSTANARM", { skip_on_cran() skip_if_offline() + skip_if_not_or_load_if_installed("logspline", "2.1.21") skip_if_not_or_load_if_installed("rstanarm") - skip_if_not_or_load_if_installed("BayesFactor") skip_if_not_or_load_if_installed("httr") - skip_if_not_or_load_if_installed("brms") fit <- suppressMessages(stan_glm(mpg ~ ., data = mtcars, refresh = 0)) @@ -72,8 +67,11 @@ test_that("bayesfactor_parameters RSTANARM", { set.seed(333) BF1 <- bayesfactor_parameters(fit, verbose = FALSE) + BF3 <- bayesfactor_parameters(insight::get_parameters(fit), insight::get_parameters(fit_p), verbose = FALSE) expect_equal(BF1, BF2) + expect_equal(BF1[["Parameter"]], BF3[["Parameter"]]) + expect_equal(BF1[["log_BF"]], BF3[["log_BF"]]) model_flat <- suppressMessages( stan_glm(extra ~ group, data = sleep, prior = NULL, refresh = 0) @@ -94,8 +92,7 @@ test_that("bayesfactor_parameters RSTANARM", { test_that("bayesfactor_parameters BRMS", { skip_if_offline() - skip_if_not_or_load_if_installed("rstanarm") - skip_if_not_or_load_if_installed("BayesFactor") + skip_if_not_or_load_if_installed("logspline", "2.1.21") skip_if_not_or_load_if_installed("httr") skip_if_not_or_load_if_installed("brms") skip_if_not_or_load_if_installed("cmdstanr")