poole_cuing_analysis.R

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....