Vocal coordination analysis

Waterrail vocal coordination

Source code and data found at

Purpose

• Test if waterrail pairs coordinate their calls

Analysis flowchart

flowchart
  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"))

source("~/Dropbox/R_package_testing/warbleR/R/test_coordination.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

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]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

[[10]]

[[11]]

[[12]]

[[13]]

[[14]]

[[15]]

[[16]]

[[17]]

[[18]]

[[19]]

[[20]]

[[21]]

[[22]]

[[23]]

[[24]]

[[25]]

[[26]]

[[27]]

[[28]]

[[29]]

[[30]]

[[31]]

[[32]]

[[33]]

[[34]]

[[35]]

[[36]]

[[37]]

[[38]]

[[39]]

[[40]]

[[41]]

[[42]]

[[43]]

[[44]]

[[45]]

[[46]]

[[47]]

3 Random assigning of IDs to test coordination while dealing with unknown ID

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

This allows to test how the coordination scores change as we vary the proportion of calls assigned to each individual (a and b) among the non-overlapping events. By running multiple iterations of this random assignment and calculating the coordination scores, we can assess whether the observed coordination is significantly different from what would be expected by chance by simulating all possible scenarios between an equal vocal activity between the two individuals (i.e. 50% of calls for each) and a scenario where one individual is more vocally active than the other (e.g. 10% of calls for one individual and 90% for the other).

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) {

    dat_list <- lapply(unique(dat$sing.event), function(y) {
        Y <- dat[dat$sing.event == y, ]
        Y <- Y[order(Y$start), ]
        Y$indiv <- NA

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

        # the rest assign them randomly
        Y$indiv[is.na(Y$overlp)] <- c("a", "b")[rbinom(sum(is.na(Y$overlp)),
            size = 1, prob = x) + 1]

        if (length(unique(Y$indiv)) == 1)
            Y <- NULL
        return(Y)
    })

    dat_sim <- do.call(rbind, dat_list)

    rep_out <- pbreplicate(n = 100, expr = test_coordination(dat_sim)$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))
        
    return(out)
    })

coor_scores_df <- do.call(rbind, coor_scores_list)

coor_scores_df$prob_label <- paste("Prop individual a:", coor_scores_df$prob)



ggplot(coor_scores_df,
aes(x = prob_label, y = coor_scores)) +
# add half-violin from {ggdist} package
ggdist::stat_halfeye(
# fill = fill_color,
alpha = 0.5,
# custom bandwidth
adjust = .5,
# adjust height
width = .6,
.width = 0,
# move geom to the cright
justification = -.2,
point_colour = NA
) +
geom_boxplot(# fill = fill_color,
width = .15,
# remove outliers
outlier.shape = NA) +
# add justified jitter from the {gghalves} package
gghalves::geom_half_point(
# color = fill_color,
# draw jitter on the left
side = "l",
# control range of jitter
range_scale = .4,
# add some transparency
alpha = .5,
transformation = ggplot2::position_jitter(height = 0)
) +
scale_color_viridis_d(option = "G", end = 0.8) +
scale_fill_viridis_d(option = "G",
end = 0.8,
alpha = 0.6) +
# ylim(c(NA, 250)) +
theme(legend.position = "none") + facet_wrap(~prob_label, scales = "free_x") +
 labs(x = "Proportion of calls for individual b", y = "Coordination score") + 
    geom_hline(yintercept = 0, col = "red2") +  theme(axis.text.x = element_blank(),  # Remove x-axis labels
        axis.ticks.x = element_blank()) + theme_classic()

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) {

    dat_list <- lapply(unique(dat$sing.event), function(y) {
        Y <- dat[dat$sing.event == y, ]
        Y <- Y[order(Y$start), ]
        Y$indiv <- NA

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

        # the rest assign them randomly
        Y$indiv[is.na(Y$overlp)] <- c("a", "b")[rbinom(sum(is.na(Y$overlp)),
            size = 1, prob = x) + 1]

        if (length(unique(Y$indiv)) == 1)
            Y <- NULL
        return(Y)
    })

    dat_sim <- do.call(rbind, dat_list)

    rep_out <- pbreplicate(n = 100, expr = test_coordination(X = dat_sim[,
        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))
        
    return(out)
    })

coor_scores_df <- do.call(rbind, coor_scores_list)

coor_scores_df$prob_label <- paste("Prop individual a:", coor_scores_df$prob)



ggplot(coor_scores_df,
aes(x = prob_label, y = coor_scores)) +
# add half-violin from {ggdist} package
ggdist::stat_halfeye(
# fill = fill_color,
alpha = 0.5,
# custom bandwidth
adjust = .5,
# adjust height
width = .6,
.width = 0,
# move geom to the cright
justification = -.2,
point_colour = NA
) +
geom_boxplot(# fill = fill_color,
width = .15,
# remove outliers
outlier.shape = NA) +
# add justified jitter from the {gghalves} package
gghalves::geom_half_point(
# color = fill_color,
# draw jitter on the left
side = "l",
# control range of jitter
range_scale = .4,
# add some transparency
alpha = .5,
transformation = ggplot2::position_jitter(height = 0)
) +
scale_color_viridis_d(option = "G", end = 0.8) +
scale_fill_viridis_d(option = "G",
end = 0.8,
alpha = 0.6) +
# ylim(c(NA, 250)) +
theme(legend.position = "none") + facet_wrap(~prob_label, scales = "free_x") +
 labs(x = "Proportion of calls for individual b", y = "Coordination score") + 
    geom_hline(yintercept = 0, col = "red2") +  theme(axis.text.x = element_blank(),  # Remove x-axis labels
        axis.ticks.x = element_blank()) + theme_classic()

Takeaways

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

---

Session information

R version 4.5.2 (2025-10-31)
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] gghalves_0.1.4     ggplot2_4.0.2      readxl_1.4.5       pbapply_1.7-4     
[5] warbleR_1.1.37     NatureSounds_1.0.5 knitr_1.51         seewave_2.2.4     
[9] tuneR_1.4.7       

loaded via a namespace (and not attached):
 [1] ggdist_3.3.3         generics_0.1.4       bitops_1.0-9        
 [4] stringi_1.8.7        digest_0.6.39        magrittr_2.0.4      
 [7] evaluate_1.0.5       grid_4.5.2           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           formatR_1.14        
[16] httr_1.4.7           scales_1.4.0         cli_3.6.5           
[19] rlang_1.1.7          crayon_1.5.3         fftw_1.0-9          
[22] withr_3.0.2          remotes_2.5.0        yaml_2.3.12         
[25] packrat_0.9.3        tools_4.5.2          parallel_4.5.2      
[28] dplyr_1.1.4          curl_7.0.0           vctrs_0.7.1         
[31] R6_2.6.1             proxy_0.4-29         lifecycle_1.0.5     
[34] dtw_1.23-1           stringr_1.6.0        htmlwidgets_1.6.4   
[37] MASS_7.3-65          pkgconfig_2.0.3      xaringanExtra_0.8.0 
[40] pillar_1.11.1        gtable_0.3.6         glue_1.8.0          
[43] Rcpp_1.1.1           tidyselect_1.2.1     tibble_3.3.0        
[46] xfun_0.56            rstudioapi_0.17.1    farver_2.1.2        
[49] rjson_0.2.23         htmltools_0.5.9      labeling_0.4.3      
[52] rmarkdown_2.30       testthat_3.2.3       signal_1.8-1        
[55] compiler_4.5.2       S7_0.2.1             distributional_0.5.0
[58] RCurl_1.98-1.17