```{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()
```


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


```{r}
# Exclude subjects who failed the control questions 
tdata <- subset(tdata, intro_check == "2: A little squid swam from the right to the left.")
tdata <- subset(tdata, timing_check_correct == "correct")
tdata <- subset(tdata, learning_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 preparations 

Check n in each condition:

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


Factorize:

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

Relabel rating columns for red and yellow entity according to "first object" (so that they are called first_obj_rating and second_obj_rating):

```{r}
left_first <- subset(tdata, first_object == "left")
colnames(left_first)[10] ="first_obj_caused_rating"
colnames(left_first)[11] ="second_obj_caused_rating"
colnames(left_first)[12] ="first_obj_maintained_rating"
colnames(left_first)[13] ="second_obj_maintained_rating"

right_first <- subset(tdata, first_object == "right")
colnames(right_first)[11] ="first_obj_caused_rating"
colnames(right_first)[10] ="second_obj_caused_rating"
colnames(right_first)[13] ="first_obj_maintained_rating"
colnames(right_first)[12] ="second_obj_maintained_rating"

right_first <- right_first[c("subj_code", "condition", "desktop_conf", "attent_conf", "intro_check", 
                             "structure", "first_object", "learn_check_left", "learn_check_right",
                             "first_obj_caused_rating", "second_obj_caused_rating", "first_obj_maintained_rating", "second_obj_maintained_rating",
                             "timing_check", "recorded_at", 
                             "age", "gender", "tech_issues", "timing_check_correct", "learning_check_left_correct", 
                             "learning_check_right_correct", "learning_check_correct")]


names(left_first) %in% names(right_first)

tdata <- rbind(left_first, right_first)
```

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","structure","first_obj_caused_rating", "second_obj_caused_rating",
                                      "first_obj_maintained_rating", "second_obj_maintained_rating"))

# into long format:
library(tidyr)
tdata_long <- gather(tdata_sub, entity, rating, 3:6)


# factorize entity 
tdata_long$entity <- factor(tdata_long$entity, levels = c("first_obj_caused_rating", "second_obj_caused_rating", 
                                                          "first_obj_maintained_rating", "second_obj_maintained_rating"), 
                            labels = c("first \ncaused", "second \ncaused","first \nmaintained", "second \nmaintained"))
```



# Graphs and summaries


```{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


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, 10.3), breaks=seq(0, 10, 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")+
  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)+
  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)

g

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

```




```{r}
library(ggridges)
g2 <- ggplot(tdata_long, 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)+
   #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))+
  annotate("text", y = 0.7, x = 0, label = c("completely disagree"), angle = 0)+
  annotate("text", y = 0.7, x = 10, label = c("completely agree"), angle = 0)

g2

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

## Mixed ANOVA 

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

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

```
Predicted interaction effect was significant-


Visualize also the obtained main effect of "structure" (reversibility)

```{r}
tdata_sub <- tdata_long


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, 10.3), breaks=seq(0, 10, 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")+
  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)+
  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)

g

```

Plot main effect of entity (causal statement):

```{r}

tdata_sub <- tdata_long


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, 10.3), breaks=seq(0, 10, 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(a1, c("structure", "entity")) 
ls2

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


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

The planned contrasts from the list above are " first..caused - second..caused" in both "irreversible" and "reversible". Both are significant.

Get effect standardized effect sizes for the contrasts

First get the SDs that need to be pooled (and also the correlations between measures for the within-comparisons)

```{r}
library(dplyr)

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


```{r}
library(dplyr)

tdata %>%
  group_by(structure) %>%
  summarize(cor=cor(first_obj_caused_rating, second_obj_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

#structure    entity             lsmean    SE df lower.CL upper.CL
# irreversible first..caused        9.64 0.178 98    9.286     9.99
# irreversible second..caused       0.92 0.408 98    0.110     1.73

#Cor
#irreversible	-0.2529154			
#reversible	-0.3145000	

stats <- d.dep.t.rm(
  m1 = 9.64,
  m2 = 0.92,
  sd1 = 0.875051,
  sd2 = 1.700420,
  n = 50,
  a = 0.05,
  r = -0.2529154
)

stats$estimate
```

```{r}
# reversible condition 

#structure    entity             lsmean    SE df lower.CL upper.CL
# reversible   first..caused        9.40 0.178 98    9.046     9.75
# reversible   second..caused       2.74 0.408 98    1.930     3.55

#Cor
#irreversible	-0.2529154			
#reversible	-0.3145000	


stats <- d.dep.t.rm(
  m1 = 9.40,
  m2 = 2.74,
  sd1 = 1.551826,
  sd2 = 3.713242,
  n = 50,
  a = 0.05,
  r = -0.3145000
)

stats$estimate
```



The final planned contrast compares the maintainer ratings for the second cause in "irreversible" and "reversible". It must be looked up from the long list of contrasts below: 

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


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

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

Get effect size: 

```{r}
#structure    entity             lsmean    SE df lower.CL upper.CL
#irreversible second..maintained   3.32 0.409 98    2.508     4.13
#reversible   second..maintained   9.36 0.409 98    8.548    10.17

#SDs
#irreversible	second \nmaintained	3.661381
#reversible	second \nmaintained	1.826785	

stats <- d.ind.t(
  m1 = 9.36,
  m2 = 3.32,
  sd1 = 1.826785,
  sd2 = 3.661381,
  n1 = 50,
  n2 = 50,
  a = 0.05
)

stats$estimate
```

