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


```{r}
tdata <- read.delim("data.txt", header=TRUE, sep="\t", na.strings="NA", dec=".", strip.white=TRUE)
```




# Results 

## Demographics

```{r}
# demographics 

tdata_age <- tdata 

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

1 = male, 2 = female, 3 = non-binary




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

# 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:7)


# 2. reshape into long format 
tdata_long <- tdata_sub %>% gather(variable, value, 5:7)


tdata_long$variable <- factor(tdata_long$variable, levels = c("single_strength_rating", "multiple_strength_rating", "no_conn_strength_rating"), 
                          labels = c("SC", "CC", "NC"))


tdata_long$value <- tdata_long$value * 0.01 # Recode DV so that values range from 0 to 1

```

# Graphs

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

tdata_sub <- tdata_long


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 

g <- ggplot(tdata_sub, aes(x=variable, y=valueJitter, group = sID)) +
  guides(fill=FALSE)+
  #facet_grid(Query_order ~ Cause_order)+
  #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 \ncause", "Common \ncause", "No \ncause")) +
  #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.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.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 = "Number Cause's Effects", 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=6,height=4.3)
#ggsave("results_lines.pdf",width=6,height=4.3)
```


# 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), stat.desc , basic = FALSE)

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



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

a1 <- aov_car(value ~ variable + Error(sID/(variable)), tdata_sub)
a1


```
Sign. Effect of "Variable" most likely due to the "non-cause" condition (see graph).


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

# follow-up analysis to get means and CIs

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

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

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

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

confint(t, level = 0.95)
```


No dilution 


Make a difference plot: 


```{r}

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

effect <- "Mdiff"
Mdiff <- 0.010
CI_low <- -0.0621
CI_up <- 0.0821

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.svg",width=2.5,height=4)
#ggsave("delta.pdf",width=2.5,height=4)
```

What Cohen's d is this? 


```{r}

dat <- tdata_sub


# since we have a repeated-meausres design, we now need the correlations of the ratings
library(dplyr) # for pipe operator
tdata -> t
r <- cor(t$single_strength_rating, t$multiple_strength_rating)
r


# now compute ES and SE and CI of it
# using the esc package because it gives SE of the ES directly
library(esc)

# get means and sds
m1 <- dat %>%
          filter(variable == "SC")%>%
          summarize(Mean1 = mean(value))

sd1 <- dat %>%
          filter(variable == "SC")%>%
          summarize(SD1 = sd(value))


m2 <- dat %>%
          filter(variable == "CC")%>%
          summarize(Mean2 = mean(value))

sd2 <- dat %>%
          filter(variable == "CC")%>%
          summarize(SD2 = sd(value))



d <- esc_mean_sd(
  grp1m = m1[,1], grp1sd = sd1[,1], grp1n = length(dat$sID)/3,
  grp2m = m2[,1], grp2sd = sd2[,1], grp2n = length(dat$sID)/3,
  r = r,
  es.type = "d"
)
d

d$ci.lo
d$ci.hi

d_ci <- (d$ci.hi - d$ci.lo)/2
d_ci
```

```{r}
# How many subjects would we need to make the ci-width small enough?

d <- esc_mean_sd(
  grp1m = m1[,1], grp1sd = sd1[,1], grp1n = 216,
  grp2m = m2[,1], grp2sd = sd2[,1], grp2n = 216,
  r = r,
  es.type = "d"
)
d

d_ci <- (d$ci.hi - d$ci.lo)/2
d_ci


```




