## Simone Antichi 27/05/2023##

## This is the script to obtain the results of the manuscript submitted to PeerJ Carlón-Beltrán et al.2023 : Whistle characterization of long-beaked common dolphin (Delphinus delphis bairdii) in La Paz Bay, Gulf of California

library(readxl)
Supplemental_data_S2 <- read_excel("C:/Users/lenovo/Desktop/UABCS/PeerJ delphinus paper second resubmission/Supplemental data S2.xlsx", 
                                   sheet = "Whistle parameters")

dataset = Supplemental_data_S2

## descriptive statistics ##


sprintf(min(dataset$`Duration (s)`), fmt = '%#.2f')
sprintf(max(dataset$`Duration (s)`), fmt = '%#.2f')
sprintf(mean(dataset$`Duration (s)`), fmt = '%#.2f')
sprintf(sd(dataset$`Duration (s)`), fmt = '%#.2f')

sprintf(min(dataset$`Maximum frequency (kHz)`), fmt = '%#.2f')
sprintf(max(dataset$`Maximum frequency (kHz)`), fmt = '%#.2f')
sprintf(mean(dataset$`Maximum frequency (kHz)`), fmt = '%#.2f')
sprintf(sd(dataset$`Maximum frequency (kHz)`), fmt = '%#.2f')

sprintf(min(dataset$`Minimum frequency (kHz)`), fmt = '%#.2f')
sprintf(max(dataset$`Minimum frequency (kHz)`), fmt = '%#.2f')
sprintf(mean(dataset$`Minimum frequency (kHz)`), fmt = '%#.2f')
sprintf(sd(dataset$`Minimum frequency (kHz)`), fmt = '%#.2f')

sprintf(min(dataset$`Frequency range (kHz)`), fmt = '%#.2f')
sprintf(max(dataset$`Frequency range (kHz)`), fmt = '%#.2f')
sprintf(mean(dataset$`Frequency range (kHz)`), fmt = '%#.2f')
sprintf(sd(dataset$`Frequency range (kHz)`), fmt = '%#.2f')

sprintf(min(dataset$`Starting frequency (kHz)`), fmt = '%#.2f')
sprintf(max(dataset$`Starting frequency (kHz)`), fmt = '%#.2f')
sprintf(mean(dataset$`Starting frequency (kHz)`), fmt = '%#.2f')
sprintf(sd(dataset$`Starting frequency (kHz)`), fmt = '%#.2f')

sprintf(min(dataset$`Ending frequency (kHz)`), fmt = '%#.2f')
sprintf(max(dataset$`Ending frequency (kHz)`), fmt = '%#.2f')
sprintf(mean(dataset$`Ending frequency (kHz)`), fmt = '%#.2f')
sprintf(sd(dataset$`Ending frequency (kHz)`), fmt = '%#.2f')

sprintf(min(dataset$`Peak frequency (kHz)`), fmt = '%#.2f')
sprintf(max(dataset$`Peak frequency (kHz)`), fmt = '%#.2f')
sprintf(mean(dataset$`Peak frequency (kHz)`), fmt = '%#.2f')
sprintf(sd(dataset$`Peak frequency (kHz)`), fmt = '%#.2f')

sprintf(min(dataset$`Inflection points`), fmt = '%#.2f')
sprintf(max(dataset$`Inflection points`), fmt = '%#.2f')
sprintf(mean(dataset$`Inflection points`), fmt = '%#.2f')
sprintf(sd(dataset$`Inflection points`), fmt = '%#.2f')

sprintf(min(dataset$Steps), fmt = '%#.2f')
sprintf(max(dataset$Steps), fmt = '%#.2f')
sprintf(mean(dataset$Steps), fmt = '%#.2f')
sprintf(sd(dataset$Steps), fmt = '%#.2f')



## histogram ##

values <- c(34.44, 28.89, 11.85, 10, 7.41, 7.41)
group <- c("upsweep","concave","downsweep","sine","convex","constant frequency") 
barplot(values, col = "#99CCFF", names.arg = group, xlab = "", ylab = "",ylim = c(0,50), xaxt = "n", yaxt = "n") 
axis(2, seq(0,50,10), cex.axis=1.4)
mtext("Frequency (%)", side=2, line=2.5, cex=2)
mtext("Whistle contour", side=3, line=-42, at=4, cex=2)
mtext("upsweep", side=3, line=-39.5, at=0.7, cex=1.8)
mtext("concave", side=3, line=-39.5, at=1.8, cex=1.8)
mtext("downsweep", side=3, line=-39.5, at=3.2, cex=1.8)
mtext("sine", side=3, line=-39.5, at=4.3, cex=1.8)
mtext("convex", side=3, line=-39.5, at=5.5, cex=1.8)
mtext("constant frequency", side=3, line=-39.5, at=6.8, cex=1.8)
abline(h=c(0),col="black")

### chi-suared test ##
contours = table(dataset$Contour)
contours = as.data.frame(contours)

frequency_contours = contours$Freq
chisq.test(contours$Freq)


### comparison with Oswald et al., 2007 ###

shapiro.test(dataset$`Duration (s)`)
shapiro.test(dataset$`Maximum frequency (kHz)`)
shapiro.test(dataset$`Minimum frequency (kHz)`)
shapiro.test(dataset$`Frequency range (kHz)`)
shapiro.test(dataset$`Starting frequency (kHz)`)
shapiro.test(dataset$`Ending frequency (kHz)`)
shapiro.test(dataset$`Peak frequency (kHz)`)
shapiro.test(dataset$`Inflection points`)
shapiro.test(dataset$Steps)


wilcox.test(dataset$`Duration (s)`, mu = 0.62, alternative = "less")
wilcox.test(dataset$`Starting frequency (kHz)`, mu = 10.87, alternative = "less")
wilcox.test(dataset$`Ending frequency (kHz)`, mu = 14.46, alternative = "less")
wilcox.test(dataset$`Minimum frequency (kHz)`, mu = 8.48, alternative = "less")
wilcox.test(dataset$`Maximum frequency (kHz)`, mu = 16.21, alternative = "less")
wilcox.test(dataset$`Inflection points`, mu = 1.59, alternative = "less")


### PERMANOVA ###
library(tidyverse)
library(vegan)
library(readxl)
library(devtools)
library(pairwiseAdonis)


Supplemental_data_S2 <- read_excel("C:/Users/lenovo/Desktop/UABCS/PeerJ delphinus paper second resubmission/Supplemental data S2.xlsx", 
                                                           sheet = "Permanova", na = "NA")
data  = Supplemental_data_S2


adonis2(log(data [,1:9] + 1) ~ data$Individual_name, method = "euclidean", center = TRUE, scale = TRUE)
pairwise.adonis(log(data [,1:9] + 1), data$Individual_name, p.adjust.m = "bonferroni")