Bayesian RM ANOVA w/ random slopes

Comments

• vandenman

  June 2022

  Hi Corey,

  
  

  Happy to hear someone read the preprint 😊

  
  

  That is correct, you need to use lmBF which only fits one model at a time (anovaBF does the wrong thing, or well, the wrong thing in this scenario). Let's look at an example in the jasp data library, namely Alcohol attitudes.

  In jasp 0.16.3 I get this output:

  

  To replicate this in R you can run this code:

  library(tidyr)
  library(BayesFactor)
  # read the data file directly from GitHub (https://github.com/jasp-stats/jasp-desktop/blob/stable/Resources/Data%20Sets/Data%20Library/3.%20ANOVA/Alcohol%20Attitudes.csv)
  dat <- read.csv("https://raw.githubusercontent.com/jasp-stats/jasp-desktop/stable/Resources/Data%20Sets/Data%20Library/3.%20ANOVA/Alcohol%20Attitudes.csv")
  
  # reshape the data to match bayes factor input
  dat_long <- dat |> pivot_longer(!participant, names_to = c("Drink", "Imagery"), names_pattern = "(beer|wine|water|)(pos|neg|neu)", values_to = "DV")
  
  # these formulas are indeed the full and null model
  full <- DV ~ Drink + Imagery + participant + Drink*Imagery + Drink*participant + Imagery*participant
  null <- DV ~                   participant +                 Drink*participant + Imagery*participant
  
  # however, there is an easier way to generate the entire modelspace
  # generate all formulas without random slopes
  formulas <- BayesFactor::enumerateGeneralModels(DV ~ Drink + Imagery + participant + Drink*Imagery, whichModels = "withmain", neverExclude = "participant")
  
  # add the random slopes
  formulas <- lapply(formulas, function(f) {
    update.formula(f, ~ . + Drink:participant + Imagery:participant)
  })
  
  # fit all models
  bfs <- lapply(formulas, lmBF, data = dat_long)
  
  # reduce the list to vector type inside BayesFactor
  bfs_c <- Reduce(c, bfs)
  
  # divide by the best performing model an sort decreasingly, as in the Jasp output
  sort(bfs_c / max(bfs_c), decreasing = TRUE)
  #> Bayes factor analysis
  #> --------------
  #> [1] Drink + Imagery + participant + Drink:Imagery + Drink:participant + Imagery:participant : 1            ±0%
  #> [2] Drink + Imagery + participant + Drink:participant + Imagery:participant                 : 6.211023e-08 ±3.21%
  #> [3] Imagery + participant + participant:Drink + Imagery:participant                         : 3.469623e-12 ±2.89%
  #> [4] Drink + participant + Drink:participant + participant:Imagery                           : 2.777261e-41 ±3.49%
  #> [5] participant + participant:Drink + participant:Imagery                                   : 6.836036e-42 ±2.72%
  #> 
  #> Against denominator:
  #>   DV ~ Drink + Imagery + participant + Drink:Imagery + Drink:participant + Imagery:participant 
  #> ---
  #> Bayes factor type: BFlinearModel, JZS
  

  The order of the models is exactly the same as in Jasp. Note that there are small differences in the values of the Bayes factor due to approximation error. If you increase the numerical accuracy in both Jasp and R then the differences should disappear, although that will (significantly) increase the running time.

  
  

  Let me know if anything is unclear!

  Cheers,

  Don
• coreyw

  June 2022 edited June 2022

  Hi Don,

  That is very helpful! Thanks for sharing the code. My next question would be how to progress from the model comparison output to the analysis of effects output...

  • If using a matched models approach would we compute the main and interaction effects as follows?

  mainEff <- c(bfs_c[1]/bfs_c[5], bfs_c[2]/bfs_c[5]) # drink, imagery
  interEff <- bfs_c[4]/bfs_c[3] # drink * imagery
  
  effects <- c("drink", "imagery", "drink * imagery")
  BFincl <- c(exp(mainEff@bayesFactor[["bf"]][1]), exp(mainEff@bayesFactor[["bf"]][2]), exp(interEff@bayesFactor[["bf"]]))
  error <- c(mainEff@bayesFactor[["error"]][1], mainEff@bayesFactor[["error"]][2], interEff@bayesFactor[["error"]])
  effTbl <- data.frame(effects,BFincl,error)
  

  • These seem to match relatively well with the output in JASP, however, some effects seem to have a large discrepancy. For example, imagery in JASP had a BFincl of 1.32e14 whereas when I ran it in R I got a BFincl of 5.02e29. I understand there is an expectation for some variation but it seems quite large.
  • How many iterations would we need to run to get a better match between numerical accuracy? I am currently using 100,000. Is set seed = 1 equivalent for R and JASP?

  Also, how would we get the Analysis of Effects table using the 'Across all models' approach? Based on my understanding, it has something to do with summing across models but not sure what the actual order of operations is.

  Thanks,

  Corey