Bayesian one proportion test: different results in JASP and BayesFactor package

Comments

• EJ

  November 2022

  Dear Lewend,

  Can you give the specific numbers for the test? And your H0 is p = 1/3, right?

  Cheers,

  E.J.
• Lewend

  November 2022

  Dear E.J.,

  Thank you so much and apologies for my late answer! Sorry, yes, H0 is p = 1/3. Attached is a table that contains the specific numbers. For instance, the second row shows that BF10 = 0.038. When using the BayesFactor package, I get BF10 = 0.30.

  Many thanks again!

  Kind regards,

  Lewend

  bayesianproptest.PNG

  
  
• EJ

  November 2022

  Hi Lewend,

  Well, these are different models. JASP uses the test proposed by Jeffreys -- under H1, theta gets a uniform prior from 0 to 1. You can change that distribution by adjusting the a and b parameters of the beta distribution. In BayesFactor, the prior is a logistic distribution on the logit of theta, centered at 0.3. I like the BayesFactor specification, but it is more difficult to communicate. I do think we ought to implement it, and in such a way that it is easy to see what the specification of the prior on the logit scale means on the probability scale...

  Cheers,

  E.J.
• Lewend

  November 2022

  Dear E.J.,

  Thank you so much for the explanation. This really helps!

  Kind regards,

  Lewend
• EJ

  November 2022

  I'll look into it a bit more.

  E.J.
• Lewend

  November 2022

  Thank you! I'll stay tuned then!

  Kind regards,

  Lewend
• EJ

  November 2022

  OK. So first I confirm that in BayesFactor yields BF01 = 3.3 (approximately)

  #R code:

  library(BayesFactor)

  proportionBF(52,173,1/3)

  Second, the default Bayesian test in JASP gives BF01 = 7.5 (approximately). See screenshot.

  https://forum.cogsci.nl/uploads/046/RWTZ1HESB118.png

  
  

  Now let's figure out what beta prior comes closest to the prior from the BayesFactor test. To do this, we use the JAGS module. Steps:

  First I'm going to deviate slightly from the BayesFactor specification and assign logit(p) a Normal distribution, centered on logit(1/3), which is log(.5). Using JAGS I'll draw from this prior and transform the samples back to the rate scale. Code:

  model{

  logit.p0 <- log(.5)

  logit.p1 ~ dnorm(logit.p0,r)

  r <- .707

  p1 <- exp(logit.p1)/(1+exp(logit.p1))

  }

  I am monitoring p1 and the MCMC samples indicate this distribution has a mean of .37 and an SD of .22:

  https://forum.cogsci.nl/uploads/530/KBZ0EU559SBD.png

  
  

  I'll match these to the moments of the beta distribution and this yields a = 1.41 and b = 2.40. So now I'll return to the default test in JASP, but instead of specifying a beta(1,1) prior I'll use a beta(1.41,2.40) prior. The results:

  https://forum.cogsci.nl/uploads/960/745BGRLDW2LD.png

  
  

  So now we have BF01 = 4.7. closer to 3.3 then before, but I'd like it to be a little closer still. Perhaps the beta does not quite fit to the logit-normal; or perhaps it is the fact that I used a normal instead of a logistic. Or perhaps I made a mistake!

  E.J.
• Lewend

  November 2022

  Hi E.J.,

  Apologies again for the late answer. This is so helpful, thank you so much! Not sure why but I feel like the test in JASP is somewhat more intuitive.

  Kind regards,

  Lewend
• Quentin

  November 2022

  Hi Lewend and E.J.,

  
  

  A couple of comments on E.J.'s implementation:

  1) logit.p1 ~ dnorm(logit.p0,r): JAGS uses the precision parameterization, not the SD parameterization.

  2) r <- .707: The default in proportionBF seems to be 0.5.

  3) If one wants to match the logistic and normal better, one could use the actual SD of the logistic which is given by: sqrt(r^2 * pi^2 / 3).

  When I adjust these things, I obtain a = 2.15 and b = 3.87 which results in BF01 = 3.729. R code below.

  
  

  library(BayesFactor)

  r <- .5

  1 / proportionBF(52,173,1/3, rscale = r)

  
  

  sd.logistic <- sqrt(r^2 * pi^2 / 3) # https://en.wikipedia.org/wiki/Logistic_distribution

  
  

  logit.p0 <- log(.5)

  logit.p1 <- rnorm(1e6, logit.p0, sd.logistic)

  p1 <- exp(logit.p1)/(1+exp(logit.p1))

  
  

  m <- mean(p1)

  s <- sd(p1)

  
  

  print(m)

  print(s)

  
  

  a <- m * (m * (1 - m) / s^2 - 1)

  b <- (1 - m) * (m * (1 - m) / s^2 - 1)

  
  

  print(a)

  print(b)

  
  

  hist(p1, probability = TRUE, breaks = 30)

  plot(function(x) dbeta(x, a, b), add = TRUE)
• EJ

  December 2022

  Thanks Quentin!

  E.J.
• Lewend

  January 2023

  Hi Quentin,

  Apologies for the late reply! Thank you so much for your help and the detailed explanation!

  Kind regards,

  Lewend