## Exp4
# Author: Hongyu Li 2020
  
library(lme4)
library(AICcmodavg)
library(tidyr)
library(ggplot2)
library(dplyr)
library(plyr)
library(imputeTS)
library(janitor)
library(Hmisc)
library(texreg)
library(varhandle)
library(knitr)
library(cowplot)
library(MASS)

# Data preparation
raw_num<-read.csv("./Exp4.csv",na.strings = "")%>%
  drop_na(SPL)%>%
  dplyr::select(SPL,stimulus,response,ï..date,duration,frogNum)%>%
  mutate(frogID=paste(ï..date,frogNum,sep=" "))

## process the raw data to get the count for number of females responded
raw<-raw_num%>%
  mutate(SPL=as.factor(SPL))

count<-count(raw, c("stimulus","SPL","response"))%>%
  drop_na()%>%
  spread(response,freq)

TOTAL<-count$YES+count$NO
count<-cbind(count,TOTAL)

count<-binconf(x=count$YES, n=count$TOTAL)%>%
  cbind(count,.)

# response rate regardless of stimulus type
count_global<-count(raw, c("SPL","response"))%>%
  drop_na()%>%
  spread(response,freq)

TOTAL_global<-count_global$YES+count_global$NO
count_global<-cbind(count_global,TOTAL_global)

count_global<-binconf(x=count_global$YES, n=count_global$TOTAL)%>%
  cbind(count_global,.)


# how does response rate change with SPL regardless of stimulus type?
# Fit model and plotting the model
signal_64_num<-subset(raw_num,SPL>=64)
glm64.1<-glmer(response~ 1+(1|frogID), data=signal_64_num, family=binomial(link = "logit"))
glm64.2<-glmer(response~ 1+(1|frogID)+SPL, data=signal_64_num, family=binomial(link = "logit"))
anova(glm64.1,glm64.2) # SPL is significant

fig_SNR<-
  ggplot(data=count_global,aes(x=SPL,y=PointEst))+
  geom_point(data=count_global,position = position_dodge(width = 0.2)) +
  geom_errorbar(data=count_global,aes(ymin=Lower, ymax=Upper,width=0.2),position = position_dodge(width = 0.2))+
  xlab("Signal to noise ratios (dB)") +
  ylab("Response proportions")+
  ylim(0:1)+
  theme_classic()+
  scale_x_discrete(labels=c("No acoustic \n signal","-12","-6","0","6","12"),breaks=c("0","64","70","76","82","88"),limits=c("0","64","70","76","82","88"))+
  theme(legend.position = c(0.13,0.9))+
  theme(legend.title=element_blank())+
  theme(axis.text=element_text(colour="black"))+
  theme(text=element_text(size=16))

## Add the regression line for response rate ~ SNR
# glm64.0<-glm(response~ 1+SPL, data=signal_64_num, family=binomial(link = "logit"))
# plotModel_64<-data.frame(SPL=seq(from = 64, to = 88, by = 0.02))%>%
#   cbind((predict(glm64.0,.,type = "response",se.fit=TRUE)))%>%
#   mutate(PointEst=fit)%>%
#   mutate(Upper=fit+1.96*se.fit,Lower=fit-1.96*se.fit)
# 
# model_64_CI_Upper<-select(plotModel_64,SPL,PointEst=Upper)
# model_64_CI_Lower<-select(plotModel_64,SPL,PointEst=Lower)
# 
# plotModel_64$SPL=as.factor(plotModel_64$SPL)
# model_64_CI_Upper$SPL=as.factor(model_64_CI_Upper$SPL)
# model_64_CI_Lower$SPL=as.factor(model_64_CI_Lower$SPL)
# 
# fig_SNR_model<-
#   ggplot(data=count_global,aes(x=SPL,y=PointEst))+
#   geom_line(data=plotModel_64,aes(group=1))+
#   geom_point(data=count_global,position = position_dodge(width = 0.2)) +
#   geom_errorbar(data=count_global,aes(ymin=Lower, ymax=Upper,width=0.2),position = position_dodge(width = 0.2))+
#   xlab("Signal to noise ratios (dB)") +
#   ylab("Response proportions")+
#   ylim(0:1)+
#   theme_classic()+
#   scale_x_discrete(labels=c("No acoustic \n signal","-12","-6","0","6","12"),breaks=c("0","64","70","76","82","88"),limits=c("0","64","70","76","82","88"))+
#   theme(legend.position = c(0.13,0.9))+
#   theme(legend.title=element_blank())+
#   theme(axis.text=element_text(colour="black"))+
#   theme(text=element_text(size=16))

save_plot("fig_SNR.png",fig_SNR,dpi=300)
save_plot("fig_SNR.eps",fig_SNR)

# Response rate ~ SPL + stimulus type
## Make graphs
count$stimulus<-factor(count$stimulus,levels(count$stimulus)[c(2,3,1)])

signal_present<-subset(count,SPL!=0)
signal_absent<-subset(count,SPL==0)

fig_SNR_stimulus<-
  ggplot(data=signal_present,aes(x=SPL,shape=stimulus,y=PointEst,group=stimulus))+
  geom_point(data=signal_present,position = position_dodge(width = 0.35)) +
  geom_point(data=signal_absent,position = position_dodge(width = 0.35)) +
  geom_errorbar(aes(ymin=Lower, ymax=Upper,width=0.5),position = position_dodge(width = 0.35))+
  geom_errorbar(data=signal_absent,aes(ymin=Lower, ymax=Upper,width=0.5),position = position_dodge(width = 0.35))+
  xlab("Signal to noise ratios (dB)") +
  ylab("Response proportions") +
  ylim(0:1)+
  theme_classic()+
  labs(shape="Stimulus types")+
  scale_shape_discrete(labels = c(expression(NM^{"+"}), expression(SM^{"+"}),expression(DM^{"+"}))) +
  scale_x_discrete(labels=c("No acoustic \n signal","-12","-6","0","6","12"),breaks=c("0","64","70","76","82","88"),limits=c("0","64","70","76","82","88"))+
  theme(legend.position = c(0.13,0.9))+
  theme(legend.title=element_blank())+
  theme(axis.text=element_text(colour="black"))+
  theme(text=element_text(size=16))

save_plot("fig_SNR_stimulus.png",fig_SNR_stimulus,dpi=300)
save_plot("fig_SNR_stimulus.eps",fig_SNR_stimulus)

# Response rates to noise as visual treatment varies
signal_absent_num<-subset(raw_num,SPL==0)
glm.0<-glmer(response~ 1+(1|frogID), data=signal_absent_num, family=binomial(link = "logit"))
glm.1<-glmer(response~ 1+(1|frogID)+stimulus, data=signal_absent_num, family=binomial(link = "logit"))
anova(glm.0,glm.1)
# The response rate to noise is not affected by stimulus types

## Full models (numerical SPL)
glm64.1<-glmer(response~ 1+(1|frogID), data=signal_64_num, family=binomial(link = "logit"))
AICc(glm64.1)
summary(glm64.1)

glm64.2<-glmer(response~ 1+SPL+(1|frogID), data=signal_64_num, family=binomial(link = "logit"))
AICc(glm64.2)
summary(glm64.2)

glm64.3<-glmer(response~ 1+SPL+(1|frogID)+stimulus, data=signal_64_num, family=binomial(link = "logit"))
summary(glm64.3)
AICc(glm64.3)

glm64.4<-glmer(response~ 1+SPL+(1|frogID)+stimulus+stimulus*SPL, data=signal_64_num, family=binomial(link = "logit"))
summary(glm64.4)
AICc(glm64.4)

# Generate a statistics table
texreg::texreg(list(glm64.1,glm64.2,glm64.3,glm64.4),digits=3,star=NULL)

# response latency
latency<-subset(raw_num,response=="YES")%>%
  drop_na(duration)
latency$stimulus<-factor(latency$stimulus,levels=c("call","SM","DM"))

latency$SPL<-as.factor(latency$SPL)

fig_SNR_latency_box<-
  ggplot(data=latency,aes(x=SPL,y=duration,fill=stimulus))+
  geom_boxplot(width=0.5)+
  xlab("Signal to noise ratios (dB)") +
  ylab("Response latency (s)")+
  scale_fill_manual(values=c("grey100","grey88","grey76"),labels = c(expression(NM^{"+"}),expression(SM^{"+"}),expression(DM^{"+"})))+
  scale_x_discrete(labels=c("No acoustic \n signal","-12","-6","0","6","12"),breaks=c("0","64","70","76","82","88"),limits=c("0","64","70","76","82","88"))+
  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"))

save_plot("fig_SNR_latency_box.png",fig_SNR_latency,base_height=3,base_width=6,dpi=300)
save_plot("fig_SNR_latency_box.eps",fig_SNR_latency,base_height=3,base_width=6)

###
#The violin plot of response latency
fig_SNR_latency_violin<-
  ggplot(data=latency,aes(x=SPL,y=duration,fill=stimulus))+
  geom_violin(scale="count")+
  xlab("SPL measured at 50cm") +
  ylab("latency (s)") +
  ggtitle("response latency")

# Combined data
dummy<-rep(1, each=153)
log<-log(latency$duration)
latency<-cbind(latency,dummy,log)
fig_SNR_latency<-
  ggplot(data=latency,aes(x=dummy,y=duration))+
  geom_violin(scale="count")+
  xlab("SPL measured at 50cm") +
  ylab("latency (s)") +
  stat_summary(fun.data="mean_sdl", fun.args=list(mult=1))+
  ggtitle("response latency")

latency_64<-subset(latency,SPL!="0")
latency_64$SPL<-as.numeric(latency_64$SPL)

# Is that Gaussian distribution?
qqnorm(latency_64$duration)

# How about log tranformation?
qqnorm(latency_64$log)

### Density plot for original and predicted values #######################
density <- density(latency_64$duration,bw=20) # returns the density data
distribution<-plot(density,col="red",ylim=c(0,.007))

# Gaussian
fitdistr(latency_64$duration, densfun = "normal")
x <- seq(-100, 400,by=1)
y <- dnorm(x, mean=133.096296, sd=73.049539)
distribution+lines(x, y)

# Lognormal
fitdistr(latency_64$duration, densfun = "log-normal")
x <- seq(-100, 400,by=1)
y <- dlnorm(x, meanlog=4.71377362, sdlog=0.63235653)
distribution+lines(x, y)

# Weibull
x <- seq(0,400,by=1)
fitdist(latency_64$duration,"weibull")
y <- dweibull(x,shape=1.922712,150.486312)
distribution+lines(x, y)

# Weibull model
library(survival)
wei.SNR.0<-survreg(Surv(duration,dummy)~frailty(frogID),latency_64,dist="weibull")
wei.SNR.1<-survreg(Surv(duration,dummy)~SPL+frailty(frogID),latency_64,dist="weibull")
wei.SNR.2<-survreg(Surv(duration,dummy)~stimulus+frailty(frogID),latency_64,dist="weibull")
wei.SNR.3<-survreg(Surv(duration,dummy)~stimulus+SPL+frailty(frogID),latency_64,dist="weibull")

anova(wei.SNR.0,wei.SNR.1) #SPL is not significant
anova(wei.SNR.0,wei.SNR.2) #stimulus is not significant
anova(wei.SNR.0,wei.SNR.3) #SPL*stimulus not significant

library(SurvRegCensCov)
ConvertWeibull(wei.SNR.0,conf.level = 0.95)

WeibullDiag(Surv(duration,dummy)~stimulus,data=latency)
WeibullDiag(Surv(duration,dummy)~SPL,data=latency)