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

# Data handling

```{r}
tdata <- read_csv("Supp_Exp_data.csv")
#tdata <- read.table("dummy_dat.csv", sep=';', header = T)

#tdata <- subset(tdata, subj_code != "RhoYo25w8umI") # This person needs to be excluded because they definitely did not understand the task

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


# Results

```{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 <- tdata
tdata_long <- tdata_sub
tdata_long$disease1 <- factor(tdata_long$disease1, levels = c("Vilosa", "Pylium"), 
                               labels = c("Vilosa", "Pylium"))

tdata_long$effects1 <- factor(tdata_long$effects1, levels = c("Lian", "Gludon"), 
                               labels = c("Lian", "Gludon"))

tdata_long$source <- factor(tdata_long$source, levels = c("no_info", "none", "expl_info"), 
                            labels = c("no info", "'don't know'", "unreadable"))

# Recode ratings so that high values represent Hn and low values represent Hw


dat_narrow_left <- subset(tdata_long, scale_orientation == "narrow_left")
dat_narrow_left$rating_rec <- (dat_narrow_left$dv_rating) - 5

dat_narrow_right <- subset(tdata_long, scale_orientation == "narrow_right")
dat_narrow_right$rating_rec <- (dat_narrow_right$dv_rating * -1) + 5

tdata_long <- rbind(dat_narrow_left, dat_narrow_right)





table(tdata_long$source, tdata_long$disease1, tdata_long$effects1)
```

## Graphs

```{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 = 30, vjust = 0.7), 
        axis.text.x = element_text(size = 14), 
        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$subj_code <- factor(tdata_sub$subj_code)

pd <- position_dodge(width = 0.3)

tdata_sub$valueJitter <- jitter(tdata_sub$rating_rec, 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 = source, y = valueJitter)) +
  guides(fill=FALSE)+
  #facet_grid( ~ Mutation1, Effects1)+
  #ggtitle("Subjects' causal srength ratings") +
  scale_y_continuous(limits = c(-5.3, 5.3), breaks=seq(-5, 5, 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 = rating_rec, group = source, fill = source), color = NA, position=position_dodge(1), alpha = 0.4)+
  #geom_line(position = pd, color = "black", size = 1, alpha=0.04) +
  geom_hline(yintercept=0, linetype="dashed", color = "black")+
  geom_jitter(aes(color = source), alpha = 0.5, width = 0.15, height = 0.2) +
  stat_summary(aes(y = rating_rec,group=1), fun.data = mean_cl_boot, geom = "errorbar", width = 0, size = 1) +
  stat_summary(aes(y = rating_rec,group=1, fill = source), fun.y=mean, geom="point", color = "black", shape = 22, size = 2, group=1, alpha = 1)+
  stat_summary(aes(y = rating_rec,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 = "Information about latent feature", y = "Explanation rating") +
  scale_color_manual(name = "Strength",values=c("#66c2a5", "#e78ac3", "#8da0cb", "#a6d854"))+
  scale_fill_manual(name = "Strength",values=c("#66c2a5", "#e78ac3", "#8da0cb", "#a6d854"))+
  annotate("text", x = 0.5, y = 3.5, label = c("broad-scope"), angle = 90)+
  annotate("text", x = 0.5, y = -3.5, label = c("narrow-scope"), angle = 90)+
  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)
g


ggsave("results_means.svg",width=5.5,height=5)
ggsave("results_means.pdf",width=5.5, height=5)
```


## Stats

Use the ggplot_build package to see a table with the means and CI values plotted in the graph: 

```{r}
values <- ggplot_build(g)$data[[4]] # values are shown in the 4th panel 
values
```


```{r}
library(rcompanion)


ci_table <- groupwiseMean(rating_rec ~ source,
              data        = tdata_long,
              traditional = FALSE,
              percentile  = TRUE)


ci_table
(ci_width <- ci_table$Percentile.upper - ci_table$Percentile.lower)


```


get group medians: 

```{r}
library(rcompanion)
# groupwiseMean(rating_rec ~ Features + Knowledge,
#               data        = tdata_long,
#               traditional = FALSE,
#               percentile  = TRUE)

groupwiseMedian(rating_rec ~ source,
                data        = tdata_long,
                bca         = FALSE,
                percentile  = TRUE,
                R           = 1000)
```



Count the different ratings observed:

```{r}
counts <- tdata_long %>%
  group_by(source, rating_rec) %>%
  summarize(n = n()) %>% 
  mutate(pct = n/sum(n),
         lbl = scales::percent(pct))

counts
```

Make a linear trend analysis to test if bias becomes smaller across conditions: 


```{r}
contrasts(tdata_long$source) <- "contr.poly" # fits a polynomial contrast 
LinearModel <- lm(rating_rec ~ source, data=tdata_long) # runs a linear model with the corresponding contrast coding
summary(LinearModel)
```

## Graph for counts

Make a plot with the different rating categories: 

```{r}

tdata_long$category[tdata_long$rating_rec < 0] <- "narrow"
tdata_long$category[tdata_long$rating_rec == 0] <- "unbiased"
tdata_long$category[tdata_long$rating_rec > 0] <- "broad"

counts2 <- tdata_long %>%
  group_by(source, category) %>%
  summarize(n = n()) %>% 
  mutate(pct = n/sum(n),
         lbl = scales::percent(pct))

counts2

counts2$category <- factor(counts2$category, levels = c("unbiased", "narrow", "broad"), labels = c("unbiased", "narrow l.s.", "broad l.s."))

```


Get proportion CIs for different conditions:

```{r}
library(PropCIs)
library(DescTools)
library(purrr)

counts_noinfo <- subset(counts2, source == "no info")
counts_dontknow <- subset(counts2, source == "'don't know'")
counts_unreadable <- subset(counts2, source == "unreadable")

(MultinomCI(counts_noinfo$n,
           conf.level=0.95,
           method="sisonglaz") -> selection_ci_noinf)
  

(MultinomCI(counts_dontknow$n,
           conf.level=0.95,
           method="sisonglaz") -> selection_ci_dontknow)

(MultinomCI(counts_unreadable$n,
           conf.level=0.95,
           method="sisonglaz") -> selection_ci_unreadable)
  



ci_low <- c(selection_ci_noinf[,2], selection_ci_dontknow[,2],  selection_ci_unreadable[,2])

ci_up <- c(selection_ci_noinf[,3], selection_ci_dontknow[,3],  selection_ci_unreadable[,3])


```


```{r}
plotdata <- counts2

plotdata$ci_low <- ci_low
plotdata$ci_up <- ci_up
```






Plot:

```{r}
plotdata <- counts2

library(scales)
theme_set(theme_light(base_size = 12, base_family = "Poppins"))

g<- ggplot(plotdata, 
       aes(x = category,
           y = pct,
           fill = category)) +
  facet_grid( ~ source)+
  geom_bar(stat = "identity",
           position = "dodge") +
  scale_y_continuous(limits = seq(0, 2),
                     breaks = seq(0, 1, .25),
                     expand = c(0,0),
                     label = percent) +
  #coord_cartesian(xlim =c(1, 7), ylim = c(0, 1.1))+
  #coord_cartesian(clip = "off")+
  geom_text(aes(label = lbl), 
            size = 3.5,
            position = position_dodge(width = 1),
            vjust = -3) +
  scale_fill_manual(name = "Strength",values=c("#66c2a5", "#e78ac3", "#8da0cb", "#a6d854"))+
  #scale_fill_brewer(palette = "Pastel1") +
  labs(y = "Percentage", 
       fill = "Explanatory preference",
       x = "Explanatory preference")+
  geom_pointrange(ymin = ci_low, ymax = ci_up, position = position_dodge(width = 0.89), shape = 22, size = 0.3)+
  #annotate(geom = "hline",yintercept = 0.5, y = 0.5, color = "black", size = 1, linetype='dotted')+
  #annotate("pointrange", x = plotdata$Transformation, y = plotdata$pct, 
   #        ymin = plotdata$ci_low, 
    #       ymax = plotdata$ci_up, 
     #      colour = "black", size = 0.8, shape = 22, fill = Transformation, fatten = 1)+
  #annotate("text", x = pvalues_x, y = Inf, label = pvalues, size = 4, vjust = 1.8)+
  theme(legend.position = "none", axis.title = element_text(size = 20), axis.text = element_text(size = 13, color = "black"),
        legend.text = element_text(size = 13),legend.title = element_text(size = 13), strip.text.x = element_text(size = 13))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

g

#ggsave("selections_between.pdf",width=6,height=5)
ggsave("categories.svg",width=10,height=5)
ggsave("categories.pdf",width=10,height=5)
```

## Proportion tests 

Test if more subjects show the bias in the "don't know" as copmared to the "unreadable" condition. 

```{r}
prop.test(x = c(plotdata$n[8], plotdata$n[5]), n = c(48, 48), alternative = "less", correct = F)
 
```

As predicted, proportion of subjects committing the narrow latent scope bias is significantly smaller in the "explanation" condition (unreadable).

Also test if more subjects show the bias in the "no information condition" than in the "unreadable condition".

```{r}
prop.test(x = c(plotdata$n[8], plotdata$n[2]), n = c(48, 48), alternative = "less", correct = F)
```



