Song geographic variation of Bearded Bellbirds

Obtaining nature media with the R package suwo

Author

Published

July 30, 2026

Source code and data found at https://github.com/maRce10/suwo-paper

Code

# options to customize chunk outputs
knitr::opts_chunk$set(
  tidy.opts = list(width.cutoff = 65), 
  tidy = TRUE,
  message = FALSE
 )

This analysis serves as the animal-behavior-research case study accompanying the manuscript “Obtaining nature media with the R package suwo”, illustrating how suwo-retrieved media can be carried into a full downstream analytical pipeline

Purpose

Obtain and curate recordings of Bearded Bellbird (Procnias averano) and their associated metadata from multiple online repositories, using the R package suwo.
Quantify acoustic structure in Bearded Bellbird “bock” songs (Snow 1970).
Infer geographic variation in vocalizations, using recordings assembled with the R package suwo and its covariation with time, geographic distance and population structure

Analysis flowchart

flowchart LR
  A[Download data] --> B(Format data) 
  B --> C(Graphs)
  C --> D('MFCCs<br/>PCA')
  D --> E(Statistical model)

style A fill:#382A5433
style B fill:#395D9C33
style C fill:#3497A933
style D fill:#60CEAC33
style E fill:#60CEAC33

Load packages

Code

# if sketchy not installed install it
if (!require("sketchy")) {
    install.packages("sketchy")
}

# install/ load packages

sketchy::load_packages(packages = c("knitr", "warbleR", "Rraven",
    "suwo", "ggplot2", "viridis", "ggtext", "brms", "geosphere", "brmsish",
    "ggcorrplot", "ggdist", "posterior", "bayestestR"))

1 Custom functions

broken_stick_k() implements the PCA retention rule used in the pipeline below.

Code

## ============================================================
## PC RETENTION: BROKEN-STICK MODEL (Frontier 1976; MacArthur
## 1957)
## ============================================================

broken_stick_k <- function(X) {

    X <- as.matrix(X)
    ev <- prcomp(X, scale. = TRUE)$sdev^2
    p <- length(ev)
    prop <- ev/sum(ev)

    # expected proportion for component k under random division
    # of variance
    bs <- vapply(seq_len(p), function(k) sum(1/(k:p))/p, numeric(1))

    above <- prop > bs
    first_fail <- which(!above)[1]
    k <- if (is.na(first_fail))
        p else first_fail - 1L

    list(k = max(1L, k), proportion = prop, expected = bs, cum_var = cumsum(prop))
}

dir_word <- function(row) if (row$Estimate > 0) "increased" else "decreased"
cred_word <- function(is_credible) if (is_credible) "a credible" else "no credible"

2 Obtain recordings

The following code queries and downloads all available Bearded Bellbird recordings from 5 repositories:

Code

# search species

# time process
init_time <- Sys.time()
options(suwo_species = "Procnias averano", suwo_format = "sound")
mac <- query_macaulay(path = tempdir())
xc <- query_xenocanto()
wa <- query_wikiaves()
inat <- query_inaturalist()
gbif <- query_gbif()
end_time <- Sys.time()

dura <- end_time - init_time


# merge data
dat <- merge_metadata(wa, inat, gbif, xc, mac)

# find duplicates and exclude them
dat <- find_duplicates(dat)

dedup_dat <- remove_duplicates(dat)


# download all files first
dm <- download_media(dedup_dat, path = "./data/raw/p_averano_recordings",
    cores = 10)

# save metadata
write.csv(dm, "./data/processed/p_averano_recordings_metadata.csv",
    row.names = FALSE)

# now look again for duplicates using file size and user name as
# criteria, and remove them
dm <- find_duplicates(metadata = dm, criteria = paste("user_name > 0.7",
    "file_size > 0.99", sep = " & "))

dedup_dm <- remove_duplicates(dm)

dm <- download_media(dedup_dat, path = "./data/raw/p_averano_recordings",
    cores = 10)

dedup_dm <- dm[!duplicated(dm$duplicate_group) | is.na(dm$duplicate_group),
    ]

duplicates <- list.files("./data/raw/p_averano_recordings")[!list.files("./data/raw/p_averano_recordings") %in%
    dedup_dm$downloaded_file_name]

unlink(file.path("./data/raw/p_averano_recordings", duplicates))

dedup_dm <- download_media(dedup_dm, path = "./data/raw/p_averano_recordings",
    cores = 23)


write.csv(dedup_dm, "./data/processed/p_averano_recordings_metadata.csv",
    row.names = FALSE)

3 Curate annotations

Code

txts <- list.files("/home/m/Insync/marceloa27@gmail.com/Google Drive/grabaciones_para_ordenar/Anotaciones/Selecciones/",
    full.names = TRUE)

# file.copy(from = txts, to =
# './data/raw/p_averano_annotations/', overwrite = TRUE)

anns <- imp_raven(path = "./data/raw/p_averano_annotations/", warbler.format = TRUE,
    name.from.file = TRUE, ext.case = "lower", all.data = TRUE)
ncol(anns)

anns$sound.files <- NULL


sfs <- gsub(".Band.Limited.Energy.Detector.wav|.Table1.wav", ".wav",
    anns$sound.files)
anns$sound.files <- sfs

# head(anns)

anns <- check_sels(anns, path = "./data/raw/p_averano_recordings",
    fix.selec = TRUE)

Code

metadata <- read.csv("./data/processed/p_averano_recordings_metadata.csv", stringsAsFactors = FALSE)

metadata <- metadata[!is.na(metadata$downloaded_file_name), ]

# convert all to .wav
metadata$downloaded_file_name <- gsub(".mp3$|.m4a$", ".wav", metadata$downloaded_file_name)


# all(unique(anns$sound.files) %in% metadata$downloaded_file_name)


metadata <- metadata[metadata$downloaded_file_name %in% unique(anns$sound.files), ]

metadata$type <- sapply(metadata$downloaded_file_name, function(x){
    
    out <- paste(unique(na.omit(anns$Tipo[anns$sound.files == x])), collapse = "-")
    
    if (length(out) == 0) {
      out <- NA
    }
    
    return(out)
})


# add lat lon
localidades <- data.frame(
  archivo = c(
    "Cumaca Valley",
    "Jaqueira/PE",
    "Guaramiranga/CE",
    "Pacatuba/CE",
    "Maranguape/CE",
    "São Luís/MA",
    "Caxias/MA",
    "FORA DO BRASIL/EX",
    "Dianópolis/TO",
    "Araguaína/TO",
    "Fazenda Nazaret, estado do Piauí"
  ),
  Latitude = c(
    10.7000,     # Cumaca Valley, Trinidad y Tobago
    -8.9967,     # Jaqueira/PE
    -4.2725,     # Guaramiranga/CE
    -3.9822,     # Pacatuba/CE
    -3.8908,     # Maranguape/CE
    -2.5297,     # São Luís/MA
    -4.8592,     # Caxias/MA
    NA,          # FORA DO BRASIL/EX - no es una localidad real
    -11.6236,    # Dianópolis/TO
    -7.1911,     # Araguaína/TO
    -9.3353      # Fazenda Nazaré, PI (confianza media)
  ),
  Longitude = c(
    -61.1500,
    -35.7014,
    -38.9328,
    -38.6206,
    -38.6844,
    -44.3028,
    -43.3567,
    NA,
    -46.8228,
    -48.2072,
    -45.5725
  )
)

for (i in which(is.na(metadata$latitude))){
    
metadata$latitude[i] <- localidades$Latitude[localidades$archivo == metadata$locality[i]]

metadata$longitude[i] <- localidades$Longitude[localidades$archivo == metadata$locality[i]]    

}


metadata$population <- "non_class"

# brazil population
metadata$population <- ifelse(grepl("razil", metadata$country) & is.na(metadata$longitude), "Brazil", metadata$population)
metadata$population <- ifelse(metadata$latitude <= -1.407 & metadata$longitude < -2.63, "Brazil", metadata$population)


# trinidad and tobago
metadata$population <- ifelse(grepl("rinidad", metadata$country), "Trinidad and Tobago", metadata$population)
metadata$population <- ifelse(metadata$latitude > 9.67 & metadata$longitude > -62.01 & metadata$population == "non_class", "Trinidad and Tobago", metadata$population)

# venezuela SE
metadata$population <- ifelse(grepl("non_class", metadata$population) & metadata$latitude <  7.56, "Venezuela SE", metadata$population)

# venezuela SE
metadata$population <- ifelse(grepl("non_class", metadata$population) & metadata$latitude >  7.56, "Venezuela N", metadata$population)

# refill brazil
metadata$population <- ifelse(is.na(metadata$population) & metadata$country ==  "Brazil", "Brazil", metadata$population)




anns$population <- sapply(anns$sound.files, function(x){
  pop <- metadata$population[metadata$downloaded_file_name == x]
  if (length(pop) == 0) {
    pop <- NA
  }
  return(pop)
})

table(anns$population, anns$Tipo)

anns$latitude <- sapply(anns$sound.files, function(x){
  lat <- metadata$latitude[metadata$downloaded_file_name == x]
  if (length(lat) == 0) {
    lat <- NA
  }
  return(lat)
})

anns$longitude <- sapply(anns$sound.files, function(x){
  lon <- metadata$longitude[metadata$downloaded_file_name == x]
  if (length(lon) == 0) {
    lon <- NA
  }
  return(lon)
})

anns$locality <- sapply(anns$sound.files, function(x){
  loc <- metadata$locality[metadata$downloaded_file_name == x]
  if (length(loc) == 0) {
    loc <- NA
  }
  return(loc)
})

anns$year <- sapply(anns$sound.files, function(x){
  yr <- metadata$date[metadata$downloaded_file_name == x]
  yr <- substring(yr, 1, 4)
  
  if (length(yr) == 0) {
    yr <- NA
  }
  return(yr)
})

unique(anns$locality[is.na(anns$latitude)])

est <- selection_table(anns, path = "./data/raw/p_averano_recordings", extended = TRUE)

est$year <-  anns$year

attr(est, "metadata") <- metadata

saveRDS(est, "./data/processed/p_averano_extended_selection_table.rds")

4 Map locations

Code

est <- readRDS("./data/processed/p_averano_extended_selection_table.rds")

est$orig.sound.files <- attr(est, "check.res")$orig.sound.files

metadata_est <- attr(est, "metadata")

map_locations(metadata_est, tags = c("latitude", "longitude", "population"),
    by = "population", palette = function(n) viridis::mako(n, begin = 0.2,
        end = 0.75))

5 Principal component analysis

Acoustic structure for “bock” songs was summarized with a set of MFCC-derived statistics (mfcc_stats(), warbleR) computed over the frequency band spanned by the annotated selections, then reduced with a principal component analysis (PCA) on the standardized statistics. Pairwise Euclidean distance between individuals in this reduced space is the acoustic-distance response used in every model below, so the number of components retained here directly determines what the response variable is, not just how it is summarized.

Why PCA, and why the number of retained components matters more than it might seem. Because PCA is an orthonormal rotation, Euclidean distance computed on all components is mathematically identical to Euclidean distance on the standardized input features — rotating axes does not change distances between points. Truncating to fewer components is therefore not a data-reduction convenience but a deliberate low-rank denoising step (Eckart and Young 1936): components beyond some point are assumed to carry more measurement noise than acoustic signal, and dropping them removes that noise from the response. The question this raises is not “what fraction of variance should I keep” (there is no principled answer to that) but “at which component does structure give way to noise.”

Why the broken-stick rule, specifically. A fixed cumulative-variance threshold (70%, 95%, …) is a common convention but has no sampling theory behind it and is explicitly discouraged as a stopping rule (Cangelosi and Goriely 2007) — the threshold itself is arbitrary. The eigenvalue-greater-than-one (“Kaiser”) criterion is the rule most frequently reported in avian bioacoustic studies, but it is documented to over-retain components substantially, and the number it retains scales with the number of input variables rather than with genuine structure in the data (Zwick and Velicer 1986; Jackson 1993) — a real liability here, since the MFCC statistic matrix has many columns. The broken-stick model instead compares each component’s observed proportion of variance against the proportion it would be expected to carry if the total variance were divided at random among all components (as if breaking a stick into random pieces); a component is retained only if it exceeds that null expectation, and retention proceeds sequentially until the first component that fails to clear its own threshold. This was the best-performing stopping rule in Jackson’s (1993, Ecology 74:2204–2214) simulation comparison across ecological datasets, and it is also the most conservative (parsimonious) of the widely used criteria — the right direction to err in here, specifically because the acoustic dissimilarity is the response variable rather than a predictor: every additional noise component retained adds variance to the response, which widens posterior intervals and reduces power to detect genuine geographic or temporal structure downstream. Retaining too few components risks discarding real signal; retaining too many risks diluting the response with noise. The broken-stick rule was chosen because, per the simulation evidence above, it is the criterion best calibrated to that specific trade-off, not because it is the most common choice in the literature.

Code

pc_rule <- "broken_stick"
fit_suffix <- if (identical(pc_rule, "broken_stick")) "" else paste0("_",
    pc_rule)
pca_file <- paste0("./data/processed/pca_song_type_1", fit_suffix,
    ".rds")

est <- readRDS("./data/processed/p_averano_extended_selection_table.rds")
est$orig.sound.files <- attr(est, "check.res")$orig.sound.files

est_tp <- est[est$Tipo == "1", ]

## ---- PCA on standardized MFCC statistics
## --------------------------------

mfcc_tp <- mfcc_stats(est_tp, bp = c(min(est_tp$bottom.freq), max(est_tp$top.freq)))
mfcc_feat <- mfcc_tp[, c(-1, -2)]

pca_tp <- prcomp(mfcc_feat, scale. = TRUE)
cum_var <- summary(pca_tp)$importance[3, ]

## PC retention via the broken-stick rule (see Custom functions
## above for the rationale). pc_rule only changes the
## cached-file path (fit_suffix) at this point, not which rule
## is applied -- broken_stick_k() runs unconditionally.
n_pcs <- broken_stick_k(mfcc_feat)$k

pcs <- paste0("PC", seq_len(n_pcs))
pc_scores <- as.data.frame(pca_tp$x[, seq_len(n_pcs), drop = FALSE])
colnames(pc_scores) <- pcs

## fix_extended_selection_table() re-attaches proper
## extended_selection_table attributes/wave-object references
## after cbind() adds plain columns -- needed so downstream
## functions that expect a genuine extended_selection_ table
## (e.g. sig2noise(), single-bracket subsetting) keep working
## correctly.
est_tp2 <- est_tp
est_tp <- cbind(est_tp, pc_scores)
est_tp <- fix_extended_selection_table(X = est_tp, Y = est_tp2)

## saved so the PCA never needs to be recomputed downstream: the
## diagnostic section below and the Model fitting section both
## read this file rather than re-running mfcc_stats()/prcomp() a
## second and third time (which previously risked the two
## silently drifting apart).
saveRDS(list(est_tp = est_tp, pcs = pcs, n_pcs = n_pcs, cum_var = cum_var),
    pca_file)

Code

pc_rule <- "broken_stick"
fit_suffix <- if (identical(pc_rule, "broken_stick")) "" else paste0("_",
    pc_rule)
pca_file <- paste0("./data/processed/pca_song_type_1", fit_suffix,
    ".rds")

pca_out <- readRDS(pca_file)

The first 14 principal components were retained, which together explain 58.4% of the variance in the standardized MFCC statistics. The retained components are used to compute pairwise Euclidean distances between individuals, which is the acoustic-distance response variable used in all models below.`

6 Check distances

Before fitting the model, the pairwise predictors are inspected directly: the shape of the geographic-distance distribution (including how the zero-distance pairs created by shared-locality coordinates are handled), and the correlation structure among predictors — particularly between same_population and geo_distance_sc, and between time_separation and geographic distance. These checks motivate two choices used throughout the modeling pipeline below: the log-transform and additive constant applied to geographic distance, and the decision to include same_population and geo_distance_sc together in the same model despite their moderate correlation.

Code

pc_rule <- "broken_stick"
fit_suffix <- if (identical(pc_rule, "broken_stick")) "" else paste0("_",
    pc_rule)

pca_file <- paste0("./data/processed/pca_song_type_1", fit_suffix,
    ".rds")
res_file <- paste0("./data/processed/results_song_type_1", fit_suffix,
    ".rds")
dist_long_file <- paste0("./data/processed/dist_long_song_type_1",
    fit_suffix, ".rds")

pca_out <- readRDS(pca_file)
est_tp <- pca_out$est_tp
pcs <- pca_out$pcs
n_pcs <- pca_out$n_pcs
cum_var <- pca_out$cum_var

## PCA recap -- printed here, immediately before any regression
## model is fitted (full detail is in the 'Principal component
## analysis' section).
cat("==== PCA recap ====\n")

==== PCA recap ====

Code

cat("PCs retained:", n_pcs, "| cumulative variance explained:", round(cum_var[n_pcs] *
    100, 1), "%\n\n")

PCs retained: 14 | cumulative variance explained: 58.4 %

Code

## ---- 1b. Signal-to-noise ratio per selection
## --------------------------- Recording quality, used below as
## a nuisance covariate (see Statistical analysis text). mar
## sets the margin adjacent to each selection over which
## background noise is sampled and must not overlap a
## neighbouring signal -- CHECK against this dataset with
## snr_spectrograms() before trusting the value below, which is
## a starting point (warbleR tutorial default), not a tuned
## setting. bp matches mfcc_stats() above so SNR is measured
## over the same frequency band as the acoustic features it is
## meant to help account for.

snr_tp <- sig2noise(est_tp, mar = 0.05, bp = c(min(est_tp$bottom.freq),
    max(est_tp$top.freq)))
est_tp$SNR <- snr_tp$SNR

## ---- 2. Aggregate by individual (mean PCs / SNR per
## orig.sound.files) --

agg_pcs <- aggregate(est_tp[, pcs, drop = FALSE], by = list(orig.sound.files = est_tp$orig.sound.files),
    FUN = mean)

agg_snr <- aggregate(est_tp$SNR, by = list(orig.sound.files = est_tp$orig.sound.files),
    FUN = mean)
names(agg_snr)[ncol(agg_snr)] <- "SNR"

agg_meta <- est_tp[!duplicated(est_tp$orig.sound.files), c("orig.sound.files",
    "population", "latitude", "longitude", "year")]

agg_indiv <- merge(agg_pcs, agg_meta, by = "orig.sound.files")
agg_indiv <- merge(agg_indiv, agg_snr, by = "orig.sound.files")
rownames(agg_indiv) <- agg_indiv$orig.sound.files

## year is built with substring() upstream and is therefore
## character; dist() needs it numeric.
agg_indiv$year <- as.numeric(as.character(agg_indiv$year))

## ---- 3. Pairwise distance matrices
## -------------------------------------

dist_acoustic_mat <- as.matrix(dist(agg_indiv[, pcs, drop = FALSE],
    method = "euclidean"))

coords <- as.matrix(agg_indiv[, c("longitude", "latitude")])
dist_geo_mat <- geosphere::distm(coords, fun = distHaversine)/1000

dist_time_mat <- as.matrix(dist(agg_indiv[, "year", drop = FALSE],
    method = "euclidean"))

## pairwise minimum SNR: a pair's acoustic dissimilarity is
## corrupted by whichever of its two recordings is noisier, not
## by their average, so the predictor takes the pair's minimum
## rather than the pair's mean.
snr_pair_mat <- outer(agg_indiv$SNR, agg_indiv$SNR, pmin)

rownames(dist_acoustic_mat) <- colnames(dist_acoustic_mat) <- agg_indiv$orig.sound.files
rownames(dist_geo_mat) <- colnames(dist_geo_mat) <- agg_indiv$orig.sound.files
rownames(dist_time_mat) <- colnames(dist_time_mat) <- agg_indiv$orig.sound.files
rownames(snr_pair_mat) <- colnames(snr_pair_mat) <- agg_indiv$orig.sound.files

## ---- 4. Long format (upper triangle)
## -----------------------------------

idx <- which(upper.tri(dist_acoustic_mat), arr.ind = TRUE)

dist_long <- data.frame(individual1 = rownames(dist_acoustic_mat)[idx[,
    1]], individual2 = colnames(dist_acoustic_mat)[idx[, 2]], acoustic_dissimilarity = dist_acoustic_mat[idx],
    geo_distance = dist_geo_mat[idx], time_separation = dist_time_mat[idx],
    min_snr = snr_pair_mat[idx])

pop_lookup <- setNames(agg_indiv$population, agg_indiv$orig.sound.files)
dist_long$population1 <- pop_lookup[dist_long$individual1]
dist_long$population2 <- pop_lookup[dist_long$individual2]

dist_long$same_population <- as.factor(as.integer(dist_long$population1 ==
    dist_long$population2))

dist_long$individual1 <- factor(dist_long$individual1)
dist_long$individual2 <- factor(dist_long$individual2)
dist_long$population1 <- factor(dist_long$population1)
dist_long$population2 <- factor(dist_long$population2)

## ---- 5. Transformations and scaling
## -------------------------------------

dist_long$acoustic_dissimilarity_sc <- scale(dist_long$acoustic_dissimilarity)[,
    1]

geo_const <- unname(quantile(dist_long$geo_distance, 0.25, na.rm = TRUE))

dist_long$log_geo_distance <- log(dist_long$geo_distance + geo_const)

## Divide by 2 SD, not 1 SD (Gelman 2008): puts this continuous
## predictor's coefficient on the same 'typical full swing'
## scale as the untouched binary same_population (0/1)
## predictor, so binary-vs- continuous effect-size comparisons
## in the results text are valid.
dist_long$geo_distance_sc <- scale(dist_long$log_geo_distance)[, 1]/2

cat("geo_const (km):", signif(geo_const, 3), "| zero-distance pairs:",
    sum(dist_long$geo_distance == 0), "(", round(100 * mean(dist_long$geo_distance ==
        0), 1), "% )\n")

geo_const (km): 529 | zero-distance pairs: NA ( NA % )

Code

## same 2-SD convention applied consistently to every continuous
## predictor
dist_long$time_separation_sc <- scale(dist_long$time_separation)[,
    1]/2

## 2-SD scaled, not log-transformed: unlike geographic distance,
## min_snr (dB) has no comparable right-skew / zero-mass issue
## here -- revisit if the observed distribution of min_snr turns
## out to be skewed once run.
dist_long$min_snr_sc <- scale(dist_long$min_snr)[, 1]/2

## ---- 5b. Filter to complete cases
## --------------------------------------- Every variable used
## in the model is filtered together here, once, so that the
## fitted model and every downstream check use exactly the same
## set of pairs.

model_vars <- c("acoustic_dissimilarity_sc", "same_population", "geo_distance_sc",
    "time_separation_sc", "min_snr_sc", "population1", "population2",
    "individual1", "individual2")

n_before <- nrow(dist_long)
dist_long <- dist_long[complete.cases(dist_long[, model_vars]), ]
n_after <- nrow(dist_long)

if (n_after < n_before) {
    cat("rows dropped due to NAs:", n_before - n_after, "of", n_before,
        "(", round(100 * (n_before - n_after)/n_before, 1), "% )\n")
}

rows dropped due to NAs: 466 of 7021 ( 6.6 % )

Code

## mm(individual1, individual2) and mm(population1, population2)
## require both columns of each pair to share exactly the same
## set of levels.

individual_levels <- union(as.character(dist_long$individual1), as.character(dist_long$individual2))
dist_long$individual1 <- factor(as.character(dist_long$individual1),
    levels = individual_levels)
dist_long$individual2 <- factor(as.character(dist_long$individual2),
    levels = individual_levels)

population_levels <- union(as.character(dist_long$population1), as.character(dist_long$population2))
dist_long$population1 <- factor(as.character(dist_long$population1),
    levels = population_levels)
dist_long$population2 <- factor(as.character(dist_long$population2),
    levels = population_levels)

dist_long$same_population <- droplevels(dist_long$same_population)

saveRDS(dist_long, dist_long_file)

6.1 Data description

Code

pc_rule <- "broken_stick"
fit_suffix <- if (identical(pc_rule, "broken_stick")) "" else paste0("_",
    pc_rule)
dist_long <- readRDS(paste0("./data/processed/dist_long_song_type_1",
    fit_suffix, ".rds"))

dist_long$geo_km_adj <- exp(dist_long$log_geo_distance)
dist_long$acoustic_prop <- dist_long$acoustic_dissimilarity/max(dist_long$acoustic_dissimilarity,
    na.rm = TRUE)

ggplot(dist_long, aes(x = geo_km_adj, fill = same_population)) + geom_histogram(bins = 50,
    position = "identity", alpha = 0.6) + scale_x_log10() + scale_fill_viridis_d(begin = 0.2,
    end = 0.75, labels = c("Different", "Same")) + labs(x = "Geographic distance (km, log scale)",
    y = "Number of pairs", fill = "Population") + theme_classic()

Code

sp <- as.numeric(as.character(dist_long$same_population))
minpos <- min(dist_long$geo_distance[dist_long$geo_distance > 0],
    na.rm = TRUE)

cat("pairs:", nrow(dist_long), "| individuals:", length(unique(c(as.character(dist_long$individual1),
    as.character(dist_long$individual2)))), "\n")

pairs: 6555 | individuals: 115

Code

cat("geo_distance == 0:", sum(dist_long$geo_distance == 0), "(", round(100 *
    mean(dist_long$geo_distance == 0), 1), "% )\n")

geo_distance == 0: 386 ( 5.9 % )

Code

cat("log_geo range:", round(range(dist_long$log_geo_distance), 2),
    "\n")

log_geo range: 6.27 8.5

Code

cat("pairs in 100-1000 km gap:", sum(dist_long$geo_distance > 100 &
    dist_long$geo_distance < 1000), "(", round(100 * mean(dist_long$geo_distance >
    100 & dist_long$geo_distance < 1000), 1), "% )\n")

pairs in 100-1000 km gap: 1858 ( 28.3 % )

Code

cat("cor(log_geo, same_population):", round(cor(dist_long$log_geo_distance,
    sp), 3), "\n")

cor(log_geo, same_population): -0.783

Code

cat("cor(time_separation, log_geo):", round(cor(dist_long$time_separation,
    dist_long$log_geo_distance), 3), "\n")

cor(time_separation, log_geo): -0.008

Code

pos <- dist_long$geo_distance[dist_long$geo_distance > 0]
consts <- c(min(pos)/10, min(pos), quantile(dist_long$geo_distance,
    0.25), 1)
names(consts) <- c("min/10 (previous)", "min", "P25 (HH20, used)",
    "1 km")

6.2 Check for collinearity

Collinearity among all predictors

Code

pred_vars <- c("min_snr_sc", "same_population", "time_separation_sc",
    "geo_distance_sc")

pred_df <- dist_long[, pred_vars]
pred_df$same_population <- as.numeric(as.character(pred_df$same_population))

corr_mat <- cor(pred_df, use = "pairwise.complete.obs")

cols <- mako(10, alpha = 0.8, begin = 0.2, end = 0.75)

ggcorrplot(corr_mat, type = "lower", lab = TRUE, lab_size = 4, colors = c(cols[1],
    "white", cols[10]), title = "Pairwise predictor correlations")

Variance inflation factors (VIFs) were computed for the fixed effects in the model to check for collinearity. The VIFs are interpreted as follows: values below 5 indicate low concern, values between 5 and 10 indicate moderate concern, and values above 10 indicate high concern. In this case, all VIFs are below 5, indicating that collinearity is not a major concern for the predictors in the model.

Code

vif_mod <- lm(acoustic_dissimilarity_sc ~ min_snr_sc + same_population +
    time_separation_sc + geo_distance_sc, data = dist_long)

vif_vals <- car::vif(vif_mod)
vif_df <- data.frame(term = names(vif_vals), vif = as.numeric(vif_vals))

## conventional rule-of-thumb thresholds: <5 low concern, 5-10
## moderate, >10 high. Some fields use stricter cutoffs (2.5 /
## 4); doesn't matter for your numbers specifically, since
## ~1.7-2.2 clears either convention.
vif_df$severity <- cut(vif_df$vif, breaks = c(0, 5, 10, Inf), labels = c("low",
    "moderate", "high"))

ggplot(vif_df, aes(x = vif, y = reorder(term, vif), color = severity)) +
    geom_vline(xintercept = c(5, 10), linetype = "dashed", color = "grey60") +
    geom_segment(aes(x = 1, xend = vif, yend = term), linewidth = 1) +
    geom_point(size = 3) + scale_color_manual(values = c(low = cols[1],
    moderate = "orange", high = cols[10])) + labs(x = "Variance Inflation Factor",
    y = NULL, color = "Collinearity") + theme_minimal(base_size = 13)

7 Statistical analysis

To evaluate whether acoustic divergence among Bearded Bellbird “bock” songs is related to geographic distance, temporal separation between recordings, and population identity, we fitted a single Bayesian mixed-effects model on pairwise acoustic dissimilarities, accounting for the non-independence inherent to pairwise distance data.

Acoustic structure for “bock” songs was summarized with a principal component analysis (PCA) on standardized MFCC-derived acoustic parameters; see the Principal component analysis section above for the full rationale behind the retention rule used. PC scores were averaged across selections belonging to the same individual, and pairwise Euclidean distances in this reduced acoustic space were used as the response variable.

Pairwise geographic distance was calculated as the great-circle (haversine) distance between individuals, log-transformed prior to standardization to correct for the right-skewed distribution typical of pairwise spatial data (many nearby pairs, few very distant pairs). We use haversine great-circle distance because latitude/longitude are angles on a sphere rather than flat coordinates, and your populations span a wide enough latitude range that treating them as flat x/y would systematically distort distances in a way that correlates with the very population layout you’re testing. Because individuals recorded at the same locality were assigned identical coordinates, a small proportion of pairs have a geographic distance of exactly zero; following Hausdorf and Hennig (2020), the 25th percentile of all pairwise distances was added before log-transformation, rather than a constant derived from the smallest observed positive distance, which would place that point mass far into the lower tail of the predictor and give it disproportionate leverage on the estimated slope (see Check distances above). Pairwise temporal distance was calculated as the absolute difference in recording year. Both distance measures, along with the acoustic response, were standardized (z-scores) prior to modeling.

Recording quality was included as a nuisance covariate. In a pairwise-distance framework, a fixed recording quality specific effect does not cancel out of a pair’s distance but enters it as the difference between the two recordings’ conditions, so uneven quality across individuals is a potential source of spurious geographic or temporal signal rather than simply extra noise – a concern that applies with particular force to temporal separation, since older recordings are close to guaranteed to have used an older generation of equipment. Signal-to-noise ratio (SNR) was measured for every selection with sig2noise() (warbleR; Araya-Salas and Smith-Vidaurre 2017), over the same frequency band used for the MFCC statistics above, and averaged across selections belonging to the same individual. The covariate entering the model for each pair was the minimum, rather than the mean, of the two individuals’ SNR: a pair’s acoustic dissimilarity is corrupted by whichever recording is noisier, not by the pair’s average quality, so the minimum better reflects the measurement error affecting that specific comparison. SNR was standardized (z-score) prior to modeling.

The model was specified as:

\[ \begin{split} \text{acoustic dissimilarity}_{ij} &\sim \text{min SNR}_{ij} + \text{same population}_{ij} \\ &\quad + \text{geographic distance}_{ij} \times \text{temporal distance}_{ij} \\ &\quad + (1 \mid \text{mm(individual}_i,\text{individual}_j)) \\ &\quad + (1 \mid \text{mm(population}_i,\text{population}_j)) \end{split} \]

where:

\(\text{acoustic dissimilarity}_{ij}\) is the Euclidean distance between the songs of individuals \(i\) and \(j\) in the reduced PCA space.
\(\text{SNR}_{ij}\) is the standardized minimum signal-to-noise ratio between individuals \(i\) and \(j\), included as a nuisance covariate for recording quality.
\(\text{same population}_{ij}\) is a binary predictor indicating whether individuals \(i\) and \(j\) belong to the same population (1) or to different populations (0).
\(\text{geographic distance}_{ij}\) is the log-standardized geographic distance between individuals \(i\) and \(j\).
\(\text{temporal distance}_{ij}\) is the standardized difference in recording year between individuals \(i\) and \(j\).
mm(individual\(_i\), individual\(_j\)) is a multi-membership random intercept accounting for the repeated use of the same individuals across pairwise comparisons.
mm(population\(_i\), population\(_j\)) is a multi-membership random intercept accounting for the repeated use of the same populations across pairwise comparisons.

The interaction between geographic and temporal distance was included to model the occurrence of vocal divergence through time even within a small geographic range.

Model specifications

The model was fitted in a Bayesian framework using the brms package.
The response variable was pairwise Euclidean acoustic dissimilarity and was modeled with a Gaussian error distribution.
Individual identity was modeled as a multi-membership random effect because each pairwise comparison involves two individuals, each of which contributes to multiple pairwise distances.
Population identity was modeled as a multi-membership random effect for the same reason.
Only unique pairwise comparisons were retained and self-comparisons were excluded.
Recording quality (minimum pairwise SNR, standardized) was included as a fixed-effect covariate, so that variance attributable to recording conditions is not misattributed to geography, time, or population identity.
Mildly regularizing priors were specified for all parameters and scaled to the standardized response: Normal(0, 1) for the intercept, Normal(0, 0.5) for fixed effects, Exponential(2) for random-effect SDs, and Exponential(1) for the residual SD. Because the response and both distance predictors are z-scored, standardized coefficients with |β| > 1 are effectively impossible and random-effect SDs cannot plausibly approach the total SD of the response; the fixed-effect and random-effect SD priors place their mass accordingly. This matters most for the population-level term, which is informed by only four populations and is therefore substantially prior-driven. These priors are checked against the model’s own implied predictions, before the data are seen, in the Prior predictive check below.
The model was fitted using Hamiltonian Monte Carlo as implemented in Stan through the cmdstanr backend, with within-chain thread parallelization, 4 chains, 4 cores, and 4,000 iterations per chain.

7.1 Model fitting

The pairwise dataset built and saved in Check distances above is loaded here, together with the priors described earlier, and used to fit two models: a prior-only version, used below for the prior predictive check, and the model fitted to the observed data. Both brm() calls are cached to disk via file =; file_refit = "on_change" ensures a stale cache is never silently reused if the priors, data, or formula change without the file name changing.

Code

## dist_long is built once, in 'Check distances' above, and
## saved to disk there; loading it here (rather than rebuilding
## it, or relying on it still being in memory from an earlier
## chunk) guarantees the model is fitted to exactly the same
## pairwise dataset the diagnostics above were run on,
## regardless of the order chunks happen to be run in.

pc_rule <- "broken_stick"
fit_suffix <- if (identical(pc_rule, "broken_stick")) "" else paste0("_",
    pc_rule)

pca_file <- paste0("./data/processed/pca_song_type_1", fit_suffix,
    ".rds")
res_file <- paste0("./data/processed/results_song_type_1", fit_suffix,
    ".rds")

pca_out <- readRDS(pca_file)
n_pcs <- pca_out$n_pcs
cum_var <- pca_out$cum_var

dist_long <- readRDS(paste0("./data/processed/dist_long_song_type_1",
    fit_suffix, ".rds"))

## geo_const was already used to build geo_distance_sc inside
## dist_long above; recomputed here only so it can be stored in
## results_type1 for reporting, not because it affects the data
## in any way.
geo_const <- unname(quantile(dist_long$geo_distance, 0.25, na.rm = TRUE))

## Priors rescaled for standardized response and predictors.
## Response and predictors are z-scored (SD = 1), so |beta| > 1
## is near-impossible and random-effect SDs cannot plausibly
## exceed the total SD of the response.  With only 4
## populations, the population-level sd is only weakly
## identified by the data and the prior does real work there.
## What these priors imply, before the data are seen, is checked
## directly in the 'Prior predictive check' section below.
priors <- c(prior(normal(0, 1), class = "Intercept"), prior(normal(0,
    0.5), class = "b"), prior(exponential(2), class = "sd"), prior(exponential(1),
    class = "sigma"))


fit_path <- function(name) paste0("./data/processed/fits/", name,
    "_song_type_1", fit_suffix)

## ============================================================
## 6. PRIOR PREDICTIVE CHECK MODEL
## ============================================================
## sample_prior = 'only' draws parameters from the priors alone
## and ignores the observed response entirely; it still needs
## the predictor columns and formula/family to build the model
## matrix. See 'Prior predictive check' below for what this is
## diagnostic of.
mod_prior_only <- brm(acoustic_dissimilarity_sc ~ min_snr_sc + same_population +
    time_separation_sc * geo_distance_sc + (1 | mm(individual1, individual2)) +
    (1 | mm(population1, population2)), data = dist_long, family = gaussian(),
    prior = priors, sample_prior = "only", chains = 4, cores = 4,
    iter = 4000, backend = "cmdstanr", threads = threading(8), control = list(adapt_delta = 0.95,
        max_treedepth = 15), file = fit_path("mod_prior_only"), file_refit = "on_change")
## ============================================================
## 7. THE MODEL
## ============================================================
## Every predictor implied by the study's hypotheses in one
## model, fitted to the observed data. See 'Statistical
## analysis' above for why this replaces a comparison across
## candidate models.

mod_interaction <- brm(acoustic_dissimilarity_sc ~ min_snr_sc + same_population +
    time_separation_sc * geo_distance_sc + (1 | mm(individual1, individual2)) +
    (1 | mm(population1, population2)), data = dist_long, family = gaussian(),
    prior = priors, chains = 4, cores = 4, iter = 4000, backend = "cmdstanr",
    threads = threading(8), control = list(adapt_delta = 0.95, max_treedepth = 15),
    file = fit_path("mod_interaction"), file_refit = "on_change")

results_type1 <- list(n_pcs = n_pcs, pc_rule = pc_rule, var_explained = cum_var[n_pcs],
    geo_const = geo_const, dist_long = dist_long, prior_check_model = mod_prior_only,
    model = mod_interaction)

saveRDS(results_type1, res_file)

7.1.1 Prior predictive check

Before looking at what the observed data say, we checked what the priors specified above imply on their own: parameter values were drawn from the priors and used to generate synthetic acoustic dissimilarities (sample_prior = "only" in brms), using the real predictor values but none of the information in the actual response. This step catches priors that look reasonable individually but combine, once passed through the full model, to imply an implausible range of outcomes – for example, if the fixed effects and random-effect SDs jointly implied standardized acoustic dissimilarities routinely several SDs from zero, that would signal the priors are too permissive before any conclusion has been drawn from real data.

One feature of this check is worth anticipating rather than mistaking for a problem. Acoustic dissimilarity is a non-negative quantity with a hard floor at zero, while the model uses a Gaussian likelihood on the z-scored response; a Gaussian prior (and posterior) predictive distribution will therefore assign some mass to negative standardized distances regardless of how the priors are tuned. This is a structural feature shared with related additive dyadic-distance models (e.g. MLPE; Clarke et al. 2002), not evidence of a misspecified prior. What the plot below is actually diagnostic of is whether the bulk of the prior predictive distribution covers a plausible range for a standardized response – roughly within a few SDs of zero – rather than being wildly over- or under-dispersed.

Code

results_type1 <- readRDS("./data/processed/results_song_type_1.rds")
mod_prior_only <- results_type1$prior_check_model

## Does the bulk of the prior-predictive distribution of
## standardized acoustic dissimilarity stay within a plausible
## range, rather than being wildly over- or under-dispersed?
## Some negative-distance mass is expected here (see prose
## above) and is not itself a problem.
pp_check(mod_prior_only, ndraws = 100) + labs(title = "Prior predictive check")

7.1.2 Posterior predictive check

With the prior predictive check above informing whether the priors themselves are reasonable, the model fitted to the observed data is next checked against the data it was fitted to, in two complementary ways. An overall density overlay checks whether the model reproduces the general shape of the observed acoustic-distance distribution, including its right tail (subject to the same floor-at-zero caveat noted above). A grouped check by population asks a more targeted question: does the model reproduce each population’s mean acoustic dissimilarity, both within that population and against every other population? This speaks directly to whether mm(population1, population2) is capturing real population-level structure rather than averaging over it.

Code

mod_interaction <- results_type1$model
dist_long <- results_type1$dist_long

## Overall shape of the observed vs. predicted acoustic-distance
## distribution.
pp_check(mod_interaction, ndraws = 100, type = "dens_overlay") + ggplot2::labs(title = "Posterior predictive check -- bock songs")

Code

## does the model reproduce each subgroup's shape, or does
## pooling same/different-population pairs together wash out
## something the model isn't capturing?
pp_check(mod_interaction, type = "dens_overlay_grouped", group = "same_population",
    ndraws = 100)

Code

## Does the model reproduce each population's mean acoustic
## dissimilarity, within-population and against every other
## population?
pp_check(mod_interaction, type = "stat_grouped", group = "population1",
    stat = "mean")

7.1.3 Results and biological interpretation

Code

pc_rule <- "broken_stick"
fit_suffix <- if (identical(pc_rule, "broken_stick")) "" else paste0("_",
    pc_rule)

results_type1 <- readRDS(paste0("./data/processed/results_song_type_1",
    fit_suffix, ".rds"))
mod_interaction <- results_type1$model
dist_long <- results_type1$dist_long

cat("\n\n==========================\n bock songs", "| PC rule:", as.character(results_type1$pc_rule),
    "| PCs used:", results_type1$n_pcs, "| Cumulative variance:",
    round(results_type1$var_explained * 100, 1), "%\n")



==========================
 bock songs | PC rule: broken_stick | PCs used: 14 | Cumulative variance: 58.4 %

Code

summ <- extended_summary(mod_interaction, highlight = TRUE, trace = FALSE,
    remove.intercepts = TRUE, gsub.pattern = c("b_min_snr_sc", "b_same_population1",
        "b_time_separation_sc", "b_geo_distance_sc", "Temporal\nseparation:geo_distance_sc"),
    gsub.replacement = c("Minimum\nSNR", "Same\npopulation", "Temporal\nseparation",
        "Geographic\ndistance", "Temporal separation ×\ngeographic distance"),
    return = TRUE)

summ$coef_table_html

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
Minimum SNR	-0.341	-0.407	-0.276	1.003	2011.256	3162.646
Same population	-0.707	-0.782	-0.633	1.001	2725.597	3980.645
Temporal separation	0.054	0.006	0.102	1.001	2896.391	4482.631
Geographic distance	0.040	-0.033	0.113	1.001	2644.687	3948.011
Temporal separation × geographic distance	-0.273	-0.328	-0.219	1.001	4939.292	5241.067

Code

## ---- gather posterior draws for the fixed effects (excluding intercept) --

draws <- as_draws_df(mod_interaction)

fx_cols <- grep("^b_", names(draws), value = TRUE)
fx_cols <- setdiff(fx_cols, "b_Intercept")

draws_long <- do.call(rbind, lapply(fx_cols, function(v) {
  data.frame(.variable = v, .value = draws[[v]])
}))

## ---- point estimate + 95% CI per parameter --------------------------------
## median_qi matches stat_halfeye()'s own default point_interval, so the
## printed numbers are exactly what the plotted point/interval represent.

summary_df <- do.call(rbind, lapply(fx_cols, function(v) {
  qs <- quantile(draws[[v]], c(0.025, 0.5, 0.975))
  data.frame(.variable = v, estimate = qs[2], lower = qs[1], upper = qs[3])
}))
rownames(summary_df) <- NULL

summary_df$label_text <- sprintf("%.2f [%.2f, %.2f]",
                                  summary_df$estimate, summary_df$lower, summary_df$upper)

## ---- order by |estimate|, high to low -------------------------------------
## A discrete y-axis reads factor levels bottom-to-top, so the LARGEST
## |estimate| has to be the LAST level to land at the top of the plot.

summary_df  <- summary_df[order(abs(summary_df$estimate)), ]
param_order <- summary_df$.variable

## ---- custom labels, keyed to the actual term names ------------------------

label_map <- setNames(
  c("Minimum\nSNR", "Same\npopulation", "Temporal\nseparation",
    "Geographic\ndistance", "Temporal\nseparation\n×\ngeographic\ndistance"),
  c("b_min_snr_sc", "b_same_population1", "b_time_separation_sc",
    "b_geo_distance_sc", "b_time_separation_sc:geo_distance_sc")
)

draws_long$.variable <- factor(draws_long$.variable, levels = param_order, labels = label_map[param_order])
summary_df$.variable <- factor(summary_df$.variable, levels = param_order, labels = label_map[param_order])

## ---- plot ------------------------------------------------------------------

summary_df$nudge_y <- c(0.45, 0.6, 0.6, 0.47, 0.43)  # adjust as needed for label placement)

cat("Probability of direction")

Probability of direction

Code

p_direction(mod_interaction, effects = "fixed")

	Parameter	pd	Effects	Component
2	b_Intercept	0.550750	fixed	conditional
3	b_min_snr_sc	1.000000	fixed	conditional
4	b_same_population1	1.000000	fixed	conditional
5	b_time_separation_sc	0.987250	fixed	conditional
1	b_geo_distance_sc	0.861375	fixed	conditional
6	b_time_separation_sc:geo_distance_sc	1.000000	fixed	conditional

Code

gg_estimates <- ggplot(draws_long, aes(x = .value, y = .variable)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray") +
  stat_halfeye(
    fill       = cols[9],
    color      = "black",
    slab_alpha = 0.6,
    point_size = 2,
    scale      = 0.6   # caps hill height, leaving guaranteed headroom for the text below
  ) +
  geom_text(
    data          = summary_df,
    aes(x = estimate, y = .variable, label = label_text,  nudge_y = nudge_y),
    inherit.aes   = FALSE,
    vjust         = 0,
   
    size          = 3
  ) +
  scale_y_discrete(expand = expansion(add = c(0.5, 0.8))) +
  labs(x = "Change in acoustic dissimilarity (SD)", y = "Predictor") +
  theme_classic(base_size = 13) + xlim(-0.9, 0.24)

gg_estimates

The figure below shows the estimated effect of geographic distance on acoustic divergence as a function of temporal separation between recordings, computed directly from the joint posterior distribution of the geographic-distance and geographic-by-temporal-distance interaction coefficients rather than from a separate model. At each value of temporal distance, the simple slope of geographic distance was calculated as the sum of the main geographic-distance coefficient and the interaction coefficient scaled by that temporal-distance value — the standard decomposition of an interaction term into its conditional slopes. The solid line shows the posterior median of this slope across the observed range of temporal separation, and the shaded band gives the corresponding 95% credible interval; the dashed horizontal line marks zero, the threshold the credible interval must clear to indicate a geographic-distance effect distinguishable from no effect at a given point along the temporal range.

Code

## ---- posterior draws for the two terms that define the
## geo-distance slope

draws <- as_draws_df(mod_interaction)
b_geo <- draws$b_geo_distance_sc
b_int <- draws$`b_time_separation_sc:geo_distance_sc`

## ---- sequence of time_separation across the OBSERVED range
## only ------------ Raw years for the x-axis (more
## interpretable than z-scores); converted to the same
## standardization used to fit the model so the slope is
## computed on the scale the coefficients actually apply to.
## Kept inside the observed range deliberately -- extrapolating
## the line past it would imply precision the data don't
## support, especially given how thin the sample gets at the
## extremes (Venezuela SE in particular).

time_mean <- mean(dist_long$time_separation)
time_sd <- sd(dist_long$time_separation)

time_seq_raw <- seq(min(dist_long$time_separation), max(dist_long$time_separation),
    length.out = 100)
time_seq_sc <- (time_seq_raw - time_mean)/time_sd

## ---- simple slope of geo_distance_sc at every point in that
## sequence ----- Same identity as the earlier p10/p50/p90
## table, just evaluated over a continuous grid instead of three
## quantiles: slope = b_geo + b_int * time.

slope_draws <- outer(b_geo, rep(1, length(time_seq_sc))) + outer(b_int,
    time_seq_sc)

slope_df <- data.frame(time_separation = time_seq_raw, estimate = apply(slope_draws,
    2, median), lower = apply(slope_draws, 2, quantile, 0.025), upper = apply(slope_draws,
    2, quantile, 0.975))

## ---- plot
## ------------------------------------------------------------------

gg_inter <- ggplot(slope_df, aes(x = time_separation, y = estimate)) +
    geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
    geom_ribbon(aes(ymin = lower, ymax = upper), fill = cols[9], alpha = 0.25) +
    geom_line(color = cols[9], linewidth = 1) + labs(x = "Temporal separation (years)",
    y = "Effect of geographic distance\non acoustic dissimilarity") +
    theme_classic(base_size = 13)

gg_inter

Code

# combine gg_estimates and gg_inter into a single figure
gg_combined <- ggpubr::ggarrange(gg_estimates, gg_inter, ncol = 2,
    nrow = 1, labels = c("A", "B"), heights = c(1, 1.2))

ggsave(filename = paste0("./output/fig_results_bock", fit_suffix,
    ".png"), plot = gg_combined, device = grDevices::png, width = 10,
    height = 5, units = "in", dpi = 300)

Code

## ---- load the fitted model
## ------------------------------------------------

pc_rule <- "broken_stick"
fit_suffix <- if (identical(pc_rule, "broken_stick")) "" else paste0("_",
    pc_rule)

results_type1 <- readRDS(paste0("./data/processed/results_song_type_1",
    fit_suffix, ".rds"))
mod_interaction <- results_type1$model

## ---- pull each predictor's row from fixef()
## -------------------------------

fx <- as.data.frame(fixef(mod_interaction))

get_row <- function(pattern) {
    hit <- fx[grepl(pattern, rownames(fx)), , drop = FALSE]
    hit[1, ]
}

snr_row <- get_row("^min_snr_sc$")
samepop_row <- get_row("^same_population")
geo_row <- get_row("^geo_distance_sc$")
time_row <- get_row("^time_separation_sc$")
int_row <- get_row("time_separation_sc:geo_distance_sc")

## ---- credible = 95% CI excludes zero
## --------------------------------------

credible <- function(row) sign(row$Q2.5) == sign(row$Q97.5)

snr_credible <- credible(snr_row)
samepop_credible <- credible(samepop_row)
geo_credible <- credible(geo_row)
time_credible <- credible(time_row)
int_credible <- credible(int_row)

## ---- wording helpers used inline in the text
## ------------------------------

dir_word <- function(row) if (row$Estimate > 0) "increased" else "decreased"
cred_word <- function(is_credible) if (is_credible) "a credible" else "no credible"

## ---- posterior probability the interaction coefficient is
## positive -------

hyp_int <- hypothesis(mod_interaction, "time_separation_sc:geo_distance_sc > 0")

interaction_sentence <- if (int_credible) {
    "The interaction was credible (its 95% credible interval excludes zero), indicating that the relationship between geographic separation and acoustic divergence itself depends on how much time has passed between the compared recordings -- the specific signature predicted by the cultural drift hypothesis, under which vocal traditions in isolated populations wander apart at a rate that compounds with both distance and time rather than either acting alone."
} else {
    "The interaction's 95% credible interval included zero, providing no support, at the statistical power available here, for acoustic divergence compounding jointly with distance and time. This does not rule out cultural drift outright: geography and/or time may still each contribute independently (see their respective coefficients above), but the data do not support the stronger claim that isolation must persist over time to produce divergence beyond what distance and time contribute on their own."
}

Overall patterns

A single, fully specified model was used for inference, rather than a comparison across candidate models (see Statistical analysis for the rationale). Every coefficient below is interpreted while holding all other predictors constant.
Recording quality (minimum pairwise SNR) had a credible association with acoustic dissimilarity (β = -0.341, 95% CrI [-0.407, -0.276]), consistent with noisier recordings inflating apparent acoustic divergence and supporting its inclusion as a covariate.
Belonging to the same population was associated with decreased acoustic dissimilarity relative to belonging to different populations (β = -0.707, 95% CrI [-0.782, -0.633]), a credible effect.
Geographic distance was associated with increased acoustic dissimilarity (β = 0.04, 95% CrI [-0.033, 0.113]), though the credible interval included zero.
Temporal distance on its own showed a credible effect on acoustic dissimilarity (β = 0.054, 95% CrI [0.006, 0.102]).
The geographic × temporal distance interaction – the coefficient that most directly tests the cultural drift hypothesis – showed a credible effect (β = -0.273, 95% CrI [-0.328, -0.219]), with a posterior probability of 0% that the coefficient is positive.

Overall biological interpretation

Acoustic divergence in Bearded Bellbird “bock” songs is shaped by a combination of geography, population identity, and recording quality, rather than by any single overwhelming factor.
Geographic distance alone did not show a credible association with acoustic divergence once the other predictors were accounted for.
Population identity contributed explanatory power beyond geographic distance alone, suggesting that population boundaries carry some additional structure – a dispersal barrier or habitat discontinuity, for instance – beyond what raw geographic distance captures.
The interaction was credible (its 95% credible interval excludes zero), indicating that the relationship between geographic separation and acoustic divergence itself depends on how much time has passed between the compared recordings – the specific signature predicted by the cultural drift hypothesis, under which vocal traditions in isolated populations wander apart at a rate that compounds with both distance and time rather than either acting alone.
Recording quality had a detectable influence on measured acoustic dissimilarity, confirming that controlling for it was a necessary methodological step rather than a precaution with no practical consequence. This is a data-quality observation rather than a biological one.
Taken together, these results indicate that acoustic divergence in this species compounds with both geographic isolation and the time elapsed since populations diverged. This case study illustrates how recordings assembled programmatically with suwo can be carried through a complete geographic-variation analysis, from raw media to a fully specified, checked, and interpreted statistical model.

References

Snow, B. K. (1970). A field study of the Bearded Bellbird in Trinidad. Ibis, 112(3), 299-329.
Araya‐Salas, M., Smith‐Vidaurre, G., & Golding, N. (2017). warbleR: an r package to streamline analysis of animal acoustic signals. Methods in Ecology & Evolution, 8(2).

Session information

Click to see

─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.6.1 (2026-06-24)
 os       Ubuntu 22.04.5 LTS
 system   x86_64, linux-gnu
 ui       X11
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       America/Costa_Rica
 date     2026-07-30
 pandoc   3.8.3 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/x86_64/ (via rmarkdown)
 quarto   1.7.31 @ /usr/local/bin/quarto

─ Packages ───────────────────────────────────────────────────────────────────
 package        * version    date (UTC) lib source
 abind            1.4-8      2024-09-12 [1] CRAN (R 4.6.1)
 ape              5.8-1      2024-12-16 [1] CRAN (R 4.6.1)
 arrayhelpers     1.1-1      2026-07-05 [1] CRAN (R 4.6.1)
 backports        1.5.1      2026-04-03 [1] CRAN (R 4.6.1)
 bayesplot        1.15.0     2025-12-12 [1] CRAN (R 4.6.1)
 bayestestR     * 0.18.1     2026-05-24 [1] CRAN (R 4.6.1)
 bitops           1.0-9      2024-10-03 [1] CRAN (R 4.6.1)
 bridgesampling   1.2-1      2025-11-19 [1] CRAN (R 4.6.1)
 brio             1.1.5      2024-04-24 [1] CRAN (R 4.6.1)
 brms           * 2.23.0     2025-09-09 [1] CRAN (R 4.6.1)
 brmsish        * 1.0.0      2026-07-14 [1] Github (maRce10/brmsish@81ab826)
 Brobdingnag      1.2-9      2022-10-19 [1] CRAN (R 4.6.1)
 broom            1.0.13     2026-05-14 [1] CRAN (R 4.6.1)
 cachem           1.1.0      2024-05-16 [1] CRAN (R 4.6.1)
 car              3.1-5      2026-02-03 [1] CRAN (R 4.6.1)
 carData          3.0-6      2026-01-30 [1] CRAN (R 4.6.1)
 checkmate        2.3.4      2026-02-03 [1] CRAN (R 4.6.1)
 cli              3.6.6      2026-04-09 [1] CRAN (R 4.6.1)
 cmdstanr         0.9.0.9001 2026-07-14 [1] Github (stan-dev/cmdstanr@2763c0f)
 coda             0.19-4.1   2024-01-31 [1] CRAN (R 4.6.1)
 codetools        0.2-20     2024-03-31 [1] CRAN (R 4.6.1)
 cowplot          1.2.0      2025-07-07 [1] CRAN (R 4.6.1)
 crayon           1.5.3      2024-06-20 [1] CRAN (R 4.6.1)
 crosstalk        1.2.2      2025-08-26 [1] CRAN (R 4.6.1)
 curl             7.1.0      2026-04-22 [1] CRAN (R 4.6.1)
 datawizard       1.3.1      2026-04-26 [1] CRAN (R 4.6.1)
 devtools         2.5.2      2026-04-30 [1] CRAN (R 4.6.1)
 digest           0.6.39     2025-11-19 [1] CRAN (R 4.6.1)
 distributional   0.8.1      2026-06-27 [1] CRAN (R 4.6.1)
 dplyr            1.2.1      2026-04-03 [1] CRAN (R 4.6.1)
 dtw              1.23-3     2026-06-09 [1] CRAN (R 4.6.1)
 ellipsis         0.3.3      2026-04-04 [1] CRAN (R 4.6.1)
 emmeans          2.0.4      2026-07-15 [1] CRAN (R 4.6.1)
 estimability     2.0.0      2026-06-26 [1] CRAN (R 4.6.1)
 evaluate         1.0.5      2025-08-27 [1] CRAN (R 4.6.1)
 farver           2.1.2      2024-05-13 [1] CRAN (R 4.6.1)
 fastmap          1.2.0      2024-05-15 [1] CRAN (R 4.6.1)
 fftw             1.0-9      2024-09-20 [1] CRAN (R 4.6.1)
 formatR          1.14       2023-01-17 [1] CRAN (R 4.6.1)
 Formula          1.2-5      2023-02-24 [1] CRAN (R 4.6.1)
 fs               2.1.0      2026-04-18 [1] CRAN (R 4.6.1)
 generics         0.1.4      2025-05-09 [1] CRAN (R 4.6.1)
 geosphere      * 1.6-8      2026-04-05 [1] CRAN (R 4.6.1)
 ggcorrplot     * 0.3.0      2026-07-24 [1] CRAN (R 4.6.1)
 ggdist         * 3.3.3      2025-04-23 [1] CRAN (R 4.6.1)
 ggplot2        * 4.0.3      2026-04-22 [1] CRAN (R 4.6.1)
 ggpubr           1.0.0      2026-07-06 [1] CRAN (R 4.6.1)
 ggsignif         0.6.4      2022-10-13 [1] CRAN (R 4.6.1)
 ggtext         * 0.1.2      2022-09-16 [1] CRAN (R 4.6.1)
 glue             1.8.1      2026-04-17 [1] CRAN (R 4.6.1)
 gridExtra        2.3.1      2026-06-25 [1] CRAN (R 4.6.1)
 gridtext         0.1.6      2026-02-19 [1] CRAN (R 4.6.1)
 gtable           0.3.6      2024-10-25 [1] CRAN (R 4.6.1)
 htmltools        0.5.9      2025-12-04 [1] CRAN (R 4.6.1)
 htmlwidgets      1.6.4      2023-12-06 [1] CRAN (R 4.6.1)
 httr             1.4.8      2026-02-13 [1] CRAN (R 4.6.1)
 inline           0.3.21     2025-01-09 [1] CRAN (R 4.6.1)
 insight          1.5.2      2026-06-28 [1] CRAN (R 4.6.1)
 jquerylib        0.1.4      2021-04-26 [1] CRAN (R 4.6.1)
 jsonlite         2.0.0      2025-03-27 [1] CRAN (R 4.6.1)
 kableExtra       1.4.1      2026-07-08 [1] CRAN (R 4.6.1)
 knitr          * 1.51       2025-12-20 [1] CRAN (R 4.6.1)
 labeling         0.4.3      2023-08-29 [1] CRAN (R 4.6.1)
 lattice          0.22-9     2026-02-09 [3] CRAN (R 4.5.2)
 leaflet          2.2.3      2025-09-04 [1] CRAN (R 4.6.1)
 lifecycle        1.0.5      2026-01-08 [1] CRAN (R 4.6.1)
 loo              2.10.0     2026-06-26 [1] CRAN (R 4.6.1)
 magrittr         2.0.5      2026-04-04 [1] CRAN (R 4.6.1)
 MASS             7.3-66     2026-07-15 [1] CRAN (R 4.6.1)
 Matrix           1.7-5      2026-03-21 [3] CRAN (R 4.5.3)
 matrixStats      1.5.0      2025-01-07 [1] CRAN (R 4.6.1)
 memoise          2.0.1      2021-11-26 [1] CRAN (R 4.6.1)
 mvtnorm          1.4-2      2026-07-12 [1] CRAN (R 4.6.1)
 NatureSounds   * 1.0.5      2025-01-17 [1] CRAN (R 4.6.1)
 nlme             3.1-169    2026-03-27 [3] CRAN (R 4.5.3)
 otel             0.2.0      2025-08-29 [1] CRAN (R 4.6.1)
 packrat          0.9.3      2025-06-16 [1] CRAN (R 4.6.1)
 pbapply          1.7-4      2025-07-20 [1] CRAN (R 4.6.1)
 pillar           1.11.1     2025-09-17 [1] CRAN (R 4.6.1)
 pkgbuild         1.4.8      2025-05-26 [1] CRAN (R 4.6.1)
 pkgconfig        2.0.3      2019-09-22 [1] CRAN (R 4.6.1)
 pkgload          1.5.3      2026-06-15 [1] CRAN (R 4.6.1)
 plyr             1.8.9      2023-10-02 [1] CRAN (R 4.6.1)
 posterior      * 1.7.0      2026-04-01 [1] CRAN (R 4.6.1)
 processx         3.9.0      2026-04-22 [1] CRAN (R 4.6.1)
 proxy            0.4-29     2025-12-29 [1] CRAN (R 4.6.1)
 ps               1.9.3      2026-04-20 [1] CRAN (R 4.6.1)
 purrr            1.2.2      2026-04-10 [1] CRAN (R 4.6.1)
 QuickJSR         1.7.0      2025-03-31 [2] CRAN (R 4.5.0)
 R6               2.6.1      2025-02-15 [1] CRAN (R 4.6.1)
 RColorBrewer     1.1-3      2022-04-03 [1] CRAN (R 4.6.1)
 Rcpp           * 1.1.2      2026-07-05 [1] CRAN (R 4.6.1)
 RcppParallel     5.1.11-2   2026-03-05 [1] CRAN (R 4.6.1)
 RCurl            1.98-1.19  2026-06-03 [1] CRAN (R 4.6.1)
 remotes          2.5.0      2024-03-17 [1] CRAN (R 4.6.1)
 reshape2         1.4.5      2025-11-12 [1] CRAN (R 4.6.1)
 rjson            0.2.23     2024-09-16 [1] CRAN (R 4.6.1)
 rlang            1.3.0      2026-07-05 [1] CRAN (R 4.6.1)
 rmarkdown        2.31       2026-03-26 [1] CRAN (R 4.6.1)
 Rraven         * 1.0.16     2017-11-17 [1] Github (maRce10/Rraven@5fea9e1)
 rstan            2.32.7     2025-03-10 [1] CRAN (R 4.6.1)
 rstantools       2.6.0      2026-01-10 [1] CRAN (R 4.6.1)
 rstatix          1.1.0      2026-07-23 [1] CRAN (R 4.6.1)
 rstudioapi       0.19.0     2026-06-11 [1] CRAN (R 4.6.1)
 S7               0.2.2      2026-04-22 [1] CRAN (R 4.6.1)
 scales           1.4.0      2025-04-24 [1] CRAN (R 4.6.1)
 seewave        * 2.2.4      2025-08-19 [1] CRAN (R 4.6.1)
 sessioninfo      1.2.4      2026-06-04 [1] CRAN (R 4.6.1)
 signal           1.8-1      2024-06-26 [1] CRAN (R 4.6.1)
 sketchy        * 1.0.5      2025-01-16 [1] CRAN (R 4.6.1)
 StanHeaders      2.32.10    2024-07-15 [1] CRAN (R 4.6.1)
 stringi          1.8.7      2025-03-27 [1] CRAN (R 4.6.1)
 stringr          1.6.0      2025-11-04 [1] CRAN (R 4.6.1)
 suwo           * 0.2.2      2026-07-23 [1] local
 svglite          2.2.2      2025-10-21 [1] CRAN (R 4.6.1)
 svUnit           1.0.8      2025-08-26 [1] CRAN (R 4.6.1)
 systemfonts      1.3.2      2026-03-05 [1] CRAN (R 4.6.1)
 tensorA          0.36.2.1   2023-12-13 [1] CRAN (R 4.6.1)
 testthat         3.3.2      2026-01-11 [1] CRAN (R 4.6.1)
 textshaping      1.0.5      2026-03-06 [1] CRAN (R 4.6.1)
 tibble           3.3.1      2026-01-11 [1] CRAN (R 4.6.1)
 tidybayes        3.0.7      2024-09-15 [1] CRAN (R 4.6.1)
 tidyr            1.3.2      2025-12-19 [1] CRAN (R 4.6.1)
 tidyselect       1.2.1      2024-03-11 [1] CRAN (R 4.6.1)
 tuneR          * 1.4.7      2024-04-17 [1] CRAN (R 4.6.1)
 usethis          3.2.1      2025-09-06 [1] CRAN (R 4.6.1)
 V8               8.2.0      2026-04-21 [1] CRAN (R 4.6.1)
 vctrs            0.7.3      2026-04-11 [1] CRAN (R 4.6.1)
 viridis        * 0.6.5      2024-01-29 [1] CRAN (R 4.6.1)
 viridisLite    * 0.4.3      2026-02-04 [1] CRAN (R 4.6.1)
 warbleR        * 1.1.37     2025-10-22 [1] CRAN (R 4.6.1)
 withr            3.0.3      2026-06-19 [1] CRAN (R 4.6.1)
 xaringanExtra    0.8.0      2024-05-19 [1] CRAN (R 4.6.1)
 xfun             0.60       2026-07-09 [1] CRAN (R 4.6.1)
 xml2             1.6.0      2026-06-22 [1] CRAN (R 4.6.1)
 xtable           1.8-8      2026-02-22 [1] CRAN (R 4.6.1)
 yaml             2.3.12     2025-12-10 [1] CRAN (R 4.6.1)

 [1] /usr/local/lib/R/site-library
 [2] /usr/lib/R/site-library
 [3] /usr/lib/R/library
 * ── Packages attached to the search path.

──────────────────────────────────────────────────────────────────────────────