bf_sims.R

library(tidyverse)
library(BFpack)

# function to simulate bayes factors in Kat's experiments
# takes:
#  - n subjects in each condition
#  - p1 for first condition success probability
#  - p2 for second condition success probability
sim_bfs <- function(n, p1, p2) {
  
  # simulate data from p1, p2
  data <- tibble(condition = c(rep("reliable", n), 
                                   rep("unreliable", n)),
                     harder = c(rbernoulli(n, p = p1), 
                                rbernoulli(n, p = p2)))

  # run a GLM and use the BFpack package to convert to a directional hypothesis test
  mod <- glm(harder ~ condition, family = "binomial", data = data)
  bfs <- BF(mod, hypothesis = "condition > 0")
  
  # get a contingency table so we can use the BayesFactor package to do a contingency table BF test
  d <- data %>% 
    group_by(condition) %>%
    summarise(sum_easier = sum(!harder), 
              sum_harder = sum(harder))
  data_mat <- as.matrix(d[,2:3])
  
  bfs2 <- BayesFactor::contingencyTableBF(x = data_mat,
                                  sampleType = "indepMulti", 
                                  fixedMargin = "rows")
  
  # extract each of these and return them as a tibble
  return(tibble(bf_bfpack = bfs$BFmatrix_confirmatory[1,2], 
                bf_bayesfactor = BayesFactor::extractBF(bfs2)$bf,
                p1_hat = sum(kat_data$harder[kat_data$condition == "reliable"]),
                p2_hat = sum(kat_data$harder[kat_data$condition == "unreliable"])))
}

sims <- expand_grid(nsim = 1:100, 
                    p1 = .5,
                    p2 = c(.5,.6,.7,.8,.9),  
                    n = c(10,20,30)) %>%
  rowwise() %>%
  mutate(bf = sim_bfs(n, p1, p2)) %>% 
  unnest(cols = c(bf))

ggplot(data = sims, 
       aes(x = log(bf_bayesfactor))) +
  geom_histogram() + 
  facet_grid(n~p2)+
  geom_vline(col = "black", xintercept = 0, lty = 2) + 
  geom_vline(col = "green", xintercept = log(10), lty = 2) +
  geom_vline(col = "red", xintercept = log(.333), lty = 2)