Degradation statistics

Sex-dependent signal propagation in L. palavanensis

Author

Marcelo Araya-Salas

Published

July 11, 2024

Source code and data found at https://github.com/maRce10/call_degradation_l_palavanensis

Purpose

Measure degradation on re-recorded files
Run stats

Analysis flowchart

flowchart
  A[Time sync files] --> B(Measure degradation) 
  B --> C(Regression models)
  C --> D(Model selection)

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("knitr", "formatR", "baRulho",
    "viridis", "warbleR", "Rraven", "brms", "ggplot2", "corrplot",
    "emmeans", "ggsignif", "lme4", "brmsish"))

Warning: replacing previous import 'brms::rstudent_t' by 'ggdist::rstudent_t'
when loading 'brmsish'

Warning: replacing previous import 'brms::dstudent_t' by 'ggdist::dstudent_t'
when loading 'brmsish'

Warning: replacing previous import 'brms::qstudent_t' by 'ggdist::qstudent_t'
when loading 'brmsish'

Warning: replacing previous import 'brms::pstudent_t' by 'ggdist::pstudent_t'
when loading 'brmsish'

1 Time sync all replicates

1.1 Females

Code

master_annotations_males <- read.csv("./data/raw/master_l.palavanensis_males_annotations.csv")

master_annotations <- rbind(master_annotations_females, master_annotations_males)


exp_raven(master_annotations, path = "./data/raw", file.name = "pooled_master_annotations",
    sound.file.path = "./data/raw/recordings")

female_files <- list.files("./data/raw/recordings", pattern = "_F_")

female_starts <- find_markers(X = master_annotations_females, markers = c("start_marker",
    "end_marker"), path = "./data/raw/recordings", cores = 3, test.files = female_files)

female_starts$start[female_starts$sound.files == "PB_Lpa_0016_T1_Ashton_F_16m.wav" &
    female_starts$marker == "start_marker"] <- 0.8659
female_starts$end[female_starts$sound.files == "PB_Lpa_0016_T1_Ashton_F_16m.wav" &
    female_starts$marker == "start_marker"] <- 1.8659
female_starts$scores[female_starts$sound.files == "PB_Lpa_0016_T1_Ashton_F_16m.wav" &
    female_starts$marker == "start_marker"] <- 1

female_starts$start[female_starts$sound.files == "PB_Lpa_0062_T1_Ashton_F_20m.wav" &
    female_starts$marker == "start_marker"] <- 1.6465
female_starts$end[female_starts$sound.files == "PB_Lpa_0062_T1_Ashton_F_20m.wav" &
    female_starts$marker == "start_marker"] <- 2.6465
female_starts$scores[female_starts$sound.files == "PB_Lpa_0062_T1_Ashton_F_20m.wav" &
    female_starts$marker == "start_marker"] <- 1


# female_starts$sound.files[(abs(female_starts$time.mismatch) >
# 1) & !is.na(female_starts$time.mismatch)]

warbleR::info_sound_files("./data/raw/recordings")

alg_females <- align_test_files(X = master_annotations_females, Y = female_starts,
    path = "./data/raw/recordings", by.song = FALSE)

alg_females$row <- 1:nrow(alg_females)

exp_raven(alg_females, path = "./data/processed", file.name = "check_alignment_females",
    sound.file.path = "./data/raw/recordings", single.file = TRUE)

# getOption('baRulho')$files_to_check_align_test_files

cs <- check_sels(alg_females, path = "./data/raw/recordings")

alg_females <- manual_realign(X = alg_females, Y = master_annotations_females,
    path = "./data/raw/recordings", flim = c(0, 6), marker = "end_marker")

1.2 Males

Code

master_annotations_males <- read.csv("./data/raw/master_l.palavanensis_males_annotations.csv")

male_files <- list.files("./data/raw/recordings", pattern = "_M_")

male_starts <- find_markers(X = master_annotations_males, markers = c("start_marker",
    "end_marker"), path = "./data/raw/recordings", cores = 3, test.files = male_files)

# male_starts$start[male_starts$sound.files ==
# 'PB_Lpa_0019_T1_Apan_M_1m.wav' & male_starts$marker ==
# 'start_marker'] <- 2.8459
# male_starts$end[male_starts$sound.files ==
# 'PB_Lpa_0019_T1_Apan_M_1m.wav' & male_starts$marker ==
# 'start_marker'] <- 3.8459
# male_starts$scores[male_starts$sound.files ==
# 'PB_Lpa_0019_T1_Apan_M_1m.wav' & male_starts$marker ==
# 'start_marker'] <- 1


# male_starts[(abs(male_starts$time.mismatch) > 1) &
# !is.na(male_starts$time.mismatch),]

alg_males <- align_test_files(X = master_annotations_males, Y = male_starts,
    path = "./data/raw/recordings", by.song = FALSE)

alg_females$row <- 1:nrow(alg_females)

exp_raven(alg_males, path = "./data/processed", file.name = "check_alignment_males",
    sound.file.path = "./data/raw/recordings", single.file = TRUE)

# getOption('baRulho')$files_to_check_align_test_files

cs <- check_sels(alg_females, path = "./data/raw/recordings")

alg_males <- manual_realign(X = alg_males, Y = master_annotations_males,
    path = "./data/raw/recordings", flim = c(0, 6), marker = "end_marker")

alg_females$sex <- "females"
alg_females$row <- NULL
alg_males$sex <- "males"

alg_anns <- rbind(alg_females, alg_males)

alg_anns$distance <- sapply(strsplit(alg_anns$sound.files, "_"), "[[",
    7)
alg_anns$distance <- gsub(".wav", "", alg_anns$distance)
alg_anns <- alg_anns[alg_anns$distance != "10m", ]

alg_anns$site <- sapply(strsplit(alg_anns$sound.files, "_"), "[[",
    5)
alg_anns$transect <- sapply(strsplit(alg_anns$sound.files, "_"), "[[",
    4)
alg_anns$transect <- sapply(strsplit(alg_anns$sound.files, "_"), "[[",
    4)

# keep only 1 reference
alg_anns <- alg_anns[!(alg_anns$distance == "1m" & alg_anns$site ==
    "Apan"), ]

table(alg_anns$distance, alg_anns$transect, alg_anns$sex)

alg_anns <- alg_anns[grep("marker", alg_anns$sound.id, invert = TRUE),
    ]

write.csv(alg_anns, "./data/processed/pooled_aligned_annotations.csv",
    row.names = FALSE)

alg_anns_est <- selection_table(alg_anns, path = "./data/raw/recordings",
    extended = TRUE)

alg_anns_est <- resample_est(alg_anns_est, samp.rate = 22.05)

saveRDS(alg_anns_est, "./data/raw/extended_selection_table_rerecorded_sounds.RDS")

2 Measure degradation

2.1 Measure

Code

alg_anns_est <- readRDS("../data/processed/extended_selection_table_rerecorded_sounds.RDS")

table(alg_anns_est$sound.id)

table(alg_anns_est$distance)

alg_anns_est$distance <- as.numeric(gsub("m", "", alg_anns_est$distance))
cores <- 3

alg_anns_est <- set_reference_sounds(alg_anns_est)


# run blur ratio
alg_anns_est <- blur_ratio(alg_anns_est, cores = cores)

# run Spectrum blur ratio
alg_anns_est <- spectrum_blur_ratio(alg_anns_est, cores = cores)

# run envelope correlation
alg_anns_est <- excess_attenuation(alg_anns_est, cores = cores)

# run envelope correlation
alg_anns_est <- envelope_correlation(alg_anns_est, cores = cores)

# run spectrum correlation
alg_anns_est <- spectrum_correlation(alg_anns_est, cores = cores)

# run signal to noise ratio
alg_anns_est <- signal_to_noise_ratio(alg_anns_est, cores = cores,
    mar = 0.03)

# run tail to noise ratio
alg_anns_est <- tail_to_signal_ratio(alg_anns_est, cores = cores,
    tsr.formula = 2, mar = 0.03)

names(alg_anns_est)[ncol(alg_anns_est)] <- "tail.to.noise.ratio"

# run tail to signal ratio
alg_anns_est <- tail_to_signal_ratio(alg_anns_est, cores = cores,
    tsr.formula = 1, mar = 0.03)

# run spcc
source("~/Dropbox/R_package_testing/baRulho/R/spcc.R")
source("~/Dropbox/R_package_testing/warbleR/R/cross_correlation.R")
source("~/Dropbox/R_package_testing/warbleR/R/internal_functions.R")

alg_anns_est <- spcc(X = alg_anns_est, cores = cores)

alg_anns_est <- alg_anns_est[alg_anns_est$distance != 1, ]

alg.tests <- as.data.frame(alg_anns_est)

write.csv(alg.tests, "./data/processed/degradation_metrics.csv", row.names = FALSE)

3 Stats

Code

degrad_dat <- read.csv("./data/processed/degradation_metrics.csv")

degrad_measures <- c("blur.ratio", "spectrum.blur.ratio", "excess.attenuation",
    "envelope.correlation", "spectrum.correlation", "signal.to.noise.ratio",
    "tail.to.noise.ratio", "tail.to.signal.ratio", "cross.correlation")

pca <- prcomp(degrad_dat[, degrad_measures], scale = TRUE)

pca
summary(pca)

degrad_dat$degrad.pc1 <- pca$x[, 1]
degrad_dat$original.sound.file <- sapply(strsplit(degrad_dat$sound.id,
    ".WAV"), "[[", 1)

degrad_dat$distance_f <- factor(degrad_dat$distance, ordered = TRUE)

library(lmerTest)
int_mod <- lmer(degrad.pc1 ~ sex * site + (1 | transect) + (1 | original.sound.file) +
    (1 | distance_f), data = degrad_dat)

summary(int_mod)

int_mod2 <- lmer(degrad.pc1 ~ sex * site + distance_f + (1 | transect) +
    (1 | original.sound.file), data = degrad_dat)

summary(int_mod2)

iter <- 5000
chains <- 4
priors <- c(prior(normal(0, 6), class = "b"))

mod <- brm(degrad.pc1 ~ sex * site + mo(distance_f) + (1 | transect) +
    (1 | original.sound.file), data = degrad_dat, prior = priors,
    iter = iter, chains = chains, cores = chains, control = list(adapt_delta = 0.99,
        max_treedepth = 15), file = "./data/processed/interaction_model.RDS",
    file_refit = "always", seed = 123)

Code

mod <- readRDS("./data/processed/interaction_model.RDS")

extended_summary(mod, trace.palette = mako, highlight = TRUE, remove.intercepts = TRUE)

3.1 mod

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 6) Intercept-student_t(3, 0.3, 2.9) sd-student_t(3, 0, 2.9) sigma-student_t(3, 0, 2.9) simo-dirichlet(1)	degrad.pc1 ~ sex * site + mo(distance_f) + (1 \| transect) + (1 \| original.sound.file)	5000	4	1	2500	0 (0%)	0	1995.689	3336.118	123

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_sexmales	-0.498	-0.795	-0.196	1	1995.689	3336.118
b_siteAshton	0.088	0.010	0.164	1	7268.361	6994.519
b_sexmales:siteAshton	0.105	-0.010	0.218	1	7194.042	7474.724

3.1.1 Descriptive stats

Code

# Total test sounds: `r nrow(degrad_dat)`

agg <- aggregate(selec ~ habitat.type + distance, degrad_dat, length)

agg$bytran <- round(agg$selec/3, 0)

names(agg) <- c("habitat type", "distance", "total test sounds", "sounds per transect")

agg

Takeaways

Session information

R version 4.3.2 (2023-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.2 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

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.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] brmsish_1.0.0      lme4_1.1-35.3      Matrix_1.6-5       ggsignif_0.6.4    
 [5] emmeans_1.9.0      corrplot_0.92      ggplot2_3.5.1      brms_2.21.0       
 [9] Rcpp_1.0.12        Rraven_1.0.13      viridis_0.6.5      viridisLite_0.4.2 
[13] baRulho_2.1.1      ohun_1.0.2         warbleR_1.1.30     NatureSounds_1.0.4
[17] seewave_2.2.3      tuneR_1.4.7        formatR_1.14       knitr_1.48        

loaded via a namespace (and not attached):
  [1] tensorA_0.36.2.1     rstudioapi_0.16.0    jsonlite_1.8.8      
  [4] magrittr_2.0.3       TH.data_1.1-2        sketchy_1.0.3       
  [7] estimability_1.4.1   farver_2.1.2         nloptr_2.0.3        
 [10] rmarkdown_2.27       vctrs_0.6.5          minqa_1.2.7         
 [13] RCurl_1.98-1.14      htmltools_0.5.8.1    distributional_0.4.0
 [16] curl_5.2.0           tidybayes_3.0.6      signal_1.8-1        
 [19] StanHeaders_2.32.9   KernSmooth_2.23-20   htmlwidgets_1.6.4   
 [22] plyr_1.8.9           sandwich_3.1-0       testthat_3.2.1.1    
 [25] zoo_1.8-12           igraph_2.0.3         lifecycle_1.0.4     
 [28] pkgconfig_2.0.3      R6_2.5.1             fastmap_1.2.0       
 [31] digest_0.6.36        colorspace_2.1-0     labeling_0.4.3      
 [34] fansi_1.0.6          abind_1.4-5          compiler_4.3.2      
 [37] proxy_0.4-27         remotes_2.5.0        withr_3.0.0         
 [40] backports_1.5.0      inline_0.3.19        DBI_1.2.3           
 [43] QuickJSR_1.2.2       pkgbuild_1.4.4       MASS_7.3-55         
 [46] rjson_0.2.21         classInt_0.4-10      loo_2.7.0           
 [49] tools_4.3.2          units_0.8-5          ape_5.8             
 [52] fftw_1.0-8           glue_1.7.0           nlme_3.1-155        
 [55] grid_4.3.2           sf_1.0-16            checkmate_2.3.1     
 [58] reshape2_1.4.4       generics_0.1.3       gtable_0.3.5        
 [61] class_7.3-20         tidyr_1.3.1          xml2_1.3.6          
 [64] Deriv_4.1.3          utf8_1.2.4           pillar_1.9.0        
 [67] ggdist_3.3.2         stringr_1.5.1        posterior_1.5.0     
 [70] splines_4.3.2        dplyr_1.1.4          lattice_0.20-45     
 [73] survival_3.2-13      tidyselect_1.2.1     pbapply_1.7-2       
 [76] arrayhelpers_1.1-0   gridExtra_2.3        V8_4.4.2            
 [79] svglite_2.1.3        stats4_4.3.2         xfun_0.45           
 [82] bridgesampling_1.1-2 brio_1.1.5           matrixStats_1.3.0   
 [85] rstan_2.32.6         stringi_1.8.4        yaml_2.3.9          
 [88] boot_1.3-28          kableExtra_1.4.0     xaringanExtra_0.7.0 
 [91] evaluate_0.24.0      codetools_0.2-18     dtw_1.23-1          
 [94] tibble_3.2.1         cli_3.6.3            RcppParallel_5.1.7  
 [97] xtable_1.8-4         Sim.DiffProc_4.9     systemfonts_1.1.0   
[100] munsell_0.5.1        coda_0.19-4.1        png_0.1-8           
[103] svUnit_1.0.6         parallel_4.3.2       rstantools_2.4.0    
[106] bayesplot_1.11.1     Brobdingnag_1.2-9    bitops_1.0-7        
[109] mvtnorm_1.2-5        scales_1.3.0         e1071_1.7-14        
[112] packrat_0.9.2        purrr_1.0.2          crayon_1.5.3        
[115] rlang_1.1.4          cowplot_1.1.3        multcomp_1.4-25