Echolocation plasticity

Author

Marcelo Araya-Salas

Published

January 17, 2026

1 Purpose

  • The first goal of this report

  • The second goal of this report

 

2 Report overview

 

3 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", "warbleR", "Rraven", "ohun", "readxl", "brms", "brmsish", "viridis", "ggplot2"))
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'
Code 
my.viridis <- function(...) mako(alpha = 0.5, begin = 0.3, end = 0.7, ...)


fill_color <- my.viridis(10)[3]

options(brms.file_refit = "always", mc.cores = 2)

theme_set(theme(text = element_text(size = 20)))

4 Clipping recordings

Code 
sound_file_folder <- "/media/marce/PBD (A)/murcielagos_la_selva_2026/4_channel_sound_files/"

# leer los archivos de raven
anns <- imp_raven(path = sound_file_folder, warbler.format = TRUE, all.data = TRUE, name.from.file = TRUE, ext.case = "lower")

feature_reference(anns, units = c("ms", "kHz"))

cut_sels(anns, path = sound_file_folder, dest.path = "./data/raw/clips_2026", mar = 0.003)

warbleR_options(wav.path = "./data/raw/clips_2026")

full_spectrograms(dest.path = "./data/processed/spectrograms_2026", flim = c(20, 120), sxrow = 0.05, rows = 8, fast.spec = FALSE, ovlp = 90, X = NULL, horizontal = TRUE, parallel = 4, collevels = seq(-100, 0, 5), suffix = "0.25s")

anns <- anns[order(anns$sound.files, anns$start), ]


est_bats <- selection_table(anns, extended = TRUE, mar = 0.002, path = sound_file_folder)

saveRDS(est_bats, "./data/processed/est_bats_2026.RDS")

5 Importing annotations

Code 
est_bats <-  readRDS("./data/processed/est_bats_2026.RDS")

sp_feat <- spectro_analysis(est_bats, ovlp = 70, wl = 200, wl.freq = 512)


est_bats$freq <- sp_feat$meanpeakf
est_bats$duration <- sp_feat$time.IQR

est_bats$pc1 <- prcomp(sp_feat[, -c(1:2)], scale. = TRUE)$x[, 1]

dat_time <- attributes(est_bats)$check.results

est_bats$gaps <- NA   
for(i in 2:nrow(est_bats)){
    if(dat_time$orig.sound.files[i] == dat_time$orig.sound.files[i-1]){
        est_bats$gaps[i] <- dat_time$orig.start[i] - dat_time$orig.start[i-1]
    } else {
        est_bats$gaps[i] <- NA
    }
}

est_bats$gaps[est_bats$gaps > 30] <- NA

est_bats$rate <- 1 / est_bats$gaps

est_bats$orig.sf <- dat_time$orig.sound.files
est_bats$orig.start <- dat_time$orig.start

metadata <- read_xlsx("./data/raw/datos_ecolocalizacion_murcis_la_selva.xlsx")

metadata <- metadata[complete.cases(metadata), ]

est_bats$sp <- sapply(est_bats$orig.sf, function(x)
    metadata$especie[metadata$`nombre archivo` == x][1])
    
est_bats$indiv <- sapply(est_bats$orig.sf, function(x)
    metadata$individuo[metadata$`nombre archivo` == x][1])

est_bats$indiv <- paste0("ind", est_bats$indiv)

est_bats <- est_bats[!is.na(est_bats$sp), ]


est_bats$treatment <- sapply(est_bats$orig.sf, function(x)
    metadata$tratamiento[metadata$`nombre archivo` == x][1])


est_bats$guild <- ifelse(est_bats$sp == "M. nigricans", "edge", "narrow")

saveRDS(est_bats, "./data/processed/acoustic_analysis_data_2026.RDS")

6 Statistical analysis

6.1 Predicciones

Code 
sim <- data.frame(
    guild = rep(c("Cerrado", "Abierto"), each = 100),
    duration = c(rnorm(100, mean = 0.03, sd = 0.01), rnorm(100, mean = 0.05, sd = 0.01)))


ggplot(sim, aes(x = guild, y = duration)) + 
    geom_boxplot(fill = mako(10, alpha = 0.3)[4]) +
    labs(x = "Ambiente", y = "Duración") +
    # remove y axis tick labels
    theme_classic(base_size = 30) + theme(axis.text.y=element_blank(),
        axis.ticks.y=element_blank() 
        )

Code 
sim2 <- data.frame(
    guild = rep(c("Cerrado", "Abierto"), each = 100),
    rate = c(rnorm(100, mean = 0.05, sd = 0.01), rnorm(100, mean = 0.032, sd = 0.01)))


ggplot(sim2, aes(x = guild, y = rate)) + 
    geom_boxplot(fill = mako(10, alpha = 0.3)[4]) +
    labs(x = "Ambiente", y = "Tasa de llamados") +
    # remove y axis tick labels
    theme_classic(base_size = 30) + theme(axis.text.y=element_blank(),
        axis.ticks.y=element_blank() 
        )

Code 
sim3 <- data.frame(
    guild = rep(c("Cerrado", "Abierto"), each = 100),
    freq = c(rnorm(100, mean = 0.05, sd = 0.01), rnorm(100, mean = 0.032, sd = 0.01)))


ggplot(sim3, aes(x = guild, y = freq)) + 
    geom_boxplot(fill = mako(10, alpha = 0.3)[4]) +
    labs(x = "Ambiente", y = "Frecuencia (kHz)") +
    # remove y axis tick labels
    theme_classic(base_size = 30) + theme(axis.text.y=element_blank(),
        axis.ticks.y=element_blank() 
        )

6.2 Macroevolution

6.2.1 By guild

Code 
est_bats <- readRDS("./data/processed/acoustic_analysis_data_2026.RDS")

iter <- 5000
chains <- 2

priors <- c(
  # Population-level effects (including intercept)
  # For log-link Gamma, intercept is on log scale
  set_prior("normal(1, 2)", class = "Intercept"),  # log-scale: exp(1) ≈ 2.7 mean duration
  set_prior("normal(0, 1)", class = "b"),
  
  # Random effect standard deviations
  set_prior("exponential(1)", class = "sd", group = "indiv"),
  set_prior("exponential(1)", class = "sd", group = "sp"),
  
  # Gamma shape parameter (nu or shape)
  # Constrained to be positive
  set_prior("exponential(1)", class = "shape")
)

dur_mod <- brm(
        formula = duration ~ guild + (1 |
                                          indiv) + (1 |
                                                        sp) + ar(p = 2, time = orig.start, gr = orig.sf) + (1 | treatment),
        family = Gamma(),
        data = est_bats,
        prior = priors,
        chains = chains,
        iter = iter,
        file = "./data/processed/duration_macro_model.RDS"
    )

priors <- c(
  # Population-level (fixed) effects
  set_prior("normal(0, 1)", class = "b"),
  
  # Intercept - centered around plausible frequency value
  # Adjust mean based on your data (e.g., kHz for bat calls)
  set_prior("normal(0, 10)", class = "Intercept"),
  
  # Random effect standard deviations
  set_prior("exponential(1)", class = "sd", group = "indiv"),
  set_prior("exponential(1)", class = "sd", group = "sp"),
  
  # Residual standard deviation
  set_prior("exponential(1)", class = "sigma")
)

freq_mod <- brm(
        formula = freq ~ guild + (1 |
                                          indiv) + (1 |
                                                        sp) + ar(p = 2, time = orig.start, gr = orig.sf) + (1 | treatment),
        family = gaussian(),
        data = est_bats,
        prior = priors,
        chains = chains,
        iter = iter,
        file = "./data/processed/freq_macro_model.RDS"
    )

priors <- c(
  # Population-level effects (including intercept)
  # For log-link Gamma, intercept is on log scale
  set_prior("normal(1, 2)", class = "Intercept"),  # log-scale: exp(1) ≈ 2.7 mean duration
  set_prior("normal(0, 1)", class = "b"),
  
  # Random effect standard deviations
  set_prior("exponential(1)", class = "sd", group = "indiv"),
  set_prior("exponential(1)", class = "sd", group = "sp"),
  
  # Gamma shape parameter (nu or shape)
  # Constrained to be positive
  set_prior("exponential(1)", class = "shape")
)

rate_mod <- brm(
        formula = rate ~ guild + (1 |
                                          indiv) + (1 |
                                                        sp) + ar(p = 2, time = orig.start, gr = orig.sf) + (1 | treatment),
        family = Gamma(),
        data = est_bats,
        prior = priors,
        chains = chains,
        iter = iter,
        file = "./data/processed/rate_macro_model.RDS"
    )


priors <- c(
  # Population-level effects
  # PCA scores usually range ~ -3 to 3
  set_prior("normal(0, 1)", class = "b"),
  
  # Intercept - centered at 0 (mean of PCA scores)
  set_prior("normal(0, 1)", class = "Intercept"),
  
  # Random effect standard deviations
  # Expect smaller variation than total variance (SD=1)
  set_prior("exponential(2)", class = "sd", group = "indiv"),
  set_prior("exponential(2)", class = "sd", group = "sp"),
  
  # Residual standard deviation
  # Should be less than 1 (total SD of PCA scores)
  set_prior("exponential(2)", class = "sigma")
)

pc1_mod <- brm(
        formula = pc1 ~ guild + (1 |
                                          indiv) + (1 |
                                                        sp) + ar(p = 2, time = orig.start, gr = orig.sf) + (1 | treatment),
        data = est_bats,
        prior = priors,
        chains = chains,
        iter = iter,
        file = "./data/processed/pc1_macro_model.RDS"
    )

Duration:

Code 
extended_summary(read.file = "./data/processed/duration_macro_model.RDS",
    n.posterior = 1000, fill = fill_color, trace.palette = my.viridis, remove.intercepts = FALSE, highlight = TRUE, print.name = FALSE)
formula family priors iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 duration ~ guild + (1 | indiv) + (1 | sp) + ar(p = 2, time = orig.start, gr = orig.sf) + (1 | treatment) gamma (inverse) b-normal(0, 1) Intercept-normal(1, 2) sd-student_t(3, 0, 327.1) sd-exponential(1) sd-exponential(1) sderr-student_t(3, 0, 327.1) shape-exponential(1) 5000 1 1 2500 627 (0.25%) 0 595.178 1246.381 225527250
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept 1.194 -3.206 5.317 1.004 595.178 1246.381
b_guildnarrow 0.002 -1.970 1.876 1.006 1615.282 1325.867

Frequency:

Code 
extended_summary(read.file = "./data/processed/freq_macro_model.RDS",
    n.posterior = 1000, fill = fill_color, trace.palette = my.viridis, remove.intercepts = FALSE, highlight = TRUE, print.name = FALSE)
formula family priors iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 freq ~ guild + (1 | indiv) + (1 | sp) + ar(p = 2, time = orig.start, gr = orig.sf) + (1 | treatment) gaussian (identity) b-normal(0, 1) Intercept-normal(0, 10) sd-student_t(3, 0, 16.6) sd-exponential(1) sd-exponential(1) sigma-exponential(1) 5000 2 1 2500 2 (4e-04%) 0 6167.505 3292.26 1360243337
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept 7.756 -12.702 28.143 1 6167.505 3903.956
b_guildnarrow -0.094 -2.025 1.838 1.001 7452.598 3292.260

Rate:

Code 
extended_summary(read.file = "./data/processed/rate_macro_model.RDS",
    n.posterior = 1000, fill = fill_color, trace.palette = my.viridis, remove.intercepts = FALSE, highlight = TRUE, print.name = FALSE)
formula family priors iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 rate ~ guild + (1 | indiv) + (1 | sp) + ar(p = 2, time = orig.start, gr = orig.sf) + (1 | treatment) gamma (inverse) b-normal(0, 1) Intercept-normal(1, 2) sd-student_t(3, 0, 2.5) sd-exponential(1) sd-exponential(1) sderr-student_t(3, 0, 2.5) shape-exponential(1) 5000 2 1 2500 1321 (0.26%) 0 456.934 139.856 1225555643
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept 0.070 -0.046 0.180 1.007 456.934 139.856
b_guildnarrow -0.002 -0.127 0.117 1.004 457.989 227.569

PC1:

Code 
extended_summary(read.file = "./data/processed/pc1_macro_model.RDS",
    n.posterior = 1000, fill = fill_color, trace.palette = my.viridis, remove.intercepts = FALSE, highlight = TRUE, print.name = FALSE)
formula family priors iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 pc1 ~ guild + (1 | indiv) + (1 | sp) + ar(p = 2, time = orig.start, gr = orig.sf) + (1 | treatment) gaussian (identity) b-normal(0, 1) Intercept-normal(0, 1) sd-student_t(3, 0, 2.8) sd-exponential(2) sd-exponential(2) sigma-exponential(2) 5000 2 1 2500 6 (0.0012%) 0 3709.555 3411.262 1688528368
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept -0.157 -2.336 2.117 1 3709.555 3411.262
b_guildnarrow 0.176 -1.612 2.004 1 4938.247 3908.156

6.2.2 By species

Code 
est_bats <- readRDS("./data/processed/acoustic_analysis_data_2026.RDS")

iter <- 5000
chains <- 2

priors <- c(
  # Population-level effects (including intercept)
  # For log-link Gamma, intercept is on log scale
  set_prior("normal(1, 2)", class = "Intercept"),  # log-scale: exp(1) ≈ 2.7 mean duration
  set_prior("normal(0, 1)", class = "b"),
  
  # Random effect standard deviations
  set_prior("exponential(1)", class = "sd", group = "indiv"),
  set_prior("exponential(1)", class = "sd", group = "sp"),
  
  # Gamma shape parameter (nu or shape)
  # Constrained to be positive
  set_prior("exponential(1)", class = "shape")
)

dur_sp_mod <- brm(
        formula = duration ~ sp + (1 |
                                          indiv) + (1 |
                                                        sp) +  (1 | treatment),
        family = Gamma(),
        data = est_bats,
        prior = priors,
        chains = chains,
        iter = iter,
        file = "./data/processed/duration_sp_macro_model.RDS"
    )

priors <- c(
  # Population-level (fixed) effects
  set_prior("normal(0, 1)", class = "b"),
  
  # Intercept - centered around plausible frequency value
  # Adjust mean based on your data (e.g., kHz for bat calls)
  set_prior("normal(0, 10)", class = "Intercept"),
  
  # Random effect standard deviations
  set_prior("exponential(1)", class = "sd", group = "indiv"),
  set_prior("exponential(1)", class = "sd", group = "sp"),
  
  # Residual standard deviation
  set_prior("exponential(1)", class = "sigma")
)

freq_mod <- brm(
        formula = freq ~ sp + (1 |
                                          indiv) + (1 |
                                                        sp) +  (1 | treatment),
        family = gaussian(),
        data = est_bats,
        prior = priors,
        chains = chains,
        iter = iter,
        file = "./data/processed/freq_sp_macro_model.RDS"
    )

priors <- c(
  # Population-level effects (including intercept)
  # For log-link Gamma, intercept is on log scale
  set_prior("normal(1, 2)", class = "Intercept"),  # log-scale: exp(1) ≈ 2.7 mean duration
  set_prior("normal(0, 1)", class = "b"),
  
  # Random effect standard deviations
  set_prior("exponential(1)", class = "sd", group = "indiv"),
  set_prior("exponential(1)", class = "sd", group = "sp"),
  
  # Gamma shape parameter (nu or shape)
  # Constrained to be positive
  set_prior("exponential(1)", class = "shape")
)

rate_mod <- brm(
        formula = rate ~ sp + (1 |
                                          indiv) + (1 |
                                                        sp) +  (1 | treatment),
        family = Gamma(),
        data = est_bats,
        prior = priors,
        chains = chains,
        iter = iter,
        file = "./data/processed/rate_sp_macro_model.RDS"
    )


priors <- c(
  # Population-level effects
  # PCA scores usually range ~ -3 to 3
  set_prior("normal(0, 1)", class = "b"),
  
  # Intercept - centered at 0 (mean of PCA scores)
  set_prior("normal(0, 1)", class = "Intercept"),
  
  # Random effect standard deviations
  # Expect smaller variation than total variance (SD=1)
  set_prior("exponential(2)", class = "sd", group = "indiv"),
  set_prior("exponential(2)", class = "sd", group = "sp"),
  
  # Residual standard deviation
  # Should be less than 1 (total SD of PCA scores)
  set_prior("exponential(2)", class = "sigma")
)

pc1_mod <- brm(
        formula = pc1 ~ sp + (1 |
                                          indiv) + (1 |
                                                        sp) + (1 | treatment),
        data = est_bats,
        prior = priors,
        chains = chains,
        iter = iter,
        file = "./data/processed/pc1_sp_macro_model.RDS"
    )

Duration:

Code 
extended_summary(read.file = "./data/processed/duration_sp_macro_model.RDS",
    n.posterior = 1000, fill = fill_color, trace.palette = my.viridis, remove.intercepts = FALSE, highlight = TRUE, print.name = FALSE)
formula family priors iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 duration ~ sp + (1 | indiv) + (1 | sp) + (1 | treatment) gamma (inverse) b-normal(0, 1) Intercept-normal(1, 2) sd-student_t(3, 0, 327.1) sd-exponential(1) sd-exponential(1) shape-exponential(1) 5000 2 1 2500 2 (4e-04%) 0 7573.91 3277.513 1638628191
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept 1.008 -3.120 5.177 1 8299.987 3570.739
b_spC.castanea 0.032 -1.896 2.023 1 8801.910 3840.320
b_spC.sowelli 0.031 -1.916 1.961 1.002 8518.652 3555.365
b_spM.nigricans -0.015 -1.983 1.958 1 7573.910 3574.675
b_spP.discolor 0.051 -1.930 2.058 1 8700.270 3277.513

Frequency:

Code 
extended_summary(read.file = "./data/processed/freq_sp_macro_model.RDS",
    n.posterior = 1000, fill = fill_color, trace.palette = my.viridis, remove.intercepts = FALSE, highlight = TRUE, print.name = FALSE)
formula family priors iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 freq ~ sp + (1 | indiv) + (1 | sp) + ar(p = 2, time = orig.start, gr = orig.sf) + (1 | treatment) gaussian (identity) b-normal(0, 1) Intercept-normal(0, 10) sd-student_t(3, 0, 16.6) sd-exponential(1) sd-exponential(1) sigma-exponential(1) 5000 2 1 2500 3 (6e-04%) 0 5849.45 3551.522 301751998
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept 7.577 -13.413 28.868 1 5849.450 3556.663
b_spC.castanea 0.127 -1.762 2.059 1.002 8521.845 3855.300
b_spC.sowelli -0.038 -1.912 1.871 1 8166.803 3551.522
b_spM.nigricans 0.106 -1.880 2.048 1.001 10965.464 3904.477
b_spP.discolor -0.351 -2.286 1.593 1.001 7615.628 3668.233

Rate:

Code 
extended_summary(read.file = "./data/processed/rate_sp_macro_model.RDS",
    n.posterior = 1000, fill = fill_color, trace.palette = my.viridis, remove.intercepts = FALSE, highlight = TRUE, print.name = FALSE)
formula family priors iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 rate ~ sp + (1 | indiv) + (1 | sp) + (1 | treatment) gamma (inverse) b-normal(0, 1) Intercept-normal(1, 2) sd-student_t(3, 0, 2.5) sd-exponential(1) sd-exponential(1) shape-exponential(1) 5000 2 1 2500 2427 (0.49%) 0 124.242 123.226 1972908817
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept 0.055 -0.340 0.382 1.02 139.228 123.226
b_spC.castanea 0.061 -0.478 0.587 1.003 192.114 173.191
b_spC.sowelli 0.007 -0.507 0.610 1.01 124.242 138.404
b_spM.nigricans 0.006 -0.582 0.526 1.007 165.680 184.588
b_spP.discolor 0.003 -0.433 0.411 1.018 151.865 137.235

PC1:

Code 
extended_summary(read.file = "./data/processed/pc1_sp_macro_model.RDS",
    n.posterior = 1000, fill = fill_color, trace.palette = my.viridis, remove.intercepts = FALSE, highlight = TRUE, print.name = FALSE)
formula family priors iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 pc1 ~ sp + (1 | indiv) + (1 | sp) + (1 | treatment) gaussian (identity) b-normal(0, 1) Intercept-normal(0, 1) sd-student_t(3, 0, 2.8) sd-exponential(2) sd-exponential(2) sigma-exponential(2) 5000 2 1 2500 13 (0.0026%) 0 3034.276 2809.475 1649834349
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept 0.034 -1.734 1.743 1.001 3882.995 2890.248
b_spC.castanea -1.102 -2.775 0.796 1.003 3593.872 2809.475
b_spC.sowelli -0.074 -1.631 1.564 1 5465.547 3489.570
b_spM.nigricans -0.300 -1.929 1.395 1.002 5860.231 3629.612
b_spP.discolor 1.583 -0.544 3.552 1 3034.276 3243.044

6.3 Plasticity

6.3.1 Only habitat

Code 
priors <- c(
  # Population-level effects (including intercept)
  # For log-link Gamma, intercept is on log scale
  set_prior("normal(1, 2)", class = "Intercept"),  # log-scale: exp(1) ≈ 2.7 mean duration
  set_prior("normal(0, 1)", class = "b"),
  
  # Random effect standard deviations
  set_prior("exponential(1)", class = "sd", group = "indiv"),
  set_prior("exponential(1)", class = "sd", group = "sp"),
  
  # Gamma shape parameter (nu or shape)
  # Constrained to be positive
  set_prior("exponential(1)", class = "shape")
)

dur_plast_mod <- brm(
        formula = duration ~ treatment + (1 |
                                          indiv) + (1 |
                                                        sp) + ar(p = 2, time = orig.start, gr = orig.sf),
        family = Gamma(),
        data = est_bats,
        prior = priors,
        chains = chains,
        iter = iter,
        file = "./data/processed/duration_plasticity_model.RDS")

priors <- c(
  # Population-level (fixed) effects
  set_prior("normal(0, 1)", class = "b"),
  
  # Intercept - centered around plausible frequency value
  # Adjust mean based on your data (e.g., kHz for bat calls)
  set_prior("normal(0, 10)", class = "Intercept"),
  
  # Random effect standard deviations
  set_prior("exponential(1)", class = "sd", group = "indiv"),
  set_prior("exponential(1)", class = "sd", group = "sp"),
  
  # Residual standard deviation
  set_prior("exponential(1)", class = "sigma")
)

freq_plast_mod <- brm(
        formula = freq ~ treatment + (1 |
                                          indiv) + (1 |
                                                        sp) + ar(p = 2, time = orig.start, gr = orig.sf),
        family = gaussian(),
        data = est_bats,
        prior = priors,
        chains = chains,
        iter = iter,
        file = "./data/processed/freq_plasticity_model.RDS"
    )

priors <- c(
  # Population-level effects (including intercept)
  # For log-link Gamma, intercept is on log scale
  set_prior("normal(1, 2)", class = "Intercept"),  # log-scale: exp(1) ≈ 2.7 mean duration
  set_prior("normal(0, 1)", class = "b"),
  
  # Random effect standard deviations
  set_prior("exponential(1)", class = "sd", group = "indiv"),
  set_prior("exponential(1)", class = "sd", group = "sp"),
  
  # Gamma shape parameter (nu or shape)
  # Constrained to be positive
  set_prior("exponential(1)", class = "shape")
)

rate_plast_mod <- brm(
        formula = rate ~ treatment + (1 |
                                          indiv) + (1 |
                                                        sp) + ar(p = 2, time = orig.start, gr = orig.sf),
        family = Gamma(),
        data = est_bats,
        prior = priors,
        chains = chains,
        iter = iter,
        file = "./data/processed/rate_plasticity_model.RDS"
    )


priors <- c(
  # Population-level effects
  # PCA scores usually range ~ -3 to 3
  set_prior("normal(0, 1)", class = "b"),
  
  # Intercept - centered at 0 (mean of PCA scores)
  set_prior("normal(0, 1)", class = "Intercept"),
  
  # Random effect standard deviations
  # Expect smaller variation than total variance (SD=1)
  set_prior("exponential(2)", class = "sd", group = "indiv"),
  set_prior("exponential(2)", class = "sd", group = "sp"),
  
  # Residual standard deviation
  # Should be less than 1 (total SD of PCA scores)
  set_prior("exponential(2)", class = "sigma")
)

pc1_plast_mod <- brm(
        formula = pc1 ~ treatment + (1 |
                                          indiv) + (1 |
                                                        sp) + ar(p = 2, time = orig.start, gr = orig.sf),
        data = est_bats,
        prior = priors,
        chains = chains,
        iter = iter,
        file = "./data/processed/pc1_plasticity_model.RDS"
    )

Duration:

Code 
extended_summary(read.file = "./data/processed/duration_plasticity_model.RDS",
    n.posterior = 1000, fill = fill_color, trace.palette = my.viridis, remove.intercepts = FALSE, highlight = TRUE, print.name = FALSE)
formula family priors iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 duration ~ treatment + (1 | indiv) + (1 | sp) + ar(p = 2, time = orig.start, gr = orig.sf) gamma (inverse) b-normal(0, 1) Intercept-normal(1, 2) sd-student_t(3, 0, 327.1) sd-exponential(1) sd-exponential(1) sderr-student_t(3, 0, 327.1) shape-exponential(1) 5000 2 1 2500 1166 (0.23%) 0 1457.818 1391.327 1813161592
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept 1.798 -2.137 5.847 1.006 1562.714 1391.327
b_treatmentafuera -0.032 -1.887 1.823 1.002 1457.818 1729.567
b_treatmentcerrado 0.014 -1.913 1.925 1.006 2239.856 1465.158

Frequency:

Code 
extended_summary(read.file = "./data/processed/freq_plasticity_model.RDS",
    n.posterior = 1000, fill = fill_color, trace.palette = my.viridis, remove.intercepts = FALSE, highlight = TRUE, print.name = FALSE)
formula family priors iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 freq ~ treatment + (1 | indiv) + (1 | sp) + ar(p = 2, time = orig.start, gr = orig.sf) gaussian (identity) b-normal(0, 1) Intercept-normal(0, 10) sd-student_t(3, 0, 16.6) sd-exponential(1) sd-exponential(1) sigma-exponential(1) 5000 2 1 2500 4 (8e-04%) 0 967.097 1297.08 1188006390
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept 56.527 34.998 72.048 1.001 967.097 1297.080
b_treatmentafuera -0.358 -2.277 1.539 1.001 5740.568 3282.569
b_treatmentcerrado -0.074 -1.899 1.757 1 5652.238 3517.020

Rate:

Code 
extended_summary(read.file = "./data/processed/rate_plasticity_model.RDS",
    n.posterior = 1000, fill = fill_color, trace.palette = my.viridis, remove.intercepts = FALSE, highlight = TRUE, print.name = FALSE)
formula family priors iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 rate ~ treatment + (1 | indiv) + (1 | sp) + ar(p = 2, time = orig.start, gr = orig.sf) gamma (inverse) b-normal(0, 1) Intercept-normal(1, 2) sd-student_t(3, 0, 2.5) sd-exponential(1) sd-exponential(1) sderr-student_t(3, 0, 2.5) shape-exponential(1) 5000 2 1 2500 1379 (0.28%) 0 1309.412 1140.78 1326725731
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept 0.072 0.033 0.111 1.001 1309.412 1140.780
b_treatmentafuera 0.004 -0.023 0.033 1 3607.585 2428.173
b_treatmentcerrado -0.007 -0.023 0.008 1 3409.088 2660.836

PC1:

Code 
extended_summary(read.file = "./data/processed/pc1_plasticity_model.RDS",
    n.posterior = 1000, fill = fill_color, trace.palette = my.viridis, remove.intercepts = FALSE, highlight = TRUE, print.name = FALSE)
formula family priors iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 pc1 ~ treatment + (1 | indiv) + (1 | sp) + ar(p = 2, time = orig.start, gr = orig.sf) gaussian (identity) b-normal(0, 1) Intercept-normal(0, 1) sd-student_t(3, 0, 2.8) sd-exponential(2) sd-exponential(2) sigma-exponential(2) 5000 2 1 2500 0 (0%) 0 3269.256 3238.144 1004510779
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept 0.086 -1.243 1.469 1 3269.256 3509.648
b_treatmentafuera 0.958 -0.403 2.241 1 4817.504 3238.144
b_treatmentcerrado -0.694 -1.520 0.071 1 5927.799 3276.075

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


# contrasts
contrasts <-contrasts(
  fit = mod,
  predictor = "treatment",
  n.posterior = 2000,
  level.sep = " VS ",
  html.table = FALSE,
  plot = TRUE,
  highlight = TRUE,
  fill = fill_color
)

6.3.2 Habitat by guild

Code 
iter <- 10000
chains <- 4

priors <- c(
  # Population-level effects (including intercept)
  # For log-link Gamma, intercept is on log scale
  set_prior("normal(1, 2)", class = "Intercept"),  # log-scale: exp(1) ≈ 2.7 mean duration
  set_prior("normal(0, 1)", class = "b"),
  
  # Random effect standard deviations
  set_prior("exponential(1)", class = "sd", group = "indiv"),
  set_prior("exponential(1)", class = "sd", group = "sp"),
  
  # Gamma shape parameter (nu or shape)
  # Constrained to be positive
  set_prior("exponential(1)", class = "shape")
)

dur_plast_mod <- brm(
        formula = duration ~ treatment * guild + (1 |
                                          indiv) + (1 |
                                                        sp) + ar(p = 2, time = orig.start, gr = orig.sf),
        family = Gamma(),
        data = est_bats,
        prior = priors,
        chains = chains,
        iter = iter,
        file = "./data/processed/duration_plasticity_interaction_model.RDS")

priors <- c(
  # Population-level (fixed) effects
  set_prior("normal(0, 1)", class = "b"),
  
  # Intercept - centered around plausible frequency value
  # Adjust mean based on your data (e.g., kHz for bat calls)
  set_prior("normal(0, 10)", class = "Intercept"),
  
  # Random effect standard deviations
  set_prior("exponential(1)", class = "sd", group = "indiv"),
  set_prior("exponential(1)", class = "sd", group = "sp"),
  
  # Residual standard deviation
  set_prior("exponential(1)", class = "sigma")
)

freq_plast_mod <- brm(
        formula = freq ~ treatment * guild + (1 |
                                          indiv) + (1 |
                                                        sp) + ar(p = 2, time = orig.start, gr = orig.sf),
        family = gaussian(),
        data = est_bats,
        prior = priors,
        chains = chains,
        iter = iter,
        file = "./data/processed/freq_plasticity_interaction_model.RDS"
    )

priors <- c(
  # Population-level effects (including intercept)
  # For log-link Gamma, intercept is on log scale
  set_prior("normal(1, 2)", class = "Intercept"),  # log-scale: exp(1) ≈ 2.7 mean duration
  set_prior("normal(0, 1)", class = "b"),
  
  # Random effect standard deviations
  set_prior("exponential(1)", class = "sd", group = "indiv"),
  set_prior("exponential(1)", class = "sd", group = "sp"),
  
  # Gamma shape parameter (nu or shape)
  # Constrained to be positive
  set_prior("exponential(1)", class = "shape")
)

rate_plast_mod <- brm(
        formula = rate ~ treatment * guild + (1 |
                                          indiv) + (1 |
                                                        sp) + ar(p = 2, time = orig.start, gr = orig.sf),
        family = Gamma(),
        data = est_bats,
        prior = priors,
        chains = chains,
        iter = iter,
        file = "./data/processed/rate_plasticity_interaction_model.RDS"
    )


priors <- c(
  # Population-level effects
  # PCA scores usually range ~ -3 to 3
  set_prior("normal(0, 1)", class = "b"),
  
  # Intercept - centered at 0 (mean of PCA scores)
  set_prior("normal(0, 1)", class = "Intercept"),
  
  # Random effect standard deviations
  # Expect smaller variation than total variance (SD=1)
  set_prior("exponential(2)", class = "sd", group = "indiv"),
  set_prior("exponential(2)", class = "sd", group = "sp"),
  
  # Residual standard deviation
  # Should be less than 1 (total SD of PCA scores)
  set_prior("exponential(2)", class = "sigma")
)

pc1_plast_mod <- brm(
        formula = pc1 ~ treatment * guild + (1 |
                                          indiv) + (1 |
                                                        sp) + ar(p = 2, time = orig.start, gr = orig.sf),
        data = est_bats,
        prior = priors,
        chains = chains,
        iter = iter,
        file = "./data/processed/pc1_plasticity_interaction_model.RDS"
    )

Duration:

Code 
extended_summary(read.file = "./data/processed/duration_plasticity_interaction_model.RDS",
    n.posterior = 1000, fill = fill_color, trace.palette = my.viridis, remove.intercepts = FALSE, highlight = TRUE, print.name = FALSE)

Frequency:

Code 
extended_summary(read.file = "./data/processed/freq_plasticity_interaction_model.RDS",
    n.posterior = 1000, fill = fill_color, trace.palette = my.viridis, remove.intercepts = FALSE, highlight = TRUE, print.name = FALSE)

Rate:

Code 
extended_summary(read.file = "./data/processed/rate_plasticity_interaction_model.RDS",
    n.posterior = 1000, fill = fill_color, trace.palette = my.viridis, remove.intercepts = FALSE, highlight = TRUE, print.name = FALSE)

PC1:

Code 
extended_summary(read.file = "./data/processed/pc1_plasticity_interaction_model.RDS",
    n.posterior = 1000, fill = fill_color, trace.palette = my.viridis, remove.intercepts = FALSE, highlight = TRUE, print.name = FALSE)

7 Plots

7.1 By guild

Code 
est_bats <- readRDS("./data/processed/acoustic_analysis_data_2026.RDS")

est_bats$guild <- ifelse(est_bats$sp == "M. nigricans", "Cerrado", "Abierto")

ggplot(est_bats, aes(x = guild, y = pc1)) + 
    geom_boxplot(fill = mako(10, alpha = 0.3)[4]) +
    labs(x = "Ambiente", y = "PC1 estructura de\nllamados")

Code 
ggplot(est_bats, aes(x = guild, y = freq)) + 
    geom_boxplot(fill = mako(10, alpha = 0.3)[4]) +
    labs(x = "Ambiente", y = "Frequencia pico (kHz)")

Code 
ggplot(est_bats, aes(x = guild, y = duration)) + 
    geom_boxplot(fill = mako(10, alpha = 0.3)[4]) +
    labs(x = "Ambiente", y = "Duración (s)")

Code 
ggplot(est_bats, aes(x = guild, y = rate)) + 
    geom_boxplot(fill = mako(10, alpha = 0.3)[4]) +
    labs(x = "Ambiente", y = "Tasa de llamados")
Warning: Removed 56 rows containing non-finite outside the scale range
(`stat_boxplot()`).

7.2 By treatment

Code 
ggplot(est_bats, aes(x = treatment, y = pc1)) + 
    geom_boxplot(fill = mako(10, alpha = 0.3)[4]) +
    labs(x = "Tratamiento", y = "PC1")

Code 
ggplot(est_bats, aes(x = treatment, y = freq)) + 
    geom_boxplot(fill = mako(10, alpha = 0.3)[4]) +
    labs(x = "Tratamiento", y = "Frecuencia pico (kHz)")

Code 
ggplot(est_bats, aes(x = treatment, y = duration)) + 
    geom_boxplot(fill = mako(10, alpha = 0.3)[4]) +
    labs(x = "Tratamiento", y = "Duración (s)")

Code 
ggplot(est_bats, aes(x = treatment, y = rate)) + 
    geom_boxplot(fill = mako(10, alpha = 0.3)[4]) +
    labs(x = "Tratamiento", y = "Tasa de llamados")
Warning: Removed 56 rows containing non-finite outside the scale range
(`stat_boxplot()`).

7.3 By species

Code 
ggplot(est_bats, aes(x = sp, y = pc1)) + 
    geom_boxplot(fill = mako(10, alpha = 0.3)[4]) +
    labs(x = "Especie", y = "PC1") +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))

Code 
ggplot(est_bats, aes(x = sp, y = freq)) + 
    geom_boxplot(fill = mako(10, alpha = 0.3)[4]) +
    labs(x = "Especie", y = "Frecuencia pico (kHz)") +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))

Code 
ggplot(est_bats, aes(x = sp, y = duration)) + 
    geom_boxplot(fill = mako(10, alpha = 0.3)[4]) +
    labs(x = "Especie", y = "Duración (s)") +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))

Code 
ggplot(est_bats, aes(x = sp, y = rate)) + 
    geom_boxplot(fill = mako(10, alpha = 0.3)[4]) +
    labs(x = "Especie", y = "Tasa de llamados") +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
Warning: Removed 56 rows containing non-finite outside the scale range
(`stat_boxplot()`).

7.4 By species and habitat

Code 
ggplot(est_bats, aes(x = sp, y = pc1, fill = treatment)) + 
    scale_fill_viridis_d(option = "G", alpha = 0.7) +
    geom_boxplot() +
    labs(x = "Especie", y = "PC1 estructura de\nllamados", fill = "Tratamiento") +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))

Code 
ggplot(est_bats, aes(x = sp, y = freq, fill = treatment)) + 
    scale_fill_viridis_d(option = "G", alpha = 0.7) +
    geom_boxplot() +
    labs(x = "Especie", y = "Frecuencia (kHz)", fill = "Tratamiento")

Code 
ggplot(est_bats, aes(x = sp, y = duration, fill = treatment)) + 
    scale_fill_viridis_d(option = "G", alpha = 0.7) +
    geom_boxplot() +
    labs(x = "Especie", y = "Duración (s)", fill = "Tratamiento") +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))

Code 
ggplot(est_bats, aes(x = sp, y = rate, fill = treatment)) + 
    scale_fill_viridis_d(option = "G", alpha = 0.7) +
    geom_boxplot() +
    labs(x = "Especie", y = "Tasa de llamados", fill = "Tratamiento") +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
Warning: Removed 56 rows containing non-finite outside the scale range
(`stat_boxplot()`).

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] ggplot2_4.0.1      viridis_0.6.5      viridisLite_0.4.2  brmsish_1.0.0     
 [5] brms_2.23.0        Rcpp_1.1.0         readxl_1.4.3       ohun_1.0.3        
 [9] Rraven_1.0.14      warbleR_1.1.36     NatureSounds_1.0.5 seewave_2.2.3     
[13] tuneR_1.4.7        formatR_1.14       knitr_1.50        

loaded via a namespace (and not attached):
  [1] DBI_1.2.3             bitops_1.0-9          pbapply_1.7-4        
  [4] gridExtra_2.3         inline_0.3.21         remotes_2.5.0        
  [7] testthat_3.2.3        sandwich_3.1-0        rlang_1.1.6          
 [10] magrittr_2.0.4        multcomp_1.4-25       matrixStats_1.5.0    
 [13] e1071_1.7-16          compiler_4.3.2        loo_2.8.0            
 [16] reshape2_1.4.5        systemfonts_1.3.1     vctrs_0.6.5          
 [19] stringr_1.6.0         arrayhelpers_1.1-0    pkgconfig_2.0.3      
 [22] crayon_1.5.3          fastmap_1.2.0         backports_1.5.0      
 [25] labeling_0.4.3        rmarkdown_2.30        purrr_1.2.0          
 [28] xfun_0.55             jsonlite_2.0.0        tidybayes_3.0.7      
 [31] parallel_4.3.2        R6_2.6.1              StanHeaders_2.32.10  
 [34] stringi_1.8.7         RColorBrewer_1.1-3    brio_1.1.5           
 [37] cellranger_1.1.0      estimability_1.4.1    rstan_2.32.7         
 [40] zoo_1.8-12            bayesplot_1.15.0      Matrix_1.6-5         
 [43] splines_4.3.2         igraph_2.1.4          tidyselect_1.2.1     
 [46] rstudioapi_0.17.1     abind_1.4-8           yaml_2.3.12          
 [49] dtw_1.23-1            codetools_0.2-18      curl_7.0.0           
 [52] pkgbuild_1.4.8        plyr_1.8.9            lattice_0.20-45      
 [55] tibble_3.3.0          withr_3.0.2           bridgesampling_1.2-1 
 [58] S7_0.2.1              posterior_1.6.1       coda_0.19-4.1        
 [61] evaluate_1.0.5        signal_1.8-1          survival_3.2-13      
 [64] sf_1.0-22             sketchy_1.0.6         units_1.0-0          
 [67] proxy_0.4-27          RcppParallel_5.1.11-1 xml2_1.5.1           
 [70] ggdist_3.3.3          pillar_1.11.1         tensorA_0.36.2.1     
 [73] packrat_0.9.2         KernSmooth_2.23-20    stats4_4.3.2         
 [76] checkmate_2.3.3       distributional_0.5.0  generics_0.1.4       
 [79] RCurl_1.98-1.17       rstantools_2.5.0      scales_1.4.0         
 [82] xtable_1.8-4          class_7.3-20          glue_1.8.0           
 [85] emmeans_1.9.0         tools_4.3.2           xaringanExtra_0.8.0  
 [88] mvtnorm_1.3-3         cowplot_1.2.0         grid_4.3.2           
 [91] tidyr_1.3.1           ape_5.8-1             QuickJSR_1.8.1       
 [94] nlme_3.1-155          cli_3.6.5             kableExtra_1.4.0     
 [97] textshaping_1.0.4     svUnit_1.0.8          svglite_2.2.2        
[100] Brobdingnag_1.2-9     dplyr_1.1.4           V8_4.4.2             
[103] gtable_0.3.6          fftw_1.0-9            digest_0.6.39        
[106] classInt_0.4-11       TH.data_1.1-2         rjson_0.2.23         
[109] htmlwidgets_1.6.4     farver_2.1.2          htmltools_0.5.9      
[112] lifecycle_1.0.4       httr_1.4.7            MASS_7.3-55