Introduction to BFpack

Bayesian t testing

ttest1 <- bain::t_test(therapeutic, alternative = "greater", mu = 5)

get_estimates(ttest1)

library(BFpack)
BF1 <- BF(ttest1)

hypothesis <- "mu = 5; mu < 5; mu > 5"
BF(ttest1, hypothesis = hypothesis)

hypothesis <- "mu = 5"
BF(ttest1, hypothesis = hypothesis)

BF(ttest1, hypothesis = "mu = 5; mu < 5; mu > 5", prior.hyp = c(.5,0,.5))

fmri1 <- cbind(fmri,ones=1)
mlm1 <- lm(cbind(Superficial,Middle,Deep) ~ -1 + ones, data = fmri1)

BFmlm1 <- BF(mlm1, hypothesis="ones_on_Superficial=ones_on_Middle=ones_on_Deep=0;
                               (ones_on_Superficial,ones_on_Middle,ones_on_Deep)>0")

Analysis of variance

First an analysis of variance (ANOVA) model is fitted using the aov fuction in R.

aov1 <- aov(price ~ anchor * motivation, data = tvprices)

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))

BF(aov1)

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:

constraints2 <- "anchorrounded > 0 & motivationlow < 0 &
  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).

Testing independent group variances

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 BF function.

bartlett1 <- bartlett_test(x = attention$accuracy, g = attention$group)

On an object of this class, by default BF executes an exploratory test of homogeneity (equality) of variances against an unconstrained (full) model

BF(bartlett1)

The group variances have names ADHD, Controls, and TS. This can be retrieved by running

get_estimates(bartlett1)

Let’s say we want to test whether a hypothesis (H1) which assumes that group variances of groups Controls and TS are equal and smaller than the group variance of the ADHD group, a hypothesis (H2) which assumes that the group variances of ADHD and 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:

hypothesis <- "Controls = TS < ADHD; Controls < TS = ADHD; Controls = TS = ADHD"
set.seed(358)
BF_var <- BF(bartlett1, hypothesis)

A comprehensive output of this analysis can be obtained by running:

summary(BF_var)

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)

Logistic regression

An example hypothesis test is considered under a logistic regression model. First a logistic regression model is fitted using the glm function

fit_glm <- glm(sent ~ ztrust + zfWHR + zAfro + glasses + attract + maturity +
               tattoos, family = binomial(), data = wilson)

By default exploratory exhaustive tests are executed of whether the separate regression coefficients are zero, negative, or positive:

BF(fit_glm)

The names of the regression coefficients on which constrained hypotheses can be formualted can be extracted using the get_estimates function.

get_estimates(fit_glm)

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_glm <- BF(fit_glm, hypothesis = "ztrust > (zfWHR, zAfro) > 0;
             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)

Correlation analysis

By default BF performs exhaustice tests of whether the separate correlations are zero, negative, or positive.

set.seed(123)
cor1 <- cor_test(memory[,1:3])
BF1 <- BF(cor1)
print(BF1)

The names of the correlations is constructed using the names of the variables separated by _with_:

get_estimates(cor1)

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.

BF2 <- BF(cor1, hypothesis = "Del_with_Im = Wmn_with_Im = Wmn_with_Del > 0")
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

memoryHC <- subset(memory,Group=="HC")[,-(4:7)]
memorySZ <- subset(memory,Group=="SZ")[,-(4:7)]
set.seed(123)
cor1 <- cor_test(memoryHC,memorySZ)

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)

get_estimates(cor1)

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 (g2) via

set.seed(123)
BF6_cor <- BF(cor1, hypothesis =
  "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")

By running 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)

Univariate/Multivariate multiple regression

For a univariate regression model, by default an exhaustive test is executed of whether an effect is zero, negative, or postive.

lm1 <- lm(Superficial ~ Face + Vehicle, data = fmri)
BF1 <- BF(lm1)
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 complement argument

BF2 <- BF(lm1, hypothesis = "Vehicle > 0 & Face < 0; Vehicle = Face = 0",
          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.

lm2 <- lm(cbind(Superficial, Middle, Deep) ~ Face + Vehicle,
              data = fmri)
constraint2 <- "Face_on_Deep = Face_on_Superficial = Face_on_Middle < 0 <
     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)
BF3 <- BF(lm2, hypothesis = constraint2)
summary(BF3)

For more information about the methodology, we refer the interested reader to Mulder and Olsson-Collentine (2019) and Mulder and Gu (2021)

Testing the effect size and beteen-study heterogeneity in a meta-analysis

For illustrative purposes we generate a hypothetical simulated dataset

set.seed(123)
tau2 <- 0.05
vi <- runif(50, min = 0.01, max = 0.2)
yi <- rnorm(50, mean = 0, sd = sqrt(vi+tau2))

where 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 BF function:

res <- metafor::rma(yi = yi, vi = vi)
BFmeta <- BF(res)

The 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:

summary(BFmeta)

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:

BFmeta$estimates

For more information about the methodology, we refer the interested reader to Van Aert and Mulder (2021)

Running BF on 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

poisson1 <- glm(formula = breaks ~ wool + tension, data = datasets::warpbreaks,
             family = poisson)

The estimates, the error covariance matrix, and the sample size are extracted from the fitted model

estimates <- poisson1$coefficients
covmatrix <- vcov(poisson1)
samplesize <- nobs(poisson1)

Constrained hypotheses on the parameters names(estimates) can then be tested as follows

BF1 <- BF(estimates, Sigma = covmatrix, n = samplesize, hypothesis = 
  "woolB > tensionM > tensionH; woolB = tensionM = tensionH")

Note that the same hypothesis test would be executed when calling

BF2 <- BF(poisson1, hypothesis = "woolB > tensionM > tensionH;
          woolB = tensionM = tensionH")

because the same Bayes factor is used when running BF on an object of class glm (see Method: Bayes factor using Gaussian approximations when calling print(BF1) and print(BF2)).

For more information about the methodology, we refer the interested reader to Gu et al. (2018)