Flight coordination and calling rate
Group flight coordination in Thyroptera
Marcelo Araya-Salas, PhD & Gloriana Chaverri, PhD
06-09-2021
Load packages
## add 'developer' to packages to be installed from github
x <- c("devtools", "maRce10/warbleR", "maRce10/Rraven", "readxl", "viridis", "ggplot2")
aa <- lapply(x, function(y) {
# get pakage name
pkg <- strsplit(y, "/")[[1]]
pkg <- pkg[length(pkg)]
# check if installed, if not then install
if (!pkg %in% installed.packages()[,"Package"]) {
if (grepl("/", y)) devtools::install_github(y, force = TRUE) else
install.packages(y)
}
# load package
a <- try(require(pkg, character.only = T), silent = T)
if (!a) remove.packages(pkg)
})Functions and global parameters
warbleR_options(wav.path = "~/Dropbox/Recordings/flight_coordination_Thyroptera/converted_sound_files_90_kHz/", wl = 300, parallel = parallel::detectCores() - 4, bp = "frange", fast = F, threshold = 15, ovlp = 20)
opts_knit$set(root.dir = "..")
# set evaluation false
opts_chunk$set(fig.width = 12, fig.height = 12, warning = FALSE, message = FALSE)
# frames per second
fps <- 120clls <- readRDS("./data/processed/curated_extended_selection_table_inquiry_calls_2020_&_2021.RDS")
metadat <- read.csv("./data/processed/metadata_inquiry_calls_2020_&_2021.csv")
# list files
video_analysis_folders <- list.dirs("/home/m/Dropbox/Salidas de videos analizados/Salidas50FramesXBicho", full.names = FALSE)
# filter
video_analysis_folders <- video_analysis_folders[video_analysis_folders != "" & !grepl("volver", video_analysis_folders)]
# get name of videos
videos <- sapply(strsplit(video_analysis_folders, "_"), "[[", 1)
videos <- as.numeric(gsub("T", "", videos))
# get name.video variable to match data in metadat and clls
years <- as.numeric(sapply(strsplit(video_analysis_folders, "_"), "[[", 3))
year.video <- paste(years, videos, sep = "-")
# put in data frame
video_df <- data.frame(directories = video_analysis_folders, video = videos, year = years, year.video = year.video)
# add metadata
metadat$year.video <- paste(metadat$Video, metadat$year, sep = "-")
video_df$year.audio <- sapply(1:nrow(video_df), function(x){
metadat$year.audio[metadat$year.video == video_df$year.video[x]][1]
})
metadat$year.video <- paste(metadat$year, metadat$Video, sep = "-")
video_df$experiment <- sapply(1:nrow(video_df), function(x){
metadat$Experimento[metadat$year.video == video_df$year.video[x]][1]
})
video_df$flight.time <- sapply(1:nrow(video_df), function(x){
metadat$Tiempo.de.vuelo[metadat$year.video == video_df$year.video[x]][1]
})
# fix flight time
video_df$flight.time <- video_df$flight.time * 1440
# get number of calls
video_df$n.calls <- sapply(1:nrow(video_df), function(x){
sum(clls$year.audio == video_df$year.audio[x])
})
# get chek.results data to get original time coordinates
attr_clls <- attr(clls, "check.results")
# convert to a selection table refering to original sound files
clls <- as.data.frame(clls)
clls$sound.files <- attr_clls$orig.sound.files
clls$start <- attr_clls$orig.start
clls$end <- attr_clls$orig.end
clls$selec <- attr_clls$orig.selec
clls <- clls[clls$year.audio %in% video_df$year.audio, ]
# check_sels(clls, parallel = 1)exp_raven(X = clls, path = "./data/processed", file.name = "audio_files_for_50_frame_video_analysis", sound.file.path = .Options$warbleR$path)
# manually annotate clapsclap_calls <- imp_raven(files = "audio_files_for_50_frame_video_analysis.txt", path = "./data/processed", all.data = TRUE, warbler.format = TRUE)
# fix time position starting from clap and remove clap row
clap_calls <- do.call(rbind, lapply(unique(clap_calls$sound.files), function(x){
X <- clap_calls[clap_calls$sound.files == x, ]
X$start <- X$start - X$start[X$class == "clap"]
# remove end just to make sure is not used
X$end <- NULL
return(X)
}))# get call rate for the first 2 minutes
video_df$n.calls.2.min <- sapply(1:nrow(video_df), function(x){
sum(clap_calls$year.audio == video_df$year.audio[x] & clap_calls$start <= 120)
})Flight pairwise distance and call rate
fls <- list.files("/home/m/Dropbox/Salidas de videos analizados/Salidas50FramesXBicho", full.names = TRUE, recursive = TRUE)
fls <- grep("unpaired", fls, value = TRUE)
fls <- grep("volver", invert = TRUE, fls, value = TRUE)
# read data with frame position for clap
frame_pos <- as.data.frame(read_excel("~/Dropbox/Salidas de videos analizados/Salidas50FramesXBicho/framesXToma50Puntos.xlsx"))
# fix names
names(frame_pos)[-1] <- paste(sapply(strsplit(names(frame_pos)[-1], "_"), "[[", 3), gsub("^T", "", sapply(strsplit(names(frame_pos)[-1], "_"), "[[", 1)), sep = "-")
# get mean distance and time of frame
video.3d.coords_l <- lapply(fls, function(x){
X <- read.csv(x)
X$video <- gsub("-unpaired-points-xyz.csv", "", basename(x))
# get coordinates in long format
coor_df <- data.frame(x = stack(X[, grep("x", names(X))])$values, y = stack(X[, grep("y", names(X))])$values, z = stack(X[, grep("z", names(X))])$values, ind = rep(1:((ncol(X) -1) / 3), each = nrow(X)), frame = 1:nrow(X))
# add year.video label
X$year.video <- paste(sapply(strsplit(X$video[1], "_"), "[[", 3), gsub("^T", "", sapply(strsplit(X$video[1], "_"), "[[", 1)), sep = "-")
# get distances among bats
dists <- sapply(unique(coor_df$frame), function(w){
W <- coor_df[coor_df$frame == w, ]
dst <- dist(W[, !names(coor_df) %in% c("ind", "frame")])
return(mean(dst, na.rm = TRUE)[1])
})
# put into a data frame
dist_df <- data.frame(year.video = X$year.video[1], frame = unique(coor_df$frame), distance = dists, n.ind = length(unique(coor_df$ind)))
# remove NA rows
dist_df <- dist_df[complete.cases(dist_df), ]
# get column with position of frames
sub_frame_pos <- frame_pos[, names(frame_pos) %in% X$year.video[1]]
sub_frame_pos <- c(na.omit(sub_frame_pos))
# make frame position relative to clap
sub_frame_pos <- sub_frame_pos - sub_frame_pos[1]
# convert to seconds
frame_time <- sub_frame_pos / fps
# remove clap
frame_time <- frame_time[-1]
# order by frame
dist_df <- dist_df[order(dist_df$frame), ]
#######CHECK####
if (nrow(dist_df) > length(frame_time))
dist_df <- dist_df[1:length(frame_time), ]
if (nrow(dist_df) < length(frame_time))
frame_time <- frame_time[1:nrow(dist_df)]
# add to output data frame
dist_df$frame_time <- frame_time
return(dist_df)
})
## add call rate 10 s before frame
mean_dist_call_rate_l <- lapply(video.3d.coords_l, function(x){
sub_clap <- clap_calls[clap_calls$year.audio == metadat$year.audio[metadat$year.video == x$year.video[1]][1],]
x$call.count.10.s <- sapply(1:nrow(x), function(y)
sum(sub_clap$start > (x$frame_time[y] - 10) & sub_clap$start < x$frame_time[y])
)
x$norm.call.count.10.s <- x$call.count.10.s / max(x$call.count.10.s)
x$call.count.20.s <- sapply(1:nrow(x), function(y)
sum(sub_clap$start > (x$frame_time[y] - 20) & sub_clap$start < x$frame_time[y])
)
x$norm.call.count.20.s <- x$call.count.20.s / max(x$call.count.20.s)
x$norm.distance <- x$distance / max(x$distance)
return(x)
})
# put in a data frame
mean_dist_call_rate <- do.call(rbind, mean_dist_call_rate_l)# for ggplot
mean_dist_call_rate_gg_l <- lapply(mean_dist_call_rate_l, function(x)
data.frame(type = rep(c("distance", "call.count.10", "call.count.20"), each = nrow(x)), value = c(x$norm.distance, x$norm.call.count.10.s, x$norm.call.count.20.s), year.video = x$year.video[1], time = rep(x$frame_time, 3))
)
mean_dist_call_rate_gg <- do.call(rbind, mean_dist_call_rate_gg_l)
ggplot(mean_dist_call_rate_gg, aes(x = time, y = value, group = type, color = type)) +
geom_point() +
geom_line() +
scale_color_viridis_d(end = 0.8, alpha = 0.6) +
facet_wrap(~ year.video, scales = "free_x", ncol = 4) +
theme_classic() 
Statistical analysis
Bayesian mixed effect models:
- distance ~ call.count.10.s + (1 | year.video)
- distance ~ call.count.20.s + (1 | year.video)
No p-value. Statistical significance can be evaluated as credibility interval (l-95% CI and u-95%) not overlapping with 0.
brm_model_10s <- brm(iter = 2500,
distance ~ call.count.10.s + (1 | year.video),
data = mean_dist_call_rate,
family = gaussian(), silent = 2,
# control = list(adapt_delta = 0.9),
chains = 2
#,
# prior = c(
# prior(normal(0, 10), "b"),
# prior(normal(0, 50), "Intercept"),
# prior(student_t(3, 0, 20), "sd"),
# prior(student_t(3, 0, 20), "sigma")
# )
)
brm_model_20s <- brm(iter = 2500,
distance ~ call.count.20.s + (1 | year.video),
data = mean_dist_call_rate,
family = gaussian(), silent = 2,
chains = 2
)
saveRDS(list(brm_model_20s = brm_model_20s, brm_model_10s = brm_model_10s), "./data/processed/brms_models_coordination_vs_call_rate.RDS")10 s model:
attach(readRDS("./data/processed/brms_models_coordination_vs_call_rate.RDS"))
summary(brm_model_10s)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: distance ~ call.count.10.s + (1 | year.video)
## Data: mean_dist_call_rate (Number of observations: 890)
## Draws: 2 chains, each with iter = 5000; warmup = 2500; thin = 1;
## total post-warmup draws = 5000
##
## Group-Level Effects:
## ~year.video (Number of levels: 18)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.55 0.11 0.39 0.80 1.00 610 1160
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 1.63 0.14 1.36 1.91 1.00 537 857
## call.count.10.s 0.01 0.01 -0.01 0.04 1.00 1856 2685
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.53 0.01 0.51 0.56 1.00 2615 3064
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).20 s model:
summary(brm_model_20s)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: distance ~ call.count.20.s + (1 | year.video)
## Data: mean_dist_call_rate (Number of observations: 890)
## Draws: 2 chains, each with iter = 2500; warmup = 1250; thin = 1;
## total post-warmup draws = 2500
##
## Group-Level Effects:
## ~year.video (Number of levels: 18)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.56 0.11 0.39 0.80 1.01 196 496
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 1.66 0.14 1.38 1.93 1.00 231 424
## call.count.20.s 0.00 0.01 -0.01 0.01 1.00 1173 1508
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.53 0.01 0.51 0.56 1.00 1083 1189
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Session information
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
##
## locale:
## [1] LC_CTYPE=pt_BR.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=es_CR.UTF-8 LC_COLLATE=pt_BR.UTF-8
## [5] LC_MONETARY=es_CR.UTF-8 LC_MESSAGES=pt_BR.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
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggplot2_3.3.5 viridis_0.6.1 viridisLite_0.4.0 readxl_1.3.1
## [5] Rraven_1.0.13 warbleR_1.1.27 NatureSounds_1.0.4 knitr_1.33
## [9] seewave_2.1.8 tuneR_1.3.3 devtools_2.4.2 usethis_2.0.1
##
## loaded via a namespace (and not attached):
## [1] backports_1.2.1 plyr_1.8.6 igraph_1.2.6
## [4] splines_4.1.0 crosstalk_1.1.1 rstantools_2.1.1
## [7] inline_0.3.19 digest_0.6.27 htmltools_0.5.1.1
## [10] rsconnect_0.8.24 fansi_0.5.0 magrittr_2.0.1
## [13] checkmate_2.0.0 memoise_2.0.0 remotes_2.4.0
## [16] brms_2.16.1 RcppParallel_5.1.4 matrixStats_0.60.1
## [19] xts_0.12.1 prettyunits_1.1.1 colorspace_2.0-2
## [22] signal_0.7-7 xfun_0.24 dplyr_1.0.7
## [25] callr_3.7.0 crayon_1.4.1 RCurl_1.98-1.3
## [28] jsonlite_1.7.2 lme4_1.1-27.1 zoo_1.8-9
## [31] glue_1.4.2 gtable_0.3.0 V8_3.4.2
## [34] distributional_0.2.2 pkgbuild_1.2.0 rstan_2.21.2
## [37] abind_1.4-5 scales_1.1.1 mvtnorm_1.1-2
## [40] DBI_1.1.1 miniUI_0.1.1.1 Rcpp_1.0.7
## [43] dtw_1.22-3 xtable_1.8-4 diffobj_0.3.4
## [46] proxy_0.4-26 StanHeaders_2.21.0-7 stats4_4.1.0
## [49] DT_0.19 htmlwidgets_1.5.3 threejs_0.3.3
## [52] posterior_1.0.1 ellipsis_0.3.2 pkgconfig_2.0.3
## [55] loo_2.4.1 farver_2.1.0 sass_0.4.0
## [58] utf8_1.2.2 tidyselect_1.1.1 labeling_0.4.2
## [61] rlang_0.4.11 reshape2_1.4.4 later_1.3.0
## [64] munsell_0.5.0 cellranger_1.1.0 tools_4.1.0
## [67] cachem_1.0.5 cli_3.0.1 generics_0.1.0
## [70] ggridges_0.5.3 evaluate_0.14 stringr_1.4.0
## [73] fastmap_1.1.0 yaml_2.2.1 processx_3.5.2
## [76] fs_1.5.0 purrr_0.3.4 pbapply_1.4-3
## [79] nlme_3.1-152 mime_0.11 projpred_2.0.2
## [82] compiler_4.1.0 bayesplot_1.8.1 shinythemes_1.2.0
## [85] rstudioapi_0.13 curl_4.3.2 gamm4_0.2-6
## [88] testthat_3.0.4 tibble_3.1.3 bslib_0.2.5.1
## [91] stringi_1.7.3 highr_0.9 ps_1.6.0
## [94] desc_1.3.0 Brobdingnag_1.2-6 lattice_0.20-44
## [97] Matrix_1.3-4 nloptr_1.2.2.2 markdown_1.1
## [100] shinyjs_2.0.0 fftw_1.0-6 tensorA_0.36.2
## [103] vctrs_0.3.8 pillar_1.6.1 lifecycle_1.0.0
## [106] jquerylib_0.1.4 bridgesampling_1.1-2 bitops_1.0-7
## [109] httpuv_1.6.2 R6_2.5.0 promises_1.2.0.1
## [112] gridExtra_2.3 codetools_0.2-18 sessioninfo_1.1.1
## [115] boot_1.3-28 colourpicker_1.1.0 MASS_7.3-54
## [118] gtools_3.9.2 assertthat_0.2.1 pkgload_1.2.1
## [121] rprojroot_2.0.2 rjson_0.2.20 withr_2.4.2
## [124] shinystan_2.5.0 mgcv_1.8-36 parallel_4.1.0
## [127] grid_4.1.0 coda_0.19-4 minqa_1.2.4
## [130] rmarkdown_2.10 shiny_1.6.0 base64enc_0.1-3
## [133] dygraphs_1.1.1.6