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


```{r, echo = FALSE}

# to create a 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_sub <- tdata_sub %>% gather(variable, value, 4:6)

tdata_sub$variable <- factor(tdata_sub$variable, levels = c("X1_strength_rating", "X3_strength_rating", "X7_strength_rating"), 
                          labels = c("one", "three", "seven"))

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 

capacity.labs <- c("Capacity: no info", "Capacity: same", "Capacity: different")
names(capacity.labs) <- c("none", "same", "different")


g <- ggplot(tdata_sub, aes(x=variable, y=valueJitter, group = sID)) +
  guides(fill=FALSE)+
  #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", "Multiple-effects \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.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=15.5,height=9)
#ggsave("results_graph_lines.pdf",width=15.5,height=9)
```





# 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 ~ Crystal_Order*Q_Order + Error(sID/(variable)), tdata_sub)
a1


```



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

# follow-up analysis 

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

Compute Cohen's d for the two dilution effects (1 vs 3 and 1 vs 7 effects)

```{r}


# one vs. three

tdata_13 <- subset(tdata_sub, variable != "seven")
tdata_13 <- subset(tdata_13, select = c(1:3,5:6))


# one vs. seven 

tdata_17 <- subset(tdata_sub, variable != "three")
tdata_17 <- subset(tdata_17, select = c(1:3,5:6))


```



```{r}

dat <- tdata_13


# since we have a repeated-meausres design, we now need the correlations of the ratings
library(dplyr) # for pipe operator
dat_wide <- spread(tdata_13, variable, value)
dat_wide -> t
r <- cor(t$one, t$three)
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 == "one")%>%
          summarize(Mean1 = mean(value))

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


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

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



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


```{r}

dat <- tdata_17


# since we have a repeated-meausres design, we now need the correlations of the ratings
library(dplyr) # for pipe operator
dat_wide <- spread(tdata_17, variable, value)
dat_wide -> t
r <- cor(t$one, t$seven)
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 == "one")%>%
          summarize(Mean1 = mean(value))

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


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

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



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













