```{r setup, include=FALSE}
knitr::opts_chunk$set(warning = FALSE, message = FALSE) 
# packages
library(ez)
library(reshape2)
library(reshape)
library(ggplot2)
library(plyr)
library(pastecs)
library(ez)
library(data.table)
library(tidyverse) 

library(showtext)
library(readr)


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

Read in the data:

```{r}
tdata <- read_csv("exp_data.csv")
```

Remove participants who failed the timing check query:

```{r}
tdata <- subset(tdata, timing_check_correct == "correct")
```


# Subject demographics

```{r}
# demographics 

min(tdata$age)
max(tdata$age)
mean(tdata$age)
sd(tdata$age)

# 1 = male, 2 = female, 3 = other
table(tdata$gender)
```

1 = male, 2 = female, 3 = non-binary, 4 = prefer not to say


# Data preparation 

Check n in each condition:

```{r}
table(tdata$structure, tdata$test_queries)
```


Factorize:

```{r}
tdata$structure <- factor(tdata$structure, levels = c("irreversible", "reversible"), labels = c("irreversible", "reversible"))
tdata$test_queries <- factor(tdata$test_queries, levels = c("original", "novel"), labels = c("original", "novel"))
```


Make a subset containing only the columns relevant for analyses and turn into long format:

```{r}
# Subset: 
tdata_sub <- subset(tdata, select = c("subj_code", "condition", "structure", "test_queries", "first_caused_rating", "second_caused_rating",
                                      "first_maintained_rating", "second_maintained_rating"))

tdata_original <- subset(tdata_sub, test_queries == "original", select = c(1:6))
tdata_novel <- subset(tdata_sub, test_queries == "novel")

# into long format:
library(tidyr)
tdata_long_original <- gather(tdata_original, entity, rating, 5:6)
tdata_long_novel <- gather(tdata_novel, entity, rating, 5:8)

# factorize entity 
tdata_long_original$entity <- factor(tdata_long_original$entity, levels = c("first_caused_rating", "second_caused_rating"), 
                            labels = c("first \ncaused", "second \ncaused"))



tdata_long_novel$entity <- factor(tdata_long_novel$entity, levels = c("first_caused_rating", "second_caused_rating", 
                                                          "first_maintained_rating", "second_maintained_rating"), 
                            labels = c("first \ncaused", "second \ncaused","first \nmaintained", "second \nmaintained"))


tdata_long_novel$rating <- as.numeric(tdata_long_novel$rating)

```


# Replication part of study 

## Graphs and summaries

Replication part of the study 

```{r}
myTheme <- theme(plot.title = element_text(face="bold", size = 22),
        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"))


tdata_sub <- tdata_long_original


library(see)
## first, turn sID into a factor
tdata_sub$subj_code <- factor(tdata_sub$subj_code)

pd <- position_dodge(width = 0.3)

tdata_sub$valueJitter <- jitter(tdata_sub$rating, factor = 0.01, amount = 0.004)

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



g <- ggplot(tdata_sub, aes(x = entity, y = valueJitter)) +
  guides(fill=FALSE)+
  facet_grid( ~ structure, labeller = label_both)+
  scale_y_continuous(limits = c(-0.3, 8.3), breaks=seq(0, 8, 1), expand = c(0,0)) +
  geom_violinhalf(aes(y = rating, group = entity, fill = structure), color = NA, 
                position=position_dodge(1), alpha = 0.4)+
  geom_jitter(aes(color = structure), alpha = 0.5, width = 0.15, height = 0.2) +
  stat_summary(aes(y = rating, group=1), fun.data = mean_cl_boot, 
               geom = "errorbar", width = 0, size = 1) +
  stat_summary(aes(y = rating, group=1, color = structure), fun.y=mean, geom="line", 
               shape = 22, size = 1.5, alpha = .7)+
  stat_summary(aes(y = rating, group=1, fill = structure), fun.y=mean, geom="point", 
               color = "black", shape = 22, size = 3, group=1, alpha = 1)+
  stat_summary(aes(y = rating,group=1), fun.y=median, geom="point", color = "black", shape = 3, size = 4, 
               group=1, alpha = 1, position = position_dodge(width = 0.5))+
  labs(x = "Causal Statement", y = "Agreement rating") +
  scale_color_manual(name = "Strength",values=c("#66c2a5", "#e78ac3", "#8da0cb", "#a6d854"))+
  scale_fill_manual(name = "Strength",values=c("#66c2a5", "#e78ac3", "#8da0cb", "#a6d854"))+
  theme(legend.position = "none")+
  annotate("text", x = 0.5, y = 2, label = c("completely disagree"), angle = 90)+
  annotate("text", x = 0.5, y = 8.2, label = c("completely agree"), angle = 90)+
  myTheme+
  theme(panel.grid.major = element_line(color = "lightgrey",
                                        size = 0.5,
                                        linetype = 'dotted'))+
  stat_summary(aes(label=round(after_stat(y),2)), fun.y=mean, geom="text", size=5,
             vjust = -2)

g

#ggsave("results_means_mainDV_original.svg",width=8,height=5)
#ggsave("results_means_mainDV_original.pdf",width=8,height=5)

```




```{r}
library(ggridges)
g2 <- ggplot(tdata_sub, aes(x = rating, y = entity, fill = entity)) +
  facet_grid( ~ structure,labeller = label_both)+
  scale_x_continuous(breaks = seq(0, 10, 1))+
  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)+
  scale_fill_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 = "Agreement rating", y = "Causal Statement") +
  scale_y_discrete(limits=rev)+
  myTheme+
  theme_light(base_family = "Poppins", base_size = 20)+
  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))

g2

#ggsave("results_dist_original.svg",width=8,height=5)
#ggsave("results_dist_original.pdf",width=8,height=5)
```


## Mixed ANOVA 

```{r}
library(afex)
library(emmeans)

a1 <- aov_car(rating ~ structure*entity + Error(subj_code/entity), tdata_long_original, 
              anova_table = list(es = "pes"))
a1

```

No significant interaction between structure and entity (causal statement), but a clear effect of entity (causal statement). This replicates the original study.


## Contrasts 

```{r}

library(lsmeans)
# means

ls2 <- lsmeans(a1, c("structure", "entity")) 
ls2

contrasts <- emmeans(a1, ~ entity|structure)
s <- pairs(contrasts, adjust = "none") 


s
confint(s, level = 0.95)
```
First cause receives higher ratings in both reversibility conditions. 

Get standardized effect sizes for the contrasts

First get the SDs that need to be pooled (and the correlation between measures for the within-tests)

```{r}
library(dplyr)

tdata_long_original %>%
  group_by(structure, entity) %>%
  summarise_at(vars(rating), list(name=sd))
```


```{r}
tdata_original %>%
  group_by(structure, test_queries) %>%
  summarize(cor=cor(first_caused_rating, second_caused_rating))
```


```{r}
# using the functions from the MOTE package (see https://matthewbjane.quarto.pub/guide-to-effect-sizes-and-confidence-intervals/Standardized-Mean-Differences.html#sec-repeated-measures-drm)

library(MOTE)

# irreversible condition
stats <- d.dep.t.rm(
  m1 = 6.80,
  m2 = 1.36,
  sd1 = 2.254625,
  sd2 = 2.23383,
  n = 25,
  a = 0.05,
  r = -0.5063090	
)

stats$estimate

```

```{r}
# reversible condition

stats <- d.dep.t.rm(
  m1 = 6.92,
  m2 = 3.24,
  sd1 = 2.379075,
  sd2 = 3.573980,
  n = 25,
  a = 0.05,
  r = -0.3798763
)

stats$estimate
```


# Novel part of the study 

```{r}

tdata_sub <- tdata_long_novel


library(see)
## first, turn sID into a factor
tdata_sub$subj_code <- factor(tdata_sub$subj_code)

pd <- position_dodge(width = 0.3)

tdata_sub$valueJitter <- jitter(tdata_sub$rating, factor = 0.01, amount = 0.004)

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



g <- ggplot(tdata_sub, aes(x = entity, y = valueJitter)) +
  guides(fill=FALSE)+
  facet_grid( ~ structure, labeller = label_both)+
  scale_y_continuous(limits = c(-0.3, 8.3), breaks=seq(0, 8, 1), expand = c(0,0)) +
  geom_jitter(aes(color = structure), alpha = 0.5, width = 0.15, height = 0.2) +
  stat_summary(aes(y = rating, group=1), fun.data = mean_cl_boot, 
               geom = "errorbar", width = 0, size = 1) +
  stat_summary(aes(y = rating, group=1, color = structure), fun.y=mean, geom="line", 
               shape = 22, size = 1.5, alpha = .7)+
  stat_summary(aes(y = rating, group=1, fill = structure), fun.y=mean, geom="point", 
               color = "black", shape = 22, size = 3, group=1, alpha = 1)+
  stat_summary(aes(y = rating,group=1), fun.y=median, geom="point", color = "black", shape = 3, size = 4, 
               group=1, alpha = 1, position = position_dodge(width = 0.5))+
  labs(x = "Causal Statement", y = "Agreement Rating") +
  scale_color_manual(name = "Strength",values=c("#66c2a5", "#e78ac3", "#8da0cb", "#a6d854"))+
  scale_fill_manual(name = "Strength",values=c("#66c2a5", "#e78ac3", "#8da0cb", "#a6d854"))+
  theme(legend.position = "none")+
  annotate("text", x = 0.5, y = 1.7, label = c("completely disagree"), angle = 90)+
  annotate("text", x = 0.5, y = 6.7, label = c("completely agree"), angle = 90)+
  myTheme+
  theme(panel.grid.major = element_line(color = "lightgrey",
                                        size = 0.5,
                                        linetype = 'dotted'))+
  stat_summary(aes(label=round(after_stat(y),2)), fun.y=mean, geom="text", size=5,
             vjust = -6)

g

#ggsave("results_means_mainDV_novel.svg",width=12,height=5)
#ggsave("results_means_mainDV_novel.pdf",width=12,height=5)

```



```{r}
library(ggridges)
g2 <- ggplot(tdata_sub, aes(x = rating, y = entity, fill = structure)) +
  facet_grid( ~ structure,labeller = label_both)+
  scale_x_continuous(breaks = seq(0, 10, 1))+
  geom_density_ridges(alpha = 0.5)+
  scale_fill_manual(values=c("#66c2a5", "#e78ac3", "#8da0cb", "#a6d854"))+
  labs(x = "Agreement Rating", y = "Causal Statement") +
  scale_y_discrete(limits=rev)+
  myTheme+
  theme_light(base_family = "Poppins", base_size = 20)+
  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", y = 0.7, x = 0, label = c("completely disagree"), angle = 0)+
  annotate("text", y = 0.7, x = 9, label = c("completely agree"), angle = 0)

g2

#ggsave("results_dist_novel.svg",width=12,height=5)
#ggsave("results_dist_novel.pdf",width=12,height=5)
```


## Mixed ANOVA 

```{r}
library(afex)
library(emmeans)

a2 <- aov_car(rating ~ structure*entity + Error(subj_code/entity), tdata_long_novel, anova_table = list(es = "pes"))
a2

```
Predicted interaction effect was obtained, but also visualize the main effects: 

Plot main effect of structure: 

```{r}

tdata_sub <- tdata_long_novel


library(see)
## first, turn sID into a factor
tdata_sub$subj_code <- factor(tdata_sub$subj_code)

pd <- position_dodge(width = 0.3)

tdata_sub$valueJitter <- jitter(tdata_sub$rating, factor = 0.01, amount = 0.004)

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


g <- ggplot(tdata_sub, aes(x = structure, y = valueJitter)) +
  guides(fill=FALSE)+
  scale_y_continuous(limits = c(-0.3, 8.3), breaks=seq(0, 8, 1), expand = c(0,0)) +
  geom_jitter(aes(color = structure), alpha = 0.5, width = 0.15, height = 0.2) +
  stat_summary(aes(y = rating, group=1), fun.data = mean_cl_boot, 
               geom = "errorbar", width = 0, size = 1) +
  stat_summary(aes(y = rating, group=1, fill = structure), fun.y=mean, geom="point", 
               color = "black", shape = 22, size = 3, group=1, alpha = 1)+
  stat_summary(aes(y = rating,group=1), fun.y=median, geom="point", color = "black", shape = 3, size = 4, 
               group=1, alpha = 1, position = position_dodge(width = 0.5))+
  labs(x = "Causal Statement", y = "Agreement Rating") +
  scale_color_manual(name = "Strength",values=c("#66c2a5", "#e78ac3", "#8da0cb", "#a6d854"))+
  scale_fill_manual(name = "Strength",values=c("#66c2a5", "#e78ac3", "#8da0cb", "#a6d854"))+
  theme(legend.position = "none")+
  annotate("text", x = 0.5, y = 1.7, label = c("completely disagree"), angle = 90)+
  annotate("text", x = 0.5, y = 6.7, label = c("completely agree"), angle = 90)+
  myTheme+
  theme(panel.grid.major = element_line(color = "lightgrey",
                                        size = 0.5,
                                        linetype = 'dotted'))+
  stat_summary(aes(label=round(after_stat(y),2)), fun.y=mean, geom="text", size=5,
             vjust = -6)

g

```

Plot main effect of entity (causal statement):

```{r}

tdata_sub <- tdata_long_novel


library(see)
## first, turn sID into a factor
tdata_sub$subj_code <- factor(tdata_sub$subj_code)

pd <- position_dodge(width = 0.3)

tdata_sub$valueJitter <- jitter(tdata_sub$rating, factor = 0.01, amount = 0.004)

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


g <- ggplot(tdata_sub, aes(x = entity, y = valueJitter)) +
  guides(fill=FALSE)+
  scale_y_continuous(limits = c(-0.3, 8.3), breaks=seq(0, 8, 1), expand = c(0,0)) +
  geom_jitter(aes(color = entity), alpha = 0.5, width = 0.15, height = 0.2) +
  stat_summary(aes(y = rating, group=1), fun.data = mean_cl_boot, 
               geom = "errorbar", width = 0, size = 1) +
  stat_summary(aes(y = rating, group=1, fill = entity), fun.y=mean, geom="point", 
               color = "black", shape = 22, size = 3, group=1, alpha = 1)+
  stat_summary(aes(y = rating,group=1), fun.y=median, geom="point", color = "black", shape = 3, size = 4, 
               group=1, alpha = 1, position = position_dodge(width = 0.5))+
  labs(x = "Causal Statement", y = "Agreement Rating") +
  scale_color_manual(name = "Strength",values=c("#66c2a5", "#e78ac3", "#8da0cb", "#a6d854"))+
  scale_fill_manual(name = "Strength",values=c("#66c2a5", "#e78ac3", "#8da0cb", "#a6d854"))+
  theme(legend.position = "none")+
  annotate("text", x = 0.5, y = 1.7, label = c("completely disagree"), angle = 90)+
  annotate("text", x = 0.5, y = 6.7, label = c("completely agree"), angle = 90)+
  myTheme+
  theme(panel.grid.major = element_line(color = "lightgrey",
                                        size = 0.5,
                                        linetype = 'dotted'))+
  stat_summary(aes(label=round(after_stat(y),2)), fun.y=mean, geom="text", size=5,
             vjust = -6)

g
```


## Contrasts 

```{r}

library(lsmeans)
# means

ls2 <- lsmeans(a2, c("structure", "entity")) 
ls2

contrasts <- emmeans(a2, ~ entity|structure)
s <- pairs(contrasts, adjust = "none")


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

The planned contrasts from this list are "first caused - second caused" in each reversibility condition. They are both significant. 


Get effect sizes for the relevant contrasts

First get the SDs that need to be pooled (and the correlation between the measures included in within-comparisons)

```{r}
library(dplyr)

tdata_long_novel %>%
  group_by(structure, entity) %>%
  summarise_at(vars(rating), list(name=sd))
```

```{r}
library(dplyr)

tdata_novel %>%
  group_by(structure, test_queries) %>%
  summarize(cor=cor(first_caused_rating, second_caused_rating))
```





```{r}
# irreversible condition
stats <- d.dep.t.rm(
  m1 = 7.70,
  m2 = 0.76,
  sd1 = 1.182353,
  sd2 = 2.036002,
  n = 50,
  a = 0.05,
  r = -0.1492076
)

stats$estimate
```

```{r}
# reversible condition
stats <- d.dep.t.rm(
  m1 = 7.44,
  m2 = 1.94,
  sd1 = 0.213 * sqrt(50),
  sd2 = 0.380 * sqrt(50),
  n = 50,
  a = 0.05,
  r = -0.4834880
)

stats$estimate
```

In the paper I also say that maintainer ratings for the first cause are higher in the reversible than in the irreversible condition. An effect size for this is given by: 

```{r}
# structure    entity             lsmean    SE df lower.CL upper.CL
# irreversible first..maintained    2.74 0.464 98  1.81965     3.66
# reversible   first..maintained    5.44 0.464 98  4.51965     6.36

stats <- d.ind.t(
  m1 = 5.44,
  m2 = 2.74,
  sd1 = 1.774651,
  sd2 = 1.182353,
  n1 = 50,
  n2 = 50,
  a = 0.05
)

stats$estimate
 
```



The other planned contrast is the one comparing the maintainer ratings between the two reversibility conditions. The following code will produce a long list of comparisons, from which the relevant one will be "second..maintained irreversible - second..maintained reversible"


```{r}
contrasts <- emmeans(a2, ~ entity*structure)
s <- pairs(contrasts, adjust = "none")


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

The relevant contrast "second..maintained irreversible - second..maintained reversible" is significant

Get a d value for this one too: 


```{r}
# structure    entity             lsmean    SE df lower.CL upper.CL
# irreversible second..maintained   2.24 0.471 98  1.30472     3.18
# reversible   second..maintained   5.12 0.471 98  4.18472     6.06


#SDs
#irreversible	first \ncaused	1.182353		
#irreversible	second \ncaused	2.036002		
#irreversible	first \nmaintained	3.427440		
#irreversible	second \nmaintained	3.365855		
#reversible	first \ncaused	1.774651		
#reversible	second \ncaused	3.203378		
#reversible	first \nmaintained	3.124361		
#reversible	second \nmaintained	3.298979


stats <- d.ind.t(
  m1 = 5.12,
  m2 = 2.24,
  sd1 = 3.298979,
  sd2 = 3.365855,
  n1 = 50,
  n2 = 50,
  a = 0.05
)

stats$estimate
```
