This is a “pro” worksheet, meaning that it assumes familiarity with
the material covered in sections up to and including Going Further
With R. For this worksheet, you’ll need to use the
rminrdata project; details of how to access this are
provided on the preprocessing worksheet.
Once you’ve opened that project, create a file anova5.R
in which to enter your R commands.
Twofactor withinsubjects
In More on Bayesian ANOVA we looked at the
simplest possible way of doing this kind of analysis within R. There is
a reasonable case to be made that the approach we took there is too
simple. In particular, it assumes that the true size of the effect is
the same in every participant. In other words, it assumes the absence of
individual differences. As discussed by van den Berg et al. (2022), this
is not only unlikely to be true, where there is more than one
withinsubjects factor, it can also lead to substantively different
conclusions to a traditional (NHST) ANOVA. This led them to propose an
alternative analysis that is more similar in concept to the traditional
analysis. In what follows we demonstrate how to do this analysis for the
same data set as used in Exercise 2 of More on
Bayesian ANOVA.
Whether emulating the models of traditional ANOVA is in itself the
best thing to do is debatable, as traditional ANOVA is potentially
problematic in a number of respects other than the use of NHST. However,
this is a discussion that we’ll leave for another time  if you’re
curious, read Rouder
et al. (2022).
Load and preprocess data
The following commands are part of the answer to Exercise 2 of More on Bayesian ANOVA and are not discussed
further here:
# Load packages
library(tidyverse)
library(BayesFactor)
# Load data, select Test phase, select appropriate columns, drop NAs
raw < read_csv('casestudies/chrismitchell/primingfaces.csv') %>%
filter(Running == "Test") %>%
select(Subject, Congruency, Load, RT) %>%
drop_na()
# Create subjectlevel summary
priming < raw %>%
mutate(Subject = factor(Subject), Congruency = factor(Congruency),
Load = factor(Load)) %>%
group_by(Subject, Congruency, Load) %>% summarise(RT = mean(RT))
Calculate some initial Bayes factors
Enter this code into your script and run it; it may take a minute or
so to produce a result:
# Calculate Bayes Factors
BFo < generalTestBF(
formula = RT ~ Congruency*Load + Subject + Subject:Congruency + Subject:Load,
data = data.frame(priming),
whichRandom = "Subject",
neverExclude = c("Subject", "Subject:Congruency", "Subj:Load"),
whichModels = "all"
)
# Display Bayes Factors
BFo
Bayes factor analysis

[1] Congruency + Subject + Congruency:Subject + Load:Subject : 1.883333e+30 ±1.15%
[2] Load + Subject + Congruency:Subject + Load:Subject : 7.871082e+33 ±0.79%
[3] Congruency:Load + Subject + Congruency:Subject + Load:Subject : 2.085742e+30 ±6.18%
[4] Congruency + Load + Subject + Congruency:Subject + Load:Subject : 2.734706e+33 ±1.83%
[5] Congruency + Congruency:Load + Subject + Congruency:Subject + Load:Subject : 6.585271e+29 ±1.56%
[6] Load + Congruency:Load + Subject + Congruency:Subject + Load:Subject : 2.940871e+33 ±2.25%
[7] Congruency + Load + Congruency:Load + Subject + Congruency:Subject + Load:Subject : 1.240869e+33 ±20.75%
[8] Subject + Congruency:Subject + Load:Subject : 5.655589e+30 ±0.38%
Against denominator:
Intercept only

Bayes factor type: BFlinearModel, JZS
Explanation of command
Much of the above code is the same as previous worksheets and is not
redescribed here. The new components are:
generalTestBF
 Where you previously used the
anovaBF
command, we now use this more powerful command,
which we need in order to do this more complex analysis.
formula = RT ~ Congruency*Load + Subject + Subject:Congruency + Subject:Load
 This specifies what you wish to calculate. It has the same form as the
simpler calculation you did in Exercise 2, but adds two new terms:
Subject:Congruency
and Subject:Load
. The
:
means “interaction”. So, the two new terms we’re adding
here are an interaction of the Subject factor with the Congruency
condition, and an interaction of the Subject factor with the Load
condition. These interactions allow for the possibility that (i) the
true size of the effect of congruency differs between individuals, and
(ii) the true size of the effect of Load differs between
individuals.
neverExclude = c("Subject", "Subject:Congruency", "Subj:Load"), whichModels = "all"
 There is nothing conceptually deep about these two lines of code, they
just tell the generalTestBF
command to generate the set of
Bayes Factors we need in the next step. When using this code with your
own data, all you need to do is ensure that the
neverExclude
terms include the Subject factor, and the
interaction of that factor with each of your two withinsubject factors
(in this case Congruency and Load)
Explanation of output
Although longer, the table is much the same as the one we saw in
Exercise 2. Each row is a Bayes Factor, and in each case that Bayes
Factor is calculated against a null model (Intercept only

the model that RT is unaffected by any of the factors, including the
participants). So, for example,
[8] Subject + Congruency:Subject + Load:Subject
has a large
BF in favour of the alternative model that RT is affected by individual
differences in overall RT and the size of the Congruency and Load
effects.
Calculate further Bayes Factors
In this (and many other) experiments, our focus is on the overall
effect of our experimental manipulation (Congruency, Load) on RT, rather
than on individual differences. We can get Bayes Factors for these
effects by using the subject model (model 8 above), rather than
the intercept, as the null model.
We do this in much the same way as we calculated the Bayes Factor for
the interaction in Exercise 2, i.e. by the /
operator. What
we want to do is apply that operator for each of the models, which we
can do with the following command:
# Recalculate Bayes Factors, taking subject model as denominator model
BF < BFo[length(BFo)] / BFo[length(BFo)]
# Display results
BF
Bayes factor analysis

[1] Congruency + Subject + Congruency:Subject + Load:Subject : 0.3330038 ±1.21%
[2] Load + Subject + Congruency:Subject + Load:Subject : 1391.735 ±0.87%
[3] Congruency:Load + Subject + Congruency:Subject + Load:Subject : 0.368793 ±6.19%
[4] Congruency + Load + Subject + Congruency:Subject + Load:Subject : 483.5404 ±1.87%
[5] Congruency + Congruency:Load + Subject + Congruency:Subject + Load:Subject : 0.1164383 ±1.61%
[6] Load + Congruency:Load + Subject + Congruency:Subject + Load:Subject : 519.9937 ±2.29%
[7] Congruency + Load + Congruency:Load + Subject + Congruency:Subject + Load:Subject : 219.4058 ±20.76%
Against denominator:
RT ~ Subject + Congruency:Subject + Load:Subject

Bayes factor type: BFlinearModel, JZS
Explanation of output
The effect of this command is to give us a Bayes Factor for each of
the other models, relative to the subject model. It achieves
this in the following way:
Explanation of command
length(BFo)
 This gives the number of rows in the table
of results (8 in this case)
BFo[length(BFo)]
 This returns the last row of the
table
BFo[length(BFo)]
 The 
operator in this
case means “remove”, so this returns all the rows apart from the last
one.
BFo[length(BFo)] / BFo[length(BFo)]
 Putting the two
previous commands together, this gives us a BF for each of the models
(except the last one), compared to the last one.
Pruning the Bayes Factors
That’s quite a lot of Bayes Factors! Some of these, we’re not going
to use. Specifically, there are Bayes Factors in that list for models
where interactions between the manipulated factors (Congruency, Load)
are assumed without the presence of corresponding main effects of those
manipulated factors. For example, model 3 is
Congruency:Load + Subject + Congruency:Subject + Load:Subject
,
which includes an interaction between Congurency and Load, but not a
main effect of either Congruency or Load. These kinds of models are
illdefined, which is sometimes considered to be a bad thing.
So, we won’t consider them here. A relatively accessible discussion of
the problems of illdefined models is provided by Rouder
et al. (2022)
We remove the models that are illdefined on the manipulated factors
by identifying them manually, and then removing them from the list. In
this case, the illdefined models are 3, 5 and 6; and we remove them
using the 
operator.
# Remove illdefined models
BFp < BF[c(3, 5, 6)]
# Display remaining set
BFp
Bayes factor analysis

[1] Congruency + Subject + Congruency:Subject + Load:Subject : 0.3330038 ±1.21%
[2] Load + Subject + Congruency:Subject + Load:Subject : 1391.735 ±0.87%
[3] Congruency + Load + Subject + Congruency:Subject + Load:Subject : 483.5404 ±1.87%
[4] Congruency + Load + Congruency:Load + Subject + Congruency:Subject + Load:Subject : 219.4058 ±20.76%
Against denominator:
RT ~ Subject + Congruency:Subject + Load:Subject

Bayes factor type: BFlinearModel, JZS
Main effects
At this point, we have an table of Bayes Factors of the same form as
Exercise 2 of the previous worksheet. Row 1 is the main effect of
Conrguency. Row 2 is the main effect of Load. The Bayes Factors are not
quite the same as in Exercise 2  we shouldn’t expect them to be, as
this is a different test. But the conclusions are the same. There is
moderate evidence for the null in the case of Congruency, and strong
evidence for an effect of Load.
Interaction
Finally, as in Exercise 2, we assess evidence for the interaction by
comparing two models  the model that contains just the two main effects
(model 3), and the model that also contains the interaction term (model
4)
# Calculate Bayes Factor for interaction
BFp[4] / BFp[3]
Bayes factor analysis

[1] Congruency + Load + Congruency:Load + Subject + Congruency:Subject + Load:Subject : 0.4537487 ±20.84%
Against denominator:
RT ~ Congruency + Load + Subject + Congruency:Subject + Load:Subject

Bayes factor type: BFlinearModel, JZS
Again, the Bayes Factor is not the same as before, because this is a
different test, but the conclusion is similar  there is moderate
evidence for the absence of an interaction.
Robustness of analysis
You might feel a bit uneasy that the results of two analyses  both
called Bayesian ANOVAs  are not exactly the same. You should not.
Analyzing data always involves making assumptions about those data and
the underlying psychological processes that generated them. For example,
here, the assumptions you make about the presence of individual
differences can matter for your conclusions about the presence of group
effects. Making reasonable assumptions is part of being a scientist, and
your conclusions can be affected by the assumptions that you make.
What we saw in this case was that Exercise 2, and the current
analysis, led to largely the same substantive conclusions, despite
starting from different assumptions. When this happens, it increases our
confidence in the conclusions we have drawn  we say the analysis is
robust across a range of assumptions. If your substantive conclusions
depend on your assumptions, and those assumptions are not well founded,
then your data may not be good enough to answer your questions.
This material is distributed under a Creative Commons licence.
CCBYSA 4.0.