Verifying JZS Bayes Factor for One-Way ANOVA
Hello, "BayesFactor" users.
I am revisiting Morey et al. (2011) to verify the mathematical expression for the Jeffreys-Zellner-Siow Bayes factor for balanced one-way ANOVA, as presented on page 374 (Equation 9). However, my results seem to differ from those produced by the anovaBF() function.
The R code I used is as follows:
JZSBF10 <- function(r=0.5, N, J, F_ratio) {
#' Input -
#' r: the prior scale on the variability of fixed effects, i.e.,
#' `rscaleFixed="medium"` or `rscaleFixed=0.5` in BayesFactor::anovaBF
#' N: sample size in a balanced group
#' J: number of groups
#' F_ratio: F-ratio
#' Output - the JZS Bayes factor for balanced one-way ANOVA in Eq. 9
integrand <- function(g) {
K <- 1 + (N * g) / (((J - 1) / (N * J - J)) * F_ratio + 1)
(r / sqrt(2 * pi)) * g^(-3/2) * exp(-r^2 / (2 * g)) * (N * g + 1)^(-(N * J - J) / 2) * K^(-(N * J - 1) / 2)
}
integrate(integrand, lower=0, upper=Inf)
}
test <- aov(weight ~ group, PlantGrowth) # one-way ANOVA with three balanced groups, n = 10
F_ratio <- summary(test)[[1]]["group", "F value"] # F-ratio
JZSBF10(N=10, J=3, F_ratio=F_ratio) # verifying Eq. 9
BayesFactor::anovaBF(weight ~ group, PlantGrowth, progress=F)
The integral returns `2.237442e-05 with absolute error < 4.1e-05`, while the anovaBF() function returns `3.896995 ±0.01%`.
Thanks for any comments.
Comments
I've done some digging and this looks like a typo in the paper. There is a stray negative exponent for (Ng+1) - if you change that to a positive exponent, you'll get the same answer as BayesFactor. Looking back at Liang et al (2008) the exponent on the (1+g) term (that becomes 1+Ng in our setup) is indeed positive, not negative. If you're super curious I could re-derive this for you (or, if you'd like to try it yourself I can help).
You can see the integrand that BayesFactor uses here: https://github.com/richarddmorey/BayesFactor/blob/master/pkg/BayesFactor/R/oneWayAOV-utility.R
I've used omega = K/(N*g+1) but you can see it is the same with a bit of algebra.
Thank you, Dr. Morey, for your prompt response.
After changing the sign of the exponent in (Ng + 1), the computation is now correct.