### Exp3
# Author: Hongyu Li

library(dplyr)
library(tidyr)
library(ggplot2)
library(lme4)
library(cowplot)
library(psych)
library(fitdistrplus)
library(survival)
library(eha)
library(MASS)
library(psych)
library(distr)

sam<-read.csv("./Exp3_threshold.csv")%>%
  dplyr::select(frogID,Stimulus,threshold,previous_tests,order,latency)%>%
  drop_na()
sam$Stimulus<-factor(sam$Stimulus,levels=c("Hch","HchSM","HchDM"))

##############################################################################################
### threshold analysis
# violin lot
label<-c(expression(NM^{"+"}),expression(SM^{"+"}),expression(DM^{"+"}))
fig_sam=
  ggplot(data=sam,aes(x=Stimulus,y=threshold)) +
  geom_violin(width=0.5)+
  stat_summary(fun.data="mean_sdl", fun.args=list(mult=1))+
  ylab("Modulation depth thresholds") +
  ylim(0:1)+
  xlab("Stimulus types")+
  theme_classic()+
  theme(text=element_text(size=16))+
  theme(axis.text=element_text(colour="black"))+
  scale_x_discrete(limits=c("Hch","HchSM","HchDM"),labels=label)

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

# SAM threshold analysis, generalized mixed-effected model
sam.summary<-describeBy(sam$threshold, sam$Stimulus,mat=TRUE)
sam.summary

sam.lm0<- lmer(threshold ~ 1+(1|frogID), data=sam,REML = FALSE)
sam.lm1<- lmer(threshold ~ 1+(1|frogID)+order, data=sam,REML = FALSE)
sam.lm2<- lmer(threshold ~ 1+(1|frogID)+Stimulus, data=sam,REML = FALSE)
anova(sam.lm0,sam.lm1)
anova(sam.lm0,sam.lm2)


##############################################################################################
# Latency analysis 
# Violin plot
label<-c(expression(NM^{"+"}),expression(SM^{"+"}),expression(DM^{"+"}))
fig_sam_latency=
  ggplot(data=sam,aes(x=Stimulus,y=latency)) +
  geom_violin(width=0.5)+
  stat_summary(fun.data="mean_sdl", fun.args=list(mult=1))+
  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(limits=c("Hch","HchSM","HchDM"),labels=label)

##############################################################################################
### models #######################################################
# non-parametric KS test, pairwise comparison
sam_latency_wide<-spread(sam, key = Stimulus, value = latency)
ks.test(sam_latency_wide$Hch,sam_latency_wide$HchSM) #d=0.146,p=0.912
ks.test(sam_latency_wide$Hch,sam_latency_wide$HchDM) #d=0.133,p=0.953
ks.test(sam_latency_wide$HchSM,sam_latency_wide$HchDM) #d=0.179,p=0.730

# Violin plot of compiled data. Check the distribution of response latency to decide models to be applied
dummy<-rep(1, each=89)
log<-log(sam$latency)
sam<-cbind(sam,dummy,log)

fig_sam_latency_global=
  ggplot(data=sam,aes(x=dummy,y=latency)) +
  geom_violin(width=0.5)+
  stat_summary(fun.data="mean_sdl", fun.args=list(mult=1))+
  ylab("response latency")

# Gaussian models ################################################
qqnorm(sam$latency)

## Linear regression models for response latency
sam.latency.lm0<- lmer(latency ~ 1+(1|frogID), data=sam,REML = FALSE)
sam.latency.lm1<- lmer(latency ~ 1+(1|frogID)+order, data=sam,REML = FALSE)
sam.latency.lm2<- lmer(latency ~ 1+(1|frogID)+previous_tests, data=sam,REML = FALSE)
anova(sam.latency.lm0,sam.latency.lm1) #order, chi(1)=3.372,p=0.066
anova(sam.latency.lm0,sam.latency.lm2) #previous_tests, chi(1)=1.937,p=0.164
# order should be included in the model to evaluate the significance of stimulus

sam.latency.lm3<- lmer(latency ~ 1+(1|frogID)+order+Stimulus, data=sam,REML = FALSE)
anova(sam.latency.lm2,sam.latency.lm4) #stimulus, chi(2)=0.552,p=0.759
# Stimulus is not significant in explaning response latency

## Log models for response latency ###############################
# qqplot log tranformation
qqnorm(sam$log)

sam.latency.log0<- lmer(log ~ 1+(1|frogID), data=sam,REML = FALSE)
sam.latency.log1<- lmer(log ~ 1+(1|frogID)+order, data=sam,REML = FALSE)
sam.latency.log2<- lmer(log ~ 1+(1|frogID)+previous_tests, data=sam,REML = FALSE)
anova(sam.latency.log0,sam.latency.log1) #order, chi(1)=0.968,p=0.325
anova(sam.latency.log0,sam.latency.log2) #previous_tests, chi(1)=0.164,p=0.685
# neither should be included in the model to evaluate the significance of stimulus

sam.latency.log3<- lmer(latency ~ 1+(1|frogID)+Stimulus, data=sam,REML = FALSE)
anova(sam.latency.log0,sam.latency.log3) #stimulus, chi(2)=0,p=1
# Stimulus is not significant in explaning response latency

# Use survival:surveg to do the same analysis (same method as weibull analysis)
log.sam.0<-survreg(Surv(latency,dummy)~frailty(frogID),sam,dist="lognormal")
log.sam.1<-survreg(Surv(latency,dummy)~Stimulus+frailty(frogID),sam,dist="lognormal")
log.sam.2<-survreg(Surv(latency,dummy)~order+frailty(frogID),sam,dist="lognormal")
log.sam.3<-survreg(Surv(latency,dummy)~previous_tests+frailty(frogID),sam,dist="lognormal")
anova(log.sam.0,log.sam.1)
anova(log.sam.0,log.sam.2)
anova(log.sam.0,log.sam.3)

## Weibull models ################################################
library(survival)
wei.sam.0<-survreg(Surv(latency,dummy)~frailty(frogID),sam,dist="weibull")
wei.sam.1<-survreg(Surv(latency,dummy)~Stimulus+frailty(frogID),sam,dist="weibull")
wei.sam.2<-survreg(Surv(latency,dummy)~order+frailty(frogID),sam,dist="weibull")
wei.sam.3<-survreg(Surv(latency,dummy)~previous_tests+frailty(frogID),sam,dist="weibull")
anova(wei.sam.0,wei.sam.1) #stimulus, chi(2.144)=2.944,p=0.253
anova(wei.sam.0,wei.sam.2) #order, chi(1.267)=1.512,p=0.287
anova(wei.sam.0,wei.sam.3) #previous tests, chi(1.080)=0.165,p=0.716

# Get the hazard parameters for the models
library(SurvRegCensCov)
ConvertWeibull(wei.sam.0,conf.level = 0.95)
ConvertWeibull(wei.sam.1,conf.level = 0.95)

# Check the cumulative hazrd funcion. Should be linear and parallel to each other
WeibullDiag(Surv(latency,dummy)~Stimulus,data=sam)

### Goodness of fit plots ##########################################################################
library(fitdistrplus)
gofstat(fitdist(sam$latency,"norm"))
gofstat(fitdist(sam$latency,"lnorm"))
gofstat(fitdist(sam$latency,"weibull"))

### Hazard function for 3 models####################################
coxreg.sam<-coxreg(Surv(latency,dummy)~1,sam)

# Gaussian
phreg.sam<-phreg(Surv(latency,dummy)~frailty(frogID),sam,dist='norm')
check.dist(coxreg.sam,phreg.sam)

# Lognormal
phreg.sam<-phreg(Surv(latency,dummy)~frailty(frogID),sam,dist='lognormal')
check.dist(coxreg.sam,phreg.sam)

# Weibull
phreg.sam<-phreg(Surv(latency,dummy)~frailty(frogID),sam,dist='weibull')
check.dist(coxreg.sam,phreg.sam)

### Density plot for original and predicted values #######################
density <- density(sam$latency,bw=20) # returns the density data

# Gaussian
fitdistr(sam$latency, densfun = "normal")
x <- seq(-100, 400,by=1)
y <- dnorm(x, mean=132.741573, sd=78.689909)
plot(density,col="red",ylim=c(0,.007))
lines(x, y)

# Lognormal
fitdistr(sam$latency, densfun = "log-normal")
x <- seq(-100, 400,by=1)
y <- dlnorm(x, meanlog=4.67549427, sdlog=0.69757346)
plot(density,col="red",ylim=c(0,.007))
lines(x, y)

# Weibull
x <- seq(0,300,by=1)
fitdist(sam$latency,"weibull")
y <- dweibull(x,shape=1.75403,149.44315)
plot(density,col="red",ylim=c(0,.007))
lines(x, y)

### Survival curve for the 3 models ######################################
#plot kaplan-meier estimate, per stimulus
s <- with(sam,Surv(latency,dummy))
sam_KM <- survfit(s ~ Stimulus,data=sam)
plot(sam_KM)

# Predicted latency
lm.sam.4<-survreg(Surv(latency,dummy)~Stimulus,sam,dist="gaussian")
lines(predict(lm.sam.4, newdata=list(Stimulus="Hch"),type="quantile",p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col="red")
lines(predict(lm.sam.4, newdata=list(Stimulus="HchSM"),type="quantile",p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col="red")
lines(predict(lm.sam.4, newdata=list(Stimulus="HchDM"),type="quantile",p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col="red")

log.sam.4<-survreg(Surv(latency,dummy)~Stimulus,sam,dist="lognormal")
lines(predict(log.sam.4, newdata=list(Stimulus="Hch"),type="quantile",p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col="red")
lines(predict(log.sam.4, newdata=list(Stimulus="HchSM"),type="quantile",p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col="red")
lines(predict(log.sam.4, newdata=list(Stimulus="HchDM"),type="quantile",p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col="red")

wei.sam.4<-survreg(Surv(latency,dummy)~Stimulus,sam,dist="weibull")
lines(predict(wei.sam.4, newdata=list(Stimulus="Hch"),type="quantile",p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col="red")
lines(predict(wei.sam.4, newdata=list(Stimulus="HchSM"),type="quantile",p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col="red")
lines(predict(wei.sam.4, newdata=list(Stimulus="HchDM"),type="quantile",p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col="red")