Skip to content
Snippets Groups Projects
Commit f86d3bbd authored by Loïc Lehnhoff's avatar Loïc Lehnhoff
Browse files

README update

parent ab21a120
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
README.md
+ 1
− 3
View file @ f86d3bbd
......@@ -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
Stats/.Rhistory
+ 3
− 474
View file @ f86d3bbd
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)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment