Verifying JZS Bayes Factor for One-Way ANOVA

sherloconan

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.

richarddmorey

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.

sherloconan

Thank you, Dr. Morey, for your prompt response.

After changing the sign of the exponent in (Ng + 1), the computation is now correct.