JimGrange · November 29, 2023 08:02
diff --git a/gistfile1.txt b/gistfile1.txt
 # load packages -----------------------------------------------------------
 library(tidyverse)
 library(brms)
 library(tidybayes)
 library(bayesplot)
 library(emmeans)
 library(bayestestR)



 # import data -------------------------------------------------------------

 #--- ...although for now we're just making data up

 # for reproducibility
 set.seed(123)  

 # set number of subjects & number of trials per cell (here set to 100 for each)
 n_subjects <- 100
 n_trials <- 100 * 4

 # simulate trial-level data
 data <- tibble(
  crossing(
    id = 1:n_subjects, 
    trial = 1:n_trials
  )
 ) |> 
  mutate(congruency = sample(c("congruent", "incongruent"), 
                             n_subjects * n_trials, 
                             replace = TRUE), 
         cue = sample(c("informative", "noninformative"), 
                      n_subjects * n_trials, 
                      replace = TRUE)) |> 
  mutate(rt = rnorm(n_subjects * n_trials, 500, 60) + 
           rexp(n_subjects * n_trials, 1 / 100))

 # calculate aggregate data
 data <- data |> 
  group_by(id, congruency, cue) |> 
  summarise(rt = mean(rt))



 # fit the bayesian regressions --------------------------------------------

 # fit the interaction model
 int_model <- brm(rt ~ congruency * cue + (1|id), 
                 data = data,
                 iter = 5000, 
                 warmup = 2000,
                 chains = 4,
                 cores = 4,
                 seed = 123,
                 save_pars = save_pars(all = TRUE))

 # model summary & model checks
 summary(int_model)
 plot(int_model)
 pp_check(int_model, ndraws = 100)

 # fit the main effects model
 main_model <- brm(rt ~ congruency + cue + (1|id), 
                  data = data, 
                  iter = 5000, 
                  warmup = 2000,
                  chains = 4,
                  cores = 4,
                  seed = 123,
                  save_pars = save_pars(all = TRUE))

 # model summary & model checks
 summary(main_model)
 plot(main_model)
 pp_check(main_model, ndraws = 100)



 # model selection & posterior summary -------------------------------------

 #--- calculate bayes factor
 bayesfactor_models(int_model, denominator = main_model)


 #--- ROPE

 # extract the posterior of the interaction term in the interaction model
 int_posterior <- int_model |> 
  spread_draws(`b_congruencyincongruent:cuenoninformative`) |> 
  dplyr::select(`b_congruencyincongruent:cuenoninformative`) |> 
  pull()

 # calculate percentage of posterior within rope
 rope_min <- -0.2
 rope_max <- 0.2
 int_in_rope <- (sum(int_posterior > rope_min & int_posterior < rope_max) / 
                  length(int_posterior)) * 100



 # post-hocs ---------------------------------------------------------------

 #--- simple effects

 # main effect of congruency
 main_eff_cong <- pairs(emmeans(int_model, ~ congruency))

 # main effect of cue
 main_eff_cue <- pairs(emmeans(int_model, ~ cue))

 # interaction
 int <- emmeans(int_model, ~ congruency * cue)
 interaction <- pairs(int, by = "cue") 



 # robustness checks -------------------------------------------------------

 # etc....
	# load packages -----------------------------------------------------------
	library(tidyverse)
	library(brms)
	library(tidybayes)
	library(bayesplot)
	library(emmeans)
	library(bayestestR)



	# import data -------------------------------------------------------------

	#--- ...although for now we're just making data up

	# for reproducibility
	set.seed(123)

	# set number of subjects & number of trials per cell (here set to 100 for each)
	n_subjects <- 100
	n_trials <- 100 * 4

	# simulate trial-level data
	data <- tibble(
	crossing(
	id = 1:n_subjects,
	trial = 1:n_trials
	)
	) \|>
	mutate(congruency = sample(c("congruent", "incongruent"),
	n_subjects * n_trials,
	replace = TRUE),
	cue = sample(c("informative", "noninformative"),
	n_subjects * n_trials,
	replace = TRUE)) \|>
	mutate(rt = rnorm(n_subjects * n_trials, 500, 60) +
	rexp(n_subjects * n_trials, 1 / 100))

	# calculate aggregate data
	data <- data \|>
	group_by(id, congruency, cue) \|>
	summarise(rt = mean(rt))



	# fit the bayesian regressions --------------------------------------------

	# fit the interaction model
	int_model <- brm(rt ~ congruency * cue + (1\|id),
	data = data,
	iter = 5000,
	warmup = 2000,
	chains = 4,
	cores = 4,
	seed = 123,
	save_pars = save_pars(all = TRUE))

	# model summary & model checks
	summary(int_model)
	plot(int_model)
	pp_check(int_model, ndraws = 100)

	# fit the main effects model
	main_model <- brm(rt ~ congruency + cue + (1\|id),
	data = data,
	iter = 5000,
	warmup = 2000,
	chains = 4,
	cores = 4,
	seed = 123,
	save_pars = save_pars(all = TRUE))

	# model summary & model checks
	summary(main_model)
	plot(main_model)
	pp_check(main_model, ndraws = 100)



	# model selection & posterior summary -------------------------------------

	#--- calculate bayes factor
	bayesfactor_models(int_model, denominator = main_model)


	#--- ROPE

	# extract the posterior of the interaction term in the interaction model
	int_posterior <- int_model \|>
	spread_draws(`b_congruencyincongruent:cuenoninformative`) \|>
	dplyr::select(`b_congruencyincongruent:cuenoninformative`) \|>
	pull()

	# calculate percentage of posterior within rope
	rope_min <- -0.2
	rope_max <- 0.2
	int_in_rope <- (sum(int_posterior > rope_min & int_posterior < rope_max) /
	length(int_posterior)) * 100



	# post-hocs ---------------------------------------------------------------

	#--- simple effects

	# main effect of congruency
	main_eff_cong <- pairs(emmeans(int_model, ~ congruency))

	# main effect of cue
	main_eff_cue <- pairs(emmeans(int_model, ~ cue))

	# interaction
	int <- emmeans(int_model, ~ congruency * cue)
	interaction <- pairs(int, by = "cue")



	# robustness checks -------------------------------------------------------

	# etc....