properly compute BF_M for non uniform priors

R/commonAnovaBayesian.R

@@ -782,26 +782,36 @@ BANOVAcomputMatchedInclusion <- function(effectNames, effects.matrix, interactio
     # no errors, proceed normally and complete the table
 
     logsumbfs <- logSumExp(logbfs + logprior)
-    internalTable[["P(M|data)"]] <-  exp(logbfs + logprior - logsumbfs)
+    logPostProbModel <- logbfs + logprior - logsumbfs
+    internalTable[["P(M|data)"]] <-  exp(logPostProbModel)
 
-    nmodels <- nrow(internalTable)
-    for (i in seq_len(nmodels)) {
-      internalTable[i, "BFM"] <- logbfs[i] - logSumExp(logbfs[-i]) + log(nmodels - 1L)
-    }
+    logNumPostOdds  <- logPostProbModel
+    logDenPostOdds  <- log1mexp(logNumPostOdds)
+    logNumPriorOdds <- logprior
+    logDenPriorOdds <- log1mexp(logNumPriorOdds)
+    internalTable[["BFM"]] <- logNumPostOdds - logDenPostOdds + logDenPriorOdds - logNumPriorOdds
 
   } else {
     # create table excluding failed models
 
     idxGood <- !is.na(logbfs)
     widxGood <- which(idxGood)
-    logsumbfs <- logSumExp(logbfs[idxGood])
-    internalTable[["P(M|data)"]] <-  exp(logbfs - logsumbfs)
 
-    nmodels <- sum(idxGood)
+    # normalize the prior w.r.t the non-failed models
+    logpriorSubset <- logprior[idxGood]
+    logpriorSubset <- logpriorSubset - logSumExp(logpriorSubset)
+
+    logsumbfs <- logSumExp(logbfs[idxGood] + logpriorSubset)
+    logPostProbModel <- logbfs[idxGood] + logpriorSubset - logsumbfs
+    internalTable[widxGood, "P(M|data)"] <- exp(logPostProbModel)
+
+    logNumPostOdds  <- logPostProbModel
+    logDenPostOdds  <- log1mexp(logNumPostOdds)
+    logNumPriorOdds <- logpriorSubset
+    logDenPriorOdds <- log1mexp(logNumPriorOdds)
+    internalTable[widxGood, "BFM"] <- logNumPostOdds - logDenPostOdds + logDenPriorOdds - logNumPriorOdds
+
     widxBad <- which(!idxGood)
-    for (i in widxGood) {
-      internalTable[["BFM"]][i] <- logbfs[i] - logSumExp(logbfs[-c(i, widxBad)]) + log(nmodels - 1L)
-    }
 
     internalTable[widxBad, "P(M|data)"] <- NaN
     internalTable[widxBad, "BFM"]       <- NaN
@@ -2832,6 +2842,17 @@ BANOVAcomputMatchedInclusion <- function(effectNames, effects.matrix, interactio
   return(interactions.matrix)
 }
 
+#' Accurately compute log(1 - exp(x))
+#'
+#' @param x numeric value or vector
+#'
+#' @details See https://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf
+#'
+log1mexp <- function(x) {
+  a0 <- -0.69314718055994528623 # log(1 / 2)
+  ifelse(x < a0, log1p(-exp(x)), log(-expm1(x)))
+}
+
 # Model prior ----
 .BANOVAcomputePriorModelProbs <- function(models, nuisance, options) {
 

tests/testthat/test-anovabayesian.R

@@ -223,6 +223,12 @@ test_that("Model prior changes posterior model probabilities", {
     } else {
       jaspTools::expect_equal_tables(tableModelComparison, reference[[modelPrior]][["modelComparison"]])
       jaspTools::expect_equal_tables(tableEffects,         reference[[modelPrior]][["posteriorSummary"]])
+
+      # internal consistency check
+      BFM              <- vapply(tableModelComparison, `[[`, "BFM",       FUN.VALUE = numeric(1))
+      prior_model_prob <- vapply(tableModelComparison, `[[`, "P(M)",      FUN.VALUE = numeric(1))
+      post_model_prob  <- vapply(tableModelComparison, `[[`, "P(M|data)", FUN.VALUE = numeric(1))
+      testthat::expect_equal(BFM, post_model_prob / (1 - post_model_prob) * (1 - prior_model_prob) / prior_model_prob)
     }
   }
   if (createReference) saveRDS(reference, referencePath)