This (mini) R package allows for further computation of BFs based on
models built with BayesFactor
(but not only).
You can install BFEffects
from
github with:
# install.packages("devtools")
devtools::install_github("mattansb/BFEffects")
The following example is a reproduction of the example given here by Richard Morey (a must read to truly understand restricted model BFs).
library(BayesFactor)
# Load the data
disgust_data <- read.table(url('http://www.learnbayes.org/disgust_example.txt'),header=TRUE)
# Compute Model
BF <- anovaBF(score ~ condition, data = disgust_data, progress = FALSE)
# Sample from the posterior distribution
samples <- posterior(BF, iterations = 10000, progress = FALSE)
# Check order constraint
consistent <-
(samples[, "condition-control"] > samples[, "condition-lemon"]) &
(samples[, "condition-sulfur"] > samples[, "condition-control"])
N_consistent <- sum(consistent)
# Compute the Bayes factor of the restriction to the full model
bf_restriction_against_full <- (N_consistent / 10000) / (1 / 6)
# Compute the Bayes factor of our restriction against the null hypothesis
bf_full_against_null <- unname(as.vector(BF))
bf_restriction_against_null <- bf_restriction_against_full * bf_full_against_null
After all that, we get:
c(BF_res.full = bf_restriction_against_full,
BF_full.0 = bf_full_against_null,
BF_res.0 = bf_restriction_against_null)
## BF_res.full BF_full.0 BF_res.0
## 4.3884000 0.7738164 3.3958161
restrict
… or, we could use restrict
!
library(BayesFactor)
library(BFEffects)
# Load the data
disgust_data <- read.table(url('http://www.learnbayes.org/disgust_example.txt'),header=TRUE)
# Compute Model
BF <- anovaBF(score ~ condition, data = disgust_data, progress = FALSE)
# Restrict the model
BF.r <- restrict(BF, lemon < control & control < sulfur)
BF.r
## Restricted Model Probabilities:
## Prior Posterior
## Probability: 0.166 0.728
##
##
## BFs:
## Restricted Full Null
## Restricted 1.000 0.228 0.294
## Full 4.391 1.000 1.292
## Null 3.398 0.774 1.000
## ---
## columns are the numerator and rows are the denominator
As we can see, we’ve reproduced the results.
This function has been moved to
bayestestR::bayesfactor_inclusion
This is an R version of JASP’s output >
effects
option:
library(BayesFactor)
library(BFEffects)
BF <- generalTestBF(len ~ supp * dose, ToothGrowth, progress = FALSE)
BF
## Bayes factor analysis
## --------------
## [1] supp : 1.198757 ±0.01%
## [2] dose : 353886181229 ±0.01%
## [3] supp + dose : 1.536665e+13 ±45.97%
## [4] supp + dose + supp:dose : 1.086319e+13 ±3.11%
##
## Against denominator:
## Intercept only
## ---
## Bayes factor type: BFlinearModel, JZS
inclusionBF(BF)
## Using 'bayestestR::bayesfactor_inclusion'
## Pr(prior) Pr(posterior) Inclusion.BF
## supp 0.6 0.987 4.94130e+01
## dose 0.6 1.000 8.06023e+12
## supp:dose 0.2 0.409 2.76400e+00
## ---
## Inclusion BFs compared among all models.
This function has been moved to
bayestestR::bayesfactor_models
If you wish to compute BFs from non-BayesFactor
models, you can do so
with the following functions:
BIC_BFs
- compute BF based on BIC approcimations (see more
here).stan_BFs
- compute BF from models fitted with brms
and
rstanarm
packages via marginal likelihoods (using the
bridgesampling
package).These functions produce a BFGrid
object that prints nicely:
library(BFEffects)
mo0 <- lm(Sepal.Length ~ 1, data = iris)
mo1 <- lm(Sepal.Length ~ Species, data = iris)
mo2 <- lm(Sepal.Length ~ Species + Petal.Length, data = iris)
mo3 <- lm(Sepal.Length ~ Species * Petal.Length, data = iris)
BIC_BFs(mo0,mo1,mo2,mo3, .den = 1) # .den determines which of the models to compare to
## Using 'bayestestR::bayesfactor_models'
## Bayes factor analysis
## --------------
##
## [2] Species 1.695852e+29
## [3] Species + Petal.Length 5.843105e+55
## [4] Species * Petal.Length 2.203243e+54
##
## Against denominator:
## [1] 1
## ---
## Bayes factor type: BIC approximation
BayesFactor
thingsinferBF
?Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.