```{r setup, include=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")
```

Exclude subjects who failed at least one control questions 

```{r}
tdata <- subset(tdata, intro_check == "2: A green ball lies on a brown wooden board. Then a yellow ball falls down onto the wooden board. The green ball then shoots into the air.")

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}
red_first <- subset(tdata, first_object == "red")
colnames(red_first)[10] ="first_caused_rating"
colnames(red_first)[11] ="second_caused_rating"
colnames(red_first)[12] = "first_maintained_rating"
colnames(red_first)[13] = "second_maintained_rating"

yellow_first <- subset(tdata, first_object == "yellow")
colnames(yellow_first)[10] ="second_caused_rating"
colnames(yellow_first)[11] ="first_caused_rating"
colnames(yellow_first)[12] = "second_maintained_rating"
colnames(yellow_first)[13] = "first_maintained_rating"


tdata <- rbind(red_first, yellow_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_caused_rating", "second_caused_rating", "first_maintained_rating", "second_maintained_rating"))

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


# factorize entity 
tdata_long$statement <- factor(tdata_long$statement, levels = c("first_caused_rating", "second_caused_rating", "first_maintained_rating", "second_maintained_rating")
                            , labels = c("first \ncaused", "second \ncaused", "first \nmaintained", "second \nmaintained"))
```


```{r}
# load stringr library
library(stringr)
 
# Split name column into firstname and last name
tdata_long[c('cause', 'role')] <- str_split_fixed(tdata_long$statement, ' ', 2)


tdata_long$role<- factor(tdata_long$role, levels = c("\ncaused", "\nmaintained")
                            , labels = c("caused", "maintained"))


tdata_long$cause<- factor(tdata_long$cause, levels = c("first", "second")
                            , labels = c("first", "second"))
```



# 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 = statement, 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"))+
  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)+
  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 = -4)

g

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

g

```




```{r}
library(ggridges)
g2 <- ggplot(tdata_long, aes(x = rating, y = statement, 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.5, x = 1, label = c("completely disagree"), angle = 0)+
  annotate("text", y = 0.5, x = 9, 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*statement + Error(subj_code/statement), tdata_long, anova_table = list(es = "pes"))
a1

```

Predicted interaction effect is significant, but other main effects should also be visualized.

Visualize Main effect of structure (causal 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 Structure", 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"))+
  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)+
  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 = -4)

g


```

Main effect of 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 = statement, 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 = statement), 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 = statement), 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"))+
  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)+
  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 = -4)

g


```

## Contrasts 

```{r}

library(lsmeans)
# means

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

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


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

Planned contrasts in the list above are "first caused - second caused" in each reversibility condition. 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, statement) %>%
  summarise_at(vars(rating), list(name=sd))
```

```{r}
library(dplyr)

tdata %>%
  group_by(structure) %>%
  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

#structure    statement          lsmean    SE df lower.CL upper.CL
# irreversible first..caused        9.62 0.131 98    9.361     9.88
# irreversible second..caused       1.16 0.420 98    0.327     1.99



stats <- d.dep.t.rm(
  m1 = 9.62,
  m2 = 1.16,
  sd1 = 1.1228608,
  sd2 = 2.4773257,
  n = 50,
  a = 0.05,
  r = -0.3371902
)

stats$estimate
```


```{r}

# reversible condition

#structure    statement          lsmean    SE df lower.CL upper.CL
# reversible   first..caused        9.80 0.131 98    9.541    10.06
# reversible   second..caused       1.86 0.420 98    1.027     2.69



stats <- d.dep.t.rm(
  m1 = 9.80,
  m2 = 1.86,
  sd1 = 0.6700594,
  sd2 = 3.3867087,
  n = 50,
  a = 0.05,
  r = -0.2194333	
)

stats$estimate
```


Now contrast testing the maintainer ratings in each reversibility condition. The code below will produce a long list from which the correct contrast must be looked up.

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


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

The planned contrast from this list is the one comparing the maintainer ratings between the reversibility conditions: "second..maintained irreversible - second..maintained reversible", which is significant.


Get effect size 

```{r}

#structure    statement          lsmean    SE df lower.CL upper.CL
#reversible   second..maintained   7.94 0.485 98    6.978     8.90
#irreversible second..maintained   3.60 0.485 98    2.638     4.56

#SDs
#reversible	second \nmaintained	3.2602116
#irreversible	second \nmaintained	3.5856858	

stats <- d.ind.t(
  m1 = 7.94,
  m2 = 3.60,
  sd1 = 3.2602116,
  sd2 = 3.5856858,
  n1 = 50,
  n2 = 50,
  a = 0.05
)

stats$estimate
```
