Acoustic analysis

Echolocation plasticity

1 Purpose

• Qunatify effect of habitat vegetation density in bat echolocation call structure

 

2 Report overview

• You can have the sections listed here, for instance:

  • Lorem ipsum
  • Takeaways

 

3 Load packages

# 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", "ggtext"))

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'

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 Custom functions

get_stable_loadings <- function(pca,
                                pcs = 4,
                                B = 1000,
                                cum_threshold = 0.5,
                                freq_threshold = 0.75,
                                seed = 123) {
  
  set.seed(seed)
  
  # Reconstruct centered data
  X <- pca$x %*% t(pca$rotation)
  
  p <- ncol(pca$rotation)
  orig_rot <- pca$rotation[, 1:pcs]
  
  boot_loadings <- array(NA, dim = c(p, pcs, B))
  
  # -----------------------------
  # Bootstrap PCA
  # -----------------------------
  
  for (b in 1:B) {
    
    idx <- sample(1:nrow(X), replace = TRUE)
    Xb <- X[idx, ]
    
    pca_b <- prcomp(Xb, scale. = TRUE)
    rot_b <- pca_b$rotation[, 1:pcs]
    
    # Align signs
    for (k in 1:pcs) {
      if (cor(rot_b[, k], orig_rot[, k]) < 0) {
        rot_b[, k] <- -rot_b[, k]
      }
    }
    
    boot_loadings[, , b] <- rot_b
  }
  
  # -----------------------------
  # Summary statistics
  # -----------------------------

  abs_boot <- abs(boot_loadings)

mean_loading <- apply(abs_boot, c(1,2), mean)
ci_lower <- apply(abs_boot, c(1,2), quantile, 0.025)
ci_upper <- apply(abs_boot, c(1,2), quantile, 0.975)
  
  # -----------------------------
  # Stability frequency
  # -----------------------------
  
  top_freq <- matrix(0, nrow = p, ncol = pcs)
  
  for (b in 1:B) {
    for (k in 1:pcs) {
      
      sq <- boot_loadings[, k, b]^2
      ord <- order(sq, decreasing = TRUE)
      cumprop <- cumsum(sq[ord]) / sum(sq)
      
      selected <- ord[cumprop <= cum_threshold]
      selected <- c(selected, ord[min(which(cumprop >= cum_threshold))])
      
      top_freq[selected, k] <- top_freq[selected, k] + 1
    }
  }
  
  top_freq <- top_freq / B
  
  stable <- (top_freq >= freq_threshold) &
            (ci_lower > 0 | ci_upper < 0)
  
  # -----------------------------
  # Return tidy dataframe
  # -----------------------------
  
  data.frame(
    variable = rep(rownames(pca$rotation), pcs),
    ind = rep(colnames(pca$rotation)[1:pcs], each = p),
    mean_loading = as.vector(mean_loading),
    ci_lower = as.vector(ci_lower),
    ci_upper = as.vector(ci_upper),
    freq = as.vector(top_freq),
    stable = as.vector(stable)
  )
}

5 Clipping recordings

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

6 Importing annotations

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

pca <- prcomp(sp_feat[, -c(1:2)], scale. = TRUE)
saveRDS(pca, "./data/processed/pca_acoustic_features_2026.RDS")
est_bats$pc1 <- pca$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)

sp_feat2 <- sp_feat[!is.na(est_bats$sp), -c(1:3)]

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

est_bats <- cbind(est_bats, sp_feat2)

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

7 Statistical analysis

7.1 PCA

pca <- readRDS("./data/processed/pca_acoustic_features_2026.RDS")

# summary(pca)

# plot rotation values by PC
pca_rot <- as.data.frame(pca$rotation[, 1:4])
pca_var <- round(summary(pca)$importance[2, ] * 100)

pca_rot_stck <- stack(pca_rot)

pca_rot_stck$variable <- rownames(pca_rot)
pca_rot_stck$values[pca_rot_stck$ind == "PC1"] <- pca_rot_stck$values[pca_rot_stck$ind ==
    "PC1"]
pca_rot_stck$Sign <- ifelse(pca_rot_stck$values > 0, "Positive", "Negative")
pca_rot_stck$rotation <- abs(pca_rot_stck$values)
pca_rot_stck$ind_var <- paste0(pca_rot_stck$ind, " (", sapply(pca_rot_stck$ind,
    function(x) pca_var[names(pca_var) == x]), "%)")


pca_rot_stck$top_vars <- ave(
  abs(pca_rot_stck$values),
  pca_rot_stck$ind,
  FUN = function(x) {
    
    # Order decreasing
    ord <- order(x, decreasing = TRUE)
    x_sorted <- x[ord]
    
    # Cumulative proportion
    cumprop <- cumsum(x_sorted) / sum(x_sorted)
    
    selected_sorted <- cumprop <= 0.50 # variables with 50% of contribution 
    selected_sorted[which(cumprop >= 0.50)[1]] <- TRUE
    
    # Return logical vector in original order
    selected <- logical(length(x))
    selected[ord] <- selected_sorted
    
    selected
  }
)


# Create facet-specific variable
pca_rot_stck$var_facet <- paste(pca_rot_stck$ind, pca_rot_stck$variable, sep = "_")

# Reorder within each ind by rotation (largest at top after coord_flip)
pca_rot_stck <- do.call(rbind, lapply(split(pca_rot_stck, pca_rot_stck$ind), function(df) {
  
  df$var_facet <- factor(
    df$var_facet,
    levels = df$var_facet[order(df$rotation)]
  )
  
  df
}))



# add which variables are stable
stable_df <- get_stable_loadings(pca, pcs = 4, B = 1000, cum_threshold = 0.5, freq_threshold = 0.5, seed = 123)

pca_rot_stck <- merge(
  pca_rot_stck,
  stable_df,
  by = c("variable", "ind"),
  all.x = TRUE
)

pca_rot_stck$top_vars <- ifelse(pca_rot_stck$stable, 0.6, 1)


# Build colored labels per row

# Colored labels
pca_rot_stck$label_col <- ifelse(
  pca_rot_stck$stable,
  paste0("<span style='color:black;'>", pca_rot_stck$variable, "</span>"),
  paste0("<span style='color:gray50;'>", pca_rot_stck$variable, "</span>")
)

# Named vector for labels
label_vec <- setNames(pca_rot_stck$label_col, pca_rot_stck$var_facet)

# absolute CI
pca_rot_stck$ci_low_plot  <- abs(pca_rot_stck$ci_lower)
pca_rot_stck$ci_high_plot <- abs(pca_rot_stck$ci_upper)


# Plot
ggplot(pca_rot_stck, aes(x = var_facet, y = rotation, fill = Sign, alpha = as.factor(top_vars))) +
  geom_col() +
      geom_errorbar(aes(ymin = ci_low_plot,
                    ymax = ci_high_plot),
                width = 0.1,
                linewidth = 0.2) +
  coord_flip() +
  scale_alpha_manual(values = pca_rot_stck$top_vars, guide = NULL) +
  scale_x_discrete(labels = label_vec, name = "Variable") +
  scale_fill_viridis_d(alpha = 0.7, begin = 0.2, end = 0.8) +
  facet_wrap(~ind_var, scales = "free_y") +
  theme_classic() +
  theme(
    axis.text.y = element_markdown()
  )

Bar plots show variable loadings for the first four principal components. Loadings are computed from a PCA performed on standardized acoustic variables (mean-centered and scaled to unit variance). Bars represent the mean absolute loading estimated from 1,000 bootstrap resamples of the original dataset. For each of 1,000 bootstrap iterations, we resampled observations (rows of the original data matrix) with replacement to generate a dataset of equal size to the original. Horizontal whiskers indicate 95% bootstrap confidence intervals (2.5–97.5 percentiles of the bootstrap distribution). Variables are ordered within each component by their contribution (squared loading). Bar color indicates the sign of the original loading (green = positive, purple = negative). Transparency highlights whether a variable consistently contributes to the component: variables are considered stable when (i) their bootstrap confidence interval does not overlap zero and (ii) they are selected among the set of variables jointly explaining at least 50% of the component variance in ≥50% of bootstrap replicates. This procedure allows assessment of both the magnitude and robustness of variable contributions to each principal component.

7.2 Predictions

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

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

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

7.3 Macroevolution

7.3.1 By guild

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 * 2,
        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:

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	473 (0.19%)	0	1220.5	988.289	177549093

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	1.006	-3.183	5.310	1.002	1220.500	988.289
b_guildnarrow	0.020	-1.882	1.864	1	1664.249	1170.528

Frequency:

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, 15.9) sd-exponential(1) sd-exponential(1) sigma-exponential(1)	5000	2	1	2500	1 (2e-04%)	0	4591.644	3685.438	2074619662

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	7.864	-12.422	28.547	1	4591.644	3685.438
b_guildnarrow	-0.107	-2.085	1.834	1	7818.703	3761.970

Rate:

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)	10000	1	1	5000	1123 (0.22%)	0	1261.86	937.138	2095919153

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	0.069	-0.052	0.190	1.001	1307.401	1288.319
b_guildnarrow	-0.001	-0.126	0.133	1	1261.860	937.138

PC1:

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.6) sd-exponential(2) sd-exponential(2) sigma-exponential(2)	5000	2	1	2500	4 (8e-04%)	0	3663.832	3436.589	539921948

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-0.276	-2.342	1.861	1.001	3663.832	3436.589
b_guildnarrow	0.249	-1.589	2.026	1	5435.457	3646.392

7.3.2 By species

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"),

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

  # Residual standard deviation
  set_prior("exponential(1)", class = "sigma")
)

freq_mod <- brm(
        formula = freq ~ sp + (1 |
                                          indiv) +  (1 | treatment) + ar(p = 2, time = orig.start, gr = orig.sf),
        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"),

  # 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 | treatment) + ar(p = 2, time = orig.start, gr = orig.sf),
        family = Gamma(),
        data = est_bats,
        prior = priors,
        chains = chains,
        iter = iter * 2,
        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"),

  # 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 | treatment) + ar(p = 2, time = orig.start, gr = orig.sf),
        data = est_bats,
        prior = priors,
        chains = chains,
        iter = iter,
        file = "./data/processed/pc1_sp_macro_model.RDS"
    )

Duration:

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 | treatment)	gamma (inverse)	b-normal(0, 1) Intercept-normal(1, 2) sd-student_t(3, 0, 327.1) sd-exponential(1) shape-exponential(1)	5000	2	1	2500	6 (0.0012%)	0	7614.584	3361.836	467736439

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	0.981	-3.017	4.972	1.001	8218.990	3661.641
b_spC.castanea	0.078	-1.845	1.946	1.001	8358.965	3833.927
b_spC.sowelli	0.019	-2.047	2.049	1	7614.584	3361.836
b_spM.nigricans	-0.008	-1.914	1.923	1	8836.799	3653.982
b_spP.discolor	0.012	-1.987	1.936	1	10002.967	3620.365

Frequency:

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 | treatment) + ar(p = 2, time = orig.start, gr = orig.sf)	gaussian (identity)	b-normal(0, 1) Intercept-normal(0, 10) sd-student_t(3, 0, 15.9) sd-exponential(1) sigma-exponential(1)	5000	2	1	2500	10 (0.002%)	0	4673.039	3051.127	646708152

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	8.237	-12.012	28.751	1.001	4673.039	3473.930
b_spC.castanea	0.316	-1.640	2.273	1	5968.647	3051.127
b_spC.sowelli	-0.035	-1.917	1.914	1	5868.048	3639.344
b_spM.nigricans	0.070	-1.893	2.050	1	6156.631	3535.837
b_spP.discolor	-0.622	-2.490	1.275	1.001	5889.512	3550.230

Rate:

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 | treatment) + 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) sderr-student_t(3, 0, 2.5) shape-exponential(1)	10000	1	1	5000	148 (0.03%)	0	1088.599	686.747	696759494

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	0.060	0.005	0.114	1	1088.599	686.747
b_spC.castanea	0.052	0.000	0.107	1.002	3065.442	2533.847
b_spC.sowelli	0.007	-0.037	0.056	1.001	2697.004	2572.708
b_spM.nigricans	0.008	-0.052	0.069	1.001	3360.607	2937.087
b_spP.discolor	-0.021	-0.064	0.025	1.004	2713.169	2223.111

PC1:

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 | treatment) + 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.6) sd-exponential(2) sigma-exponential(2)	5000	2	1	2500	21 (0.0042%)	0	1199.795	708.66	432212221

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-0.207	-1.622	1.127	1.002	1319.042	708.660
b_spC.castanea	-1.996	-3.200	-0.627	1.006	1712.827	1849.149
b_spC.sowelli	-0.434	-1.565	0.779	1.002	3350.068	2630.150
b_spM.nigricans	-0.694	-2.104	0.814	1.002	2675.381	2592.366
b_spP.discolor	3.235	1.592	4.501	1.005	1199.795	1515.622

PC1 contrasts:

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


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

7.4 Plasticity

7.4.1 Only habitat

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),
        family = Gamma(),
        data = est_bats,
        prior = priors,
        chains = chains,
        iter = iter * 2,
        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:

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)	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)	10000	2	1	5000	0 (0%)	0	14329.29	6800.629	964104920

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	1.693	-2.287	5.612	1.001	16258.15	7335.222
b_treatmentafuera	-0.016	-1.984	1.903	1	14956.31	6800.629
b_treatmentcerrado	0.039	-1.896	2.024	1	14329.29	7078.364

Frequency:

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, 15.9) sd-exponential(1) sd-exponential(1) sigma-exponential(1)	5000	2	1	2500	0 (0%)	0	1266.387	2306.524	2034681173

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	57.413	37.429	72.090	1	1266.387	2306.524
b_treatmentafuera	-0.638	-2.570	1.283	1	6118.916	3561.138
b_treatmentcerrado	0.064	-1.787	1.911	1.001	6273.085	3613.961

Rate:

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	1275 (0.26%)	0	1515.517	1205.672	300310293

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	0.071	0.033	0.109	1.001	1515.517	1205.672
b_treatmentafuera	0.001	-0.022	0.027	1	3013.530	2772.100
b_treatmentcerrado	-0.007	-0.023	0.008	1	3341.345	2445.124

PC1:

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.6) sd-exponential(2) sd-exponential(2) sigma-exponential(2)	5000	2	1	2500	0 (0%)	0	2294.978	3156.907	936212012

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	0.007	-1.423	1.408	1.001	2294.978	3156.907
b_treatmentafuera	0.492	-0.644	1.582	1	4738.195	3698.792
b_treatmentcerrado	-0.662	-1.503	0.118	1.001	5191.332	3588.190

PC1 contrasts:

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
)

7.4.2 Habitat by guild

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:

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:

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:

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:

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)

8 Plots

8.1 By guild

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

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

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

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 58 rows containing non-finite outside the scale range
(`stat_boxplot()`).

8.2 By treatment

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

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

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

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 58 rows containing non-finite outside the scale range
(`stat_boxplot()`).

8.3 By species

Predicted variables

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

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

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

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 58 rows containing non-finite outside the scale range
(`stat_boxplot()`).

Additional variables

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

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

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

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

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

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

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

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

8.4 By species and habitat

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

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

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

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 58 rows containing non-finite outside the scale range
(`stat_boxplot()`).

Takeaways

 

---

 

Session information

Warning in system("timedatectl", intern = TRUE): running command 'timedatectl'
had status 1

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] ggtext_0.1.2       ggplot2_4.0.2      viridis_0.6.5      viridisLite_0.4.3 
 [5] brmsish_1.0.0      brms_2.23.0        Rcpp_1.1.1         readxl_1.4.5      
 [9] ohun_1.0.5         Rraven_1.0.16      warbleR_1.1.37     NatureSounds_1.0.5
[13] seewave_2.2.4      tuneR_1.4.7        formatR_1.14       knitr_1.51        

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-1        rlang_1.1.7          
 [10] magrittr_2.0.4        multcomp_1.4-28       matrixStats_1.5.0    
 [13] e1071_1.7-17          compiler_4.5.2        loo_2.9.0            
 [16] reshape2_1.4.5        systemfonts_1.3.1     vctrs_0.7.1          
 [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        markdown_2.0         
 [28] purrr_1.2.0           xfun_0.56             litedown_0.7         
 [31] jsonlite_2.0.0        tidybayes_3.0.7       parallel_4.5.2       
 [34] R6_2.6.1              StanHeaders_2.32.10   stringi_1.8.7        
 [37] RColorBrewer_1.1-3    brio_1.1.5            cellranger_1.1.0     
 [40] estimability_1.5.1    rstan_2.32.7          zoo_1.8-14           
 [43] bayesplot_1.15.0      Matrix_1.7-4          splines_4.5.2        
 [46] igraph_2.2.1          tidyselect_1.2.1      rstudioapi_0.17.1    
 [49] abind_1.4-8           yaml_2.3.12           dtw_1.23-1           
 [52] codetools_0.2-20      curl_7.0.0            pkgbuild_1.4.8       
 [55] plyr_1.8.9            lattice_0.22-7        tibble_3.3.0         
 [58] withr_3.0.2           bridgesampling_1.2-1  S7_0.2.1             
 [61] posterior_1.6.1       coda_0.19-4.1         evaluate_1.0.5       
 [64] signal_1.8-1          survival_3.8-3        sf_1.0-24            
 [67] sketchy_1.0.7         units_1.0-0           proxy_0.4-29         
 [70] RcppParallel_5.1.11-1 xml2_1.5.2            ggdist_3.3.3         
 [73] pillar_1.11.1         tensorA_0.36.2.1      packrat_0.9.3        
 [76] KernSmooth_2.23-26    stats4_4.5.2          checkmate_2.3.4      
 [79] distributional_0.5.0  generics_0.1.4        RCurl_1.98-1.17      
 [82] commonmark_1.9.5      rstantools_2.5.0      scales_1.4.0         
 [85] xtable_1.8-4          class_7.3-23          glue_1.8.0           
 [88] emmeans_1.11.1        tools_4.5.2           xaringanExtra_0.8.0  
 [91] mvtnorm_1.3-3         cowplot_1.2.0         grid_4.5.2           
 [94] tidyr_1.3.2           ape_5.8-1             QuickJSR_1.8.1       
 [97] nlme_3.1-168          cli_3.6.5             kableExtra_1.4.0     
[100] textshaping_1.0.4     svUnit_1.0.8          svglite_2.2.2        
[103] Brobdingnag_1.2-9     dplyr_1.1.4           V8_6.0.4             
[106] gtable_0.3.6          fftw_1.0-9            digest_0.6.39        
[109] classInt_0.4-11       TH.data_1.1-3         rjson_0.2.23         
[112] htmlwidgets_1.6.4     farver_2.1.2          htmltools_0.5.9      
[115] lifecycle_1.0.5       httr_1.4.7            gridtext_0.1.5       
[118] MASS_7.3-65