---
title: "Analisis paper TOL"
output:
  html_document:
    toc: yes
    df_print: paged
  html_notebook:
    toc: yes
    encoding: "UTF-8"
editor_options:
  chunk_output_type: inline
---

# Tower of London: Problem Parameters and Error Analysis in 5-year-old children from Poor Homes.

```{r setup, include=FALSE} 

knitr::opts_chunk$set(warning = FALSE, message = FALSE) 
```

## 0.  Packages and data loading

```{r}

rm(list = ls())

library(flextable)
library(tidyverse)
library(ggpubr)
library(multcomp)
library(rstatix)
library(rcompanion)

options(scipen = 999)

```



```{r}

## Data loading

data = read.csv("Giovannettietal2024JCCS.csv")  # Main data

errores = read.csv("error_types.csv")  # Error trajectories data


```

## 1. Figure 2. Trial-by-trial attrition in TOL task

```{r, fig.width = 10, fig.height= 5}


data %>% 
  filter(Movimientos < 6)%>%
  ungroup() %>% 
  mutate(n = length(unique(username))) %>% 
  
  group_by(trial_number, Movimientos) %>%
  
  summarise(jugados = (sum(!is.na(trial_result))/n)*100)%>%
  
  ggplot(aes(x=factor(trial_number) , y=jugados)) +
  geom_line(aes(group = 1), size = 0.8, alpha = 0.2)+
  geom_line(aes(group = 1), size = 0.8, alpha = 0.2)+
  geom_point(aes(y=jugados, colour = factor(Movimientos)), size = 4)+
  
  xlab("\nTrial number")+
  ylab("Participants\n atrittion (%)\n")+
  
  scale_color_discrete(name = "Required movements")+
  coord_cartesian(ylim = c(50,100))+
  theme_minimal(base_size = 18)-> p

ggarrange(p, ncol = 1, nrow = 1, common.legend = T, legend = "bottom") -> figure2

figure2

# ggsave("Final_figures/Figure 2.png", figure2, dpi = 300, height = 5, width = 10, bg = "white")





```

# 2. Final sample construction

Final sample is composed by all children who at least got to trial number 14. This final sample is extracted from the `data` object and filtered from reaction times higher or lower than 3 SDs.

```{r}

final_data = data%>%
  filter(username %in% errores$username) %>%
  group_by(username)%>%
  filter(trial_number <= 14)%>%
  mutate(jugados = n()) %>%
  filter(jugados == 14)%>%
  filter(Movimientos >= 3)%>%
  ungroup()%>%
  
  # RT's filtering

  mutate(media_tiempo_mov_poblacional = mean(Movement_time_correct,
                                             na.rm = T),
         sd_tiempo_mov_poblacional = sd(Movement_time_correct,
                                        na.rm = T))%>%
  mutate(media_tiempo_planif_poblacional = mean(Planning_time_correct,
                                                na.rm = T),
         sd_tiempo_planif_poblacional = sd(Planning_time_correct,
                                           na.rm = T))%>%

  filter(is.na(Movement_time_correct)|Movement_time_correct < (media_tiempo_mov_poblacional + sd_tiempo_mov_poblacional*3))%>%
  filter(is.na(Planning_time_correct)|Planning_time_correct < 60)%>%
  mutate(tipo_ensayo = paste(Search.Depth, Goal.Hierarchy, sep = "_"))

final_sample_size = length(unique(final_data$username))



```


Final sample consist of a total of `r final_sample_size`

## 2.1. Table S1. Sample description.

```{r}

UBN = final_data%>%
  filter(trial_number == 14)%>%
  filter(!is.na(NBI))%>%
  rename("UBN" = NBI) %>% 
  summarise(UBN = round((sum(UBN)/n())*100,
            digits = 2))
  

SEX = final_data %>% 
  filter(trial_number == 14)%>%
  group_by(SEXO) %>% 
  summarise(freq = n())  


AGE = final_data %>% 
  filter(trial_number == 14) %>% 
  summarise(Mean = mean(EDAD_INICIO, na.rm = T),
            SD = sd(EDAD_INICIO, na.rm = T)) %>% 
  mutate(across(where(is.numeric), ~round(.x, digits = 2))) %>% 
  summarise(Age = paste0(Mean, " (", SD, ")")) 


data.frame(Variable = c("n","Age (years)", "Sex (female)", "UBN (%)"),
           Value = c(sum(SEX$freq), AGE$Age, SEX$freq[1], UBN$UBN)) %>% 
  flextable() %>% 
  set_caption(
    caption = as_paragraph(
      as_chunk("Table S1. Sample description", props = fp_text_default(bold = TRUE))
    ),
    align_with_table = F
  ) %>% 
  autofit()






```




# 3. Table 2. Descriptive and normality analyses for cognitive measures.

## 3.1 Descriptive analyses
 

```{r}


# Descrptive analyses
 
tabla_descriptivos <- final_data%>%
  mutate(trialresult = trial_result,
         tiempomov = Movement_time_correct,
         mediatiempoplan = Planning_time_correct,
         rtplanerror = Planning_time_error,
         rtmoverror = Movement_time_error)%>%
  dplyr::select(trialresult, mediatiempoplan, tiempomov,rtplanerror,rtmoverror)%>%
  summarise(across(where(is.numeric), list(Mean = ~ mean(.x, na.rm = T),
                                           SD = ~ sd(.x, na.rm = T),
                                           Median = ~ median(.x, na.rm = T),
                                           Q25 = ~ quantile(.x, 0.25, na.rm = T),
                                           Q75 = ~ quantile(.x, 0.75, na.rm = T)),
                  .names = "{.col}.{.fn}"))%>%
  pivot_longer(cols = c(starts_with("trialresult"), starts_with("mediatiempoplan"), starts_with("tiempomov"), starts_with("rtplanerror"), starts_with("rtmoverror")),
               names_to = "Variable")%>%
  extract(Variable, c("Variables", "B"),regex = "([[:alnum:]]+).([[:alnum:]]+)")%>%
  pivot_wider(names_from = B, values_from = value) %>% 
    mutate(Variables = case_when(
    Variables == "mediatiempoplan" ~ "Planning time",
    Variables == "trialresult" ~ "Trial result",
    Variables == "tiempomov" ~ "Movement time",
    Variables == "rtplanerror" ~ "Error planning time",
    Variables == "rtmoverror" ~ "Error movement time"
  ))


```


## 3.2. Normality analyses

```{r}



tabla_normalidad <- as.data.frame(matrix(NA, 5,2))

colnames(tabla_normalidad) <- c("Variables", "Normality")

tabla_normalidad[,1] <- c("Trial result", "Planning time", "Movement time", "Error planning time", "Error movement time")

tabla_normalidad[1, 2] <- NA
tabla_normalidad$Normality[2] <- as.numeric(shapiro.test(final_data$Planning_time_correct)[2])
tabla_normalidad$Normality[3] <- as.numeric(shapiro.test(final_data$Movement_time_correct)[2])
tabla_normalidad$Normality[4] <- as.numeric(shapiro.test(final_data$Planning_time_error)[2])
tabla_normalidad$Normality[5] <- as.numeric(shapiro.test(final_data$Movement_time_error)[2])






```

## 3.3. Table Formatting

```{r}


Table_2 = tabla_descriptivos %>% 
  left_join(tabla_normalidad) %>% 
  mutate('Shapiro Wilk (p)' = ifelse(Normality < 0.001, "< 0.001", signif(Normality, 3)), .keep = "unused")%>%
  
  
  
  flextable()%>%
  colformat_double(digits = 2)%>%
  merge_v(c( "Variables"))  %>%
  set_caption(
    caption = as_paragraph(
      as_chunk("Table 2. Descriptive and normality analyses for cognitive measures.", props = fp_text_default(bold = TRUE))
    ),
    align_with_table = F
  ) %>% autofit()

Table_2



```

# 4. Correct trial analyses

## 4.1. Performance between three- and four-movement trials

### 4.1.1. Table S2. GLM's for performance between and three- and four-movement trials

```{r}

# GLM for trial_result between three- and four-movement trials

binom = glm(trial_result ~ Movimientos, family = "binomial", data = final_data)


table_binom = binom %>%
  tidy() %>% filter(term != "(Intercept)") %>%
  rename("Contrast" = term) %>% 
  mutate(Contrast = c("Movements")) %>% 
     rename("Estimate" = estimate) %>% 
  dplyr::select(-statistic)

table_binom%>% 
  
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", signif(p.value, 3)), .keep = "unused")%>% 
  mutate(across(where(is.numeric),~round(.x, digits = 2))) %>%
  flextable() %>%
  set_caption(
    caption = as_paragraph(
      as_chunk("Table S2. GLM's for performance between and three- and four-movement trials", props = fp_text_default(bold = TRUE))
    ),
    align_with_table = F
  ) %>%
  autofit()

```


### 4.1.2. Table S3. Planning and movement times between three- and four-movement trials


```{r}

Variables = c("Planning time", "Movement time")
p.value <- rep(NA, 2)
effect.size <- rep(NA, 2)

table_times_results <- data.frame(Variables, p.value, effect.size)

# Planning times

    table_times_results$p.value[1] <- as.numeric(wilcox.test(Planning_time_correct ~ Movimientos, 
                                                        data = final_data)[3])
    
    table_times_results$effect.size[1] <- as.numeric(wilcoxonR(final_data$Planning_time_correct,
                                                          final_data$Movimientos, 
                                                          data = final_data))
    
# Movement times
    
    
        table_times_results$p.value[2] <- as.numeric(wilcox.test(Movement_time_correct ~ Movimientos, 
                                                        data = final_data)[3])
    
    table_times_results$effect.size[2] <- as.numeric(wilcoxonR(final_data$Movement_time_correct,
                                                          final_data$Movimientos, 
                                                          data = final_data))
    
    
    

    
table_times_results  %>% 
  
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", signif(p.value, 3)), .keep = "unused")%>% 
  mutate(across(where(is.numeric),~round(.x, digits = 2))) %>%
  flextable() %>%
  set_caption(
    caption = as_paragraph(
      as_chunk("Table S3. Planning and movement times between three- and four-movement trials", props = fp_text_default(bold = TRUE))
    ),
    align_with_table = F
  ) %>%
  autofit()
    
    



```


### 4.1.3. Figure 4. Correct trial analyses between three- and four-movement trials.

```{r, fig.height= 5, fig.width= 10}


final_data %>% 
  group_by(Movimientos) %>%
  summarise(correctos = sum(trial_result, na.rm = T), jugados = sum(!is.na(trial_result)),
            proporcion_correctos = correctos/jugados,
            proporcion_errores = 1- proporcion_correctos,
            sd_correctos = (sd(trial_result)/sqrt(n())))%>%
  ggplot(aes(x=factor(Movimientos) , y=proporcion_correctos)) +
  
  geom_point(size = 4)+
  geom_errorbar(aes(ymin = (proporcion_correctos - sd_correctos), ymax = (proporcion_correctos + sd_correctos)),
                width = 0.2, size = 1)+
  theme_minimal(base_size = 18)+
  ylim(0,1)+
  ylab("Proportion of correct trials\n")+
  xlab("\nRequired movements")+
  theme(plot.margin = margin(1, 0, 0, 0, "cm"),
        strip.text = element_text(vjust = 15, margin = margin(20,1,1,1)))-> p1


final_data %>% 
  group_by(Movimientos) %>%

  ggplot(aes(x=factor(Movimientos) , y= Planning_time_correct)) +
  geom_boxplot()+
  theme_minimal(base_size = 18)+
  ylim(0,50)+
  ylab("Planning time\n")+
  xlab("\nRequired movements")+
  theme(plot.margin = margin(1, 0, 0, 0, "cm"),
        strip.text = element_text(vjust = 15, margin = margin(20,1,1,1))) -> p2


final_data %>% 
  group_by(Movimientos) %>%

  ggplot(aes(x=factor(Movimientos) , y= Movement_time_correct)) +
  geom_boxplot()+
  theme_minimal(base_size = 18)+
  ylim(0,20)+
  ylab("Movement time\n")+
  xlab("\nRequired movements")+
  theme(plot.margin = margin(1, 0, 0, 0, "cm"),
        strip.text = element_text(vjust = 15, margin = margin(20,1,1,1)))-> p3


ggarrange(p1, p2, p3, ncol = 3, nrow = 1, labels =   c("A", "B", "C"), font.label = list(size = 18))-> figure4


figure4


# ggsave("Final_figures/Figure 4.pdf", figure4, dpi = 300, height = 5, width = 10, bg = "white")


```


## 4.2. Performance within three-movement trials

### 4.2.1. Table S5. GLM's for performance within three-movement trials


```{r}

# GLM for trial_result within three-movement trials

final_data %>%                                                                                                  # Data subsetting
  mutate(tipo_ensayo = paste(Search.Depth, Goal.Hierarchy, sep = "_"))%>%
  mutate(tipo_ensayo = factor(tipo_ensayo, levels = c("0_tower", "0_partial","1_partial", "1_plane")))%>%
  filter(Movimientos == 3)->  final_data_3


binom_3 <- glm(trial_result ~ tipo_ensayo, family = "binomial", data = final_data_3)                            # Model implementation
# summary(binom_3)


contrasts_binom_3 = glht(binom_3, linfct = mcp(tipo_ensayo = "Tukey")) %>%                                      # Contrasts extracted from the model
   summary(test = adjusted("holm")) %>% 
   tidy() %>% 
   dplyr::select(contrast, estimate,std.error, adj.p.value) %>% 
   rename("Contrast" = contrast,
          "Estimate" = estimate,
          "p.value" = adj.p.value) 




contrasts_binom_3 %>% 
  
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", signif(p.value, 3)), .keep = "unused")%>% 
  mutate(across(where(is.numeric),~round(.x, digits = 2))) %>%
  flextable() %>%
  set_caption(
    caption = as_paragraph(
      as_chunk("Table S5. GLM's for performance within three-movement trials", props = fp_text_default(bold = TRUE))
    ),
    align_with_table = F
  ) %>%
  autofit()








```


### 4.2.2. Table S6. Planning and movement times within three-movement trials (Kruskal-Wallis)

```{r}

final_data_3_long = final_data_3 %>% 
  dplyr::select(tipo_ensayo, Planning_time_correct, Movement_time_correct) %>% 
  pivot_longer(-c(tipo_ensayo), names_to = "Variables", values_to = "values")


final_data_3_long %>%
  summarise(kruskal_test = tidy(kruskal.test(values ~ tipo_ensayo)), .by = "Variables") %>% 
  unnest(kruskal_test) %>% 
  dplyr::select(Variables, statistic, p.value)%>% 
  
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", round(p.value, 3)), .keep = "unused")%>% 
  mutate(across(where(is.numeric),~round(.x, digits = 2))) %>%
  flextable() %>%
  set_caption(
    caption = as_paragraph(
      as_chunk("Table S6. Planning and movement times within three-movement trials (Kruskal-Wallis)", props = fp_text_default(bold = TRUE))
    ),
    align_with_table = F
  ) %>% 
  autofit()

```

### 4.2.3.  Table S7. Planning and movement times within three-movement trials (contrasts)

```{r}

# P.values


time_contrasts_3 = final_data_3_long %>% 
  reframe(contrasts = tidy(
      
      pairwise.wilcox.test(values,
                         tipo_ensayo,
                         paired = F, 
                         p.adjust.method = "holm")),
      .by = Variables) %>% 
  
  unnest(contrasts)
  

# Effect sizes

time_effsize_3 <- final_data_3_long %>%
  group_by(Variables) %>%
  do(effsizes = wilcox_effsize(values ~ tipo_ensayo, data = .)) %>% 
  unnest(effsizes) %>% 
  rename("group1" = group2,
         "group2" = group1) %>% 
  dplyr::select(Variables, group1, group2, effsize) 


# Complete table



time_contrasts_3 %>%
  left_join(time_effsize_3) %>% 
  mutate(Contrast = paste(group1, "-", group2), .keep = "unused") %>% 
  dplyr::select(Variables, Contrast, everything())%>% 
  
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", round(p.value, 3)), .keep = "unused")%>% 
  mutate(across(where(is.numeric),~round(.x, digits = 2))) %>%
  flextable() %>%
  set_caption(
    caption = as_paragraph(
      as_chunk("Table S7. Planning and movement times within three-movement trials (contrasts)", props = fp_text_default(bold = TRUE))
    ),
    align_with_table = F
  ) %>%
  merge_v(j = "Variables") %>%
  valign(j = "Variables",valign = "top") %>% 
  autofit() %>% 
  fix_border_issues()
  


```

### 4.2.4. Figure 5. Correct trial analyses within three-movement trials

```{r, fig.height= 5, fig.width= 10}

## Plot theming

customize = function(myplot){
  
  myplot +
    theme_minimal(base_size = 15)+
    theme(plot.margin = margin(0, 0.25, 0.75, 0.25, "cm"),
          strip.text = element_text(vjust = 22, margin = margin(70,1,1,1)),
          axis.title.x = element_text(vjust = - 1)) 
  
  
}


## Plot data agreggation

fig5_plot_data = final_data %>% 
  mutate(Goal.Hierarchy = case_when(
  Goal.Hierarchy == "tower" ~ "Tower\n hierarchy",
  Goal.Hierarchy == "partial" ~ "Partial\n hierarchy",
  Goal.Hierarchy == "plane" ~ "Plane\n hierarchy"))%>%
  mutate(Search.Depth = case_when(
  Search.Depth == "0" ~ "Low",
  Search.Depth == "1" ~ "High"))%>%
  mutate(Goal.Hierarchy = factor(Goal.Hierarchy, levels = c("Tower\n hierarchy", "Partial\n hierarchy", "Plane\n hierarchy")))%>%
  mutate(Search.Depth = factor(Search.Depth, levels = c("Low", "High"))) %>% 
    group_by(Search.Depth, Goal.Hierarchy) %>%
  filter(Movimientos == 3) 


## Trial_result within three-movement trials


fig5_plot_data%>%

  summarise(correctos = sum(trial_result, na.rm = T), jugados = sum(!is.na(trial_result)),
            proporcion_correctos = correctos/jugados,
            proporcion_errores = 1- proporcion_correctos,
            sd_correctos = (sd(trial_result)/sqrt(n())))%>%
  ggplot(aes(x=factor(Search.Depth) , y=proporcion_correctos)) +
  
  geom_point(size = 4)+
  geom_errorbar(aes(ymin = (proporcion_correctos - sd_correctos), ymax = (proporcion_correctos + sd_correctos)),
                width = 0.3, size = 1)+
  facet_grid(. ~ Goal.Hierarchy) +
  theme_minimal()%>%
  customize()+

  ylim(0,1)+
  ylab("Proportion of correct trials")+
  xlab("Search depth") -> p1


## Planning_time_correct within three-movement trials


fig5_plot_data %>%

  ggplot(aes(x=factor(Search.Depth) , y=Planning_time_correct))+
  geom_boxplot()+

  facet_grid(. ~ Goal.Hierarchy) +
  theme_minimal()%>% 
  customize()+

  ylim(0,50)+
  ylab("Planning time")+
  xlab("Search depth")  -> p2


## Movement_time_correct within three-movement trials

fig5_plot_data%>%

  ggplot(aes(x=factor(Search.Depth) , y=Movement_time_correct))+
  geom_boxplot()+

  facet_grid(. ~ Goal.Hierarchy)+
  theme_minimal()%>% 
  customize()+

 
  scale_y_continuous(limits = c(0,10), breaks = seq(0,10, by = 2))+
  ylab("Movement time")+
  xlab("Search depth")  -> p3



ggarrange(p1, p2, p3, ncol = 3, nrow = 1, labels =   c("A", "B", "C"), font.label = list(size = 18), widths = c(1,1.1))-> figure5

figure5
  
# 
# ggsave("Final_figures/Figure 5.pdf", figure5, dpi = 300, height = 5, width = 10, bg = "white")



```


## 4.3. Trial result within four-movement trials

### 4.3.1. Table S8. GLM's for performance within four-movement trials

```{r}

# GLM for trial_result within three-movement trials

final_data %>%                                                                                                  # Data subsetting
  mutate(tipo_ensayo = paste(Search.Depth, Goal.Hierarchy, sep = "_"))%>%
  mutate(tipo_ensayo = factor(tipo_ensayo, levels = c("1_tower", "0_partial","1_partial", "0_plane")))%>%
  filter(Movimientos == 4)->  final_data_4

binom_4 <- glm(trial_result ~ tipo_ensayo, family = "binomial", data = final_data_4)                            # Model implementation
# summary(binom_3)


contrasts_binom_4 = glht(binom_4, linfct = mcp(tipo_ensayo = "Tukey")) %>%                                      # Contrasts extracted from the model
   summary(test = adjusted("holm")) %>% 
   tidy() %>% 
   dplyr::select(contrast, estimate,std.error, adj.p.value) %>% 
   rename("Contrast" = contrast,
          "Estimate" = estimate,
          "p.value" = adj.p.value) 




contrasts_binom_4 %>% 
  
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", signif(p.value, 3)), .keep = "unused")%>% 
  mutate(across(where(is.numeric),~round(.x, digits = 2))) %>%
  flextable() %>%
  set_caption(
    caption = as_paragraph(
      as_chunk("Table S8. GLM's for performance within four-movement trials", props = fp_text_default(bold = TRUE))
    ),
    align_with_table = F
  ) %>%
  autofit()







```

### 4.3.2. Table S9. Planning and movement times within four-movement trials (Kruskal-Wallis)

```{r}

final_data_4_long = final_data_4 %>% 
  dplyr::select(tipo_ensayo, Planning_time_correct, Movement_time_correct) %>% 
  pivot_longer(-c(tipo_ensayo), names_to = "Variables", values_to = "values")


final_data_4_long %>%
  summarise(kruskal_test = tidy(kruskal.test(values ~ tipo_ensayo)), .by = "Variables") %>% 
  unnest(kruskal_test) %>% 
  dplyr::select(Variables, statistic, p.value)%>% 
  
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", round(p.value, 3)), .keep = "unused")%>% 
  mutate(across(where(is.numeric),~round(.x, digits = 2))) %>%
  flextable() %>%
  set_caption(
    caption = as_paragraph(
      as_chunk("Table S9. Planning and movement times within four-movement trials (Kruskal-Wallis)", props = fp_text_default(bold = TRUE))
    ),
    align_with_table = F
  ) %>% 
  autofit()

```


### 4.3.3. Table S10. Planning and movement times within four-movement trials (contrasts)

```{r}

# P.values


time_contrasts_4 = final_data_4_long %>% 
  reframe(contrasts = tidy(
      
      pairwise.wilcox.test(values,
                         tipo_ensayo,
                         paired = F, 
                         p.adjust.method = "holm")),
      .by = Variables) %>% 
  
  unnest(contrasts)
  

# Effect sizes

time_effsize_4 <- final_data_4_long %>%
  group_by(Variables) %>%
  do(effsizes = wilcox_effsize(values ~ tipo_ensayo, data = .)) %>% 
  unnest(effsizes) %>% 
  rename("group1" = group2,
         "group2" = group1) %>% 
  dplyr::select(Variables, group1, group2, effsize) 


# Complete table



time_contrasts_4 %>%
  left_join(time_effsize_4) %>% 
  mutate(Contrast = paste(group1, "-", group2), .keep = "unused") %>% 
  dplyr::select(Variables, Contrast, everything())%>% 
  
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", round(p.value, 3)), .keep = "unused")%>% 
  mutate(across(where(is.numeric),~round(.x, digits = 2))) %>%
  flextable() %>%
  set_caption(
    caption = as_paragraph(
      as_chunk("Table S10. Planning and movement times within four-movement trials (contrasts)", props = fp_text_default(bold = TRUE))
    ),
    align_with_table = F
  ) %>%
  merge_v(j = "Variables") %>%
  valign(j = "Variables",valign = "top") %>% 
  autofit() %>% 
  fix_border_issues()
  


```


### 4.3.4. Figure 6. Correct trial analyses within four-movement trials

```{r, fig.height= 5, fig.width= 10}

## Plot theming

customize = function(myplot){
  
  myplot +
    theme_minimal(base_size = 15)+
    theme(plot.margin = margin(0, 0.25, 0.75, 0.25, "cm"),
          strip.text = element_text(vjust = 22, margin = margin(70,1,1,1)),
          axis.title.x = element_text(vjust = - 1)) 
  
  
}


## Plot data agreggation

fig5_plot_data = final_data %>% 
  mutate(Goal.Hierarchy = case_when(
  Goal.Hierarchy == "tower" ~ "Tower\n hierarchy",
  Goal.Hierarchy == "partial" ~ "Partial\n hierarchy",
  Goal.Hierarchy == "plane" ~ "Plane\n hierarchy"))%>%
  mutate(Search.Depth = case_when(
  Search.Depth == "0" ~ "Low",
  Search.Depth == "1" ~ "High"))%>%
  mutate(Goal.Hierarchy = factor(Goal.Hierarchy, levels = c("Tower\n hierarchy", "Partial\n hierarchy", "Plane\n hierarchy")))%>%
  mutate(Search.Depth = factor(Search.Depth, levels = c("Low", "High"))) %>% 
    group_by(Search.Depth, Goal.Hierarchy) %>%
  filter(Movimientos == 4) 


## Trial_result within three-movement trials


fig5_plot_data%>%

  summarise(correctos = sum(trial_result, na.rm = T), jugados = sum(!is.na(trial_result)),
            proporcion_correctos = correctos/jugados,
            proporcion_errores = 1- proporcion_correctos,
            sd_correctos = (sd(trial_result)/sqrt(n())))%>%
  ggplot(aes(x=factor(Search.Depth) , y=proporcion_correctos)) +
  
  geom_point(size = 4)+
  geom_errorbar(aes(ymin = (proporcion_correctos - sd_correctos), ymax = (proporcion_correctos + sd_correctos)),
                width = 0.3, size = 1)+
  facet_grid(. ~ Goal.Hierarchy) +
  theme_minimal()%>%
  customize()+

  ylim(0,1)+
  ylab("Proportion of correct trials")+
  xlab("Search depth") -> p1


## Planning_time_correct within three-movement trials


fig5_plot_data %>%

  ggplot(aes(x=factor(Search.Depth) , y=Planning_time_correct))+
  geom_boxplot()+

  facet_grid(. ~ Goal.Hierarchy) +
  theme_minimal()%>% 
  customize()+

  ylim(0,50)+
  ylab("Planning time")+
  xlab("Search depth")  -> p2


## Movement_time_correct within three-movement trials

fig5_plot_data%>%

  ggplot(aes(x=factor(Search.Depth) , y=Movement_time_correct))+
  geom_boxplot()+

  facet_grid(. ~ Goal.Hierarchy)+
  theme_minimal()%>% 
  customize()+

 
  scale_y_continuous(limits = c(0,10), breaks = seq(0,10, by = 2))+
  ylab("Movement time")+
  xlab("Search depth")  -> p3



ggarrange(p1, p2, p3, ncol = 3, nrow = 1, labels =   c("A", "B", "C"), font.label = list(size = 18), widths = c(1,1.1))-> figure6

figure6
  

# ggsave("Final_figures/Figure 6.pdf", figure6, dpi = 300, height = 5, width = 10, bg = "white")





```

# 5. Error analyses

## 5.1. Table S11. Error response times between three- and four-movement trials


```{r}

# I subset the data to retain trials with equal or +2 movements

final_data%>%
  filter(total_movements <= Movimientos+2)-> final_data_equal_movements



Variables = c("Planning time", "Movement time")
p.value <- rep(NA, 2)
effect.size <- rep(NA, 2)

table_times_results <- data.frame(Variables, p.value, effect.size)

# Planning times

    table_times_results$p.value[1] <- as.numeric(wilcox.test(Planning_time_error ~ Movimientos, 
                                                        data = final_data)[3])
    
    table_times_results$effect.size[1] <- as.numeric(wilcoxonR(final_data$Planning_time_error,
                                                          final_data$Movimientos, 
                                                          data = final_data))
    
# Movement times
    
    
        table_times_results$p.value[2] <- as.numeric(wilcox.test(Movement_time_error ~ Movimientos, 
                                                        data = final_data)[3])
    
    table_times_results$effect.size[2] <- as.numeric(wilcoxonR(final_data$Movement_time_error,
                                                          final_data$Movimientos, 
                                                          data = final_data))
    
    
    

    
table_times_results  %>% 
  
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", signif(p.value, 3)), .keep = "unused")%>% 
  mutate(across(where(is.numeric),~round(.x, digits = 2))) %>%
  flextable() %>%
  set_caption(
    caption = as_paragraph(
      as_chunk("Table S11. Error response times between three- and four-movement trials", props = fp_text_default(bold = TRUE))
    ),
    align_with_table = F
  ) %>%
  autofit()
    
    



```

## 5.2. Error response times within three-movement trials

```{r}

final_data_equal_movements %>%                                                                            # Data subsetting
  mutate(tipo_ensayo = paste(Search.Depth, Goal.Hierarchy, sep = "_"))%>%
  mutate(tipo_ensayo = factor(tipo_ensayo, levels = c("0_tower", "0_partial","1_partial", "1_plane")))%>%
  filter(Movimientos == 3)->  final_data_eq_3



```


### 5.2.1. Table S12. Error response times within three-movement trials (Kruskal-Wallis)

```{r}

final_data_3_long = final_data_eq_3 %>% 
  dplyr::select(tipo_ensayo, Planning_time_error, Movement_time_error) %>% 
  pivot_longer(-c(tipo_ensayo), names_to = "Variables", values_to = "values")


final_data_3_long %>%
  summarise(kruskal_test = tidy(kruskal.test(values ~ tipo_ensayo)), .by = "Variables") %>% 
  unnest(kruskal_test) %>% 
  dplyr::select(Variables, statistic, p.value)%>% 
  
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", round(p.value, 3)), .keep = "unused")%>% 
  mutate(across(where(is.numeric),~round(.x, digits = 2))) %>%
  flextable() %>%
  set_caption(
    caption = as_paragraph(
      as_chunk("Table S12. Error response times within three-movement trials (Kruskal-Wallis)", props = fp_text_default(bold = TRUE))
    ),
    align_with_table = F
  ) %>% 
  autofit()

```

### 5.2.2. Table S13. Error response times within three-movement trials (contrasts)

```{r}

# P.values


time_contrasts_3 = final_data_3_long %>% 
  reframe(contrasts = tidy(
      
      pairwise.wilcox.test(values,
                         tipo_ensayo,
                         paired = F, 
                         p.adjust.method = "holm")),
      .by = Variables) %>% 
  
  unnest(contrasts)
  

# Effect sizes

time_effsize_3 <- final_data_3_long %>%
  group_by(Variables) %>%
  do(effsizes = wilcox_effsize(values ~ tipo_ensayo, data = .)) %>% 
  unnest(effsizes) %>% 
  rename("group1" = group2,
         "group2" = group1) %>% 
  dplyr::select(Variables, group1, group2, effsize) 


# Complete table



time_contrasts_3 %>%
  left_join(time_effsize_3) %>% 
  mutate(Contrast = paste(group1, "-", group2), .keep = "unused") %>% 
  dplyr::select(Variables, Contrast, everything())%>% 
  
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", round(p.value, 3)), .keep = "unused")%>% 
  mutate(across(where(is.numeric),~round(.x, digits = 2))) %>%
  flextable() %>%
  set_caption(
    caption = as_paragraph(
      as_chunk("Table S13. Error response times within three-movement trials (contrasts)", props = fp_text_default(bold = TRUE))
    ),
    align_with_table = F
  ) %>%
  merge_v(j = "Variables") %>%
  valign(j = "Variables",valign = "top") %>% 
  autofit() %>% 
  fix_border_issues()
  


```


### 5.2.3. Figure 7. Error response times within three-movement trials.

```{r, fig.height= 5, fig.width= 10}

## Plot theming

customize = function(myplot){
  
  myplot +
    theme_minimal(base_size = 15)+
    theme(plot.margin = margin(0, 0.25, 0.75, 0.25, "cm"),
          strip.text = element_text(vjust = 22, margin = margin(70,1,1,1)),
          axis.title.x = element_text(vjust = - 1)) 
  
  
}


## Plot data agreggation

fig5_plot_data = final_data_equal_movements %>% 
  mutate(Goal.Hierarchy = case_when(
  Goal.Hierarchy == "tower" ~ "Tower\n hierarchy",
  Goal.Hierarchy == "partial" ~ "Partial\n hierarchy",
  Goal.Hierarchy == "plane" ~ "Plane\n hierarchy"))%>%
  mutate(Search.Depth = case_when(
  Search.Depth == "0" ~ "Low",
  Search.Depth == "1" ~ "High"))%>%
  mutate(Goal.Hierarchy = factor(Goal.Hierarchy, levels = c("Tower\n hierarchy", "Partial\n hierarchy", "Plane\n hierarchy")))%>%
  mutate(Search.Depth = factor(Search.Depth, levels = c("Low", "High"))) %>% 
    group_by(Search.Depth, Goal.Hierarchy) %>%
  filter(Movimientos == 3) 




## Planning_time_correct within three-movement trials


fig5_plot_data %>%

  ggplot(aes(x=factor(Search.Depth) , y=Planning_time_error))+
  geom_boxplot()+

  facet_grid(. ~ Goal.Hierarchy) +
  theme_minimal()%>% 
  customize()+

  ylim(0,50)+
  ylab("Planning error time")+
  xlab("Search depth")  -> p2


## Movement_time_correct within three-movement trials

fig5_plot_data%>%

  ggplot(aes(x=factor(Search.Depth) , y=Movement_time_error))+
  geom_boxplot()+

  facet_grid(. ~ Goal.Hierarchy)+
  theme_minimal()%>% 
  customize()+

 
  scale_y_continuous(limits = c(0,10), breaks = seq(0,10, by = 2))+
  ylab("Movement error time")+
  xlab("Search depth")  -> p3



ggarrange(p2, p3, ncol = 2, nrow = 1, labels =   c("A", "B"), font.label = list(size = 18), widths = c(1,1))-> figure7

figure7
  

# ggsave("Final_figures/Figure 7.pdf", figure7, dpi = 300, height = 5, width = 10, bg = "white")




```


## 5.3. Error response times within four-movement trials

```{r}

final_data_equal_movements %>%                                                                            # Data subsetting
  mutate(tipo_ensayo = paste(Search.Depth, Goal.Hierarchy, sep = "_"))%>%
  mutate(tipo_ensayo = factor(tipo_ensayo, levels = c("1_tower", "0_partial","1_partial", "0_plane")))%>%
  filter(Movimientos == 4)->  final_data_eq_4



```


### 5.3.1. Table S14. Error response times within four-movement trials (Kruskal-Wallis)

```{r}

final_data_4_long = final_data_eq_4 %>% 
  dplyr::select(tipo_ensayo, Planning_time_error, Movement_time_error) %>% 
  pivot_longer(-c(tipo_ensayo), names_to = "Variables", values_to = "values")


final_data_4_long %>%
  summarise(kruskal_test = tidy(kruskal.test(values ~ tipo_ensayo)), .by = "Variables") %>% 
  unnest(kruskal_test) %>% 
  dplyr::select(Variables, statistic, p.value)%>% 
  
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", round(p.value, 3)), .keep = "unused")%>% 
  mutate(across(where(is.numeric),~round(.x, digits = 2))) %>%
  flextable() %>%
  set_caption(
    caption = as_paragraph(
      as_chunk("Table S14. Error response times within three-movement trials (Kruskal-Wallis)", props = fp_text_default(bold = TRUE))
    ),
    align_with_table = F
  ) %>% 
  autofit()

```


### 5.3.2. Table S15. Error response times within four-movement trials (contrasts)

```{r}

# P.values


time_contrasts_4 = final_data_4_long %>% 
  reframe(contrasts = tidy(
      
      pairwise.wilcox.test(values,
                         tipo_ensayo,
                         paired = F, 
                         p.adjust.method = "holm")),
      .by = Variables) %>% 
  
  unnest(contrasts)
  

# Effect sizes

time_effsize_4 <- final_data_4_long %>%
  group_by(Variables) %>%
  do(effsizes = wilcox_effsize(values ~ tipo_ensayo, data = .)) %>% 
  unnest(effsizes) %>% 
  rename("group1" = group2,
         "group2" = group1) %>% 
  dplyr::select(Variables, group1, group2, effsize) 


# Complete table



time_contrasts_4 %>%
  left_join(time_effsize_4) %>% 
  mutate(Contrast = paste(group1, "-", group2), .keep = "unused") %>% 
  dplyr::select(Variables, Contrast, everything())%>% 
  
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", round(p.value, 3)), .keep = "unused")%>% 
  mutate(across(where(is.numeric),~round(.x, digits = 2))) %>%
  flextable() %>%
  set_caption(
    caption = as_paragraph(
      as_chunk("Table S15. Error response times within three-movement trials (contrasts)", props = fp_text_default(bold = TRUE))
    ),
    align_with_table = F
  ) %>%
  merge_v(j = "Variables") %>%
  valign(j = "Variables",valign = "top") %>% 
  autofit() %>% 
  fix_border_issues()
  


```


### 5.3.3. Figure 8. Error response times within four-movement trials.

```{r, fig.height= 5, fig.width= 10}

## Plot theming

customize = function(myplot){
  
  myplot +
    theme_minimal(base_size = 15)+
    theme(plot.margin = margin(0, 0.25, 0.75, 0.25, "cm"),
          strip.text = element_text(vjust = 22, margin = margin(70,1,1,1)),
          axis.title.x = element_text(vjust = - 1)) 
  
  
}


## Plot data agreggation

fig5_plot_data = final_data_equal_movements %>% 
  mutate(Goal.Hierarchy = case_when(
  Goal.Hierarchy == "tower" ~ "Tower\n hierarchy",
  Goal.Hierarchy == "partial" ~ "Partial\n hierarchy",
  Goal.Hierarchy == "plane" ~ "Plane\n hierarchy"))%>%
  mutate(Search.Depth = case_when(
  Search.Depth == "0" ~ "Low",
  Search.Depth == "1" ~ "High"))%>%
  mutate(Goal.Hierarchy = factor(Goal.Hierarchy, levels = c("Tower\n hierarchy", "Partial\n hierarchy", "Plane\n hierarchy")))%>%
  mutate(Search.Depth = factor(Search.Depth, levels = c("Low", "High"))) %>% 
    group_by(Search.Depth, Goal.Hierarchy) %>%
  filter(Movimientos == 4) 




## Planning_time_correct within three-movement trials


fig5_plot_data %>%

  ggplot(aes(x=factor(Search.Depth) , y=Planning_time_error))+
  geom_boxplot()+

  facet_grid(. ~ Goal.Hierarchy) +
  theme_minimal()%>% 
  customize()+

  ylim(0,50)+
  ylab("Planning error time")+
  xlab("Search depth")  -> p2


## Movement_time_correct within three-movement trials

fig5_plot_data%>%

  ggplot(aes(x=factor(Search.Depth) , y=Movement_time_error))+
  geom_boxplot()+

  facet_grid(. ~ Goal.Hierarchy)+
  theme_minimal()%>% 
  customize()+

 
  scale_y_continuous(limits = c(0,10), breaks = seq(0,10, by = 2))+
  ylab("Movement error time")+
  xlab("Search depth")  -> p3



ggarrange(p2, p3, ncol = 2, nrow = 1, labels =   c("A", "B"), font.label = list(size = 18), widths = c(1,1))-> figure8

figure8
  

# ggsave("Final_figures/Figure 8.pdf", figure8, dpi = 300, height = 5, width = 10, bg = "white")




```


## 5.4. Error types

### 5.4.1. Preparing data

```{r}

# Loading error trajectories 

tray_error = left_join(final_data, errores, by = c("username", "trial_number"))


# Subseting and labeling data

tray_error_completa = tray_error%>%
  
  filter(trial_result == 0)%>%                                    # Select only error trials
  filter(!is.na(tipo_error))%>%                                   # Select trials with available trajectories
  filter(Movimientos == 3 | Movimientos == 4)%>%                    # Select only 3 and 4 movement trials
  mutate(Cantidad_mov = case_when(                                  # Cantidad_mov refers to amount of moments (if they exceded or not)
    total_movements <= Movimientos ~ "Required movements",
    total_movements > Movimientos & total_movements <= Movimientos*3 ~ "More movements than required",
    total_movements > Movimientos*3 ~ "Movement excedence"))%>%
  mutate(Config_final = case_when(                                  # Labeling error types
    tipo_error == 1 ~ "Required",
    tipo_error == 2 ~ "Required",
    tipo_error == 3 ~ "Required",
    tipo_error == 4 ~ "Permuted",
    tipo_error == 5 ~ "Permuted",
    tipo_error == 6 ~ "Not required",
    tipo_error == 7 ~ "Not played",
    tipo_error == 8 ~ "Not required"))%>%
  mutate(trayec_prueba = paste(Primer_mov_codif, Config_final, sep = " "))


# Aggregating error types data

tray_error_completa%>%
  filter(tipo_error != 7)%>%
  group_by(Movimientos, Config_final)%>%
  summarise(cant_errores = n())%>%
  group_by(Movimientos)%>%
  mutate(total_errores = sum(cant_errores))%>%
  ungroup()%>%
  mutate(proporcion_errores = cant_errores/total_errores) -> tray_error_SD_GH


```
### 5.4.2. Table S16. Error types between three- and four-movement trials (Chisq. test).

```{r}

tray_error_SD_GH%>%
  
  dplyr::select(Config_final, cant_errores, Movimientos)%>%
  pivot_wider(names_from = "Config_final", values_from = cant_errores, values_fill = 0) %>% 
  as.data.frame() %>% 
  column_to_rownames("Movimientos") %>% 
  mutate(tidy(chisq.test(.))) %>% 
  dplyr::select(-parameter, -method) %>% 
  rownames_to_column("Movements") %>% 
  flextable() %>% 
    colformat_double(digits = 2)%>%
  merge_v(c("statistic" ,"p.value"))  %>%
  set_caption(
    caption = as_paragraph(
      as_chunk("Table S16. Error types between three- and four-movement trials (Chisq. test).", props = fp_text_default(bold = TRUE))
    ),
    align_with_table = F
  ) %>% 
  fix_border_issues() %>%
  autofit()
  

```

### 5.4.3. Figure S1. Error types between three- and four-movement trials.

```{r}


tray_error_SD_GH%>%
  mutate(Config_final = factor(Config_final, levels = c("Not required", "Permuted", "Required")))%>%
  ggplot(aes(x = factor(Movimientos), y = proporcion_errores, fill = Config_final))+
  geom_bar(stat = "identity", position = "stack")+
  geom_text(aes(label=ifelse(proporcion_errores >= 0.05, paste0(trunc(proporcion_errores*100^2)/10^2,"%"),"")),
            position=position_stack(vjust=0.5), colour="white")+


  coord_flip() +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(y="\nPercentaje of errors", x="Movements\n", fill = "Error type")+
  theme(panel.spacing=unit(-1, "lines"))+
  theme_minimal() ->p1

p1

```


### 5.4.4. Table S17. Error types between three-movement trials (Fisher test).

```{r}

# Subsetting data

tray_error_completa_3 = tray_error_completa%>%
  filter(Movimientos == 3)


# Aggregating data


tray_error_completa_3%>%
  filter(trial_result == 0)%>%
  filter(!is.na(tipo_error))%>%
  filter(tipo_error != 7)%>%
  group_by(Movimientos, Search.Depth, Goal.Hierarchy ,Config_final)%>%
  summarise(cant_errores = n())%>%
  group_by(Movimientos, Search.Depth, Goal.Hierarchy)%>%
  mutate(total_errores = sum(cant_errores))%>%
  ungroup()%>%
  mutate(proporcion_errores = cant_errores/total_errores)%>%
  mutate(Search.Depth = case_when(
    Search.Depth == 1 ~ "High",
    Search.Depth == 0 ~ "Low"
  ))%>%
  mutate(Goal.Hierarchy = factor(Goal.Hierarchy, levels = c("tower", "partial", "plane")))%>%
  mutate(Trial_type = interaction(Search.Depth, Goal.Hierarchy))-> tray_error_SD_GH_3


## Fisher test

tray_error_SD_GH_3%>%
  dplyr::select(Config_final, cant_errores, Trial_type)%>%
  pivot_wider(names_from = "Config_final", values_from = cant_errores, values_fill = 0) %>% 
  

  as.data.frame() %>% 
  column_to_rownames("Trial_type") %>%
  mutate(tidy(fisher.test
              (., simulate.p.value)))%>%
  dplyr::select(-alternative, -method) %>%
  rownames_to_column("Trial_type") %>%
  flextable() %>%
    colformat_double(digits = 3)%>%
  merge_v(c("p.value"))  %>%
  set_caption(
    caption = as_paragraph(
      as_chunk("Table S17. Error types between four-movement trials (Fisher test).", props = fp_text_default(bold = TRUE))
    ),
    align_with_table = F
  ) %>%
  fix_border_issues() %>%
  autofit()
  

```




### 5.4.5. Table S18. Error types between four-movement trials (Fisher test).

```{r}

# Subsetting data

tray_error_completa_4 = tray_error_completa%>%
  filter(Movimientos == 4)


# Aggregating data


tray_error_completa_4%>%
  filter(trial_result == 0)%>%
  filter(!is.na(tipo_error))%>%
  filter(tipo_error != 7)%>%
  group_by(Movimientos, Search.Depth, Goal.Hierarchy ,Config_final)%>%
  summarise(cant_errores = n())%>%
  group_by(Movimientos, Search.Depth, Goal.Hierarchy)%>%
  mutate(total_errores = sum(cant_errores))%>%
  ungroup()%>%
  mutate(proporcion_errores = cant_errores/total_errores)%>%
  mutate(Search.Depth = case_when(
    Search.Depth == 1 ~ "High",
    Search.Depth == 0 ~ "Low"
  ))%>%
  mutate(Goal.Hierarchy = factor(Goal.Hierarchy, levels = c("tower", "partial", "plane")))%>%
  mutate(Trial_type = interaction(Search.Depth, Goal.Hierarchy))-> tray_error_SD_GH_4


## Fisher test

tray_error_SD_GH_4%>%
  dplyr::select(Config_final, cant_errores, Trial_type)%>%
  pivot_wider(names_from = "Config_final", values_from = cant_errores, values_fill = 0) %>% 
  

  as.data.frame() %>% 
  column_to_rownames("Trial_type") %>%
  mutate(tidy(fisher.test
              (., simulate.p.value)))%>%
  dplyr::select(-alternative, -method) %>%
  rownames_to_column("Trial_type") %>%
  flextable() %>%
    colformat_double(digits = 3)%>%
  merge_v(c("p.value"))  %>%
  set_caption(
    caption = as_paragraph(
      as_chunk("Table S18. Error types between and four-movement trials (Fisher test).", props = fp_text_default(bold = TRUE))
    ),
    align_with_table = F
  ) %>%
  fix_border_issues() %>%
  autofit()
  

```


### 5.4.6. Figure 9. Error types in three- and four-movement trials

```{r, fig.height= 5, fig.width= 10}


tray_error_SD_GH_3%>%
  mutate(Config_final = factor(Config_final, levels = c("Not required", "Permuted", "Required")))%>%
  ggplot(aes(x = factor(Search.Depth), y = proporcion_errores, fill = Config_final))+
  geom_bar(stat = "identity", position = "stack")+
  geom_text(aes(label=ifelse(proporcion_errores >= 0.05, paste0(trunc(proporcion_errores*100^2)/10^2,"%"),"")),
            position=position_stack(vjust=0.5), colour="white")+


  coord_flip() +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(y="\nPercentaje of errors", x="Search Depth\n", fill = "Error type")+
  facet_grid(Goal.Hierarchy ~ .)+
  theme(panel.spacing=unit(-1, "lines"))+
  theme_minimal() ->p1

# p1




tray_error_SD_GH_4%>%
    mutate(Config_final = factor(Config_final, levels = c("Not required", "Permuted", "Required")))%>%
  ggplot(aes(x = factor(Search.Depth), y = proporcion_errores, fill = Config_final))+
  geom_bar(stat = "identity", position = "stack")+
  geom_text(aes(label=ifelse(proporcion_errores >= 0.05, paste0(trunc(proporcion_errores*100^2)/10^2,"%"),"")),
            position=position_stack(vjust=0.5), colour="white")+


  coord_flip() +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(y="\nPercentaje of errors", x="Search Depth\n", fill = "Error type")+
  facet_grid(Goal.Hierarchy ~ .)+
  theme(panel.spacing=unit(-1, "lines"))+
  theme_minimal() -> p2

# p2



ggarrange(p1, p2, ncol = 2, nrow = 1, labels =   c("A", "B", "C"), common.legend = T, legend = "bottom") -> figure9

figure9

# ggsave("Final_figures/Figure 9.pdf", figure9, dpi = 300, height = 5, width = 10, bg = "white")
# ggsave("Final_figures/Figure 9.png", figure9, dpi = 300, height = 5, width = 10, bg = "white")

```


```{r}

tray_error_completa %>%
  filter(trial_result == 0)%>%
  filter(!is.na(tipo_error))%>%
  filter(tipo_error != 7)%>%
  filter(Movimientos == 4) %>%
  # filter(total_movements <= 5)%>%
  mutate(Config_final = case_when(
    Config_final == "Not required" ~ "Not required/permuted",
    Config_final == "Permuted" ~ "Not required/permuted",
    Config_final == "Required" ~ "Required"
  ))%>%
  # mutate(Config_final = factor(Config_final, levels = c("No requerida", "Bola o torre permutada", "Requerida")))%>%
  mutate(Goal.Hierarchy = factor(Goal.Hierarchy, levels = c("tower", "partial", "plane")))%>%
  mutate(Config_final = factor(Config_final, levels = c("Not required/permuted", "Required")))%>%
  ggplot(aes(x=factor(Search.Depth) , y=Movement_time_error, fill = Config_final))+
  geom_boxplot()+
  facet_grid(. ~ Goal.Hierarchy)+

  theme_minimal()+
  theme(plot.margin = margin(1, 0, 0, 0, "cm"),
        strip.text = element_text(vjust = 15, margin = margin(20,1,1,1)))



tray_error_completa_4%>%
  filter(trial_result == 0)%>%
  filter(!is.na(tipo_error))%>%
  filter(tipo_error != 7)%>%
  # filter(total_movements <= 5)%>%
  mutate(Config_final = case_when(
    Config_final == "Not required" ~ "Not required/permuted",
    Config_final == "Permuted" ~ "Not required/permuted",
    Config_final == "Required" ~ "Required"
  ))%>%
  # mutate(Config_final = factor(Config_final, levels = c("No requerida", "Bola o torre permutada", "Requerida")))%>%
  mutate(Goal.Hierarchy = factor(Goal.Hierarchy, levels = c("tower", "partial", "plane")))%>%
  mutate(Config_final = factor(Config_final, levels = c("Not required/permuted", "Required")))%>%
  ggplot(aes(x=factor(Search.Depth) , y=Movement_time_error, fill = Config_final))+
  geom_boxplot()+
  facet_grid(. ~ Goal.Hierarchy)+

  theme_minimal()+
  theme(plot.margin = margin(1, 0, 0, 0, "cm"),
        strip.text = element_text(vjust = 15, margin = margin(20,1,1,1)))
  # # group_by(Movimientos, Search.Depth, Goal.Hierarchy ,Config_final)%>%
  # summarise(cant_errores = n())%>%
  # group_by(Movimientos, Search.Depth, Goal.Hierarchy)%>%
  # mutate(total_errores = sum(cant_errores))%>%
  # ungroup()%>%
  # mutate(proporcion_errores = cant_errores/total_errores)%>%
  # mutate(Search.Depth = case_when(
  #   Search.Depth == 1 ~ "High",
  #   Search.Depth == 0 ~ "Low"
  # ))%>%
  # mutate(Goal.Hierarchy = factor(Goal.Hierarchy, levels = c("tower", "partial", "plane")))%>%
  # mutate(Trial_type = interaction(Search.Depth, Goal.Hierarchy))-> tray_error_SD_GH_3

```


### 5.4.7 Table S19. Error response times between different error types (three-movement trials)

```{r}

tray_error_completa_3 %>%
 filter(trial_result == 0)%>%
  filter(!is.na(tipo_error))%>%
  filter(tipo_error != 7)%>%
  filter(Movimientos == 3) %>%
  mutate(Config_final = case_when(
    Config_final == "Not required" ~ "Not required/permuted",
    Config_final == "Permuted" ~ "Not required/permuted",
    Config_final == "Required" ~ "Required"
  ))%>%
  mutate(Goal.Hierarchy = factor(Goal.Hierarchy, levels = c("tower", "partial", "plane")))%>%
  mutate(Config_final = factor(Config_final, levels = c("Not required/permuted", "Required")))->final_data_3


pvalues_planning <- final_data_3 %>%
  group_by(tipo_ensayo) %>%
  group_map(~ {
    test_result <- wilcox.test(Planning_time_error ~ Config_final, data = .x)
    tidy(test_result)
  }) %>%
  bind_rows(.id = "tipo_ensayo") %>% 
  select(tipo_ensayo, p.value) %>% 
  mutate(tipo_ensayo = sort(unique(final_data_3$tipo_ensayo)),
         Variable = c(rep("Planning_time_error", 4 ))) %>% 
  select(Variable, everything())


pvalues_movement <- final_data_3 %>%
  group_by(tipo_ensayo) %>%
  group_map(~ {
    test_result <- wilcox.test(Movement_time_error ~ Config_final, data = .x)
    tidy(test_result)
  }) %>%
  bind_rows(.id = "tipo_ensayo") %>% 
  select(tipo_ensayo, p.value) %>% 
  mutate(tipo_ensayo = sort(unique(final_data_3$tipo_ensayo)),
         Variable = c(rep("Movement_time_error", 4 ))) %>% 
  select(Variable, everything())


wilcox_results = rbind(pvalues_planning, pvalues_movement)


effsizes_planning <- final_data_3 %>%
  group_by(tipo_ensayo) %>%
  group_map(~ {
    test_result <- wilcoxonR(.x$Planning_time_error,
                             .x$Config_final, data = .x)
    tidy(test_result)
  }) %>%
  bind_rows(.id = "tipo_ensayo")%>% 
  select(tipo_ensayo, x) %>% 
  rename("Effect_size" = x) %>% 
  mutate(tipo_ensayo = sort(unique(final_data_3$tipo_ensayo)),
         Variable = c(rep("Planning_time_error", 4 ))) %>% 
  select(Variable, everything())


effsizes_movement <- final_data_3 %>%
  group_by(tipo_ensayo) %>%
  group_map(~ {
    test_result <- wilcoxonR(.x$Movement_time_error,
                             .x$Config_final, data = .x)
    tidy(test_result)
  }) %>%
  bind_rows(.id = "tipo_ensayo")%>% 
  select(tipo_ensayo, x) %>% 
  rename("Effect_size" = x) %>% 
  mutate(tipo_ensayo = sort(unique(final_data_3$tipo_ensayo)),
         Variable = c(rep("Movement_time_error", 4 ))) %>% 
  select(Variable, everything())


effsizes_results = rbind(effsizes_planning, effsizes_movement)


wilcox_results = wilcox_results %>% left_join(effsizes_results)
  
wilcox_results %>% 
  flextable() %>%
    colformat_double(digits = 3)%>%
  merge_v(c("p.value"))  %>%
  set_caption(
    caption = as_paragraph(
      as_chunk("Table S19. Error response times between different error types (three-movement trials)", props = fp_text_default(bold = TRUE))
    ),
    align_with_table = F
  ) %>%
  fix_border_issues() %>%
  autofit()



```



### 5.4.8. Table S20. Error response times between different error types (four-movement trials)

```{r}

tray_error_completa_4 %>%
 filter(trial_result == 0)%>%
  filter(!is.na(tipo_error))%>%
  filter(tipo_error != 7)%>%
  filter(Movimientos == 4) %>%
  mutate(Config_final = case_when(
    Config_final == "Not required" ~ "Not required/permuted",
    Config_final == "Permuted" ~ "Not required/permuted",
    Config_final == "Required" ~ "Required"
  ))%>%
  mutate(Goal.Hierarchy = factor(Goal.Hierarchy, levels = c("tower", "partial", "plane")))%>%
  mutate(Config_final = factor(Config_final, levels = c("Not required/permuted", "Required")))->final_data_4


pvalues_planning <- final_data_4 %>%
  group_by(tipo_ensayo) %>%
  group_map(~ {
    test_result <- wilcox.test(Planning_time_error ~ Config_final, data = .x)
    tidy(test_result)
  }) %>%
  bind_rows(.id = "tipo_ensayo") %>% 
  select(tipo_ensayo, p.value) %>% 
  mutate(tipo_ensayo = sort(unique(final_data_4$tipo_ensayo)),
         Variable = c(rep("Planning_time_error", 4 ))) %>% 
  select(Variable, everything())


pvalues_movement <- final_data_4 %>%
  group_by(tipo_ensayo) %>%
  group_map(~ {
    test_result <- wilcox.test(Movement_time_error ~ Config_final, data = .x)
    tidy(test_result)
  }) %>%
  bind_rows(.id = "tipo_ensayo") %>% 
  select(tipo_ensayo, p.value) %>% 
  mutate(tipo_ensayo = sort(unique(final_data_4$tipo_ensayo)),
         Variable = c(rep("Movement_time_error", 4 ))) %>% 
  select(Variable, everything())


wilcox_results = rbind(pvalues_planning, pvalues_movement)


effsizes_planning <- final_data_4 %>%
  group_by(tipo_ensayo) %>%
  group_map(~ {
    test_result <- wilcoxonR(.x$Planning_time_error,
                             .x$Config_final, data = .x)
    tidy(test_result)
  }) %>%
  bind_rows(.id = "tipo_ensayo")%>% 
  select(tipo_ensayo, x) %>% 
  rename("Effect_size" = x) %>% 
  mutate(tipo_ensayo = sort(unique(final_data_4$tipo_ensayo)),
         Variable = c(rep("Planning_time_error", 4 ))) %>% 
  select(Variable, everything())


effsizes_movement <- final_data_4 %>%
  group_by(tipo_ensayo) %>%
  group_map(~ {
    test_result <- wilcoxonR(.x$Movement_time_error,
                             .x$Config_final, data = .x)
    tidy(test_result)
  }) %>%
  bind_rows(.id = "tipo_ensayo")%>% 
  select(tipo_ensayo, x) %>% 
  rename("Effect_size" = x) %>% 
  mutate(tipo_ensayo = sort(unique(final_data_4$tipo_ensayo)),
         Variable = c(rep("Movement_time_error", 4 ))) %>% 
  select(Variable, everything())


effsizes_results = rbind(effsizes_planning, effsizes_movement)


wilcox_results = wilcox_results %>% left_join(effsizes_results)
  
wilcox_results %>% 
  flextable() %>%
    colformat_double(digits = 3)%>%
  merge_v(c("p.value"))  %>%
  set_caption(
    caption = as_paragraph(
      as_chunk("Table S20. Error response times between different error types (four-movement trials)", props = fp_text_default(bold = TRUE))
    ),
    align_with_table = F
  ) %>%
  fix_border_issues() %>%
  autofit()



```

### 5.4.9. Figure 10. Error planning and movement times between error types in four-movement trials.

```{r, fig.height= 8, fig.width= 10}


## Plot theming

customize = function(myplot){
  
  myplot +
    theme_minimal(base_size = 15)+
    theme(plot.margin = margin(0, 0.25, 0.75, 0.25, "cm"),
          strip.text = element_text(vjust = -2, margin = margin(0,1,1,1)),
          axis.title.x = element_text(vjust = - 1)) 
  
  
}

tray_error_completa %>%
  filter(trial_result == 0)%>%
  filter(!is.na(tipo_error))%>%
  filter(tipo_error != 7)%>%
  filter(Movimientos == 4) %>%
  
  filter(Cantidad_mov != "Exceso movimientos")%>%
 mutate(Config_final = case_when(
    Config_final == "Not required" ~ "Not required/permuted",
    Config_final == "Permuted" ~ "Not required/permuted",
    Config_final == "Required" ~ "Required"
  ))%>%
    mutate(Goal.Hierarchy = case_when(
  Goal.Hierarchy == "tower" ~ "Tower\n hierarchy\n ",
  Goal.Hierarchy == "partial" ~ "Partial\n hierarchy\n ",
  Goal.Hierarchy == "plane" ~ "Plane\n hierarchy\n "))%>%
  mutate(Search.Depth = case_when(
  Search.Depth == "0" ~ "Low",
  Search.Depth == "1" ~ "High"))%>%
  mutate(Goal.Hierarchy = factor(Goal.Hierarchy, levels = c("Tower\n hierarchy\n ", "Partial\n hierarchy\n ", "Plane\n hierarchy\n ")))%>%
  mutate(Search.Depth = factor(Search.Depth, levels = c("Low", "High")))%>%

  mutate(Config_final = factor(Config_final, levels = c("Not required/permuted", "Required")))%>%


  ggplot(aes(x=factor(Search.Depth) , y=Planning_time_error, fill = Config_final))+
  geom_boxplot()+
  facet_grid(. ~ Goal.Hierarchy)+

  theme_minimal() %>% 
  customize()+
  labs(y= "\nPlanning error times", x= "Search Depth\n", fill = "Error type")->p1





tray_error_completa %>%
  filter(trial_result == 0)%>%
  filter(!is.na(tipo_error))%>%
  filter(tipo_error != 7)%>%
  filter(Movimientos == 4) %>%
  
  filter(Cantidad_mov != "Exceso movimientos")%>%
 mutate(Config_final = case_when(
    Config_final == "Not required" ~ "Not required/permuted",
    Config_final == "Permuted" ~ "Not required/permuted",
    Config_final == "Required" ~ "Required"
  ))%>%
  mutate(Goal.Hierarchy = case_when(
  Goal.Hierarchy == "tower" ~ "Tower\n hierarchy\n ",
  Goal.Hierarchy == "partial" ~ "Partial\n hierarchy\n ",
  Goal.Hierarchy == "plane" ~ "Plane\n hierarchy\n "))%>%
  mutate(Search.Depth = case_when(
  Search.Depth == "0" ~ "Low",
  Search.Depth == "1" ~ "High"))%>%
  mutate(Goal.Hierarchy = factor(Goal.Hierarchy, levels = c("Tower\n hierarchy\n ", "Partial\n hierarchy\n ", "Plane\n hierarchy\n ")))%>%
  mutate(Search.Depth = factor(Search.Depth, levels = c("Low", "High")))%>%


  ggplot(aes(x=factor(Search.Depth) , y=Movement_time_error, fill = Config_final))+
  geom_boxplot()+
  facet_grid(. ~ Goal.Hierarchy)+

  theme_minimal() %>% 
  customize()+
  labs(y= "\nMovement error times", x= "Search Depth\n", fill = "Error type")   -> p2



ggarrange(p1, p2, ncol = 1, nrow = 2, labels =   c("A", "B"), common.legend = T, legend = "bottom", font.label = list(size = 18))-> figure10

figure10

# ggsave("Final_figures/Figure 10.pdf", figure10, dpi = 300, height = 8, width = 10, bg = "white")

```

