# Multicopter-assisted acoustic surveys compared with ground-based transects
# Analysis for the data file in this repository. Running the whole script
# reproduces Table 1 and the three GLMM results as reported in the manuscript.
# R >= 4.4

library(readr)
library(dplyr)
library(glmmTMB)
library(emmeans)
library(performance)
library(car)


#>>>>>>>>> settings <<<<<<<<<<<

# the multicopter flew at 120 m and carried the recorder 16 m below it, so
# everything aloft is reported at 104 m. the data file already says 104m.
group_levels  <- c("Myotis", "Nyctaloid", "Pipistrelloid")
habitat_levels <- c("forest", "meadow", "waterbody")

# Myotis is almost absent at 104 m (9 files in total), the model does not
# converge sensibly on that, so the aerial habitat model runs on the other two.
groups_aloft <- c("Nyctaloid", "Pipistrelloid")


#>>>>>>>>>>>> load the data <<<<<<<<<<<<<<<<

# one row per transect survey and species, bat_activity = number of one-second
# files holding at least one identified call. zeros are kept, they are real.
# the river Neckar site and the two species outside the three echolocation
# groups are already dropped in this file, nothing extra filtered here.
bats <- read_tsv("bat_activity_transects_EN.txt", show_col_types = FALSE) %>%
  mutate(
    echolocation_group = factor(echolocation_group, levels = group_levels),
    habitat            = factor(habitat, levels = habitat_levels),
    recording_height   = factor(recording_height, levels = c("ground", "104m")),
    multicopter        = factor(multicopter),
    # one transect survey: a habitat walked once in one night at one site with
    # one multicopter. ground and 104 m share it, that is the pairing unit.
    survey_id  = paste(survey_night, study_site, multicopter, habitat, transect_replicate),
    # the night block a survey sits in, used for the nested random effect
    site_night = paste(survey_night, study_site, multicopter)
  )


#>>>>>> Table 1: activity per group and height <<<<<<<<<

table_1 <- bats %>%
  group_by(echolocation_group, recording_height) %>%
  summarise(one_second_files = sum(bat_activity), .groups = "drop")
print(as.data.frame(table_1))
cat("total:", sum(table_1$one_second_files), "\n")


#>>>>>>>>>> aggregation for the models <<<<<<<<<<<<<<

# activity summed per group over the given keys
aggregate_by <- function(df, keys) {
  df %>% group_by(across(all_of(keys))) %>%
    summarise(activity_sum = sum(bat_activity), .groups = "drop")
}

# Q1 works per survey, all 87 of them
q1_data <- aggregate_by(bats,
  c("survey_id", "echolocation_group", "recording_height", "multicopter"))

# the two habitat models work per habitat and night block, which pools the
# repeat visits of 5 of the 79 units. dropping those 5 leaves every contrast
# on the same side of its threshold.
ground_data <- bats %>% filter(recording_height == "ground") %>%
  aggregate_by(c("site_night", "habitat", "echolocation_group", "multicopter"))

aloft_data <- bats %>%
  filter(recording_height == "104m", echolocation_group %in% groups_aloft) %>%
  mutate(echolocation_group = droplevels(echolocation_group)) %>%
  aggregate_by(c("site_night", "habitat", "echolocation_group", "multicopter"))


#>>>>>>> why negative binomial and not Poisson <<<<<<<<<<<<

# Pearson dispersion of the same model fitted as Poisson
poisson_dispersion <- function(formula, data) {
  m <- glmmTMB(formula, family = poisson, data = data)
  sum(residuals(m, type = "pearson")^2) / df.residual(m)
}


#>>>>>>>>>>>>> the three models <<<<<<<<<<

# Q1, does activity differ between ground and 104 m per group. the random
# intercept is the survey, which is what pairs the two heights.
model_height <- glmmTMB(
  activity_sum ~ echolocation_group * recording_height * multicopter + (1 | survey_id),
  family = nbinom2, data = q1_data)

# Q2 on the ground, does habitat matter. habitat sits inside the night block,
# so the random part is nested. this also covers the three habitat transects of
# one night not being independent of each other.
model_habitat_ground <- glmmTMB(
  activity_sum ~ echolocation_group * habitat * multicopter + (1 | site_night/habitat),
  family = nbinom2, data = ground_data)

# Q2 at 104 m, same structure without Myotis
model_habitat_aloft <- glmmTMB(
  activity_sum ~ echolocation_group * habitat * multicopter + (1 | site_night/habitat),
  family = nbinom2, data = aloft_data)


#>>>>>> post hoc contrasts <<<<<<<<<<<<<<<

# Tukey within each echolocation group, so the correction runs over the three
# contrasts that are actually of interest and not over all fifteen.
height_contrasts <- pairs(
  emmeans(model_height, ~ recording_height | echolocation_group, type = "response"),
  adjust = "tukey")

ground_contrasts <- pairs(
  emmeans(model_habitat_ground, ~ habitat | echolocation_group, type = "response"),
  adjust = "tukey")

aloft_contrasts <- pairs(
  emmeans(model_habitat_aloft, ~ habitat | echolocation_group, type = "response"),
  adjust = "tukey")

# Q3, the two multicopters against each other at each height. taken from the
# habitat models, averaged over habitat and group.
multicopter_ground <- pairs(emmeans(model_habitat_ground, ~ multicopter))
multicopter_aloft  <- pairs(emmeans(model_habitat_aloft,  ~ multicopter))


#>>>>>>>>>>> run <<<<<<<<<<<<<<<<<

# dispersion of the three models under Poisson, for comparison
cat("\n>>> Poisson dispersion of the same models <<<\n")
cat("height model      :", round(poisson_dispersion(
  activity_sum ~ echolocation_group * recording_height * multicopter + (1 | survey_id), q1_data), 1), "\n")
cat("habitat, ground   :", round(poisson_dispersion(
  activity_sum ~ echolocation_group * habitat * multicopter + (1 | site_night/habitat), ground_data), 1), "\n")
cat("habitat, 104 m    :", round(poisson_dispersion(
  activity_sum ~ echolocation_group * habitat * multicopter + (1 | site_night/habitat), aloft_data), 1), "\n")

# Q1 contrasts, model fit and the type III omnibus tests
cat("\n>>> Q1 ground versus 104 m <<<\n")
print(as.data.frame(height_contrasts))
cat("conditional R2:", round(r2_nakagawa(model_height)$R2_conditional, 3), "\n")
print(Anova(model_height, type = "III"))

# Q2 contrasts and model fit, ground level
cat("\n>>> Q2 habitat on the ground <<<\n")
print(as.data.frame(ground_contrasts))
cat("conditional R2:", round(r2_nakagawa(model_habitat_ground)$R2_conditional, 3), "\n")

# Q2 contrasts and model fit, 104 m
cat("\n>>> Q2 habitat at 104 m <<<\n")
print(as.data.frame(aloft_contrasts))
cat("conditional R2:", round(r2_nakagawa(model_habitat_aloft)$R2_conditional, 3), "\n")

# Q3, the two multicopter types at each height
cat("\n>>> Q3 multicopter type <<<\n")
cat("ground:\n"); print(as.data.frame(multicopter_ground))
cat("104 m:\n");  print(as.data.frame(multicopter_aloft))