Skip to content
Snippets Groups Projects
Select Git revision
  • 4a92670e67b110d81d84b95095afa657815ccf5a
  • master default
2 results

__init__.py

Blame
  • .Rhistory 24.19 KiB
    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  ####
    # whistles
    table <- cld(emmeans(mod.whi, pairwise~fishing_net, adjust="tukey"), Letters = letters)
    myletters_df <- data.frame(fishing_net=table$fishing_net,
    letter = trimws(table$.group))
    barPlot(computeStats(acoustic.dta, fishing_net, whistling_time_per_dolphin/n_bins),
    myletters_df, fishing_net,
    old_names = c("SSF","F"), new_names = c("Absent", "Present"),
    xname="Presence/Asence of fishing net", height=1,
    ytitle="Mean whistling time per dolphin per min (in sec)")
    # BBP
    table <- cld(emmeans(mod.bbp, pairwise~fishing_net, adjust="tukey"), Letters = letters)
    myletters_df <- data.frame(fishing_net=table$fishing_net,
    letter = trimws(table$.group))
    barPlot(computeStats(acoustic.dta, fishing_net, BBPs_per_dolphin),
    myletters_df, fishing_net,
    old_names = c("SSF","F"), new_names = c("Absent", "Present"),
    xname="Presence/Asence of fishing net", height=.6,
    ytitle="Mean number of BBPs per dolphin per min")
    # Clicks
    table <- cld(emmeans(mod.cli, pairwise~fishing_net, adjust="tukey"), Letters = letters)
    myletters_df <- data.frame(fishing_net=table$fishing_net,
    letter = trimws(table$.group))
    barPlot(computeStats(acoustic.dta, fishing_net, clicks_per_dolphin),
    myletters_df, fishing_net,
    old_names = c("SSF","F"), new_names = c("Absent", "Present"),
    xname="Presence/Asence of fishing net", height=100,
    ytitle="Mean number of clicks per dolphin per min")
    #### Acoustic plots  ####
    # Whistles
    table <- cld(emmeans(mod.whi, pairwise~acoustic, adjust="tukey"), Letters = letters)
    myletters_df <- data.frame(acoustic=table$acoustic,letter = trimws(table$.group))
    barPlot(computeStats(acoustic.dta, acoustic, whistling_time_per_dolphin/n_bins),
    myletters_df, acoustic, height=1.3, ytitle="Mean whistling time per dolphin per min (in sec)",
    old_names = c("AV","AV+D","D","D+AP","AP"),
    new_names = c("BEF","BEF+DUR","DUR", "DUR+AFT", "AFT"),
    xname="Activation sequence")
    # BBPs
    table <- cld(emmeans(mod.bbp, pairwise~acoustic, adjust="tukey"), Letters = letters)
    myletters_df <- data.frame(acoustic=table$acoustic,letter = trimws(table$.group))
    barPlot(computeStats(acoustic.dta, acoustic, BBPs_per_dolphin),
    myletters_df, acoustic, height=1.2, ytitle="Mean number of BBPs per dolphin per min",
    old_names = c("AV","AV+D","D","D+AP","AP"),
    new_names = c("BEF","BEF+DUR","DUR", "DUR+AFT", "AFT"),
    xname="Activation sequence")
    # Clicks
    table <- cld(emmeans(mod.cli, pairwise~acoustic, adjust="tukey"), Letters = letters)
    myletters_df <- data.frame(acoustic=table$acoustic,letter = trimws(table$.group))
    barPlot(computeStats(acoustic.dta, acoustic, clicks_per_dolphin),
    myletters_df, acoustic, height=155, ytitle="Mean number of clicks per dolphin per min",
    old_names = c("AV","AV+D","D","D+AP","AP"),
    new_names = c("BEF","BEF+DUR","DUR", "DUR+AFT", "AFT"),
    xname="Activation sequence")
    #### Interaction fishing_net:acoustic plots  ####
    # Whistles
    letters_df <- computeLetters(emmeans(mod.whi, pairwise~fishing_net:acoustic, adjust="tukey"),
    "fishing_net")
    letters_df$acoustic <- computeLetters(emmeans(mod.whi, 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, whistling_time_per_dolphin/n_bins, two=acoustic),
    NULL, acoustic, fill=fishing_net,
    old_names = c("AV","AV+D","D","D+AP","AP"), ytitle="Mean whistling time per dolphin per min (in sec)",
    new_names = c("BEF","BEF+DUR","DUR", "DUR+AFT", "AFT"),
    xname="Activation sequence", height=c(.95,.95,.95,1,.95,1,.95,1,1,1),
    colours=c("#E69F00","#999999"), size=5,
    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")
    View(acoustic.dta)
    sum(acoustic.dta$date=="09/07/2021")
    56/361
    ########################################################################
    #    STATISTICS
    #    Author : Loic LEHNHOFF
    #    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
    n_bins = 187.5 # number of bins per sec for spectrograms (whistles)
    ################# 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)
    # 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.
    View(acoustic.dta)
    sum(acoustic.dta$date=="09/07/2021")
    sum(acoustic.dta$acoustic=="T")