Weissman-replication-analysis.Rmd

---
title: "Weissman-replication-analysis"
author: "Marton Kovacs"
date: '`r format(Sys.time(), "%d %B, %Y")`'
output: html_document
---

# Load packages

```{r load packages, warning = FALSE, message = FALSE}
if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse, osfr, lme4, BayesFactor, DescTools, sjstats, car, broom, papaja, cowplot)
```

# Load helper functions

```{r load helper functions}
source("utils.R")
```

# Download data from OSF
## OSF auth (until project is public)

Read/write token should be created on OSF and it should be named as "osf_token_write"

```{r osf authentication}
osf_auth(token = read_lines(list.files(pattern ="osf_token_write_.*.txt")))
```

## Connect to data OSF folder
```{r osf connect to OSF folder, message = FALSE}
data_guid <- "9knds"

weissman_project <- osf_retrieve_node(data_guid)
```

## Download data locally

```{r osf download data, message = FALSE, echo = T, results = "hide"}
local_data_pth <- file.path("data","Processed")

create_local_structure(local_data_pth)

data_files <-
  weissman_project %>%
  osf_ls_files() %>%
  filter(name == "Processed") %>%
  osf_ls_files()

data_files %>%
  group_by(name) %>% # for each experiment type
  do(download_files(.,local_data_pth))

# uncomment following line to remove the data
# remove_local_data(local_data_pth)
```

# Import data

```{r, message = FALSE, warning = FALSE}
# Saving subfolder names
subfolder <- list("primeprobe", "flanker", "stroop", "simon")

# Reading data
processed <-
  tibble(task = subfolder,
         rt_data = map(subfolder,
                       ~ read_plus(subfolder_name = .x,
                                   pattern = ".tsv$",
                                   path = "data/Processed/",
                                   sep = "\t",
                                   include = "_Rt_")),
         acc_data = map(subfolder,
                        ~ read_plus(subfolder_name = .x,
                                    pattern = ".tsv$",
                                    path = "data/Processed/",
                                    sep = "\t",
                                    include = "_Acc_")),
         acc_exp_data = map(subfolder,
                        ~ read_plus(subfolder_name = .x,
                                    pattern = ".tsv$",
                                    path = "data/Processed/",
                                    sep = "\t",
                                    include = "Acc_Exploratory_")))

# Reading original study effect sizes
original <- read_tsv("original_study_data.tsv")
```

# Set task names

```{r}
processed$rt_data <- set_names(processed$rt_data, processed$task)

processed$acc_data <- set_names(processed$acc_data, processed$task)

processed$acc_exp_data <- set_names(processed$acc_exp_data, processed$task)
```

# Descriptive statistics of the sample size
## Number of participants left after exclusion

The number of participants is the same for the reaction time, the accuracy, and the exploratory analyses as well

```{r}
map(processed$rt_data,
    . %>% 
      distinct(participant_id) %>% 
      count())
```

## Number of responses after exclusion
### For the reaction time analysis

```{r}
map(processed$rt_data,
    . %>% 
      count())
```

### For the accuracy analysis

```{r}
map(processed$acc_data,
    . %>% 
      count())
```

### For the accuracy exploratory analysis

```{r}
map(processed$acc_exp_data,
    . %>% 
      count())
```

## Age

We calculated the following statistics based on the reaction time data file, but because it contains the same participants as the accuracy data, and the exploratory data files there should be no difference between the results.

```{r}
map_df(processed$rt_data,
    . %>% 
      distinct(participant_id, .keep_all = T) %>% 
      group_by(task) %>% 
      summarize(median_age = median(age, na.rm = T),
                min_age = min(age, na.rm = T),
                max_age = max(age, na.rm = T),
                quart1_age = quantile(age, probs = 0.25, na.rm = T),
                quart3_age = quantile(age, probs = 0.75, na.rm = T))) %>% 
      
      knitr::kable(caption = "Age demographics")
```

## Gender

```{r}
map_df(processed$rt_data,
    . %>% 
      distinct(participant_id, .keep_all = T) %>% 
      group_by(task,gender) %>% 
      count() %>% 
      ungroup() %>% 
      filter(gender %in% c("female", "male")) %>% 
      mutate(prop = n / sum(n) * 100)) %>% 
      knitr::kable(caption = "Gender demographics")
```

## Language

```{r}
map_df(processed$rt_data,
    . %>% 
      distinct(participant_id, .keep_all = T) %>% 
      mutate(loc = case_when(loc == "HUN" ~ "Hungarian",
                             loc == "CZ" ~ "Czech")) %>% 
      group_by(task,loc) %>% 
      count() %>% 
      ungroup() %>% 
      mutate(prop = n / sum(n) * 100)) %>% 
      knitr::kable(caption = "Language demographics")
```

## Education demographics

```{r}
map_df(processed$rt_data,
    . %>% 
      distinct(participant_id, .keep_all = T) %>% 
      group_by(task,education) %>% 
      count() %>% 
      ungroup() %>% 
      filter(education %ni% c("blank", "refused")) %>% 
      mutate(prop = n / sum(n) * 100)) %>% 
      knitr::kable(caption = "Education demographics")
```

# Figures of CSE
## For the reaction time data

```{r}
processed <-
  processed %>% 
  mutate(cse_plot_rt_data = map(rt_data,
                                . %>% 
                                  mutate(isPrevCongruent = case_when(isPrevCongruent ==  0L ~ "Incongruent",
                                                                     isPrevCongruent ==  1L ~ "Congruent"),
                                         isCongruent = case_when(isCongruent ==  0L ~ "Incongruent",
                                                                 isCongruent ==  1L ~ "Congruent")) %>%
                                  group_by(participant_id, isPrevCongruent, isCongruent) %>% 
                                  summarise(participant_mean_rt = mean(rt, na.rm = T)) %>% 
                                  group_by(isPrevCongruent, isCongruent) %>% 
                                  summarise(N = n(),
                                            mean_rt = mean(participant_mean_rt, na.rm = T),
                                            sd_rt = sd(participant_mean_rt, na.rm = T),
                                            se_rt = sd_rt / sqrt(N))),
         cse_plot_rt = map2(cse_plot_rt_data, task, 
                           ~ ggplot(.x, aes(x = isPrevCongruent, y = mean_rt, shape = isCongruent, group = isCongruent)) +
                             geom_path() +
                             geom_point(size = 2) +
                             geom_errorbar(aes(ymin = mean_rt - se_rt, ymax = mean_rt + se_rt), width=.1) +
                             scale_shape_manual(values=c(4, 16)) +
                             scale_y_continuous(limits = c(550, 1000)) +
                             scale_x_discrete(expand=c(1, 0)) +
                             ggtitle(stringr::str_to_title(.y)) +
                             xlab("Congruency of the previous trial")+
                             ylab("Reaction time") +
                             guides(shape = guide_legend(title="Congruency of \n the current trial")) +
                             papaja::theme_apa() +
                             theme(legend.position = c(0.85, 0.5),
                                   axis.line = element_line(color = "black"))))
```

## Creating figure for the paper

```{r}
cowplot::plot_grid(plotlist = processed$cse_plot_rt)
```

## Saving the figure

```{r}
ggsave("figures/weissman_replication_rt_cse.png", width = 14.4, height = 8, plot = last_plot())
```

## For the accuracy data

```{r}
processed <-
  processed %>% 
  mutate(cse_plot_acc_data = map(acc_data,
                                . %>% 
                                  mutate(isPrevCongruent = case_when(isPrevCongruent ==  0L ~ "Incongruent",
                                                                     isPrevCongruent ==  1L ~ "Congruent"),
                                         isCongruent = case_when(isCongruent ==  0L ~ "Incongruent",
                                                                 isCongruent ==  1L ~ "Congruent")) %>%
                                  group_by(participant_id, isPrevCongruent, isCongruent) %>% 
                                  summarise(participant_mean_acc = mean(isCorrect, na.rm = T)) %>% 
                                  group_by(isPrevCongruent, isCongruent) %>% 
                                  summarise(N = n(),
                                            mean_acc = mean(participant_mean_acc, na.rm = T),
                                            sd_acc = sd(participant_mean_acc, na.rm = T),
                                            se_acc = sd_acc / sqrt(N))),
         cse_plot_acc = map2(cse_plot_acc_data, task,
                           ~ ggplot(.x, aes(x = isPrevCongruent, y = mean_acc, shape = isCongruent, group = isCongruent)) +
                             geom_path() +
                             geom_point(size = 2) +
                             geom_errorbar(aes(ymin = mean_acc - se_acc, ymax = mean_acc + se_acc), width=.1) +
                             scale_shape_manual(values=c(4, 16)) + 
                             scale_y_continuous(limits = c(0.93, 1)) +
                             scale_x_discrete(expand=c(1, 0)) +
                             ggtitle(stringr::str_to_title(.y)) +
                             xlab("Congruency of the previous trial")+
                             ylab("Accuracy") +
                             guides(shape = guide_legend(title="Congruency of \n the current trial")) +
                             papaja::theme_apa() +
                             theme(legend.position = c(0.85, 0.5),
                                   axis.line = element_line(color = "black"))))
```

## Creating figure for the paper

```{r}
cowplot::plot_grid(plotlist = processed$cse_plot_acc)
```

## Saving the figure

```{r}
ggsave("figures/weissman_replication_acc_cse.png", width = 14.4, height = 8, plot = last_plot())
```

## For the exploratory accuracy data

```{r}
processed <-
  processed %>% 
  mutate(cse_plot_acc_exp_data = map(acc_exp_data,
                                . %>% 
                                  mutate(isPrevCongruent = case_when(isPrevCongruent ==  0L ~ "Incongruent",
                                                                     isPrevCongruent ==  1L ~ "Congruent"),
                                         isCongruent = case_when(isCongruent ==  0L ~ "Incongruent",
                                                                 isCongruent ==  1L ~ "Congruent")) %>%
                                  group_by(participant_id, isPrevCongruent, isCongruent) %>% 
                                  summarise(participant_mean_acc = mean(isCorrect, na.rm = T)) %>% 
                                  group_by(isPrevCongruent, isCongruent) %>% 
                                  summarise(N = n(),
                                            mean_acc = mean(participant_mean_acc, na.rm = T),
                                            sd_acc = sd(participant_mean_acc, na.rm = T),
                                            se_acc = sd_acc / sqrt(N))),
         cse_plot_acc_exp = map2(cse_plot_acc_exp_data, task,
                           ~ ggplot(.x, aes(x = isPrevCongruent, y = mean_acc, shape = isCongruent, group = isCongruent)) +
                             geom_path() +
                             geom_point(size = 2) +
                             geom_errorbar(aes(ymin = mean_acc - se_acc, ymax = mean_acc + se_acc), width=.1) +
                             scale_shape_manual(values=c(4, 16)) + 
                             scale_y_continuous(limits = c(0.93, 1)) +
                             scale_x_discrete(expand=c(1, 0)) +
                             ggtitle(stringr::str_to_title(.y)) +
                             xlab("Congruency of the previous trial")+
                             ylab("Accuracy") +
                             guides(shape = guide_legend(title="Congruency of \n the current trial")) +
                             papaja::theme_apa() +
                             theme(legend.position = c(0.85, 0.5),
                                   axis.line = element_line(color = "black"))))
```

## Creating figure for the paper

```{r}
cowplot::plot_grid(plotlist = processed$cse_plot_acc_exp)
```

## Saving the figure

```{r}
ggsave("figures/weissman_replication_acc_exp_cse.png", width = 14.4, height = 8, plot = last_plot())
```

# Reaction time analyis
## I. analysis: Testing whether the CSE was present in the different tasks
### Preprocessing data

We are calculating mean reaction time for the analysis for each participant in each condition for each task.

```{r calculating mean dev rt}
processed <-
  processed %>% 
  mutate(rt_data = map(rt_data,
                       . %>% 
                         group_by(participant_id, isCongruent, isPrevCongruent) %>% 
                         summarise(rtConditionalMean = mean(rt, na.rm = T)) %>%
                         ungroup() %>% 
                         mutate(isCongruent = as.factor(isCongruent),
                                isPrevCongruent = as.factor(isPrevCongruent),
                                participant_id = as.factor(participant_id))))
```

We are calculating the main congruency effect of the current (n) trial, the main congruency effect of the previous trial (n-1) and the interaction effect (CSE) between them for each participant.

```{r}
processed <-
  processed %>% 
  mutate(rt_effect_data = map(rt_data,
                              . %>% 
                                ungroup() %>% 
                                mutate(condition = case_when(isPrevCongruent == 0L & isCongruent == 0L ~ "ii",
                                                             isPrevCongruent == 0L & isCongruent == 1L ~ "ic",
                                                             isPrevCongruent == 1L & isCongruent == 0L ~ "ci",
                                                             isPrevCongruent == 1L & isCongruent == 1L ~ "cc",
                                                             TRUE ~ NA_character_)) %>% 
                                select(-isPrevCongruent, -isCongruent) %>% 
                                spread(key = condition, value = rtConditionalMean) %>% 
                                mutate(congruencyEffect = ((ci + ii) / 2) - ((cc + ic) / 2),
                                       prevCongruencyEffect = ((ci + cc) / 2) - ((ii + ic) / 2),
                                       cseEffect = (ci - cc) - (ii - ic))))
```

### Running the analysis

We are running a 2 × 2 repeated-measures ANOVA per task with mean RT as dependent variable. The two factors of the ANOVA are Previous Trial Congruency (congruent, incongruent) and Current Trial Congruency (congruent, incongruent).

```{r}
processed <-
  processed %>% 
  mutate(anova_rt = map(rt_data,
                        ~ aov(rtConditionalMean ~ isCongruent * isPrevCongruent + Error(participant_id / (isCongruent * isPrevCongruent)), data = .x)))

# Print results with names
processed$anova_rt %>% 
  set_names(., processed$task) %>% 
  map(., summary)
```

The interaction was only significant for the first task, the primeprobe task.

We are calculating the partial eta square as an effect size estimate for the ANOVA.

```{r}
processed <-
  processed %>% 
  mutate(anova_rt_eta = map(anova_rt,
                            ~ eta_sq(.x, partial = TRUE)))

# Print results with names
processed$anova_rt_eta %>% 
  set_names(., processed$task) %>% 
  map(., summary)
```

### Calculating Bayes factors
#### Preprocessing data for the Bayes factor analysis

We are substracting the F values for each effect (main congruency effect of current trial, main congruency effect of previous trial and interaction effect) from the ANOVA.

```{r}
processed <-
  processed %>% 
  mutate(anova_rt_f = map(anova_rt,
                          . %>%
                            broom::tidy() %>% 
                            select(term, statistic) %>% 
                            transmute(term = case_when(term == "isCongruent" ~ "mainEffectIsCongruent",
                                                       term == "isPrevCongruent" ~ "mainEffectIsPrevCongruent",
                                                       term == "isCongruent:isPrevCongruent" ~ "interactionEffect",
                                                       term == "Residuals" ~ term),
                                      fValue = statistic) %>% 
                            filter(!is.na(fValue))))
```

We are calculating the mean effects (main congruency effect of current trial, main congruency effect of previous trial and interaction effect) for each task.

```{r}
processed <-
  processed %>% 
  mutate(anova_rt_raw_effect = map(rt_effect_data, 
                                   . %>% 
                                     summarise(mainEffectIsCongruent = mean(congruencyEffect, na.rm = T),
                                               mainEffectIsPrevCongruent = mean(prevCongruencyEffect, na.rm = T),
                                               interactionEffect = mean(cseEffect, na.rm = T),
                                               nParticipant = n()) %>%
                                     gather(key = "term", value = "rawEffect", -nParticipant)))
```

We are adding the calucalted F-values for each effect to the calculated mean effects.

```{r}
processed <-
  processed %>% 
  mutate(anova_rt_raw_effect = map2(anova_rt_raw_effect, anova_rt_f,
                                    ~ inner_join(.x, .y, by = "term")))
```

Calculating the SEs.

```{r}
processed <- 
  processed %>% 
  mutate(anova_rt_raw_effect = map(anova_rt_raw_effect,
                                   . %>% 
                                     mutate(se = rawEffect / sqrt(fValue))))
```

Preprocessing effect sizes that Weissman et al. (2014) found for the H1 model of the Bayes factor analysis.

```{r}
original_rt_effect <-
  original %>% 
  filter(measure == "rt") %>% 
  select(task, congruency, cse) %>% 
  gather(key = "term", value = "originalRawEffect", -task) %>%
  mutate(term = case_when(term == "congruency" ~ "mainEffectIsCongruent",
                          term == "cse" ~ "interactionEffect")) %>% 
  group_by(task) %>% 
  nest() %>% 
  rename(original_effect_reaction_time = data)
```

Joining the effect sizes of the replicated study with the effect sizes of the current study.

```{r}
processed <-
  processed %>%
  mutate(task = as.character(task)) %>% # vector needs to be character class (instead of a list) for joining
  left_join(., original_rt_effect, by = "task") %>% 
  mutate(anova_rt_bf_effect = map2(anova_rt_raw_effect, original_effect_reaction_time,
                              ~ left_join(.x, .y, by = "term")))
```

#### Calculating the Bayes factor
##### For the interaction effect

A half-normal distribution is used with a mode of 0 and a SD equal to the half of the congruency effect (in ms) for the particular task in the original data set by Weissman et al. (2014) to model the prediction of the interaction effect. 

```{r}
processed <-
  processed %>%
  mutate(bf_interaction_rt = map(anova_rt_bf_effect,
                                 ~ Bf(sd = .x$se[3],
                                      obtained = .x$rawEffect[3],
                                      dfdata = .x$nParticipant[1] - 1,
                                      meanoftheory = 0,
                                      sdtheory = .x$originalRawEffect[1] / 2,
                                      dftheory = 10^10,
                                      tail = 1)))

# Print results with names
processed$bf_interaction_rt %>% 
  set_names(., processed$task)
```

###### Finding the robustness region for the interaction effect

```{r}
# Finding the lower boundary of the robustness region
## Primeprobe
primeprobe_rt <- processed$anova_rt_bf_effect[[1]]

Bf(sd = primeprobe_rt$se[3],
   obtained = primeprobe_rt$rawEffect[3],
   dfdata = primeprobe_rt$nParticipant[1] - 1,
   meanoftheory = 0,
   sdtheory = 0.60,
   dftheory = 10^10,
   tail = 1)

## Flanker
flanker_rt <- processed$anova_rt_bf_effect[[2]]

Bf(sd = flanker_rt$se[3],
   obtained = flanker_rt$rawEffect[3],
   dfdata = flanker_rt$nParticipant[1] - 1,
   meanoftheory = 0,
   sdtheory = 14.6,
   dftheory = 10^10,
   tail = 1)

## Stroop
stroop_rt <- processed$anova_rt_bf_effect[[3]]

### The lower limit is 0, because the Bf is inconslusive

## Simon
simon_rt <- processed$anova_rt_bf_effect[[4]]

## The lower limit is 0, because the Bf is inconslusive

# Finding the upper boundary of the robustness region
## Primeprobe
Bf(sd = primeprobe_rt$se[3],
   obtained = primeprobe_rt$rawEffect[3],
   dfdata = primeprobe_rt$nParticipant[1] - 1,
   meanoftheory = 0,
   sdtheory = 6423,
   dftheory = 10^10,
   tail = 1)

## Flanker

### The upper boundary is infinite, because the Bf supports the null model

## Stroop
Bf(sd = stroop_rt$se[3],
   obtained = stroop_rt$rawEffect[3],
   dfdata = stroop_rt$nParticipant[1] - 1,
   meanoftheory = 0,
   sdtheory = 222,
   dftheory = 10^10,
   tail = 1)

## Simon
Bf(sd = simon_rt$se[3],
   obtained = simon_rt$rawEffect[3],
   dfdata = simon_rt$nParticipant[1] - 1,
   meanoftheory = 0,
   sdtheory = 120,
   dftheory = 10^10,
   tail = 1)
```

##### For the main congruency effect of the current trial

To test the congruency main effect, the SD of the H1 model is set to the congruency effect reported by Weissman et al. (2014) for the given task. The other settings stay the same.

```{r}
processed <-
  processed %>%
  mutate(bf_congruency_rt = map(anova_rt_bf_effect,
                                ~ Bf(sd = .x$se[1],
                                     obtained = .x$rawEffect[1],
                                     dfdata = .x$nParticipant[1] - 1,
                                     meanoftheory = 0,
                                     sdtheory = .x$originalRawEffect[1],
                                     dftheory = 10^10,
                                     tail = 1)))

# Print results with names
processed$bf_congruency_rt %>% 
  set_names(., processed$task)
```

###### Finding the robustness region for the interaction effect

```{r}
# Finding the lower boundary of the robustness region
## Primeprobe
Bf(sd = primeprobe_rt$se[1],
   obtained = primeprobe_rt$rawEffect[1],
   dfdata = primeprobe_rt$nParticipant[1] - 1,
   meanoftheory = 0,
   sdtheory = 0.74,
   dftheory = 10^10,
   tail = 1)

## Flanker
Bf(sd = flanker_rt$se[1],
   obtained = flanker_rt$rawEffect[1],
   dfdata = flanker_rt$nParticipant[1] - 1,
   meanoftheory = 0,
   sdtheory = 1,
   dftheory = 10^10,
   tail = 1)

## Stroop
Bf(sd = stroop_rt$se[1],
   obtained = stroop_rt$rawEffect[1],
   dfdata = stroop_rt$nParticipant[1] - 1,
   meanoftheory = 0,
   sdtheory = 1.75,
   dftheory = 10^10,
   tail = 1)

## Simon
Bf(sd = simon_rt$se[1],
   obtained = simon_rt$rawEffect[1],
   dfdata = simon_rt$nParticipant[1] - 1,
   meanoftheory = 0,
   sdtheory = 0.74,
   dftheory = 10^10,
   tail = 1)

# Finding the upper boundary of the robustness region
## Primeprobe
Bf(sd = primeprobe_rt$se[1],
   obtained = primeprobe_rt$rawEffect[1],
   dfdata = primeprobe_rt$nParticipant[1] - 1,
   meanoftheory = 0,
   sdtheory = 26512,
   dftheory = 10^10,
   tail = 1)

## Flanker
Bf(sd = flanker_rt$se[1],
   obtained = flanker_rt$rawEffect[1],
   dfdata = flanker_rt$nParticipant[1] - 1,
   meanoftheory = 0,
   sdtheory = 17181,
   dftheory = 10^10,
   tail = 1)

## Stroop
Bf(sd = stroop_rt$se[1],
   obtained = stroop_rt$rawEffect[1],
   dfdata = stroop_rt$nParticipant[1] - 1,
   meanoftheory = 0,
   sdtheory = 47403,
   dftheory = 10^10,
   tail = 1)

## Simon
Bf(sd = simon_rt$se[1],
   obtained = simon_rt$rawEffect[1],
   dfdata = simon_rt$nParticipant[1] - 1,
   meanoftheory = 0,
   sdtheory = 15000,
   dftheory = 10^10,
   tail = 1)
```

### Follow up simple effects analysis
#### Preprocessing data

Saving the data of the task (primeprobe) with significant interaction between current trial congruency and previous trial congruency.

```{r}
primeprobe <- processed$rt_effect_data[[1]]
```

#### Checking the number of responses

```{r}
primeprobe %>% 
  count()
```

#### Running the analysis

If the Bayes factor is greater than the preset threshold of B = 3 for the interaction effect in a given task we run a simple effects analyses, contrasting current-congruent congruent (cC) and current-incongruent congruent (iC) trials.

```{r current congruent follow up comparison for rt}
current_congruent <- 
  t.test(primeprobe$ic, primeprobe$cc, paired = TRUE) %>% 
  broom::tidy()

current_congruent %>% knitr::kable()
```

current-congruent incongruent (cI) and current-incongruent incongruent (iI) trials in another.

```{r current incongruent follow up comparison for rt}
current_incongruent <- 
  t.test(primeprobe$ci, primeprobe$ii, paired = TRUE) %>% 
  broom::tidy()

current_incongruent %>% knitr::kable()
```

#### Calculating the Bayes factor for the follow up analysis

The prior H1 model is a half-normal distribution with a mode of 0 and an SD equal to half of the CSE estimate for that given task, i.e., the originally reported congruency effect divided by four.

##### Current congruent follow up analysis

```{r}
Bf(sd = current_congruent$estimate / current_congruent$statistic,
   obtained = current_congruent$estimate,
   dfdata = current_congruent$parameter,
   meanoftheory = 0,
   sdtheory = processed$anova_rt_bf_effect[[1]]$originalRawEffect[1] / 4,
   dftheory = 10^10,
   tail = 1)
```

###### Finding the robustness region for the current congruent follow up analysis

```{r}
# Finding the lower boundary of the robustness region
Bf(sd = current_congruent$estimate / current_congruent$statistic,
   obtained = current_congruent$estimate,
   dfdata = current_congruent$parameter,
   meanoftheory = 0,
   sdtheory = 0.525,
   dftheory = 10^10,
   tail = 1)
```

```{r}
# Finding the upper boundary of the robustness region
Bf(sd = current_congruent$estimate / current_congruent$statistic,
   obtained = current_congruent$estimate,
   dfdata = current_congruent$parameter,
   meanoftheory = 0,
   sdtheory = 2713,
   dftheory = 10^10,
   tail = 1)
```

##### Current incongruent follow up analysis

```{r}
Bf(sd = current_incongruent$estimate / current_incongruent$statistic,
   obtained = current_incongruent$estimate,
   dfdata = current_incongruent$parameter,
   meanoftheory = 0,
   sdtheory = processed$anova_rt_bf_effect[[1]]$originalRawEffect[1] / 4,
   dftheory = 10^10,
   tail = 1)
```

###### Finding the robustness region for the current incongruent follow up analysis

```{r}
# Finding the lower boundary of the robustness region
Bf(sd = current_incongruent$estimate / current_incongruent$statistic,
   obtained = current_incongruent$estimate,
   dfdata = current_incongruent$parameter,
   meanoftheory = 0,
   sdtheory = 0.53,
   dftheory = 10^10,
   tail = 1)
```

```{r}
# Finding the upper boundary of the robustness region
Bf(sd = current_incongruent$estimate / current_incongruent$statistic,
   obtained = current_incongruent$estimate,
   dfdata = current_incongruent$parameter,
   meanoftheory = 0,
   sdtheory = 2973,
   dftheory = 10^10,
   tail = 1)
```

## II. analysis: Testing whether the size of the CSE varies across tasks
### Creating the datatable

We are combining the effects (main congruency effect of current trial, main congruency effect of previous trial and interaction effect) for each task (primeprobe, Flanker, Stroop, Simon) in one datatable.

```{r, warning = FALSE}
interaction_rt_data <- 
  processed %>% 
  select(task, rt_effect_data) %>% 
  unnest(rt_effect_data) %>% 
  mutate(task = as.factor(task))
```

### Testing the assumptions

We are testing the homogenity of variances.

```{r}
car::leveneTest(cseEffect ~ task, data = interaction_rt_data)
```

The assumption of homogenity of variances is violated. Therefore, we are running a Kruskal Wallis H test instead of an ANOVA. 

### Running the analysis

We are running a  Kruskal-Wallis H test with task (levels: Prime-Probe, Flanker, Stroop, Simon) as the single factor, and the CSE index as the dependent variable to examine if the size of the effect differs across tasks.

```{r}
kruskal.test(cseEffect ~ task, data = interaction_rt_data)
```

### Calculating pairwise comparisons

A Dunn’s test is used for pairwise comparisons.

```{r}
dunn_test <- DescTools::DunnTest(cseEffect ~ task, data = interaction_rt_data)

dunn_test
```

Creating a tidy datatable from the Dunn test output.

```{r}
dunn_test <-
  dunn_test[[1]] %>% 
  as_tibble(rownames = "comparison") %>% 
  janitor::clean_names()
```

#### Bf analysis of pairwise comparisons
##### Preprocessing data for the Bf analyses

```{r}
# Stroop and Primeprobe comparison
stroop_primeprobe_comparison <- 
  t.test(cseEffect ~ task,
         data = filter(interaction_rt_data,
                       task %in% c("stroop", "primeprobe"))) %>% 
  broom::tidy()

# Stroop and Simon comparison
stroop_simon_comparison <- 
  t.test(cseEffect ~ task,
         data = filter(interaction_rt_data,
                       task %in% c("stroop", "simon"))) %>% 
  broom::tidy()

# Stroop and Flanker comparison
stroop_flanker_comparison <- 
  t.test(cseEffect ~ task,
         data = filter(interaction_rt_data,
                       task %in% c("stroop", "flanker"))) %>% 
  broom::tidy()

# Primeprobe and Flanker comparison
primeprobe_flanker_comparison <- 
  t.test(cseEffect ~ task,
         data = filter(interaction_rt_data,
                       task %in% c("primeprobe", "flanker"))) %>% 
  broom::tidy()

# Primeprobe and Simon comparison
primeprobe_simon_comparison <- 
  t.test(cseEffect ~ task,
         data = filter(interaction_rt_data,
                       task %in% c("primeprobe", "simon"))) %>% 
  broom::tidy()

# Flanker and Simon comparison
flanker_simon_comparison <- 
  t.test(cseEffect ~ task,
         data = filter(interaction_rt_data,
                       task %in% c("flanker", "simon"))) %>% 
  broom::tidy()
```

##### Running the analysis

A half-normal distribution is used with a mode of 0 and SD equal to half of the greater congruency effect of the two in any given pair as reported by Weissman et al. (2014).

Tasks in a descending order based on their congruency effect according to the original study:
Stroop > Primeprobe > Flanker > Simon

We are using a two-tailed test for the analysis.

###### Stroop and Primeprobe comparison

```{r}
# Bf Analysis
Bf(obtained = stroop_primeprobe_comparison$estimate,
   sd = stroop_primeprobe_comparison$estimate / stroop_primeprobe_comparison$statistic,
   dfdata = stroop_primeprobe_comparison$parameter,
   meanoftheory = 0,
   sdtheory = original$congruency[1] / 2,
   dftheory = 10^10,
   tail = 2)

# Finding the robustness region for the current incongruent follow up analysis
## Finding the lower boundary of the robustness region
Bf(obtained = stroop_primeprobe_comparison$estimate,
   sd = stroop_primeprobe_comparison$estimate / stroop_primeprobe_comparison$statistic,
   dfdata = stroop_primeprobe_comparison$parameter,
   meanoftheory = 0,
   sdtheory = 28.8,
   dftheory = 10^10,
   tail = 2)

## Finding the upper boundary of the robustness region
### The upper boundary is infinite, because the Bf supports the null model
```

###### Stroop and Simon comparison

```{r}
# Bf Analysis
Bf(obtained = stroop_simon_comparison$estimate,
   sd = stroop_simon_comparison$estimate / stroop_simon_comparison$statistic,
   dfdata = stroop_simon_comparison$parameter,
   meanoftheory = 0,
   sdtheory = original$congruency[1] / 2,
   dftheory = 10^10,
   tail = 2)

# Finding the robustness region for the current incongruent follow up analysis
## Finding the lower boundary of the robustness region
Bf(obtained = stroop_simon_comparison$estimate,
   sd = stroop_simon_comparison$estimate / stroop_simon_comparison$statistic,
   dfdata = stroop_simon_comparison$parameter,
   meanoftheory = 0,
     sdtheory = 37.7,
   dftheory = 10^10,
   tail = 2)

## Finding the upper boundary of the robustness region
## The upper boundary is infinite, because the Bf supports the null model
```

###### Stroop and Flanker comparison

```{r}
# Bf Analysis
Bf(obtained = stroop_flanker_comparison$estimate,
   sd = stroop_flanker_comparison$estimate / stroop_flanker_comparison$statistic,
   dfdata = stroop_flanker_comparison$parameter,
   meanoftheory = 0,
   sdtheory = original$congruency[1] / 2,
   dftheory = 10^10,
   tail = 2)

# Finding the robustness region for the current incongruent follow up analysis
## Finding the lower boundary of the robustness region
### The lower limit is 0, because the Bf is inconslusive

## Finding the upper boundary of the robustness region
Bf(obtained = stroop_flanker_comparison$estimate,
   sd = stroop_flanker_comparison$estimate / stroop_flanker_comparison$statistic,
   dfdata = stroop_flanker_comparison$parameter,
   meanoftheory = 0,
   sdtheory = 140,
   dftheory = 10^10,
   tail = 2)
```

###### Primeprobe and Flanker comparison

```{r}
# Bf Analysis
Bf(obtained = primeprobe_flanker_comparison$estimate,
   sd = primeprobe_flanker_comparison$estimate / primeprobe_flanker_comparison$statistic,
   dfdata = primeprobe_flanker_comparison$parameter,
   meanoftheory = 0,
   sdtheory = original$congruency[4] / 2,
   dftheory = 10^10,
   tail = 2)

# Finding the robustness region for the current incongruent follow up analysis
## Finding the lower boundary of the robustness region
Bf(obtained = primeprobe_flanker_comparison$estimate,
   sd = primeprobe_flanker_comparison$estimate / primeprobe_flanker_comparison$statistic,
   dfdata = primeprobe_flanker_comparison$parameter,
   meanoftheory = 0,
   sdtheory = 6.3,
   dftheory = 10^10,
   tail = 2)

## Finding the upper boundary of the robustness region
Bf(obtained = primeprobe_flanker_comparison$estimate,
   sd = primeprobe_flanker_comparison$estimate / primeprobe_flanker_comparison$statistic,
   dfdata = primeprobe_flanker_comparison$parameter,
   meanoftheory = 0,
   sdtheory = 101,
   dftheory = 10^10,
   tail = 2)
```

###### Primeprobe and Simon comparison

```{r}
# Bf Analysis
Bf(obtained = primeprobe_simon_comparison$estimate,
   sd = primeprobe_simon_comparison$estimate / primeprobe_simon_comparison$statistic,
   dfdata = primeprobe_simon_comparison$parameter,
   meanoftheory = 0,
   sdtheory = original$congruency[4] / 2,
   dftheory = 10^10,
   tail = 2)

# Finding the robustness region for the current incongruent follow up analysis
## Finding the lower boundary of the robustness region
Bf(obtained = primeprobe_simon_comparison$estimate,
   sd = primeprobe_simon_comparison$estimate / primeprobe_simon_comparison$statistic,
   dfdata = primeprobe_simon_comparison$parameter,
   meanoftheory = 0,
   sdtheory = 35,
   dftheory = 10^10,
   tail = 2)

## Finding the upper boundary of the robustness region
## The upper boundary is infinite, because the Bf supports the null model
```

###### Flanker and Simon comparison

```{r}
# Bf Analysis
Bf(obtained = flanker_simon_comparison$estimate,
   sd = flanker_simon_comparison$estimate / flanker_simon_comparison$statistic,
   dfdata = flanker_simon_comparison$parameter,
   meanoftheory = 0,
   sdtheory = original$congruency[3] / 2,
   dftheory = 10^10,
   tail = 2)

# Finding the robustness region for the current incongruent follow up analysis
## Finding the lower boundary of the robustness region
### The lower limit is 0, because the Bf is inconslusive

## Finding the upper boundary of the robustness region
Bf(obtained = flanker_simon_comparison$estimate,
   sd = flanker_simon_comparison$estimate / flanker_simon_comparison$statistic,
   dfdata = flanker_simon_comparison$parameter,
   meanoftheory = 0,
   sdtheory = 83,
   dftheory = 10^10,
   tail = 2)
```

## III. analysis: Testing whether the size of the CSE changes as a function of the magnitude of interference

We are calculating correlations between the CSE index and the congruency effect for each task.

```{r}
processed <-
  processed %>% 
  mutate(corr = map(rt_effect_data,
                    ~ cor.test(.x$congruencyEffect, .x$cseEffect, method = "pearson") %>% 
                      broom::tidy()))

# Print the results with names
processed$corr %>% 
  set_names(., processed$task) %>% 
  knitr::kable()
```

### Calculating Bayes factors
#### Preprocessing data for the Bayes factor analysis

Preprocessing correlation coefficients that Weissman et al. (2014) found for the H1 model of the Bayes factor analysis.

```{r}
original_corr_effect <- 
  original %>% 
  filter(measure == "rt") %>% 
  select(task, corr) %>% 
  rename(original_corr_estimate = corr)
```

We are Fisher's z transforming the correlation coefficient and calculating the number of participants and the SE for the specification of the Bf model.

```{r}
processed <-
  processed %>% 
  mutate(corr_bf_effect = map(corr,
                              ~ tibble(estimateFisherZ = DescTools::FisherZ(.x$estimate),
                                       df = .x$parameter,
                                       se = 1 / sqrt(df + 1))))
```

We are joining the original preprocessed effects to the effects of the current study.

```{r}
processed <-
  processed %>% 
  left_join(., original_corr_effect, by = "task") %>% 
  mutate(corr_bf_effect = map2(corr_bf_effect, original_corr_estimate,
                               ~ mutate(.x, originalEstimateFisherZ = DescTools::FisherZ(.y))))
```

#### Running the Bayes factor analysis

A two-tailed normal distribution is used as a prior, with a mode of 0 and a SD of 0.549 (corresponding to an r of 0.5).

```{r}
processed <-
  processed %>% 
  mutate(bf_corr = map(corr_bf_effect,
                     ~ Bf(sd = .x$se,
                          obtained = .x$estimateFisherZ,
                          dfdata = .x$df,
                          meanoftheory = 0,
                          sdtheory = .x$originalEstimateFisherZ,
                          dftheory = 10^10,
                          tail = 1)))

# Print the results with names
processed$bf_corr %>% 
  set_names(., processed$task)
```

###### Finding the robustness region for the correlation analysis

```{r}
# Finding the lower boundary of the robustness region
## Primeprobe
primeprobe_corr <- processed$corr_bf_effect[[1]]

### The lower limit is 0, because the Bf is inconslusive

## Flanker
flanker_corr <- processed$corr_bf_effect[[2]]

### The lower limit is 0, because the Bf is inconslusive

## Stroop
stroop_corr <- processed$corr_bf_effect[[3]]

### The lower limit is 0, because the Bf is inconslusive

## Simon
simon_corr <- processed$corr_bf_effect[[4]]

Bf(sd = simon_corr$se,
   obtained = simon_corr$estimateFisherZ,
   dfdata = simon_corr$df,
   meanoftheory = 0,
   sdtheory = 0.054,
   dftheory = 10^10,
   tail = 1)
```

```{r}
# Finding the upper boundary of the robustness region
## Primeprobe
Bf(sd = primeprobe_corr$se,
   obtained = primeprobe_corr$estimateFisherZ,
   dfdata = primeprobe_corr$df,
   meanoftheory = 0,
   sdtheory = 1.3,
   dftheory = 10^10,
   tail = 1)

## Flanker
Bf(sd = flanker_corr$se,
   obtained = flanker_corr$estimateFisherZ,
   dfdata = flanker_corr$df,
   meanoftheory = 0,
   sdtheory = 0.89,
   dftheory = 10^10,
   tail = 1)

## Stroop
Bf(sd = stroop_corr$se,
   obtained = stroop_corr$estimateFisherZ,
   dfdata = stroop_corr$df,
   meanoftheory = 0,
   sdtheory = 0.45,
   dftheory = 10^10,
   tail = 1)

## Simon
Bf(sd = simon_corr$se,
   obtained = simon_corr$estimateFisherZ,
   dfdata = simon_corr$df,
   meanoftheory = 0,
   sdtheory = 0.99,
   dftheory = 10^10,
   tail = 1)
```

# Accuracy analysis
## I. analysis: Testing whether the CSE was present in the different tasks
### Preprocessing data

We are calculating the accuracy for the analysis for each participant in each condition for each task.

```{r calculating mean dev acc}
processed <-
  processed %>% 
  mutate(acc_data = map(acc_data,
                        . %>% 
                          group_by(participant_id, isCongruent, isPrevCongruent) %>% 
                          summarise(AccConditionalMean = mean(isCorrect, na.rm = T)) %>% 
                          ungroup() %>% 
                          mutate(isCongruent = as_factor(isCongruent),
                                 isPrevCongruent = as_factor(isPrevCongruent),
                                 participant_id = as_factor(participant_id))))
```

We are calculating the main congruency effect of the current (n) trial, the main congruency effect of the previous trial (n-1) and the interaction effect (CSE) between them for each participant.

```{r}
processed <-
  processed %>% 
  mutate(acc_effect_data = map(acc_data,
                              . %>% 
                                ungroup() %>% 
                                mutate(condition = case_when(isPrevCongruent == 0L & isCongruent == 0L ~ "ii",
                                                             isPrevCongruent == 0L & isCongruent == 1L ~ "ic",
                                                             isPrevCongruent == 1L & isCongruent == 0L ~ "ci",
                                                             isPrevCongruent == 1L & isCongruent == 1L ~ "cc",
                                                             TRUE ~ NA_character_)) %>% 
                                select(-isPrevCongruent, -isCongruent) %>% 
                                spread(key = condition, value = AccConditionalMean) %>% 
                                mutate(congruencyEffect = ((ci + ii) / 2) - ((cc + ic) / 2),
                                       prevCongruencyEffect = ((ci + cc) / 2) - ((ii + ic) / 2),
                                       cseEffect = (ci - cc) - (ii - ic))))
```

### Running the analysis

We are running a 2 × 2 repeated-measures ANOVA per task with accuracy as dependent variable. The two factors of the ANOVA are Previous Trial Congruency (congruent, incongruent) and Current Trial Congruency (congruent, incongruent).

```{r}
processed <-
  processed %>% 
  mutate(anova_acc = map(acc_data,
                        ~ aov(AccConditionalMean ~ isCongruent * isPrevCongruent + Error(participant_id / (isCongruent * isPrevCongruent)), data = .x)))

# Print the results with names
processed$anova_acc %>% 
  set_names(., processed$task) %>% 
  map(., summary)
```

None of the interaction effects are significant, therefore we will not run any follow up analysis for the accuracy analysis.

We are calculating the partial eta square as an effect size estimate for the ANOVA.

```{r}
processed <-
  processed %>% 
  mutate(anova_acc_eta = map(anova_acc,
                            ~ eta_sq(.x, partial = TRUE) %>% 
                              filter(stratum == "participant_id:isCongruent:isPrevCongruent")))

# Print results with names
processed$anova_acc_eta %>% 
  set_names(., processed$task) %>% 
  knitr::kable()
```

### Calculating Bayes factors
#### Preprocessing data for the Bayes factor analysis

We are substracting the F values for each effect (main congruency effect of current trial, main congruency effect of previous trial and interaction effect) from the ANOVA.

```{r}
processed <-
  processed %>% 
  mutate(anova_acc_f = map(anova_acc,
                          . %>%
                            broom::tidy() %>% 
                            select(term, statistic) %>% 
                            transmute(term = case_when(term == "isCongruent" ~ "mainEffectIsCongruent",
                                                       term == "isPrevCongruent" ~ "mainEffectIsPrevCongruent",
                                                       term == "isCongruent:isPrevCongruent" ~ "interactionEffect",
                                                       term == "Residuals" ~ term),
                                      fValue = statistic)))
```

We are calculating the mean effects (main congruency effect of current trial, main congruency effect of previous trial and interaction effect) for each task.

```{r}
processed <-
  processed %>% 
  mutate(anova_acc_raw_effect = map(acc_effect_data,
                                   . %>% 
                                     summarise(mainEffectIsCongruent = mean(congruencyEffect, na.rm = T),
                                               mainEffectIsPrevCongruent = mean(prevCongruencyEffect, na.rm = T),
                                               interactionEffect = mean(cseEffect, na.rm = T),
                                               nParticipant = n()) %>%
                                     gather(key = "term", value = "rawEffect", -nParticipant)))
```

We are adding the calucalted F-values for each effect to the calculated mean effects.

```{r}
processed <-
  processed %>% 
  mutate(anova_acc_raw_effect = map2(anova_acc_raw_effect, anova_acc_f,
                                    ~ inner_join(.x, .y, by = "term") %>% 
                                      mutate(se = rawEffect / sqrt(fValue))))
```

Preprocessing effect sizes that Weissman et al. (2014) found for the H1 model of the Bayes factor analysis.

```{r}
original_acc_effect <-
  original %>% 
  filter(measure == "acc") %>% 
  select(task, congruency, cse) %>% 
  gather(key = "term", value = "originalRawEffect", -task) %>%
  mutate(term = case_when(term == "congruency" ~ "mainEffectIsCongruent",
                          term == "cse" ~ "interactionEffect")) %>% 
  group_by(task) %>% 
  nest() %>% 
  rename(original_effect_accuracy = data)
```

Joining the effect sizes of the replicated study with the effect sizes of the current study.

```{r}
processed <-
  processed %>%
  mutate(task = as.character(task)) %>% # vector needs to be character class (instead of a list) for joining
  left_join(., original_acc_effect, by = "task") %>% 
  mutate(anova_acc_bf_effect = map2(anova_acc_raw_effect, original_effect_accuracy,
                              ~ left_join(.x, .y, by = "term")))
```

#### Calculating the Bayes factor
##### For the interaction effect

A half-normal distribution is used with a mode of 0 and a SD equal to the half of the congruency effect (in ms) for the particular task in the original data set by Weissman et al. (2014) to model the prediction of the interaction effect. 

```{r}
processed <-
  processed %>% 
  mutate(bf_interaction_acc = map(anova_acc_bf_effect,
                                 ~ Bf(sd = .x$se[3],
                                      obtained = -(.x$rawEffect[3]),
                                      dfdata = .x$nParticipant[1] - 1,
                                      meanoftheory = 0,
                                      sdtheory = .x$originalRawEffect[1] / 2,
                                      dftheory = 10^10,
                                      tail = 1)))

# Print results with names
processed$bf_interaction_acc %>% 
  set_names(., processed$task)
```

###### Finding the robustness region for the interaction effect

```{r}
# Finding the lower boundary of the robustness region
## Primeprobe
primeprobe_acc <- processed$anova_acc_bf_effect[[1]]

### The lower limit is 0, because the Bf is inconslusive

## Flanker
flanker_acc <- processed$anova_acc_bf_effect[[2]]

### The lower limit is 0, because the Bf is inconslusive

## Stroop
stroop_acc <- processed$anova_acc_bf_effect[[3]]

### The lower limit is 0, because the Bf is inconslusive

## Simon
simon_acc <- processed$anova_acc_bf_effect[[4]]

Bf(sd = simon_acc$se[3],
   obtained = -(simon_acc$rawEffect[3]),
   dfdata = simon_acc$nParticipant[1] - 1,
   meanoftheory = 0,
   sdtheory = 0.0021, 
   dftheory = 10^10,
   tail = 1)

# Finding the upper boundary of the robustness region
## Primeprobe

Bf(sd = primeprobe_acc$se[3],
   obtained = -(primeprobe_acc$rawEffect[3]),
   dfdata = primeprobe_acc$nParticipant[1] - 1,
   meanoftheory = 0,
   sdtheory = 0.0113,
   dftheory = 10^10,
   tail = 1)

## Flanker
Bf(sd = flanker_acc$se[3],
   obtained = -(flanker_acc$rawEffect[3]),
   dfdata = flanker_acc$nParticipant[1] - 1,
   meanoftheory = 0,
   sdtheory = 0.0062,
   dftheory = 10^10,
   tail = 1)

## Stroop
Bf(sd = stroop_acc$se[3],
   obtained = -(stroop_acc$rawEffect[3]),
   dfdata = stroop_acc$nParticipant[1] - 1,
   meanoftheory = 0,
   sdtheory = 0.063,
   dftheory = 10^10,
   tail = 1)

## Simon
Bf(sd = simon_acc$se[3],
   obtained = -(simon_acc$rawEffect[3]),
   dfdata = simon_acc$nParticipant[1] - 1,
   meanoftheory = 0,
   sdtheory = 0.085,
   dftheory = 10^10,
   tail = 1)
```

##### For the main congruency effect of the current trial

To test the congruency main effect, the SD of the H1 model is set to the congruency effect reported by Weissman et al. (2014) for the given task. The other settings stay the same.

```{r}
processed <-
  processed %>% 
  mutate(bf_congruency_acc = map(anova_acc_bf_effect,
                                ~ Bf(sd = .x$se[1],
                                     obtained = -(.x$rawEffect[1]),
                                     dfdata = .x$nParticipant[1] - 1,
                                     meanoftheory = 0,
                                     sdtheory = .x$originalRawEffect[1],
                                     dftheory = 10^10,
                                     tail = 1))) 

# Print results with names
processed$bf_congruency_acc %>% 
  set_names(., processed$task)
```

###### Finding the robustness region for the main congruency effect

```{r}
# Finding the lower boundary of the robustness region
## Primeprobe
Bf(sd = primeprobe_acc$se[1],
   obtained = -(primeprobe_acc$rawEffect[1]),
   dfdata = primeprobe_acc$nParticipant[1] - 1,
   meanoftheory = 0,
   sdtheory = 0.00049,
   dftheory = 10^10,
   tail = 1)

## Flanker

### The lower limit is 0, because the Bf is inconslusive

## Stroop
Bf(sd = stroop_acc$se[1],
   obtained = -(stroop_acc$rawEffect[1]),
   dfdata = stroop_acc$nParticipant[1] - 1,
   meanoftheory = 0,
   sdtheory = 0.00049,
   dftheory = 10^10,
   tail = 1)

## Simon
Bf(sd = simon_acc$se[1],
   obtained = -(simon_acc$rawEffect[1]),
   dfdata = simon_acc$nParticipant[1] - 1,
   meanoftheory = 0,
   sdtheory = 0.00082,
   dftheory = 10^10,
   tail = 1)

# Finding the upper boundary of the robustness region
## Primeprobe
Bf(sd = primeprobe_acc$se[1],
   obtained = -(primeprobe_acc$rawEffect[1]),
   dfdata = primeprobe_acc$nParticipant[1] - 1,
   meanoftheory = 0,
   sdtheory = 3.32,
   dftheory = 10^10,
   tail = 1)

## Flanker
Bf(sd = flanker_acc$se[1],
   obtained = -(flanker_acc$rawEffect[1]),
   dfdata = flanker_acc$nParticipant[1] - 1,
   meanoftheory = 0,
   sdtheory = 0.0152,
   dftheory = 10^10,
   tail = 1)

## Stroop
Bf(sd = stroop_acc$se[1],
   obtained = -(stroop_acc$rawEffect[1]),
   dfdata = stroop_acc$nParticipant[1] - 1,
   meanoftheory = 0,
   sdtheory = 4.449,
   dftheory = 10^10,
   tail = 1)

## Simon
Bf(sd = simon_acc$se[1],
   obtained = -(simon_acc$rawEffect[1]),
   dfdata = simon_acc$nParticipant[1] - 1,
   meanoftheory = 0,
   sdtheory = 10.04,
   dftheory = 10^10,
   tail = 1)
```

### Follow up simple effects analysis
#### Preprocessing data

Saving the data of the task (primeprobe) with significant interaction between current trial congruency and previous trial congruency.

```{r}
simon <- processed$acc_effect_data[[4]]
```

#### Running the analysis

If the Bayes factor is greater than the preset threshold of B = 3 for the interaction effect in a given task we run a simple effects analyses, contrasting current-congruent congruent (cC) and current-incongruent congruent (iC) trials.

```{r current congruent follow up comparison for acc}
current_congruent <- 
  t.test(simon$cc, simon$ic, paired = TRUE) %>% 
  broom::tidy()

current_congruent %>% knitr::kable()
```

current-congruent incongruent (cI) and current-incongruent incongruent (iI) trials in another.

```{r current incongruent follow up comparison for acc}
current_incongruent <- 
  t.test(simon$ci, simon$ii, paired = TRUE) %>% 
  broom::tidy()

current_incongruent %>% knitr::kable()
```

#### Calculating the Bayes factor for the follow up analysis

The prior H1 model is a half-normal distribution with a mode of 0 and an SD equal to half of the CSE estimate for that given task, i.e., the originally reported congruency effect divided by four.

#### Current congruent follow-up comparison

```{r}
Bf(sd = current_congruent$estimate / current_congruent$statistic,
   obtained = current_congruent$estimate,
   dfdata = current_congruent$parameter,
   meanoftheory = 0,
   sdtheory = processed$anova_acc_bf_effect[[4]]$originalRawEffect[1] / 4,
   dftheory = 10^10,
   tail = 1)
```

##### Finding the robustness region for the current congruent follow up analysis

```{r}
# Finding the lower boundary of the robustness region
Bf(sd = current_congruent$estimate / current_congruent$statistic,
   obtained = current_congruent$estimate,
   dfdata = current_congruent$parameter,
   meanoftheory = 0,
   sdtheory = 0.00145,
   dftheory = 10^10,
   tail = 1)

# Finding the upper boundary of the robustness region
Bf(sd = current_congruent$estimate / current_congruent$statistic,
   obtained = current_congruent$estimate,
   dfdata = current_congruent$parameter,
   meanoftheory = 0,
   sdtheory = 0.026,
   dftheory = 10^10,
   tail = 1)
```

We found moderate evidence for the H1 (t = 2.39, bf = 6.12, RR[0.00083, 0.015]).

#### Current icongruent follow-up comparison

```{r}
Bf(sd = current_incongruent$estimate / current_incongruent$statistic,
   obtained = -(current_incongruent$estimate),
   dfdata = current_incongruent$parameter,
   meanoftheory = 0,
   sdtheory = processed$anova_acc_bf_effect[[4]]$originalRawEffect[1] / 4,
   dftheory = 10^10,
   tail = 1)
```

##### Finding the robustness region for the current incongruent follow up analysis

```{r}
# Finding the lower boundary of the robustness region

### The lower limit is 0, because the Bf is inconslusive

# Finding the upper boundary of the robustness region
Bf(sd = current_incongruent$estimate / current_incongruent$statistic,
   obtained = -(current_incongruent$estimate),
   dfdata = current_incongruent$parameter,
   meanoftheory = 0,
   sdtheory = 0.062,
   dftheory = 10^10,
   tail = 1)
```

# Exploratory analyses
## I. analysis: Testing whether the CSE was present in the different tasks based on error rates

One of the reviewers requested to compute the CSE based on the error rates after excluding the post error trials, therefore we account for the possible confounding effect of post-error processing.

### Preprocessing data

We are calculating the accuracy for the analysis for each participant in each condition for each task.

```{r calculating mean dev acc exp}
processed <-
  processed %>% 
  mutate(acc_exp_data = map(acc_exp_data,
                        . %>% 
                          group_by(participant_id, isCongruent, isPrevCongruent) %>% 
                          summarise(AccConditionalMean = mean(isCorrect, na.rm = T)) %>% 
                          ungroup() %>% 
                          mutate(isCongruent = as_factor(isCongruent),
                                 isPrevCongruent = as_factor(isPrevCongruent),
                                 participant_id = as_factor(participant_id))))
```

We are calculating the main congruency effect of the current (n) trial, the main congruency effect of the previous trial (n-1) and the interaction effect (CSE) between them for each participant.

```{r}
processed <-
  processed %>% 
  mutate(acc_exp_effect_data = map(acc_exp_data,
                              . %>% 
                                ungroup() %>% 
                                mutate(condition = case_when(isPrevCongruent == 0L & isCongruent == 0L ~ "ii",
                                                             isPrevCongruent == 0L & isCongruent == 1L ~ "ic",
                                                             isPrevCongruent == 1L & isCongruent == 0L ~ "ci",
                                                             isPrevCongruent == 1L & isCongruent == 1L ~ "cc",
                                                             TRUE ~ NA_character_)) %>% 
                                select(-isPrevCongruent, -isCongruent) %>% 
                                spread(key = condition, value = AccConditionalMean) %>% 
                                mutate(congruencyEffect = ((ci + ii) / 2) - ((cc + ic) / 2),
                                       prevCongruencyEffect = ((ci + cc) / 2) - ((ii + ic) / 2),
                                       cseEffect = (ci - cc) - (ii - ic))))
```

### Running the analysis

We are running a 2 × 2 repeated-measures ANOVA per task with accuracy as dependent variable. The two factors of the ANOVA are Previous Trial Congruency (congruent, incongruent) and Current Trial Congruency (congruent, incongruent).

```{r}
processed <-
  processed %>% 
  mutate(anova_acc_exp = map(acc_exp_data,
                        ~ aov(AccConditionalMean ~ isCongruent * isPrevCongruent + Error(participant_id / (isCongruent * isPrevCongruent)), data = .x)))

# Print the results with names
processed$anova_acc_exp %>% 
  set_names(., processed$task) %>% 
  map(., summary)
```

None of the interaction effects are significant, therefore we will not run any follow up analysis for the accuracy analysis.

We are calculating the partial eta square as an effect size estimate for the ANOVA.

```{r}
processed <-
  processed %>% 
  mutate(anova_acc_exp_eta = map(anova_acc_exp,
                            ~ eta_sq(.x, partial = TRUE) %>% 
                              filter(stratum == "participant_id:isCongruent:isPrevCongruent")))

# Print results with names
processed$anova_acc_exp_eta %>% 
  set_names(., processed$task) %>% 
  knitr::kable()
```

### Calculating Bayes factors
#### Preprocessing data for the Bayes factor analysis

We are substracting the F values for each effect (main congruency effect of current trial, main congruency effect of previous trial and interaction effect) from the ANOVA.

```{r}
processed <-
  processed %>% 
  mutate(anova_acc_exp_f = map(anova_acc_exp,
                          . %>%
                            broom::tidy() %>% 
                            select(term, statistic) %>% 
                            transmute(term = case_when(term == "isCongruent" ~ "mainEffectIsCongruent",
                                                       term == "isPrevCongruent" ~ "mainEffectIsPrevCongruent",
                                                       term == "isCongruent:isPrevCongruent" ~ "interactionEffect",
                                                       term == "Residuals" ~ term),
                                      fValue = statistic)))
```

We are calculating the mean effects (main congruency effect of current trial, main congruency effect of previous trial and interaction effect) for each task.

```{r}
processed <-
  processed %>% 
  mutate(anova_acc_exp_raw_effect = map(acc_exp_effect_data,
                                   . %>% 
                                     summarise(mainEffectIsCongruent = mean(congruencyEffect, na.rm = T),
                                               mainEffectIsPrevCongruent = mean(prevCongruencyEffect, na.rm = T),
                                               interactionEffect = mean(cseEffect, na.rm = T),
                                               nParticipant = n()) %>%
                                     gather(key = "term", value = "rawEffect", -nParticipant)))
```

We are adding the calucalted F-values for each effect to the calculated mean effects.

```{r}
processed <-
  processed %>% 
  mutate(anova_acc_exp_raw_effect = map2(anova_acc_exp_raw_effect, anova_acc_exp_f,
                                    ~ inner_join(.x, .y, by = "term") %>% 
                                      mutate(se = rawEffect / sqrt(fValue))))
```

We are using the same preprocessed original effect sizes as we used for the conformatory accuracy analysis saved in the object named: original_acc_effect.

Joining the effect sizes of the replicated study with the effect sizes of the current study.

```{r}
processed <-
  processed %>%
  mutate(task = as.character(task)) %>% # vector needs to be character class (instead of a list) for joining
  mutate(anova_acc_exp_bf_effect = map2(anova_acc_exp_raw_effect, original_effect_accuracy,
                              ~ left_join(.x, .y, by = "term")))
```

#### Calculating the Bayes factor
##### For the interaction effect

A half-normal distribution is used with a mode of 0 and a SD equal to the half of the congruency effect (in ms) for the particular task in the original data set by Weissman et al. (2014) to model the prediction of the interaction effect. 

```{r}
processed <-
  processed %>% 
  mutate(bf_interaction_acc_exp = map(anova_acc_exp_bf_effect,
                                 ~ Bf(sd = .x$se[3],
                                      obtained = -(.x$rawEffect[3]),
                                      dfdata = .x$nParticipant[1] - 1,
                                      meanoftheory = 0,
                                      sdtheory = .x$originalRawEffect[1] / 2,
                                      dftheory = 10^10,
                                      tail = 1)))

# Print results with names
processed$bf_interaction_acc_exp %>% 
  set_names(., processed$task)
```

###### Finding the robustness region for the interaction effect

```{r}
# Finding the lower boundary of the robustness region
## Primeprobe
primeprobe_acc_exp <- processed$anova_acc_exp_bf_effect[[1]]

### The lower limit is 0, because the Bf is inconslusive

## Flanker
flanker_acc_exp <- processed$anova_acc_exp_bf_effect[[2]]

### The lower limit is 0, because the Bf is inconslusive

## Stroop
stroop_acc_exp <- processed$anova_acc_exp_bf_effect[[3]]

### The lower limit is 0, because the Bf is inconslusive

## Simon
simon_acc_exp <- processed$anova_acc_exp_bf_effect[[4]]

# Bf(sd = simon_acc_exp$se[3],
#    obtained = -(simon_acc_exp$rawEffect[3]),
#    dfdata = simon_acc_exp$nParticipant[1] - 1,
#    meanoftheory = 0,
#    sdtheory = , 
#    dftheory = 10^10,
#    tail = 1)

# Finding the upper boundary of the robustness region
## Primeprobe
# Bf(sd = primeprobe_acc_exp$se[3],
#    obtained = -(primeprobe_acc_exp$rawEffect[3]),
#    dfdata = primeprobe_acc_exp$nParticipant[1] - 1,
#    meanoftheory = 0,
#    sdtheory = ,
#    dftheory = 10^10,
#    tail = 1)

## Flanker
# Bf(sd = flanker_acc_exp$se[3],
#    obtained = -(flanker_acc_exp$rawEffect[3]),
#    dfdata = flanker_acc_exp$nParticipant[1] - 1,
#    meanoftheory = 0,
#    sdtheory = ,
#    dftheory = 10^10,
#    tail = 1)

## Stroop
# Bf(sd = stroop_acc_exp$se[3],
#    obtained = -(stroop_acc_exp$rawEffect[3]),
#    dfdata = stroop_acc_exp$nParticipant[1] - 1,
#    meanoftheory = 0,
#    sdtheory = ,
#    dftheory = 10^10,
#    tail = 1)

## Simon
# Bf(sd = simon_acc_exp$se[3],
#    obtained = -(simon_acc_exp$rawEffect[3]),
#    dfdata = simon_acc_exp$nParticipant[1] - 1,
#    meanoftheory = 0,
#    sdtheory = ,
#    dftheory = 10^10,
#    tail = 1)
```

##### For the main congruency effect of the current trial

To test the congruency main effect, the SD of the H1 model is set to the congruency effect reported by Weissman et al. (2014) for the given task. The other settings stay the same.

```{r}
processed <-
  processed %>% 
  mutate(bf_congruency_acc_exp = map(anova_acc_exp_bf_effect,
                                ~ Bf(sd = .x$se[1],
                                     obtained = -(.x$rawEffect[1]),
                                     dfdata = .x$nParticipant[1] - 1,
                                     meanoftheory = 0,
                                     sdtheory = .x$originalRawEffect[1],
                                     dftheory = 10^10,
                                     tail = 1))) 

# Print results with names
processed$bf_congruency_acc_exp %>% 
  set_names(., processed$task)
```

###### Finding the robustness region for the main congruency effect

```{r}
# Finding the lower boundary of the robustness region
## Primeprobe
# Bf(sd = primeprobe_acc_exp$se[1],
#    obtained = -(primeprobe_acc_exp$rawEffect[1]),
#    dfdata = primeprobe_acc_exp$nParticipant[1] - 1,
#    meanoftheory = 0,
#    sdtheory = ,
#    dftheory = 10^10,
#    tail = 1)

## Flanker

### The lower limit is 0, because the Bf is inconslusive

## Stroop
# Bf(sd = stroop_acc_exp$se[1],
#    obtained = -(stroop_acc_exp$rawEffect[1]),
#    dfdata = stroop_acc_exp$nParticipant[1] - 1,
#    meanoftheory = 0,
#    sdtheory = ,
#    dftheory = 10^10,
#    tail = 1)

## Simon
# Bf(sd = simon_acc_exp$se[1],
#    obtained = -(simon_acc_exp$rawEffect[1]),
#    dfdata = simon_acc_exp$nParticipant[1] - 1,
#    meanoftheory = 0,
#    sdtheory = ,
#    dftheory = 10^10,
#    tail = 1)

# Finding the upper boundary of the robustness region
## Primeprobe
# Bf(sd = primeprobe_acc_exp$se[1],
#    obtained = -(primeprobe_acc_exp$rawEffect[1]),
#    dfdata = primeprobe_acc_exp$nParticipant[1] - 1,
#    meanoftheory = 0,
#    sdtheory = ,
#    dftheory = 10^10,
#    tail = 1)

## Flanker
# Bf(sd = flanker_acc_exp$se[1],
#    obtained = -(flanker_acc_exp$rawEffect[1]),
#    dfdata = flanker_acc_exp$nParticipant[1] - 1,
#    meanoftheory = 0,
#    sdtheory = ,
#    dftheory = 10^10,
#    tail = 1)

## Stroop
# Bf(sd = stroop_acc_exp$se[1],
#    obtained = -(stroop_acc_exp$rawEffect[1]),
#    dfdata = stroop_acc_exp$nParticipant[1] - 1,
#    meanoftheory = 0,
#    sdtheory = ,
#    dftheory = 10^10,
#    tail = 1)

## Simon
# Bf(sd = simon_acc_exp$se[1],
#    obtained = -(simon_acc_exp$rawEffect[1]),
#    dfdata = simon_acc_exp$nParticipant[1] - 1,
#    meanoftheory = 0,
#    sdtheory = ,
#    dftheory = 10^10,
#    tail = 1)
```

### Follow up simple effects analysis
#### Preprocessing data

Saving the data of the task (primeprobe) with significant interaction between current trial congruency and previous trial congruency.

```{r}
simon <- processed$acc_exp_effect_data[[4]]
```

#### Running the analysis

If the Bayes factor is greater than the preset threshold of B = 3 for the interaction effect in a given task we run a simple effects analyses, contrasting current-congruent congruent (cC) and current-incongruent congruent (iC) trials.

```{r current congruent follow up comparison for acc exp}
current_congruent <- 
  t.test(simon$cc, simon$ic, paired = TRUE) %>% 
  broom::tidy()

current_congruent %>% knitr::kable()
```

current-congruent incongruent (cI) and current-incongruent incongruent (iI) trials in another.

```{r current incongruent follow up comparison for acc exp}
current_incongruent <- 
  t.test(simon$ci, simon$ii, paired = TRUE) %>% 
  broom::tidy()

current_incongruent %>% knitr::kable()
```

#### Calculating the Bayes factor for the follow up analysis

The prior H1 model is a half-normal distribution with a mode of 0 and an SD equal to half of the CSE estimate for that given task, i.e., the originally reported congruency effect divided by four.

#### Current congruent follow-up comparison

```{r}
Bf(sd = current_congruent$estimate / current_congruent$statistic,
   obtained = current_congruent$estimate,
   dfdata = current_congruent$parameter,
   meanoftheory = 0,
   sdtheory = processed$anova_acc_exp_bf_effect[[4]]$originalRawEffect[1] / 4,
   dftheory = 10^10,
   tail = 1)
```

##### Finding the robustness region for the current congruent follow up analysis

```{r}
# Finding the lower boundary of the robustness region
# Bf(sd = current_congruent$estimate / current_congruent$statistic,
#    obtained = current_congruent$estimate,
#    dfdata = current_congruent$parameter,
#    meanoftheory = 0,
#    sdtheory = ,
#    dftheory = 10^10,
#    tail = 1)

# Finding the upper boundary of the robustness region
# Bf(sd = current_congruent$estimate / current_congruent$statistic,
#    obtained = current_congruent$estimate,
#    dfdata = current_congruent$parameter,
#    meanoftheory = 0,
#    sdtheory = ,
#    dftheory = 10^10,
#    tail = 1)
```

We found moderate evidence for the H1 (t = 2.39, bf = 6.12, RR[0.00083, 0.015]).

#### Current icongruent follow-up comparison

```{r}
Bf(sd = current_incongruent$estimate / current_incongruent$statistic,
   obtained = -(current_incongruent$estimate),
   dfdata = current_incongruent$parameter,
   meanoftheory = 0,
   sdtheory = processed$anova_acc_exp_bf_effect[[4]]$originalRawEffect[1] / 4,
   dftheory = 10^10,
   tail = 1)
```

##### Finding the robustness region for the current incongruent follow up analysis

```{r}
# Finding the lower boundary of the robustness region

### The lower limit is 0, because the Bf is inconslusive

# Finding the upper boundary of the robustness region
# Bf(sd = current_incongruent$estimate / current_incongruent$statistic,
#    obtained = -(current_incongruent$estimate),
#    dfdata = current_incongruent$parameter,
#    meanoftheory = 0,
#    sdtheory = ,
#    dftheory = 10^10,
#    tail = 1)
```

## II. analysis: Testing whether the CSE was present with BIS

One of the reviewers requested a test with Balanced Integration Scores (BIS, Liesefeldt & Janczyk, 2019, Behavior Research Methods, 51, 40–60).

### Joining the mean accuracy and mean rt data per condition per participant

```{r}
processed <-
  processed %>% 
  mutate(bis_data = map2(rt_data, acc_data,
                         ~ left_join(.x, .y, by = c("participant_id", "isCongruent", "isPrevCongruent"))))
```

### Calculating BIS score

```{r}
processed <-
  processed %>% 
  mutate(bis_data = map(bis_data,
                        . %>% 
                          mutate(N = length(participant_id),
                                 bis_score = ((AccConditionalMean - mean(AccConditionalMean)) /  sqrt( ((N - 1) / N) * var(AccConditionalMean))) - ((rtConditionalMean - mean(rtConditionalMean)) /  sqrt( ((N - 1) / N) * var(rtConditionalMean))))))
```

### Running analysis

```{r}
processed <-
  processed %>%
  mutate(anova_bis = map(bis_data,
                        ~ aov(bis_score ~ isCongruent * isPrevCongruent + Error(participant_id / (isCongruent * isPrevCongruent)), data = .x)))

# Print the results with names
processed$anova_bis %>%
  set_names(., processed$task) %>%
  map(., summary)
```