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
Vocal coordination analysis
Waterrail vocal coordination
Purpose
- Test if waterrail pairs coordinate their calls
Analysis flowchart
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