diff --git a/.gitignore b/.gitignore
index 29f7900d3395f00378168dbcc04e6d0a0cd2c3ff..61ac5e194db331708d2d46e12ae0ea12c4d44392 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,2 +1,3 @@
 
 *.sav
+Stats/.Rhistory
diff --git a/Stats/.Rhistory b/Stats/.Rhistory
index be7778ea363a655ee8271b09c294e72990357753..363f710925e38fb4086501d3a4a42c1e6636ad78 100644
--- a/Stats/.Rhistory
+++ b/Stats/.Rhistory
@@ -1,68 +1,3 @@
-library(MASS)
-library(lmtest)
-library(multcomp)
-library(emmeans)
-library(dplyr)        # "%>%" function
-library(forcats)      # "fct_relevel" function
-library(stringr)      # "gsub" function
-library(rcompanion)   # "fullPTable" function
-library(multcompView) # "multcompLetters" function
-library(ggplot2)
-library(pgirmess)
-#library(tidyquant)    # geom_ma() if rolling average needed
-################# DATASET IMPORTS #####################################
-folder <- './../'
-whistles.dta <-read.table(file=paste0(folder,
-'Whistles/Evaluation/whistles_durations.csv'),
-sep = ',', header=TRUE)
-whistles.dta <- whistles.dta[order(whistles.dta$audio_names),]
-bbp.dta <-read.table(file=paste0(folder,
-'BBPs/Results/16-06-22_14h00_number_of_BBP.csv'),
-sep = ',', header=TRUE)
-bbp.dta <- bbp.dta[order(bbp.dta$audio_names),]
-clicks.dta <-read.table(file=paste0(folder,
-'Clicks/Results/projection_updated_number_of_clicks_02052022.csv'), #number_of_clicks_02052022.csv
-sep = ',', header=TRUE)
-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)
-# 3.6%, 53.7% & 24.7% 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 bytheir 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.
@@ -510,3 +445,68 @@ print( posthocKW(data_test$BBPs_per_dolphin, data_test$ID))
 print( posthocKW(data_test$clicks_per_dolphin, data_test$ID))
 hist(acoustic.dta$ID)
 hist(acoustic.dta$number)
+Adapted from Yannick OUTREMAN
+#    Agrocampus Ouest - 2020
+#######################################################################
+library(pscl)
+library(MASS)
+library(lmtest)
+library(multcomp)
+library(emmeans)
+library(dplyr)        # "%>%" function
+library(forcats)      # "fct_relevel" function
+library(stringr)      # "gsub" function
+library(rcompanion)   # "fullPTable" function
+library(multcompView) # "multcompLetters" function
+library(ggplot2)
+library(pgirmess)
+library(postHoc)
+#library(tidyquant)    # geom_ma() if rolling average needed
+################# DATASET IMPORTS #####################################
+folder <- './../'
+whistles.dta <-read.table(file=paste0(folder,
+'Whistles/Evaluation/whistles_durations.csv'),
+sep = ',', header=TRUE)
+whistles.dta <- whistles.dta[order(whistles.dta$audio_names),]
+bbp.dta <-read.table(file=paste0(folder,
+'BBPs/Results/16-06-22_14h00_number_of_BBP.csv'),
+sep = ',', header=TRUE)
+bbp.dta <- bbp.dta[order(bbp.dta$audio_names),]
+clicks.dta <-read.table(file=paste0(folder,
+'Clicks/Results/projection_updated_number_of_clicks_02052022.csv'), #number_of_clicks_02052022.csv
+sep = ',', header=TRUE)
+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)
+# 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)
diff --git a/Stats/BBP-click-whistles_3models.R b/Stats/BBP-click-whistles_3models.R
index 9d674b6e804cd97e7023279ff4f226c9d09e3dc8..0744f046b3664c566e051d6c79a80f7d30009270 100644
--- a/Stats/BBP-click-whistles_3models.R
+++ b/Stats/BBP-click-whistles_3models.R
@@ -76,10 +76,10 @@ plot(acoustic.dta) # nothing that we can see
 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)
-# 3.6%, 53.7% & 24.7% of our data are zeros. Will have to be dealt with.
+# 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 bytheir behavior, the emission of a pinger or a fishing net.
+# 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.
@@ -174,7 +174,6 @@ nb.whi <- glm.nb(total_whistles_duration ~
 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)
-shapiro.test(residuals(mod.whi)) # H0 : normality -> not rejected if p>0.05
 dwtest(mod.whi) # H0 -> independent if p>0.05 (autocorrelation if p<0.05)
 bptest(mod.whi) # H0 -> homoscedasticity if p<0.05
 # No normality but we do not need it
@@ -190,7 +189,6 @@ mod.bbp <- glm.nb(number_of_bbp ~ acoustic + fishing_net + behavior
 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
-# Normality not needed in GLM, hypotheses verified !
 mod.bbp$deviance/mod.bbp$df.residual 
 # slight underdispersion, not improved with ZINB so we keep this
 
@@ -203,7 +201,6 @@ car::Anova(mod.cli, type=3)
 shapiro.test(residuals(mod.cli)) # H0 : normality -> cannot be rejected if p > 0.05
 dwtest(mod.cli) # H0 -> independent if p>0.05 (autocorrelation if p<0.05)
 bptest(mod.cli) # H0 -> homoscedasticity if p<0.05
-# Normality not needed in GLM, hypotheses verified !
 mod.cli$deviance/mod.cli$df.residual
 # slight overdispersion. (ZINB does not clearly improve results so we keep this)