############################# Complete Analyses and graphs 
############################# Article F-D calls in Parids 
############################# SALIS Ambre, 2023, France 


################################ Packages #####
library(plyr)
library(dplyr)
library(tidyr)
library(ggplot2)
library(Rmisc)
library(cowplot)
library(mgcv)
library(lme4)
library(RVAideMemoire)
library(emmeans)
library(MASS)
library(MuMIn)
library(rsq)
library(glmmTMB)
library(DHARMa)
library(factoextra)
library(cluster)
library(nlme)
################################ DATA ######
#### One example of the raw data used to obtain the SYNTAX dataframe: 
RAW<-read.table("ExempleRaw_BT1.txt", h=T)

#### The acoustic data : 
ACOUSTIC<-read.table("Acoustic.txt", h=T)

#### The data of the calls : 
SYNTAX<-read.table("Syntax.txt", h=T)

#### The raw data of simulation:
RAW_SIMUL<-read.table("SIMUL_TOTAL.txt", h=T)

#### The complete final data: 
FINAL<-read.table("Recordings.txt", h=T)



#####
#####


################################ Acoustic variations between species ####

#####
######################### Graphs #### 
#####
############### Bandwidth ####
ACOUSTIC$NOTE <- factor(ACOUSTIC$NOTE, levels = c("FME", "D"))

tgc <- summarySE(ACOUSTIC, measurevar="band_mean", groupvars=c("SP", "NOTE"))
tgc

p<-ggplot(ACOUSTIC, aes(x=SP, y=band_mean, color = NOTE)) +
  geom_violin(aes(x=SP,y=band_mean, fill = NOTE, color=NA),alpha=0.25, trim=T, position = position_dodge(width = 0.3))+
  geom_errorbar(data= tgc,size=1.1, aes(ymin=band_mean-ci, ymax=band_mean+ci, group=NOTE), colour="grey30", width=.05, position = position_dodge(.3))+
  geom_point(data= tgc,size=4, alpha=0.9, position = position_dodge(width = 0.3)) +
  scale_y_log10() + theme_minimal()
p1<-p+scale_fill_manual(values=c("seagreen", "violetred4"))+scale_color_manual(values=c("violetred4", "seagreen"))+
  ggtitle("Bandwidth") +
  xlab("Species") + ylab("Bandwidth (Hz)")+theme(plot.title = element_text(hjust = 0.5))
p1<-p1+ theme(
  plot.title = element_text(size=20, face="bold"),
  axis.title.x = element_text(size=20, face="bold"),
  axis.title.y = element_text(size=20, face="bold"), 
  axis.text = element_text(size=20, face = "bold"),
  legend.title = element_text(size = 20, face="bold"),
  legend.text = element_text(size = 20, face="bold"))




############### Peak Frequency ####

tgc <- summarySE(ACOUSTIC, measurevar="peakfreq_mean", groupvars=c("SP", "NOTE"))
tgc



p<-ggplot(ACOUSTIC, aes(x=SP, y=peakfreq_mean, color = NOTE)) +
  geom_violin(aes(x=SP,y=peakfreq_mean, fill = NOTE, color=NA),alpha=0.25, trim=T, position = position_dodge(width = 0.3))+
  geom_errorbar(data= tgc,size=1.1, aes(ymin=peakfreq_mean-ci, ymax=peakfreq_mean+ci, group=NOTE), colour="grey30", width=.05, position = position_dodge(.3))+
  geom_point(data= tgc,size=4, alpha=0.9, position = position_dodge(width = 0.3)) +
  scale_y_log10() + theme_minimal()
p1<-p+scale_fill_manual(values=c("seagreen", "violetred4"))+scale_color_manual(values=c("violetred4", "seagreen"))+
  ggtitle("Peak Frequency") +
  xlab("Species") + ylab("Peak Frequency (Hz)")+theme(plot.title = element_text(hjust = 0.5))
p2<-p1+ theme(
  plot.title = element_text(size=20, face="bold"),
  axis.title.x = element_text(size=20, face="bold"),
  axis.title.y = element_text(size=20, face="bold"), 
  axis.text = element_text(size=20, face = "bold"),
  legend.title = element_text(size = 20, face="bold"),
  legend.text = element_text(size = 20, face="bold"))

############### Duration #### 

tgc <- summarySE(ACOUSTIC, measurevar="Duration", groupvars=c("SP", "NOTE"))
tgc



p<-ggplot(ACOUSTIC, aes(x=SP, y=Duration, color = NOTE)) +
  geom_violin(aes(x=SP, y=Duration, fill = NOTE, color=NA), alpha=0.25, trim=T, position = position_dodge(width = 0.3))+
  geom_errorbar(data= tgc,size=1.1, aes(ymin=Duration-ci, ymax=Duration+ci, group=NOTE), colour="grey30", width=.05, position = position_dodge(.3))+
  geom_point(data= tgc,size=4,alpha=0.9, position = position_dodge(width = 0.3)) +
  scale_y_log10() + theme_minimal()
p3<-p+scale_fill_manual(values=c("seagreen", "violetred4"))+scale_color_manual(values=c("violetred4", "seagreen"))+
  ggtitle("Duration") +
  xlab("Species") + ylab("Duration (sec)")+theme(plot.title = element_text(hjust = 0.5))
p3<-p3+ theme(
  plot.title = element_text(size=20, face="bold"),
  axis.title.x = element_text(size=20, face="bold"),
  axis.title.y = element_text(size=20, face="bold"), 
  axis.text = element_text(size=20, face = "bold"),
  legend.title = element_text(size = 20, face="bold"),
  legend.text = element_text(size = 20, face="bold"))





############### Modulation #### 
ACOUSTIC$modulation<-abs(ACOUSTIC$peakfreq_end- ACOUSTIC$PeakFreq_start)
tgc <- summarySE(ACOUSTIC, measurevar="modulation", groupvars=c("SP", "NOTE"))
tgc



p<-ggplot(ACOUSTIC, aes(x=SP, y=modulation, color = NOTE)) +
  geom_violin(aes(x=SP, y=modulation, fill = NOTE, color=NA),alpha=0.25, trim=T, position = position_dodge(width = 0.3))+
  geom_errorbar(data= tgc,size=1.1, aes(ymin=modulation-ci, ymax=modulation+ci, group=NOTE), colour="grey30", width=.05, position = position_dodge(.3))+
  geom_point(data= tgc,size=4,alpha=0.9, position = position_dodge(width = 0.3)) +
  theme_minimal()
p4<-p+scale_fill_manual(values=c("seagreen", "violetred4"))+scale_color_manual(values=c("violetred4", "seagreen"))+
  ggtitle("Modulation") +
  xlab("Species") + ylab("Modulation (Hz)")+theme(plot.title = element_text(hjust = 0.5))
p4<-p4+ theme(
  plot.title = element_text(size=20, face="bold"),
  axis.title.x = element_text(size=20, face="bold"),
  axis.title.y = element_text(size=20, face="bold"), 
  axis.text = element_text(size=20, face = "bold"),
  legend.title = element_text(size = 20, face="bold"),
  legend.text = element_text(size = 20, face="bold"))



# plot final
plot_grid(p1, p2, p3, p4, labels=c("A", "B", "C", "D"), ncol = 2, nrow = 2, label_size = 30)

#####
######################### Models ####

# 1 modèle per sp :
AA_GT<-ACOUSTIC[ACOUSTIC$SP=="Parusmajor",]
AA_BT<-ACOUSTIC[ACOUSTIC$SP=="Paruscaeruleus",]
AA_CT<-ACOUSTIC[ACOUSTIC$SP=="Parusater",]
AA_TUTI<-ACOUSTIC[ACOUSTIC$SP=="Baeolophusbicolor",]
AA_CACH<-ACOUSTIC[ACOUSTIC$SP=="Paruscarolinensis",]
AA_NONN<-ACOUSTIC[ACOUSTIC$SP=="Paruspalustris",]
AA_HUPP<-ACOUSTIC[ACOUSTIC$SP=="Paruscristatus",]


#####
############### Bandwidth ####

##### GT #####

model1 <- glmmTMB(band_mean ~ (1|File) + NOTE ,
                  data = AA_GT,
                  family = truncated_nbinom1)
summary(testResiduals(model1))



summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")



#### BT #### ok


model1 <- glmmTMB(band_mean ~ (1|File) + NOTE ,
                  data = AA_BT,
                  family = truncated_nbinom1)
summary(testResiduals(model1))



summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")



#### CT #### 

model1 <- glmmTMB(band_mean ~ (1|File) + NOTE ,
                  data = AA_CT,
                  family = truncated_nbinom2)
summary(testResiduals(model1))





summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")





#### NONN #### ok
model1 <- glmmTMB(band_mean ~ (1|File) + NOTE ,
                  data = AA_NONN,
                  family = truncated_nbinom2)
summary(testResiduals(model1))


summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")


#### HUPP ####
model1 <- glmmTMB(band_mean ~ (1|File) + NOTE ,
                  data = AA_HUPP,
                  family = truncated_nbinom1)
summary(testResiduals(model1))



summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")



#### CACH #### OK
model1 <- glmmTMB(band_mean ~ (1|File) + NOTE ,
                  data = AA_CACH,
                  family = truncated_nbinom1)
summary(testResiduals(model1))


summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")


#### TUTI ####


model1 <- glmmTMB(band_mean ~ (1|File) + NOTE ,
                  data = AA_TUTI,
                  family = truncated_nbinom1)
summary(testResiduals(model1))


summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")













############### Peak frequency ####

#### GT ####
model1 <- glmmTMB(peakfreq_mean ~ (1|File) + NOTE ,
                  data = AA_GT,
                  family = truncated_nbinom2)

summary(testResiduals(model1))


summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")




#### BT #### 

model1 <- glmmTMB(peakfreq_mean ~ (1|File) + NOTE ,
                  data = AA_BT,
                  family = truncated_nbinom1)
summary(testResiduals(model1))

summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")


#### CT ####  


model1 <- glmmTMB(peakfreq_mean ~ (1|File) + NOTE ,
                  data = AA_CT,
                  family = truncated_nbinom1)
summary(testResiduals(model1))

summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")


#### NONN ####

model1 <- glmmTMB(peakfreq_mean ~ (1|File) + NOTE ,
                  data = AA_NONN,
                  family = truncated_nbinom1)
summary(testResiduals(model1))

summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")



#### HUPP ####


model1 <- glmmTMB(peakfreq_mean ~ (1|File) + NOTE ,
                  data = AA_HUPP,
                  family = truncated_nbinom1)
summary(testResiduals(model1))

summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")



#### CACH #### 


model1 <- glmmTMB(peakfreq_mean ~ (1|File) + NOTE ,
                  data = AA_CACH,
                  family = truncated_nbinom1)
summary(testResiduals(model1))

summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")


#### TUTI ####
model1 <- glmmTMB(peakfreq_mean ~ (1|File) + NOTE ,
                  data = AA_TUTI,
                  family = truncated_nbinom1)
summary(testResiduals(model1))

summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")






############### Duration#### 

#### GT ####


model1<-lmer(log(AA_GT$Duration) ~ AA_GT$NOTE + (1|AA_GT$File))

plot(model1) # check residuals
plot(random.effects(model1)) # check des random effects
summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")




#### BT ####
model1<-lmer(log(AA_BT$Duration*100) ~ AA_BT$NOTE + (1|AA_BT$File))

plotresid(model1) # check residuals
plot(random.effects(model1)) # check des random effects
summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")


#### CT ####


model1<-lmer((AA_CT$Duration*100) ~ AA_CT$NOTE + (1|AA_CT$File))


plotresid(model1) # check residuals
plot(random.effects(model1)) # check des random effects
r.squaredGLMM(model1)
summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")


#### NONN ####

model1<-lmer(log(AA_NONN$Duration*100) ~ AA_NONN$NOTE + (1|AA_NONN$File))


plotresid(model1) # check residuals
plot(random.effects(model1)) # check des random effects

r.squaredGLMM(model1)
summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")


#### HUPP ####

model1<-lmer(log(AA_HUPP$Duration*100) ~ AA_HUPP$NOTE + (1|AA_HUPP$File))


plotresid(model1) # check residuals
plot(random.effects(model1)) # check des random effects
r.squaredGLMM(model1)
summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")


#### CACH ####
model1<-lmer((AA_CACH$Duration*100) ~ AA_CACH$NOTE + (1|AA_CACH$File))


plotresid(model1) # check residuals
plot(random.effects(model1)) # check des random effects
r.squaredGLMM(model1)
summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")


#### TUTI ####


model1<-lmer(log(AA_TUTI$Duration*100) ~ AA_TUTI$NOTE + (1|AA_TUTI$File))

plotresid(model1) # check residuals
plot(random.effects(model1)) # check des random effects
r.squaredGLMM(model1)
summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")




############### Modulation ####

#### GT ####

model1 <- glmmTMB(modulation ~ (1|File) + NOTE ,
                 data = AA_GT,
                 ziformula = ~1,
                 family = nbinom2) # nbinom2 gives the best AIcc
summary(testResiduals(model1))

summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")


#### BT ####
model1 <- glmmTMB(modulation ~ (1|File) + NOTE ,
                  data = AA_BT,
                  ziformula = ~1,
                  family = nbinom1) # nbinom1 gives the best AIcc
summary(testResiduals(model1))

summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")


#### CT ####
model1 <- glmmTMB(modulation ~ (1|File) + NOTE ,
                  data = AA_CT,
                  ziformula = ~1,
                  family = nbinom2) # nbinom2 gives the best AIcc

summary(testResiduals(model1))

summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")


#### NONN ####
model1 <- glmmTMB(modulation ~ (1|File) + NOTE ,
                  data = AA_NONN,
                  ziformula = ~1,
                  family = nbinom1) # nbinom1 gives the best AIcc
summary(testResiduals(model1))
summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")


#### HUPP ####
model1 <- glmmTMB(modulation ~ (1|File) + NOTE ,
                  data = AA_HUPP,
                  ziformula = ~1,
                  family = nbinom2) # nbinom2 gives the best AIcc

summary(testResiduals(model1))
summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")


#### CACH ####
model1 <- glmmTMB(modulation ~ (1|File) + NOTE ,
                  data = AA_CACH,
                  ziformula = ~1,
                  family = nbinom2) # nbinom2 gives the best AIcc
summary(testResiduals(model1))

summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")


#### TUTI ####

model1 <- glmmTMB(modulation ~ (1|File) + NOTE ,
                  data = AA_TUTI,
                  ziformula = ~1,
                  family = nbinom2) # nbinom2 gives the best AIcc
summary(testResiduals(model1))


summary(model1)
emmeans(model1, list(pairwise ~ NOTE), adjust = "tukey", type= "response")









#####
################################# Clustering F from D notes #####



#####
################ GT  ####

# Preparing a dataframe with only the four acoustic variables
AA_2<-scale(AA_GT[, c(6,9,10,11)])


# Testing which clustering linkage method is the most efficient:
m <- c( "average", "single", "complete", "ward")
names(m) <- c( "average", "single", "complete", "ward")

ac <- function(x) {
  agnes(AA_2, method = x)$ac
}

sapply(m, ac) # WARD is the best 


# clustering :
clust <- agnes(AA_2, method = "ward")

#produce dendrogram
pltree(clust, cex = 0.6, hang = -1, main = "Dendrogram")
plot(agnes(AA_2, method = "ward"))

#calculate gap statistic for each number of clusters
gap_stat <- clusGap(AA_2, FUN = hcut, nstart = 25, K.max = 15, B = 100)

#produce plot of clusters vs. gap statistic
fviz_gap_stat(gap_stat) # 10 is the best



#compute distance matrix
d <- dist(AA_2, method = "euclidean")

#perform hierarchical clustering using Ward's method
final_clust <- hclust(d, method = "ward.D2" )

#cut the dendrogram into 2 clusters: 
groups <- cutree(final_clust, k=2)

#find number of observations in each cluster
table(groups)

final_data <- cbind(AA_GT, cluster = groups)
table(final_data$cluster, final_data$NOTE) # check if our classification is congruent with the clustering method


#cut the dendrogram into 10 clusters: 
groups <- cutree(final_clust, k=10)

#find number of observations in each cluster
table(groups)
final_data <- cbind(AA_GT, cluster = groups)
table(final_data$cluster, final_data$NOTE)




################ BLUE TIT ####
AA_2<-scale(AA_BT[, c(6,9,10,11)])


# Testing which clustering linkage method is the most efficient:
m <- c( "average", "single", "complete", "ward")
names(m) <- c( "average", "single", "complete", "ward")

ac <- function(x) {
  agnes(AA_2, method = x)$ac
}

sapply(m, ac) # WARD is the best 


# clustering :
clust <- agnes(AA_2, method = "ward")

#produce dendrogram
pltree(clust, cex = 0.6, hang = -1, main = "Dendrogram")
plot(agnes(AA_2, method = "ward"))

#calculate gap statistic for each number of clusters
gap_stat <- clusGap(AA_2, FUN = hcut, nstart = 25, K.max = 15, B = 100)

#produce plot of clusters vs. gap statistic
fviz_gap_stat(gap_stat) # 8 is the best



#compute distance matrix
d <- dist(AA_2, method = "euclidean")

#perform hierarchical clustering using Ward's method
final_clust <- hclust(d, method = "ward.D2" )

#cut the dendrogram into 2 clusters: 
groups <- cutree(final_clust, k=2)

#find number of observations in each cluster
table(groups)

final_data <- cbind(AA_BT, cluster = groups)
table(final_data$cluster, final_data$NOTE) # check if our classification is congruent with the clustering method


#cut the dendrogram into 10 clusters: 
groups <- cutree(final_clust, k=8)

#find number of observations in each cluster
table(groups)
final_data <- cbind(AA_BT, cluster = groups)
table(final_data$cluster, final_data$NOTE)

################ COAL TIT ####

AA_2<-scale(AA_CT[, c(6,9,10,11)])


# Testing which clustering linkage method is the most efficient:
m <- c( "average", "single", "complete", "ward")
names(m) <- c( "average", "single", "complete", "ward")

ac <- function(x) {
  agnes(AA_2, method = x)$ac
}

sapply(m, ac) # WARD is the best 


# clustering :
clust <- agnes(AA_2, method = "ward")

#produce dendrogram
pltree(clust, cex = 0.6, hang = -1, main = "Dendrogram")
plot(agnes(AA_2, method = "ward"))

#calculate gap statistic for each number of clusters
gap_stat <- clusGap(AA_2, FUN = hcut, nstart =25, K.max = 15, B = 100)

#produce plot of clusters vs. gap statistic
fviz_gap_stat(gap_stat) # 9 is the best



#compute distance matrix
d <- dist(AA_2, method = "euclidean")

#perform hierarchical clustering using Ward's method
final_clust <- hclust(d, method = "ward.D2" )

#cut the dendrogram into 2 clusters: 
groups <- cutree(final_clust, k=2)

#find number of observations in each cluster
table(groups)

final_data <- cbind(AA_CT, cluster = groups)
table(final_data$cluster, final_data$NOTE) # check if our classification is congruent with the clustering method


#cut the dendrogram into 9 clusters: 
groups <- cutree(final_clust, k=9)

#find number of observations in each cluster
table(groups)
final_data <- cbind(AA_CT, cluster = groups)
table(final_data$cluster, final_data$NOTE)




################ NONNETTE ####
AA_2<-scale(AA_NONN[, c(6,9,10,11)])


# Testing which clustering linkage method is the most efficient:
m <- c( "average", "single", "complete", "ward")
names(m) <- c( "average", "single", "complete", "ward")

ac <- function(x) {
  agnes(AA_2, method = x)$ac
}

sapply(m, ac) # WARD is the best 


# clustering :
clust <- agnes(AA_2, method = "ward")

#produce dendrogram
pltree(clust, cex = 0.6, hang = -1, main = "Dendrogram")
plot(agnes(AA_2, method = "ward"))

#calculate gap statistic for each number of clusters
gap_stat <- clusGap(AA_2, FUN = hcut, nstart = 25, K.max = 15, B = 100)

#produce plot of clusters vs. gap statistic
fviz_gap_stat(gap_stat) # 5 is the best



#compute distance matrix
d <- dist(AA_2, method = "euclidean")

#perform hierarchical clustering using Ward's method
final_clust <- hclust(d, method = "ward.D2" )

#cut the dendrogram into 2 clusters: 
groups <- cutree(final_clust, k=2)

#find number of observations in each cluster
table(groups)

final_data <- cbind(AA_NONN, cluster = groups)
table(final_data$cluster, final_data$NOTE) # check if our classification is congruent with the clustering method


#cut the dendrogram into 5 clusters: 
groups <- cutree(final_clust, k=5)

#find number of observations in each cluster
table(groups)
final_data <- cbind(AA_NONN, cluster = groups)
table(final_data$cluster, final_data$NOTE)




################ HUPP   ####

AA_2<-scale(AA_HUPP[, c(6,9,10,11)])


# Testing which clustering linkage method is the most efficient:
m <- c( "average", "single", "complete", "ward")
names(m) <- c( "average", "single", "complete", "ward")

ac <- function(x) {
  agnes(AA_2, method = x)$ac
}

sapply(m, ac) # WARD is the best 


# clustering :
clust <- agnes(AA_2, method = "ward")

#produce dendrogram
pltree(clust, cex = 0.6, hang = -1, main = "Dendrogram")
plot(agnes(AA_2, method = "ward"))

#calculate gap statistic for each number of clusters
gap_stat <- clusGap(AA_2, FUN = hcut, nstart = 25, K.max = 15, B = 100)

#produce plot of clusters vs. gap statistic
fviz_gap_stat(gap_stat) # 8 is the best



#compute distance matrix
d <- dist(AA_2, method = "euclidean")

#perform hierarchical clustering using Ward's method
final_clust <- hclust(d, method = "ward.D2" )

#cut the dendrogram into 2 clusters: 
groups <- cutree(final_clust, k=2)

#find number of observations in each cluster
table(groups)

final_data <- cbind(AA_HUPP, cluster = groups)
table(final_data$cluster, final_data$NOTE) # check if our classification is congruent with the clustering method


#cut the dendrogram into 8 clusters: 
groups <- cutree(final_clust, k=8)

#find number of observations in each cluster
table(groups)
final_data <- cbind(AA_HUPP, cluster = groups)
table(final_data$cluster, final_data$NOTE)

################ CACH ####
AA_2<-scale(AA_CACH[, c(6,9,10,11)])


# Testing which clustering linkage method is the most efficient:
m <- c( "average", "single", "complete", "ward")
names(m) <- c( "average", "single", "complete", "ward")

ac <- function(x) {
  agnes(AA_2, method = x)$ac
}

sapply(m, ac) # WARD is the best 


# clustering :
clust <- agnes(AA_2, method = "ward")

#produce dendrogram
pltree(clust, cex = 0.6, hang = -1, main = "Dendrogram")
plot(agnes(AA_2, method = "ward"))

#calculate gap statistic for each number of clusters
gap_stat <- clusGap(AA_2, FUN = hcut, nstart = 25, K.max = 15, B = 100)

#produce plot of clusters vs. gap statistic
fviz_gap_stat(gap_stat) # 8 is the best



#compute distance matrix
d <- dist(AA_2, method = "euclidean")

#perform hierarchical clustering using Ward's method
final_clust <- hclust(d, method = "ward.D2" )

#cut the dendrogram into 2 clusters: 
groups <- cutree(final_clust, k=2)

#find number of observations in each cluster
table(groups)

final_data <- cbind(AA_CACH, cluster = groups)
table(final_data$cluster, final_data$NOTE) # check if our classification is congruent with the clustering method


#cut the dendrogram into 8 clusters: 
groups <- cutree(final_clust, k=8)

#find number of observations in each cluster
table(groups)
final_data <- cbind(AA_CACH, cluster = groups)
table(final_data$cluster, final_data$NOTE)



################ TUTI ####

AA_2<-scale(AA_TUTI[, c(6,9,10,11)])


# Testing which clustering linkage method is the most efficient:
m <- c( "average", "single", "complete", "ward")
names(m) <- c( "average", "single", "complete", "ward")

ac <- function(x) {
  agnes(AA_2, method = x)$ac
}

sapply(m, ac) # WARD is the best 


# clustering :
clust <- agnes(AA_2, method = "ward")

#produce dendrogram
pltree(clust, cex = 0.6, hang = -1, main = "Dendrogram")
plot(agnes(AA_2, method = "ward"))

#calculate gap statistic for each number of clusters
gap_stat <- clusGap(AA_2, FUN = hcut, nstart = 25, K.max = 10, B = 100)

#produce plot of clusters vs. gap statistic
fviz_gap_stat(gap_stat) # 6 is the best



#compute distance matrix
d <- dist(AA_2, method = "euclidean")

#perform hierarchical clustering using Ward's method
final_clust <- hclust(d, method = "ward.D2" )

#cut the dendrogram into 2 clusters: 
groups <- cutree(final_clust, k=2)

#find number of observations in each cluster
table(groups)

final_data <- cbind(AA_TUTI, cluster = groups)

table(final_data$cluster, final_data$NOTE) # check if our classification is congruent with the clustering method


#cut the dendrogram into 6 clusters: 
groups <- cutree(final_clust, k=6)

#find number of observations in each cluster
table(groups)
final_data <- cbind(AA_TUTI, cluster = groups)
table(final_data$cluster, final_data$NOTE)




#####

################################# Analysis of syntax ####
#####
########## try the script only once to be sure  ######
# Tableau pour machine learning 


d <- read.table("ExempleRaw_BT1.txt")




#Changement nom des colonnes
colnames(d) <- c("file","notes", "start","end","score")

# Création colonne avec temps de chaque note 
d <- data.frame(d, lag=d$end-d$start)
hist(d$lag[d$notes=="S"], breaks = 50) # déterminer espace entre cris 



# Détermination de quelles lignes sont des silences plus petits que 0.15 sec et suppression
todelete <- which(d$note=="S" & d$lag<0.15)
d <- d[-todelete,]


# On crée une colonne qui sépare les clusters de notes 
d <- data.frame(d, seqid=-1)
current.seqid <- 1
for (i in 1:nrow(d)){
  d$seqid[i] <- current.seqid
  if(d$note[i]=="S") current.seqid <- current.seqid+1
}

# On supprime définitivement les Silences restants
issilence <- which(d$note=="S")
d <- d[-issilence,]
d
# On crée une liste qui sépare nos clusters 
splittedd <- split(d, d$seqid)

# on regroupe et compte les différents types de clusters 
A<- sapply(splittedd, function(seq){
  paste(seq$note,collapse="")
})


# Crée un tableau avec les clusters col1, les fréquences col2: 
Tableau<-data.frame(A)
colnames(Tableau) <- c("Clusters")

# Créer une nouvelle colonne qui compte le nombre de lettres dans la colonne cluster

library(stringr)
Tableau$NbFME<-str_count(Tableau$Clusters, regex("FME"))
Tableau$NbD<-str_count(Tableau$Clusters, regex("D"))
Tableau$Nbtotal<-Tableau$NbFME+Tableau$NbD

# Donner la première lettre des clusters #
Tableau$First<-gsub('\\b(\\pL)\\pL{0,}|.','\\U\\1',Tableau$Clusters,perl = TRUE)
Tableau

# ajouter colonne avec le nb de transitions A-D ou D-A : 

k=1
i=1
Transii=0
library(stringi)

for (k in 1:nrow(Tableau)) { 
  
  if(grepl("D", Tableau[k,]$Clusters, fixed = TRUE) == "TRUE"){
    Dpos<-stri_locate_all(pattern = 'D', Tableau[k,]$Clusters, fixed = TRUE)
    Dtable<-data.frame(Dpos)
    Dtable<-data.frame(Dtable, Note="D") } 
  else { Dtable<-data.frame(start=NULL, end=NULL)}
  
  if(grepl("FME", Tableau[k,]$Clusters, fixed = TRUE)== "TRUE"){
    Apos<-stri_locate_all(pattern = 'FME', Tableau[k,]$Clusters, fixed = TRUE)
    Atable<-data.frame(Apos)
    Atable<-data.frame(Atable, Note="FME")  }
  else { Atable<-data.frame(start=NULL, end=NULL)}
  
  Combine<- rbind(Atable, Dtable)
  Combine<- Combine[order(Combine[,2]), ]
  
  if(nrow(Combine) != 1 ) {
    for (i in 1:(nrow(Combine)-1)) 
    { if(Combine[i,]$Note != Combine[i+1,]$Note) Transii<-Transii+1 } 
  }else {Transii<- 0} 
  
  
  
  Tableau$Transition[k]<-Transii
  Transii<-0}  


# On rajoute le nom du fichier 
Tableau$Fichier<-"1"



# On rajoute le nom de l'espèce
Tableau$Espece<-"BT"

nomtableau<-paste0("BT", "1")

Tableau
assign(nomtableau, Tableau)






#####
########################### Simulation of F-D calls ####
#####
####### Show the script for  simulation #####





library(stringi)
library(stringr)




# This function creates N calls (N= number of calls in the natural sample for a given species)
# The size is randomly picked into the natural distribution of length in the natural sample
# The Probability of F and D notes is also picked given the natural probabilities found in our natural data
# There is not rule for the order 
# The function also gives different information on the calls created (first note, size of call, etc...)


F_VS_D<- function(SP, prob_bass, prob_haute,nboot) {
  boot<-numeric(nboot) 
  for(i in 1:nboot) {
    aa<-sample(c("F", "D"), sample(SYNTAX$Nbtotal[SYNTAX$Species==SP], 1),  prob = c(prob_bass, prob_haute), replace = T)
    boot[i] <-paste(aa, collapse = "")
    Simulation<- as.data.frame(boot)
    
    Simulation$NbFME<-str_count(Simulation$boot, regex("F"))
    Simulation$NbD<-str_count(Simulation$boot, regex("D"))
    Simulation$Nbtotal<-Simulation$NbFME+Simulation$NbD
    Simulation$First<-gsub('\\b(\\pL)\\pL{0,}|.','\\U\\1',Simulation$boot,perl = TRUE)
  }
  
  for (k in 1:nrow(Simulation)) { 
    
    if(grepl("D", Simulation[k,]$boot, fixed = TRUE) == "TRUE"){
      Dpos<-stri_locate_all(pattern = 'D', Simulation[k,]$boot, fixed = TRUE)
      Dtable<-data.frame(Dpos)
      Dtable<-data.frame(Dtable, Note="D") } 
    else { Dtable<-data.frame(start=NULL, end=NULL)}
    
    if(grepl("F", Simulation[k,]$boot, fixed = TRUE)== "TRUE"){
      Apos<-stri_locate_all(pattern = 'F', Simulation[k,]$boot, fixed = TRUE)
      Atable<-data.frame(Apos)
      Atable<-data.frame(Atable, Note="F")  }
    else { Atable<-data.frame(start=NULL, end=NULL)}
    
    Combine<- rbind(Atable, Dtable)
    Combine<- Combine[order(Combine[,2]), ]
    
    if(nrow(Combine) != 1 ) {
      for (i in 1:(nrow(Combine)-1)) 
      { if(Combine[i,]$Note != Combine[i+1,]$Note) Transii<-Transii+1 } 
    }else {Transii<- 0} 
    
    
    Simulation$Transition[k]<-Transii
    Transii<-0}
  
  
  Simulation$TransiYN<- ifelse(Simulation$Transition==1, 1,0)
  Simulation$FirstF<- ifelse(Simulation$First=="F", 1,0)
  Simulation$sumFMED<- Simulation$FirstF + Simulation$TransiYN
  Simulation$FMEDYN<- ifelse(Simulation$sumFMED==2, 1,0)
  
  
  return(Simulation)
} 




#### Example :
layout(1)
essai<-F_VS_D("Hupp", 0.13, 0.87, 658)

# confirmation qu'on a bien des cris de la même taille:
hist(essai$Nbtotal)
hist(SYNTAX_hupp$Nbtotal)

# confirmation qu'on a bien le bon % de notes F et D: 
sum(SYNTAX$NbFME[SYNTAX$Species=="Hupp"])/ sum(SYNTAX$Nbtotal[SYNTAX$Species=="Hupp"])
sum(essai$NbFME)/ sum(essai$Nbtotal)




#####
########### For each species, we repeat the simulation and test with wilcoxon test the differences with our data ######




########### Simulations and graphs ######

#### For the following function to work, we need :
# Pourcentage de notes F of the given species:
sum(SYNTAX$NbFME[SYNTAX$Species=="CACH"])/ sum(SYNTAX$Nbtotal[SYNTAX$Species=="CACH"])

#Number of calls of the given species:
sum(SYNTAX$Species=="CACH")


#### Function that produces 1000 dataset with same properties as the natural ones, 
# but with no rules regarding order of the notes: 

coucou3<- function(nboot2 =1000) {
  boot2<-numeric(nboot2) 
  for(i in 1:nboot2) {
    
    Simul1<- as.data.frame(F_VS_D("Nonn", 0.53, 0.47,398)) #NB: for each species we have to adapt this part, with the species name, the proportion of D notes, F notes and the total number of calls
    boot2[i] <- sum(Simul1$FMEDYN)/398
    
    
  }
  
  return(boot2 )
} 


### COAL TIT
COALSIMUL<-coucou3()
head(COALSIMUL)
COALSIMUL<-as.data.frame(COALSIMUL)
sum(COALSIMUL$COALSIMUL>0.0534)/1000


SYNTAX_CT<-SYNTAX[SYNTAX$Species=="CT",]
sum(SYNTAX_CT$FMED)/617

PLOTCT<-ggplot(COALSIMUL, aes(x=COALSIMUL)) +
  geom_histogram(fill="black", position="dodge", bins = 20)+
  geom_vline(xintercept = 0.053,
             linetype="dashed", color="red", size=1 )+
  labs(y= "Frequency", x = "% FME-D calls", title= "P. ater")+
  theme_bw(base_size = 22)



### GREAT TIT
GREATSIMUL<-coucou3()
head(GREATSIMUL)
GREATSIMUL<-as.data.frame(GREATSIMUL)

SYNTAX_GT<-SYNTAX[SYNTAX$Species=="GT",]
sum(SYNTAX_GT$FMED)/334

PLOTGT<-ggplot(GREATSIMUL, aes(x=GREATSIMUL)) +
  geom_histogram(fill="black", position="dodge", bins = 40)+
  geom_vline(xintercept = 0.47,
             linetype="dashed", color="red", size=1)+
  labs(y= "Frequency", x = "% FME-D calls", title= "P. major")+
  theme_bw(base_size = 22)




### BLUE TIT
BLUESIMUL<-coucou3()
head(BLUESIMUL)
BLUESIMUL<-as.data.frame(BLUESIMUL)

SYNTAX_BT<-SYNTAX[SYNTAX$Species=="BT",]
sum(SYNTAX_BT$FMED)/330

PLOTBT<-ggplot(BLUESIMUL, aes(x=BLUESIMUL)) +
  geom_histogram(fill="black", position="dodge", bins = 40)+
  geom_vline(xintercept = 0.25,
             linetype="dashed", color="red", size=1)+
  labs(y= "Frequency", x = "% FME-D calls", title= "C. caeruleus")+
  theme_bw(base_size = 22)



### HUPP TIT
HUPPSIMUL<-coucou3()
head(HUPPSIMUL)
HUPPSIMUL<-as.data.frame(HUPPSIMUL)

SYNTAX_hupp<-SYNTAX[SYNTAX$Species=="Hupp",]
sum(SYNTAX_hupp$FMED)/658

PLOThupp<-ggplot(HUPPSIMUL, aes(x=HUPPSIMUL)) +
  geom_histogram(fill="black", position="dodge", bins = 40)+
  geom_vline(xintercept = 0.31,
             linetype="dashed", color="red", size=1)+
  labs(y= "Frequency", x = "% FME-D calls", title= "L. cristatus")+
  theme_bw(base_size = 22)

### nonn  : 

NONNSIMUL<-coucou3()
head(NONNSIMUL)
NONNSIMUL<-as.data.frame(NONNSIMUL)

SYNTAX_nonn<-SYNTAX[SYNTAX$Species=="Nonn",]
sum(SYNTAX_nonn$FMED)/398

PLOTnonn<-ggplot(NONNSIMUL, aes(x=NONNSIMUL)) +
  geom_histogram(fill="black", position="dodge", bins = 40)+
  geom_vline(xintercept = 0.35,
             linetype="dashed", color="red", size=1)+
  labs(y= "Frequency", x = "% FME-D calls", title= "P. palustris")+
  theme_bw(base_size = 22)


### TUTI

NONNTUTI<-coucou3()
head(NONNTUTI)
NONNTUTI<-as.data.frame(NONNTUTI)

SYNTAX_TUTI<-SYNTAX[SYNTAX$Species=="TUTI",]
sum(SYNTAX_TUTI$FMED)/359

PLOTTUTI<-ggplot(NONNTUTI, aes(x=NONNTUTI)) +
  geom_histogram(fill="black", position="dodge", bins = 30)+
  geom_vline(xintercept = 0.23,
             linetype="dashed", color="red", size=1)+
  labs(y= "Frequency", x = "% FME-D calls", title= "B. bicolor")+
  theme_bw(base_size = 22)



#### CACH

SIMULCACH<-coucou3()
head(SIMULCACH)
SIMULCACH<-as.data.frame(SIMULCACH)

SYNTAX_CACH<-SYNTAX[SYNTAX$Species=="CACH",]
sum(SYNTAX_CACH$FMED)/263

PLOTCACH<-ggplot(SIMULCACH, aes(x=SIMULCACH)) +
  geom_histogram(fill="black", position="dodge", bins = 30)+
  geom_vline(xintercept = 0.65,
             linetype="dashed", color="red", size=1)+  
  labs(y= "Frequency", x = "% FME-D calls", title= "P. carolinensis")+
  theme_bw(base_size = 22)


  ### Graph final 

plot_grid(PLOTTUTI, PLOTCT, PLOTBT,PLOTCACH, PLOThupp,PLOTGT, PLOTnonn,labels=c("A", "B", "C", "D", "E", "F", "G"), ncol = 2, nrow = 4, label_size = 30)






#####
########################### Number of F and D notes in F-D calls #####
coucou<-read.table("SYNTAX_var.txt", h=T)
coucou$ID<-as.factor(coucou$ID)
coucouGT<-coucou[coucou$Espece=="GT",]
coucouCT<-coucou[coucou$Espece=="CT",]
coucouBT<-coucou[coucou$Espece=="BT",]
coucouNONN<-coucou[coucou$Espece=="Nonn",]
coucouTUTI<-coucou[coucou$Espece=="TUTI",]
coucouHUPP<-coucou[coucou$Espece=="Hupp",]
coucouCACH<-coucou[coucou$Espece=="CACH",]


############ GRAPH ##############



coucou$Note2 = factor(coucou$Note, levels = c("F", "D"), ordered = TRUE)
coucou$Espece = factor(coucou$Espece, levels = c("TUTI", "CT", "BT", "CACH", "Hupp", "GT", "Nonn"), ordered = TRUE)


p<-ggplot(coucou, aes(x = Note2, y = Nb)) + 
  geom_line(aes(group = ID), color = "grey") + 
  geom_point(size = 2, color= "gray") +
  geom_boxplot(aes(fill = Note2), alpha = 0.8, fill=c("seagreen", "violetred4","seagreen", "violetred4","seagreen", "violetred4","seagreen", "violetred4","seagreen", "violetred4","seagreen", "violetred4","seagreen", "violetred4"))+
  facet_grid(~ Espece)
p +xlab("Type of Note")+ylab("Number of Notes")+  theme_minimal()+
  theme(strip.background = element_rect(colour="black", fill="white", size=1, linetype="solid"), 
        axis.title.x = element_text(color="black", size=22),
        axis.title.y = element_text(color="black", size=22), 
        legend.text = element_text(color = "black", size = 22),
        axis.text =element_text(color = "black", size = 22),
        strip.text.x = element_text(size = 22))






######## MODELS ######



# GT

model1<-glmmTMB(Nb ~ Note + (1|Fichier)+ (1|ID), zi= ~Note, data=coucouGT, family=truncated_nbinom2,  control = glmmTMBControl(optimizer = optim, optArgs = list(method="BFGS")))
summary(testResiduals(model1))

summary(model1)
emmeans(model1, list(pairwise ~ Note), adjust = "tukey", type= "response")

# BT

model1<-glmmTMB(Nb ~ Note + (1|Fichier)+ (1|ID), zi= ~Note, data=coucouBT, family=truncated_nbinom2)
summary(testResiduals(model1))


summary(model1)
emmeans(model1, list(pairwise ~ Note), adjust = "tukey", type= "response")




# CT

model1<-glmmTMB(Nb ~ Note + (1|Fichier)+ (1|ID), zi= ~Note, data=coucouCT, family=truncated_nbinom2)
summary(testResiduals(model1))


summary(model1)
emmeans(model1, list(pairwise ~ Note), adjust = "tukey", type= "response")


# hupp

model1<-glmmTMB(Nb ~ Note + (1|Fichier)+ (1|ID), zi= ~Note , data=coucouHUPP, family=truncated_nbinom1,  control = glmmTMBControl(optimizer = optim, optArgs = list(method="BFGS")))
summary(testResiduals(model1))


summary(model1)
emmeans(model1, list(pairwise ~ Note), adjust = "tukey", type= "response")


# nonn
model1<-glmmTMB(Nb ~ Note + (1|Fichier)+ (1|ID), zi=~ Note, data=coucouHUPP, family=truncated_nbinom1,  control = glmmTMBControl(optimizer = optim, optArgs = list(method="BFGS")))
summary(testResiduals(model1))

summary(model1)
emmeans(model1, list(pairwise ~ Note), adjust = "tukey", type= "response")


# tuti
model1<-glmmTMB(Nb ~ Note + (1|Fichier)+ (1|ID), zi=~ Note, data=coucouTUTI, family=truncated_nbinom1,  control = glmmTMBControl(optimizer = optim, optArgs = list(method="BFGS")))
summary(testResiduals(model1))

summary(model1)
emmeans(model1, list(pairwise ~ Note), adjust = "tukey", type= "response")

# cach
model1<-glmmTMB(Nb ~ Note + (1|Fichier)+ (1|ID), zi=~ Note, data=coucouCACH, family=truncated_nbinom1,  control = glmmTMBControl(optimizer = optim, optArgs = list(method="BFGS")))
summary(testResiduals(model1))

summary(model1)
emmeans(model1, list(pairwise ~ Note), adjust = "tukey", type= "response")