9 Bayesian Inference for 1-way ANOVAs

There are several ways to approach Bayesian Inference, even for a simple design such as a 1way ANOVA. Options exist both within and outside of the R ecosystem. This chapter does not intend to be comprehensive in that regard. One approach that has gained some currency in psychological research is the use of bayes factors as an element of evaluating hypotheses. That is the focus taken here, using the BayesFactor package.

The emphasis on bayes factors here is not to be taken as a strong recommendation that Bayesian inference for ANOVA should proceed that way. There are many opposing voices. bayes factor usage has been criticised on a variety of grounds and it is beyond the scope of this tutorial document to explore that discussion. One point of distinction is the characterization of a bayes factor approach as too focused on hypotheses, rather than estimation. In that regard alternatives exist, for example:

1. The bayesanova package (Kelter, 2022) avoids bayes factors and uses a different kind of posterior index.
2. Fuloria has provided a tutorial on use of the bayesanova package.
3. Golicher has provided a tutorial on use of the rjags package for Bayesian Anova.

9.1 A BayesFactor approach to Oneway ANOVA

This section provides simple code to extract a bayes factor for the omnibus model and then attempts to demonstrate a way to evaluate contrasts within that model.

Initially, using the BayesFactor package of Morey et. al. (Morey & Rouder, 2018), we can use their default anovaBF approach. In order to set up the analysis for contrasts that are done below, the contrast for “factora” is re-specified here as the standard orthogonal set used above. This BayesFactor approach is the alternative to the “fit.1” model used above.

The ‘anovaBF’ default approach is to compare the full model to the intercept-only model. A Jeffries prior is used for in the default approach of the anovaBF function for hypothetical population means and variances. Some control over the characteristics of the prior distribution is possible in anovaBF, but that is beyond the scope of this document. Readers are urged to begin with the help pages for the BayesFactor package and the anovaBF function.

#library(BayesFactor)
#set contrasts for factora to the orthogonal set
contrasts.factora <- matrix(c(2,-1,-1, 0,1,-1),ncol=2)
contrasts(hays$factora) <- contrasts.factora
#contrasts(factora)
BayesFactor::anovaBF(dv~factora, data=hays)

Bayes factor analysis
--------------
[1] factora : 6.926757 ±0.01%

Against denominator:
  Intercept only 
---
Bayes factor type: BFlinearModel, JZS

This indicates some evidence against the null, roughly 7 times the evidence for the null. It is not a strikingly strong degree of evidence.

9.2 Analytical Contrasts and BayesFactor

A direct BF approach to individual contrasts with the BayesFactor doesn’t appear to be available. An alternative would be to use the ‘lmBF’ function on the analogous “lm” model. However, ‘lmBF’ won’t work with IVs that are factors, so we will try to trick it. If we manually code the contrasts as new variables in the data frame we can name those as IVs. I’ve created a modified .csv file that has these new variables (ac1 and ac2). Read the new file and look at the first few and last few lines:

hays3 <- read.csv("data/hays3.csv", stringsAsFactors=T)
knitr::kable(psych::headTail(hays3))

	factora	dv	ac1	ac2
1	control	27	2	0
2	control	28	2	0
3	control	33	2	0
4	control	19	2	0
…	NA	…	…	…
27	slow	22	-1	-1
28	slow	17	-1	-1
29	slow	20	-1	-1
30	slow	23	-1	-1

The lmBF function allows us to evaluate those models and to compare them. The substantially larger Bayes Factor for the ac1-only model indicates it is a better fit. This is because (as we know), the means of the 2nd and 3rd groups are so close so we don’t really need that second vector in the model.

lmbf.full <- BayesFactor::lmBF(dv~ac1+ac2, data=hays3)
lmbf.ac1 <-  BayesFactor::lmBF(dv~ac1, data=hays3)
#lmbf.ac2 <-  BayesFactor::lmBF(dv~ac2, data=hays3)
#lmbf.full
#lmbf.ac1
ful_vs_ac1 <- c(lmbf.full,lmbf.ac1)
ful_vs_ac1

Bayes factor analysis
--------------
[1] ac1 + ac2 : 6.959137 ±0%
[2] ac1       : 20.4168  ±0%

Against denominator:
  Intercept only 
---
Bayes factor type: BFlinearModel, JZS

A different function from the same package permits direct comparison of all regression models for a given set of IVs. It gives results that mirror what we saw above.

BayesFactor::generalTestBF(dv~ac1+ac2, data=hays3, whichModels="all")

Bayes factor analysis
--------------
[1] ac1       : 20.4168   ±0%
[2] ac2       : 0.3530796 ±0%
[3] ac1 + ac2 : 6.959137  ±0%

Against denominator:
  Intercept only 
---
Bayes factor type: BFlinearModel, JZS

To use this approach, it is recommended that the student explore a good bit more background on use of bayes factor approaches. The detailed treatment of this method is beyond the bounds of the 511 class.

One can also examine considerable online help from Morey and others on the BayesFactor package:

Morey on BayesFactor

Allaire, J., Xie, Y., McPherson, J., Luraschi, J., Ushey, K., Atkins, A., … Iannone, R. (2018). Rmarkdown: Dynamic documents for r. Retrieved from https://CRAN.R-project.org/package=rmarkdown

Howell, D. C. (2013). Statistical methods for psychology (8th ed., pp. xix, 770 p.). Book, Belmont, CA: Wadsworth Cengage Learning.

Kelter, R. (2022). Bayesanova: An r package for bayesian inference in the analysis of variance via markov chain monte carlo in gaussian mixture models. The R Journal, 14, 54–78.

Morey, R. D., & Rouder, J. N. (2018). BayesFactor: Computation of bayes factors for common designs. Retrieved from https://CRAN.R-project.org/package=BayesFactor