diff --git a/README.md b/README.md
index 2a6ec2616f34c24db69703e6f1fdbbcdc2cf391c..f1f626d37d7168537773e0623c3fa51b96f3edd1 100644
--- a/README.md
+++ b/README.md
@@ -23,9 +23,7 @@ $ pip install name-package
 .matlab files can be executed in matlab or octave.
 
 ## Contact
-Feel free to contact me if you have questions, tips or anything else to say. I'd really appreciate it!
-
-Loïc Lehnhoff - <loic.lehnhoff@gmail.com>
+Please contact [Loïc Lehnhoff](mailto:loic.lehnhoff@gmail.com), [Bastien Mérigot](mailto:bastien.merigot@umontpellier.fr) or [Hervé Glotin](mailto:herve.glotin@lis-lab.fr) for any question related to this repository.
 
 ## Related to
 Results and scripts are presented in details in: Lehnhoff, L.; Glotin, H.; Bernard, S.; Dabin, W.; Le Gall, Y.; Menut, E.; Meheust, E.; Peltier, H.; Pochat, A.; Pochat, K.; Rimaud, T.; Sourget, Q.; Spitz, J.; Van Canneyt, O.; Mérigot, B. Behavioural Responses of Common Dolphins Delphinus delphis to a Bio-Inspired Acoustic Device for Limiting Fishery By-Catch. Sustainability 2022, 14, 13186. https://doi.org/10.3390/su142013186
\ No newline at end of file
diff --git a/Stats/.Rhistory b/Stats/.Rhistory
index 2d7f5bdf7fc481d7600c8b044f236b475791a9cb..711e717faf579bdc3bd57c3e030de797b9010243 100644
--- a/Stats/.Rhistory
+++ b/Stats/.Rhistory
@@ -1,257 +1,3 @@
-legend_title="Fishing net", legend_labs=c("Present", "Absent"))
-# BBPs
-letters_df <- computeLetters(emmeans(mod.bbp, pairwise~fishing_net:acoustic, adjust="tukey"),
-"fishing_net")
-letters_df$acoustic <- computeLetters(emmeans(mod.bbp, pairwise~fishing_net:acoustic, adjust="tukey"),
-"acoustic")$acoustic
-letters_df <- letters_df[, c("acoustic","fishing_net","letter")]
-letters_df$letter <- gsub(" ", "", letters_df$letter)
-barPlot(computeStats(acoustic.dta, fishing_net, BBPs_per_dolphin, two=acoustic),
-NULL, acoustic, fill=fishing_net,
-old_names = c("AV","AV+D","D","D+AP","AP"), ytitle="Mean number of BBPs per dolphin per min",
-new_names = c("BEF","BEF+DUR","DUR", "DUR+AFT", "AFT"),
-xname="Activation sequence", height=c(1.65,1.65,1.72,1.65,1.72,1.65,1.65,1.72,1.72,1.72),
-colours=c("#E69F00","#999999"), size=5,
-legend_title="Fishing net", legend_labs=c("Present", "Absent"))
-# Clicks
-letters_df <- computeLetters(emmeans(mod.cli, pairwise~fishing_net:acoustic, adjust="tukey"),
-"fishing_net")
-letters_df$acoustic <- computeLetters(emmeans(mod.cli, pairwise~fishing_net:acoustic, adjust="tukey"),
-"acoustic")$acoustic
-letters_df <- letters_df[, c("acoustic","fishing_net","letter")]
-letters_df$letter <- gsub(" ", "", letters_df$letter)
-barPlot(computeStats(acoustic.dta, fishing_net, clicks_per_dolphin, two=acoustic),
-NULL, acoustic, fill=fishing_net,
-old_names = c("AV","AV+D","D","D+AP","AP"), ytitle="Mean number of clicks per dolphin per min",
-new_names = c("BEF","BEF+DUR","DUR", "DUR+AFT", "AFT"),
-xname="Activation sequence", height=c(180,180,187,187,180,187,180,180,187,187),
-colours=c("#E69F00","#999999"), size=5,
-legend_title="Fishing net", legend_labs=c("Present", "Absent"))
-#### Behaviour plots ####
-# Whistles
-table <- cld(emmeans(mod.whi, pairwise~behavior, adjust="tukey"), Letters = letters)
-myletters_df <- data.frame(behavior=table$behavior,letter = trimws(table$.group))
-barPlot(computeStats(acoustic.dta, behavior, whistling_time_per_dolphin/n_bins),
-myletters_df, behavior, height=1.5, ytitle="Mean whistling time per dolphin per min (in sec)",
-old_names = c("CHAS", "DEPL", "SOCI"),
-new_names = c("Foraging", "Travelling", "Socialising"),
-xname="Behaviours of dolphins")
-# BBPs
-# real effect measured in model
-table <- cld(emmeans(mod.bbp, pairwise~behavior, adjust="tukey"), Letters = letters)
-myletters_df <- data.frame(acoustic=table$behavior,letter = trimws(table$.group))
-barPlot(computeStats(acoustic.dta, behavior, BBPs_per_dolphin),
-myletters_df, behavior, height=1.2, ytitle="Mean number of BBPs per dolphin per min",
-old_names = c("CHAS", "DEPL", "SOCI"),
-new_names = c("Foraging", "Travelling", "Socialising"),
-xname="Behaviours of dolphins")
-# Clicks
-# no significant effect in click statistical model so all the same letters
-myletters_df <- data.frame(behavior=unique(acoustic.dta$behavior),
-letter = rep("a",length(unique(acoustic.dta$behavior))))
-barPlot(computeStats(acoustic.dta, behavior, clicks_per_dolphin),
-myletters_df,
-behavior, old_names = c("CHAS", "DEPL", "SOCI"),
-new_names = c("Foraging", "Travelling", "Socialising"),
-xname="Behaviours of dolphins", height=150,
-ytitle="Mean number of clicks per dolphin per min")
-#### Nets plots + KW analysis ####
-# Whistles
-#KW test
-kruskal.test(acoustic.dta$whistling_time_per_dolphin ~ acoustic.dta$net)
-# p<0.05 so post-hoc
-kruskalmc(acoustic.dta$whistling_time_per_dolphin, acoustic.dta$net)
-# DIY : letters
-myletters_df <- data.frame(net=c("SSF", "chalut_blanc", "chalut_vert", "tremail", "grand_filet"),
-letter = c("a","ad","bd","cd","a"))
-barPlot(computeStats(acoustic.dta, net, whistling_time_per_dolphin/n_bins),
-NULL,
-net, old_names = c("SSF", "chalut_blanc", "chalut_vert", "tremail", "grand_filet"),
-new_names = c("Absent", "Nylon trawl net", "PE trawl net", "Nylon gill net", "Long nylon gill net"),
-xname="Fishing nets", height=.6,
-ytitle="Mean whistling time per dolphin per min (in sec)")+
-theme(axis.text.x=element_text(size=8.5))
-# BBPs
-#KW test
-kruskal.test(acoustic.dta$BBPs_per_dolphin ~ acoustic.dta$net)
-# p<0.05 so post-hoc
-kruskalmc(acoustic.dta$BBPs_per_dolphin, acoustic.dta$net)
-# DIY : letters
-myletters_df <- data.frame(net=c("SSF", "chalut_blanc", "chalut_vert", "tremail", "grand_filet"),
-letter = c("a","a","a","a","a"))
-barPlot(computeStats(acoustic.dta, net, BBPs_per_dolphin),
-NULL,
-net, old_names = c("SSF", "chalut_blanc", "chalut_vert", "tremail", "grand_filet"),
-new_names = c("Absent", "Nylon trawl net", "PE trawl net", "Nylon gill net", "Long nylon gill net"),
-xname="Fishing nets", height=.8,
-ytitle="Mean number of BBPs per dolphin per min")+
-theme(axis.text.x=element_text(size=8.5))
-# Clicks
-#KW test
-kruskal.test(acoustic.dta$clicks_per_dolphin ~ acoustic.dta$net)
-# p<0.05 so post-hoc
-kruskalmc(acoustic.dta$clicks_per_dolphin, acoustic.dta$net)
-# DIY : letters
-myletters_df <- data.frame(net=c("SSF", "chalut_blanc", "chalut_vert", "tremail", "grand_filet"),
-letter = c("ae","ad","bd","cd","e"))
-barPlot(computeStats(acoustic.dta, net, clicks_per_dolphin),
-NULL,
-net, old_names = c("SSF", "chalut_blanc", "chalut_vert", "tremail", "grand_filet"),
-new_names = c("Absent", "Nylon trawl net", "PE trawl net", "Nylon gill net", "Long nylon gill net"),
-xname="Fishing nets", height=120,
-ytitle="Mean number of clicks per dolphin per min")+
-theme(axis.text.x=element_text(size=8.5))
-#### Beacon plots + KW analysis (letters not shown for readability) ####
-# Whistles
-#KW test
-kruskal.test(acoustic.dta$whistling_time_per_dolphin ~ acoustic.dta$beacon)
-names = computeStats(acoustic.dta, beacon, whistling_time_per_dolphin/n_bins)["beacon"]
-barPlot(computeStats(acoustic.dta, beacon, whistling_time_per_dolphin/n_bins),
-NULL,
-beacon, old_names = unlist(names), new_names = unlist(names),
-xname="Signals from bio-inspired beacon", height=0.9, size=3,
-ytitle="Mean whistling time per dolphin per min (in sec)")+
-theme(axis.text.x=element_text(size=8))+
-scale_x_discrete(guide=guide_axis(n.dodge = 2))
-# NC stands for "Unknown". Corresponding to categories where the beacon was not turned on yet ('BEF')
-# BBPs
-#KW test
-kruskal.test(acoustic.dta$BBPs_per_dolphin ~ acoustic.dta$beacon)
-names = computeStats(acoustic.dta, beacon, BBPs_per_dolphin)["beacon"]
-barPlot(computeStats(acoustic.dta, beacon, BBPs_per_dolphin),
-NULL,
-beacon, old_names = unlist(names), new_names = unlist(names),
-xname="Signals from bio-inspired beacon", height=0.5, size=3,
-ytitle="Mean number of BBPs per dolphin per min")+
-theme(axis.text.x=element_text(size=8))+
-scale_x_discrete(guide=guide_axis(n.dodge = 2))
-# NC stands for "Unknown". Corresponding to categories where the beacon was not turned on yet ('BEF')
-# Clicks
-#KW test
-kruskal.test(acoustic.dta$clicks_per_dolphin ~ acoustic.dta$beacon)
-names = computeStats(acoustic.dta, beacon, clicks_per_dolphin)["beacon"]
-barPlot(computeStats(acoustic.dta, beacon, clicks_per_dolphin),
-NULL,
-beacon, old_names = unlist(names), unlist(names),
-xname="Signals from bio-inspired beacon", height=150, size=3,
-ytitle="Mean number of clicks per dolphin per min")+
-theme(axis.text.x=element_text(size=8))+
-scale_x_discrete(guide=guide_axis(n.dodge = 2))
-# NC stands for "Unknown". Corresponding to categories where the beacon was not turned on yet ('BEF')
-#### Plots by number of dolphins ####
-# Whistles
-numb_stats_w <- computeStats(acoustic.dta, number, total_whistles_duration/n_bins)
-numb_stats_w[is.na(numb_stats_w)] <- 0
-numb_stats_w$number <- as.factor(numb_stats_w$number)
-numb_stats_w %>%
-ggplot(aes(x=number, y=mean, group=1)) +
-geom_errorbar(aes(x=number, ymin=mean-ic, ymax=mean+ic),
-color="red", width=.1, show.legend = FALSE)+
-geom_point() + geom_line() +
-theme_classic() + theme(text=element_text(size=12)) +
-ylab("Mean whistling time per dolphin per min (in sec)")+
-xlab("Number of dolphins in group")
-# BBPs
-numb_stats_b <- computeStats(acoustic.dta, number, number_of_bbp)
-numb_stats_b[is.na(numb_stats_b)] <- 0
-numb_stats_b$number <- as.factor(numb_stats_b$number)
-numb_stats_b %>%
-ggplot(aes(x=number, y=mean, group=1)) +
-geom_errorbar(aes(x=number, ymin=mean-ic, ymax=mean+ic),
-color="red", width=.1, show.legend = FALSE)+
-geom_point() + geom_line() +
-theme_classic() + theme(text=element_text(size=12)) +
-ylab("Number of BBPs per dolphin per min")+
-xlab("Number of dolphins in group")
-# Clicks
-numb_stats_c <- computeStats(acoustic.dta, number, number_of_clicks)
-numb_stats_c[is.na(numb_stats_c)] <- 0
-numb_stats_c$number <- as.factor(numb_stats_c$number)
-numb_stats_c %>%
-ggplot(aes(x=number, y=mean, group=1)) +
-geom_errorbar(aes(x=number, ymin=mean-ic, ymax=mean+ic),
-color="red", width=.1)+
-geom_point() + geom_line() +
-theme_classic() + theme(text=element_text(size=12)) +
-ylab("Mean number of clicks per dolphin per min")+
-xlab("Number of echolocation clicks in group")
-#### Plots by Group ID ####
-# Whistles
-numb_stats_w <- computeStats(acoustic.dta, ID, whistling_time_per_dolphin/n_bins)
-numb_stats_w[is.na(numb_stats_w)] <- 0
-numb_stats_w$ID <- as.factor(numb_stats_w$ID)
-numb_stats_w %>%
-ggplot(aes(x=ID, y=mean, group=1)) +
-geom_errorbar(aes(x=ID, ymin=mean-sd, ymax=mean+ci),
-color="red", width=.1, show.legend = FALSE)+
-geom_point() + scale_x_discrete(guide = guide_axis(n.dodge = 2))+
-theme_light() + theme(text=element_text(size=12)) +
-ylab("Mean whistling time per dolphin per min (in sec)")+
-xlab("ID of dolphins group")
-# BBPs
-numb_stats_b <- computeStats(acoustic.dta, ID, BBPs_per_dolphin)
-numb_stats_b[is.na(numb_stats_b)] <- 0
-numb_stats_b$ID <- as.factor(numb_stats_b$ID)
-numb_stats_b %>%
-ggplot(aes(x=ID, y=mean, group=1)) +
-geom_errorbar(aes(x=ID, ymin=mean-sd, ymax=mean+ci),
-color="red", width=.1, show.legend = FALSE)+
-geom_point() + scale_x_discrete(guide = guide_axis(n.dodge = 2))+
-theme_light() + theme(text=element_text(size=12)) +
-ylab("Number of BBPs per dolphin per min")+
-xlab("ID of dolphins group")
-# Clicks
-numb_stats_c <- computeStats(acoustic.dta, ID, clicks_per_dolphin)
-numb_stats_c[is.na(numb_stats_c)] <- 0
-numb_stats_c$ID <- as.factor(numb_stats_c$ID)
-numb_stats_c %>%
-ggplot(aes(x=ID, y=mean, group=1)) +
-geom_errorbar(aes(x=ID, ymin=mean-sd, ymax=mean+ci),
-color="red", width=.1)+
-geom_point() + scale_x_discrete(guide = guide_axis(n.dodge = 2))+
-theme_light() + theme(text=element_text(size=12)) +
-ylab("Mean number of clicks per dolphin per min")+
-xlab("ID of dolphin group")
-# KW test on IDs
-# whistles (excluding groups "2" because 1 sample and groups because no whistles recorded )
-data_test <- acoustic.dta[acoustic.dta$ID!="2",]
-print( posthocKW(data_test$whistling_time_per_dolphin, data_test$ID))
-print( posthocKW(data_test$BBPs_per_dolphin, data_test$ID))
-print( posthocKW(data_test$clicks_per_dolphin, data_test$ID))
-#### Plots by days ####
-# clicks
-numb_stats_c <- computeStats(acoustic.dta, date, clicks_per_dolphin)
-numb_stats_c %>%
-ggplot(aes(x=date, y=mean, group=1)) +
-geom_errorbar(aes(x=date, ymin=mean-sd, ymax=mean+ci),
-color="black", width=.1)+
-geom_point() + scale_x_discrete(guide = guide_axis(n.dodge=2)) +
-theme_classic() + theme(text=element_text(size=12)) +
-ylab("Mean number of clicks per dolphin per min")+
-xlab("Days of recording")
-# BBPs
-numb_stats_c <- computeStats(acoustic.dta, date, BBPs_per_dolphin)
-numb_stats_c %>%
-ggplot(aes(x=date, y=mean, group=1)) +
-geom_errorbar(aes(x=date, ymin=mean-sd, ymax=mean+ci),
-color="black", width=.1)+
-geom_point() + scale_x_discrete(guide = guide_axis(n.dodge=2)) +
-theme_classic() + theme(text=element_text(size=12)) +
-ylab("Number of BBPs per dolphin per min")+
-xlab("Days of recording")
-# Whistles
-numb_stats_c <- computeStats(acoustic.dta, date, whistling_time_per_dolphin/n_bins)
-numb_stats_c %>%
-ggplot(aes(x=date, y=mean, group=1)) +
-geom_errorbar(aes(x=date, ymin=mean-sd, ymax=mean+ci),
-color="black", width=.1)+
-geom_point() + scale_x_discrete(guide = guide_axis(n.dodge=2)) +
-theme_classic() + theme(text=element_text(size=12)) +
-ylab("Mean whistling time per dolphin per min (in sec)")+
-xlab("Days of recording")
-computeStats(acoustic.dta, acoustic, clicks_per_dolphin)
-108/43.8
-88.2/43.8
 ########################################################################
 #    STATISTICS
 #    Author : Loic LEHNHOFF
@@ -290,223 +36,6 @@ clicks.dta <- clicks.dta[order(clicks.dta$audio_names),]
 # Merge files into 1 dataset
 acoustic.dta <- clicks.dta
 acoustic.dta$number_of_bbp <- bbp.dta$number_of_BBP
-acoustic.dta$total_whistles_duration <- whistles.dta$total_whistles_duration
-rm(whistles.dta, bbp.dta, clicks.dta)
-# add group IDs
-id2020 <- read.table(file=paste0(folder, 'CSV_data/Audio_Data_2020.csv'),
-sep = ',', header=TRUE)[1:396,]
-id2021 <- read.table(file=paste0(folder, 'CSV_data/Audio_Data_2021.csv'),
-sep = ',', header=TRUE)[1:96,]
-id2021$ID <- id2021$ID+max(id2020$ID)
-id2021$Seq <- id2021$Seq+max(id2020$Seq)
-id.dta <- rbind(id2020, id2021)
-id.dta$Fichier.Audio <-str_sub(id.dta$Fichier.Audio, -27, -5)
-acoustic.dta$ID <- rep(-1, 490)
-for (name in acoustic.dta$audio_names){
-acoustic.dta$ID[match(name, acoustic.dta$audio_names)] <- id.dta$ID[match(name, id.dta$Fichier.Audio)]
-}
-acoustic.dta$ID <- as.factor(acoustic.dta$ID)
-rm(id2020, id2021, id.dta)
-# suppress "T" acoustic data (other groups not tested on our variables)
-acoustic.dta <- acoustic.dta[acoustic.dta$acoustic!="T",]
-# shuffle dataframe
-acoustic.dta <- acoustic.dta[sample(1:nrow(acoustic.dta)), ]
-acoustic.dta$acoustic <- factor(acoustic.dta$acoustic)
-#################### DATA INSPECTION  #################################
-# Data description
-names(acoustic.dta)
-# self explenatory except acoustic : correspond to the activation sequence.
-# Look for obvious correlations
-plot(acoustic.dta) # nothing that we can see
-# Look for zero-inflation
-100*sum(acoustic.dta$number_of_clicks == 0)/nrow(acoustic.dta)
-100*sum(acoustic.dta$number_of_bbp == 0)/nrow(acoustic.dta)
-100*sum(acoustic.dta$total_whistles_duration == 0)/nrow(acoustic.dta)
-# 5.8%, 53.7% & 27.1% of our data are zeros. Will have to be dealt with.
-# QUESTION: This study is aimed at understanding if dolphin's acoustic activity
-# is influenced by their behavior, the emission of a pinger or a fishing net.
-# Dependent variables (Y): number_of_clicks, number_of_bbp, total_whistles_duration.
-# Explanatory variables (X): acoustic, fishing_net, behavior, beacon, net, number.
-# What are the H0/ H1 hypotheses ?
-# H0 : No influence of any of the explanatory variables on a dependant one.
-# H1 : Influence of an explanatory variable on a dependent one.
-##################### DATA EXPLORATION ################################
-# Y Outlier detection
-par(mfrow=c(2,3))
-boxplot(acoustic.dta$total_whistles_duration, col='red',
-ylab='total_whistles_duration')
-boxplot(acoustic.dta$number_of_bbp, col='red',
-ylab='number_of_bbp')
-boxplot(acoustic.dta$number_of_clicks, col='red',
-ylab='number_of_clicks')
-dotchart(acoustic.dta$total_whistles_duration, pch=16,
-xlab='total_whistles_duration', col='red')
-dotchart(acoustic.dta$number_of_bbp, pch=16,
-xlab='number_of_bbp', col='red')
-dotchart(acoustic.dta$number_of_clicks, pch=16,
-xlab='number_of_clicks', col='red')
-# Y distribution
-par(mfrow=c(2,3))
-hist(acoustic.dta$total_whistles_duration, col='red', breaks=8,
-xlab='total_whistles_duration', ylab='number')
-hist(acoustic.dta$number_of_bbp, col='red', breaks=8,
-xlab='number_of_bbp', ylab='number')
-hist(acoustic.dta$number_of_clicks, col='red', breaks=8,
-xlab='number_of_clicks', ylab='number')
-qqnorm(acoustic.dta$total_whistles_duration, col='red', pch=16)
-qqline(acoustic.dta$total_whistles_duration)
-qqnorm(acoustic.dta$number_of_bbp, col='red', pch=16)
-qqline(acoustic.dta$number_of_bbp)
-qqnorm(acoustic.dta$number_of_clicks, col='red', pch=16)
-qqline(acoustic.dta$number_of_clicks)
-shapiro.test(acoustic.dta$total_whistles_duration)
-shapiro.test(acoustic.dta$number_of_bbp)
-shapiro.test(acoustic.dta$number_of_clicks)
-# p-values are significant => they do not follow normal distributions
-# will need a transformation or the use of a glim model
-# X Number of individuals per level
-summary(factor(acoustic.dta$acoustic))
-summary(factor(acoustic.dta$fishing_net))
-summary(factor(acoustic.dta$behavior))
-summary(factor(acoustic.dta$beacon))
-summary(factor(acoustic.dta$net))
-table(factor(acoustic.dta$acoustic),factor(acoustic.dta$fishing_net))
-table(factor(acoustic.dta$acoustic),factor(acoustic.dta$behavior))
-table(factor(acoustic.dta$behavior),factor(acoustic.dta$acoustic))
-ftable(factor(acoustic.dta$fishing_net), factor(acoustic.dta$behavior), factor(acoustic.dta$acoustic))
-# => unbalanced, no big deal but will need more work (no orthogonality):
-# Effects can depend on the order of the variables
-# => Beacon and net have modalities with <10 individuals => analysis impossible
-# => They will be treated apart from the rest as they are likely to be biased
-##################### STATISTICAL MODELLING ###########################
-### Model tested
-# GLM: General linear model (residual hypothesis: normality, homoscedasticity, independant)
-# GLIM: Generalized linear model (residual hypothesis: uncorrelated residuals)
-# NB : Negative Binomial model (residual hypothesis: independantM)
-# ZINB: Zero inflated negative binomial model (residual hypothesis: independant)
-# We are using number as an offset (more dolphins => more signals)
-# beacon and net explanatory variables could not be tested in models
-# as they contain information already present in "fishing_net" which is more
-# interesting to keep for our study. They will be treated after
-# (using kruskall-Wallis non-parametric test)
-# fishing_net, behavior and acoustic where tested with their interactions.
-# If a variable is it in a model, it is because it had no significant effect.
-par(mfrow=c(1,1))
-### Model for whistles
-# Residual hypotheses not verified for GLM
-# Overdipsersion when using GLIM (negative binomial)
-# Using ZINB:
-zero.whi <- zeroinfl(total_whistles_duration ~
-acoustic + fishing_net + behavior + offset(log(number)),
-data=acoustic.dta, dist='negbin')
-nb.whi <- glm.nb(total_whistles_duration ~
-acoustic + fishing_net + behavior + offset(log(number)),
-data=acoustic.dta)
-# comparison ZINB VS NB model
-vuong(zero.whi, nb.whi)  #(if p-value<0.05 then first model in comparison is better)
-mod.whi <- zero.whi # => zeroinflated model is indeed better suited
-car::Anova(mod.whi, type=3)
-dwtest(mod.whi) # H0 -> independent if p>0.05 (autocorrelation if p<0.05)
-bptest(mod.whi) # H0 -> homoscedasticity if p<0.05
-mod.whi$df.null/mod.whi$df.residual
-# no dispersion, perfect
-### Model for BBP
-# No normality of residuals for GLM
-# overdispersion with GLIM quasipoisson
-#try with glim NB:
-mod.bbp <- glm.nb(number_of_bbp ~ acoustic + fishing_net + behavior
-+ offset(log(number)),
-data=acoustic.dta)
-car::Anova(mod.bbp, type=3)
-dwtest(mod.bbp) # H0 -> independent if p>0.05 (autocorrelation if p<0.05)
-bptest(mod.bbp) # H0 -> homoscedasticity if p<0.05
-mod.bbp$deviance/mod.bbp$df.residual
-# slight underdispersion, not improved with ZINB so we keep this
-### Model for clicks
-# Using NB model:
-mod.cli <- glm.nb(number_of_clicks ~ acoustic + fishing_net + acoustic:fishing_net + offset(log(number)),
-data=acoustic.dta)
-car::Anova(mod.cli, type=3)
-dwtest(mod.cli) # H0 -> independent if p>0.05 (autocorrelation if p<0.05)
-bptest(mod.cli) # H0 -> homoscedasticity if p<0.05
-mod.cli$deviance/mod.cli$df.residual
-# slight overdispersion. (ZINB does not clearly improve results so we keep this)
-# FYI1: Comparison of combination of explanatory variables between models
-# were compared based on BIC criterion.
-# The model with the lowest BIC was kept (and is the one shown)
-# FYI2: log(number of dolphin per group) does have an effect on data but we have
-# no interest in investigating it, that is why we use it as an offset.
-##################### Boxplots and comparisons #####################
-### Functions to compute stats
-computeLetters <- function(temp, category) {
-test <- multcomp::cld(object = temp$emmeans,
-Letters = letters)
-myletters_df <- data.frame(category=test[,category],
-letter=trimws(test$.group))
-colnames(myletters_df)[1] <- category
-return(myletters_df)
-}
-computeStats <- function(data, category, values, two=NULL, three=NULL) {
-my_sum <- data %>%
-group_by({{category}}, {{two}}, {{three}}) %>%
-summarise(
-n=n(),
-mean=mean({{values}}),
-sd=sd({{values}})
-) %>%
-mutate( se=sd/sqrt(n))  %>%
-mutate( ic=se * qt((1-0.05)/2 + .5, n-1))
-return(my_sum)
-}
-barPlot <- function(dta, signif, category, old_names, new_names, fill=NULL, size=5,
-height, xname="", colours="black", legend_title="", legend_labs="",ytitle=""){
-if (!is.null(signif)){colnames(signif)[1] <- "use"}
-dta %>%
-mutate(use=fct_relevel({{category}}, old_names)) %>%
-ggplot(aes(x=use, y=mean, group={{fill}}, fill={{fill}},color={{fill}})) +
-{if(length(colours)==1)geom_point(color=colours, position=position_dodge(.5))}+
-{if(length(colours)==2)geom_point(position=position_dodge(.5), show.legend = FALSE)}+
-{if(length(colours)==2)scale_color_manual(values=colours, name=legend_title, labels=legend_labs)}+
-scale_x_discrete(breaks=old_names,
-labels=new_names)+
-ylab(ytitle)+
-xlab(xname)+
-theme_classic()+ theme(text=element_text(size=12))+
-{if(!is.null(signif))geom_text(data=signif, aes(label=letter, y=height), size=size,
-colour="black", position=position_dodge(.5))}+
-geom_errorbar(aes(x=use, ymin=mean-ic, ymax=mean+ic), position=position_dodge(.5), width=.1, show.legend = FALSE)
-}
-####Introducing variables averaged per dolphins ####
-# since we introduced an offset, variables can be divided by the number of dolphins
-acoustic.dta$whistling_time_per_dolphin <- acoustic.dta$total_whistles_duration/acoustic.dta$number
-acoustic.dta$BBPs_per_dolphin <- acoustic.dta$number_of_bbp/acoustic.dta$number
-acoustic.dta$clicks_per_dolphin <- acoustic.dta$number_of_clicks/acoustic.dta$number
-#### Fishing net  ####
-#### Behaviour plots ####
-# Whistles
-table <- cld(emmeans(mod.whi, pairwise~behavior, adjust="tukey"), Letters = letters)
-myletters_df <- data.frame(behavior=table$behavior,letter = trimws(table$.group))
-barPlot(computeStats(acoustic.dta, behavior, whistling_time_per_dolphin/n_bins),
-myletters_df, behavior, height=1.5, ytitle="Mean whistling time per dolphin per min (in sec)",
-old_names = c("CHAS", "DEPL", "SOCI"),
-new_names = c("Foraging", "Travelling", "Socialising"),
-xname="Behaviours of dolphins")
-# BBPs
-# real effect measured in model
-table <- cld(emmeans(mod.bbp, pairwise~behavior, adjust="tukey"), Letters = letters)
-myletters_df <- data.frame(acoustic=table$behavior,letter = trimws(table$.group))
-barPlot(computeStats(acoustic.dta, behavior, BBPs_per_dolphin),
-myletters_df, behavior, height=1.2, ytitle="Mean number of BBPs per dolphin per min",
-old_names = c("CHAS", "DEPL", "SOCI"),
-new_names = c("Foraging", "Travelling", "Socialising"),
-xname="Behaviours of dolphins")
-# Clicks
-# no significant effect in click statistical model so all the same letters
-myletters_df <- data.frame(behavior=unique(acoustic.dta$behavior),
-letter = rep("a",length(unique(acoustic.dta$behavior))))
-barPlot(computeStats(acoustic.dta, behavior, clicks_per_dolphin),
-myletters_df,
-behavior, old_names = c("CHAS", "DEPL", "SOCI"),
-new_names = c("Foraging", "Travelling", "Socialising"),
-xname="Behaviours of dolphins", height=150,
-ytitle="Mean number of clicks per dolphin per min")
+acoustic.dta
+acoustic.dta$number_of_bbp
+mean(acoustic.dta$number_of_bbp)