Vocal coordination analysis

Waterrail vocal coordination

Source code and data found at

Purpose

• Test if waterrail pairs coordinate their calls

Analysis flowchart

flowchart LR
  A[Read data] --> B(Format data) 
  B --> C(Simulate random IDs)
  C --> D(Run coordination)
  D --> E(Check coordination scores) 

style A fill:#44015466
style B fill:#3E4A894D
style C fill:#26828E4D
style D fill:#6DCD594D

Load packages

# knitr is require for creating html/pdf/word reports
# formatR is used for soft-wrapping code

# install/ load packages
sketchy::load_packages(
  packages = c("warbleR", "pbapply", "readxl", "ggplot2", "gghalves", "formatR", "cowplot")
)

# custom function based on gghalves::geom_half_point
geom_half_point <- function(
  mapping = NULL,
  data = NULL,
  side = "l",
  range_scale = 0.4,
  alpha = 1,
  pch = 16,
  size = 1.5,
  transformation = NULL,   # <- absorb but ignore
  ...
) {

  ggplot2::geom_point(
    mapping = mapping,
    data = data,
    alpha = alpha,
    shape = pch,
    size = size,
    position = position_half_jitter(
      width = 0.1,
      height = 0,
      side = side,
      range_scale = range_scale
    ),
    ...
  )
}

position_half_jitter <- function(width = 0.1,
                                 height = 0,
                                 side = "l",
                                 range_scale = 0.4) {

  ggplot2::ggproto(
    NULL, ggplot2::Position,

    compute_panel = function(self, data, params, scales) {

      # apply regular jitter first
      data$x <- data$x + stats::runif(
        nrow(data),
        -width / 2,
        width / 2
      )

      # shift to one side
      shift <- if (side %in% c("l", "left")) {
        -range_scale / 2
      } else {
        range_scale / 2
      }

      data$x <- data$x + shift

      data
    }
  )
}

source("./scripts/test_coordination_impute.R")
source("~/Dropbox/R_package_testing/warbleR/R/internal_functions.R")

1 Read and format data

dat <-readxl::read_excel("./data/raw/Water_rail.xlsx")
dat <-readxl::read_excel("./data/raw/Water rail database v4.xlsx")
dat <- as.data.frame(dat)
dat$sing.event <- paste(dat$No., dat$`Before dawn (0), after dawn (1)`, sep = "-")
dat$start <- as.numeric(dat$`Call series begin time (s)`)
dat$end <- as.numeric(dat$`Call series end time (s)`)
dat$overlp <- ifelse(dat$`Overlapping series (no-0/yes-1)` == 0, NA, 1)

dat$indiv <- NA

# set those that overlap
dat$indiv[!is.na(dat$overlp)] <- rep(c("a", "b"), sum(!is.na(dat$overlp)) / 2)

dat$indiv[is.na(dat$overlp)] <- c("a", "b")[rbinom(sum(is.na(dat$overlp)), 1, 0.5) + 1]

2 Explore data

Here we print 5 of the singing events to visualize the coordination of calls. The plot_coordination() function from the warbleR package is used to create a plot of the vocal events, showing their timing and overlap.

agg_ind <- aggregate(indiv ~ sing.event, dat, function(x) length(unique(x)))

dat2 <- dat[dat$sing.event %in% agg_ind$sing.event[agg_ind$indiv == 2],]

a <- warbleR::plot_coordination(dat2, img = FALSE, ovlp = TRUE, only.coor = FALSE)

Warning: At least one singing event with less than 10 vocalizations

print(a[1:5])

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

3 Random assigning of IDs to test coordination

The routine below:

• Assigns different IDs to the vocal events that overlap (i.e. those that are coordinated)
• Assigns random IDs to the vocal events that do not overlap and assigns
• An unlabelled signal that overlaps nothing gets its own independent binomial draw. Filling in all missing IDs this way is repeated 100 times, giving a set of complete datasets that reflect the uncertainty about who produced the unlabelled signals. * Each imputed dataset is then treated as one observed dataset and randomized 1000 times while keeping its imputed labels. Its observed coordination statistic is compared only with its own null distribution.
• Finally, the observed values, expected values, coordination scores and p-values are averaged across imputations. The averaged p-value equals the proportion of all randomized datasets, pooled across imputations, that are at least as extreme as their own imputed observed data.

To quantify coordination we use a modified version of the coordination score (sensu Araya-Salas et al. 2017) calculated as:

\[ \text{bounded score} = \frac{\text{obs} - \text{mean.random}} {\text{obs} + \text{mean.random}} \]

Positive values indicate a tendency to overlap while negative values indicate a tendency to alternate. NA values will be returned when events cannot be randomized (e.g. too few signals).

3.1 Time lapse among vocal units

# probabilities controling the proportion of calls for the 2
# individuals
probs <- seq(0.1, 0.5, by = 0.05)

# set global options (this can also be set within the function
# call)
warbleR_options(iterations = 1000, pb = FALSE, ovlp.method = "time.closest")

# run over different probability values
out <- lapply(probs, function(x) {

    rep_out <- pbreplicate(n = 100, expr = test_coordination_impute(dat)$coor.score.bounded, cl = 20)

    return(rep_out)
})


saveRDS(out, "./data/processed/coordination_scores_by_id_proportion_time_lapse.RDS")

coor_scores_list <- readRDS("./data/processed/coordination_scores_by_id_proportion_time_lapse.RDS")

# probabilities controling the proportion of calls for the 2 individuals
probs <- seq(0.1, 0.5, by = 0.05)

coor_scores_list <- lapply(seq_along(coor_scores_list), function(x) {
    X <- coor_scores_list[[x]]
    out <- data.frame(prob = probs[x], coor_scores = c(X))
    out$proportion <- sum(out$coor_scores < 0) / nrow(out)
    
    return(out)
})

coor_scores_df <- do.call(rbind, coor_scores_list)


coor_scores_df$prob_label <- paste0(
    coor_scores_df$prob,
    "/",
    1 - coor_scores_df$prob,
    " (",
    round(coor_scores_df$proportion * 100, 0),
    "% < 0)"
)


labels_vec <- c(
  "0.1/0.9", "0.15/0.85", "0.2/0.8", "0.25/0.75",
  "0.3/0.7", "0.35/0.65", "0.4/0.6", "0.45/0.55",
  "0.5/0.5"
)

coor_scores_df$type <- "Time lapse"

gg_time <- ggplot(coor_scores_df, aes(x = factor(prob), y = coor_scores)) +
      geom_hline(yintercept = 0, col = "red2", lty = 2, lwd = 1) +
  ggdist::stat_halfeye(
    fill = "#A0DFB9FF",
    color = NA,
    alpha = 0.5,
    adjust = .5,
    width = .6,
    .width = 0,
    justification = -.2,
    point_colour = NA
  ) +

  geom_boxplot(
    fill = "#A0DFB9FF",
    color = "black",
    width = .15,
    outlier.shape = NA
  ) +

  geom_half_point(
    color = "#A0DFB9FF",
    side = "l",
    range_scale = .4,
    alpha = .7,
    transformation = ggplot2::position_jitter(height = 0)
  ) +

  scale_x_discrete(labels = labels_vec) +

  labs(
    x = "Proportion of calls (individual 1 / individual 2)",
    y = "Coordination score"
  ) +

  theme_classic() +
  theme(
    legend.position = "none"
  ) +
    facet_wrap(~ type, nrow = 1)

gg_time

Average proportion of coordination scores below 0 (i.e. alternating) per probability

aggregate(proportion ~ prob, coor_scores_df, mean)

prob	proportion
0.10	0.7513953
0.15	0.7495349
0.20	0.7527907
0.25	0.7511628
0.30	0.7520930
0.35	0.7509302
0.40	0.7523256
0.45	0.7511628
0.50	0.7516279

Average proportion of coordination scores below 0 (i.e. alternating) across all probabilities

mean(coor_scores_df$proportion)

[1] 0.751447

and its range:

range(coor_scores_df$proportion)

[1] 0.7495349 0.7527907

3.2 Overlap counts

# probabilities controling the proportion of calls for the 2
# individuals
probs <- seq(0.1, 0.5, by = 0.05)

# set global options (this can also be set within the function
# call)
warbleR_options(iterations = 100, pb = FALSE, ovlp.method = "count")

# run over different probability values
out <- lapply(probs, function(x) {

    rep_out <- pbreplicate(n = 100, expr = test_coordination_impute(X = dat[,
        c("sing.event", "indiv", "start", "end")])$coor.score.bounded, cl = 20)

    return(rep_out)
})


saveRDS(out, "./data/processed/coordination_scores_by_id_proportion_counts.RDS")

coor_scores_list <- readRDS("./data/processed/coordination_scores_by_id_proportion_counts.RDS")

# probabilities controling the proportion of calls for the 2 individuals
probs <- seq(0.1, 0.5, by = 0.05)

coor_scores_list <- lapply(seq_along(coor_scores_list), function(x) {
        X <- coor_scores_list[[x]]
    out <- data.frame(prob = probs[x], coor_scores = c(X))
    
    out$proportion <- sum(out$coor_scores > 0) / nrow(out)
        
    return(out)
    })

coor_scores_df <- do.call(rbind, coor_scores_list)


labels_vec <- c(
  "0.1/0.9", "0.15/0.85", "0.2/0.8", "0.25/0.75",
  "0.3/0.7", "0.35/0.65", "0.4/0.6", "0.45/0.55",
  "0.5/0.5"
)

coor_scores_df$type <- "Overlap counts"

gg_count <- ggplot(coor_scores_df, aes(x = factor(prob), y = coor_scores)) +
      geom_hline(yintercept = 0, col = "red2", lty = 2, lwd = 1) +
  ggdist::stat_halfeye(
    fill = "#A0DFB9FF",
    color = NA,
    alpha = 0.5,
    adjust = .5,
    width = .6,
    .width = 0,
    justification = -.2,
    point_colour = NA
  ) +

  geom_boxplot(
    fill = "#A0DFB9FF",
    color = "black",
    width = .15,
    outlier.shape = NA
  ) +

  geom_half_point(
    color = "#A0DFB9FF",
    side = "l",
    range_scale = .4,
    alpha = .7,
    transformation = ggplot2::position_jitter(height = 0)
  ) +

  scale_x_discrete(labels = labels_vec) +
  labs(
    x = "Proportion of calls (individual 1 / individual 2)",
    y = "Coordination score"
  ) +
  theme_classic() +
  theme(
    legend.position = "none"
  ) +
    facet_wrap(~ type, nrow = 1)

gg_count

Average proportion of coordination scores above 0 (i.e. synchronizing) per probability

aggregate(proportion ~ prob, coor_scores_df, mean)

prob	proportion
0.10	0.8604651
0.15	0.8604651
0.20	0.8604651
0.25	0.8604651
0.30	0.8604651
0.35	0.8604651
0.40	0.8604651
0.45	0.8604651
0.50	0.8604651

Average proportion of coordination scores above 0 (i.e. synchronizing) across all iterations and probabilities

mean(coor_scores_df$proportion)

[1] 0.8604651

and its range:

range(coor_scores_df$proportion)

[1] 0.8604651 0.8604651

# combine plots
cowplot::plot_grid(gg_time + ylim(c(-1, 1)), gg_count + ylim(c(-1, 1)), nrow = 2)

Takeaways

• Waterrail pairs seem to sing synchronously and overlap their calls, as indicated by the positive coordination scores

---

Session information

R version 4.6.1 (2026-06-24)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.5 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0  LAPACK version 3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       

time zone: America/Costa_Rica
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] cowplot_1.2.0      formatR_1.14       gghalves_0.1.4     ggplot2_4.0.3     
 [5] readxl_1.5.0       pbapply_1.7-4      warbleR_1.1.37     NatureSounds_1.0.5
 [9] knitr_1.52         seewave_2.2.4      tuneR_1.4.7       

loaded via a namespace (and not attached):
 [1] ggdist_3.3.3         generics_0.1.4       bitops_1.1-0        
 [4] stringi_1.8.9        digest_0.6.39        magrittr_2.0.5      
 [7] evaluate_1.0.5       grid_4.6.1           RColorBrewer_1.1-3  
[10] fastmap_1.2.0        sketchy_1.0.7        cellranger_1.1.0    
[13] jsonlite_2.0.0       brio_1.1.5           httr_1.4.9          
[16] scales_1.4.0         cli_3.6.6            rlang_1.3.0         
[19] crayon_1.5.3         fftw_1.0-9           withr_3.0.3         
[22] yaml_2.3.12          otel_0.2.0           tools_4.6.1         
[25] parallel_4.6.1       dplyr_1.2.1          curl_8.0.0          
[28] vctrs_0.7.3          R6_2.6.1             proxy_0.4-29        
[31] lifecycle_1.0.5      dtw_1.23-3           stringr_1.6.0       
[34] htmlwidgets_1.6.4    MASS_7.3-66          pkgconfig_2.0.3     
[37] pillar_1.11.1        gtable_0.3.6         glue_1.8.1          
[40] Rcpp_1.1.2           tidyselect_1.2.1     tibble_3.3.1        
[43] xfun_0.61            rstudioapi_0.19.0    farver_2.1.2        
[46] rjson_0.2.23         htmltools_0.5.9      labeling_0.4.3      
[49] rmarkdown_2.32       testthat_3.3.2       signal_1.8-1        
[52] compiler_4.6.1       S7_0.2.2             distributional_0.8.1
[55] RCurl_1.98-1.19