BFpack contains a collection of functions for Bayesian hypothesis testing using Bayes factors and posterior probabilities in R. The main function
BF needs a fitted model
x as input argument. Depending on the class of the fitted model, a standard hypothesis test is executed by default. For example, if
x is a fitted regression model of class
lm then posterior probabilities are computed of whether each separate coefficient is zero, negative, or positive (assuming equal prior probabilities). If one has specific hypotheses with equality and/or order constraints on the parameters under the fitted model
x then these can be formulated using the
hypothesis argument (a character string), possibly together prior probabilities for the hypotheses via the
prior.hyp argument (default all hypotheses are equally likely a priori), and the
complement argument which is a logical stating whether the complement hypotheses should be included in the case (
TRUE by default).
Alternatively, when the model of interest is not of a class that is currently supported,
x can also be a named numeric vector containing the estimates of the model parameters of interest, together with the error covariance matrix in the argument
Sigma, and the sample size used to obtain the estimates, to perform an approximate Bayes factor test using large sample theory.
The key references for the package are
Mulder, J., Williams, D. R., Gu, X., Tomarken, A., Boeing-Messing, F., Olsson-Collentine, A., Meijerink, M., Menke, J., van Aert, R., Fox, J.-P., Hoijtink, H., Rosseel, Y., Wagenmakers, E.-J., and van Lissa, C. (to appear). BFpack: Flexible Bayes Factor Testing of Scientific Theories in R. Journal of Statistical Software.
Mulder, J., van Lissa, C., Gu, X., Olsson-Collentine, A., Boeing-Messing, F., Williams, D. R., Fox, J.-P., Menke, J., et al. (2021). BFpack: Flexible Bayes Factor Testing of Scientific Expectations. (Version 0.3.2) https://CRAN.R-project.org/package=BFpack
BF(x, hypothesis, prior.hyp = NULL, complement = TRUE, ...)
x, a fitted model object that is obtained using a R-function. The object can be obtained via the following R functions:
t_testfor t testing,
bartlett_testfor testing independent group variances,
aovfor AN(C)OVA testing,
manovafor MAN(C)OVA testing,
lmfor linear regresssion analysis,
cor_testfor correlation analysis,
lmercurrently for testing intraclass correlations in random intercept models,
glmfor generalized linear models,
coxphfor survival analysis,
survregfor survival analysis,
polrfor ordinal regression,
zeroinflfor zero-inflated regression,
xcan also be a named vector with estimates of the key parameters.
hypothesis, a character string specifying the hypotheses with equality and/or order constraints on the key parameters of interest.
hypothesis = NULLwhich executes exploratory hypothesis tests (examples below).
model1is a fitted model object.
&. Hypotheses are separated using a semi-colon
;. For example
hypothesis = "weight > height & height > 0; weight = height = 0"implies that the first hypothesis assumes that the parameter
weightis larger than the parameter
heightand that the parameter
heightis positive, and the second hypothesis assumes that the two parameters are equal to zero. Note that the first hypothesis could equivalently have been written as
weight > height > 0.
prior.hyp, a numeric vector specifying the prior probabilities of the hypotheses of the
hypothesisargument. The default setting is
prior.hyp = NULLwhich sets equal prior probabilities.
complement, a logical value which specified if a complement hypothesis is included in the tested hypotheses specified under
hypothesis. The default setting is
TRUE. The complement hypothesis covers the remaining parameters space that is not covered by the constrained hypotheses. For example, if an equality hypothesis and an order hypothesis are formulated, say,
hypothesis = "weight = height = length; weight > height > length", the complement hypothesis covers the remaining subspace where neither
"weight = height = length"holds, nor
"weight > height > length"holds.
Alternatively if one is interesting in testing hypotheses on a model object which has a class that is currently not supported, an approximate Bayesian test can be executed with the following (additional) arguments
x, a named numeric vector of the estimates (e.g., MLE) of the parameters of interest where the labels are equal to the names of the parameters which are used for the
Sigma, the approximate posterior covariance matrix (e.g,. error covariance matrix) of the parameters of interest.
n, the sample size that was used to acquire the estimates and covariance matrix.
The output is of class
BF. By running the
BF object, a short overview of the results are presented. By running the
summary function on the
BF object, a comprehensive overview of the results are presented.
First a classical one sample t test is executed for the test value \(\mu = 5\) on the therapeutic data
<- bain::t_test(therapeutic, alternative = "greater", mu = 5)ttest1
t_test function is part of the bain package. The function is equivalent to the standard
t.test function with the addition that the returned object contains additional output than the standard
To see which parameters can be tested on this object run
which shows that the only parameter that can be tested is the population mean which has name
To perform an exploratory Bayesian t test of whether the population mean is equal to, smaller than, or larger than the null value (which is
5 here, as specified when defining the
ttest1 object), one needs to run
BF function on the object.
library(BFpack) <- BF(ttest1)BF1
This executes an exploratory (‘exhaustive’) test around the null value:
H1: mu = 5 versus
H2: mu < 5 versus
H3: mu > 5 assuming equal prior probabilities for
H3 of 1/3. The output presents the posterior probabilities for the three hypotheses.
The same test would be executed when the same hypotheses are explicitly specified using the
<- "mu = 5; mu < 5; mu > 5" hypothesis BF(ttest1, hypothesis = hypothesis)
In the above test the complement hypothesis is excluded automatically as the formualted hypothesis under the
hypothesis argument cover the complete parameter space. Furthermore, when testing hypotheses via the
hypothesis argument, the output also presents an
Evidence matrix containing the Bayes factors between the hypotheses formulated in the
A standard two-sided test around the null value
mu is executed by setting the hypothesis argument equal to the precise null hypothesis so that the complement hypothesis (which is included by default) corresponds to the hypothesis that assumes that the population mean is anything but the null value
<- "mu = 5" hypothesis BF(ttest1, hypothesis = hypothesis)
prior.hyp can be used to specify different prior probabilities for the hypotheses. For example, when the left one-tailed hypothesis is not possible based on prior considerations (e.g., see Mulder et al. (2021, Section 4.1, preprint)) while the precise (null) hypothesis and the right one-tailed hypothesis are equally likely, the argument
prior.hyp should be a vector specifying the prior probabilities of the respective hypotheses
BF(ttest1, hypothesis = "mu = 5; mu < 5; mu > 5", prior.hyp = c(.5,0,.5))
For more information about the methodology, we refer the interested reader to Mulder et al. (2021) and Mulder and Gu (2021).
First an analysis of variance (ANOVA) model is fitted using the
aov fuction in
<- aov(price ~ anchor * motivation, data = tvprices)aov1
Next a Bayesian test can be performed on the fitted object. By default exploratory tests are executed of whether the individual main and interaction effects are zero or not (corresponding to the full model) (see Mulder et al. (2021, Section 4.2))
One can also test for specific equal/order hypotheses based on scientific expectations such as whether
anchorrounded is positive,
motivationlow is negative, and the interaction effect
anchorrounded:motivationlow is negative (see Mulder et al. (2021, Section 4.2)) versus null hypothesis versus the complement hypothesis (which assumes that the constraints of neither two hypotheses hold). This test can be executed as follows:
<- "anchorrounded > 0 & motivationlow < 0 & constraints2 anchorrounded:motivationlow < 0; anchorrounded = 0 & motivationlow = 0 & anchorrounded:motivationlow = 0" set.seed(1234) BF(aov1, hypothesis = constraints2)
For more information about the methodology, we refer the interested reader to Mulder et al. (2021) and Mulder and Gu (2021).
First a classical significance test is executed using the
bartlett_test function, which is part of the BFpack package. The function is equivalent to the standard
bartlett.test function with the addition that the returned object contains additional output needed for the test using the
<- bartlett_test(x = attention$accuracy, g = attention$group)bartlett1
On an object of this class, by default
BF executes an exploratory test of homogeneity (equality) of variances against an unconstrained (full) model
The group variances have names
TS. This can be retrieved by running
Let’s say we want to test whether a hypothesis (H1) which assumes that group variances of groups
TS are equal and smaller than the group variance of the
ADHD group, a hypothesis (H2) which assumes that the group variances of
TS are equal and larger than the
Controls group, a hypothesis (H3) which assumes all group variances are equal, and a complement hypothesis (H4). To do this we run the following:
<- "Controls = TS < ADHD; Controls < TS = ADHD; Controls = TS = ADHD" hypothesis set.seed(358) <- BF(bartlett1, hypothesis)BF_var
A comprehensive output of this analysis can be obtained by running:
which presents the results of an exploratory analysis and the results of a confirmatory analysis (based on the hypotheses formulated under the
hypothesis argument). The exploratory analysis tests a hypothesis which assumes that the variances are equal across groups (homogeneity of variances) versus an alternative unrestricted hypothesis. The output shows that the posterior probabilities of these two hypotheses are approximately 0.803 and 0.197 (assuming equal priori probabilities). Note that the p value in the classical Bartlett test for these data equals 0.1638 which implies that the hypothesis of homogeneity of variances cannot be rejected using common significance levels, such as 0.05 or 0.01. Note however that this p value cannot be used as a measure for the evidence in the data in favor of homogeneity of group variances. This can be done using the proposed Bayes factor test which shows that the probability that the variances are equal is approximately 0.803. Also note that the exploratory test could equivalently tested via the
hypothesis argument by running
BF(bartlett1, "Controls = TS = ADHD").
The confirmatory test shows that H1 receives strongest support from the data, but H2 and H3 are viable competitors. It appears that even the complement H3 cannot be ruled out entirely given a posterior prob- ability of 0.058. To conclude, the results indicate that TS population are as heterogeneous in their attentional performances as the healthy control population in this specific task, but further research would be required to obtain more conclusive evidence.
For more information about the methodology, we refer the interested reader to Boing-Messing et al. (2017)
An example hypothesis test is considered under a logistic regression model. First a logistic regression model is fitted using the
<- glm(sent ~ ztrust + zfWHR + zAfro + glasses + attract + maturity + fit_glm family = binomial(), data = wilson) tattoos,
By default exploratory exhaustive tests are executed of whether the separate regression coefficients are zero, negative, or positive:
The names of the regression coefficients on which constrained hypotheses can be formualted can be extracted using the
Two different hypotheses are formulated with competing equality and/or order constraints on the regression coefficients of interest Mulder et al. (2021, Section 4.4)
<- BF(fit_glm, hypothesis = "ztrust > (zfWHR, zAfro) > 0; BF_glm ztrust > zfWHR = zAfro = 0") summary(BF_glm)
By calling the
summary function on the output object of class
BF, the results of the exploratory tests are presented of whether each separate parameter is zero, negative, or positive, and the results of the confirmatory test of the hypotheses under the
hypothesis argument are presented. When the hypotheses do not cover the complete parameter space, by default the complement hypothesis is added which covers the remaining parameter space that is not covered by the constraints under the hypotheses of interest. In the above example, the complement hypothesis covers the parameter space where neither
"ztrust > (zfWHR, zAfro) > 0" holds, nor where
"ztrust > zfWHR = zAfro = 0" holds.
For more information about the methodology, we refer the interested reader to Gu et al. (2018) and Mulder et al. (2021)
BF performs exhaustice tests of whether the separate correlations are zero, negative, or positive.
set.seed(123) <- cor_test(memory[,1:3]) cor1 <- BF(cor1) BF1 print(BF1)
The names of the correlations is constructed using the names of the variables separated by
Specific hypotheses based on prior/theoretical considerations can be tested using the
hypothesis argument. As an example we show here how to test whether all correlations are equal and positive versus its complement.
<- BF(cor1, hypothesis = "Del_with_Im = Wmn_with_Im = Wmn_with_Del > 0") BF2 print(BF2)
We can also test equality and order constraints on correlations across different groups. As the seventh column of the
memory object is a group indicator, let us first create different objects for the two different groups, and perform Bayesian estimation on the correlation matrices of the two different groups
<- subset(memory,Group=="HC")[,-(4:7)] memoryHC <- subset(memory,Group=="SZ")[,-(4:7)] memorySZ set.seed(123) <- cor_test(memoryHC,memorySZ)cor1
In this case with multiple groups by default exploratory tests are executed of whether the correlations are zero, negative, or positive for each separate group (e.g., correlations in group
gr1 are denoted by
_in_gr1 at the end of the name)
Next we test the one-sided hypothesis that the respective correlations in the first group (
g1) are larger than the correlations in the second group (
set.seed(123) <- BF(cor1, hypothesis = BF6_cor "Del_with_Im_in_g1 > Del_with_Im_in_g2 & Del_with_Wmn_in_g1 > Del_with_Wmn_in_g2 & Im_with_Wmn_in_g1 > Im_with_Wmn_in_g2")
print(BF6_cor), the output shows that the one-sided hypothesis received a posterior probability of 0.991 and the alternative received a posterior probability of .009 (assuming equal prior probabilities).
For more information about the methodology, we refer the interested reader to Mulder (2016) and Mulder and Gelissen (2019)
For a univariate regression model, by default an exhaustive test is executed of whether an effect is zero, negative, or postive.
<- lm(Superficial ~ Face + Vehicle, data = fmri) lm1 <- BF(lm1) BF1 print(BF1)
Hypotheses can be tested with equality and/or order constraints on the effects of interest. If prefered the complement hypothesis can be omitted using the
<- BF(lm1, hypothesis = "Vehicle > 0 & Face < 0; Vehicle = Face = 0", BF2 complement = FALSE) print(BF2)
In a multivariate regression model hypotheses can be tested on the effects on the same dependent variable, and on effects across different dependent variables. The name of an effect is constructed as the name of the predictor variable and the dependent variable separated by
_on_. Testing hypotheses with both constraints within a dependent variable and across dependent variables makes use of a Monte Carlo estimate which may take a few seconds.
<- lm(cbind(Superficial, Middle, Deep) ~ Face + Vehicle, lm2 data = fmri) <- "Face_on_Deep = Face_on_Superficial = Face_on_Middle < 0 < constraint2 Vehicle_on_Deep = Vehicle_on_Superficial = Vehicle_on_Middle; Face_on_Deep < Face_on_Superficial = Face_on_Middle < 0 < Vehicle_on_Deep = Vehicle_on_Superficial = Vehicle_on_Middle" set.seed(123) <- BF(lm2, hypothesis = constraint2) BF3 summary(BF3)
For more information about the methodology, we refer the interested reader to Mulder and Olsson-Collentine (2019) and Mulder and Gu (2021)
For illustrative purposes we generate a hypothetical simulated dataset
set.seed(123) <- 0.05 tau2 <- runif(50, min = 0.01, max = 0.2) vi <- rnorm(50, mean = 0, sd = sqrt(vi+tau2))yi
tau2 denotes the true between-study heterogeneity,
vi is a vector containing the squared standard errors of 50 studies, and
yi is a vector containing the estimated effects sizes in the 50 studies. To test the overall effect size and the between-study heterogeneity using
BFpack, an initial meta-analysis needs to be executed using the
metafor package. Subsequently the output is plugged into the
<- metafor::rma(yi = yi, vi = vi) res <- BF(res)BFmeta
summary output gives the posterior probabilities for a zero, negative, and positive between-study heterogeneity
I^2 and overall effect size
mu assuming equal prior probabilities:
The results indicate evidence for positive between-study heterogeneity (suggesting that a random effects meta-analysis model is appropriate) and for a zero overall effect size.
The unconstrained estimates (posterior mean and median) and the lower and upper bound of the 95% Bayesian credible intervals can be obtained by calling:
For more information about the methodology, we refer the interested reader to Van Aert and Mulder (2021)
BFon a named vector
The input for the
BF function can also be a named vector containing the estimates of the parameters of interest. In this case the error covariance matrix of the estimates is also needed via the
Sigma argument, as well as the sample size that was used for obtaining the estimates via the
n argument. Bayes factors are then computed using Gaussian approximations of the likelihood (and posterior), similar as in classical Wald test.
We illustrate this for a Poisson regression model
<- glm(formula = breaks ~ wool + tension, data = datasets::warpbreaks, poisson1 family = poisson)
The estimates, the error covariance matrix, and the sample size are extracted from the fitted model
<- poisson1$coefficients estimates <- vcov(poisson1) covmatrix <- nobs(poisson1)samplesize
Constrained hypotheses on the parameters
names(estimates) can then be tested as follows
<- BF(estimates, Sigma = covmatrix, n = samplesize, hypothesis = BF1 "woolB > tensionM > tensionH; woolB = tensionM = tensionH")
Note that the same hypothesis test would be executed when calling
<- BF(poisson1, hypothesis = "woolB > tensionM > tensionH; BF2 woolB = tensionM = tensionH")
because the same Bayes factor is used when running
BF on an object of class
Method: Bayes factor using Gaussian approximations when calling
For more information about the methodology, we refer the interested reader to Gu et al. (2018)