```{r, echo=FALSE, message=FALSE, warning=FALSE}
# packages
library(ez)
library(reshape2)
library(reshape)
library(ggplot2)
library(pastecs)
library(plyr)
library(ez)
library(data.table)
library(tidyverse)
library(readxl)
library(readr)
library(lme4)
library(RColorBrewer)
library(lattice)
library(survival)
library(Formula)
library(ggpubr)
library(scales)
library(lme4)
library(multcomp)
library(zoo)
library(car)
library(emmeans)

library(showtext)
font_add_google("Poppins", "Poppins")
font_add_google("Roboto Mono", "Roboto Mono")
showtext_auto()
```


# Data preparation

```{r, message=FALSE, warning=FALSE}

data <- read_csv("Exp_data.csv")
```

```{r}
# delete trials 1 to 4 

data <- subset(data, trialNr > 4)


#new column with correct answer text, trial type, likelihood ratio, and strategies used
data <- data %>% 
  mutate(correctAns2 = if_else(correctAns == 1, "Correct", "Wrong"),
         urnSize = ifelse(trialNr %% 2 == 0, "Big", "Small"),
         LR = if_else(trialNr<7, "Infinite",
                         if_else(trialNr<15 & trialNr>6, "LR 4","LR 1.5")),
         trialType = if_else(trialNr<7, "Control",
                            if_else(trialNr>6 & trialNr<9,"A", 
                                    if_else(trialNr>14 & trialNr<17, "A", 
                                            if_else(trialNr>8 & trialNr<11, "B",
                                                    if_else(trialNr>16 & trialNr<19, "B", 
                                                            if_else(trialNr>10 & trialNr<13, "C", 
                                                                    if_else(trialNr>18 & trialNr<21, "C", "D"))))))))


data <- subset(data, select = c(sID_new, age, sex, trialNr, correctAns, correctAns2, urnSize, LR, trialType))

data_withC <- data

data <- subset(data, trialType != "Control")

data$trialOrd <- rep(1:16,90)
data$trialOrd2 <- rep(c("First_half","First_half","First_half","First_half","First_half","First_half","First_half","First_half","Second_half","Second_half","Second_half","Second_half","Second_half","Second_half","Second_half","Second_half"),90)



```


```{r}
data$age <- factor(data$age, levels = c(3:5), labels = c("3 years", "4 years", "5 years"))

data$correctAns2 <- factor(data$correctAns2, levels = c("Wrong", "Correct"))

data$trialType <- factor(data$trialType, levels = c("A", "B", "C", "D"))

data$LR <- factor(data$LR, levels = c("LR 1.5", "LR 4"))

data$trialOrd2 <- factor(data$trialOrd2, levels = c("First_half", "Second_half"))

```


# Initial Step: Check performance in the control trials

```{r}
data_control <- subset(data_withC, trialType == "Control")


proptest_data <- data_control %>%
  group_by(age, correctAns2) %>%
  summarize(n = n()) %>% 
  mutate(pct = n/sum(n),
         lbl = scales::percent(pct))

```

Against chance-level performance:

```{r}
age3 <- prop.test(proptest_data$n[1], proptest_data$n[1] + proptest_data$n[2], p = 0.5)
age3

age4 <- prop.test(proptest_data$n[3], proptest_data$n[3] + proptest_data$n[4], p = 0.5)
age4

age5 <- prop.test(proptest_data$n[5], proptest_data$n[5] + proptest_data$n[6], p = 0.5)
age5
```

Very high accuracy in all age groups in the control trials. 


# Analyses of test trials


## Graphs


Define a theme:

```{r}
myTheme <- theme(plot.title = element_text(face="bold", size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20),
        axis.text.x = element_text(size = 14, angle = 0), 
        axis.text.y = element_text(size = 16, angle = 0),
        legend.text = element_text(size = 18),
        legend.title = element_text(face = "bold", size = 18),
        strip.text.x = element_text(size = 18),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(), 
        axis.line.x = element_line(colour = "black"), 
        axis.line.y = element_line(colour = "black"),
        axis.text = element_text(colour ="black"), 
        axis.ticks = element_line(colour ="black"))

```



Make the plots:

1. Age overall

```{r, warning=FALSE}
# create a summary dataset that also contains the percentages
# AGE

plotdata <- data %>%
  group_by(sID_new, age, correctAns2) %>%
  summarize(n = n()) #%>% 
  #mutate(pct = n/sum(n),
   #      lbl = scales::percent(pct))

plotdata_age_correct <- subset(plotdata, correctAns2 == "Correct")
plotdata_age_wrong <- subset(plotdata, correctAns2 == "Wrong" & n == 16)

plotdata_age_wrong[plotdata_age_wrong=="Wrong"] <- "Correct"
plotdata_age_wrong[plotdata_age_wrong==16] <- 0

plotdata_age <- rbind(plotdata_age_correct, plotdata_age_wrong)
```


```{r, warning=FALSE}

plotdat <- plotdata_age

theme_set(theme_light(base_size = 20, base_family = "Poppins"))


library(ggdist)
library(ggridges)

g_age <- ggplot(plotdat, aes(x = n, y = age, fill = age)) +
  #facet_grid( ~ LR,labeller = label_both)+
  ggtitle("Performance by age") +
  scale_x_continuous(limits = c(-1, 17), breaks=seq(0, 16, 1), expand = c(0,0))+
  stat_slab(aes(thickness = after_stat(pdf*n)), scale = 0.8, alpha = 0.7, fill = "#66c2a5") +
  stat_dotsinterval(side = "bottom", scale = 0.7, slab_linewidth = NA, point_interval = "mean_qi", color = "black", interval_alpha = 0, 
                    slab_alpha = 0.7, point_alpha = 0, fill = "#66c2a5") +
  #stat_dotsinterval()+
  #geom_density_ridges(alpha = 0.5)+
   #stat_summary(aes(x = rating_rec), fun.x=mean, geom="point", 
  #             color = "black", shape = 22, size = 2, group=1, alpha = 1)+
  stat_summary(aes(y = age, color = age), fun.data = mean_cl_boot, 
               geom = "errorbar", width = 0, size = 1, color = "black") +
  stat_summary(aes(y = age, fill = age), fun.y=mean, geom="point", 
               color = "black", shape = 22, size = 3, group=1, alpha = 1, fill = "black")+
  scale_fill_manual(values=c("#66c2a5", "#e78ac3", "#8da0cb", "#a6d854"))+
  scale_color_manual(values=c("#66c2a5", "#e78ac3", "#8da0cb", "#a6d854"))+
  #scale_fill_viridis_c(name = "Explanation \nRating", option = "C", breaks=c(-5,0,5), labels=c("narrow scope", "no preference", "broad scope"))+
  labs(x = "Number of correct choices", y = "Age group") +
  stat_summary(aes(label=round(after_stat(x),2)), fun.y=mean, geom="text", size=5.5,
             vjust = -1)+
  #scale_y_discrete(limits=rev)+
  myTheme+
  theme_light(base_family = "Poppins", base_size = 15)+
  theme(panel.grid = element_blank(), axis.text = element_text(colour ="black"))+
  theme(legend.position="none",
        legend.title=element_blank(),legend.key.width = unit(1.95, 'cm'))+
  theme(axis.text.y = element_text(size = 14, angle = 0))+
  annotate('text', x = 10, y = 3.5, label = '80.4%',size = 5.5)+
  annotate('text', x = 7, y = 2.5, label = '70.4%',size = 5.5)+
  annotate('text', x = 4, y = 1.3, label = '63.6%',size = 5.5)
  #+
   #annotate('curve', x = 13, y = 3.4, yend = 3.5, xend = 11.3, linewidth = 0.8, curvature = 0.2, arrow = arrow(length = unit(0.2, 'cm')))
  

g_age

ggsave("PScores_dist_age_bw.svg",width=6.5,height=5)
ggsave("PScores_dist_age_bw.pdf",width=6.5,height=5)
```



2. LR

```{r, warning=FALSE}

plotdata <- data %>%
  group_by(sID_new, age, LR, correctAns2) %>%
  summarize(n = n()) #%>% 
  #mutate(pct = n/sum(n),
   #      lbl = scales::percent(pct))

plotdata_LR_correct <- subset(plotdata, correctAns2 == "Correct")
plotdata_LR_wrong <- subset(plotdata, correctAns2 == "Wrong" & n == 8)

plotdata_LR_wrong[plotdata_LR_wrong=="Wrong"] <- "Correct"
plotdata_LR_wrong[plotdata_LR_wrong==8] <- 0

plotdata_LR <- rbind(plotdata_LR_correct, plotdata_LR_wrong)
```

```{r, warning=FALSE}

plotdat <- plotdata_LR

theme_set(theme_light(base_size = 20, base_family = "Poppins"))


library(ggdist)
library(ggridges)
g_LR <- ggplot(plotdat, aes(x = n, y = LR)) +
  ggtitle("Performance by LR") +
  #facet_grid( ~ age, labeller = label_both)+
  scale_x_continuous(limits = c(-1, 9), breaks=seq(0, 8, 1), expand = c(0,0))+
  stat_slab(aes(thickness = after_stat(pdf*n)), scale = 0.8, alpha = 0.7, fill = "#e78ac3") +
  stat_dotsinterval(side = "bottom", scale = 0.7, slab_linewidth = NA, point_interval = "mean_qi", color = "black", interval_alpha = 0, 
                    slab_alpha = 0.7, point_alpha = 0, fill = "#e78ac3") +
  #stat_dotsinterval()+
  #geom_density_ridges(alpha = 0.5)+
   #stat_summary(aes(x = rating_rec), fun.x=mean, geom="point", 
  #             color = "black", shape = 22, size = 2, group=1, alpha = 1)+
  stat_summary(aes(y = LR), fun.data = mean_cl_boot, 
               geom = "errorbar", width = 0, size = 1, color = "black") +
  stat_summary(aes(y = LR), fun.y=mean, geom="point", 
               color = "black", shape = 22, size = 3, group=1, alpha = 1, fill = "black")+
  scale_fill_manual(values=c("#66c2a5", "#e78ac3", "#8da0cb", "#a6d854"))+
  scale_color_manual(values=c("#66c2a5", "#e78ac3", "#8da0cb", "#a6d854"))+
  #scale_fill_viridis_c(name = "Explanation \nRating", option = "C", breaks=c(-5,0,5), labels=c("narrow scope", "no preference", "broad scope"))+
  labs(x = "Number of correct choices", y = "Likelihood ratio") +
  stat_summary(aes(label=round(after_stat(x),2)), fun.y=mean, geom="text", size=5.5,
             vjust = -1)+
  #scale_y_discrete(limits=rev)+
  myTheme+
  theme_light(base_family = "Poppins", base_size = 15)+
  theme(panel.grid = element_blank(), axis.text = element_text(colour ="black"))+
  theme(legend.position="none",
        legend.title=element_blank(),legend.key.width = unit(1.95, 'cm'))+
  theme(axis.text.y = element_text(size = 14, angle = 0))+
  annotate('text', x = 1, y = 2.5, label = '71.6%',size = 5.5)+
  annotate('text', x = 1, y = 1.3, label = '71.3%',size = 5.5)
  
  

g_LR

ggsave("PScores_dist_LR.svg",width=6,height=4)
ggsave("PScores_dist_LR.pdf",width=6,height=4)
```


3. Trial Type

```{r, warning=FALSE}
plotdata <- data %>%
  group_by(sID_new, age, trialType, correctAns2) %>%
  summarize(n = n()) #%>% 
  #mutate(pct = n/sum(n),
   #      lbl = scales::percent(pct))

plotdata_trial_correct <- subset(plotdata, correctAns2 == "Correct")
plotdata_trial_wrong <- subset(plotdata, correctAns2 == "Wrong" & n == 4)

plotdata_trial_wrong[plotdata_trial_wrong=="Wrong"] <- "Correct"
plotdata_trial_wrong[plotdata_trial_wrong==4] <- 0

plotdata_trial <- rbind(plotdata_trial_correct, plotdata_trial_wrong)
```


```{r, warning=FALSE}
plotdat <- plotdata_trial

theme_set(theme_light(base_size = 20, base_family = "Poppins"))

colnames(plotdat) <- c("sID_new", "age", "Trial type", "correctAns2", "n") 


library(ggdist)
library(ggridges)
g_TT <- ggplot(plotdat, aes(x = n, y = `Trial type`)) +
  ggtitle("Performance by trial type")+
  #facet_grid( ~ age, labeller = label_both)+
  scale_x_continuous(limits = c(-1, 5), breaks=seq(0, 4, 1), expand = c(0,0))+
  stat_slab(aes(thickness = after_stat(pdf*n)), scale = 0.8, alpha = 0.7, fill = "#8da0cb") +
  stat_dotsinterval(side = "bottom", scale = 0.7, slab_linewidth = NA, point_interval = "mean_qi", color = "black", interval_alpha = 0, 
                    slab_alpha = 0.7, point_alpha = 0, fill = "#8da0cb") +
  #stat_dotsinterval()+
  #geom_density_ridges(alpha = 0.5)+
   #stat_summary(aes(x = rating_rec), fun.x=mean, geom="point", 
  #             color = "black", shape = 22, size = 2, group=1, alpha = 1)+
  stat_summary(aes(y = `Trial type`), fun.data = mean_cl_boot, 
               geom = "errorbar", width = 0, size = 1, color = "black") +
  stat_summary(aes(y = `Trial type`), fun.y=mean, geom="point", 
               color = "black", shape = 22, size = 3, group=1, alpha = 1, fill ="black")+
  scale_fill_manual(values=c("#66c2a5", "#e78ac3", "#8da0cb", "#a6d854"))+
  scale_color_manual(values=c("#66c2a5", "#e78ac3", "#8da0cb", "#a6d854"))+
  #scale_fill_viridis_c(name = "Explanation \nRating", option = "C", breaks=c(-5,0,5), labels=c("narrow scope", "no preference", "broad scope"))+
  labs(x = "Number of correct choices", y = "Trial type") +
  stat_summary(aes(label=round(after_stat(x),2)), fun.y=mean, geom="text", size=5.5,
             vjust = -1)+
  scale_y_discrete(limits=rev)+
  myTheme+
  theme_light(base_family = "Poppins", base_size = 15)+
  theme(panel.grid = element_blank(), axis.text = element_text(colour ="black"))+
  theme(legend.position="none",
        legend.title=element_blank(),legend.key.width = unit(1.95, 'cm'))+
  theme(axis.text.y = element_text(size = 14, angle = 0))+  
  annotate('text', x = 0.5, y = 4.3, label = '68.3%',size = 5.5)+
  annotate('text', x = 0.5, y = 3.3, label = '75.5%',size = 5.5)+
  annotate('text', x = 0.5, y = 2.3, label = '66.8%',size = 5.5)+
  annotate('text', x = 0.5, y = 1.3, label = '75.25%',size = 5.5)
  

g_TT

ggsave("PScores_dist_trialType.svg",width=6,height=4)
ggsave("PScores_dist_trialType.pdf",width=6,height=4)
```


4. Combine the plots


```{r, warning=FALSE, message= FALSE}
library(ggpubr)

figure1 <- ggarrange(g_age, 
                    ggarrange(g_LR, g_TT, labels = c("B", "C"), ncol=2, font.label = list(size = 20, color = "black")),
                    nrow = 2, 
                    labels = "A",
                    font.label = list(size = 20, color = "black"))
figure1

ggsave("ExpRes_comb_bw.svg",width=10,height=10)
ggsave("ExpRes_comb_bw.pdf",width=10,height=10)


```

## Statistical Analyses 

### Overall effect of age

ANOVA:

```{r, warning=FALSE}
library(afex)
library(emmeans)

a1 <- aov_car(n ~ age + Error(sID_new/1), plotdata_age, 
              anova_table = list(es = "pes"))
a1

```

Polynomial contrast to test the trends:

```{r, warning=FALSE}
contrasts(plotdata_age$age) <- "contr.poly"

LinearModel.3 <- lm(n ~ age, data=plotdata_age)
summary(LinearModel.3)
```

Use emmeans to get the individual means and their CIs:

```{r, warning=FALSE}

library(lsmeans)
# means

ls2 <- lsmeans(a1, c("age")) 
ls2

contrasts <- emmeans(a1, ~ age)
s <- pairs(contrasts, adjust = "none") 


s
confint(s, level = 0.95)
```

## Proportion tests 


```{r, warning=FALSE}
proptest_data <- data %>%
  group_by(age, correctAns2) %>%
  summarize(n = n()) %>% 
  mutate(pct = n/sum(n),
         lbl = scales::percent(pct))

```

Against chance-level performance:

```{r, warning=FALSE}
age3 <- prop.test(proptest_data$n[2], proptest_data$n[1] + proptest_data$n[2], p = 0.5)
age3

age4 <- prop.test(proptest_data$n[4], proptest_data$n[3] + proptest_data$n[4], p = 0.5)
age4

age5 <- prop.test(proptest_data$n[6], proptest_data$n[5] + proptest_data$n[6], p = 0.5)
age5
```

Against the best heuristic strategy (More good and Less bad predict 75% correct urns): 

3- and 4-year-olds have a performance below 0.75, but the 5-year-olds are better. Is it significantly better?

```{r, warning=FALSE}
age5 <- prop.test(proptest_data$n[6], proptest_data$n[5] + proptest_data$n[6], p = 0.75, alternative = "greater")
age5
```

## GLMM 

to test also the other factors included in the experiment. 

```{r, warning=FALSE}

data$correctAns <- as.numeric(data$correctAns)

contr=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun=100000)) # settings needed for model to converge


# null model
null = glmer(correctAns ~ 1 + (1|sID_new), data=data, family="binomial", control = contr)

# add age
age = glmer(correctAns ~ age + (1|sID_new), data=data, family="binomial", control = contr)

# add urnsize
urnsize = glmer(correctAns ~ age + urnSize + (1|sID_new), data=data, family="binomial", control = contr)

# add LR 
LR = glmer(correctAns ~ age + urnSize + LR + (1|sID_new), data=data, family="binomial", control = contr)

# add trial type
trialType = glmer(correctAns ~ age + urnSize + LR + trialType + (1|sID_new), data=data, family="binomial", control = contr)


# add interaction between LR and trial type
interTTLR = glmer(correctAns ~ age + urnSize + LR + trialType + trialType*LR + (1|sID_new), data=data, family="binomial", control = contr)

# add interaction between LR and age
interLRage = glmer(correctAns ~ age + urnSize + LR + trialType + trialType*LR + LR*age + (1|sID_new), data=data, family="binomial", control = contr)

# add interaction between TT and age
interTTage = glmer(correctAns ~ age + urnSize + LR + trialType + trialType*LR + LR*age + trialType*age + (1|sID_new), data=data, family="binomial", control = contr)



# add interaction between LR and trial type and age
interTTLRage = glmer(correctAns ~ age + urnSize + LR + trialType + trialType*LR + LR*age + trialType*LR*age + (1|sID_new), data=data, family="binomial", control = contr)



# test models
anova(null, age, urnsize, LR, trialType, interTTLR, interLRage, interTTage, interTTLRage)
```

In addition to age, trial type, the interaction between LR and trial type is significant. 

Follow up on the interaction by looking at a graph:

```{r, warning=FALSE}
plotdata <- data %>%
  group_by(sID_new, LR, trialType, correctAns2) %>%
  summarize(n = n()) #%>% 
  #mutate(pct = n/sum(n),
   #      lbl = scales::percent(pct))


plotdata_trialLR_correct <- subset(plotdata, correctAns2 == "Correct")
plotdata_trialLR_wrong <- subset(plotdata, correctAns2 == "Wrong" & n == 2)

plotdata_trialLR_wrong[plotdata_trialLR_wrong=="Wrong"] <- "Correct"
plotdata_trialLR_wrong[plotdata_trialLR_wrong==2] <- 0

plotdata_trialLR <- rbind(plotdata_trialLR_correct, plotdata_trialLR_wrong)
```




```{r, warning=FALSE}
plotdat <- plotdata_trialLR

theme_set(theme_light(base_size = 20, base_family = "Poppins"))

colnames(plotdat) <- c("sID_new", "LR", "Trial type", "correctAns2", "n") 



# annotations LR 1.5  
ann_text_LR15_A <- data.frame(n = 2.5,`Trial type`=4.4,lab = "Text",
                       LR = factor("LR 1.5",levels = c("LR 1.5", "LR 4")))
colnames(ann_text_LR15_A)[2] <- "Trial type"

ann_text_LR15_B <- data.frame(n = 2.5,`Trial type`=3.4,lab = "Text",
                       LR = factor("LR 1.5",levels = c("LR 1.5", "LR 4")))
colnames(ann_text_LR15_B)[2] <- "Trial type"

ann_text_LR15_C <- data.frame(n = 2.5,`Trial type`=2.4,lab = "Text",
                       LR = factor("LR 1.5",levels = c("LR 1.5", "LR 4")))
colnames(ann_text_LR15_C)[2] <- "Trial type"

ann_text_LR15_D <- data.frame(n = 2.5,`Trial type`=1.4,lab = "Text",
                       LR = factor("LR 1.5",levels = c("LR 1.5", "LR 4")))
colnames(ann_text_LR15_D)[2] <- "Trial type"


# annotations LR 4 
ann_text_LR4_A <- data.frame(n = 2.5,`Trial type`=4.4,lab = "Text",
                       LR = factor("LR 4",levels = c("LR 1.5", "LR 4")))
colnames(ann_text_LR4_A)[2] <- "Trial type"

ann_text_LR4_B <- data.frame(n = 2.5,`Trial type`=3.4,lab = "Text",
                       LR = factor("LR 4",levels = c("LR 1.5", "LR 4")))
colnames(ann_text_LR4_B)[2] <- "Trial type"

ann_text_LR4_C <- data.frame(n = 2.5,`Trial type`=2.4,lab = "Text",
                       LR = factor("LR 4",levels = c("LR 1.5", "LR 4")))
colnames(ann_text_LR4_C)[2] <- "Trial type"

ann_text_LR4_D <- data.frame(n = 2.5,`Trial type`=1.4,lab = "Text",
                       LR = factor("LR 4",levels = c("LR 1.5", "LR 4")))
colnames(ann_text_LR4_D)[2] <- "Trial type"



library(ggdist)
library(ggridges)
g_TTLR <- ggplot(plotdat, aes(x = n, y = `Trial type`, fill = LR)) +
  ggtitle("Performance by trial type and LR")+
  facet_grid( ~ LR)+
  scale_x_continuous(limits = c(-1, 3), breaks=seq(0, 2, 1), expand = c(0,0))+
  stat_slab(aes(thickness = after_stat(pdf*n)), scale = 0.8, alpha = 0.7) +
  stat_dotsinterval(side = "bottom", scale = 0.7, slab_linewidth = NA, point_interval = "mean_qi", color = "black", interval_alpha = 0, 
                    slab_alpha = 0.7, point_alpha = 0) +
  #stat_dotsinterval()+
  #geom_density_ridges(alpha = 0.5)+
   #stat_summary(aes(x = rating_rec), fun.x=mean, geom="point", 
  #             color = "black", shape = 22, size = 2, group=1, alpha = 1)+
  stat_summary(aes(y = `Trial type`, color = LR), fun.data = mean_cl_boot, 
               geom = "errorbar", width = 0, size = 1, color = "black") +
  stat_summary(aes(y = `Trial type`, fill = LR), fun.y=mean, geom="point", 
               color = "black", shape = 22, size = 3, group=1, alpha = 1, fill = "black")+
  scale_fill_manual(values=c("#8da0cb", "#8da0cb", "#8da0cb", "#a6d854"))+
  scale_color_manual(values=c("#8da0cb", "#8da0cb", "#8da0cb", "#a6d854"))+
  #scale_fill_viridis_c(name = "Explanation \nRating", option = "C", breaks=c(-5,0,5), labels=c("narrow scope", "no preference", "broad scope"))+
  labs(x = "Number of correct choices", y = "Trial type") +
  stat_summary(aes(label=round(after_stat(x),2)), fun.y=mean, geom="text", size=5.5,
             vjust = -1)+
  scale_y_discrete(limits=rev)+
  myTheme+
  theme_light(base_family = "Poppins", base_size = 15)+
  theme(panel.grid = element_blank(), axis.text = element_text(colour ="black"))+
  theme(legend.position="none",
        legend.title=element_blank(),legend.key.width = unit(1.95, 'cm'))+
  theme(axis.text.y = element_text(size = 14, angle = 0))+
  geom_text(data = ann_text_LR15_A,label = "62.0%", size = 3.5)+ # LR 15
  geom_text(data = ann_text_LR15_B,label = "77.0%", size = 3.5)+
  geom_text(data = ann_text_LR15_C,label = "66.0%", size = 3.5)+
  geom_text(data = ann_text_LR15_D,label = "79.5%", size = 3.5)+
  geom_text(data = ann_text_LR4_A,label = "74.5%", size = 3.5)+ # LR 4
  geom_text(data = ann_text_LR4_B,label = "74.0%", size = 3.5)+
  geom_text(data = ann_text_LR4_C,label = "67.0%", size = 3.5)+
  geom_text(data = ann_text_LR4_D,label = "71.0%", size = 3.5)

g_TTLR

ggsave("PScores_dist_interTTLR.svg",width=6,height=4)
ggsave("PScores_dist_interTTLR.pdf",width=6,height=4)
```


A and C trials were harder, but only for the more difficult LR


Do the tests (contrasts): 

```{r, warning=FALSE}

library(lsmeans)
# means

ls2 <- lsmeans(interTTLR, c("age")) 
ls2

contrasts <- emmeans(interTTLR, ~ trialType|LR)
s <- pairs(contrasts, adjust = "none") 


s
confint(s, level = 0.95)
```




## Exploratory analyses 

Might there be a learning effect in the different age groups? 

```{r, warning=FALSE}

trialOrder = glmer(correctAns ~ age + urnSize + LR + trialType + trialType*LR + trialOrd + (1|sID_new), data=data, family="binomial", control = contr) # take out the interactions that were not a sig. improvement 

Inter_trialOrderAge = glmer(correctAns ~ age + urnSize + LR + trialType + trialType*LR + trialOrd + trialOrd*age + (1|sID_new), data=data, family="binomial", control = contr) # take out the interactions that were not a sig. improvement 

trialOrder2 = glmer(correctAns ~ age + urnSize + LR + trialType + trialType*LR + trialOrd2 + (1|sID_new), data=data, family="binomial", control = contr) # take out the interactions that were not a sig. improvement 

Inter_trialOrder2Age = glmer(correctAns ~ age + urnSize + LR + trialType + trialType*LR + trialOrd2 + trialOrd2*age + (1|sID_new), data=data, family="binomial", control = contr) # take out the interactions that were not a sig. improvement 


# test models
anova(null, age, urnsize, LR, trialType, interTTLR, trialOrder)

anova(null, age, urnsize, LR, trialType, interTTLR, trialOrder, Inter_trialOrderAge)

anova(null, age, urnsize, LR, trialType, interTTLR, trialOrder2)

anova(null, age, urnsize, LR, trialType, interTTLR, trialOrder2, Inter_trialOrder2Age)

gender = glmer(correctAns ~ age + urnSize + LR + trialType + trialType*LR + trialOrd + sex + (1|sID_new), data=data, family="binomial", control = contr) # take out the interactions that were not a sign. improvement


anova(null, age, urnsize, LR, trialType, interTTLR, trialOrder, gender)

```

```{r, warning=FALSE}
#learning effect in 8 first trials to compare results with Gualtieri et al.
data_8trials <- subset(data, trialOrd <9)

# null model
null_8trials = glmer(correctAns ~ 1 + (1|sID_new), data=data_8trials, family="binomial", control = contr)

# add age
age_8trials = glmer(correctAns ~ age + (1|sID_new), data=data_8trials, family="binomial", control = contr)

# add urnsize
urnsize_8trials = glmer(correctAns ~ age + urnSize + (1|sID_new), data=data_8trials, family="binomial", control = contr)

# add LR 
LR_8trials = glmer(correctAns ~ age + urnSize + LR + (1|sID_new), data=data_8trials, family="binomial", control = contr)

# add trial type
trialType_8trials = glmer(correctAns ~ age + urnSize + LR + trialType + (1|sID_new), data=data_8trials, family="binomial", control = contr)


# add interaction between LR and trial type
interTTLR_8trials = glmer(correctAns ~ age + urnSize + LR + trialType + trialType*LR + (1|sID_new), data=data_8trials, family="binomial", control = contr)


#add trial order
trialOrder_8trials = glmer(correctAns ~ age + urnSize + LR + trialType + trialType*LR + trialOrd + (1|sID_new), data=data_8trials, family="binomial", control = contr)

anova(null_8trials, age_8trials, urnsize_8trials, LR_8trials, trialType_8trials, interTTLR_8trials, trialOrder_8trials)

```
no effect of learning in the first 8 trials



```{r, warning=FALSE}
##Learning effect

data$trialOrd <- as.numeric(data$trialOrd)

plotdata_trialOrd <- data %>%
  group_by(age, trialOrd, correctAns2) %>%
  summarize(n = n()) %>% 
  mutate(pct = n/sum(n),
         lbl = scales::percent(pct))

plotdata_trialOrd <- subset(plotdata_trialOrd, correctAns2 == "Correct")

theme_set(theme_light(base_size = 20, base_family = "Poppins"))


reg <- ggplot(plotdata_trialOrd, aes(x=trialOrd, y=pct)) +
  facet_grid(~ age) +
  geom_point(shape=1) +    
  geom_smooth(method=lm,   
              se=FALSE, level = .95, colour = "#045a8d")+
  myTheme+
  theme(axis.title=element_text(size=18),axis.text.x = element_text(size = 12),axis.text.y = element_text(size = 18),
        legend.position = "none")+
  labs(x = "Trial order", y = "Proportion of correct choices") +
  scale_y_continuous(breaks = c(0,0.25,0.50,0.75,1), limits = c(0,1))+
  scale_x_continuous(limits = c(0, 17), breaks=seq(1, 16, 1), expand = c(0,0))+
  stat_cor(method = "pearson", label.x = 1, label.y = 0.2, size = 5, digits = 2, alternative = "two.sided")

reg

ggsave("learning_effect.pdf",width=10,height=5)
ggsave("learning_effect.svg",width=10,height=5)

```

```{r, warning=FALSE}
# compare if correlation in 5-year-olds is higher from those in 4- and 3-year-olds
library(psych)
paired.r(0.77,0.57,n=16, twotailed=T) # 5 vs. 4 years
paired.r(0.72,0.16,n=16, twotailed=T) # 5 vs. 3 years
```

## Learning Effect due to trial-error-learning?

```{r}
# delete trials 1 to 4 
data_learning <- read_csv("Exp_data.csv")

data_learning <- subset(data_learning, trialNr > 4)

data_learning$trialOrd <- rep(1:18,90)

data_learning <- subset(data_learning, select = c(sID_new, age, correctAns, trialOrd, feedBack))





data_learning$age <- factor(data_learning$age, levels = c(3:5), labels = c("3 years", "4 years", "5 years"))

data_learning$feedBack <- factor(data_learning$feedBack, levels = c("positive","negative"))

data_learning$trialOrd <- as.numeric(data_learning$trialOrd)

data_learning$correctAns <- as.numeric(data_learning$correctAns)


```


```{r}
# Define lag (number of trials to look back for feedback)
lag <- 1  # Look back at feedback from the previous trial

data_grouped <- data_learning %>%
  group_by(sID_new, age) %>%  # Group by individual ID
  mutate(correctAns.prev = lag(correctAns, 1),
         feedBack.prev = lag(feedBack, 1)) %>%  
  ungroup() 


data_grouped <- subset(data_grouped, trialOrd != 1)

data_grouped <- data_grouped %>% 
  mutate(correctAns_prev = if_else(correctAns.prev == 1, "Correct", "Wrong"))



data_grouped_5years <- subset(data_grouped, age == "5 years")
data_grouped_4years <- subset(data_grouped, age == "4 years")
data_grouped_3years <- subset(data_grouped, age == "3 years")

```

```{r}
 

# Logistic regression model (correctness as outcome)
#glm_0 <- glm(correctAns ~ 1 + (1|sID_new), data = data_grouped, family = binomial)

glm_learning_5years <- glm(correctAns ~ feedBack.prev * correctAns_prev + (1|sID_new), data = data_grouped_5years, family = binomial)
glm_learning_4years <- glm(correctAns ~ feedBack.prev * correctAns_prev + (1|sID_new), data = data_grouped_4years, family = binomial)
glm_learning_3years <- glm(correctAns ~ feedBack.prev * correctAns_prev + (1|sID_new), data = data_grouped_3years, family = binomial)
# Summary of the model
summary(glm_learning_5years)
summary(glm_learning_4years)
summary(glm_learning_3years)

#anova(glm_0,glm_learning)
```
No significant predictors. Hence, we don't expect there to be a learning effect that is based on trial-error learning. 


## Individual-level analyses


```{r message=FALSE, include=FALSE}

data <- read_csv("Exp_data.csv")
```


```{r}
# delete trials 1 to 4 

data <- subset(data, trialNr > 4)

data <- subset(data, select = c(sID_new,age,sex,trialNr,selectedUrn,correctAns))

#new column with correct answer text
data <- data %>% 
  mutate(
         given_answer = if_else(selectedUrn == 2, 0, 1))

data <- subset(data, trialNr > 6) 
```


```{r}
data$age <- factor(data$age, levels = c(3:5), labels = c("3 years", "4 years", "5 years"))
data$trialNr <- factor(data$trialNr, levels = c(7:22)) 

```



```{r}
##new data frame with strategies
#create new data frame with the urn that should be chose according to different strategies
sID_new <- rep(1:90, each = 16) 
trialNr <- rep(7:22,90) 
strat_ratio <- rep(1:0,720) 
strat_more_good <- rep(c(1,0,1,0,0,1,1,0,1,0,1,0,0,1,1,0),90)
strat_less_bad <- rep(c(1,0,1,0,1,0,0,1,1,0,1,0,1,0,0,1),90) 
strat_random <- rep(0.5,1440) 
strat_good_urn <- rep(c(1,0,0.5,0.5,0.5,0.5,0.5,0.5,1,0,0.5,0.5,0.5,0.5,0.5,0.5),90) 

strategies <- data.frame(sID_new, trialNr, strat_random, strat_ratio, strat_more_good, strat_less_bad, strat_good_urn)

strategies$trialNr <- factor(strategies$trialNr, levels = c(7:22)) 

```

```{r}
# merge both data frames
data_strategies<- merge(data, strategies, by = c("sID_new", "trialNr"))
data_strategies <- left_join(data, strategies)

data_strategies <- subset(data_strategies, select = c(sID_new,age,sex,trialNr,given_answer,strat_random,strat_ratio,strat_more_good,strat_less_bad,strat_good_urn))
```

```{r}
# calculate absolute difference between given answer and different strategies

data_strategies <- data_strategies %>%
  mutate(diff_strat_random = abs(given_answer - strat_random))
data_strategies <- data_strategies %>%
  mutate(diff_strat_ratio = abs(given_answer - strat_ratio))
data_strategies <- data_strategies %>%
  mutate(diff_strat_more_good = abs(given_answer - strat_more_good))
data_strategies <- data_strategies %>%
  mutate(diff_strat_less_bad = abs(given_answer - strat_less_bad))
data_strategies <- data_strategies %>%
  mutate(diff_strat_good_urn = abs(given_answer - strat_good_urn))
```

```{r}
#summarize data to get the sum of new columns
data_strategies_sum <- data_strategies %>%
  group_by(sID_new,age) %>%
  summarise(diff_strat_random = sum(diff_strat_random),
            diff_strat_ratio = sum(diff_strat_ratio),
            diff_strat_more_good = sum(diff_strat_more_good),
            diff_strat_less_bad = sum(diff_strat_less_bad),
            diff_strat_good_urn = sum(diff_strat_good_urn))

data_strategies_final <- data_strategies_sum %>%
   mutate(random = (16 - diff_strat_random) / 16,
          ratio = (16 - diff_strat_ratio) / 16,
          more_good = (16 - diff_strat_more_good) / 16,
          less_bad = (16 - diff_strat_less_bad) / 16,
          good_urn = (16 - diff_strat_good_urn) / 16)

data_strategies_final <- subset(data_strategies_final, select = c(sID_new,age,random,ratio,more_good,less_bad,good_urn))

```

```{r}
#final step, qualitative choice
data_strategies_final <- data_strategies_final %>%
   mutate(random_qual = ifelse(random == max(c(random, ratio, more_good, less_bad, good_urn)),1, 0),
          ratio_qual = ifelse(ratio == max(c(random, ratio, more_good, less_bad, good_urn)),1, 0),
          more_good_qual = ifelse(more_good == max(c(random, ratio, more_good, less_bad, good_urn)),1,0),
          less_bad_qual = ifelse(less_bad == max(c(random, ratio, more_good, less_bad, good_urn)),1,0),
          good_urn_qual = ifelse(good_urn == max(c(random, ratio, more_good, less_bad, good_urn)),1,0))

#replace by 0 all 1 if random =1
data_strategies_final <- data_strategies_final %>%
  mutate(
    across(c(ratio_qual, more_good_qual, less_bad_qual, good_urn_qual), ~ ifelse(.x == 1 & random_qual == 1, 0, .x)),
    .after = random_qual # Ensure the new column comes after random_qual
  )

```

```{r}
#weigh 
#divide by the sum of 1s
data_strategies_final <- data_strategies_final %>%
  mutate(random_qual_weighed = random_qual / sum(random_qual,ratio_qual, more_good_qual, less_bad_qual, good_urn_qual),
         ratio_qual_weighed = ratio_qual / sum(random_qual,ratio_qual, more_good_qual, less_bad_qual, good_urn_qual),
         more_good_qual_weighed = more_good_qual / sum(random_qual,ratio_qual, more_good_qual, less_bad_qual, good_urn_qual),
         less_bad_qual_weighed = less_bad_qual / sum(random_qual,ratio_qual, more_good_qual, less_bad_qual, good_urn_qual),
         good_urn_qual_weighed = good_urn_qual / sum(random_qual,ratio_qual, more_good_qual, less_bad_qual, good_urn_qual))
```


```{r}
#final data frame with weighed values
plotdata_strategies <- subset(data_strategies_final, select = c(sID_new,age,random_qual_weighed,ratio_qual_weighed,more_good_qual_weighed,less_bad_qual_weighed,good_urn_qual_weighed))



plotdata_strategies <- rename(plotdata_strategies, random = random_qual_weighed)
plotdata_strategies <- rename(plotdata_strategies, ratio = ratio_qual_weighed)
plotdata_strategies <- rename(plotdata_strategies, more_good = more_good_qual_weighed)
plotdata_strategies <- rename(plotdata_strategies, less_bad = less_bad_qual_weighed)
plotdata_strategies <- rename(plotdata_strategies, good_urn = good_urn_qual_weighed)


#transform in long format
plotdata_strategies_long <- plotdata_strategies %>%
  gather(key = "strategy", value = "count", -sID_new, -age)


```

```{r}
#plot data
plotdata <- plotdata_strategies_long

age_strategy_summary <- plotdata %>%
  group_by(age, strategy) %>%
  summarize(average_count = mean(count)*100)

round(age_strategy_summary$average_count,3)
```



```{r}
g <- ggplot(age_strategy_summary, aes(x = age, y = average_count, fill = factor(strategy, levels = c("good_urn", "random", "less_bad", "more_good", "ratio")))) +
  geom_bar(stat = "identity") +
  labs( x = "Age", y = "Percentage") +  # Add labels and title (corrected "Percentage" to "Average Count")
  coord_cartesian(ylim = c(0, 100)) +  # Set y-axis limits from 0 to 100 
  theme_classic() +
  scale_fill_brewer(palette = "Pastel1", labels = c("Good urn", "Random", "Less bad", "More good", "Ratio"), name = "Strategy") 


g

ggsave("ind_strategies_bars.pdf",width=6,height=5)

ggsave("ind_strategies_bars.svg",width=6,height=5)

```