# Exp2
# Author: Hongyu Li 2020

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

# Define 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, ...))}

######################################################################
# Data organization
preference_18<-read.csv("./Exp2_2018.csv")%>%
  dplyr::select(ï..date,frogID,test_type,choice,latency)
count_18<-plyr::count(preference_18,c("test_type","choice"))%>%
  spread(choice,freq)

preference_19<-read.csv("./Exp2_2019.csv")%>%
  dplyr::select(ï..date,frogID,SPL,choice,latency)
preference_19$SPL<-as.factor(preference_19$SPL)
count_19<-plyr::count(preference_19,c("SPL","choice"))%>%
  tidyr::drop_na()%>%
  spread(choice,freq)


#########################################################################
# acoustic vs acoustic+SM, acoustic vs acoustic+DM, acoustic+SM vs acoustic+DM at 88dB
# combine data from 2018 & 2019
db88_summary<-count_18
db88_summary[2,3]<-count_18[2,3]+count_19[3,2]
db88_summary[2,4]<-count_18[2,4]+count_19[3,3]

db88_summary[is.na(db88_summary)] <- 0
TOTAL<-db88_summary[,2]+db88_summary[,3]+db88_summary[,4]
db88_summary<-cbind(db88_summary,TOTAL)%>%
  gather(choice,freq,CALL:SM)

db88_summary[db88_summary==0]<-NA
db88_summary<-drop_na(db88_summary)%>%
  cbind(binconf(x=.$freq,n=.$TOTAL))

bt<-function(x,n,p = 0.5){binom.test(x,n,0.5)$p.value}
p.value<-mapply(bt,db88_summary$freq,db88_summary$TOTAL)
db88_summary<-cbind(db88_summary,p.value)

db88_summary$test_type<-factor(db88_summary$test_type,levels=c("DM vs call","SM vs call","DM vs SM"))

fig_preference_88dB<-
  ggplot(data=db88_summary[c(3,4,6),],aes(x=test_type,y=PointEst))+
  geom_point(position = position_dodge(width = 0.2)) +
  geom_errorbar(position = position_dodge(width = 0.2),aes(ymin=Lower, ymax=Upper,width=0.2))+
  geom_hline(yintercept=0.5,linetype="dashed") +
  xlab("Test types") +
  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(expression(DM^{"+"}~" vs "~NM^{"+"}),expression(SM^{"+"}~" vs "~NM^{"+"}),expression(DM^{"+"}~" vs "~SM^{"+"})))

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

## latency
preference_19_88db<-subset(preference_19,SPL=="88")
colnames(preference_19_88db)[3]<-"test_type"
db88_latency<-rbind(preference_18,preference_19_88db)
db88_latency$test_type[db88_latency$test_type=="88"]<-"DM vs SM"
db88_latency$test_type<-droplevels(db88_latency$test_type)
db88_latency$frogID<-paste(db88_latency$ï..date,db88_latency$frogID, sep="-")

fig_preference_88_latency_violin<-
  ggplot(data=db88_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("Stimulus types")+
  theme_classic()

db88_latency$test_type<-factor(db88_latency$test_type,levels=c("DM vs call","SM vs call","DM vs SM"))
db88_latency$choice<-factor(db88_latency$choice,levels=c("DM","SM","CALL"))

fig_preference_88_latency_box<-
  ggplot(data=db88_latency,aes(x=test_type,y=latency,fill=choice)) +
  geom_boxplot(width=0.5)+
  ylab("Response latency (s)") +
  xlab("Alternatives provided in choice tests")+
  scale_fill_manual(values=c("grey64","grey88","grey100"),labels = c(expression(DM^{"+"}),expression(SM^{"+"}),expression(NM^{"+"})))+
  theme_classic()+
  theme(legend.position="top", 
        legend.justification="left",
        legend.direction='horizontal',
        legend.box.margin=margin(c(0.5,0.5,0.5)),
        legend.title=element_blank())+
  theme(text=element_text(size=16))+
  theme(axis.text=element_text(colour="black"))+
  scale_x_discrete(labels = c(expression(DM^{"+"}~" vs "~NM^{"+"}),expression(SM^{"+"}~" vs "~NM^{"+"}),expression(DM^{"+"}~" vs "~SM^{"+"})))

save_plot("fig_preference_88_latency_box.png",fig_preference_88_latency,base_height=3,base_width=6,dpi = 300)
save_plot("fig_preference_88_latency_box.eps",fig_preference_88_latency,base_height=3,base_width=6)

###
# Aderson Darling test
latency_SMcall<-subset(db88_latency,test_type=="SM vs call")
latency_SMcall$choice<-droplevels(latency_SMcall$choice)

latency_DMcall<-subset(db88_latency,test_type=="DM vs call")
latency_DMcall$choice<-droplevels(latency_DMcall$choice)

latency_DMSM<-subset(db88_latency,test_type=="DM vs SM")
latency_DMSM$choice<-droplevels(latency_DMSM$choice)

library(twosamples)
ad_test(latency_DMcall[latency_DMcall$choice=="DM",]$latency,latency_DMcall[latency_DMcall$choice=="CALL",]$latency)
ad_test(latency_SMcall[latency_SMcall$choice=="SM",]$latency,latency_SMcall[latency_SMcall$choice=="CALL",]$latency)
ad_test(latency_DMSM[latency_DMSM$choice=="SM",]$latency,latency_DMSM[latency_DMSM$choice=="DM",]$latency)


#########################################################################
# call+SM vs call+DM with different SPLs
## proportion
# combine data from 2018 & 20 19
SPL_summary<-count_19
SPL_summary[3,2]<-count_19[3,2]+count_18[2,3]
SPL_summary[3,3]<-count_19[3,3]+count_18[2,4]

TOTAL<-SPL_summary$DM+SPL_summary$SM
SPL_summary<-cbind(SPL_summary,TOTAL)
SPL_summary<-cbind(SPL_summary,binconf(x=SPL_summary$DM,n=SPL_summary$TOTAL))

bt<-function(x,n,p = 0.5){binom.test(x,n,0.5)$p.value}
p.value<-mapply(bt,SPL_summary$DM,SPL_summary$TOTAL)
SPL_summary<-cbind(SPL_summary,p.value)

fig_preference_SPL<-
  ggplot(data=SPL_summary,aes(x=SPL,y=PointEst))+
  geom_point(position = position_dodge(width = 0.2)) +
  geom_errorbar(position = position_dodge(width = 0.2),aes(ymin=Lower, ymax=Upper,width=0.2))+
  geom_hline(yintercept=0.5,linetype="dashed") +
  xlab("SPLs of acoustic signals at 50cm (dB)") +
  ylab(expression("Proportions of subjects")) +
  ylim(0:1)+
  theme_classic()+
  theme(text=element_text(size=16))+
  theme(axis.text=element_text(colour="black"))

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

# Models
preference_18_88db<-subset(preference_18,test_type=="DM vs SM")
colnames(preference_18_88db)[3]<-"SPL"
preference_18_88db$SPL<- "88"
preference_SPL_latency<-rbind(preference_19,preference_18_88db)

preference_SPL_lm.0<-glmer(choice~1+(1|frogID),data=preference_SPL_latency,family=binomial(link = "logit"))
preference_SPL_lm.1<-glmer(choice~1+(1|frogID)+SPL,data=preference_SPL_latency,family=binomial(link = "logit"))
anova(preference_SPL_lm.0,preference_SPL_lm.1)

#####################
## latency
preference_18_88db<-subset(preference_18,test_type=="DM vs SM")
colnames(preference_18_88db)[3]<-"SPL"
preference_18_88db$SPL<-"88"
preference_SPL_latency<-rbind(preference_19,preference_18_88db)

preference_SPL_latency$choice<-factor(preference_SPL_latency$choice,levels=c("DM","SM"))

fig_preference_SPL_latency_box<-
  ggplot(data=preference_SPL_latency,aes(x=SPL,y=latency,fill=choice)) +
  geom_boxplot(width=0.5)+
  ylab("Response latency (s)") +
  xlab("SPLs of acoustic signals at 50cm (dB)")+
  theme_classic()+
  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())+
  theme(text=element_text(size=16))+
  theme(axis.text=element_text(colour="black"))+
  scale_fill_manual(values=c("grey76","grey88"),labels = c(expression(DM^{"+"}),expression(SM^{"+"})))

fig_preference_SPL_latency_violin<-
  ggplot(data=preference_SPL_latency,aes(x=SPL,y=latency,fill=choice)) +
  geom_split_violin(scale="count")+
  # stat_summary(fun.data="mean_sdl", fun.args=list(mult=1))+
  ylab("Response latency") +
  xlab("Stimulus types")+
  theme_classic()

###
# Anderson Darling test
DMSM76_latency<-subset(preference_SPL_latency,SPL=="76")
DMSM82_latency<-subset(preference_SPL_latency,SPL=="82")

ad_test(DMSM76_latency[DMSM76_latency$choice=="DM",]$latency,DMSM76_latency[DMSM76_latency$choice=="SM",]$latency,nboots=5000)
ad_test(DMSM82_latency[DMSM82_latency$choice=="DM",]$latency,DMSM82_latency[DMSM82_latency$choice=="SM",]$latency)
ad_test(preference_SPL_latency[preference_SPL_latency$choice=="DM",]$latency,preference_SPL_latency[preference_SPL_latency$choice=="SM",]$latency)