Vocal coordination analysis

Waterrail vocal coordination

Author
Published

February 23, 2026

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

Code
# 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"))

1 Read and format data

Code
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

Check if all events have 2 individuals

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

agg_ind
sing.event indiv
1-0 2
1-1 2
10-0 2
10-1 2
11-0 1
11-1 2
12-0 2
12-1 2
13-0 2
13-1 2
14-0 2
14-1 2
15-0 2
15-1 2
16-0 2
16-1 2
17-0 2
17-1 2
18-0 2
18-1 2
19-0 2
19-1 2
2-0 2
2-1 2
20-0 2
20-1 2
21-0 2
21-1 2
22-0 2
22-1 2
23-0 2
23-1 2
24-0 2
24-1 2
25-0 2
25-1 2
3-0 2
3-1 1
4-0 2
4-1 1
5-0 2
5-1 2
6-0 2
6-1 2
7-0 2
7-1 2
8-0 2
8-1 2
9-0 2
9-1 2
Code
dat <- dat[dat$sing.event %in% agg_ind$sing.event[agg_ind$indiv ==
    2], ]

a <- warbleR::plot_coordination(dat, img = FALSE, ovlp = TRUE, only.coor = FALSE)
Warning: At least one singing event with less than 10 vocalizations
Code
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 Run coordination test

coordination score (sensu Araya-Salas et al. 2017), calculated as:

  • (obs.overlap−mean.random.ovlp)/mean.random.ovlp

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

Code
# 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 = "count")

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

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

    return(rep_out)
})


saveRDS(out, "./data/processed/coordination_scores_by_id_proportion.RDS")
Code
coor_scores_list <- readRDS("./data/processed/coordination_scores_by_id_proportion.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


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