Cape parrot vocal dialects

Author

Marcelo Araya-Salas

Published

October 22, 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_USE_aug2026 - consolidated_sound_files_CPV_contact_calls.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_USE_aug2026 - consolidated_sound_files_CPV_contact_calls.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), ]


# table(dat$Delete.)
  • 3959 calls
  • 13 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 352
central subB Hoha Forest-iGxalingenwa Nature Reserve-Marutswa Forest-Polela Sawmill-Salt Spring Farm 481
northern Amorentia 529
southern Alice Pecan Orchard-Baddaford Farm-Hogsback-King William’s Town-Schwarzwald Forest-Stutterheim 2597
  • 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
5 central subB Hoha Forest
6 central subB iGxalingenwa Nature Reserve
9 central subB Marutswa Forest
10 central subB Polela Sawmill
11 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
12 southern Schwarzwald Forest
13 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")

## -----------------------------------------------------------------------
## 7a. region + distance, additive (no interaction, no
## locality): the cleanest direct test of whether region
## explains anything beyond distance alone, isolated from both
## the interaction term (which mechanically inflates VIFs, see
## mod_region_geo) and locality
## -----------------------------------------------------------------------

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

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

## -----------------------------------------------------------------------
## 7c. between-region subset model: the mirror of 7b above.
## same_region and geo_distance_sc are highly collinear in
## mod_region_geo (r = -0.90, see caveat in 'Region x distance
## interaction' below), so the interaction model's
## between-region slope is not a clean, region-independent
## estimate. Restricting to different-region pairs only removes
## the confound the same way the within-region subset does -
## region is held constant (always 'different'), so any
## remaining distance effect is unconfounded by region identity.
## -----------------------------------------------------------------------

dist_dat_diff_region <- dist_dat[dist_dat$same_region == "different",
    ]

mod_geodist_between_region <- brm(bf(dissim_beta ~ geo_distance_sc +
    (1 | mm(rec1, rec2))), data = dist_dat_diff_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_between_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_region_plus_geodist <- add_criterion(mod_region_plus_geodist,
    "loo")
mod_all <- add_criterion(mod_all, "loo")

# fit on subsets (same-region-only / different-region-only), so
# neither is comparable via loo_compare() to the models above -
# reported separately
mod_geodist_within_region <- add_criterion(mod_geodist_within_region,
    "loo")
mod_geodist_between_region <- add_criterion(mod_geodist_between_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_region_plus_geodist <- readRDS("./data/processed/brms_model_region_plus_geodist.rds")
mod_geodist_within_region <- readRDS("./data/processed/brms_model_geographic_distance_within_region.rds")
mod_geodist_between_region <- readRDS("./data/processed/brms_model_geographic_distance_between_region.rds")

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

kable(as.data.frame(loo_comp), digits = 2)
model elpd_diff se_diff p_worse diag_diff diag_elpd elpd_loo se_elpd_loo p_loo se_p_loo looic se_looic
mod_all 0.00 0.00 NA 10375.48 66.54 105.42 2.68 -20750.97 133.09
mod_region_geo -1.67 7.12 0.59 |elpd_diff| < 4 10373.81 66.10 103.84 2.60 -20747.62 132.20
mod_region_plus_geodist -15.87 6.89 0.99 10359.62 65.85 103.91 2.59 -20719.23 131.70
mod_both -16.52 6.14 1.00 10358.96 66.14 104.37 2.64 -20717.93 132.28
mod_region -16.90 6.81 0.99 10358.58 65.94 103.00 2.57 -20717.17 131.88
mod_geodist -43.52 11.02 1.00 10331.96 65.39 103.24 2.57 -20663.92 130.78
mod_locality -153.53 19.89 1.00 10221.96 65.31 104.35 2.61 -20443.92 130.63

mod_all and mod_region_geo are statistically indistinguishable at the top of the ranking (Δ = -1.67, SE 7.12, well under 1 SE), and both clearly outperform the rest: mod_region_geo is ahead of mod_region_plus_geodist (Δ = -14.19, SE 4.50, ~3.2 SEs), mod_both (Δ = -14.84, SE 5.63, ~2.6 SEs), mod_region (Δ = -15.23, SE 4.87, ~3.1 SEs), mod_geodist alone (Δ = -41.85, SE 9.10, ~4.6 SEs), and mod_locality (Δ = -151.85, SE 19.17, ~7.9 SEs). mod_all (all three predictors, no interaction) ties with mod_region_geo on predictive fit but doesn’t isolate the interaction the way mod_region_geo does, so mod_region_geo remains the more directly interpretable model for the question of whether distance matters differently within vs. between regions. mod_region_plus_geodist - the additive region + distance model, with no interaction and no locality - sits statistically indistinguishable from mod_both and mod_region (differences well under 1 SE from each other), but credibly behind mod_region_geo and mod_all: letting the distance slope change between same- and different-region pairs, or adding locality, both capture structure that a flat regional offset on top of a linear distance effect does not.

4.1.1 Does region explain anything beyond distance alone?

The LOO table above already hints at an answer - mod_geodist (distance only) sits far below every model that includes region in some form. To test this directly and cleanly, mod_region_plus_geodist adds same_region to geo_distance_sc as a flat, additive effect - no interaction, no locality - so it isolates the question from both the interaction term (which mechanically inflates VIFs, see below) and the locality confound.

Code 
extended_summary(mod_region_plus_geodist, 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 228.947 558.626 123
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept 0.274 0.208 0.339 1.02 228.947 558.626
b_same_regionsame -0.131 -0.165 -0.097 1.001 4561.496 6630.507
b_geo_distance_sc 0.028 -0.001 0.058 1.001 5179.590 7465.584

Code 
kable(as.data.frame(loo_compare(mod_region_plus_geodist, mod_geodist)),
    digits = 2)
model elpd_diff se_diff p_worse diag_diff diag_elpd elpd_loo se_elpd_loo p_loo se_p_loo looic se_looic
mod_region_plus_geodist 0.00 0.00 NA 10359.62 65.85 103.91 2.59 -20719.23 131.70
mod_geodist -27.66 7.35 1 10331.96 65.39 103.24 2.57 -20663.92 130.78

Adding same_region improves fit over distance alone by Δ = 27.66 (SE 7.35, ~3.8 SEs) - a credible, non-trivial improvement. This shows that a linear effect of (log) distance does not capture the structure in the data, but it does not show that region matters independently of distance (see Region and distance cannot be separated below). same_region’s own coefficient stays clearly credible (-0.131, CI -0.165 to -0.097) in this specification. Once same_region is in the model, though, geo_distance_sc’s coefficient shrinks from 0.124 (CI 0.108-0.139, mod_geodist alone) down to a value that barely clears zero (0.028, CI -0.001-0.058) - the two predictors are splitting credit for shared variance, not adding independent information.

Code 
check_collinearity(mod_region_plus_geodist)
Term VIF VIF_CI_low VIF_CI_high SE_factor Tolerance Tolerance_CI_low Tolerance_CI_high
same_region 3.499212 3.352116 3.655507 1.870618 0.2857786 0.2735599 0.298319
geo_distance_sc 3.499212 3.352116 3.655507 1.870618 0.2857786 0.2735599 0.298319
Code 
post <- as_draws_df(mod_region_plus_geodist)
cor(post$b_same_regionsame, post$b_geo_distance_sc)
[1] 0.8451162

Collinearity here (VIF = 3.50 for both terms) is more moderate than in mod_region_geo (VIF = 76.67/30.59) or even mod_all (VIF = 6.51/10.73) - removing the interaction term removes most of its mechanical VIF inflation. But the posterior correlation between the two coefficients is still high (r = 0.85, close to the raw r = -0.90 between the predictors themselves), so the two effects are still substantially entangled even in this simplest, cleanest specification.

This additive model is itself credibly worse than mod_region_geo (Δ = -14.19, SE 4.50, ~3.2 SEs) and mod_all (Δ = -15.87, SE 6.89, ~2.3 SEs) - a flat regional offset on top of distance doesn’t capture as much as letting the distance slope vary by region (the interaction) or adding locality.

Region and distance cannot be separated with this sampling. Same-region and different-region pairs do not overlap in geographic distance at all: every same-region pair has geo_distance_sc ≤ 0.12 and every different-region pair has geo_distance_sc ≥ 0.32 (see the range table under Region x distance interaction). In this dataset same_region is therefore an exact function of distance - equivalent to “distance below some cutoff between 0.12 and 0.32” - and not just a correlated predictor. Any LOO gain from adding same_region can be equally explained by a nonlinear, threshold-like effect of distance located in that gap, and no model fit to these data can tell the two apart (the region x distance interaction model is likewise equivalent to a broken-stick distance effect with its break inside the gap). What the data do support is that dissimilarity does not increase with distance within regions (see Within-region distance effect) and rises sharply at the spatial scale separating them. Separating a regional effect from a distance threshold would require recordings from nearby sites in different regions, or distant sites in the same region.

4.1.2 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.36) is low concern, same_region (VIF = 6.51) is moderate
## concern, and geo_distance_sc (VIF = 10.73) 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 6.509847 6.214338 6.822103 2.551440 0.1536134 0.1465824 0.1609182
same_locality 3.356687 3.216648 3.505574 1.832126 0.2979128 0.2852600 0.3108826
geo_distance_sc 10.728885 10.225427 11.259818 3.275498 0.0932063 0.0888114 0.0977954

Collinearity results

  • geo_distance_sc shows severe collinearity (VIF = 10.73, tolerance = 0.09) — over 90% of its variance is explained by same_region and same_locality combined.
  • same_region also shows concerning collinearity (VIF = 6.51, tolerance = 0.15).
  • same_locality is the least collinear of the three (VIF = 3.36, tolerance = 0.30) and falls within the conventional low-concern range (VIF < 5), though it is not fully independent of the other two.
  • 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.3 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 378.3 1172.638 123
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept 0.048 -0.054 0.153 1.004 378.300 1172.638
b_same_regionsame 0.091 0.007 0.176 1.001 1184.339 2687.030
b_geo_distance_sc 0.405 0.271 0.541 1.002 1152.951 2619.029
b_same_regionsame:geo_distance_sc -0.394 -0.532 -0.258 1.002 1207.548 2755.544

The interaction is the key result: geo_distance_sc (0.405, CI 0.271-0.541) is the distance slope for different-region pairs - dissimilarity rises steeply with distance. The interaction term same_regionsame:geo_distance_sc (-0.394, CI -0.532 to -0.258) cancels almost all of that slope for same-region pairs: the implied within-region slope (0.405 - 0.394 = 0.011, 95% CI -0.017 to 0.041) is not credibly different from zero. Geographic distance predicts acoustic dissimilarity strongly between regions, and has no detectable effect within them.

Caveat: same_region and geo_distance_sc are themselves highly collinear.

Code 
dist_dat <- readRDS("./data/processed/dist_dat.RDS")
cor(as.numeric(dist_dat$same_region == "same"), dist_dat$geo_distance_sc)
[1] -0.8995788
Code 
check_collinearity(mod_region_geo)
Term VIF VIF_CI_low VIF_CI_high SE_factor Tolerance Tolerance_CI_low Tolerance_CI_high
same_region 76.67036 72.91710 80.61949 8.756161 0.0130429 0.0124039 0.0137142
geo_distance_sc 30.59326 29.11081 32.15388 5.531117 0.0326869 0.0311004 0.0343515

This model’s two base predictors correlate at r = -0.90 - same_region is pair-coded (0/1) and geo_distance_sc is the pairwise distance it is almost standing in for. check_collinearity() reports VIF = 76.67 for same_region and 30.59 for geo_distance_sc, but flags that this is inflated by the presence of the same_region x geo_distance_sc interaction term itself rather than a literal severity score (centering or dropping the interaction would shrink it); the raw bivariate correlation above is the more interpretable diagnostic. Practically, this means the model has little power to distinguish “these two recordings are in the same region” from “these two recordings are geographically close” - they are almost the same piece of information in this dataset - so the strong between-region distance slope reported above should be read as a property of the combined region/distance structure, not as cleanly isolating a pure, region-independent distance effect.

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.32 0.61 0.87
same -0.67 -0.33 0.12

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), while the different-region line rises from ~0.54 to ~0.60 across its own range. 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.4 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 281.032 602.054 123
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept 0.139 0.07 0.212 1.027 281.032 602.054
b_geo_distance_sc 0.034 0.00 0.070 1.001 6628.191 8650.446

4.1.5 Between-region distance effect

mod_region_geo’s between-region distance slope is estimated from a model where same_region and geo_distance_sc are themselves highly collinear (r = -0.90, see caveat above), so it isn’t a clean, region-independent estimate. As the mirror of the within-region subset above, mod_geodist_between_region was fit on different-region pairs only - region is held constant (always “different”), which removes that confound the same way the within-region subset removes it for same-region pairs:

Code 
extended_summary(mod_geodist_between_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 1658.725 3286.347 123
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept -0.035 -0.207 0.138 1.003 1658.725 3286.347
b_geo_distance_sc 0.495 0.238 0.749 1.002 2162.476 3793.423

This unconfounded estimate (geo_distance_sc = 0.495, CI 0.238-0.749) is credible and, if anything, slightly larger than mod_region_geo’s confounded between-region slope (0.405, CI 0.271-0.541) - so the interaction model isn’t overstating the between-region distance effect. Between the two subset models, the pattern is clear and consistent: distance predicts dissimilarity credibly among different-region pairs, and not at all among same-region pairs.

4.1.6 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 243.638 660.826 123
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept 0.206 0.134 0.278 1.003 243.638 660.826
b_same_regionsame -0.039 -0.085 0.008 1.001 3251.354 4808.902
b_same_localitysame 0.071 0.047 0.096 1.001 3909.844 6813.102
b_geo_distance_sc 0.152 0.101 0.204 1.001 3254.262 4910.639

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, 0.5) Intercept-normal(0, 1) phi-exponential(1) sd-exponential(2) dissim_beta ~ same_region + same_locality + (1 | mm(rec1, rec2)) 4000 4 1 1000 0 (0%) 0 144.882 384.133 123
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_Intercept 0.297 0.231 0.371 1.027 144.882 384.133
b_same_regionsame -0.164 -0.183 -0.146 1.002 4066.575 6433.716
b_same_localitysame 0.012 -0.002 0.026 1 7673.612 8455.440

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 and positive (0.071, CI 0.047-0.096) - same-locality pairs more dissimilar - whereas in mod_both the same positive direction is small and not credible (0.012, CI -0.002-0.026). Becoming credible only once distance is added 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.

Can mod_all actually tell these three effects apart? Besides the VIFs already reported, the posterior itself shows how entangled the three coefficients are - the correlation between their posterior draws is a direct measure of how much the model can trade one effect off against another, independent of any VIF formula:

Code 
post <- as_draws_df(mod_all)
kable(data.frame(pair = c("same_region vs. geo_distance_sc", "same_locality vs. geo_distance_sc",
    "same_region vs. same_locality"), posterior_r = c(cor(post$b_same_regionsame,
    post$b_geo_distance_sc), cor(post$b_same_localitysame, post$b_geo_distance_sc),
    cor(post$b_same_regionsame, post$b_same_localitysame))), digits = 2)
pair posterior_r
same_region vs. geo_distance_sc 0.91
same_locality vs. geo_distance_sc 0.82
same_region vs. same_locality 0.68

b_same_regionsame and b_geo_distance_sc trade off almost one-for-one (posterior r = 0.91) - the data can shift credit between “same region” and “distance” with little cost to fit, which is exactly what the VIFs predicted. b_same_localitysame is similarly entangled with geo_distance_sc (r = 0.82) and with b_same_regionsame (r = 0.68). None of these three coefficients is estimated as if the other two didn’t exist, so mod_all’s individual estimates - credible or not - describe how the three predictors divide up shared variance in this particular model, not independent, “all-else-equal” effects. The subset models above (within- vs. between-region) sidestep this by removing the confound at the data level instead of trying to partition it statistically.

Note: b_Intercept shows reduced ESS in both mod_all (244 bulk / 661 tail) and mod_both (145 bulk / 384 tail) relative to the other coefficients, with Rhat also slightly elevated in mod_both (1.027); mod_all’s intercept Rhat is within range (1.003). The fixed-effect and interaction coefficients themselves are all well converged (Rhat ~1.00). Worth increasing iter or inspecting trace plots for the intercept before treating these as final.

4.1.7 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 and distance cannot be separated with the current sampling. Same-region and different-region pairs do not overlap in geographic distance, so same_region is equivalent to “below a certain distance” in this dataset. Adding region to a distance-only model improves fit (mod_region_plus_geodist vs. mod_geodist: Δ = 27.66, SE 7.35, ~3.8 SEs; same_region = -0.131, CI -0.165 to -0.097), which shows that the relationship is not a simple linear function of (log) distance - but it does not show that region matters independently of distance. A true regional effect and a distance threshold that coincides with region boundaries are indistinguishable here.

  • The best-supported models let the distance effect differ between same- and different-region pairs. mod_region_geo is statistically tied with mod_all for the best out-of-sample fit (Δ = -1.67, SE 7.12) and clearly ahead of the simpler models; since mod_all doesn’t include the interaction term, mod_region_geo is the model that directly addresses whether distance operates differently within vs. between regions. It shows a steep distance effect between regions and no detectable effect within them - though same_region and geo_distance_sc are themselves highly collinear (r = -0.90, see caveat above), so this isn’t a clean separation of two independent effects.

  • Within a region, geographic distance has no detectable effect on acoustic dissimilarity. The within-region-only model (mod_geodist_within_region) finds a small positive slope that is not credibly different from zero (0.034, CI -0.0004-0.070), and the implied within-region slope from the interaction model (mod_region_geo: 0.011, CI -0.017-0.041) agrees: once pairs are restricted to the same region, distance stops mattering.

  • Between regions, geographic distance predicts dissimilarity strongly (0.405, CI 0.271-0.541 in mod_region_geo), an effect that stands in clear contrast to the null within-region one. Because same_region and geo_distance_sc are collinear in mod_region_geo (r = -0.90), this slope was re-estimated on different-region pairs only (mod_geodist_between_region), where region is held constant and the confound is removed by design: the unconfounded estimate (0.495, CI 0.238-0.749) agrees, and is if anything larger - the between-region distance effect is not an artifact of that collinearity.

  • mod_all cannot cleanly separate region, locality, and distance. Beyond the VIFs, the posterior draws of its three coefficients are highly correlated with each other (same_region-geo_distance_sc: r = 0.91; same_locality-geo_distance_sc: r = 0.82; same_region-same_locality: r = 0.68) - the model can trade credit between them at little cost to fit. This is why the within-/between-region subset models above, which remove the confound at the data level, are the more trustworthy source for the distance effect specifically.

  • 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 6.51 (region), 10.73 (distance), and 3.36 (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: calls are cohesive within regions regardless of distance, and dissimilarity rises sharply at the spatial scale that separates regions - consistent with regional dialects rather than gradual change with distance at finer scales. There are many variants within regions, but they remain more similar to each other than to variants from other regions, with no additional pull from distance within the region itself (that slope is not statistically distinguishable from zero). With the current sampling, however, a true regional effect cannot be distinguished from a distance threshold that coincides with region boundaries, since same- and different-region pairs do not overlap in distance.

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)
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 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) 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)) 4000 4 1 1000 0 (0%) 0 2576.34 4389.989 123
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_duration_Intercept -0.730 -0.935 -0.528 1.001 2801.270 4658.267
b_meanpeakf_Intercept 0.807 0.341 1.279 1.001 3951.186 5959.060
b_spent_Intercept 2.324 2.172 2.472 1.001 3816.104 5879.942
b_duration_regioncentralsubB -0.353 -0.568 -0.131 1.001 2654.846 4569.530
b_duration_regionnorthern -0.545 -0.850 -0.233 1 3506.345 5287.439
b_duration_regionsouthern -0.474 -0.683 -0.264 1.001 2576.340 4389.989
b_meanpeakf_regioncentralsubB 0.190 -0.304 0.684 1.001 4292.430 6213.882
b_meanpeakf_regionnorthern -0.456 -1.143 0.215 1 6297.589 8204.507
b_meanpeakf_regionsouthern 0.353 -0.122 0.826 1.001 3963.790 6101.431
b_spent_regioncentralsubB -0.575 -0.738 -0.410 1.001 3976.678 6723.934
b_spent_regionnorthern -0.207 -0.436 0.020 1.001 4409.911 6945.017
b_spent_regionsouthern -0.087 -0.239 0.071 1.001 3733.371 6100.988

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.136-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.487, 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.6.1 (2026-06-24)
 os       Ubuntu 22.04.5 LTS
 system   x86_64, linux-gnu
 ui       X11
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       America/Costa_Rica
 date     2026-10-22
 pandoc   3.8.3 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/x86_64/ (via rmarkdown)
 quarto   1.7.31 @ /usr/local/bin/quarto

─ Packages ───────────────────────────────────────────────────────────────────
 package          * version    date (UTC) lib source
 abind              1.4-8      2024-09-12 [2] CRAN (R 4.6.1)
 ape                5.8-1      2024-12-16 [1] CRAN (R 4.6.1)
 arrayhelpers       1.1-2      2026-07-24 [2] CRAN (R 4.6.1)
 askpass            1.2.1      2024-10-04 [1] CRAN (R 4.6.1)
 backports          1.5.1      2026-04-03 [2] CRAN (R 4.6.1)
 bayesplot          1.15.0     2025-12-12 [2] CRAN (R 4.6.1)
 bitops             1.1-0      2026-07-30 [2] CRAN (R 4.6.1)
 boot               1.3-32     2025-08-29 [4] CRAN (R 4.5.1)
 bridgesampling     1.2-1      2025-11-19 [2] CRAN (R 4.6.1)
 brio               1.1.5      2024-04-24 [2] CRAN (R 4.6.1)
 brms             * 2.23.0     2025-09-09 [2] CRAN (R 4.6.1)
 brmsish          * 1.0.0      2026-07-14 [2] Github (maRce10/brmsish@81ab826)
 Brobdingnag        1.2-9      2022-10-19 [2] CRAN (R 4.6.1)
 cachem             1.1.0      2024-05-16 [1] CRAN (R 4.6.1)
 car                3.1-5      2026-02-03 [2] CRAN (R 4.6.1)
 carData            3.0-6      2026-01-30 [2] CRAN (R 4.6.1)
 caret            * 7.0-1      2024-12-10 [2] CRAN (R 4.6.1)
 checkmate          2.3.4      2026-02-03 [2] CRAN (R 4.6.1)
 class              7.3-24     2026-08-03 [1] CRAN (R 4.6.1)
 classInt           0.4-11     2025-01-08 [2] CRAN (R 4.6.1)
 cli                3.6.6      2026-04-09 [1] CRAN (R 4.6.1)
 cluster            2.1.8.3    2026-07-30 [1] CRAN (R 4.6.1)
 cmdstanr           0.9.0.9001 2026-07-14 [2] Github (stan-dev/cmdstanr@2763c0f)
 coda               0.19-4.1   2024-01-31 [1] CRAN (R 4.6.1)
 codetools          0.2-20     2024-03-31 [2] CRAN (R 4.6.1)
 cowplot            1.2.0      2025-07-07 [2] CRAN (R 4.6.1)
 crayon             1.5.3      2024-06-20 [2] CRAN (R 4.6.1)
 curl               8.0.0      2026-08-25 [1] CRAN (R 4.6.1)
 data.table         1.18.4     2026-05-06 [2] CRAN (R 4.6.1)
 DBI                1.3.0      2026-02-25 [2] CRAN (R 4.6.1)
 deldir             2.0-4      2024-02-28 [2] CRAN (R 4.6.1)
 devtools           2.5.2      2026-04-30 [2] CRAN (R 4.6.1)
 digest             0.6.39     2025-11-19 [1] CRAN (R 4.6.1)
 distributional     0.8.1      2026-06-27 [2] CRAN (R 4.6.1)
 dplyr              1.2.1      2026-04-03 [1] CRAN (R 4.6.1)
 dtw                1.23-3     2026-06-09 [2] CRAN (R 4.6.1)
 e1071              1.7-17     2025-12-18 [2] CRAN (R 4.6.1)
 ellipsis           0.3.3      2026-04-04 [2] CRAN (R 4.6.1)
 emmeans          * 2.0.4      2026-07-15 [2] CRAN (R 4.6.1)
 energy             1.7-12     2024-08-24 [2] CRAN (R 4.6.1)
 estimability       2.0.0      2026-06-26 [2] CRAN (R 4.6.1)
 evaluate           1.0.5      2025-08-27 [1] CRAN (R 4.6.1)
 farver             2.1.2      2024-05-13 [1] CRAN (R 4.6.1)
 fastmap            1.2.0      2024-05-15 [1] CRAN (R 4.6.1)
 fftw               1.0-9      2024-09-20 [2] CRAN (R 4.6.1)
 foreach            1.5.2      2022-02-02 [1] CRAN (R 4.6.1)
 formatR          * 1.14       2023-01-17 [2] CRAN (R 4.6.1)
 Formula            1.2-6      2026-08-03 [2] CRAN (R 4.6.1)
 fs                 2.1.0      2026-04-18 [1] CRAN (R 4.6.1)
 future             1.75.0     2026-07-20 [2] CRAN (R 4.6.1)
 future.apply       1.20.2     2026-02-20 [2] CRAN (R 4.6.1)
 generics           0.1.4      2025-05-09 [1] CRAN (R 4.6.1)
 geosphere        * 1.6-8      2026-04-05 [2] CRAN (R 4.6.1)
 ggcorrplot       * 0.3.0      2026-07-24 [2] CRAN (R 4.6.1)
 ggdist             3.3.3      2025-04-23 [2] CRAN (R 4.6.1)
 ggplot2          * 4.0.3      2026-04-22 [1] CRAN (R 4.6.1)
 globals            0.19.1     2026-03-13 [2] CRAN (R 4.6.1)
 glue               1.8.1      2026-04-17 [1] CRAN (R 4.6.1)
 goftest            1.2-3      2021-10-07 [2] CRAN (R 4.6.1)
 gower              1.0.2      2024-12-17 [2] CRAN (R 4.6.1)
 gridExtra          2.3.1      2026-06-25 [2] CRAN (R 4.6.1)
 gsl                2.1-9      2025-11-10 [2] CRAN (R 4.6.1)
 gtable             0.3.6      2024-10-25 [1] CRAN (R 4.6.1)
 hardhat            1.4.3      2026-04-04 [2] CRAN (R 4.6.1)
 htmltools          0.5.9      2025-12-04 [1] CRAN (R 4.6.1)
 htmlwidgets        1.6.4      2023-12-06 [1] CRAN (R 4.6.1)
 httr               1.4.9      2026-09-01 [1] CRAN (R 4.6.1)
 igraph             2.3.4      2026-09-30 [1] CRAN (R 4.6.1)
 inline             0.3.21     2025-01-09 [2] CRAN (R 4.6.1)
 insight            1.5.2      2026-06-28 [2] CRAN (R 4.6.1)
 ipred              0.9-15     2024-07-18 [2] CRAN (R 4.6.1)
 iterators          1.0.14     2022-02-05 [1] CRAN (R 4.6.1)
 jsonlite           2.0.0      2025-03-27 [1] CRAN (R 4.6.1)
 kableExtra         1.4.1      2026-07-08 [2] CRAN (R 4.6.1)
 KernSmooth         2.23-26    2025-01-01 [4] CRAN (R 4.4.2)
 knitr            * 1.52       2026-09-06 [1] CRAN (R 4.6.1)
 labeling           0.4.3      2023-08-29 [1] CRAN (R 4.6.1)
 lattice          * 0.22-9     2026-02-09 [4] CRAN (R 4.5.2)
 lava               1.9.2      2026-06-30 [2] CRAN (R 4.6.1)
 lifecycle          1.0.5      2026-01-08 [1] CRAN (R 4.6.1)
 listenv            1.0.0      2026-06-22 [2] CRAN (R 4.6.1)
 loo                2.10.1     2026-07-24 [2] CRAN (R 4.6.1)
 lubridate          1.9.5      2026-02-04 [1] CRAN (R 4.6.1)
 magrittr           2.0.5      2026-04-04 [1] CRAN (R 4.6.1)
 MASS               7.3-66     2026-07-15 [2] CRAN (R 4.6.1)
 Matrix             1.7-6      2026-07-25 [1] CRAN (R 4.6.1)
 matrixStats        1.5.0      2025-01-07 [2] CRAN (R 4.6.1)
 memoise            2.0.1      2021-11-26 [1] CRAN (R 4.6.1)
 mgcv               1.9-4      2025-11-07 [4] CRAN (R 4.5.2)
 mlbench          * 2.1-11     2026-08-14 [1] CRAN (R 4.6.1)
 ModelMetrics       1.2.2.2    2020-03-17 [2] CRAN (R 4.6.1)
 mvtnorm            1.4-2      2026-07-12 [1] CRAN (R 4.6.1)
 NatureSounds     * 1.0.5      2025-01-17 [2] CRAN (R 4.6.1)
 nicheROVER         1.1.2      2023-10-13 [2] CRAN (R 4.6.1)
 nlme               3.1-170    2026-07-15 [1] CRAN (R 4.6.1)
 nnet               7.3-21     2026-08-03 [1] CRAN (R 4.6.1)
 ohun             * 1.0.4      2025-10-22 [2] CRAN (R 4.6.1)
 openssl            2.4.2      2026-06-09 [1] CRAN (R 4.6.1)
 otel               0.2.0      2025-08-29 [1] CRAN (R 4.6.1)
 parallelly         1.48.0     2026-06-29 [2] CRAN (R 4.6.1)
 patchwork        * 1.3.2      2025-08-25 [1] CRAN (R 4.6.1)
 pbapply          * 1.7-4      2025-07-20 [2] CRAN (R 4.6.1)
 performance      * 0.17.1     2026-06-30 [2] CRAN (R 4.6.1)
 permute          * 0.9-10     2026-02-06 [2] CRAN (R 4.6.1)
 PhenotypeSpace   * 0.1.1      2026-07-24 [2] CRAN (R 4.6.1)
 pillar             1.11.1     2025-09-17 [1] CRAN (R 4.6.1)
 pkgbuild           1.4.8      2025-05-26 [2] CRAN (R 4.6.1)
 pkgconfig          2.0.3      2019-09-22 [1] CRAN (R 4.6.1)
 pkgload            1.5.3      2026-06-15 [2] CRAN (R 4.6.1)
 plyr               1.8.9      2023-10-02 [1] CRAN (R 4.6.1)
 png                0.1-9      2026-03-15 [2] CRAN (R 4.6.1)
 polyclip           1.10-7     2024-07-23 [2] CRAN (R 4.6.1)
 posterior          1.7.0      2026-04-01 [2] CRAN (R 4.6.1)
 pROC               1.19.0.1   2025-07-31 [2] CRAN (R 4.6.1)
 processx           3.9.0      2026-04-22 [2] CRAN (R 4.6.1)
 prodlim            2026.03.11 2026-03-11 [2] CRAN (R 4.6.1)
 proxy              0.4-29     2025-12-29 [2] CRAN (R 4.6.1)
 ps                 1.9.3      2026-04-20 [2] CRAN (R 4.6.1)
 purrr              1.2.2      2026-04-10 [1] CRAN (R 4.6.1)
 QuickJSR           1.10.0     2026-05-17 [3] CRAN (R 4.6.1)
 R6                 2.6.1      2025-02-15 [1] CRAN (R 4.6.1)
 randomForest     * 4.7-1.2    2024-09-22 [2] CRAN (R 4.6.1)
 raster             3.6-32     2025-03-28 [2] CRAN (R 4.6.1)
 rbibutils          2.4.1      2026-01-21 [2] CRAN (R 4.6.1)
 RColorBrewer       1.1-3      2022-04-03 [1] CRAN (R 4.6.1)
 Rcpp             * 1.1.2      2026-07-05 [1] CRAN (R 4.6.1)
 RcppParallel       6.2.0      2026-07-30 [2] CRAN (R 4.6.1)
 RCurl              1.98-1.19  2026-06-03 [2] CRAN (R 4.6.1)
 Rdpack             2.6.6      2026-02-08 [2] CRAN (R 4.6.1)
 recipes            1.3.3      2026-05-30 [2] CRAN (R 4.6.1)
 reformulas         0.4.4      2026-02-02 [2] CRAN (R 4.6.1)
 reshape2           1.4.5      2025-11-12 [2] CRAN (R 4.6.1)
 reticulate         1.46.0     2026-04-09 [1] CRAN (R 4.6.1)
 rjson              0.2.23     2024-09-16 [2] CRAN (R 4.6.1)
 rlang              1.3.0      2026-07-05 [1] CRAN (R 4.6.1)
 rmarkdown          2.32       2026-09-01 [1] CRAN (R 4.6.1)
 rpart              4.1.27     2026-03-27 [4] CRAN (R 4.5.3)
 RSpectra           0.16-2     2024-07-18 [1] CRAN (R 4.6.1)
 rstan              2.32.7     2025-03-10 [2] CRAN (R 4.6.1)
 rstantools         2.7.0      2026-07-26 [2] CRAN (R 4.6.1)
 rstudioapi         0.19.0     2026-06-11 [1] CRAN (R 4.6.1)
 S7                 0.2.2      2026-04-22 [1] CRAN (R 4.6.1)
 scales             1.4.0      2025-04-24 [1] CRAN (R 4.6.1)
 seewave          * 2.2.4      2025-08-19 [2] CRAN (R 4.6.1)
 sessioninfo        1.2.4      2026-06-04 [2] CRAN (R 4.6.1)
 sf                 1.1-2      2026-07-23 [2] CRAN (R 4.6.1)
 signal             1.8-1      2024-06-26 [2] CRAN (R 4.6.1)
 sketchy            1.0.7      2026-10-07 [1] CRANs (R 4.6.1)
 sp                 2.2-3      2026-07-19 [2] CRAN (R 4.6.1)
 spatstat.data      3.1-9      2025-10-18 [2] CRAN (R 4.6.1)
 spatstat.explore   3.8-2      2026-07-27 [2] CRAN (R 4.6.1)
 spatstat.geom      3.8-2      2026-07-24 [2] CRAN (R 4.6.1)
 spatstat.random    3.5-1      2026-07-27 [2] CRAN (R 4.6.1)
 spatstat.sparse    3.2-0      2026-05-21 [2] CRAN (R 4.6.1)
 spatstat.univar    3.2-0      2026-05-18 [2] CRAN (R 4.6.1)
 spatstat.utils     3.2-4      2026-07-16 [2] CRAN (R 4.6.1)
 StanHeaders        2.32.10    2024-07-15 [2] CRAN (R 4.6.1)
 stringi            1.8.9      2026-08-04 [1] CRAN (R 4.6.1)
 stringr            1.6.0      2025-11-04 [1] CRAN (R 4.6.1)
 survival           3.8-9      2026-07-08 [2] CRAN (R 4.6.1)
 svglite            2.2.2      2025-10-21 [2] CRAN (R 4.6.1)
 svUnit             1.0.8      2025-08-26 [2] CRAN (R 4.6.1)
 systemfonts        1.3.2      2026-03-05 [1] CRAN (R 4.6.1)
 T4transport        0.1.8      2026-01-11 [2] CRAN (R 4.6.1)
 tensor             1.5.1      2025-06-17 [2] CRAN (R 4.6.1)
 tensorA            0.36.2.1   2023-12-13 [2] CRAN (R 4.6.1)
 terra              1.9-34     2026-06-19 [2] CRAN (R 4.6.1)
 testthat           3.3.2      2026-01-11 [2] CRAN (R 4.6.1)
 textshaping        1.0.5      2026-03-06 [2] CRAN (R 4.6.1)
 tibble             3.3.1      2026-01-11 [1] CRAN (R 4.6.1)
 tidybayes          3.0.7      2024-09-15 [2] CRAN (R 4.6.1)
 tidyr              1.3.2      2025-12-19 [1] CRAN (R 4.6.1)
 tidyselect         1.2.1      2024-03-11 [1] CRAN (R 4.6.1)
 timechange         0.4.0      2026-01-29 [1] CRAN (R 4.6.1)
 timeDate           4052.112   2026-01-28 [2] CRAN (R 4.6.1)
 tuneR            * 1.4.7      2024-04-17 [2] CRAN (R 4.6.1)
 umap             * 0.2.10.0   2023-02-01 [1] CRAN (R 4.6.1)
 units              1.0-1      2026-03-11 [2] CRAN (R 4.6.1)
 usethis            3.2.1      2025-09-06 [2] CRAN (R 4.6.1)
 V8                 8.2.0      2026-04-21 [2] CRAN (R 4.6.1)
 vctrs              0.7.3      2026-04-11 [1] CRAN (R 4.6.1)
 vegan            * 2.7-5      2026-05-25 [2] CRAN (R 4.6.1)
 viridis          * 0.6.5      2024-01-29 [2] CRAN (R 4.6.1)
 viridisLite      * 0.4.3      2026-02-04 [1] CRAN (R 4.6.1)
 warbleR          * 1.1.37     2025-10-22 [2] CRAN (R 4.6.1)
 withr              3.0.3      2026-06-19 [1] CRAN (R 4.6.1)
 xfun               0.61       2026-09-16 [1] CRAN (R 4.6.1)
 xml2               1.6.0      2026-06-22 [1] CRAN (R 4.6.1)
 xtable             1.8-8      2026-02-22 [2] CRAN (R 4.6.1)
 yaml               2.3.12     2025-12-10 [1] CRAN (R 4.6.1)

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

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