```{r setup, include=FALSE}
# packages
library(ez)
library(reshape2)
library(reshape)
library(ggplot2)
library(plyr)
library(pastecs)
library(ez)
library(data.table)
#library(Rcmdr)
#Commander()

library(showtext)

font_add_google("Poppins", "Poppins")
font_add_google("Roboto Mono", "Roboto Mono")
showtext_auto()
##########################################Data Managment############################
# read in the data file 

#tdata.wide <- read.delim("dummy_data.txt", header=TRUE, sep="\t", na.strings="NA", dec=".", strip.white=TRUE)
tdata <- read.delim("Data_Exp_ProbQuery_Simul_vs_Single_CausePresentation.txt", header=TRUE, sep="\t", na.strings="NA", dec=".", strip.white=TRUE)

```



```{r}
# demographics 
min(tdata$Age)
max(tdata$Age)
mean(tdata$Age)
sd(tdata$Age)
```

```{r}
# 1 = male, 2 = female, 3 = other
count(tdata$Sex)
```


```{r, echo = FALSE}
# reorder factor 
tdata$learning_format <- factor(tdata$learning_format, levels = c("joint", "separate"), 
                          labels = c("joint presentation", "separate presentation"))

# reorder factor 
tdata$target_cause <- factor(tdata$target_cause, levels = c("single", "multi"), 
                            labels = c("single-effect cause", "multiple-effects cause"))

# reorder factor 
tdata$target_effect <- factor(tdata$target_effect, levels = c("first", "second", "third"), 
                                     labels = c("first", "second", "third"))


# to create chart, the data must be in long format and only contain the relevant dependent variables

# 1. make a subset with only the relevant dvs 

tdata_sub <- subset(tdata, select = 1:6)


# recode dependent variables 

tdata_sub$strength_rating <- (tdata_sub$strength_rating - 1) * 0.1
```


# Graphs

## Interaction between number of effects and similarity of effects

The bars represent means and the error bars depict 95% bootstrapped CIs.

```{r, echo = FALSE}

myTheme <- theme(plot.title = element_text(face="bold", size = 22),
        axis.title.x = element_text(face = "bold", size = 20),
        axis.title.y = element_text(face = "bold", size = 20),
        axis.text.x = element_text(size = 18, angle = 0), 
        axis.text.y = element_text(size = 14, 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"))


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

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

pd <- position_dodge(width = 0.3)

tdata_sub$valueJitter <- jitter(tdata_sub$strength_rating, factor = 1, amount = 0.04)

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

# new labes for the facets 
effects.labs <- c("Knowledge of: both causes", "Knowledge of: only one cause")
names(effects.labs) <- c("joint presentation", "separate presentation")

g <- ggplot(tdata_sub, aes(x=target_cause, y=valueJitter, color = target_cause, fill = target_cause)) +
  guides(fill=FALSE)+
  facet_grid( ~ learning_format, labeller = labeller(learning_format =effects.labs))+
  #ggtitle("Subjects' causal srength ratings") +
  scale_y_continuous(limits = c(-0.05, 1.05), breaks=seq(0, 1, 0.1), expand = c(0,0)) +
  scale_x_discrete(labels=c("single-effect \n cause", "common \n cause")) +
  geom_violinhalf(position=position_dodge(1), alpha = 0.2, color = NA)+
  geom_point(position=position_jitterdodge(dodge.width=0.5), alpha = 0.2) +
  stat_summary(fun.y=mean, colour="black", geom="line", group = 1, size = 1.5, linetype = "solid", alpha = 1)+
  stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = 0, size = 1, position = position_dodge(width = 0.5), color = "black") +
  stat_summary(fun.y=mean, geom="point", color = "black", shape = 22, size = 4, group=1, alpha = 1, position = position_dodge(width = 0.5))+
  stat_summary(fun.y=median, geom="point", color = "black", shape = 3, size = 4, group=1, alpha = 1, position = position_dodge(width = 0.5))+
  labs(x = "Target Cause", y = "Causal Strength Rating") +
  scale_color_manual(name = "Entity",values=c("#fc9272","#3182bd"))+
  scale_fill_manual(name = "Entity",values=c("#fc9272","#3182bd"))+
  theme(legend.position = "none")+
  myTheme
g


#ggsave("results_points.svg",width=11.5,height=6)
#ggsave("results_point.pdf",width=11.5,height=6)
```



# Descriptive Stats

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

################################################################################################################
##################################### Statistical Analyses #####################################################
################################################################################################################
library(pastecs)
library(lme4)
library(nlme)
library(ez)


##################################### Create descriptive stats ###################################

by(tdata_sub$strength_rating, list(tdata_sub$target_cause, tdata_sub$learning_format), stat.desc , basic = FALSE)

##################################### Model ######################################################
```

# Test


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

a1 <- aov_car(strength_rating ~ learning_format*target_cause + Error(sID), tdata_sub)
a1
```



```{r}
############### 
# a conditional analysis 

# Means for the different conditions

ls2 <- lsmeans(a1, c("target_cause"), by = "learning_format") # group means by between-condition
ls2

# Tets

# simple main effects 
pairs(ls2) # compares rep-measure differences separately for each between-factor level

# interaction contrast 
pairs(pairs(ls2), by = NULL)

#test(pairs(pairs(ls2), by = NULL), joint = TRUE) # This reproduces the F-Value of the ANOVA interaction


#lsmip(a1, High_Strength_Component ~ variable) # lsemans can also produce graphs
```


```{r}
# compute the confidence interval for the singular causation differences in each between-subject condition

# joint presentation

t <- qt(0.975, 716, lower.tail = TRUE, log.p = FALSE)
#t

effect <- "Mdiff"
Mdiff <- 0.1700
SE <- 0.0243
CI <- SE*t
CI_low <- Mdiff - CI
CI_up <- Mdiff + CI

Mdiff
CI_low
CI_up



# Plot 

myTheme <- theme(plot.title = element_text(face="bold", size = 22),
        axis.title.x = element_text(face = "bold", size = 20),
        axis.title.y = element_text(face = "bold", size = 20),
        axis.text.x = element_text(size = 18, angle = 0), 
        axis.text.y = element_text(size = 25, 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"))

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



barchart <- ggplot()+
  myTheme+
  #guides(fill=FALSE)+
  #facet_wrap(~Latency + SampleSize, ncol=2)+
  #ggtitle("Mean difference (95% CI)") +
  #coord_cartesian(ylim=c(-0.1,2)) + 
  scale_y_continuous(limits = c(-0.1, 0.5), breaks=seq(-0.1, 0.5, 0.1), expand = c(0,0)) +
  scale_x_discrete(labels=c("r")) +
  #annotate("rect", xmin=1.7, xmax=2.3, ymin=0.95, ymax=1.05, color="#31a354", fill = "white", size = 1) +
  #stat_summary(fun.y=mean, colour="grey20", geom="point", shape = 21, size = 3)+
  #stat_summary(fun.y = mean, geom = "bar", position = "dodge", colour = "black")+
  #stat_summary(fun.data = mean_cl_boot, geom = "errorbar", position = position_dodge(width = 0.90), width = 0.2) +
  #geom_jitter(width = 0.3, height = 0.02, alpha = 0.6, colour = "red") +
  #ggtitle("Means (95% bootstr. CIs)") +
  #theme(axis.text.x = element_text(size = 10, angle = 0, hjust = 0.5))+
  labs(x= "", y = "Mean difference") +
  #scale_color_manual(values=c("#005083",  "#f0b64d"))# +
  #scale_fill_manual(values=c("#969696",  "#969696"))
  #annotate("point", x = 1, y = 100, colour = "firebrick", size = 2)+
  #annotate(xmin = -Inf, xmax = Inf, ymin = 4.77-1.96*0.297, ymax = 4.77+1.96*0.297, geom = "rect", alpha = 0.2, fill = "firebrick")+
  #annotate(xmin = -Inf, xmax = Inf, ymin = 5.02-1.96*0.372, ymax = 5.02+1.96*0.372, geom = "rect", alpha = 0.2, fill = "blue")+
  #annotate(geom = "hline",yintercept = 100, y = 100, color = "red")+
  annotate("pointrange", x = 1, y = Mdiff, ymin = CI_low, ymax = CI_up, colour = "black", size = 1.5, shape = 24, fill = "darkgrey")+
  #annotate("pointrange", x = 2, y = 5.02, ymin = 5.02-1.96*0.372, ymax = 5.02+1.96*0.372, colour = "blue", size = 0.8, shape = 15)+
  #annotate("text", x = 0.5, y = 2.6, family = "Poppins", size = 6, color = "gray20", label = "Impfeffekt")+
  #geom_curve(aes(x = 0.5, y = 3, xend = 0.9, yend = 4),arrow = arrow(length = unit(0.03, "npc")),color = "gray20", curvature = +0.2)+
  #annotate("text", x = 1.8, y = 2.6, family = "Poppins", size = 6, color = "gray20", label = "Dosierungseffekt")+
  #geom_curve(aes(x = 1.8, y = 3, xend = 2, yend = 4),arrow = arrow(length = unit(0.03, "npc")),color = "gray20", curvature = +0.2)+
  annotate(geom = "hline",yintercept = 0, y = 0, color = "red", size = 1.2)+
  theme(plot.background = element_rect(
    fill = "white",
    colour = "black",
    size = 1
  ))


barchart

#ggsave("delta_joint.svg",width=2.5,height=4)
#ggsave("delta_joint.pdf",width=2.5,height=4)



```



```{r}
# compute the confidence interval for the singular causation differences in each between-subject condition

# separate presentation
t <- qt(0.975, 716, lower.tail = TRUE, log.p = FALSE)
#t

effect <- "Mdiff"
Mdiff <- 0.0661 
SE <- 0.0243
CI <- SE*t
CI_low <- Mdiff - CI
CI_up <- Mdiff + CI

Mdiff
CI_low
CI_up



# Plot 

myTheme <- theme(plot.title = element_text(face="bold", size = 22),
        axis.title.x = element_text(face = "bold", size = 20),
        axis.title.y = element_text(face = "bold", size = 20),
        axis.text.x = element_text(size = 18, angle = 0), 
        axis.text.y = element_text(size = 25, 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"))

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



barchart <- ggplot()+
  myTheme+
  #guides(fill=FALSE)+
  #facet_wrap(~Latency + SampleSize, ncol=2)+
  #ggtitle("Mean difference (95% CI)") +
  #coord_cartesian(ylim=c(-0.1,2)) + 
  scale_y_continuous(limits = c(-0.1, 0.5), breaks=seq(-0.1, 0.5, 0.1), expand = c(0,0)) +
  scale_x_discrete(labels=c("r")) +
  #annotate("rect", xmin=1.7, xmax=2.3, ymin=0.95, ymax=1.05, color="#31a354", fill = "white", size = 1) +
  #stat_summary(fun.y=mean, colour="grey20", geom="point", shape = 21, size = 3)+
  #stat_summary(fun.y = mean, geom = "bar", position = "dodge", colour = "black")+
  #stat_summary(fun.data = mean_cl_boot, geom = "errorbar", position = position_dodge(width = 0.90), width = 0.2) +
  #geom_jitter(width = 0.3, height = 0.02, alpha = 0.6, colour = "red") +
  #ggtitle("Means (95% bootstr. CIs)") +
  #theme(axis.text.x = element_text(size = 10, angle = 0, hjust = 0.5))+
  labs(x= "", y = "Mean difference") +
  #scale_color_manual(values=c("#005083",  "#f0b64d"))# +
  #scale_fill_manual(values=c("#969696",  "#969696"))
  #annotate("point", x = 1, y = 100, colour = "firebrick", size = 2)+
  #annotate(xmin = -Inf, xmax = Inf, ymin = 4.77-1.96*0.297, ymax = 4.77+1.96*0.297, geom = "rect", alpha = 0.2, fill = "firebrick")+
  #annotate(xmin = -Inf, xmax = Inf, ymin = 5.02-1.96*0.372, ymax = 5.02+1.96*0.372, geom = "rect", alpha = 0.2, fill = "blue")+
  #annotate(geom = "hline",yintercept = 100, y = 100, color = "red")+
  annotate("pointrange", x = 1, y = Mdiff, ymin = CI_low, ymax = CI_up, colour = "black", size = 1.5, shape = 24, fill = "darkgrey")+
  #annotate("pointrange", x = 2, y = 5.02, ymin = 5.02-1.96*0.372, ymax = 5.02+1.96*0.372, colour = "blue", size = 0.8, shape = 15)+
  #annotate("text", x = 0.5, y = 2.6, family = "Poppins", size = 6, color = "gray20", label = "Impfeffekt")+
  #geom_curve(aes(x = 0.5, y = 3, xend = 0.9, yend = 4),arrow = arrow(length = unit(0.03, "npc")),color = "gray20", curvature = +0.2)+
  #annotate("text", x = 1.8, y = 2.6, family = "Poppins", size = 6, color = "gray20", label = "Dosierungseffekt")+
  #geom_curve(aes(x = 1.8, y = 3, xend = 2, yend = 4),arrow = arrow(length = unit(0.03, "npc")),color = "gray20", curvature = +0.2)+
  annotate(geom = "hline",yintercept = 0, y = 0, color = "red", size = 1.2)+
  theme(plot.background = element_rect(
    fill = "white",
    colour = "black",
    size = 1
  ))


barchart

#ggsave("delta_sep.svg",width=2.5,height=4)
#ggsave("delta_sep.pdf",width=2.5,height=4)



```




```{r}
# compute the confidence interval for the singular causation differences in each between-subject condition

# interaction

t <- qt(0.975, 716, lower.tail = TRUE, log.p = FALSE)
#t

effect <- "Mdiff"
Mdiff <- 0.104 
SE <- 0.0344
CI <- SE*t
CI_low <- Mdiff - CI
CI_up <- Mdiff + CI

Mdiff
CI_low
CI_up



# Plot 

myTheme <- theme(plot.title = element_text(face="bold", size = 22),
        axis.title.x = element_text(face = "bold", size = 20),
        axis.title.y = element_text(face = "bold", size = 20),
        axis.text.x = element_text(size = 18, angle = 0), 
        axis.text.y = element_text(size = 25, 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"))

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



barchart <- ggplot()+
  myTheme+
  #guides(fill=FALSE)+
  #facet_wrap(~Latency + SampleSize, ncol=2)+
  #ggtitle("Mean difference (95% CI)") +
  #coord_cartesian(ylim=c(-0.1,2)) + 
  scale_y_continuous(limits = c(-0.1, 0.5), breaks=seq(-0.1, 0.5, 0.1), expand = c(0,0)) +
  scale_x_discrete(labels=c("r")) +
  #annotate("rect", xmin=1.7, xmax=2.3, ymin=0.95, ymax=1.05, color="#31a354", fill = "white", size = 1) +
  #stat_summary(fun.y=mean, colour="grey20", geom="point", shape = 21, size = 3)+
  #stat_summary(fun.y = mean, geom = "bar", position = "dodge", colour = "black")+
  #stat_summary(fun.data = mean_cl_boot, geom = "errorbar", position = position_dodge(width = 0.90), width = 0.2) +
  #geom_jitter(width = 0.3, height = 0.02, alpha = 0.6, colour = "red") +
  #ggtitle("Means (95% bootstr. CIs)") +
  #theme(axis.text.x = element_text(size = 10, angle = 0, hjust = 0.5))+
  labs(x= "", y = "Delta Mean difference") +
  #scale_color_manual(values=c("#005083",  "#f0b64d"))# +
  #scale_fill_manual(values=c("#969696",  "#969696"))
  #annotate("point", x = 1, y = 100, colour = "firebrick", size = 2)+
  #annotate(xmin = -Inf, xmax = Inf, ymin = 4.77-1.96*0.297, ymax = 4.77+1.96*0.297, geom = "rect", alpha = 0.2, fill = "firebrick")+
  #annotate(xmin = -Inf, xmax = Inf, ymin = 5.02-1.96*0.372, ymax = 5.02+1.96*0.372, geom = "rect", alpha = 0.2, fill = "blue")+
  #annotate(geom = "hline",yintercept = 100, y = 100, color = "red")+
  annotate("pointrange", x = 1, y = Mdiff, ymin = CI_low, ymax = CI_up, colour = "black", size = 1.5, shape = 23, fill = "darkgrey")+
  #annotate("pointrange", x = 2, y = 5.02, ymin = 5.02-1.96*0.372, ymax = 5.02+1.96*0.372, colour = "blue", size = 0.8, shape = 15)+
  #annotate("text", x = 0.5, y = 2.6, family = "Poppins", size = 6, color = "gray20", label = "Impfeffekt")+
  #geom_curve(aes(x = 0.5, y = 3, xend = 0.9, yend = 4),arrow = arrow(length = unit(0.03, "npc")),color = "gray20", curvature = +0.2)+
  #annotate("text", x = 1.8, y = 2.6, family = "Poppins", size = 6, color = "gray20", label = "Dosierungseffekt")+
  #geom_curve(aes(x = 1.8, y = 3, xend = 2, yend = 4),arrow = arrow(length = unit(0.03, "npc")),color = "gray20", curvature = +0.2)+
  annotate(geom = "hline",yintercept = 0, y = 0, color = "red", size = 1.2)+
  theme(plot.background = element_rect(
    fill = "white",
    colour = "black",
    size = 1
  ))


barchart

#ggsave("interaction.svg",width=2.5,height=4)
#ggsave("interaction.pdf",width=2.5,height=4)



```




```{r}
# compute Cohen's d

library(rstatix)

# subset for the four panels shown in the figure
joint_pres <- subset(tdata_sub, learning_format == "joint presentation")
sep_pres <- subset(tdata_sub, learning_format == "separate presentation")



joint_pres %>% cohens_d(strength_rating ~ target_cause, paired = FALSE)

sep_pres %>% cohens_d(strength_rating ~ target_cause, paired = FALSE)




d_joint <- 0.7051869	
d_sep <- 0.3008925	




```


Compute Cohen's ds and their SE, CI. 

```{r}
library(dplyr) # for pipe operator
# using the esc package because it gives SE of the ES directly
library(esc)

# condition: joint presentation of causes

joint <- subset(tdata_sub, learning_format == "joint presentation")



m1 <- joint %>%
          filter(target_cause == "single-effect cause")%>%
          summarize(Mean1 = mean(strength_rating))

sd1 <- joint %>%
          filter(target_cause == "single-effect cause")%>%
          summarize(SD1 = sd(strength_rating))


m2 <- joint %>%
          filter(target_cause == "multiple-effects cause")%>%
          summarize(Mean2 = mean(strength_rating))

sd2 <- joint %>%
          filter(target_cause == "multiple-effects cause")%>%
          summarize(SD2 = sd(strength_rating))



esc_mean_sd(
  grp1m = m1[,1], grp1sd = sd1[,1], grp1n = 360/2,
  grp2m = m2[,1], grp2sd = sd2[,1], grp2n = 360/2,
  es.type = "d"
)
```

```{r}
library(dplyr) # for pipe operator
# using the esc package because it gives SE of the ES directly
library(esc)



# condition: joint presentation of causes

sep <- subset(tdata_sub, learning_format == "separate presentation")



m1 <- sep %>%
          filter(target_cause == "single-effect cause")%>%
          summarize(Mean1 = mean(strength_rating))

sd1 <- sep %>%
          filter(target_cause == "single-effect cause")%>%
          summarize(SD1 = sd(strength_rating))


m2 <- sep %>%
          filter(target_cause == "multiple-effects cause")%>%
          summarize(Mean2 = mean(strength_rating))

sd2 <- sep %>%
          filter(target_cause == "multiple-effects cause")%>%
          summarize(SD2 = sd(strength_rating))



esc_mean_sd(
  grp1m = m1[,1], grp1sd = sd1[,1], grp1n = 360/2,
  grp2m = m2[,1], grp2sd = sd2[,1], grp2n = 360/2,
  es.type = "d"
)
```




