use new dlogspline

fix #624

DESCRIPTION

@@ -92,7 +92,7 @@ Suggests:
     knitr,
     lavaan,
     lme4,
-    logspline,
+    logspline (>= 2.1.21),
     MASS,
     mclust,
     mediation,

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

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 -------------------------------------------------------------

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)

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")