diff --git a/R/commonAnovaBayesian.R b/R/commonAnovaBayesian.R
index 1550d660..095654c1 100644
--- a/R/commonAnovaBayesian.R
+++ b/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) {
diff --git a/tests/testthat/BANOVA_modelPriors.rds b/tests/testthat/BANOVA_modelPriors.rds
index 76cc934c..6531b50b 100644
Binary files a/tests/testthat/BANOVA_modelPriors.rds and b/tests/testthat/BANOVA_modelPriors.rds differ
diff --git a/tests/testthat/_snaps/anova/raincloud-plots-horizontal.svg b/tests/testthat/_snaps/anova/raincloud-plots-horizontal.svg
index 2206b22b..21990290 100644
--- a/tests/testthat/_snaps/anova/raincloud-plots-horizontal.svg
+++ b/tests/testthat/_snaps/anova/raincloud-plots-horizontal.svg
@@ -21,77 +21,77 @@
-
-
+
+
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-1
-0
-
-
-
-
-
-
-
-
-
-
-
--4
--3
--2
--1
-0
-1
-2
-3
-contNormal
-contBinom
-raincloud-plots-horizontal
+
+1
+0
+
+
+
+
+
+
+
+
+
+
+
+-4
+-3
+-2
+-1
+0
+1
+2
+3
+contNormal
+contBinom
+raincloud-plots-horizontal
diff --git a/tests/testthat/_snaps/anova/raincloud-plots-vertical.svg b/tests/testthat/_snaps/anova/raincloud-plots-vertical.svg
index 63f66743..679848cd 100644
--- a/tests/testthat/_snaps/anova/raincloud-plots-vertical.svg
+++ b/tests/testthat/_snaps/anova/raincloud-plots-vertical.svg
@@ -21,77 +21,77 @@
-
-
+
+
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
--4
--3
--2
--1
-0
-1
-2
-3
-
-
-
-
-
-
-
-
-
-
-
-1
-0
-contBinom
-contNormal
-raincloud-plots-vertical
+
+-4
+-3
+-2
+-1
+0
+1
+2
+3
+
+
+
+
+
+
+
+
+
+
+
+1
+0
+contBinom
+contNormal
+raincloud-plots-vertical
diff --git a/tests/testthat/_snaps/anovarepeatedmeasures/raincloud-plots-between-subjects.svg b/tests/testthat/_snaps/anovarepeatedmeasures/raincloud-plots-between-subjects.svg
index 594103ac..271abe51 100644
--- a/tests/testthat/_snaps/anovarepeatedmeasures/raincloud-plots-between-subjects.svg
+++ b/tests/testthat/_snaps/anovarepeatedmeasures/raincloud-plots-between-subjects.svg
@@ -21,114 +21,114 @@
-
-
+
+
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-20
-40
-60
-80
-100
-120
-140
-
-
-
-
-
-
-
-
-
-
-Male
-Female
-gender
+
+20
+40
+60
+80
+100
+120
+140
+
+
+
+
+
+
+
+
+
+
+Male
+Female
+gender
raincloud-plots-between-subjects
diff --git a/tests/testthat/_snaps/anovarepeatedmeasures/raincloud-plots.svg b/tests/testthat/_snaps/anovarepeatedmeasures/raincloud-plots.svg
index ff198058..3ab58155 100644
--- a/tests/testthat/_snaps/anovarepeatedmeasures/raincloud-plots.svg
+++ b/tests/testthat/_snaps/anovarepeatedmeasures/raincloud-plots.svg
@@ -21,182 +21,182 @@
-
-
+
+
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-20
-40
-60
-80
-100
-120
-
-
-
-
-
-
-
-
-
-
-High
-Some
-None
-Charisma
+
+20
+40
+60
+80
+100
+120
+
+
+
+
+
+
+
+
+
+
+High
+Some
+None
+Charisma
raincloud-plots
diff --git a/tests/testthat/test-anovabayesian.R b/tests/testthat/test-anovabayesian.R
index d76fdc25..325a76bf 100644
--- a/tests/testthat/test-anovabayesian.R
+++ b/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)