Cape parrot vocal dialects

Author

Marcelo Araya-Salas

Published

August 26, 2026

Code 
# set working directory as project directory or one directory above,
knitr::opts_knit$set(root.dir = "..")

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

 

1 Purpose

  • Measure acoustic structure of cape parrot contact calls

  • Compare acoustic dissimilarity between individuals from different localities and regions

 

Load packages

Code 
# knitr is require for creating html/pdf/word reports formatR is
# used for soft-wrapping code

# install/ load packages
sketchy::load_packages(packages = c("knitr", "formatR", "viridis",
    "warbleR", github = "maRce10/PhenotypeSpace", "ggplot2", "randomForest",
    "mlbench", "caret", "pbapply", "vegan", "umap", "brms", "brmsish",
    "maRce10/ohun", "geosphere", "performance", "ggcorrplot", "emmeans",
    "patchwork"))
Warning: replacing previous import 'brms::rstudent_t' by 'ggdist::rstudent_t'
when loading 'brmsish'
Warning: replacing previous import 'brms::dstudent_t' by 'ggdist::dstudent_t'
when loading 'brmsish'
Warning: replacing previous import 'brms::qstudent_t' by 'ggdist::qstudent_t'
when loading 'brmsish'
Warning: replacing previous import 'brms::pstudent_t' by 'ggdist::pstudent_t'
when loading 'brmsish'

2 Acoustic analysis

Code 
## Format data
dat <- read.csv("./data/raw/consolidated_sound_files_CPV_contact_calls_USEaug2026 - UPDATED_USE for analyses.csv")

warbleR_options(path = "./data/raw/consolidated_files")
st <- selection_table(whole.recs = TRUE)

st$region <- sapply(st$sound.files, function(x) dat$Regions..4.[dat$New_Name ==
    x][1])


st <- st[complete.cases(st$region), ]

# make a folder for each region and save the spectrograms in the
# corresponding folder

for (i in unique(st$region)) {
    dir.create(paste0("./data/processed/spectros_by_region/", i),
        showWarnings = FALSE)
    spectrograms(st[st$region == i, ], wl = 512, flim = c(0, 10),
        dest.path = paste0("./data/processed/spectros_by_region/",
            i), pal = viridis, collevels = seq(-100, 0, 5))
}

# st$sorted <- sapply(st$sound.files, function(x)
# dat$Sorted[dat$New_Name == x][1]) table(st$sorted)

# spectrograms(st, wl = 512, flim = c(0, 10), dest.path =
# './data/processed/spectrograms', pal = viridis, collevels =
# seq(-100, 0, 5)) spectrograms(st[st$sorted == 'unsorted', ],
# wl = 512, flim = c(0, 10), dest.path =
# './data/processed/unsorted_spectrograms', pal = viridis,
# collevels = seq(-100, 0, 5)) tailor_sels(st, auto.next = TRUE,
# flim = c(0, 8), collevels = seq(-100, 0, 5))

2.1 Make selection table

Code 
sel_tab <- selection_table(path = "./data/raw/consolidated_files/",
    whole.recs = TRUE)

tailored <- read.csv("./data/raw/consolidated_files/seltailor_output.csv")

tailored <- tailored[tailored$tailored == "y", ]


non_tailored <- sel_tab[!sel_tab$sound.files %in% tailored$sound.files,
    ]
non_tailored$tailored <- "n"


tailored$top.freq[is.na(tailored$bottom.freq)] <- non_tailored$bottom.freq <- min(tailored$bottom.freq,
    na.rm = TRUE)
tailored$top.freq[is.na(tailored$top.freq)] <- non_tailored$top.freq <- max(tailored$top.freq,
    na.rm = TRUE)

comm_names <- intersect(names(tailored), names(non_tailored))

all_sels <- rbind(tailored[, comm_names], non_tailored[, comm_names])

write.csv(all_sels, "./data/processed/selection_table_entire_sound_files.csv",
    row.names = FALSE)

2.2 Run cross-correlation

Code 
sel_tab <- read.csv("./data/processed/selection_table_entire_sound_files.csv")

xcorr <- cross_correlation(X = sel_tab, path = "./data/raw/consolidated_files/",
    method = 2, parallel = 1)

rownames(xcorr) <- gsub("-1$", "", rownames(xcorr))

colnames(xcorr) <- gsub("-1$", "", colnames(xcorr))

saveRDS(xcorr, "./data/processed/cross_correlation_matrix.RDS")

# less than 0.1% were undefined
sum(is.infinite(xcorr))/length(xcorr)

# convert infinite to mean xcorr
xcorr[is.infinite(xcorr)] <- mean(xcorr[!is.infinite(xcorr) & xcorr <
    1])

xcorr_mds <- cmdscale(d = as.dist(xcorr), k = 2)

rownames(xcorr_mds) <- gsub("-1$", "", rownames(xcorr_mds))

saveRDS(xcorr_mds, "./data/processed/cross_correlation_MDS.RDS")

# measure acoustic features
sel_tab$bottom.freq[is.na(sel_tab$bottom.freq)] <- min(sel_tab$bottom.freq,
    na.rm = TRUE)

sel_tab$top.freq <- max(sel_tab$top.freq)

acous_feat <- spectro_analysis(X = sel_tab, path = "./data/raw/consolidated_files/",
    parallel = 30)

acous_feat$region <- sapply(acous_feat$sound.files, function(x) dat$Regions4[dat$New_Name ==
    x][1])

acous_feat$recording <- substr(acous_feat$sound.files, 0, 4)

saveRDS(acous_feat, "./data/processed/acoustic_features.RDS")

3 Data description

Code 
# add data from second location
dat <- read.csv("./data/raw/consolidated_sound_files_CPV_contact_calls_USEaug2026 - UPDATED_USE for analyses.csv")


names(dat)[grep("Regions..4.", names(dat))] <- "Regions4"

names(dat) <- gsub("..cluster.", ".for.cluster", names(dat))

dat <- dat[!is.na(dat$Location.for.cluster) & !is.na(dat$Longitude.for.cluster) &
    !is.na(dat$Latitude.for.cluster), ]
  • 6287 calls
  • 14 localities
  • 4 regions
  • Number of localities per region:
Code 
agg <- aggregate(Location.for.cluster ~ Regions4, dat, function(x) length(unique(x)))

names(agg) <- c("region", "localities")

agg$calls <- aggregate(Location.for.cluster ~ Regions4, dat, length)[,
    2]

agg$localities <- aggregate(Location.for.cluster ~ Regions4, dat,
    function(x) paste(unique(x), collapse = "-"))[, 2]

agg
region localities calls
central subA Jetty River Lodge-Ntafufu Ecolodge 687
central subB Polela Sawmill-iGxalingenwa Nature Reserve-Marutswa Forest-Salt Spring Farm-Hoha Forest 1082
northern Amorentia 717
southern Hogsback-Schwarzwald Forest-Alice Pecan Orchard-Stutterheim-King William’s Town-Baddaford Farm 3801
  • Region of each locality:
Code 
agg <- aggregate(Regions4 ~ Location.for.cluster, dat, function(x) paste(unique(x),
    collapse = "-"))

names(agg) <- c("localities", "region")

agg[order(agg$region), 2:1]
region localities
7 central subA Jetty River Lodge
10 central subA Ntafufu Ecolodge
5 central subB Hoha Forest
6 central subB iGxalingenwa Nature Reserve
9 central subB Marutswa Forest
11 central subB Polela Sawmill
12 central subB Salt Spring Farm
2 northern Amorentia
1 southern Alice Pecan Orchard
3 southern Baddaford Farm
4 southern Hogsback
8 southern King William’s Town
13 southern Schwarzwald Forest
14 southern Stutterheim

4 Statistical analysis

4.1 Geographical levels of variation

To evaluate whether Cape Parrot contact call structure is better explained by region, locality, geographic distance, or some combination of these, we fitted a set of competing Bayesian mixed-effects models on pairwise acoustic dissimilarities between recordings, accounting for the non-independence inherent to pairwise distance data.

Model specifications

Six competing models were specified, each sharing the same multi-membership random-effect structure but differing in which fixed effect(s) they include, to test via leave-one-out cross-validation (LOO) whether region, locality, geographic distance, or some combination best explains acoustic dissimilarity between recordings:

\[ \text{mod\_region}: \quad \text{acoustic dissimilarity}_{ij} \sim \text{same region}_{ij} + (1 \mid \text{mm(recording}_i, \text{recording}_j)) \]

\[ \text{mod\_locality}: \quad \text{acoustic dissimilarity}_{ij} \sim \text{same locality}_{ij} + (1 \mid \text{mm(recording}_i, \text{recording}_j)) \]

\[ \text{mod\_geodist}: \quad \text{acoustic dissimilarity}_{ij} \sim \text{geographic distance}_{ij} + (1 \mid \text{mm(recording}_i, \text{recording}_j)) \]

\[ \text{mod\_region\_geo}: \quad \text{acoustic dissimilarity}_{ij} \sim \text{same region}_{ij} \times \text{geographic distance}_{ij} + (1 \mid \text{mm(recording}_i, \text{recording}_j)) \]

\[ \text{mod\_both}: \quad \text{acoustic dissimilarity}_{ij} \sim \text{same region}_{ij} + \text{same locality}_{ij} + (1 \mid \text{mm(recording}_i, \text{recording}_j)) \]

\[ \text{mod\_all}: \quad \text{acoustic dissimilarity}_{ij} \sim \text{same region}_{ij} + \text{same locality}_{ij} + \text{geographic distance}_{ij} + (1 \mid \text{mm(recording}_i, \text{recording}_j)) \]

where:

  • \(\text{acoustic dissimilarity}_{ij}\) is the mean pairwise cross-correlation-based dissimilarity (\(1 - r\)) between all calls of recording \(i\) and all calls of recording \(j\).

  • \(\text{same region}_{ij}\) and \(\text{same locality}_{ij}\) are binary predictors indicating whether recordings \(i\) and \(j\) come from the same region/locality (1) or different ones (0).

  • \(\text{geographic distance}_{ij}\) is the great-circle (haversine) distance between recordings \(i\) and \(j\), log-transformed (25th-percentile offset added before log to avoid zero-distance pairs at the same locality) and standardized by 2 SD (Gelman 2008), so its coefficient is on a comparable “typical full swing” scale to the binary predictors.

  • mm(recording\(_i\), recording\(_j\)) is a multi-membership random intercept accounting for the repeated use of the same recordings across pairwise comparisons.

To isolate whether geographic distance matters specifically within regions - as opposed to only driving the separation between regions - a seventh model was fit on the subset of same-region pairs only:

\[ \text{mod\_geodist\_within\_region}: \quad \text{acoustic dissimilarity}_{ij} \sim \text{geographic distance}_{ij} + (1 \mid \text{mm(recording}_i, \text{recording}_j)), \quad \text{same region}_{ij} = \text{same} \]

This avoids the extrapolation problem inherent to reading mod_region_geo‘s main effect: same-region and different-region pairs occupy almost non-overlapping ranges of geographic distance (regions are, by construction, spatially clustered), so any single-point comparison between the two groups’ fitted lines extrapolates at least one of them well outside its supported data range. A model fit on same-region pairs alone estimates the within-region distance slope directly, with no such extrapolation.

  • The response variable was mean pairwise acoustic dissimilarity (1 - cross-correlation) between recordings, rescaled to the open (0, 1) interval and modeled with a Beta error distribution.
  • The unit of analysis is the recording pair, not the individual call pair: region, locality, and geographic distance are properties of the recording, and acoustic dissimilarity was averaged across all call pairs belonging to each pair of recordings.
  • Recording identity was modeled as a multi-membership random intercept, (1 | mm(rec1, rec2)), because each recording contributes to multiple pairwise comparisons and its pairs are not independent.
  • Predictor collinearity was checked (correlation matrix and VIF) before interpreting the multi-predictor models; see Check predictor collinearity below.
  • Mildly regularizing priors were used: Normal(0, 1) for the intercept, Normal(0, 0.5) for fixed effects, Exponential(2) for the random-effect SD, and Exponential(1) for the Beta precision parameter (phi).
  • Models were 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.
Code 
xcorr <- readRDS("./data/processed/cross_correlation_matrix.RDS")

dat$recording <- substr(dat$Old_Name, 0, 4)

## ---------------------------------------------------------------
## Region vs. locality as predictors of acoustic (dis)similarity
## brms multi-membership model, recording as the mm() grouping
## factor every pairwise combination of the 116 recordings,
## filled with mean acoustic distance across their calls
## ---------------------------------------------------------------
## 2. similarity -> dissimilarity (full matrix, keep symmetric,
## no NAs)
## -----------------------------------------------------------------------

dissim <- 1 - xcorr

## -----------------------------------------------------------------------
## 3. all pairwise combinations of the 116 recordings, filled
## with the mean call-level dissimilarity between the two
## recordings
## -----------------------------------------------------------------------

rec_ids <- dat$recording  # recording id per call, same order as dissim


recordings <- sort(unique(rec_ids))  # 116 unique recordings
n_rec <- length(recordings)


pairs_idx <- t(combn(n_rec, 2))  # 6670 combinations, 116 choose 2

dist_dat <- data.frame(recording.1 = recordings[pairs_idx[, 1]], recording.2 = recordings[pairs_idx[,
    2]], stringsAsFactors = FALSE)



dist_dat$mean_dissim <- sapply(seq_len(nrow(pairs_idx)), function(p) {
    f1 <- dat$New_Name[dat$recording == dist_dat$recording.1[p]]
    f2 <- dat$New_Name[dat$recording == dist_dat$recording.2[p]]

    dists <- as.vector(dissim[rownames(dissim) %in% f1, colnames(dissim) %in%
        f2])

    dists <- dists[!is.infinite(dists)]

    mean(dists, na.rm = TRUE)
})

## -----------------------------------------------------------------------
## 4. attach one locality / region per recording (assumes a
## recording was made at a single site - check this holds in
## your data)
## -----------------------------------------------------------------------

dist_dat$locality.1 <- sapply(dist_dat$recording.1, function(x) unique(dat$Location.for.cluster[dat$recording ==
    x])[1])

dist_dat$locality.2 <- sapply(dist_dat$recording.2, function(x) unique(dat$Location.for.cluster[dat$recording ==
    x])[1])

dist_dat$region.1 <- sapply(dist_dat$recording.1, function(x) unique(dat$Regions4[dat$recording ==
    x])[1])

dist_dat$region.2 <- sapply(dist_dat$recording.2, function(x) unique(dat$Regions4[dat$recording ==
    x])[1])


# measure geographic distance
dist_dat$geo_distance <- sapply(seq_len(nrow(dist_dat)), function(x) {
    coords1 <- dat[dat$recording == dist_dat$recording.1[x], c("Longitude",
        "Latitude")][1, ]
    coords2 <- dat[dat$recording == dist_dat$recording.2[x], c("Longitude",
        "Latitude")][1, ]
    coords <- as.matrix(rbind(coords1, coords2))
    dist_geo_mat <- distm(coords, fun = distHaversine)/1000
    return(dist_geo_mat[1, 2])
})


geo_const <- unname(quantile(dist_dat$geo_distance, 0.25, na.rm = TRUE))
dist_dat$log_geo_distance <- log(dist_dat$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_dat$geo_distance_sc <- scale(dist_dat$log_geo_distance)[, 1]/2


dist_dat$same_locality <- factor(ifelse(dist_dat$locality.1 == dist_dat$locality.2,
    "same", "different"), levels = c("different", "same"))

dist_dat$same_region <- factor(ifelse(dist_dat$region.1 == dist_dat$region.2,
    "same", "different"), levels = c("different", "same"))

names(dist_dat)[names(dist_dat) %in% c("recording.1", "recording.2")] <- c("rec1",
    "rec2")
# dist_dat$rec1 <- factor(dist_dat$rec1) dist_dat$rec2 <-
# factor(dist_dat$rec2)

## -----------------------------------------------------------------------
## 5. dissimilarity to open (0,1) interval for beta family
## -----------------------------------------------------------------------

rescale01 <- function(x, eps = 1e-04) (x * (1 - 2 * eps)) + eps
dist_dat$dissim_beta <- rescale01(dist_dat$mean_dissim)

## -----------------------------------------------------------------------
## 5b. check collinearity among predictors before fitting
## multi-predictor models - locality is nested within region and
## geo_distance is partly confounded with both, so this is worth
## checking explicitly
## -----------------------------------------------------------------------

saveRDS(dist_dat, "./data/processed/dist_dat.RDS")

# correlation matrix (binary predictors coded 0/1 for
# correlation purposes)
predictor_cor <- data.frame(same_region_num = as.numeric(dist_dat$same_region ==
    "same"), same_locality_num = as.numeric(dist_dat$same_locality ==
    "same"), geo_distance_sc = dist_dat$geo_distance_sc)

cor_mat <- cor(predictor_cor, use = "complete.obs")
kable(cor_mat, digits = 2)

saveRDS(cor_mat, "./data/processed/predictor_correlation_matrix.RDS")

# VIF (generalized VIF for factors) on the fullest fixed-effect
# specification - run once mod_all is fitted, see below
# performance::check_collinearity(mod_all)

## -----------------------------------------------------------------------
## 6. priors
## -----------------------------------------------------------------------

priors <- c(prior(normal(0, 1), class = "Intercept"), prior(normal(0,
    0.5), class = "b"), prior(exponential(2), class = "sd"), prior(exponential(1),
    class = "phi"))

## -----------------------------------------------------------------------
## 7. competing models: locality vs region, same mm() structure
## -----------------------------------------------------------------------

mod_region <- brm(bf(dissim_beta ~ same_region + (1 | mm(rec1, rec2))),
    data = dist_dat, family = Beta(), prior = priors, chains = 4,
    cores = 4, iter = 4000, warmup = 1000, backend = "cmdstanr", threads = threading(8),
    control = list(adapt_delta = 0.95), seed = 123, file = "./data/processed/brms_model_region")

mod_locality <- brm(bf(dissim_beta ~ same_locality + (1 | mm(rec1,
    rec2))), data = dist_dat, family = Beta(), prior = priors, backend = "cmdstanr",
    threads = threading(8), chains = 4, cores = 4, iter = 4000, warmup = 1000,
    control = list(adapt_delta = 0.95), seed = 123, file = "./data/processed/brms_model_locality")


mod_geodist <- brm(bf(dissim_beta ~ geo_distance_sc + (1 | mm(rec1,
    rec2))), data = dist_dat, family = Beta(), prior = priors, backend = "cmdstanr",
    threads = threading(8), chains = 4, cores = 4, iter = 4000, warmup = 1000,
    control = list(adapt_delta = 0.95), seed = 123, file = "./data/processed/brms_model_geographic_distance")


mod_geodist_region <- brm(bf(dissim_beta ~ same_region * geo_distance_sc +
    (1 | mm(rec1, rec2))), data = dist_dat, family = Beta(), prior = priors,
    backend = "cmdstanr", threads = threading(8), chains = 4, cores = 4,
    iter = 4000, warmup = 1000, control = list(adapt_delta = 0.95),
    seed = 123, file = "./data/processed/brms_model_geographic_distance_region_interaction")


# full model with both terms - check for
# collinearity/confounding between region and locality before
# trusting this one

mod_both <- brm(bf(dissim_beta ~ same_region + same_locality + (1 |
    mm(rec1, rec2))), data = dist_dat, family = Beta(), prior = priors,
    backend = "cmdstanr", threads = threading(8), chains = 4, cores = 4,
    iter = 4000, warmup = 1000, control = list(adapt_delta = 0.95),
    seed = 123, file = "./data/processed/brms_model_both")

mod_all <- brm(bf(dissim_beta ~ same_region + same_locality + geo_distance_sc +
    (1 | mm(rec1, rec2))), data = dist_dat, family = Beta(), prior = priors,
    backend = "cmdstanr", threads = threading(8), chains = 4, cores = 4,
    iter = 4000, warmup = 1000, control = list(adapt_delta = 0.95),
    seed = 123, file = "./data/processed/brms_model_all")

# VIF on the full fixed-effect specification - run now that
# mod_all exists
performance::check_collinearity(mod_all)

## -----------------------------------------------------------------------
## 7b. within-region subset model: isolates the
## geographic-distance effect at the scale where it can actually
## be estimated without extrapolating across the near-disjoint
## distance ranges of same-region vs. different-region pairs
## (see caveat in Model specifications above)
## -----------------------------------------------------------------------

dist_dat_same_region <- dist_dat[dist_dat$same_region == "same", ]

mod_geodist_within_region <- brm(bf(dissim_beta ~ geo_distance_sc +
    (1 | mm(rec1, rec2))), data = dist_dat_same_region, family = Beta(),
    prior = priors, backend = "cmdstanr", threads = threading(8),
    chains = 4, cores = 4, iter = 4000, warmup = 1000, control = list(adapt_delta = 0.95),
    seed = 123, file = "./data/processed/brms_model_geographic_distance_within_region")



## -----------------------------------------------------------------------
## 8. model comparison
## -----------------------------------------------------------------------

mod_region <- add_criterion(mod_region, "loo")
mod_locality <- add_criterion(mod_locality, "loo")
mod_both <- add_criterion(mod_both, "loo")
mod_geodist <- add_criterion(mod_geodist, "loo")
mod_geodist_region <- add_criterion(mod_geodist_region, "loo")
mod_all <- add_criterion(mod_all, "loo")

# fit on a different (smaller, same-region-only) dataset, so
# this is NOT comparable via loo_compare() to the models above -
# reported separately
mod_geodist_within_region <- add_criterion(mod_geodist_within_region,
    "loo")

The following table ranks models by out-of-sample predictive fit (elpd_loo), with the top model set to elpd_diff = 0 and all others shown as the gap below it, plus a standard error (se_diff) for that gap.

Code 
mod_region <- readRDS("./data/processed/brms_model_region.rds")
mod_locality <- readRDS("./data/processed/brms_model_locality.rds")
mod_both <- readRDS("./data/processed/brms_model_both.rds")
mod_region_geo <- readRDS("./data/processed/brms_model_geographic_distance_region_interaction.rds")
mod_all <- readRDS("./data/processed/brms_model_all.rds")
mod_geodist <- readRDS("./data/processed/brms_model_geographic_distance.rds")
mod_geodist_within_region <- readRDS("./data/processed/brms_model_geographic_distance_within_region.rds")

loo_comp <- loo_compare(mod_region, mod_locality, mod_both, mod_geodist,
    mod_region_geo, mod_all)

kable(as.data.frame(loo_comp), digits = 2)
elpd_diff se_diff elpd_loo se_elpd_loo p_loo se_p_loo looic se_looic
mod_region_geo 0.00 0.00 11132.90 88.52 114.27 3.39 -22265.80 177.03
mod_all -8.47 6.78 11124.43 88.85 114.67 3.39 -22248.87 177.70
mod_geodist -17.87 6.19 11115.03 88.59 113.57 3.37 -22230.06 177.17
mod_both -25.89 8.23 11107.01 91.98 118.83 3.52 -22214.01 183.95
mod_region -27.08 7.72 11105.82 88.90 113.22 3.36 -22211.65 177.80
mod_locality -77.07 15.11 11055.82 92.22 118.18 3.53 -22111.65 184.44

mod_region_geo (the region x distance interaction) has the best expected out-of-sample fit of all six models - clearly ahead of mod_geodist alone (Δ = -17.87, SE 6.19, ~3 SEs), and far ahead of mod_region, mod_both, and mod_locality. mod_all (all three predictors, no interaction) is statistically indistinguishable from mod_region_geo (Δ = -8.47, SE 6.78, well under 1.5 SEs) but doesn’t isolate the interaction the way mod_region_geo does, so mod_region_geo is the more directly interpretable model for the question of whether distance matters differently within vs. between regions.

4.1.1 Check predictor collinearity

same_locality is nested within same_region (every same-locality pair is automatically a same-region pair) and geo_distance_sc is partly confounded with both (regions are, by construction, spatially clustered), so collinearity was checked before trusting coefficients from the multi-predictor models.

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

pred_vars <- c("same_region", "same_locality", "geo_distance_sc")
pred_df <- dist_dat[, pred_vars]
pred_df$same_region <- as.numeric(pred_df$same_region == "same")
pred_df$same_locality <- as.numeric(pred_df$same_locality == "same")

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

Code 
## Variance inflation factors (VIFs) were computed for the fixed
## effects to check for collinearity. VIFs are interpreted
## conventionally: values below 5 indicate low concern, values
## between 5 and 10 indicate moderate concern, and values above
## 10 indicate high concern. In this case, same_locality (VIF =
## 3.1) is low concern, same_region (VIF = 7.4) is moderate
## concern, and geo_distance_sc (VIF = 11.3) is high concern -
## region, locality, and geographic distance all partly encode
## the same underlying spatial structure, so mod_all's
## individual coefficients should not be read as independent
## effects.

vif_mod <- lm(dissim_beta ~ same_region + same_locality + geo_distance_sc,
    data = dist_dat)
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.
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)

Code 
check_collinearity(mod_all)
Term VIF VIF_CI_low VIF_CI_high SE_factor Tolerance Tolerance_CI_low Tolerance_CI_high
same_region 7.443463 7.124530 7.779005 2.728271 0.1343461 0.1285511 0.1403601
same_locality 3.060321 2.943171 3.184535 1.749377 0.3267631 0.3140176 0.3397696
geo_distance_sc 11.321602 10.824196 11.844192 3.364759 0.0883267 0.0844296 0.0923856

Collinearity results

  • geo_distance_sc shows severe collinearity (VIF = 11.3, tolerance = 0.09) — over 91% of its variance is explained by same_region and same_locality combined.
  • same_region also shows concerning collinearity (VIF = 7.4, tolerance = 0.13).
  • same_locality is the least collinear of the three (VIF = 3.1, tolerance = 0.33) but still above the conventional VIF = 5 caution threshold.
  • This isn’t a data-quality problem — it reflects the actual structure of the sampling design: regions are spatially clustered, so region identity, locality identity, and geographic distance all partly encode the same underlying geographic information.
  • Practical consequence: individual coefficients in mod_all (which includes all three predictors together) can’t be cleanly interpreted as independent, “holding-everything-else-constant” effects — they’re each absorbing overlapping variance. mod_all remains valid for LOO/predictive comparison, but its coefficient table shouldn’t be over-interpreted term-by-term.

4.1.2 Region x distance interaction

Code 
extended_summary(mod_region_geo, highlight = TRUE, print.name = FALSE)
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 0.5) Intercept-normal(0, 1) phi-exponential(1) sd-exponential(2) dissim_beta ~ same_region * geo_distance_sc + (1 | mm(rec1, rec2)) 4000 4 1 1000 0 (0%) 0 194.019 716.161 123
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept -0.028 -0.139 0.086 1.012 194.019 716.161
b_same_regionsame 0.197 0.113 0.280 1.001 1394.714 2366.607
b_geo_distance_sc 0.481 0.345 0.617 1.001 1344.062 2342.787
b_same_regionsame:geo_distance_sc -0.427 -0.566 -0.290 1.001 1425.763 2571.987

The interaction is the key result: geo_distance_sc (0.481, CI 0.345-0.617) is the distance slope for different-region pairs - dissimilarity rises steeply with distance. The interaction term same_regionsame:geo_distance_sc (-0.427, CI -0.566 to -0.290) nearly cancels that slope for same-region pairs, leaving an effective within-region slope close to flat (≈ 0.05). Geographic distance predicts acoustic dissimilarity strongly between regions, and barely at all within them.

Same-region and different-region pairs occupy almost non-overlapping ranges of geo_distance_sc, since regions are spatially clustered:

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

kable(data.frame(same_region = c("different", "same"), min = tapply(dist_dat$geo_distance_sc,
    dist_dat$same_region, min), mean = tapply(dist_dat$geo_distance_sc,
    dist_dat$same_region, mean), max = tapply(dist_dat$geo_distance_sc,
    dist_dat$same_region, max)), digits = 2, row.names = FALSE)
same_region min mean max
different 0.31 0.59 0.87
same -0.66 -0.35 0.11

Because of this near-disjoint range, the model’s main effect for same_region (evaluated at geo_distance_sc = 0, i.e. the overall mean distance) extrapolates outside the real data range for at least one of the two groups, and should not be read as “same-region pairs are more/less dissimilar than different-region pairs” at some typical distance - no such shared typical distance exists in the data. The plot below restricts each group’s fitted line to its own observed distance range to avoid this:

Code 
library(ggplot2)

region_ranges <- tapply(dist_dat$geo_distance_sc, dist_dat$same_region,
    range)

newdat <- do.call(rbind, lapply(names(region_ranges), function(r) {
    data.frame(same_region = r, geo_distance_sc = seq(region_ranges[[r]][1],
        region_ranges[[r]][2], length.out = 100))
}))

preds <- fitted(mod_region_geo, newdata = newdat, re_formula = NA)
newdat <- cbind(newdat, preds)

ggplot(newdat, aes(x = geo_distance_sc, y = Estimate, color = same_region,
    fill = same_region)) + geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5),
    alpha = 0.3, color = NA) + geom_line(linewidth = 1) + labs(x = "geo_distance_sc",
    y = "dissim_beta", color = "same_region", fill = "same_region") +
    theme_minimal()

Restricted to where each group actually has data: the same-region line is flat across its whole observed range (~0.53-0.55), while the different-region line rises steeply across its own range (~0.42-0.55). The two lines are not meaningfully compared at a shared distance value, since they don’t share a common range - only the difference in slopes is well supported.

4.1.3 Within-region distance effect

To estimate the within-region distance slope directly, without extrapolating across groups, mod_geodist_within_region was fit on same-region pairs only:

Code 
extended_summary(mod_geodist_within_region, highlight = TRUE, print.name = FALSE)
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 0.5) Intercept-normal(0, 1) phi-exponential(1) sd-exponential(2) dissim_beta ~ geo_distance_sc + (1 | mm(rec1, rec2)) 4000 4 1 1000 0 (0%) 0 270.76 727.738 123
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept 0.166 0.093 0.239 1.015 270.76 727.738
b_geo_distance_sc 0.070 0.033 0.106 1.001 7078.57 7956.087

4.1.4 Full model (region, locality, distance together)

Code 
extended_summary(mod_all, highlight = TRUE, print.name = FALSE)
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 0.5) Intercept-normal(0, 1) phi-exponential(1) sd-exponential(2) dissim_beta ~ same_region + same_locality + geo_distance_sc + (1 | mm(rec1, rec2)) 4000 4 1 1000 0 (0%) 0 202.281 316.368 123
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept 0.152 0.063 0.237 1.03 202.281 316.368
b_same_regionsame 0.039 -0.012 0.089 1 3412.251 5030.299
b_same_localitysame 0.053 0.028 0.078 1 4185.206 6139.623
b_geo_distance_sc 0.171 0.115 0.227 1.001 3225.264 4702.258

Code 
extended_summary(mod_both, highlight = TRUE, print.name = FALSE)
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 1) Intercept-normal(0, 2) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 2.5) dissim_beta ~ same_region + same_locality + (1 | mm(rec1, rec2)) 4000 4 1 1000 0 (0%) 0 215.76 463.087 123
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept 0.262 0.193 0.331 1.041 215.760 463.087
b_same_regionsame -0.104 -0.123 -0.085 1 5266.276 7306.824
b_same_localitysame -0.008 -0.022 0.007 1 10168.787 8880.007

mod_all’s same_region coefficient is no longer credible once locality and distance are both included (CI crosses zero) - its signal is carried through the interaction with distance (above), not as a distance-independent main effect, so this isn’t a contradiction of the interaction result. mod_all’s same_locality coefficient is credible but positive - the opposite sign from its (non-credible) direction in mod_both - which is a collinearity artifact rather than a real effect: same-locality pairs have geo_distance_sc close to zero by construction, so once distance enters the model this coefficient partly re-absorbs “distance ≈ 0” rather than a clean locality effect. mod_both (no distance) reproduces the original, simpler finding: same_region is credibly negative (pairs from the same region are more similar) and same_locality is not credible - locality adds nothing once region is in the model.

Note: b_Intercept shows Rhat slightly above 1.01 with reduced ESS in mod_all and mod_both (the fixed-effect and interaction coefficients are all well converged, Rhat ~1.00). Worth increasing iter or inspecting trace plots for the intercept before treating these as final.

4.1.5 Region and locality only (no distance)

Code 
# effect on the response scale
conditional_effects(mod_region, effects = "same_region")

Code 
conditional_effects(mod_locality, effects = "same_locality")

Takeaways

  • Region structures acoustic similarity, but not through a distance-independent effect - it works specifically through its interaction with geographic distance. The best-fitting model (mod_region_geo) shows a steep distance effect between regions and a much weaker one within regions.

  • Within a region, geographic distance has only a small effect on acoustic dissimilarity - the within-region-only model (mod_geodist_within_region) finds a credible but small positive slope (0.070, CI 0.033-0.106), and the effective within-region slope in mod_region_geo (≈0.05) is consistent with this: distance matters far less within a region than between regions, but not literally zero.

  • Between regions, geographic distance predicts dissimilarity strongly (0.481, CI 0.345-0.617 in mod_region_geo) - more distant regions are more acoustically distinct, an effect several times larger than the within-region one.

  • Locality adds nothing beyond region once distance is properly accounted for. mod_both (region + locality, no distance) shows same_locality is not credible - locality contributes nothing once region is included. mod_all (region + locality + distance together) shows both same_locality and geo_distance_sc as credible, but VIF values of 7.4 (region), 11.3 (distance), and 3.1 (locality) indicate severe collinearity among all three predictors when combined - mod_all’s individual coefficients should not be read as independent effects, since region, locality, and geographic distance all partly encode the same underlying spatial structure.

  • Caution on same-region vs. different-region pairs is needed at any single point. These two groups occupy near-disjoint ranges of geographic distance (regions are spatially clustered), so comparing their fitted lines at one distance value extrapolates outside the real data for at least one group. Only the difference in slopes is well supported; a difference in intercepts is not.

  • Practical conclusion: this argues for regional acoustic cohesion that holds fairly independently of within-region distance, alongside genuine divergence between regions that scales strongly with how far apart they are - consistent with strong differences between regions without unique, geography-independent dialects at finer scales. There are many variants within regions, but they remain more similar to each other than to variants from other regions, with only a modest additional pull from distance within the region itself.

5 Region-based feature variation

Code 
acous_feat <- readRDS("./data/processed/acoustic_features.RDS")


acous_feat$region <- as.factor(acous_feat$region)
acous_feat$recording <- as.factor(acous_feat$recording)

## sp.ent is a proportion (0-1) - rescale to the open interval for the
## beta family, and rename without a "." (same lesson as mm() earlier:
## brms's multivariate resp identifiers are built by stripping "." and
## "_" from the column name, so "sp.ent" and "sp_ent" would collide;
## use a clean name up front to avoid ambiguity)
rescale01 <- function(x, eps = 1e-4) (x * (1 - 2 * eps)) + eps
acous_feat$spent <- rescale01(acous_feat$sp.ent)

## -----------------------------------------------------------------------
## three response equations, one family each, sharing recording as a
## correlated random intercept: (1 | p | recording) uses the same ID
## letter "p" across all three so brms estimates a correlation matrix
## between the recording-level intercepts of duration, meanpeakf, and
## spectral entropy - i.e. do "noisier"/atypical recordings shift all
## three acoustic measures together, not just one at a time
## -----------------------------------------------------------------------

bf_duration  <- bf(duration  ~ region + (1 | p | recording), family = lognormal())
bf_meanpeakf <- bf(meanpeakf ~ region + (1 | p | recording), family = gaussian())
bf_spent     <- bf(spent     ~ region + (1 | p | recording), family = Beta())

priors <- c(
  prior(normal(0, 1),   class = "Intercept", resp = "duration"),
  prior(normal(0, 0.5), class = "b",         resp = "duration"),
  prior(exponential(2), class = "sd",        resp = "duration"),
  prior(exponential(1), class = "sigma",     resp = "duration"),

  prior(normal(0, 1),   class = "Intercept", resp = "meanpeakf"),
  prior(normal(0, 0.5), class = "b",         resp = "meanpeakf"),
  prior(exponential(2), class = "sd",        resp = "meanpeakf"),
  prior(exponential(1), class = "sigma",     resp = "meanpeakf"),

  prior(normal(0, 1),   class = "Intercept", resp = "spent"),
  prior(normal(0, 0.5), class = "b",         resp = "spent"),
  prior(exponential(2), class = "sd",        resp = "spent")
)

mod_mv_region <- brm(
  bf_duration + bf_meanpeakf + bf_spent,
  data = acous_feat,
  prior = priors,
  backend = "cmdstanr", threads = threading(8),
  chains = 4, cores = 4, iter = 4000, warmup = 1000,
  control = list(adapt_delta = 0.95),
  seed = 123,
  file = "./data/processed/brms_model_multivariate_region"
)

extended_summary(mod_mv_region, highlight = TRUE, print.name = FALSE)
<table class="table table-striped table-hover table-condensed table-responsive" style="font-size: 12px; width: auto !important; margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:left;">   </th>
   <th style="text-align:left;"> priors </th>
   <th style="text-align:left;"> formula </th>
   <th style="text-align:right;"> iterations </th>
   <th style="text-align:right;"> chains </th>
   <th style="text-align:right;"> thinning </th>
   <th style="text-align:right;"> warmup </th>
   <th style="text-align:left;"> diverg_transitions </th>
   <th style="text-align:left;"> rhats &gt; 1.05 </th>
   <th style="text-align:right;"> min_bulk_ESS </th>
   <th style="text-align:right;"> min_tail_ESS </th>
   <th style="text-align:left;"> seed </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> 1 </td>
   <td style="text-align:left;"> b-normal(0, 0.5)
b-normal(0, 0.5)
b-normal(0, 0.5)
Intercept-normal(0, 1)
Intercept-normal(0, 1)
Intercept-normal(0, 1)
L-lkj_corr_cholesky(1)
phi-gamma(0.01, 0.01)
sd-exponential(2)
sd-exponential(2)
sd-exponential(2)
sigma-exponential(1)
sigma-exponential(1) </td>
   <td style="text-align:left;"> list(duration = list(formula = duration ~ region + (1 | p | recording), pforms = list(), pfix = list(), family = list(family = "lognormal", link = "identity", linkfun = function (mu) 
link(mu, link = slink), linkinv = function (eta) 
inv_link(eta, link = slink), dpars = c("mu", "sigma"), type = "real", ybounds = c(0, Inf), closed = c(FALSE, NA), ad = c("weights", "subset", "cens", "trunc", "mi", "index"), specials = "logscale", link_sigma = "log"), resp = "duration", mecor = TRUE), meanpeakf = list(formula = meanpeakf ~ region + (1 | p | recording), pforms = list(), pfix = list(), family = list(family = "gaussian", link = "identity", linkfun = function (mu) 
link(mu, link = slink), linkinv = function (eta) 
inv_link(eta, link = slink), dpars = c("mu", "sigma"), type = "real", ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c("weights", "subset", "se", "cens", "trunc", "mi", "index"), normalized = c("_time_hom", "_time_het", "_lagsar", "_errorsar", "_fcor"), specials = c("residuals", "rescor")), resp = "meanpeakf", mecor = TRUE), spent = list(formula = spent ~ region + (1 | p | recording), pforms = list(), pfix = list(), family = list(family = "beta", link = "logit", linkfun = function (mu) 
link(mu, link = slink), linkinv = function (eta) 
inv_link(eta, link = slink), dpars = c("mu", "phi"), type = "real", ybounds = c(0, 1), closed = c(FALSE, FALSE), ad = c("weights", "subset", "cens", "trunc", "mi", "index"), link_phi = "log"), resp = "spent", mecor = TRUE)) </td>
   <td style="text-align:right;"> 4000 </td>
   <td style="text-align:right;"> 4 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 1000 </td>
   <td style="text-align:left;"> <span style="     ">0 (0%)</span> </td>
   <td style="text-align:left;"> <span style="     ">0</span> </td>
   <td style="text-align:right;"> 2576.34 </td>
   <td style="text-align:right;"> 4389.989 </td>
   <td style="text-align:left;"> 123 </td>
  </tr>
</tbody>
</table><table class="table table-striped table-hover table-condensed table-responsive" style="font-size: 12px; width: auto !important; margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:left;">   </th>
   <th style="text-align:right;"> Estimate </th>
   <th style="text-align:right;"> l-95% CI </th>
   <th style="text-align:right;"> u-95% CI </th>
   <th style="text-align:left;"> Rhat </th>
   <th style="text-align:right;"> Bulk_ESS </th>
   <th style="text-align:right;"> Tail_ESS </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;background-color: rgba(109, 205, 89, 77) !important;"> b_duration_Intercept </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> -0.730 </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> -0.935 </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> -0.528 </td>
   <td style="text-align:left;background-color: rgba(109, 205, 89, 77) !important;"> <span style="     ">1.001</span> </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> 2801.270 </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> 4658.267 </td>
  </tr>
  <tr>
   <td style="text-align:left;background-color: rgba(109, 205, 89, 77) !important;"> b_meanpeakf_Intercept </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> 0.807 </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> 0.341 </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> 1.279 </td>
   <td style="text-align:left;background-color: rgba(109, 205, 89, 77) !important;"> <span style="     ">1.001</span> </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> 3951.186 </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> 5959.060 </td>
  </tr>
  <tr>
   <td style="text-align:left;background-color: rgba(109, 205, 89, 77) !important;"> b_spent_Intercept </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> 2.324 </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> 2.172 </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> 2.472 </td>
   <td style="text-align:left;background-color: rgba(109, 205, 89, 77) !important;"> <span style="     ">1.001</span> </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> 3816.104 </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> 5879.942 </td>
  </tr>
  <tr>
   <td style="text-align:left;background-color: rgba(109, 205, 89, 77) !important;"> b_duration_regioncentralsubB </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> -0.353 </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> -0.568 </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> -0.131 </td>
   <td style="text-align:left;background-color: rgba(109, 205, 89, 77) !important;"> <span style="     ">1.001</span> </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> 2654.846 </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> 4569.530 </td>
  </tr>
  <tr>
   <td style="text-align:left;background-color: rgba(109, 205, 89, 77) !important;"> b_duration_regionnorthern </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> -0.545 </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> -0.850 </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> -0.233 </td>
   <td style="text-align:left;background-color: rgba(109, 205, 89, 77) !important;"> <span style="     ">1</span> </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> 3506.345 </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> 5287.439 </td>
  </tr>
  <tr>
   <td style="text-align:left;background-color: rgba(109, 205, 89, 77) !important;"> b_duration_regionsouthern </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> -0.474 </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> -0.683 </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> -0.264 </td>
   <td style="text-align:left;background-color: rgba(109, 205, 89, 77) !important;"> <span style="     ">1.001</span> </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> 2576.340 </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> 4389.989 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> b_meanpeakf_regioncentralsubB </td>
   <td style="text-align:right;"> 0.190 </td>
   <td style="text-align:right;"> -0.304 </td>
   <td style="text-align:right;"> 0.684 </td>
   <td style="text-align:left;"> <span style="     ">1.001</span> </td>
   <td style="text-align:right;"> 4292.430 </td>
   <td style="text-align:right;"> 6213.882 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> b_meanpeakf_regionnorthern </td>
   <td style="text-align:right;"> -0.456 </td>
   <td style="text-align:right;"> -1.143 </td>
   <td style="text-align:right;"> 0.215 </td>
   <td style="text-align:left;"> <span style="     ">1</span> </td>
   <td style="text-align:right;"> 6297.589 </td>
   <td style="text-align:right;"> 8204.507 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> b_meanpeakf_regionsouthern </td>
   <td style="text-align:right;"> 0.353 </td>
   <td style="text-align:right;"> -0.122 </td>
   <td style="text-align:right;"> 0.826 </td>
   <td style="text-align:left;"> <span style="     ">1.001</span> </td>
   <td style="text-align:right;"> 3963.790 </td>
   <td style="text-align:right;"> 6101.431 </td>
  </tr>
  <tr>
   <td style="text-align:left;background-color: rgba(109, 205, 89, 77) !important;"> b_spent_regioncentralsubB </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> -0.575 </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> -0.738 </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> -0.410 </td>
   <td style="text-align:left;background-color: rgba(109, 205, 89, 77) !important;"> <span style="     ">1.001</span> </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> 3976.678 </td>
   <td style="text-align:right;background-color: rgba(109, 205, 89, 77) !important;"> 6723.934 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> b_spent_regionnorthern </td>
   <td style="text-align:right;"> -0.207 </td>
   <td style="text-align:right;"> -0.436 </td>
   <td style="text-align:right;"> 0.020 </td>
   <td style="text-align:left;"> <span style="     ">1.001</span> </td>
   <td style="text-align:right;"> 4409.911 </td>
   <td style="text-align:right;"> 6945.017 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> b_spent_regionsouthern </td>
   <td style="text-align:right;"> -0.087 </td>
   <td style="text-align:right;"> -0.239 </td>
   <td style="text-align:right;"> 0.071 </td>
   <td style="text-align:left;"> <span style="     ">1.001</span> </td>
   <td style="text-align:right;"> 3733.371 </td>
   <td style="text-align:right;"> 6100.988 </td>
  </tr>
</tbody>
</table>

Code 
conditional_effects(mod_mv_region, effects = "region", resp = "duration")

Code 
conditional_effects(mod_mv_region, effects = "region", resp = "meanpeakf")

Code 
conditional_effects(mod_mv_region, effects = "region", resp = "spent")

Code 
## all pairwise region contrasts, one response at a time

emm_duration <- emmeans(mod_mv_region, ~region, resp = "duration")
emm_meanpeakf <- emmeans(mod_mv_region, ~region, resp = "meanpeakf")
emm_spent <- emmeans(mod_mv_region, ~region, resp = "spent")

contrasts_duration <- pairs(emm_duration)
contrasts_meanpeakf <- pairs(emm_meanpeakf)
contrasts_spent <- pairs(emm_spent)

# each of these returns a full posterior for every pairwise
# difference; summary() with HPD intervals is the Bayesian
# analogue of a CI table
summary(contrasts_duration, point.est = median)
contrast estimate lower.HPD upper.HPD
central subA - central subB 0.3541985 0.136489 0.572487
central subA - northern 0.5467765 0.230542 0.845921
central subA - southern 0.4748645 0.272601 0.690762
central subB - northern 0.1903080 -0.064800 0.474995
central subB - southern 0.1210120 0.046371 0.194343
northern - southern -0.0693550 -0.339819 0.184583
Code 
summary(contrasts_meanpeakf, point.est = median)
contrast estimate lower.HPD upper.HPD
central subA - central subB -0.1889290 -0.669702 0.3129130
central subA - northern 0.4546965 -0.182350 1.1739400
central subA - southern -0.3538705 -0.827519 0.1171150
central subB - northern 0.6443407 -0.024983 1.3629080
central subB - southern -0.1617815 -0.365023 0.0228440
northern - southern -0.8104185 -1.502559 -0.1593541
Code 
summary(contrasts_spent, point.est = median)
contrast estimate lower.HPD upper.HPD
central subA - central subB 0.5755130 0.4129400 0.7395750
central subA - northern 0.2073405 -0.0233093 0.4310660
central subA - southern 0.0878271 -0.0628694 0.2455020
central subB - northern -0.3686090 -0.5543530 -0.1792850
central subB - southern -0.4874686 -0.5502953 -0.4283280
northern - southern -0.1191689 -0.3001820 0.0621004
Code 
## regional marginal means, one panel per response

plot_emm <- function(emm_obj, title) {
    emm_df <- as.data.frame(emm_obj)
    ggplot(emm_df, aes(x = region, y = emmean, color = region)) +
        geom_pointrange(aes(ymin = lower.HPD, ymax = upper.HPD), linewidth = 0.8,
            size = 0.6) + scale_color_viridis_d(option = "mako", begin = 0.2,
        end = 0.8) + labs(title = title, x = NULL, y = "Estimated marginal mean") +
        theme_minimal(base_size = 12) + theme(legend.position = "none",
        axis.text.x = element_text(angle = 30, hjust = 1))
}

p_emm_duration <- plot_emm(emm_duration, "Duration")
p_emm_meanpeakf <- plot_emm(emm_meanpeakf, "Peak frequency")
p_emm_spent <- plot_emm(emm_spent, "Spectral entropy")

p_emm_duration + p_emm_meanpeakf + p_emm_spent

Code 
## -----------------------------------------------------------------------
## pairwise contrasts, one panel per response - colored by
## whether the 95% HPD interval excludes zero (credible) or not
## -----------------------------------------------------------------------

plot_contrasts <- function(contrast_obj, title) {
    contr_df <- as.data.frame(contrast_obj)
    contr_df$credible <- contr_df$lower.HPD > 0 | contr_df$upper.HPD <
        0

    ggplot(contr_df, aes(x = estimate, y = reorder(contrast, estimate),
        color = credible)) + geom_vline(xintercept = 0, linetype = "dashed",
        color = "grey50") + geom_pointrange(aes(xmin = lower.HPD,
        xmax = upper.HPD), linewidth = 0.8, size = 0.5) + scale_color_manual(values = c(`TRUE` = "#357BA2",
        `FALSE` = "grey70"), labels = c(`TRUE` = "credible", `FALSE` = "not credible")) +
        labs(title = title, x = "Estimated difference", y = NULL,
            color = NULL) + theme_minimal(base_size = 12)
}

p_contr_duration <- plot_contrasts(contrasts_duration, "Duration: pairwise contrasts")
p_contr_meanpeakf <- plot_contrasts(contrasts_meanpeakf, "Peak frequency: pairwise contrasts")
p_contr_spent <- plot_contrasts(contrasts_spent, "Spectral entropy: pairwise contrasts")

p_contr_duration/p_contr_meanpeakf/p_contr_spent

Takeaways

  • Duration: central subA is credibly longer than all three other regions (central subB: 0.354, CI 0.137-0.572; northern: 0.547, CI 0.231-0.846; southern: 0.475, CI 0.273-0.691). Among the rest, central subB is credibly longer than southern (0.121, CI 0.046-0.194), but not credibly different from northern (0.190, CI -0.065-0.475). northern vs. southern is not credible. Overall: central subA stands apart with the longest calls; the other three form a less clearly separated group.

  • Peak frequency: the only credible pairwise contrast is northern - southern (-0.810, CI -1.503 to -0.159) - northern has lower peak frequency than southern. This didn’t show up in the original vs.-reference table, since neither region differed credibly from central subA individually; the full pairwise comparison reveals a real difference the reference-coded coefficients masked. All other pairs (including both vs. central subA) are not credible.

  • Spectral entropy (log-odds scale): central subB is the outlier region here, credibly lower than every other region - vs. central subA (0.576, CI 0.413-0.740), vs. northern (-0.369, CI -0.554 to -0.179), vs. southern (-0.488, CI -0.550 to -0.428). No other pair is credible (central subA - northern and northern - southern both border zero but don’t clear it).

  • Putting the three features together: each region stands out on a different axis - central subA on duration (longest calls), central subB on spectral entropy (most tonal/least noisy), and the northern/southern split shows up only in peak frequency. No single region is acoustically distinct across all three features simultaneously; the regional structure is spread across different acoustic dimensions.

  • Practical note: the pairwise contrasts add real information beyond the original coefficient table specifically for meanpeakf - that response looked like “no regional signal” from the vs.-reference table alone, but the full 6-pair comparison shows northern and southern are genuinely different from each other, just not from the reference region.

─ 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-08-26
 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)
 askpass            1.2.1      2024-10-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)
 bitops             1.0-9      2024-10-03 [1] CRAN (R 4.5.2)
 boot               1.3-32     2025-08-29 [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)
 caret            * 7.0-1      2024-12-10 [1] CRAN (R 4.5.2)
 checkmate          2.3.4      2026-02-03 [1] CRAN (R 4.5.2)
 class              7.3-23     2025-01-01 [1] CRAN (R 4.5.2)
 classInt           0.4-11     2025-01-08 [1] CRAN (R 4.5.2)
 cli                3.6.6      2026-04-09 [1] CRAN (R 4.5.2)
 cluster            2.1.8.2    2026-02-05 [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)
 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)
 DBI                1.3.0      2026-02-25 [1] CRAN (R 4.5.2)
 deldir             2.0-4      2024-02-28 [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)
 e1071              1.7-17     2025-12-18 [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)
 energy             1.7-12     2024-08-24 [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)
 foreach            1.5.2      2022-02-02 [3] CRAN (R 4.1.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)
 future             1.70.0     2026-03-14 [1] CRAN (R 4.5.2)
 future.apply       1.20.2     2026-02-20 [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)
 globals            0.19.1     2026-03-13 [1] CRAN (R 4.5.2)
 glue               1.8.1      2026-04-17 [1] CRAN (R 4.5.2)
 goftest            1.2-3      2021-10-07 [3] CRAN (R 4.1.1)
 gower              1.0.2      2024-12-17 [1] CRAN (R 4.5.2)
 gridExtra          2.3.1      2026-06-25 [1] CRAN (R 4.5.2)
 gsl                2.1-9      2025-11-10 [1] CRAN (R 4.5.2)
 gtable             0.3.6      2024-10-25 [1] CRAN (R 4.5.2)
 hardhat            1.4.2      2025-08-20 [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)
 igraph             2.3.3      2026-06-26 [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)
 ipred              0.9-15     2024-07-18 [1] CRAN (R 4.5.2)
 iterators          1.0.14     2022-02-05 [3] CRAN (R 4.1.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)
 KernSmooth         2.23-26    2025-01-01 [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)
 lava               1.9.0      2026-04-05 [1] CRAN (R 4.5.2)
 lifecycle          1.0.5      2026-01-08 [1] CRAN (R 4.5.2)
 listenv            0.10.1     2026-03-10 [1] CRAN (R 4.5.2)
 loo                2.9.0      2025-12-23 [1] CRAN (R 4.5.2)
 lubridate          1.9.5      2026-02-04 [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)
 mgcv               1.9-4      2025-11-07 [1] CRAN (R 4.5.2)
 mlbench          * 2.1-7      2026-02-18 [1] CRAN (R 4.5.2)
 ModelMetrics       1.2.2.2    2020-03-17 [3] CRAN (R 4.0.1)
 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)
 nicheROVER         1.1.2      2023-10-13 [1] CRAN (R 4.5.2)
 nlme               3.1-168    2025-03-31 [1] CRAN (R 4.5.2)
 nnet               7.3-20     2025-01-01 [1] CRAN (R 4.5.2)
 ohun             * 1.0.4      2025-10-22 [1] CRAN (R 4.5.2)
 openssl            2.4.2      2026-06-09 [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)
 parallelly         1.46.1     2026-01-08 [1] CRAN (R 4.5.2)
 patchwork        * 1.3.2      2025-08-25 [1] CRAN (R 4.5.2)
 pbapply          * 1.7-4      2025-07-20 [1] CRAN (R 4.5.2)
 performance      * 0.17.1     2026-06-30 [1] CRAN (R 4.5.2)
 permute          * 0.9-10     2026-02-06 [1] CRAN (R 4.5.2)
 PhenotypeSpace   * 0.1.1      2026-08-17 [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)
 png                0.1-9      2026-03-15 [1] CRAN (R 4.5.2)
 polyclip           1.10-7     2024-07-23 [1] CRAN (R 4.5.2)
 posterior          1.6.1      2025-02-27 [1] CRAN (R 4.5.2)
 pROC               1.19.0.1   2025-07-31 [1] CRAN (R 4.5.2)
 processx           3.9.0      2026-04-22 [1] CRAN (R 4.5.2)
 prodlim            2026.03.11 2026-03-11 [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)
 randomForest     * 4.7-1.2    2024-09-22 [1] CRAN (R 4.5.2)
 raster             3.6-32     2025-03-28 [1] CRAN (R 4.5.2)
 rbibutils          2.4.1      2026-01-21 [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)
 Rdpack             2.6.6      2026-02-08 [1] CRAN (R 4.5.2)
 recipes            1.3.1      2025-05-21 [1] CRAN (R 4.5.2)
 reformulas         0.4.4      2026-02-02 [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)
 reticulate         1.45.0     2026-02-13 [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)
 rpart              4.1.24     2025-01-07 [1] CRAN (R 4.5.2)
 RSpectra           0.16-2     2024-07-18 [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)
 sf                 1.1-2      2026-07-23 [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)
 spatstat.data      3.1-9      2025-10-18 [1] CRAN (R 4.5.2)
 spatstat.explore   3.8-2      2026-07-27 [1] CRAN (R 4.5.2)
 spatstat.geom      3.8-2      2026-07-24 [1] CRAN (R 4.5.2)
 spatstat.random    3.5-1      2026-07-27 [1] CRAN (R 4.5.2)
 spatstat.sparse    3.2-0      2026-05-21 [1] CRAN (R 4.5.2)
 spatstat.univar    3.2-0      2026-05-18 [1] CRAN (R 4.5.2)
 spatstat.utils     3.2-4      2026-07-16 [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)
 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)
 T4transport        0.1.8      2026-01-11 [1] CRAN (R 4.5.2)
 tensor             1.5.1      2025-06-17 [1] CRAN (R 4.5.2)
 tensorA            0.36.2.1   2023-12-13 [1] CRAN (R 4.5.2)
 terra              1.9-11     2026-03-26 [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)
 timechange         0.4.0      2026-01-29 [1] CRAN (R 4.5.2)
 timeDate           4052.112   2026-01-28 [1] CRAN (R 4.5.2)
 tuneR            * 1.4.7      2024-04-17 [1] CRAN (R 4.5.2)
 umap             * 0.2.10.0   2023-02-01 [1] CRAN (R 4.5.2)
 units              1.0-1      2026-03-11 [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)
 vegan            * 2.7-5      2026-05-25 [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.

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