use new dlogspline

DESCRIPTION

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

R/bayesfactor_models.R

@@ -177,7 +177,7 @@
   cl$...$estimator <- cl$...$check_response <- NULL
 
 
   names(mods) <- sapply(cl$`...`, insight::safe_deparse)
   names(denominator) <- insight::safe_deparse(cl$denominator)
 
   mods <- .cleanup_BF_models(mods, denominator, cl)
@@ -196,7 +196,7 @@
     supported_models[!has_terms] <- FALSE
   }
 
   objects <- .safe(do.call(insight::ellipsis_info, c(mods, verbose = FALSE)))
   if (!is.null(objects)) {
     were_checked <- inherits(objects, "ListModels")
 
@@ -333,33 +333,20 @@
 
 #' @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)
   names(mods) <- sapply(cl$`...`, insight::safe_deparse)
   names(denominator) <- insight::safe_deparse(cl$denominator)
 
   mods <- .cleanup_BF_models(mods, denominator, cl)
@@ -370,22 +357,11 @@
 
 
 #' @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
@@ -397,7 +373,7 @@
   mBFs <- c(0, BayesFactor::extractBF(models, TRUE, TRUE))
   mforms <- sapply(c(models@denominator, models@numerator), function(x) x@shortName)
 
   if (!inherits(models@denominator, "BFlinearModel")) {
     mforms <- .clean_non_linBF_mods(mforms)
   } else {
     mforms[mforms == "Intercept only"] <- "1"
@@ -458,7 +434,7 @@
   out <- -outer(x$log_BF, x$log_BF, FUN = "-")
   rownames(out) <- colnames(out) <- x$Model
 
   # out <- exp(out)
 
   class(out) <- c("bayesfactor_models_matrix", class(out))
   out
@@ -471,8 +447,8 @@
   if (length(mods) == 1 && inherits(mods[[1]], "list")) {
     mods <- mods[[1]]
     mod_names <- tryCatch(
       {
         sapply(cl$`...`[[1]][-1], insight::safe_deparse)
       },
       error = function(e) {
         NULL
@@ -483,7 +459,7 @@
     }
   }
 
   if (!is.numeric(denominator[[1]])) {
     denominator_model <- which(names(mods) == names(denominator))
 
     if (length(denominator_model) == 0) {
@@ -591,7 +567,7 @@
       return(m_txt)
     },
     error = function(e) {
       return(m_names)
     }
   )
 }

R/bayesfactor_parameters.R

@@ -262,7 +262,7 @@
     prior <- posterior
     if (verbose) {
       insight::format_warning(
         "Prior not specified! Please specify a prior (in the form 'prior = distribution_normal(1000, 0, 1)') to get meaningful results."
       )
     }
   }
@@ -417,9 +417,9 @@
   }
 
 
-  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,
@@ -430,7 +430,7 @@
 
   bf_val <- data.frame(
     Parameter = colnames(posterior),
-    log_BF = log(sdbf),
+    log_BF = sdlogbf,
     stringsAsFactors = FALSE
   )
 
@@ -470,23 +470,23 @@
 
 
 #' @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)
@@ -496,36 +496,29 @@
         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")