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


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




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

```




# Results 

## Demographics

```{r}
# demographics 

# one participant indicated 3 for age. Needs to be excluded for the age analysis.
tdata_age <- subset(tdata, Age > 3) # one participant erronously reported 3 

min(tdata_age$Age)
max(tdata_age$Age)
mean(tdata_age$Age)
sd(tdata_age$Age)

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


```{r, echo = FALSE}
# reorder factor 

tdata$Order <- factor(tdata$Order, levels = c("multiple_first", "single_first"), labels = c("multiple_first", "single_first"))

tdata$Multiple_Effects <- factor(tdata$Multiple_Effects, levels = c("blue", "red"), labels = c("blue", "red"))

tdata$Process <- factor(tdata$Process, levels = c("generative", "preventive"), labels = c("generative", "preventive"))

tdata$Valence <- factor(tdata$Valence, levels = c("positive", "negative"), labels = c("positive", "negative"))

tdata$Target <- factor(tdata$Target, 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:8)


# 2. reshape into long format 

tdata_sub <- tdata_sub %>% gather(variable, value, 7:8)

tdata_sub$variable <- factor(tdata_sub$variable, levels = c("single_strength_rating", "multiple_strength_rating"), 
                          labels = c("single", "multiple"))

tdata_sub$value <- tdata_sub$value * 0.01


```

# Graphs

```{r}
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"))


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$value, factor = 1, amount = 0.04)

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

# new labes for the facets 
process.labs <- c("Process: Generative", "Process: Preventive")
names(process.labs) <- c("generative", "preventive")

valence.labs <- c("Effects: Positive", "Effects: Negative")
names(valence.labs) <- c("positive", "negative")


g <- ggplot(tdata_sub, aes(x=variable, y=valueJitter, group = sID)) +
  guides(fill=FALSE)+
  facet_grid(Valence ~ Process, labeller = labeller(Process = process.labs, Valence = valence.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")) +
  #stat_summary(fun.y = mean, geom = "bar", position = "dodge", colour = "black", alpha =0.5) +
  geom_violinhalf(aes(y = value, group = variable, fill = variable), color = NA, position=position_dodge(1), alpha = 0.2)+
  geom_line(position = pd, color = "black", size = 1, alpha=0.07) +
  geom_point(aes(color = variable), position = pd, alpha = 0.4) +
  stat_summary(aes(y = value,group=1), fun.data = mean_cl_boot, geom = "errorbar", width = 0, size = 1) +
  stat_summary(aes(y = value,group=1), fun.y=mean, colour="black", geom="line",group=1, size = 1.5, linetype = "solid", alpha = 1)+
  stat_summary(aes(y = value,group=1, fill = variable), fun.y=mean, geom="point", color = "black", shape = 22, size = 5, group=1, alpha = 1)+
  stat_summary(aes(y = value,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 = "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_lines.svg",width=15.5,height=9)
#ggsave("Exp_Supp1_res_rev.pdf",width=8.3,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$value, list(tdata_sub$variable, tdata_sub$Process, tdata_sub$Valence), stat.desc , basic = FALSE)

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



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

a1 <- aov_car(value ~ Order*Process*Valence*Multiple_Effects*Target*variable + Error(sID/(variable)), tdata_sub)
a1


```



```{r}
# same ANOVA as before
lmeModel <- lmer(value ~ Process*variable + (1|sID), data=tdata_sub)

# follow-up analysis 

ls1 <- lsmeans(a1, c("variable", "Process", "Valence")) # joint evaluation (basically gives the same table)
ls1
```



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

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

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


ls3 <- lsmeans(a1, c("Process")) # group means by between-condition
ls3

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

```



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

ls4 <- lsmeans(a1, c("variable"), by = c("Process", "Valence")) # group means by between-condition
ls4

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

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

test(pairs(pairs(ls4), 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 Cohen's d

library(rstatix)

# subset for the four panels shown in the figure
gen_pos <- subset(tdata_sub, Process == "generative" & Valence == "positive")
gen_neg <- subset(tdata_sub, Process == "generative" & Valence == "negative")

prev_pos <- subset(tdata_sub, Process == "preventive" & Valence == "positive")
prev_neg <- subset(tdata_sub, Process == "preventive" & Valence == "negative")


gen_pos %>% cohens_d(value ~ variable, paired = TRUE)
gen_neg %>% cohens_d(value ~ variable, paired = TRUE)

prev_pos %>% cohens_d(value ~ variable, paired = TRUE)
prev_neg %>% cohens_d(value ~ variable, paired = TRUE)


d_gen_pos <- 0.9161883
d_gen_neg <- 1.138706

d_prev_pos <- 0.5932288
d_prev_neg <- 0.5132523

```



```{r}
# get confidence intervals for d 

# 1) compute correlations for the ratings

gen_pos_mult <- subset(gen_pos, variable == "multiple")
gen_pos_sing <- subset(gen_pos, variable == "single")

cor.test(gen_pos_sing$value, gen_pos_mult$value)

cor_pos_gen <- -0.1982053


prev_pos_mult <- subset(prev_pos, variable == "multiple")
prev_pos_sing <- subset(prev_pos, variable == "single")

cor.test(prev_pos_sing$value, prev_pos_mult$value)

cor_pos_prev <- 0.4342252


gen_neg_mult <- subset(gen_neg, variable == "multiple")
gen_neg_sing <- subset(gen_neg, variable == "single")

cor.test(gen_neg_sing$value, gen_neg_mult$value)

cor_neg_gen <- 0.006839188


prev_neg_mult <- subset(prev_neg, variable == "multiple")
prev_neg_sing <- subset(prev_neg, variable == "single")

cor.test(prev_neg_sing$value, prev_neg_mult$value)

cor_neg_prev <-0.004405479



# 2) Now compute SE for d

n <- 120


# formula: Sqrt((1/n + d^2/n)*2*(1-r))

SEd_pos_gen <- sqrt((1/n + d_gen_pos^2/n)*2*(1-cor_pos_gen))
SEd_pos_gen

SEd_neg_gen <- sqrt((1/n + d_gen_neg^2/n)*2*(1-cor_neg_gen))
SEd_neg_gen

SEd_pos_prev <- sqrt((1/n + d_prev_pos^2/n)*2*(1-cor_pos_prev))
SEd_pos_prev

SEd_neg_prev <- sqrt((1/n + d_prev_neg^2/n)*2*(1-cor_neg_prev))
SEd_neg_prev

```
```{r}
# Confidence intervalls for d


("d gen pos")
round(d_gen_pos,2)
round((d_gen_pos - 1.96*SEd_pos_gen),2)
round((d_gen_pos + 1.96*SEd_pos_gen),2)


("d prev pos")
round(d_prev_pos,2)
round((d_prev_pos - 1.96*SEd_pos_prev),2)
round((d_prev_pos + 1.96*SEd_pos_prev),2)


("d gen neg")
round(d_gen_neg,2)
round((d_gen_neg - 1.96*SEd_neg_gen),2)
round((d_gen_neg + 1.96*SEd_neg_gen),2)


("d prev neg")
round(d_prev_neg,2)
round((d_prev_neg - 1.96*SEd_neg_prev),2)
round((d_prev_neg + 1.96*SEd_neg_prev),2)


```



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

#Process = generative, Valence = positive:
# contrast          estimate     SE  df t.ratio p.value
# single - multiple    0.315 0.0295 432 10.673  <.0001 

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

effect <- "Mdiff"
Mdiff <- 0.315
SE <- 0.0295
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_blank(),
        axis.text.x = element_text(size = 18, angle = 0), 
        axis.text.y = element_text(size = 40, angle = 0),
        legend.text = element_text(size = 18),
        legend.title = element_text(size = 22),
        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 change") +
  #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 = 2, 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 = "white",
    size = 1
  ))


barchart

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



```


```{r}
# Process = preventive, Valence = positive:
# contrast          estimate     SE  df t.ratio p.value
# single - multiple    0.156 0.0295 432  5.274  <.0001 

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

effect <- "Mdiff"
Mdiff <- 0.156
SE <- 0.0295
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_blank(),
        axis.text.x = element_text(size = 18, angle = 0), 
        axis.text.y = element_text(size = 40, angle = 0),
        legend.text = element_text(size = 18),
        legend.title = element_text(size = 22),
        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 change") +
  #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 = 2, 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 = "white",
    size = 1
  ))


barchart

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



```


```{r}
#Process = generative, Valence = negative:
# contrast          estimate     SE  df t.ratio p.value
# single - multiple    0.354 0.0295 432 11.982  <.0001 

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

effect <- "Mdiff"
Mdiff <- 0.354
SE <- 0.0295
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_blank(),
        axis.text.x = element_text(size = 18, angle = 0), 
        axis.text.y = element_text(size = 40, angle = 0),
        legend.text = element_text(size = 18),
        legend.title = element_text(size = 22),
        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 change") +
  #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 = 2, 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 = "white",
    size = 1
  ))


barchart

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



```

```{r}
#Process = preventive, Valence = negative:
# contrast          estimate     SE  df t.ratio p.value
# single - multiple    0.182 0.0295 432  6.163  <.0001

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

effect <- "Mdiff"
Mdiff <- 0.182
SE <- 0.0295
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_blank(),
        axis.text.x = element_text(size = 18, angle = 0), 
        axis.text.y = element_text(size = 40, angle = 0),
        legend.text = element_text(size = 18),
        legend.title = element_text(size = 22),
        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 change") +
  #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 = 2, 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 = "white",
    size = 1
  ))


barchart

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



```



```{r}
# contrast                                                                          estimate     SE  df t.ratio p.value
# (single - multiple generative positive) - (single - multiple preventive positive)   0.1594 0.0418 432  3.817  0.0009 

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

effect <- "Mdiff"
Mdiff <- 0.1594 
SE <- 0.0418
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_blank(),
        axis.text.x = element_text(size = 18, angle = 0), 
        axis.text.y = element_text(size = 40, angle = 0),
        legend.text = element_text(size = 18),
        legend.title = element_text(size = 22),
        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 = 2, 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 = "white",
    size = 1
  ))


barchart

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



```

```{r}
# contrast                                                                          estimate     SE  df t.ratio p.value
#(single - multiple generative negative) - (single - multiple preventive negative)   0.1718 0.0418 432  4.115  0.0003

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

effect <- "Mdiff"
Mdiff <- 0.1718
SE <- 0.0418
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_blank(),
        axis.text.x = element_text(size = 18, angle = 0), 
        axis.text.y = element_text(size = 40, angle = 0),
        legend.text = element_text(size = 18),
        legend.title = element_text(size = 22),
        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 = 2, 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 = "white",
    size = 1
  ))


barchart

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



```


# Cluster Analysis

Goal: Identify and mark subjects who perceived dilution, who perceived strengthening, and who didn't discriminate between the cases.

```{r}

data_cluster <- tdata[-c(9:15)]

# append absolut deviations as new columns
data_cluster$Delta_Rating <- data_cluster$single_strength_rating - data_cluster$multiple_strength_rating



crit = 3 # Threshold criterion (needs to be exceeded to leave the invariance cluster)

data_cluster$Change[data_cluster$Delta_Rating > crit] <- "Dilution"
data_cluster$Change[data_cluster$Delta_Rating >= -crit & data_cluster$Delta_Rating <= crit] <- "Invariance"
data_cluster$Change[data_cluster$Delta_Rating < -crit] <- "Strengthening"



data_cluster %>% count(Change)

# now append the clustering results to the main data frame 
data_cluster <- subset(data_cluster, select = c(1,10))
tdata_sub <- merge(tdata_sub, data_cluster, by = c("sID"))


```

```{r}
Dilution <- prop.test(313,480,correct=FALSE)
Dilution

Invariance <- prop.test(124,480,correct=FALSE)
Invariance

Strengthening <- prop.test(43,480,correct=FALSE)
Strengthening
```


```{r, echo = FALSE}
myTheme <- theme(plot.title = element_text(face="bold", size = 22),
        axis.title.x = element_blank(),
        axis.title.y = element_text(face = "bold", size = 20),
        axis.text.x = element_text(size = 18, 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"))


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$value, factor = 1, amount = 0.04)

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

# new labes for the facets 
effects.labs <- c("Multiple effects from same domain", "Multiple effects from different domains")
names(effects.labs) <- c("same domain", "different domains")

palette <- c("#FFAE49", "#8cb495", "#12658b")



# new labes for the facets 
process.labs <- c("Process: Generative", "Process: Preventive")
names(process.labs) <- c("generative", "preventive")

valence.labs <- c("Effects: Positive", "Effects: Negative")
names(valence.labs) <- c("positive", "negative")


g <- ggplot(tdata_sub, aes(x=variable, y=valueJitter, group = sID, color = Change)) +
  guides(fill=FALSE)+
  facet_grid(Valence ~ Process, labeller = labeller(Process = process.labs, Valence = valence.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")) +
  #stat_summary(fun.y = mean, geom = "bar", position = "dodge", colour = "black", alpha =0.5) +
  geom_violinhalf(aes(y = value, group = variable), fill = "black", color = NA, position=position_dodge(1), alpha = 0.2)+
  geom_line(position = pd, size = 1, alpha=0.2) +
  geom_point(position = pd, alpha = 0.4, size = 1.5) +
  #geom_line(position = pd, color = "black", size = 1, alpha=0.04) +
  #geom_point(aes(color = variable), position = pd, alpha = 0.2) +
  stat_summary(aes(y = value,group=1), fun.data = mean_cl_boot, geom = "errorbar", width = 0, size = 1) +
  #stat_summary(aes(y = value,group=1), fun.data = mean_cl_boot, geom = "errorbar", width = 0, size = 1) +
  #stat_summary(aes(y = value,group=1), fun.y=mean, colour="black", geom="line",group=1, size = 1.5, linetype = "solid", alpha = 0.5)+
  stat_summary(aes(y = value,group=1, fill = variable), fun.y=mean, geom="point", color = "black", fill = "black", shape = 22, size = 4, 
               group=1, alpha = 0.5)+
  stat_summary(aes(y = value,group=1), fun.y=median, geom="point", color = "black", fill = "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 = "Model",values=palette, labels = c("Dilution: 65%", "Invariance: 26%", "Strengthening: 09%"))+
  scale_fill_manual(name = "Model",values=palette)+
  theme(legend.position = "bottom")+
  myTheme
g


#ggsave("Clusters.svg",width=14.5,height=8.3)
#ggsave("Clusters.pdf",width=10.5,height=7.5)
```

```{r}
# Generative vs. Preventive 
gen <- subset(tdata_sub, Process == "generative" | Process == "preventive")
gen %>% count(Change, by = Process)
# keep in mind that data are in long format i.e. each subject has to rows
```


```{r}
# More dilution if process is generative?
prop.test(c(368/2,257/2), c(240,240), p = NULL, alternative = "two.sided",
          correct = TRUE)
```






```{r}
# models in each of the four conditions 

# generative cause leading to positive effects
posgen <- subset(tdata_sub, Process == "generative" & Valence == "positive")
posgen %>% count(Change)


```

```{r}

Dilution <- prop.test(186/2,120,correct=FALSE)
Dilution

Invariance <- prop.test(40/2,120,correct=FALSE)
Invariance

Strengthening <- prop.test(14/2,120,correct=FALSE)
Strengthening
```







```{r}
# models in each of the four conditions 

# preventive cause leading to positive effects
posprev <- subset(tdata_sub, Process == "preventive" & Valence == "positive")
posprev %>% count(Change)


```

```{r}
Dilution <- prop.test(120,240,correct=FALSE)
Dilution

Invariance <- prop.test(100,240,correct=FALSE)
Invariance

Strengthening <- prop.test(20,240,correct=FALSE)
Strengthening
```

```{r}
# models in each of the four conditions 

# generative cause leading to negative effects
posprev <- subset(tdata_sub, Process == "generative" & Valence == "negative")
posprev %>% count(Change)

```


```{r}
Dilution <- prop.test(182,240,correct=FALSE)
Dilution

Invariance <- prop.test(44,240,correct=FALSE)
Invariance

Strengthening <- prop.test(14,240,correct=FALSE)
Strengthening
```


```{r}
# models in each of the four conditions 

# preventive cause leading to negative effects
posprev <- subset(tdata_sub, Process == "preventive" & Valence == "negative")
posprev %>% count(Change)

```


```{r}
Dilution <- prop.test(138,240,correct=FALSE)
Dilution

Invariance <- prop.test(64,240,correct=FALSE)
Invariance

Strengthening <- prop.test(38,240,correct=FALSE)
Strengthening
```

```{r}
# check if the proportion of strengtheners is higher in the prev_neg condition compared to all the others 

strengthening_prev_neg <- 27
strengthening_rest <- 45
n_prev_neg <- 120
n_rest <- 120*3

x <- c(strengthening_prev_neg,strengthening_rest)
n <- c(n_prev_neg,n_rest <- 120*3)

prop.test(x, n, p = NULL, alternative = "two.sided",
          correct = TRUE)


```











