# Recorder distances from multicopters to avoid bat disturbance
# Analysis for the two data files in this repository. Running the whole script
# reproduces the LMER interaction results and the paired Wilcoxon tests as
# reported in the manuscript. R >= 4.2.1

library(readr)
library(dplyr)
library(tidyr)
library(lmerTest)
library(performance)


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

# recorder height on the 12 m rod -> recording distance in m
height_to_distance <- c("12m" = 10.0, "9.8m" = 12.2, "7.2m" = 14.8,
                        "4.6m" = 17.4, "2m" = 20.0)

species_levels <- c("Myotini", "Nyctaloid", "Pipistrelloid")

# Five distances are tested per drone and group, so the Bonferroni-corrected
# level is 0.05/5 = 0.01. The p-values themselves stay unadjusted, only the
# threshold we judge them against changes.
alpha_bonferroni <- 0.05 / 5


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

# ConVecDro only holds the complete sequences (the 9 incomplete ones were
# already dropped, so 31 are left), Mavic keeps all of them because there would
# otherwise be too few data. That is already baked into the files, so nothing
# extra gets filtered here.
load_drone <- function(file) {
  df <- read_tsv(file, show_col_types = FALSE) %>%
    mutate(
      Distance     = as.numeric(height_to_distance[Height]),
      Distance_chr = sprintf("%.1f", Distance),
      Phase        = trimws(Phase),
      ID2          = paste(Survey_Night, Site, sep = "_"),
      bat_activity = as.numeric(bat_activity)
    )
  if (any(is.na(df$Distance))) {
    stop("unmapped Height value(s): ",
         paste(unique(df$Height[is.na(df$Distance)]), collapse = ", "))
  }
  df
}


#>>>>>> LMER: phase x distance interaction <<<<<<<<<<

# Fit bat_activity ~ Phase * Distance + (1 | ID2) for a set of distances and
# take the interaction p-value at each non-reference distance (the reference is
# the first distance level). R2 is the conditional Nakagawa R2.
fit_interaction <- function(df_species, dist_levels_chr, left_label, right_labels_chr) {
  d <- df_species %>%
    filter(Distance_chr %in% dist_levels_chr) %>%
    mutate(
      Distance = factor(Distance_chr, levels = dist_levels_chr),
      Phase    = factor(Phase, levels = c("Control", "Multicopter"))
    )

  out <- data.frame(
    Comparison     = paste0(left_label, " vs ", right_labels_chr),
    p_value        = NA_real_,
    R2_conditional = NA_real_,
    stringsAsFactors = FALSE
  )
  if (nrow(d) == 0) return(out)

  tryCatch({
    m   <- lmer(bat_activity ~ Phase * Distance + (1 | ID2), data = d)
    cf  <- coef(summary(m))
    int <- cf[grep("^PhaseMulticopter:Distance", rownames(cf)), , drop = FALSE]
    if (nrow(int) > 0) {
      dist <- sprintf("%.1f", as.numeric(sub("^PhaseMulticopter:Distance", "", rownames(int))))
      pmap <- setNames(as.numeric(int[, "Pr(>|t|)"]), dist)
      out$p_value <- unname(pmap[right_labels_chr])
    }
    out$R2_conditional <- as.numeric(r2_nakagawa(m)$R2_conditional)
    out
  }, error = function(e) out)
}

# all interaction comparisons for one drone (over the three groups)
analyze_interactions <- function(df, drone_name) {
  bind_rows(lapply(species_levels, function(sp) {
    df_sp <- df %>% filter(EcholocationGroup == sp)
    # reference 12.2 m against the upper four devices
    outA <- fit_interaction(df_sp, c("12.2", "14.8", "17.4", "20.0"),
                            "12.2", c("14.8", "17.4", "20.0"))
    # orientation pair: 10.0 m downward vs 12.2 m upward
    outB <- fit_interaction(df_sp, c("10.0", "12.2"), "10.0", c("12.2"))
    res <- rbind(outA, outB)
    res$Drone   <- drone_name
    res$Species <- sp
    res[, c("Drone", "Species", "Comparison", "p_value", "R2_conditional")]
  }))
}


#>>>>>>> paired Wilcoxon: control vs multicopter per distance <<<<<<<

# Phases are compared per distance with a paired Wilcoxon signed-rank test. The
# non-parametric test is used because the paired differences are clearly not
# normal (see shapiro_diff_p in the output), so a paired t-test would not be
# appropriate. Only complete Control/Multicopter pairs (same ID2) are kept,
# effect size is r = z / sqrt(n), and significant_bonf flags the ones below the
# Bonferroni level set above.
paired_wilcoxon <- function(df, drone_name, min_pairs = 3) {
  groups <- setdiff(unique(df$EcholocationGroup), "Noise")

  bind_rows(lapply(groups, function(group) {
    wide <- df %>%
      filter(EcholocationGroup == group) %>%
      group_by(ID2, Distance, Phase) %>%
      summarise(bat_activity = sum(bat_activity, na.rm = TRUE), .groups = "drop") %>%
      pivot_wider(names_from = Phase, values_from = bat_activity)

    if (!("Control"     %in% names(wide))) wide$Control     <- NA_real_
    if (!("Multicopter" %in% names(wide))) wide$Multicopter <- NA_real_
    wide <- wide %>% filter(!is.na(Control) & !is.na(Multicopter))
    if (nrow(wide) == 0) return(NULL)

    bind_rows(lapply(sort(unique(wide$Distance)), function(dist) {
      dd <- wide %>% filter(Distance == dist)
      n  <- nrow(dd)
      if (n < min_pairs) {
        return(tibble(EcholocationGroup = group, Distance = dist, n_pairs = n,
                      p_value = NA_real_, significant_bonf = NA,
                      effect_r = NA_real_, shapiro_diff_p = NA_real_,
                      Drone = drone_name))
      }
      te   <- suppressWarnings(wilcox.test(dd$Control, dd$Multicopter, paired = TRUE, exact = FALSE))
      z    <- abs(qnorm(te$p.value / 2))
      r    <- z / sqrt(n)
      diff <- dd$Control - dd$Multicopter
      sh_p <- if (n >= 3 && n <= 5000) shapiro.test(diff)$p.value else NA_real_
      tibble(EcholocationGroup = group, Distance = dist, n_pairs = n,
             p_value = as.numeric(te$p.value),
             significant_bonf = as.numeric(te$p.value) < alpha_bonferroni,
             effect_r = as.numeric(r), shapiro_diff_p = sh_p,
             Drone = drone_name)
    }))
  }))
}


#>>>>> assumption checks for one LMER <<<<<<<<<

# Shapiro-Wilk on the residuals plus a Breusch-Pagan test for homoscedasticity.
# check_heteroscedasticity() runs the BP test on the fitted mixed model.
check_assumptions <- function(model) {
  list(
    shapiro_residuals = shapiro.test(residuals(model)),
    breusch_pagan     = check_heteroscedasticity(model)
  )
}


#>>>>>>>> sensitivity of the interaction <<<<<<<<<<<<

# The raw-count residuals break normality and homoscedasticity, so this refits
# Model 1 on the raw, sqrt- and log1p-transformed response and reports both the
# assumption checks and the interaction p-values. log1p removes the
# heteroscedasticity in every group while the ns / * / ** pattern stays the
# same, which is the robustness check mentioned in the manuscript.
sensitivity_interaction <- function(df, drone_name, group) {
  d <- df %>%
    filter(EcholocationGroup == group, Distance %in% c(12.2, 14.8, 17.4, 20.0)) %>%
    mutate(Distance = factor(Distance_chr, levels = c("12.2", "14.8", "17.4", "20.0")),
           Phase    = factor(Phase, levels = c("Control", "Multicopter")))
  bind_rows(lapply(c("raw", "sqrt", "log1p"), function(tr) {
    d$Y  <- switch(tr, raw = d$bat_activity, sqrt = sqrt(d$bat_activity), log1p = log1p(d$bat_activity))
    m    <- lmer(Y ~ Phase * Distance + (1 | ID2), data = d)
    cf   <- coef(summary(m))
    rows <- grep("^PhaseMulticopter:Distance", rownames(cf))
    ip   <- cf[rows, "Pr(>|t|)"]
    names(ip) <- sub("^PhaseMulticopter:Distance", "", rownames(cf)[rows])
    tibble(
      Drone = drone_name, Group = group, transform = tr,
      shapiro_p = shapiro.test(residuals(m))$p.value,
      breusch_pagan_p = tryCatch(as.numeric(check_heteroscedasticity(m)), error = function(e) NA_real_),
      p_12.2_vs_14.8 = ip["14.8"], p_12.2_vs_17.4 = ip["17.4"], p_12.2_vs_20.0 = ip["20.0"]
    )
  }))
}


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

convecdro <- load_drone("bat_activity_ConVecDro_EN.txt")
mavic     <- load_drone("bat_activity_Mavic_EN.txt")

# LMER interaction (matches the "Interaktion" sheet)
interaction_results <- bind_rows(
  analyze_interactions(convecdro, "ConVecDro"),
  analyze_interactions(mavic,     "DJI Mavic 2 Pro")
)
print(as.data.frame(interaction_results))

# paired Wilcoxon tests
wilcoxon_results <- bind_rows(
  paired_wilcoxon(convecdro, "ConVecDro"),
  paired_wilcoxon(mavic,     "Mavic")
) %>% arrange(Drone, EcholocationGroup, Distance)
print(as.data.frame(wilcoxon_results))

# assumption checks and the sensitivity analysis are run per model when needed:
# sensitivity_interaction(convecdro, "ConVecDro", "Pipistrelloid")