# Exp1_control
# Author: Hongyu Li 2020

library(plyr)
library(dplyr)
library(Hmisc)
library(tidyverse)
library(cowplot)
library(scales)
library(ggplot2)
library(lme4)
library(twosamples)


##################################
# Define the funciton for split violin plot (to be used later)
GeomSplitViolin <- ggproto("GeomSplitViolin", GeomViolin, 
                           draw_group = function(self, data, ..., draw_quantiles = NULL) {
                             data <- transform(data, xminv = x - violinwidth * (x - xmin), xmaxv = x + violinwidth * (xmax - x))
                             grp <- data[1, "group"]
                             newdata <- plyr::arrange(transform(data, x = if (grp %% 2 == 1) xminv else xmaxv), if (grp %% 2 == 1) y else -y)
                             newdata <- rbind(newdata[1, ], newdata, newdata[nrow(newdata), ], newdata[1, ])
                             newdata[c(1, nrow(newdata) - 1, nrow(newdata)), "x"] <- round(newdata[1, "x"])
                             
                             if (length(draw_quantiles) > 0 & !scales::zero_range(range(data$y))) {
                               stopifnot(all(draw_quantiles >= 0), all(draw_quantiles <=
                                                                         1))
                               quantiles <- ggplot2:::create_quantile_segment_frame(data, draw_quantiles)
                               aesthetics <- data[rep(1, nrow(quantiles)), setdiff(names(data), c("x", "y")), drop = FALSE]
                               aesthetics$alpha <- rep(1, nrow(quantiles))
                               both <- cbind(quantiles, aesthetics)
                               quantile_grob <- GeomPath$draw_panel(both, ...)
                               ggplot2:::ggname("geom_split_violin", grid::grobTree(GeomPolygon$draw_panel(newdata, ...), quantile_grob))
                             }
                             else {
                               ggplot2:::ggname("geom_split_violin", GeomPolygon$draw_panel(newdata, ...))
                             }
                           })

geom_split_violin <- function(mapping = NULL, data = NULL, stat = "ydensity", position = "identity", ..., 
                              draw_quantiles = NULL, trim = TRUE, scale = "area", na.rm = FALSE, 
                              show.legend = NA, inherit.aes = TRUE) {
  layer(data = data, mapping = mapping, stat = stat, geom = GeomSplitViolin, 
        position = position, show.legend = show.legend, inherit.aes = inherit.aes, 
        params = list(trim = trim, scale = scale, draw_quantiles = draw_quantiles, na.rm = na.rm, ...))}

  
  
##########################################  
# Test outcomes (Proportions of choices)

control<-read.csv("./Exp1_proportion.csv")%>%
  dplyr::mutate(proportion=choice_correct/n)

control$ï..test_type<-factor(control$ï..test_type,levels=c("Negative control","Species","Pulse number"))

fig_control<-
  ggplot(data=control,aes(x=ï..test_type,y=proportion))+
  geom_point(position = position_dodge(width = 0.2)) +
  geom_errorbar(position = position_dodge(width = 0.2),aes(ymin=CI_low, ymax=CI_high,width=0.2))+
  geom_hline(yintercept=0.5,linetype="dashed") +
  xlab("Choice tests") +
  ylab("Proportions of subjects") +
  ylim(0:1)+
  theme_classic()+
  theme(text=element_text(size=16))+
  theme(axis.text=element_text(colour="black"))+
  scale_x_discrete(labels =c("Negative control","Species \n discrimination","Pulse number \n discrimination"))

save_plot("fig_control.png",fig_control,base_width=6,base_height=4, dpi = 300)
save_plot("fig_control.eps",fig_control,base_width=6,base_height=4)

# Test whether there's a side bias
raw<-read.csv("./Exp1_NegativeControl.csv",na.strings = "")
raw$SPK_left<-as.factor(raw$SPK_left)

control_negative<-dplyr::select(raw,SPK_left,choice_side,order_left)%>%
  drop_na()

summary(control_negative)
binom.test(x=13,n=34,0.5)

glm.1<-glm(choice_side~ 1, data=control_negative, family=binomial(link = "logit"))
glm.2<-glm(choice_side~ 1+order_left, data=control_negative, family=binomial(link = "logit"))
glm.3<-glm(choice_side~ 1+SPK_left, data=control_negative, family=binomial(link = "logit"))
lmtest::lrtest(glm.1,glm.2) #order not significant
lmtest::lrtest(glm.1,glm.3) #Speaker not significant

# Generate a statistical table
texreg::texreg(list(glm.1,glm.2,glm.3),digits=3,stars=NULL)
# Use AICc to replace AIC
AICc(glm.1)
AICc(glm.2)
AICc(glm.3)

#####################################################
# Latency

control_latency<-read.csv("./Exp1_latency.csv")%>%
  tidyr::drop_na(latency)

control_latency$test.type<-factor(control_latency$test.type,levels=c("negative_control","Hch_Hver","34_26"))
control_latency$choice<-factor(control_latency$choice,levels=c("L","R","34","Hch","26"))

fig_control_latency_box<-
  ggplot(data=control_latency,aes(x=test.type,y=latency,fill=choice)) +
  geom_boxplot(width=0.5,positio = position_dodge2(width =0.5, preserve = "single"))+
  ylab("Response latency (s)") +
  xlab("Stimulus types")+
  theme_classic()+
  theme(text=element_text(size=16))+
  theme(axis.text=element_text(colour="black"))+
  scale_x_discrete(labels =c("Negative control","Species \n discrimination","Pulse number \n discrimination"))+
  theme(legend.position="top", 
        legend.justification="left",
        legend.direction='horizontal',
        legend.box.margin=margin(c(0.5,0.5,0.5,0.5)),
        legend.title=element_blank())+
  scale_fill_manual(breaks=c("34","Hch","26"),values=c("grey88","grey88","grey67","grey88","grey100"),labels = c("34 pulses","H. chrysoscelis","26 pulses"))

fig_control_latency_violin<-
  ggplot(data=control_latency,aes(x=test.type,y=latency,fill=choice)) +
  geom_split_violin(scale="count")+
  # stat_summary(fun.data="mean_sdl", fun.args=list(mult=1))+
  ylab("Response latency") +
  xlab("Tesy types")+
  theme_classic()

save_plot("fig_control_latency_box.png",fig_control_latency,base_width=6,base_height=3, dpi = 300)
save_plot("fig_control_latency_box.eps",fig_control_latency,base_width=6,base_height=3)

# non-parametric Aderson Darling test
# AD-tests involve boot strapping process. The P-values will vary each time running the test.
pulsenum_latency<-subset(control_latency,test.type=="34_26")
pulsenum_latency$choice<-droplevels(pulsenum_latency$choice)
ad_test(pulsenum_latency[pulsenum_latency$choice=="34",]$latency,pulsenum_latency[pulsenum_latency$choice=="26",]$latency)

NegCont_latency<-subset(control_latency,test.type=="negative_control")
NegCont_latency$choice<-droplevels(NegCont_latency$choice)
ad_test(NegCont_latency[NegCont_latency$choice=="L",]$latency,NegCont_latency[NegCont_latency$choice=="R",]$latency)