1 – Introduction

1Ubiquitous across the empirical sciences, analysis of variance (ANOVA) allows researchers to assess the effects of categorical predictors on a continuous outcome variable. Consider for instance an experiment by Strack, Martin, & Stepper (1988) designed to test the facial feedback hypothesis, that is, the hypothesis that people’s affective responses can be influenced by their own facial expression. Participants were randomly assigned to one of three conditions. In the lips condition, participants were instructed to hold a pen with their lips, inducing a pout. In the teeth condition, participants were instructed to hold a pen between their teeth, inducing a smile. In the control condition, participants were told to hold a pen in their nondominant hand. With the pen in the instructed position, each participant then rated four cartoons for funniness. The outcome variable was the average funniness rating across the four cartoons. The ANOVA procedure may be used to test the null hypothesis that the pen position does not result in different funniness ratings.

2ANOVAs are typically conducted using frequentist statistics, where p-values decide statistical significance in an all-or-none manner: if p < .05, the result is deemed statistically significant and the null hypothesis is rejected; if p > .05, the result is deemed statistically nonsignificant, and the null hypothesis is retained. Such binary thinking has been critiqued extensively (e.g., Amrhein, Greenland, & McShane, 2019; Cohen, 1994; Rouder, Engelhardt, McCabe, & Morey, 2016), and some perceive it as a cause of the reproducibility crisis in psychology (Cumming, 2014; but see Savalei & Dunn, 2015). In recent years, several alternatives to p-values have been suggested, for example reporting confidence intervals (Cumming, 2014; Gardner & Altman, 1986) or abandoning null hypothesis testing altogether (B. B. McShane, Gal, Gelman, Robert, & Tackett, 2019).

3Here we focus on another alternative: Bayesian inference. In the Bayesian framework, knowledge about parameters and hypotheses is updated as a function of predictive success—hypotheses that predicted the observed data relatively well receive a boost in credibility, whereas hypotheses that predicted the data relatively poorly suffer a decline (Wagenmakers, Morey, & Lee, 2016). A series of recent articles show how the Bayesian framework can supplement or supplant the frequentist p-value (e.g, Burton, Gurrin, & Campbell, 1998; Dienes & McLatchie, 2018; Jarosz & Wiley, 2014; Masson, 2011; Nathoo & Masson, 2016; Rouder et al., 2016).

4The advantages of the Bayesian paradigm over the frequentist p-value are well documented (e.g., Wagenmakers et al., 2018); for instance, with Bayesian inference researchers can incorporate prior knowledge and quantify support, both in favor and against the null-hypothesis; furthermore, this support may be monitored as the data accumulate (Stefan, Gronau, Schönbrodt, & Wagenmakers, 2019). Despite these and other advantages, Bayesian analyses are still used only sparingly in the social sciences (van der Schoot, Winter, Ryan, Zondervan–Zwijnenburg, & Depaoli, 2017). The broad adoption of Bayesian statistics—and Bayesian ANOVA in particular—is hindered by the fact that Bayesian concepts are rarely taught in applied statistics courses. Consequently, practitioners may be unsure of how to conduct a Bayesian ANOVA and interpret the results.

5To help familiarize researchers with Bayesian inference for common experimental designs, this article provides a guide for conducting and interpreting a Bayesian ANOVA with JASP (JASP Team, 2019). JASP is a free, open-source statistical software program with a graphical user interface that offers both Bayesian and frequentist analyses. Below, we first provide a brief introduction to Bayesian statistics. Subsequently, we use two data examples to explain the key concepts of ANOVA.

2 – Bayesian foundations

6This section explains some of the fundamentals of Bayesian inference. We focus on interpretation rather than mathematical detail; see the special issue on Bayesian inference by Vandekerckhove, Rouder, & Kruschke (2018) for a set of comprehensive, low-level introductions to Bayesian inference. The central goal of Bayesian inference is learning, that is, using observations to update knowledge. In an ANOVA we want to learn about the candidate models M and their condition-effect parameters β. Returning to the example of the facial feedback experiment, we commonly specify two models. The null model describes the funniness ratings using a single grand average across all three conditions, effectively stating that there is no effect of pen position. The parameters of the null model are thus the average test score and the error variance. The alternative model describes the funniness ratings using an overall average and the effect of pen position; in other words, the means of the three condition are allowed to differ. Therefore, the alternative model has five parameters: the average funniness ratings across participants, the error variance, and for each of the three pen positions the magnitude of the effect. [1]

7To start the learning process, we need to specify prior beliefs about the plausibility of each model, p(M), and about the plausible parameter values β within each model, p(β|M). [2] These prior beliefs are represented by prior distributions. Observing data D drives an update of beliefs, transforming the prior distribution over models and parameters to a joint posterior distribution, denoted p(β, M|D). [3] The updating factor—the change from prior to posterior beliefs—is determined by relative predictive performance for the observed data (Wagenmakers et al., 2016). As shown in Figure 1, the knowledge updating process forms a learning cycle, such that the posterior distribution after the first batch of data becomes the prior distribution for the next batch.

Figure 1. Bayesian learning can be conceptualized as a cyclical process of up-dating knowledge in response to prediction errors. The prediction step is deductive, and the updating step is inductive. For a detailed account see Jevons (1874/1913, Chapters XI and XII). Figure available at BayesianSpectacles.org under a CC-BY license

Figure 1. L’apprentissage bayésien peut être conceptualisé comme un processus cyclique d’actualisation des connaissances en réponse à des erreurs de prédiction. L’étape de prédiction est déductive, et l’étape de mise à jour est inductive. Pour un compte rendu détaillé, voir Jevons (1874/1913, chapitres XI et XII). Figure disponible sur BayesianSpectacles.org sous une licence CC-BY

Figure 1. Bayesian learning can be conceptualized as a cyclical process of up-dating knowledge in response to prediction errors. The prediction step is deductive, and the updating step is inductive. For a detailed account see Jevons (1874/1913, Chapters XI and XII). Figure available at BayesianSpectacles.org under a CC-BY license

Figure 1. L’apprentissage bayésien peut être conceptualisé comme un processus cyclique d’actualisation des connaissances en réponse à des erreurs de prédiction. L’étape de prédiction est déductive, et l’étape de mise à jour est inductive. Pour un compte rendu détaillé, voir Jevons (1874/1913, chapitres XI et XII). Figure disponible sur BayesianSpectacles.org sous une licence CC-BY

8Mathematically, the updating process is given by Bayes’ rule:

9

10This rule stipulates how knowledge about the relative plausibility of both models and parameters ought to be updated in light of the observed data. When the focus is on the comparison of two rival models, one generally considers only the model updating term. This term, commonly known as the Bayes factor, quantifies the relative predictive performance of the rival models, that is, the change in relative model plausibility that is brought about by the data (Etz & Wagenmakers, 2017; Jeffreys, 1939; Kass & Raftery, 1995; Wrinch & Jeffreys, 1921):

11

12When the Bayes factor BF₁₀ equals 20, the observed data are twenty times more likely to occur under M₁ than under M₀ (i.e., support for M₁ versus M₀); when the Bayes factor BF₁₀ equals 1/20, the observed data are twenty times more likely to occur under M₀ than under M₁ (i.e., support for M₀ versus M₁); when the Bayes factor BF₁₀ equals 1, the observed data are equally likely to occur under both models (i.e., neither model is supported over the other). Note that the Bayes factor is a comparison of two models and hence it is always a relative measure of evidence, that is, it quantifies the performance of one model relative to another. [4]

13Since the prior and posterior odds are both ratios of probabilities, the Bayes factor can be seen as an odds ratio that quantifies the change in belief from prior odds to posterior odds. The Bayes factor can be presented as BF_10, p(D|M₁), divided by p(D|M₀), or as its reciprocal BF₀₁, p(D|M₀) over p(D|M₁). Typically, BF₁₀ > 1 is used to quantify evidence in favor of the alternative hypothesis, whereas BF₀₁ > 1 is used to quantify evidence in favor of the null hypothesis. For instance, BF₁₀ > 1/3 can be interpreted as “the data are 1/3 times more likely under M₁ than under M₀”, but for a Bayes factor lower than 1 it is more intuitive to switch numerator and denominator and instead report the results as BF₀₁ = 3, that is, “the data are 3 times more likely under M₀ than under M₁”.

14The Bayesian paradigm differs from the frequentist paradigm in at least four key aspects. First, evidence in favor of a particular model, quantified by a Bayes factor, is a continuous measure of support. Unlike the frequentist Neyman-Pearson decision rule (usually p < 0.05), there is no need to impose all-or-none Bayes factor cut-offs for accepting or rejecting a particular model. Moreover, the Bayes factor can discriminate between “absence of evidence” (i.e., nondiagnostic data that are predicted about equally well under both models, such that the Bayes factor is close to 1) and “evidence of absence” (i.e., diagnostic data that support the null hypothesis over the alternative hypothesis).

15A second difference is that, in the Bayesian paradigm, knowledge about models M and parameters β is updated simultaneously. Consequently, it is natural to account for model uncertainty by considering all models, but assigning more weight to those models that predicted the data relatively well. This procedure is known as Bayesian model averaging (BMA; Hinne, Gronau, van den Bergh, & Wagenmakers, 2019; Hoeting, Madigan, Raftery, & Volinsky, 1999; Jevons, 1874/1913; p. 296 Jeffreys, 1939; p. 365 Jeffreys, 1961). In contrast, many frequentist analyses first select a ‘best’ model and subsequently estimate its parameters, thereby neglecting model uncertainty and producing overconfident conclusions (Claeskens & Hjort, 2008, Ch 7.4). [5] Another benefit of BMA is that point estimates and uncertainty intervals can be derived without conditioning on a specific model. This way, model uncertainty is accounted for in point estimates and uncertainty intervals.

16A third difference is that the Bayesian posterior distributions allow for direct probabilistic statements about parameters. For example, based on the posterior distribution of β we can state that we are 95% confident that the parameter lies between x and y. This range of parameter values is commonly known as a 95% credible interval. [6] Similarly, we can consider any interval from a to b and quantify our confidence that the parameter falls in that specific range. A fourth difference is that Bayesian inference automatically penalizes for complexity and thus favors parsimony (e.g., Berger & Jefferys, 1992; Jeffreys, 1961; Myung & Pitt, 1997). For instance, a model with a redundant covariate will make poor predictions. Consequently, the Bayes factor, which compares the relative predictive performance of two models, will favor the model without the redundant predictor over the model with the redundant predictor. Key is that, as the words suggest, the predictive performance is assessed using parameter values that are drawn from the prior distributions.

3 – ANOVA

17Traditionally, analysis of variance involves—as the name suggests—a comparison of variances. In the frequentist framework, the variance between each level of the categorical predictor is compared to the variance within the levels of the categorical predictor. When the categorical predictor has no effect, the population variances between the levels equals the population variances within the levels, and the sample ratio of these variances is distributed according to a central F-distribution. Under the assumption that the null hypothesis is true, we may then calculate the probability of encountering a sample ratio of variances that is at least as large as the one observed—this then yields the much-maligned yet omnipresent p-value.

18Instead, the Bayesian ANOVA contrasts the predictive performance of competing models (Rouder et al., 2016). In order to make predictions the model parameters need to be assigned prior distributions. These prior distributions could in principle be specified from subjective background knowledge, but here we follow Rouder, Morey, Speckman, & Province (2012) and use a default specification inspired by linear regression models, designed to meet general desiderata such as consistency and scale invariance (i.e., it does not matter whether the outcome variable is measured in seconds or milliseconds; see also Bayarri, Berger, Forte, & García-Donato, 2012; Liang, Paulo, Molina, Clyde, & Berger, 2008).

4 – Assumptions

19Before interpreting the results from an ANOVA, it is prudent to assess whether its main assumption holds, namely that the residuals are normally distributed. A common tool to assess the normality of the residuals is a Q-Q plot, which visualizes the quantiles of the observed residuals against the quantiles expected from a standard normal distribution. If the residuals are normally distributed, then all the points in a Q-Q plot fall on the red line in Figure 2. In contrast to a frequentist ANOVA, where the residuals are point estimates, a Bayesian ANOVA provides a probability distribution for each residual. The uncertainty in the residuals can thus be summarized by 95% credible intervals. The left panel of Figure 2 shows an example where the larger quantiles lie away from the red line, displaying a substantial deviation from normality. The right panel of Figure 2 shows residuals that are more consistent with what is expected under a normal distribution.

Figure 2. Q-Q plots of non-normally distributed residuals (left) and approximately normally distributed residuals (right). The vertical bars through each point represent the 95% central credible interval. If the data are perfectly normally distributed, all points fall on the red line. Note that the y-axis of the two panels has a different scale

Figure 2. Tracés Q-Q des résidus non distribués normalement (à gauche) et des résidus distribués à peu près normalement (à droite). Les barres verticales à travers chaque point représentent l’intervalle central de crédibilité à 95%. Si les données sont parfaitement distribuées normalement, tous les points tombent sur la ligne rouge. Notez que l’axe des y des deux tracés a une échelle différente

Figure 2. Q-Q plots of non-normally distributed residuals (left) and approximately normally distributed residuals (right). The vertical bars through each point represent the 95% central credible interval. If the data are perfectly normally distributed, all points fall on the red line. Note that the y-axis of the two panels has a different scale

Figure 2. Tracés Q-Q des résidus non distribués normalement (à gauche) et des résidus distribués à peu près normalement (à droite). Les barres verticales à travers chaque point représentent l’intervalle central de crédibilité à 95%. Si les données sont parfaitement distribuées normalement, tous les points tombent sur la ligne rouge. Notez que l’axe des y des deux tracés a une échelle différente

20Introductory texts discuss additional ANOVA assumptions, most of which follow directly from the normality of the residuals. For some of these assumptions, violations can be difficult to detect visually in a Q-Q plot. An example is sphericity, which is specific to repeated measures ANOVA. One definition of sphericity is that the variance of all pairwise difference scores is equal. In the frequentist paradigm, this assumption is usually assessed using Mauchly’s test (but see Tijmstra, 2018). Another example is homogeneity of variances, which implies that the residual variance is equal across all levels of the predictors. Homogeneity of variances can be assessed using Levene’s test (Levene, 1961).

21The following sections illustrate how to conduct and interpret a Bayesian ANOVA with JASP. JASP can be freely downloaded from https://jasp-stats.org/download/. Annotated.jasp files of the discussed analyses, data sets, and a step-by-step guide on conducting a Bayesian ANOVA in JASP are available at https://osf.io/f8krs/. We should stress that the current implementation of Bayesian ANOVA in JASP is based on the R package BayesFactor (Morey & Rouder, 2015) which is itself based on the statistical work by Rouder et al. (2012).

5 – Example I: A robot’s social skills

22Do people take longer to switch off a robot when it displays social skills? This question was studied by Horstmann et al. (2018, see their publication for online access to the complete data set) and we use their data to illustrate the key concepts of a Bayesian ANOVA. In the Horstmann et al. (2018) study, 85 participants interacted with a robot. Participants were told that the purpose of their interaction with the robot was to test a new algorithm. After two dummy tasks were completed, the instructor told the participants that they could switch off the robot if they wanted. The outcome variable was the time it took participants to switch off the robot. Here we analyze the log-transformed switch-off times since the Q-Q plot of the raw switch off times showed a violation of normality. Horstmann et al. (2018) manipulated two variables in a between-subjects design. First, they manipulated the robots’ verbal responses to be either social (e.g., “Oh yes, pizza is great. One time I ate a pizza as big as me.”) or functional (e.g., “You prefer pizza. This worked well. Let us continue.”). Second, either the robot protested to being turned off (e.g., “No! Please do not switch me off! I am scared that it [sic] will not brighten up again!”) or it did not. Therefore, the design of this study is a 2x2 between-subjects ANOVA. The data are shown in Figure 3.

Figure 3. Observed log switch-off times for the data of Horstmann et al. (2018)

Figure 3. Log des temps mis pour éteindre le robot, à partir des données de Horstmann et al. (2018)

Figure 3. Observed log switch-off times for the data of Horstmann et al. (2018)

Figure 3. Log des temps mis pour éteindre le robot, à partir des données de Horstmann et al. (2018)

6 – Example II: Post hoc tests on the houses of Hogwarts

30After executing an ANOVA and finding strong evidence that a particular predictor relates to the outcome variable, a common question arises: “Which levels of the predictor deviate from one another?”. As an illustration, consider the data from Jakob, Garcia-Garzon, Jarke, & Dablander (2019) where 847 participants filled out a ‘sorting hat’ questionnaire that determined their assignment to one of the four Houses of Hogwarts from the Harry Potter books: Gryffindor, Hufflepuff, Ravenclaw, or Slytherin. [11] Subsequently, participants filled out the dark triad questionnaire (Jones & Paulhus, 2014) that was used to derive the outcome variable: Machiavellism.

31In this example, there is only one categorical predictor: The House of Hogwarts a participant was assigned to. If we compare the model with this predictor to the null model, we find overwhelming evidence for the alternative (BF₁₀ = 6.632 × 10¹⁸). This is a clear indication that Machiavellism differs between the members of the four houses. However, this result does not indicate the houses responsible for the difference. To address that question, we need a post hoc test.

Figure 4. Observed Machiavellism scores for each of the four Houses of Hogwarts

Figure 4. Scores observés de Machiavélisme pour chacune des quatre maisons de Poudlard

Figure 4. Observed Machiavellism scores for each of the four Houses of Hogwarts

Figure 4. Scores observés de Machiavélisme pour chacune des quatre maisons de Poudlard

32For ANOVA models, the main component of a post hoc test is a t-test on all pairwise combinations of a predictor’s levels. For a Bayesian ANOVA, the main component is the Bayesian t-test. Table 5 shows the Bayesian post hoc tests for the sorting hat data. As with frequentist inference, Bayesian post hoc tests are subject to a multiple comparison problem. To control for multiplicity, we follow the approach discussed in Westfall (1997) which is an extension of the approach of Jeffreys (1938); for an overview of Bayesian methods correcting for multiplicity see for instance de Jong (2019). Westfall’s approach relates the overall null hypothesis p(H₀) that all condition means are equal to each comparison between two condition means. That way, the prior probability of the overall null hypothesis can be adjusted to correct for multiplicity, and this influences each individual comparison. The procedure to relate the overall null hypothesis to each comparison is described below.

33A condition mean μ_i is either equal to the grand mean μ with probability τ, or μ_i is drawn from a continuous distribution with probability 1 – τ. It is key that this distribution is continuous because two values drawn from a continuous distribution are never exactly equal. Thus, the probability that two condition means μ_i and μ_j are equal is p(μ_i = μ_j) = p(μ_i = μ) × p(μ_j = μ) = τ². From this, the probability of the null hypothesis that all J condition means are equal follows: p(H₀) = p(µ₁ = µ₂ = … = µ_J) = p(µ₁ = µ) × p(µ₂ = µ) × … × p(µ_J = µ) = τ ^J. Solving for τ, we obtain τ = p. Thus, the prior probability that two specific magnitudes are equal can be expressed in terms of the prior probability that all magnitudes are equal, that is p(µ_i = µ_j) = τ² = p. For example, imagine there are four conditions (J = 4) and the prior probability that all condition means are equal is 0.5. Then, the prior probability that two conditions means are equal is: . The prior odds are then . In sum, the Westfall approach involves, as a first step, Bayesian t-tests for all pairwise comparisons, which provides the unadjusted Bayes factors. In the next step, the prior model odds are adjusted by fixing the overall probability of no effect to 0.5. The adjusted prior odds and the Bayes factor are then used to calculate the adjusted posterior odds.

34Table 5 shows the results for the post hoc tests of the sorting hat example. The adjusted posterior odds show (1) evidence (i.e., odds of about 16) that Machiavellism differs between Hufflepuff and Ravenclaw; (2) evidence (i.e., odds of about 27) that Machiavellism differs between Gryffindor and Hufflepuff; (3) overwhelming evidence (i.e., odds of about 1.04 × 10⁹, 5.43 × 10¹⁶, and 5.30 × 10⁹) that Machiavellism differs between Gryffindor and Slytherin, between Hufflepuff and Slytherin, and between Ravenclaw and Slytherin, respectively; (4) evidence (i.e, odds of 1/0.0432 ≈ 23) that Machiavellism of Gryffindor and Ravenclaw is the same.

35Now that we know which Houses differ, the next step is to assess the magnitude of each House of Hogwarts on Machiavellism score. Rather than examining a table that summarizes the marginal posteriors, we plot the model averaged posteriors for each house in Figure 5. Clearly, Slytherin scores higher on Machiavellism than the other Houses whereas Hufflepuff scores lower on Machiavellism than the other Houses. Table 6 in the appendix shows the parameters estimates of the marginal posterior effects for each house.

Figure 5. Posterior distributions of the effect of each House of Hogwarts on Machiavellism. Slytherin scores higher on Machiavellism than the other Houses whereas Hufflepuff scores lower on Machiavellism than the other Houses. The horizontal error bars above each density represent 95% credible intervals

Figure 5. Distributions postérieures de l’effet de chaque maison de Poudlard sur le Machiavélisme. Serpentard obtient le plus haut score de Machiavélisme par rapport aux autres maisons, alors que Poufsouffle obtient le score de Machiavélisme le plus bas. Les barres d’erreur horizontales au-dessus de chaque densité représentent des intervalles de crédibilité à 95%

Figure 5. Posterior distributions of the effect of each House of Hogwarts on Machiavellism. Slytherin scores higher on Machiavellism than the other Houses whereas Hufflepuff scores lower on Machiavellism than the other Houses. The horizontal error bars above each density represent 95% credible intervals

Figure 5. Distributions postérieures de l’effet de chaque maison de Poudlard sur le Machiavélisme. Serpentard obtient le plus haut score de Machiavélisme par rapport aux autres maisons, alors que Poufsouffle obtient le score de Machiavélisme le plus bas. Les barres d’erreur horizontales au-dessus de chaque densité représentent des intervalles de crédibilité à 95%

7 – Concluding comments

36The goal of this paper was to provide guidance for practitioners on to conduct a Bayesian ANOVA in JASP and interpret the results. Although the focus was on ANOVAs with categorical predictors, JASP can also handle ANOVAs with additional continuous predictors. The appropriate analysis then becomes an analysis of covariance (ANCOVA) and all concepts explained here still apply. For a general guide on reporting Bayesian analyses see van Doorn et al. (2019).

37As with all statistical methods, the Bayesian ANOVA comes with limitations and caveats. For instance, when the model is severely misspecified and the residuals are non-normally distributed, the results from a standard ANOVA—whether Bayesian or frequentist—are potentially misleading and should be interpreted with care. In such cases, at least two alternatives may be considered. The first alternative is to consider a rank-based ANOVA such as the Kruskal–Wallis test (Kruskal & Wallis, 1952). This test depends only on the ordinal information in the data and hence does not make strong assumptions on how the data ought to be distributed. The second alternative is to specify a different distribution for the residuals. Using software for general Bayesian inference such as Stan (Carpenter et al., 2017) or JAGS (Plummer, 2003), it is relatively straightforward to specify any distribution for the residuals. However, this approach requires knowledge about programming and statistical modeling and is likely to be computationally intensive. Another limitation of the Bayesian ANOVA is that, especially in more complicated designs, it is not straightforward to intuit what knowledge the prior distributions represent.

38Some limitations are specific to JASP. Currently, it is not possible to use post hoc tests to examine whether the contribution of a level differs from zero, that is, to test whether a specific level deviates from the grand mean. It is also not possible to handle missing values in any other way than list-wise deletion. Another limitation relates to sample size planning. The typical planning process involves a frequentist power analysis which provides the sample size needed to achieve a certain rate of correctly detecting a true effect of a prespecified magnitude. A Bayesian equivalent of power analysis is Bayes factor design analysis (BFDA; e.g., Schönbrodt & Wagenmakers, 2018). In a sequential design, BFDA produces the expected sample sizes required to reach a target level of evidence (i.e., a target Bayes factor). In a fixed-n design, BFDA produces the expected levels of evidence, given a specification of the magnitude of the effect. At the moment of writing BFDA has not been implemented in JASP; an accessible tutorial and a Shiny app are provided by Stefan et al. (2019).

39We believe that the Bayesian ANOVA provides a perspective on the analysis of factorial designs that can fruitfully supplement or even supplant the currently dominant frequentist ANOVA. The epistemic advantages of the Bayesian paradigm are well known (e.g., Jeffreys, 1961; Wagenmakers et al., 2018) but in order to be adopted in research practice it is essential for the methodology to be implemented in an easy-to-use software package such as JASP. In addition to the software, however, practitioners also require guidance on how to interpret the results, which was the main purpose of this paper. In general, we hope that the increased use of the Bayesian ANOVA will stimulate the methodological diversity in the field, and that it will become more standard to examine the robustness of frequentist conclusions by comparing them to the Bayesian alternative.