Testing Models’ Parameters with Bayes Factors
A Bayes factor for a single parameter can be used to answer the question:
“Given the observed data, has the null hypothesis of an absence of an effect become more or less credible?”
Bayesian analysis of the Students’ (1908) Sleep data set.
Let’s use the Students’ (1908) Sleep data set (data("sleep")). The data comes from a study in which participants were administered a drug and the researchers assessed the extra hours of sleep that participants slept afterwards. We will try answering the following research question using Bayes factors:
Given the observed data, has the hypothesis that the drug (the effect of
group) has no effect on the numbers of hours of extra sleep (variableextra) become more of less credible?
The boxplot suggests that the second group has a higher number of hours of extra sleep. By how much?
Let’s fit a simple Bayesian linear model, with a prior of \(b_{group} \sim N(0, 3)\) (i.e. the prior follows a Gaussian/normal distribution with \(mean = 0\) and \(SD = 3\)), using rstanarm package:
set.seed(123)
library(rstanarm)
model <- stan_glm(
formula = extra ~ group,
data = sleep,
prior = normal(0, 3, autoscale = FALSE)
)Testing against a null-region
One way of operationlizing the null-hypothesis is by setting a null region, such that an effect that falls within this interval would be practically equivalent to the null (Kruschke, 2010). In our case, that means defining a range of effects we would consider equal to the drug having no effect at all. We can then compute the prior probability of the drug’s effect falling within this null-region, and the prior probability of the drug’s effect falling outside the null-region to get our prior odds. Say any effect smaller than an hour of extra sleep is practically equivalent to being no effect at all, we would define our prior odds as:
\[ \frac {P(b_{drug} \in [-1, 1])} {P(b_{drug} \notin [-1, 1])} \]
Given our prior has a normal distribution centered at 0 hours with a scale (an SD) of 3 hours, our priors would look like this:
and the prior odds would be 2.2.
By looking at the posterior distribution, can now compute the posterior probability of the drug’s effect falling within the null-region, and the posterior probability of the drug’s effect falling outside the null-region to get our posterior odds:
\[ \frac {P(b_{drug} \in [-1,1] | Data)} {P(b_{drug} \notin [-1,1] | Data)} \]
We can see that the center of the posterior distribution has shifted away from 0 (to ~1.5). Likewise, the posterior odds are 2, which seems to favor the effect being non-null. But, does this mean the data support the alternative over the null? Hard to say, since even before the data were observed, the priors already favored the alternative - so we need to take our priors into account here!
Let’s compute the Bayes factor as the change from the prior odds to the posterior odds: \(BF_{10} = Odds_{posterior} / Odds_{prior} = 0.9\)! This BF indicates that the data provide 1/0.9 = 1.1 times more evidence for the effect of the drug being practically nothing than it does for the drug having some clinically significant effect. Thus, although the center of distribution has shifted away from 0, and the posterior distribution seems to favor a non-null effect of the drug, it seems that given the observed data, the probability mass has overall shifted closer to the null interval, making the values in the null interval more probable! (see Non-overlapping Hypotheses in Morey & Rouder, 2011)
All of this can be achieved with the function bayesfactor_parameters(), which computes a Bayes factor for each of the model’s parameters:
My_first_BF <- bayesfactor_parameters(model, null = c(-1, 1))
My_first_BF> Bayes Factor (Null-Interval)
>
> Parameter | BF
> -------------------
> (Intercept) | 0.103
> group2 | 0.883
>
> * Evidence Against The Null: [-1.000, 1.000]We can also plot using the see package:
library(see)
plot(My_first_BF)
Note that interpretation guides for Bayes factors can be found in the effectsize package:
effectsize::interpret_bf(exp(My_first_BF$log_BF[2]), include_value = TRUE)> [1] "anecdotal evidence (BF = 1/1.13) against"
> (Rules: jeffreys1961)Testing against the point-null (0)
What if we don’t know what region would be practically equivalent to 0?
Or if we just want the null to be exactly zero? Not a problem - as the width of null region shrinks to a point, the change from the prior probability to the posterior probability of the null can be estimated by comparing the density of the null value between the two distributions.1 This ratio is called the Savage-Dickey ratio, and has the added benefit of also being an approximation of a Bayes factor comparing the estimated model against the a model in which the parameter of interest has been restricted to a point-null:
“[…] the Bayes factor for \(H_0\) versus \(H_1\) could be obtained by analytically integrating out the model parameter \(\theta\). However, the Bayes factor may likewise be obtained by only considering \(H_1\), and dividing the height of the posterior for \(\theta\) by the height of the prior for \(\theta\), at the point of interest.” (Wagenmakers, Lodewyckx, Kuriyal, & Grasman, 2010)
My_second_BF <- bayesfactor_parameters(model, null = 0)
My_second_BF> Bayes Factor (Savage-Dickey density ratio)
>
> Parameter | BF
> ----------------
> group2 | 1.24
>
> * Evidence Against The Null: 0plot(My_second_BF)
Directional hypotheses
We can also compute Bayes factors for directional hypotheses (“one sided”), if we have a prior hypotheses about the direction of the effect. This can be done by setting an order restriction on the prior distribution (which results in an order restriction on the posterior distribution) of the alternative (Morey & Wagenmakers, 2014). For example, if we have a prior hypothesis that the drug has a positive effect on the number of sleep hours, the alternative will be restricted to the region to the right of the null (point or interval):
test_group2_right <- bayesfactor_parameters(model, direction = ">")
test_group2_right> Bayes Factor (Savage-Dickey density ratio)
>
> Parameter | BF
> ----------------
> group2 | 2.37
>
> * Evidence Against The Null: 0
> * Direction: Right-Sided testplot(test_group2_right)
As we can see, given that we have an a priori assumption about the direction of the effect (that the effect is positive), the presence of an effect is 2.8 times more likely than the absence of an effect, given the observed data (or that the data are 2.8 time more probable under \(H_1\) than \(H_0\)). This indicates that, given the observed data, and a priori hypothesis, the posterior mass has shifted away from the null value, giving some evidence against the null (note that a Bayes factor of 2.8 is still considered quite weak evidence).
Thanks to the flexibility of Bayesian framework, it is also possible to compute a Bayes factor for dividing hypotheses - that is, for a null and alternative that are complementary, opposing one-sided hypotheses (Morey & Wagenmakers, 2014).
For example, above we compared an alternative of \(H_A\): the drug has a positive effects to the null \(H_0\): the drug has no effect. But we can also compare instead the same alternative to its complementary hypothesis: \(H_{-A}\): the drug has a negative effects.
test_group2_dividing <- bayesfactor_parameters(model, null = c(-Inf, 0))
test_group2_dividing> Bayes Factor (Null-Interval)
>
> Parameter | BF
> -----------------
> group2 | 20.53
>
> * Evidence Against The Null: [-Inf, 0.000]plot(test_group2_dividing)
We can see that this test produces even stronger (more conclusive) evidence than the one-sided vs. point-null test! And indeed, as a rule of thumb, the more specific the two hypotheses are, and the more distinct they are from one another, the more power our Bayes factor has!2
Thanks to the transitivity of Bayes factors, we can also use bayesfactor_parameters() to compare even more types of hypotheses, with some trickery. For example:
\[ \underbrace{BF_{0<b<1\text{ vs. }b=0}}_{\text{range vs. point}} = \underbrace{BF_{b<0\text{ vs. }b=0}}_{\text{directional vs. point}} / \underbrace{BF_{b<0\text{ vs. }0<b<1}}_{\text{directional vs. range}} \]
NOTE: See the Testing Contrasts appendix below.
Support intervals
So far we’ve seen that Bayes factors quantify relative support between competing hypotheses. However, we can also ask:
Upon observing the data, the credibility of which of the parameter’s values has increased (or decreased)?
For example, we’ve seen that the point null has become somewhat less credible after observing the data, but we might also ask which values have gained credibility given the observed data?. The resulting range of values is called the support interval as it indicates which values are supported by the data (Wagenmakers, Gronau, Dablander, & Etz, 2018). We can do this by once again comparing the prior and posterior distributions and checking where the posterior densities are higher than the prior densities.
In bayestestR, this can be achieved with the si() function:
my_first_si <- si(
posterior = data.frame(group2 = posterior),
prior = data.frame(group2 = prior),
BF = 1
)
print(my_first_si)> Support Interval
>
> Parameter | BF = 1 SI
> ------------------------
> group2 | [0.15, 3.04]The argument BF = 1 indicates that we want the interval to contain values that have gained support by a factor of at least 1 (that is, any support at all).
Visually, we can see that the credibility of all the values within this interval has increased (and likewise the credibility of all the values outside this interval has decreased):
plot(my_first_si)
We can also see the this support interval (just barely) excludes the point null (0) - whose credibility we’ve already seen has decreased by the observed data. This emphasizes the relationship between the support interval and the Bayes factor:
“The interpretation of such intervals would be analogous to how a frequentist confidence interval contains all the parameter values that would not have been rejected if tested at level \(\alpha\). For instance, a BF = 1/3 support interval encloses all values of theta for which the updating factor is not stronger than 3 against.” (Wagenmakers et al., 2018)
Thus, the choice of BF (the level of support the interval should indicate) depends on what we want our interval to represent:
- A \(BF = 1\) contains values whose credibility has merely not decreased by observing the data.
- A \(BF > 1\) contains values who received more impressive support from the data.
- A \(BF < 1\) contains values whose credibility has not been impressively decreased by observing the data. Testing against values outside this interval will produce a Bayes factor larger than \(1/BF\) in support of the alternative.
