Song geographic variation of Bearded Bellbirds

Obtaining nature media with the R package suwo

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

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

• Explore variation in “bock” songs 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

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

# install/ load packages

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

1 Custom functions

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

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

## ============================================================
## build_dist_long(): construct the pairwise dataset from a
## selection table that already carries PC scores and SNR
## ============================================================
## Which INDIVIDUALS/RECORDINGS go in is controlled entirely by
## what's passed in as bock_est, not by any argument to this
## function itself.  Assumes bock_est already has PC score
## columns (from PCA) and an SNR column (computed earlier in the
## pipeline) -- neither is recomputed here.  Arguments: bock_est
## extended_selection_table with PC scores and SNR attached pcs
## character vector of PC column names (from pca_out$pcs)
## out_file optional path to saveRDS() the result to; NULL skips
## saving Returns: the pairwise dist_long data.frame.

build_dist_long <- function(bock_est, pcs, out_file = NULL) {

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

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

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

    agg_meta <- as.data.frame(bock_est[!duplicated(bock_est$orig.sound.files),
        c("orig.sound.files", "population", "latitude", "longitude",
            "date")])

    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

    ## Convert explicitly to numeric BEFORE dist()
    agg_indiv$full_date <- as.Date(agg_indiv$date)

    ## Fractional years, not raw days: keeps this on the same
    ## 'years' scale as every existing description of temporal
    ## separation in the manuscript, while replacing the coarse
    ## integer-year-difference with real day-level precision.
    agg_indiv$year_precise <- as.numeric(agg_indiv$full_date)/365.25

    cat("Individuals going into this build:", nrow(agg_indiv), "\n")
    print(table(agg_indiv$population))

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

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

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

    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 -- revisit if skewed once inspected.
    dist_long$min_snr_sc <- scale(dist_long$min_snr)[, 1]/2

    ## ---- Filter to complete cases
    ## -----------------------------------------

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

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

    cat("Final pairs in this build:", nrow(dist_long), "\n\n")

    if (!is.null(out_file))
        saveRDS(dist_long, out_file)

    dist_long
}

2 Obtain recordings

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

# 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

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)

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

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

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

# remove those without a date
anns <- anns[!is.na(anns$date), ]

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

est$year <-  anns$year
est$date <-  anns$date

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

attr(est, "metadata") <- metadata

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

bock_est <- est[est$Tipo == 1, ]

bock_est$Tipo <- NULL

# measure SNR
bock_est <- sig2noise(bock_est, mar = 0.05,
                     bp = c(min(bock_est$bottom.freq), max(bock_est$top.freq)))


attr(bock_est, "metadata") <- metadata[metadata$downloaded_file_name %in% bock_est$orig.sound.files, ]

saveRDS(bock_est, "./data/processed/p_averano_bock_extended_selection_table.rds")

4 Remove low SNR songs

Exclude songs with SNR < 1, which are likely to be too noisy to provide reliable acoustic measurements.

bock_est <- readRDS("./data/processed/p_averano_bock_extended_selection_table.rds")

bock_metadata <- attr(bock_est, "metadata")

# make ggplot histogram
ggplot(bock_est, aes(x = SNR)) + geom_histogram(bins = 50, fill = viridis::mako(10,
    alpha = 0.5)[2], color = "black") + labs(x = "Signal-to-noise ratio (dB)",
    y = "Count") + geom_vline(xintercept = 1, color = "red", linetype = "dashed") +
    theme_classic()

8 out of 781 songs have SNR < 1 and were excluded from the analysis.

bock_est <- bock_est[bock_est$SNR >= 1, ]

5 Keep the highest SNR per recording

## ============================================================
## CAP SELECTIONS PER RECORDING (keep up to 5 highest-SNR songs)
## ============================================================
## A recording with many annotated songs would otherwise
## dominate that individual's mean PC/SNR values relative to a
## recording with only one or two -- capping at the 5 best-SNR
## selections keeps every recording's contribution to the
## per-individual average roughly comparable, while still
## preferring its cleanest songs over noisier ones from the same
## recording. Computed on plain vectors (ave/rank), then applied
## as a single row-index subset -- not split()/rbind() -- since
## this is the indexing form already established to preserve
## extended_selection_table attributes correctly (see
## fix_extended_selection_table() in the PCA chunk).

snr_rank_in_recording <- ave(-bock_est$SNR, bock_est$orig.sound.files,
    FUN = function(x) rank(x, ties.method = "first"))

n_recordings_capped <- sum(table(bock_est$orig.sound.files) > 5)
n_selections_before <- nrow(bock_est)

bock_est <- bock_est[snr_rank_in_recording <= 5, ]

cat("Recordings with >5 songs (capped):", n_recordings_capped, "\n")

Recordings with >5 songs (capped): 52 

cat("Selections before cap:", n_selections_before, "| after cap:",
    nrow(bock_est), "\n")

Selections before cap: 773 | after cap: 461 

6 Remove possible duplicated recordings from the same individual

Recordings were screened for likely duplicates of the same individual before principal component analysis:

• Recordings made within 1 km and within 7 days of one another were flagged as possible duplicates (for example, a territorial male re-recorded on a later visit to the same lek or perch), to avoid treating repeat recordings of one bird as independent individuals in the pairwise dataset.
• Distances and date separations were computed independently within each of the four populations rather than across the full dataset, since two recordings from different allopatric populations can never represent the same individual.
• Recordings linked by chains of pairwise matches (e.g. three recordings in close succession at one site, where the first and third do not themselves meet the criteria but each links through the second) were grouped as a single set of potential duplicates. Only the highest-mean-SNR recording was retained from each group; the rest, and their selections, were excluded.
• This is distinct from the pairwise minimum-SNR covariate used later in the statistical model: here SNR chooses among recordings believed to represent the same individual, whereas the covariate accounts for residual recording-quality variation among the individuals that remain.
• This procedure identified 19 groups of two or more recordings (43 of 115 total recordings), reducing the retained recordings from 115 to 91. The two largest groups (four and five recordings) were inspected manually and confirmed to be genuine repeat recordings of the same individual rather than multiple individuals recorded in close proximity.

## Recordings within km_threshold km and week_threshold days of
## each other are treated as repeat recordings of the same bird
## (e.g. a territorial male re-recorded on a later visit to the
## same lek/perch), rather than independent individuals.
## Distances are computed WITHIN each population separately,
## never across populations -- both cheaper than one large
## cross-population distance matrix, and correct, since two
## recordings from different allopatric populations can never be
## the same individual regardless of any coordinate arithmetic.
## Chains of recordings that consecutively meet the criteria
## (A-B, B-C) are grouped together even if A and C do not
## themselves meet the criteria, as long as each links through
## an intermediate recording.  Only the highest-SNR recording is
## kept from each group.

km_threshold <- 1  # km
week_threshold <- 7  # days


## one row per recording, with its location, date, population,
## and mean SNR across its selections -- used only to decide
## which RECORDINGS to keep
rec_meta <- bock_est[!duplicated(bock_est$orig.sound.files), c("orig.sound.files",
    "population", "latitude", "longitude", "date")]
rec_meta$full_date <- as.Date(rec_meta[["date"]])

rec_snr <- aggregate(bock_est$SNR, by = list(orig.sound.files = bock_est$orig.sound.files),
    FUN = mean)
names(rec_snr)[2] <- "SNR"

rec_meta <- merge(rec_meta, rec_snr, by = "orig.sound.files")

## connected components via a small base-R union-find: chains of
## recordings that consecutively meet the criteria collapse into
## one group even if the first and last don't meet the criteria
## directly.
union_find_components <- function(n, edges) {
    parent <- seq_len(n)

    find <- function(x, parent) {
        while (parent[x] != x) x <- parent[x]
        x
    }

    for (i in seq_len(nrow(edges))) {
        ra <- find(edges[i, 1], parent)
        rb <- find(edges[i, 2], parent)
        if (ra != rb)
            parent[ra] <- rb
    }

    sapply(seq_len(n), find, parent = parent)
}

## ---- loop over population: pairwise distances computed within
## each ---- population only, never across populations

rec_meta$dup_group <- NA_character_

for (pop in unique(rec_meta$population)) {

    pop_idx <- which(rec_meta$population == pop)
    pop_sub <- rec_meta[pop_idx, ]
    n_pop <- nrow(pop_sub)

    if (n_pop > 1) {
        coords_pop <- as.matrix(pop_sub[, c("longitude", "latitude")])
        geo_km_pop <- geosphere::distm(coords_pop, fun = distHaversine)/1000
        days_pop <- as.matrix(dist(as.numeric(pop_sub$full_date)))

        idx_dup_pop <- which(geo_km_pop <= km_threshold & days_pop <=
            week_threshold & upper.tri(geo_km_pop), arr.ind = TRUE)

        local_group <- union_find_components(n_pop, idx_dup_pop)
    } else {
        local_group <- 1L
    }

    ## prefix with population name so group IDs are globally
    ## unique once every population's results are combined back
    ## into one table
    rec_meta$dup_group[pop_idx] <- paste(pop, local_group, sep = "_")
}

## sanity check before collapsing: how many recordings are
## affected, and does any group span an implausibly long date
## range (a sign of chaining rather than genuine same-individual
## repeats)?
group_sizes <- table(rec_meta$dup_group)
cat("Groups with >1 recording:", sum(group_sizes > 1), "| recordings affected:",
    sum(group_sizes[group_sizes > 1]), "\n")

Groups with >1 recording: 20 | recordings affected: 45 

group_span <- aggregate(full_date ~ dup_group, data = rec_meta, FUN = function(x) as.numeric(diff(range(x))))
print(group_span[group_span$full_date > week_threshold, ])  # inspect any of these

[1] dup_group full_date
<0 rows> (or 0-length row.names)

## select the highest-SNR recording from each pool of potential
## duplicates
rec_meta <- rec_meta[order(-rec_meta$SNR), ]
recordings_to_keep <- rec_meta$orig.sound.files[!duplicated(rec_meta$dup_group)]

n_before_dedup <- length(unique(bock_est$orig.sound.files))
bock_est <- bock_est[bock_est$orig.sound.files %in% recordings_to_keep,
    ]

attr(bock_est, "metadata") <- bock_metadata[bock_metadata$downloaded_file_name %in%
    bock_est$orig.sound.files, ]

cat("Recordings before dedup:", n_before_dedup, "| after dedup:",
    length(unique(bock_est$orig.sound.files)), "\n")

Recordings before dedup: 120 | after dedup: 95 

7 Map locations

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

8 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, its bounds taken as the lowest and highest frequency limits across all annotation boxes in the dataset, then reduced with a principal component analysis (PCA) on the standardized statistics.
• MFCCs were used because they summarize the overall spectral envelope rather than tracking individual frequency components, making them well suited to the harmonically rich structure of “bock” songs, in which energy is distributed across multiple harmonics.
• Pairwise Euclidean distance between individuals in this reduced space is the acoustic-distance response used in the model below, so the number of components retained determines what the response variable is, not just how it is summarized.
• Because PCA is an orthonormal rotation, truncating to fewer components is 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.
• Components were retained using the broken-stick model: each component’s observed proportion of variance is compared against the proportion it would be expected to carry if the total variance were divided at random among all components, and a component is retained only until the first one that fails to clear that expectation.
• This was the best-performing stopping rule in Jackson’s (1993, Ecology 74:2204–2214) simulation comparison across ecological datasets, and it is the conservative choice for this analysis: acoustic dissimilarity is the model’s response variable rather than a predictor, so every additional noise component retained widens posterior intervals and reduces power to detect genuine geographic or temporal structure downstream.

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

## band taken from the annotations themselves, on both ends, so
## it matches the band already used for sig2noise() upstream
mfcc_tp <- mfcc_stats(bock_est, bp = c(min(bock_est$bottom.freq),
    max(bock_est$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
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.
bock_est2 <- bock_est
bock_est <- cbind(bock_est, pc_scores)
bock_est <- fix_extended_selection_table(X = bock_est, Y = bock_est2)

## 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(bock_est = bock_est, pcs = pcs, n_pcs = n_pcs, cum_var = cum_var),
    "./data/processed/pca_bock_broken_stick.rds")

pca_out <- readRDS("./data/processed/pca_bock_broken_stick.rds")

• The first 15 principal components were retained, together explaining 58.1% of the variance in the standardized MFCC statistics.
• These components are used to compute pairwise Euclidean distances between individuals, the acoustic-distance response variable used in the model below.

9 Check distances

• Before fitting the model, the pairwise predictors are inspected directly: the shape of the geographic-distance distribution (including how 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.

bock_est <- pca_out$bock_est
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 ====

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

PCs retained: 15 | cumulative variance explained: 58.1 %

## SNR was already measured earlier in the pipeline -- bock_est
## already carries an SNR column at this point (see the
## SNR-filtering step upstream, before PCA), so it is not
## recomputed here.

## build_dist_long() (defined once, in Custom functions)
## encapsulates aggregation by individual, all pairwise
## distances, transformations, 2-SD scaling, and NA filtering --
## see its definition for the full step-by-step logic.

dist_long <- build_dist_long(bock_est, pcs, out_file = "./data/processed/dist_long_bock.rds")

Individuals going into this build: 95 

             Brazil Trinidad and Tobago         Venezuela N        Venezuela SE 
                 40                  34                  15                   6 
geo_const (km): 530 | zero-distance pairs: 182 ( 4.1 % )
Final pairs in this build: 4465 

9.1 Data description

dist_long <- readRDS("./data/processed/dist_long_bock.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()

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: 4465 | individuals: 95 

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

geo_distance == 0: 182 ( 4.1 % )

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

log_geo range: 6.27 8.5 

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: 1281 ( 28.7 % )

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

cor(log_geo, same_population): -0.77 

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.02 

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

n_pairs <- nrow(dist_long)
n_indiv <- length(unique(c(as.character(dist_long$individual1), as.character(dist_long$individual2))))

• The pairwise dataset comprised 4465 comparisons among 95 individuals.

9.2 Check for collinearity

Collinearity among all predictors

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.
• Interpreted as: below 5, low concern; 5 to 10, moderate concern; above 10, high concern.
• All VIFs fall below 5, so collinearity is not a major concern for the predictors in this model.

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)

## agg_indiv is local to build_dist_long() and not available
## here, so a lightweight per-individual aggregation is
## recomputed directly from bock_est for this check alone.

agg_snr_check <- aggregate(bock_est$SNR, by = list(orig.sound.files = bock_est$orig.sound.files),
    FUN = mean)
names(agg_snr_check)[2] <- "SNR"

agg_year_check <- bock_est[!duplicated(bock_est$orig.sound.files),
    c("orig.sound.files", "date")]
agg_year_check$year <- as.integer(format(as.Date(agg_year_check$date),
    "%Y"))

snr_year_check <- merge(agg_snr_check, agg_year_check, by = "orig.sound.files")
snr_year_cor <- cor(snr_year_check$SNR, snr_year_check$year, use = "pairwise.complete.obs")

• The correlation between signal-to-noise ratio (SNR) and recording year was also checked, to rule out systematic drift in recording quality over time that could confound the analysis of temporal effects on acoustic divergence.
• The correlation was low (-0.34); SNR does not vary systematically with recording year.

10 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 was summarized with a PCA on standardized MFCC-derived acoustic parameters (see Principal component analysis above for the retention rule). 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.

Predictor construction

• Geographic distance was calculated as the great-circle (haversine) distance between individuals, rather than treating latitude/longitude as flat coordinates: because these populations span a wide latitude range, a flat-coordinate calculation would systematically distort distances in a way that correlates with the population layout itself.
• Geographic distance was log-transformed prior to standardization, correcting for the right-skewed distribution typical of pairwise spatial data (many nearby pairs, few very distant pairs). 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).
• Temporal separation was calculated as the absolute difference in recording date, expressed in fractional years.
• All continuous predictors (geographic distance, temporal separation, minimum SNR) were standardized by dividing by two, rather than one, standard deviation (Gelman 2008). This puts each predictor’s coefficient on the same footing as the binary same_population predictor, which was left as an untransformed 0/1 indicator: under this convention, a coefficient represents the estimated effect of a typical full swing in that predictor, matching what a coefficient on a two-category binary predictor already represents, rather than the effect of a half-swing under the more common one-SD standardization, allowing effect sizes to be compared meaningfully across the continuous and binary predictors in the model.
• 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. 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.

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 separation}_{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 separation}_{ij}\) is the standardized difference in recording date 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 distance and temporal separation was included to model the possibility that vocal divergence accumulates 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, 2-SD 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, scaled to the standardized response and predictors: 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 is z-scored and each continuous predictor is scaled to represent a two-SD swing (see above), 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 11,000 iterations per chain (1,000 warmup).
• Convergence was assessed with the potential scale reduction factor and effective sample sizes; all parameters reached R-hat = 1.00 with bulk and tail effective sample sizes above 12,000.

10.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.

pca_out <- readRDS("./data/processed/pca_bock_broken_stick.rds")
n_pcs <- pca_out$n_pcs
cum_var <- pca_out$cum_var

dist_long <- readRDS("./data/processed/dist_long_bock.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_bock 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 the standardized response and predictors
## (response z-scored, SD = 1; continuous predictors scaled to a
## 2-SD swing, see Statistical analysis above), 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,
    "_bock")

iterations <- 11000
warmup <- 1000
## ============================================================
## 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 = iterations, warmup = warmup, 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 = iterations, warmup = warmup,
    backend = "cmdstanr", threads = threading(8), control = list(adapt_delta = 0.95,
        max_treedepth = 15), file = fit_path("mod_interaction"), file_refit = "on_change")

results_bock <- list(n_pcs = n_pcs, 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_bock, "./data/processed/results_bock.rds")

10.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, so a Gaussian prior (and posterior) predictive distribution will 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.

results_bock <- readRDS("./data/processed/results_bock.rds")

## 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(results_bock$prior_check_model, ndraws = 100) + labs(title = "Prior predictive check")

10.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.

mod_interaction <- results_bock$model
dist_long <- results_bock$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")

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

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

10.1.3 Results and biological interpretation

summ <- extended_summary(mod_interaction, highlight = FALSE, 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.189	-0.269	-0.110	1	12149.33	19646.20
Same population	-0.759	-0.856	-0.663	1	15698.53	21897.35
Temporal separation	0.046	-0.014	0.106	1	18783.25	23938.73
Geographic distance	0.038	-0.054	0.131	1	15783.74	20975.14
Temporal separation × geographic distance	-0.229	-0.296	-0.162	1	25912.65	26334.40

describe_posterior(mod_interaction, test = c("pd", "rope", "bf"),
    rope_range = c(-0.1, 0.1))

	Parameter	Median	CI	CI_low	CI_high	pd	ROPE_CI	ROPE_low	ROPE_high	ROPE_Percentage	log_BF	Rhat	ESS_tail
2	b_Intercept	0.0838196	0.95	-0.6818204	0.8306799	0.599175	0.95	-0.1	0.1	0.2352368	-1.070789	1.000885	12396
3	b_min_snr_sc	-0.1892340	0.95	-0.2686283	-0.1103955	1.000000	0.95	-0.1	0.1	0.0000000	8.342134	1.000179	19646
4	b_same_population1	-0.7594415	0.95	-0.8562604	-0.6631937	1.000000	0.95	-0.1	0.1	0.0000000	35.248516	1.000169	21897
5	b_time_separation_sc	0.0464148	0.95	-0.0138683	0.1064004	0.934400	0.95	-0.1	0.1	0.9841842	-1.651079	1.000154	23939
1	b_geo_distance_sc	0.0381715	0.95	-0.0543412	0.1305531	0.788900	0.95	-0.1	0.1	0.9230000	-2.039086	1.000161	20975
6	b_time_separation_sc:geo_distance_sc	-0.2291295	0.95	-0.2957460	-0.1617190	1.000000	0.95	-0.1	0.1	0.0000000	10.678659	1.000102	26334

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

cat("Probability of direction")

Probability of direction

p_direction(mod_interaction, effects = "fixed")

	Parameter	pd	Effects	Component
2	b_Intercept	0.599175	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.934400	fixed	conditional
1	b_geo_distance_sc	0.788900	fixed	conditional
6	b_time_separation_sc:geo_distance_sc	1.000000	fixed	conditional

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), inherit.aes = FALSE, nudge_y = 0.45,
        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)

gg_estimates

ggsave(filename = "./output/fig_results_bock.png", plot = gg_estimates,
    device = grDevices::png, width = 6, height = 5, units = "in",
    dpi = 300)

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

results_bock <- readRDS("./data/processed/results_bock.rds")
mod_interaction <- results_bock$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
## negative -------

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

## ---- where does the temporal-separation slope cross zero?
## ---------------- Determines whether the negative interaction
## means the time effect merely flattens as geographic distance
## grows, or actually reverses sign within the range of
## geographic distances actually observed. The wording in the
## interpretation below branches on this rather than asserting
## either.

fx_est <- fixef(mod_interaction)[, "Estimate"]
cross <- unname(-fx_est["time_separation_sc"]/fx_est["time_separation_sc:geo_distance_sc"])
geo_rng <- range(dist_long$geo_distance_sc, na.rm = TRUE)

int_reverses <- cross > geo_rng[1] && cross < geo_rng[2]

cat("Temporal-separation slope crosses zero at geo_distance_sc =",
    round(cross, 2), "\nObserved geo_distance_sc range:", round(geo_rng,
        2), "\nReverses within observed range:", int_reverses, "\n")

Temporal-separation slope crosses zero at geo_distance_sc = 0.2 
Observed geo_distance_sc range: -0.85 0.7 
Reverses within observed range: TRUE 

Overall patterns

• Recording quality (minimum pairwise SNR) had a credible association with acoustic dissimilarity (β = -0.189, 95% CrI [-0.269, -0.11]).
• Belonging to the same population was associated with decreased acoustic dissimilarity relative to belonging to different populations (β = -0.759, 95% CrI [-0.856, -0.663]), a credible effect.
• Geographic distance was associated with increased acoustic dissimilarity (β = 0.038, 95% CrI [-0.054, 0.131]), though the credible interval included zero.
• Temporal separation was associated with increased acoustic dissimilarity (β = 0.046, 95% CrI [-0.014, 0.106]), though the credible interval included zero.
• The temporal separation × geographic distance interaction showed a credible effect (β = -0.229, 95% CrI [-0.296, -0.162]), with a posterior probability of 100% that the coefficient is negative.

TipOverall biological interpretation

• Acoustic divergence in Bearded Bellbird “bock” songs is structured primarily by population identity, with geography and time contributing jointly rather than as independent main effects.
• Population identity showed the strongest effect of any predictor, indicating that variation in “bock” songs is structured among discrete populations rather than as a smooth function of distance.
• Geographic distance alone did not show a credible association with acoustic dissimilarity once the other predictors were accounted for.
• Temporal separation alone did not show a credible association with acoustic dissimilarity once the other predictors were accounted for.
• The temporal-separation-by-geographic-distance interaction showed a credible effect, indicating that the rate at which songs diverge over time depends on how far apart the compared recordings were made.
• The negative coefficient means that temporal divergence is strongest among geographically proximate recordings, and within the observed range of geographic distance the relationship reverses, with widely separated recordings showing no further increase, or a slight decrease, in dissimilarity with elapsed time. A pattern of this kind, detectable over only a few decades, is itself notable: comparable genetically-driven phenotypic divergence typically requires many more generations, absent unusually strong selection, whereas culturally transmitted traits such as song can shift within only a few generations of vocal tutoring, making cultural evolution mediated by vocal production learning a plausible, though unconfirmed, explanation. Vocal production learning has been documented behaviorally in the congeneric Three-wattled Bellbird (Procnias tricarunculatus), further supported by evidence that its geographic song variation is decoupled from genetic population structure (Saranathan et al. 2007; Kroodsma et al. 2013), and Snow (1970) reported apparent temporal change in Bearded Bellbird song in Trinidad; together, these findings raise, without resolving, the possibility that a comparable process occurs in P. averano.
• Recording quality had a detectable influence on measured acoustic dissimilarity, confirming that controlling for it was a necessary methodological step. This is a data-quality observation rather than a biological one.
• This case study illustrates how recordings assembled programmatically with suwo can be carried through a complete geographic-variation analysis. Beyond simplifying acquisition, this is also what made the analyses above possible in the first place: harmonized location, date, and recordist metadata across repositories is what permitted the geographic-population assignments and pairwise temporal-separation calculations used throughout, inferences that would be difficult to draw reliably from any single repository’s metadata conventions alone, let alone from several reconciled by hand.

---

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. (2017). warbleR: an r package to streamline analysis of animal acoustic signals. Methods in Ecology & Evolution, 8(2), 184-191.

Session information

Click to see

─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.5.2 (2025-10-31)
 os       Ubuntu 22.04.4 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-09-04
 pandoc   3.6.3 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/x86_64/ (via rmarkdown)
 quarto   1.8.25 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/quarto

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

 [1] /home/m/R/x86_64-pc-linux-gnu-library/4.5
 [2] /usr/local/lib/R/site-library
 [3] /usr/lib/R/site-library
 [4] /usr/lib/R/library
 * ── Packages attached to the search path.

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