Skip to content

Instantly share code, notes, and snippets.

@achetverikov
Created October 21, 2021 19:44
Show Gist options
  • Save achetverikov/f8ad3c2cd2cd0c693faa5cf6ed9a99b7 to your computer and use it in GitHub Desktop.
Save achetverikov/f8ad3c2cd2cd0c693faa5cf6ed9a99b7 to your computer and use it in GitHub Desktop.
simulations with BayesFactors
library(BayesFactor)
library(future.apply)
library(future.batchtools)
# for the simulations, I used our cluster that runs PBS Torque
plan(batchtools_torque, resources = list(walltime = '01:00:00', memory = '2Gb')) # specify that the PBS Torque cluster is used and the resources we want
# create conditions combinations
conditions <- expand.grid(d = seq(0, 0.36, 0.04), N = seq(5, 65, 4))
# a function that gets BFs nrep times based on the conditions
get_bf <- function (curr_conditions, nrep = 10000){
# in each replication, we generate a random normal sample with a mean d and SD = 1 and compute the BF using default settings for one-sample t-test
bf <- replicate(nrep,
ttestBF(rnorm(curr_conditions[['N']], curr_conditions[['d']]))@bayesFactor$bf)
data.frame(N = curr_conditions[['N']], d = curr_conditions[['d']], bf = bf)
}
# run the simulations and combine the results into data.table
res <- rbindlist(future_apply(conditions, 1, get_bf, future.seed = T,
future.packages = c('data.table','BayesFactor')))
mean_bf <- res[,.(mean_log_bf = mean(bf), mean_bf = mean(exp(bf))), by = .(N,d)]
color_cuts <- c(0,.01,0.1,0.33,3,10,100,Inf)
ggplot(mean_bf, aes(x = d, y = N, fill = mean_bf))+geom_tile()+scale_fill_viridis_c(trans='log')+labs(fill = 'Mean logBF')
ggplot(mean_bf, aes(x = d, y = N, fill= cut2(mean_bf, color_cuts)))+geom_tile()+scale_fill_viridis_d()+labs(fill = 'Mean BF')
ggplot(mean_bf, aes(x = d, y = N, fill = mean_log_bf))+geom_tile()+scale_fill_viridis_c()+labs(fill = 'Mean logBF')
ggplot(mean_bf, aes(x = d, y = N, fill = cut2(mean_log_bf, cuts = log(color_cuts))))+geom_tile()+scale_fill_viridis_d()+labs(fill = 'Mean logBF')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment