```{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, desktop_conf == "1: I confirm" & attent_conf == "1: I will take it seriously")
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$condition)
table(tdata$structure, tdata$test_statement_type)
```


Factorize:

```{r}
tdata$structure <- factor(tdata$structure, levels = c("irreversible", "reversible"), labels = c("irreversible", "reversible"))
tdata$test_statement_type <- factor(tdata$test_statement_type, levels = c("trigger_first", "trigger_second", "maintainer_first", "maintainer_second"),
                                    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_long <- tdata
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$statement_rating, factor = 0.01, amount = 0.004)

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


g <- ggplot(tdata_sub, aes(x = test_statement_type, 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 = statement_rating, group=1), fun.data = mean_cl_boot, 
               geom = "errorbar", width = 0, size = 1) +
  stat_summary(aes(y = statement_rating, group=1, color = structure), fun.y=mean, geom="line", 
               shape = 22, size = 1.5, alpha = .7)+
  stat_summary(aes(y = statement_rating, group=1, fill = structure), fun.y=mean, geom="point", 
               color = "black", shape = 22, size = 3, group=1, alpha = 1)+
  stat_summary(aes(y = statement_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 = statement_rating, y = test_statement_type, 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)
```


## Between-subejcts ANOVA 

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

a1 <- aov_car(statement_rating ~ structure*test_statement_type + Error(subj_code), 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$statement_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 = statement_rating, group=1), fun.data = mean_cl_boot, 
               geom = "errorbar", width = 0, size = 1) +
  stat_summary(aes(y = statement_rating, group=1, fill = structure), fun.y=mean, geom="point", 
               color = "black", shape = 22, size = 3, group=1, alpha = 1)+
  stat_summary(aes(y = statement_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 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$statement_rating, factor = 0.01, amount = 0.004)

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

g <- ggplot(tdata_sub, aes(x = test_statement_type, 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 = test_statement_type), alpha = 0.5, width = 0.15, height = 0.2) +
  stat_summary(aes(y = statement_rating, group=1), fun.data = mean_cl_boot, 
               geom = "errorbar", width = 0, size = 1) +
  stat_summary(aes(y = statement_rating, group=1, fill = test_statement_type), fun.y=mean, geom="point", 
               color = "black", shape = 22, size = 3, group=1, alpha = 1)+
  stat_summary(aes(y = statement_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", "test_statement_type")) 
ls2

contrasts <- emmeans(a1, ~ test_statement_type|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.

```{r}
library(dplyr)

tdata_long %>%
  group_by(structure, test_statement_type) %>%
  summarise_at(vars(statement_rating), list(name=sd))
```

```{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-two-ind-group-pooled

library(MOTE)


# irreversible condition

stats <- d.ind.t(
  m1 = 9.63,
  m2 = 0.30,
  sd1 = 1.6709141,
  sd2 = 0.6512587,
  n1 = 30,
  n2 = 30,
  a = 0.05
)

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

stats <- d.ind.t(
  m1 = 8.70,
  m2 = 3.67,
  sd1 = 2.5072309,
  sd2 = 3.2412570,
  n1 = 30,
  n2 = 30,
  a = 0.05
)

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, ~ test_statement_type*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
# reversible   second \nmaintained   8.47 0.529 232    7.424     9.51
# irreversible second \nmaintained   3.27 0.529 232    2.224     4.31


#SDs
# reversible	second \nmaintained	2.5961953	
# irreversible	second \nmaintained	4.1015837		


stats <- d.ind.t(
  m1 = 8.47,
  m2 = 3.27,
  sd1 = 2.5961953,
  sd2 = 4.1015837,
  n1 = 30,
  n2 = 30,
  a = 0.05
)

stats$estimate
```

