###data management

combined<-read.csv("combined.csv")

#remove data with 'raining' and 'heavy rain'
#to check whether missing out any other raining hours
combined$bout<-as.numeric(ifelse(combined$remarks=="raining","",combined$bout))
combined$bout<-as.numeric(ifelse(combined$remarks=="heavy rain","",combined$bout))
combined$bout<-as.numeric(ifelse(combined$remarks=="Gymnogryllus too loud","",combined$bout))
combined$bout<-as.numeric(ifelse(combined$remarks=="raining+Phyllomimus","",combined$bout))

#remove randoms with completely no call
combined<-subset(combined,random!="Hindhede27-2")
combined<-subset(combined,random!="Hindhede27-5")
combined<-subset(combined,random!="Labrador89-5")
combined<-subset(combined,random!="Sentosa89-1")

#select field data for only day 1
combined<-subset(combined,day!="2")

#remove lab data without calls
combined<-subset(combined,remarks!="REPEATED")


#make factors
combined$locality<-as.factor(combined$locality)
combined$random<-as.factor(paste(combined$random,combined$date))
combined$date<-as.factor(combined$date)
combined$type<-as.factor(combined$type)
summary(combined)


#temperature statistics
mean(subset(combined,type=="field")$temp)
sd(subset(combined,type=="field")$temp)

mean(subset(combined,type=="lab")$temp)
sd(subset(combined,type=="lab")$temp)


#centering/ standardising bout and temp
fcombined<-subset(combined,type=="field")
fcombined$stdbout<-scale(fcombined$bout)
fcombined$stdtemp<-scale(fcombined$temp)
hist(fcombined$stdtemp)

lcombined<-subset(combined,type=="lab")
lcombined$stdbout<-scale(lcombined$bout)
lcombined$stdtemp<-scale(lcombined$temp)
hist(lcombined$stdtemp)

combined$stdbout<-rbind(fcombined$stdbout,lcombined$stdbout)
combined$stdbout<-scale(combined$stdbout)

combined$temp<-rbind(fcombined$stdtemp,lcombined$stdtemp)
combined$temp<-scale(combined$temp)

head(combined)




###model 1: circadian rhythm
library(mgcv)

#test between lab and field
#using locality only as random effect
m2<-gam(stdbout~type+ 	#test that intercept is different
s(hour,by=type,bs="cc")+ 	#test that the curve is different from 0 for different data type
s(type,hour,bs="fs")+ 	#test that the curve is different from each other
s(locality,bs="re")+	#random effect
scale(temp),	#effect of temp
data=combined,method="REML")

summary(m2)

gam.check(m2)

plot.gam(m2,pages=1)


#using random as random effect
m3<-gam(stdbout~type+ 	#test that intercept is different
s(hour,by=type,bs="cc")+ 	#test that the curve is different from 0 for different data type
s(type,hour,bs="fs")+ 	#test that the curve is different from each other
s(random,bs="re")+	#random effect
scale(temp),	#effect of temp
data=combined,method="REML")

summary(m3)

gam.check(m3)

plot.gam(m3,pages=1)


#m3 has higher R2 and better diagnostic than m2

par(mar = c(4.5,4.5,1,1))
plot(stdbout~hour,data=combined, type="n",
xlab="Time of the day",ylab="Standardised number of echemes",
xlim=c(0,23),ylim=c(-0.5,1.5),cex.lab=1.2)
polygon(x=c(-1,6.85,6.85,-1),y=c(-1.5,-1.5,51.5,51.5),col="#CFCFDE")  
polygon(x=c(18.95,24,24,18.95),y=c(-1.5,-1.5,51.5,51.5),col="#CFCFDE")
polygon(x=c(18.95,6.85,6.85,18.95),y=c(-1.5,-1.5,51.5,51.5),col="white",border="#FDE1C4")

new.hour<-seq(0,23,length.out=408)

newdata1<- with(combined, expand.grid(hour=new.hour,type="field",temp=mean(temp,na.rm=T),random=""))
new<-predict.gam(m3,se.fit=T,re.form=NA,newdata1,type="response") 
polygon(c(new.hour,rev(new.hour)),
    c(new$fit-(1.96*new$se.fit),rev(new$fit+(1.96*new$se.fit))), 
    col = "darkolivegreen1", border = NA)
lines(new$fit~new.hour,col="darkolivegreen4",lwd=5)



newdata1<- with(combined, expand.grid(hour=new.hour,type="lab",temp=mean(temp,na.rm=T),random=""))
new<-predict.gam(m3,se.fit=T,re.form=NA,newdata1,type="response") 
polygon(c(new.hour,rev(new.hour)),
    c(new$fit-(1.96*new$se.fit),rev(new$fit+(1.96*new$se.fit))), 
    col = rgb(0,0,1,0.3), border = NA)
lines(new$fit~new.hour,col="darkblue",lwd=5)



legend("topleft",legend=c("Field data","Laboratory data"),
lwd=5,lty=1,text.col="black",cex=1,
col=c("darkolivegreen4","darkblue"),bty="n")




#for statistics of each dataset
m3<-gam(bout~ 	#test that intercept is different
s(hour,bs="cc")+ 
s(random,bs="re")+	#random effect
scale(temp),	#effect of temp
family=quasipoisson,
data=subset(combined,type=="field"),method="REML")

summary(m3)

new.hour<-seq(0,23,length.out=408)
newdata2<- with(combined, expand.grid(hour=new.hour,type="field",temp=mean(temp,na.rm=T),random=""))

pre<-predict.gam(m3,se.fit=T,re.form=NA,newdata2,type="response") ###change model
data<-as.data.frame(cbind(call=round(pre$fit,0),hour=round(newdata1$hour,0)))
head(data)

aggregate(call~hour,data=data,FUN=mean)
aggregate(call~hour,data=data,FUN=sd)


m4<-gam(bout~ 	#test that intercept is different
s(hour,bs="cc")+ 
s(random,bs="re")+	#random effect
scale(temp),	#effect of temp
family=quasipoisson,
data=subset(combined,type=="lab"),method="REML")

summary(m4)

newdata3<- with(combined, expand.grid(hour=new.hour,type="lab",temp=mean(temp,na.rm=T),random=""))

pre1<-predict.gam(m4,se.fit=T,re.form=NA,newdata3,type="response") ###change model
data1<-as.data.frame(cbind(call=round(pre1$fit,0),hour=round(newdata1$hour,0)))
data1

aggregate(call~hour,data=data1,FUN=mean)
aggregate(call~hour,data=data1,FUN=sd)
























dawn<-read.csv("dawn.csv")

#remove data with no call at all
dawn<-subset(dawn,dawn$remarks!="no call")
dawn<-subset(dawn,dawn$remarks!="raining")
dawn<-subset(dawn,dawn$day=="1")

library(lubridate)
dawn$time<-hour(hms(dawn$time))+(minute(hms(dawn$time))/60)

dawn$locality<-as.factor(dawn$locality)
dawn$random<-as.factor(paste(dawn$random,dawn$date))
dawn$date<-as.factor(dawn$date)

summary(dawn)
unique(dawn$random)


#centering/ standardising bout and temp
fdawn<-subset(dawn,type=="field")
fdawn$stdbout<-scale(fdawn$bout)
fdawn$stdtemp<-scale(fdawn$temp)
hist(fdawn$stdtemp)

ldawn<-subset(dawn,type=="lab")
ldawn$stdbout<-scale(ldawn$bout)
ldawn$stdtemp<-scale(ldawn$temp)
hist(ldawn$stdtemp)

dawn$stdbout<-rbind(fdawn$stdbout,ldawn$stdbout)
dawn$stdbout<-scale(dawn$stdbout)[,1]

dawn$temp<-rbind(fdawn$stdtemp,ldawn$stdtemp)
dawn$temp<-scale(dawn$temp)[,1]

dawn$type<-as.factor(dawn$type)

head(dawn)



library(mgcv)

m2<-gam(bout~ 	#test that intercept is different
s(time)+ 
s(random,bs="re")+	#random effect
scale(temp),	#effect of temp
family=quasipoisson,
data=dawn,method="REML")

summary(m2)
gam.check(m2)




#test between lab and field
#using random as random effect
m3<-gam(stdbout~type+ 	#test that intercept is different
s(time,by=type,bs="fs")+ 	#test that the curve is different from 0 for different data type
s(type,time,bs="fs")+ 	#test that the curve is different from each other
s(random,bs="re")+	#random effect
scale(temp),	#effect of temp
data=dawn,method="REML")

summary(m3)
gam.check(m3)
plot.gam(m3,pages=1)


plot(stdbout~time,data=dawn, type="n",xlab="Time of the day",ylab="Standardised number of echemes",
ylim=c(-2,3),xlim=c(6,7.6))

cicada<-read.csv("cicada.csv")

cicada$sunrise<-hour(hms(cicada$sunrise))+(minute(hms(cicada$sunrise))/60)
cicada$Cstart<-hour(hms(cicada$Cstart))+(minute(hms(cicada$Cstart))/60)
cicada$C.start<-hour(hms(cicada$C.start))+(minute(hms(cicada$C.start))/60)
cicada$C.end<-hour(hms(cicada$C.end))+(minute(hms(cicada$C.end))/60)

polygon(x=c(5.5,mean(cicada$sunrise,na.rm=T),mean(cicada$sunrise,na.rm=T),5.5),
y=c(-3,-3,32,32),
col="#CFCFDE", border = NA)
polygon(x=c(mean(cicada$sunrise,na.rm=T),7.7,7.7,mean(cicada$sunrise,na.rm=T)),
y=c(-3,-3,32,32),
col="white", border = NA)


new.time<-seq(6,7.6,length.out=696)
newdata1<- with(dawn, expand.grid(time=new.time,temp=mean(temp,na.rm=T),random="",type="field"))
new<-predict.gam(m3,se.fit=T,re.form=NA,newdata1,type="response") 
polygon(c(new.time,rev(new.time)),
    c(new$fit-(1.96*new$se.fit),rev(new$fit+(1.96*new$se.fit))), 
    col="darkolivegreen1", border = NA)
lines(new$fit~new.time,col="darkolivegreen4",lwd=5,lty=1)


new.time<-seq(6,7.6,length.out=696)
newdata1<- with(dawn, expand.grid(time=new.time,temp=mean(temp,na.rm=T),random="",type="lab"))
new<-predict.gam(m3,se.fit=T,re.form=NA,newdata1,type="response") 
polygon(c(new.time,rev(new.time)),
    c(new$fit-(1.96*new$se.fit),rev(new$fit+(1.96*new$se.fit))), 
    col=rgb(0,0,1,0.3), border = NA)
lines(new$fit~new.time,col="darkblue",lwd=5,lty=1)

abline(v=mean(cicada$Cstart,na.rm=T),col="#2b2220",lwd=2.5,lty=2)
arrows(mean(cicada$C.start,na.rm=T),3,mean(cicada$C.end,na.rm=T),3,col="#72401d",lwd=2.5,code=3,length=0.1)
abline(v=mean(cicada$C.start,na.rm=T),col="#72401d",lwd=2.5)
abline(v=mean(cicada$C.end,na.rm=T),col="#72401d",lwd=2.5)

















#ANH and dominant frequency
freq<-read.csv("frequency.csv")
head(freq)
freq<-subset(freq,freq$day=="1")
freq$call<-as.factor(freq$call)

kruskal.test(domf~call,data=freq)

pairwise.wilcox.test(freq$domf,freq$call,p.adj="bonf")

library(ggplot2)
ggplot(freq, aes(x=call, y=domf,color=call)) + 
  geom_boxplot(color=c("black","#72401d","darkolivegreen4"))+
geom_jitter(shape=16, position=position_jitter(0.2))+
scale_x_discrete(limits=c("purana", "cicada^", "lebinthus"))+
labs(x="Species", y = "Dominant frequency (kHz)")+
scale_color_manual(values=c("#72401d","darkolivegreen4","black"),guide="none")+
theme_classic()


mean(subset(freq,freq$call=="cicada^")$domf,na.rm=T)
sd(subset(freq,freq$call=="cicada^")$domf,na.rm=T)

mean(subset(freq,freq$call=="purana")$domf,na.rm=T)
sd(subset(freq,freq$call=="purana")$domf,na.rm=T)





#cicada start and end statistics
cicada<-read.csv("cicada.csv")
head(cicada)

cicada$C.delta<-hour(hms(cicada$C.delta))+(minute(hms(cicada$C.delta))/60)
mean(cicada$C.delta,na.rm=T)*60
sd(cicada$C.delta,na.rm=T)*60

cicada$Cstart<-hour(hms(cicada$Cstart))+(minute(hms(cicada$Cstart))/60)
sd(cicada$Cstart,na.rm=T)*60

cicada$C.start<-hour(hms(cicada$C.start))+(minute(hms(cicada$C.start))/60)
sd(cicada$C.start,na.rm=T)*60

cicada$C.end<-hour(hms(cicada$C.end))+(minute(hms(cicada$C.end))/60)
sd(cicada$C.end,na.rm=T)*60


















# Lebinthus dawn chorus figure 2
# Hannah ter Hofstede
# Last updated Jul. 4, 2022

# clear workspace
rm(list=ls(all=TRUE)) 
# Set working directory
setwd("~/Documents/Hannah/Publications/2022 Tan Lebinthus CR/Figures/")

# Open libraries "seewave" and "tuneR"
library("seewave")
library("tuneR")
# library("viridis")   

# Load recording
Sen <- readWave("Sentosa_20210916_063402.wav")

# reduced sampling rate for 1 hour spectro
Sen.r <- resamp(Sen, 96000, 44100, output="Wave")

# 1 hour figure: in 15 minutes segments
tiff("1 hour spectro 1-15 min.tiff")
spectro(Sen.r, wl = 1024, grid = FALSE, tlim=c(0,900), collevels=seq(-60,-5,5))
dev.off()

tiff("1 hour spectro 15-30 min.tiff")
spectro(Sen.r, wl = 1024, grid = FALSE, tlim=c(900,1800), collevels=seq(-60,-5,5))
dev.off()

tiff("1 hour spectro 30-45 min.tiff")
spectro(Sen.r, wl = 1024, grid = FALSE, tlim=c(1800,2700), collevels=seq(-60,-5,5))
dev.off()

tiff("1 hour spectro.tiff 45-57 min")
spectro(Sen.r, wl = 1024, grid = FALSE, tlim=c(2700,3597), collevels=seq(-70,-5,5))
dev.off()


# Nocturnal crickets
tiff("nocturnal cricket spectro wl 512 ovlp 85.tiff", width = 1800, height = 1200, res = 200)
spectro(Sen, wl = 512, grid = FALSE, tlim=c(410.4,420.4), ovlp = 85,
        flim=c(0,22), collevels=seq(-60,0,5), trel = FALSE)
dev.off()

# ^ Cicada
tiff("Cicada ^ spectro wl 512 ovlp 85.tiff", width = 1800, height = 1200, res = 200)
spectro(Sen, wl = 512, grid = FALSE, tlim=c(1867.4,1897.4), ovlp = 85,
        flim=c(0,22), collevels=seq(-60,0,5), trel = FALSE)
dev.off()

# Purana nr. usnani cicada
tiff("Cicada Purana spectro wl 512 ovlp 85.tiff", width = 1800, height = 1200, res = 200)
spectro(Sen, wl = 512, grid = FALSE, tlim=c(2743,2746), ovlp = 85,
        flim=c(0,22), collevels=seq(-60,0,5), trel = FALSE)
dev.off()

# Lebinthus
tiff("Lebinthus 1 wl 512 ovlp 85.tiff", width = 1800, height = 1200, res = 200)
spectro(Sen, wl = 512, grid = FALSE, tlim=c(3533.4,3536.4), ovlp = 85,
        flim=c(0,22), collevels=seq(-60,0,5), trel = FALSE)
dev.off()