Code
# options to customize chunk outputs
knitr::opts_chunk$set(
message = FALSE
)Statistical analysis
Vocal character displacement in southern capuchinos
Purpose
- Analyze acoustic divergence and the role of sympatry in simple songs of Southern capuchinos
Load packages and custom functions
Code
# install knitr package if not installed
if (!requireNamespace("sketchy", quietly = TRUE)) {
install.packages("sketchy")
}
packages <- c(
"knitr",
"dplyr",
"tidyverse",
"vegan",
"geosphere",
"ecodist",
"brms",
"brmsish",
"PhenotypeSpace",
"kableExtra",
"ggplot2",
"ggtext"
)
# install/load packages
sketchy::load_packages(packages = packages)
options(knitr.kable.NA = '')
print <- function(x, row.names = FALSE) {
kb <- kable(x, row.names = row.names, digits = 4, "html")
kb <- kable_styling(kb,
bootstrap_options = c("striped", "hover", "condensed", "responsive"))
scroll_box(kb, width = "100%")
}
# set theme globally
theme_set(theme_classic(base_size = 20))
get_stable_loadings <- function(pca, pcs = 4, B = 1000, cum_threshold = 0.5,
freq_threshold = 0.75, seed = 123) {
set.seed(seed)
# Reconstruct centered data
X <- pca$x %*% t(pca$rotation)
p <- ncol(pca$rotation)
orig_rot <- pca$rotation[, 1:pcs]
boot_loadings <- array(NA, dim = c(p, pcs, B))
# ----------------------------- Bootstrap PCA
# -----------------------------
for (b in 1:B) {
idx <- sample(1:nrow(X), replace = TRUE)
Xb <- X[idx, ]
pca_b <- prcomp(Xb, scale. = TRUE)
rot_b <- pca_b$rotation[, 1:pcs]
# Align signs
for (k in 1:pcs) {
if (cor(rot_b[, k], orig_rot[, k]) < 0) {
rot_b[, k] <- -rot_b[, k]
}
}
boot_loadings[, , b] <- rot_b
}
# ----------------------------- Summary statistics
# -----------------------------
abs_boot <- abs(boot_loadings)
mean_loading <- apply(abs_boot, c(1, 2), mean)
ci_lower <- apply(abs_boot, c(1, 2), quantile, 0.025)
ci_upper <- apply(abs_boot, c(1, 2), quantile, 0.975)
# ----------------------------- Stability frequency
# -----------------------------
top_freq <- matrix(0, nrow = p, ncol = pcs)
for (b in 1:B) {
for (k in 1:pcs) {
sq <- boot_loadings[, k, b]^2
ord <- order(sq, decreasing = TRUE)
cumprop <- cumsum(sq[ord])/sum(sq)
selected <- ord[cumprop <= cum_threshold]
selected <- c(selected, ord[min(which(cumprop >= cum_threshold))])
top_freq[selected, k] <- top_freq[selected, k] + 1
}
}
top_freq <- top_freq/B
stable <- (top_freq >= freq_threshold) & (ci_lower > 0 | ci_upper <
0)
# ----------------------------- Return tidy dataframe
# -----------------------------
data.frame(variable = rep(rownames(pca$rotation), pcs), ind = rep(colnames(pca$rotation)[1:pcs],
each = p), mean_loading = as.vector(mean_loading), ci_lower = as.vector(ci_lower),
ci_upper = as.vector(ci_upper), freq = as.vector(top_freq),
stable = as.vector(stable))
}
# higliht significant rows
highlight <- function(x, estimate_col = "Estimate", lower_col = "Q2.5",
upper_col = "Q97.5", strong_fill = "#E8602DFF", moderate_fill = "#FAC127FF",
weak_fill = "#FCFFA4FF", alpha = 0.5, digits = 3) {
## -------------------------------------------------- Row
## groups --------------------------------------------------
strong_rows <- which(x$pd > 0.95)
moderate_rows <- which(x$pd > 0.9 & x$pd <= 0.95)
weak_rows <- which(x$pd > 0.8 & x$pd <= 0.9)
## -------------------------------------------------- Build
## kable --------------------------------------------------
x_kbl <- kableExtra::kbl(x, row.names = TRUE, escape = FALSE,
format = "html", digits = digits)
## -------------------------------------------------- Apply
## row highlighting
## --------------------------------------------------
if (length(strong_rows) > 0) {
x_kbl <- kableExtra::row_spec(x_kbl, row = strong_rows, background = grDevices::adjustcolor(strong_fill,
alpha.f = alpha))
}
if (length(moderate_rows) > 0) {
x_kbl <- kableExtra::row_spec(x_kbl, row = moderate_rows,
background = grDevices::adjustcolor(moderate_fill, alpha.f = alpha))
}
if (length(weak_rows) > 0) {
x_kbl <- kableExtra::row_spec(x_kbl, row = weak_rows, background = grDevices::adjustcolor(weak_fill,
alpha.f = alpha))
}
## --------------------------------------------------
## Styling
## --------------------------------------------------
x_kbl <- kableExtra::kable_styling(x_kbl, bootstrap_options = c("striped",
"hover", "condensed", "responsive"), full_width = FALSE, font_size = 12)
return(x_kbl)
}1 Simple songs
1.1 Song level
1.1.1 Read and prepare data
Code
simple_songs <- read.csv("./data/raw/Sporophila song data_song-level plus MCP MST and PCA_simple songs.csv")
individual_coord <- read.csv("./data/raw/Coordenadas_individuos_simple_songs.csv")
# assign coordinates to each individual
simple_songs_lat_long <- simple_songs |>
left_join(individual_coord, by = "Individual")
# names(simple_songs_lat_long)1.1.2 Select variables and filter populations
Code
# select variables that are not highly correlated
simple_song_dat <- simple_songs_lat_long[, c("Individual", "Lat", "Lon", "species", "Population", "location",
"meanpeakf", "num.elms", "elm.duration", "freq.range",
"song.duration", "song.rate", "gap.duration",
"freq.range.Min5toMax95", "mst")]
# remove individuals from underrepresented populations
simple_song_dat <- filter(simple_song_dat, Population != "Pal_ER")
simple_song_dat <- filter(simple_song_dat, Population != "NA")
simple_song_dat$Individual <- factor(simple_song_dat$Individual)
simple_song_dat$species <- factor(simple_song_dat$species)
simple_song_dat$Population <- factor(simple_song_dat$Population)
simple_song_dat$location <- factor(simple_song_dat$location)1.1.3 Principal Component Analysis
We used the first 5 principal components (PCs) for subsequent analyses, which together explained 93% of the variance in the data.
Code
# run PCA on acoustic variables
pca <- prcomp(simple_song_dat[, c("meanpeakf", "num.elms", "elm.duration", "freq.range",
"song.duration", "song.rate", "gap.duration",
"freq.range.Min5toMax95", "mst")], center = TRUE, scale. = TRUE)
## Run PCA Inspect variance explained summary(pca)
# plot rotation values by PC
pca_rot <- as.data.frame(pca$rotation[, 1:5])
pca_var <- round(summary(pca)$importance[2, ] * 100)
pca_rot_stck <- stack(pca_rot)
pca_rot_stck$variable <- rownames(pca_rot)
pca_rot_stck$values[pca_rot_stck$ind == "PC1"] <- pca_rot_stck$values[pca_rot_stck$ind ==
"PC1"]
pca_rot_stck$Sign <- ifelse(pca_rot_stck$values > 0, "Positive", "Negative")
pca_rot_stck$rotation <- abs(pca_rot_stck$values)
pca_rot_stck$ind_var <- paste0(pca_rot_stck$ind, " (", sapply(pca_rot_stck$ind,
function(x) pca_var[names(pca_var) == x]), "%)")
pca_rot_stck$top_vars <- ave(abs(pca_rot_stck$values), pca_rot_stck$ind,
FUN = function(x) {
# Order decreasing
ord <- order(x, decreasing = TRUE)
x_sorted <- x[ord]
# Cumulative proportion
cumprop <- cumsum(x_sorted)/sum(x_sorted)
selected_sorted <- cumprop <= 0.5 # variables with 50% of contribution
selected_sorted[which(cumprop >= 0.5)[1]] <- TRUE
# Return logical vector in original order
selected <- logical(length(x))
selected[ord] <- selected_sorted
selected
})
# Create facet-specific variable
pca_rot_stck$var_facet <- paste(pca_rot_stck$ind, pca_rot_stck$variable,
sep = "_")
# Reorder within each ind by rotation (largest at top after
# coord_flip)
pca_rot_stck <- do.call(rbind, lapply(split(pca_rot_stck, pca_rot_stck$ind),
function(df) {
df$var_facet <- factor(df$var_facet, levels = df$var_facet[order(df$rotation)])
df
}))
# add which variables are stable
stable_df <- get_stable_loadings(pca, pcs = 5, B = 1000, cum_threshold = 0.5,
freq_threshold = 0.5, seed = 123)
pca_rot_stck <- merge(pca_rot_stck, stable_df, by = c("variable",
"ind"), all.x = TRUE)
pca_rot_stck$top_vars <- ifelse(pca_rot_stck$stable, 1, 1)
# Build colored labels per row
# Colored labels
pca_rot_stck$label_col <- ifelse(pca_rot_stck$stable < 1.1, paste0("<span style='color:black;'>",
pca_rot_stck$variable, "</span>"), paste0("<span style='color:gray50;'>",
pca_rot_stck$variable, "</span>"))
# Named vector for labels
label_vec <- setNames(pca_rot_stck$label_col, pca_rot_stck$var_facet)
# absolute CI
pca_rot_stck$ci_low_plot <- abs(pca_rot_stck$ci_lower)
pca_rot_stck$ci_high_plot <- abs(pca_rot_stck$ci_upper)
# Plot
ggplot(pca_rot_stck, aes(x = var_facet, y = rotation, fill = Sign,
alpha = as.factor(top_vars))) + geom_col() + coord_flip() + scale_alpha_manual(values = pca_rot_stck$top_vars,
guide = NULL) + scale_x_discrete(labels = label_vec, name = "Variable") +
labs(x = "Rotation") + scale_fill_viridis_d(alpha = 0.7, begin = 0.2,
end = 0.8) + facet_wrap(~ind_var, scales = "free_y", nrow = 3) + theme_classic() +
theme(axis.text.y = element_markdown())
Code
# bind PCA scores with metadata
pca_scores <- cbind(
pca$x,
simple_song_dat[,c("Population", "species", "location", "Individual", "Lat", "Lon")]
)
summ_pca <- summary(pca)
# select the PCs explaining at least 90%
pcs <- grep("^PC", names(pca_scores), value = TRUE)[1:min(which(summ_pca$importance[3,] > 0.9))]
pca_scores <- pca_scores |>
mutate(
across(all_of(pcs), ~ as.numeric(.x)),
Lat = as.numeric(Lat),
Lon = as.numeric(Lon),
species = as.character(species),
location = tryCatch(
iconv(location, from = "", to = "UTF-8"),
error = function(e) location
)
)1.1.4 Acoustic and geographic distances
1.1.4.1 Build pairwise distance dataset
Code
df_clean <- pca_scores |>
mutate(
across(all_of(pcs), as.numeric),
Lat = as.numeric(Lat),
Lon = as.numeric(Lon),
species = as.character(species),
Individual = as.character(Individual)
) |>
filter(complete.cases(across(all_of(c(pcs, "Lat", "Lon", "species", "Individual")))))
# table(df_clean$Population, df_clean$species)
# assign song IDs
df_clean$song_id <- 1
for (i in 2:nrow(df_clean)) {
if (df_clean$Individual[i] == df_clean$Individual[i - 1] &&
df_clean$species[i] == df_clean$species[i - 1]) {
df_clean$song_id[i] <- df_clean$song_id[i - 1]
} else {
df_clean$song_id[i] <- df_clean$song_id[i - 1] + 1
}
}
df_clean$song_id <- paste(
df_clean$species,
sapply(strsplit(as.character(df_clean$Population), "_"), "[[", 2),
df_clean$Individual,
df_clean$song_id,
sep = "-"
)1.1.4.2 Convert to long table
Code
# acoustic distance between songs (Euclidean multivariate, unscaled PCs)
dist_acoustic <- dist(
df_clean[, pcs],
method = "euclidean"
)
dist_acoustic_mat <- as.matrix(dist_acoustic)
# geographic distance (Haversine, km) between songs
coords <- as.matrix(df_clean[, c("Lon", "Lat")]) # Lon first, Lat second
D_geo_km <- geosphere::distm(coords, fun = distHaversine) / 1000
dist_geo <- as.dist(D_geo_km)
dist_geo_mat <- as.matrix(dist_geo)
rownames(dist_acoustic_mat) <- colnames(dist_acoustic_mat) <- df_clean$song_id
rownames(dist_geo_mat) <- colnames(dist_geo_mat) <- df_clean$song_id
# use only upper triangle
idx <- which(upper.tri(dist_acoustic_mat), arr.ind = TRUE)
dist_acoustic_long <- data.frame(
id1 = rownames(dist_acoustic_mat)[idx[, 1]],
id2 = colnames(dist_acoustic_mat)[idx[, 2]],
acoustic_distance = dist_acoustic_mat[idx],
geo_distance = dist_geo_mat[idx]
)
# parse species, population, and individual from song IDs
parts1 <- strsplit(as.character(dist_acoustic_long$id1), "-")
parts2 <- strsplit(as.character(dist_acoustic_long$id2), "-")
dist_acoustic_long$species1 <- sapply(parts1, `[`, 1)
dist_acoustic_long$population1 <- sapply(parts1, `[`, 2)
dist_acoustic_long$individual1 <- sapply(parts1, `[`, 3)
dist_acoustic_long$species2 <- sapply(parts2, `[`, 1)
dist_acoustic_long$population2 <- sapply(parts2, `[`, 2)
dist_acoustic_long$individual2 <- sapply(parts2, `[`, 3)
# keep only between-species comparisons
dist_acoustic_long <- dist_acoustic_long[
dist_acoustic_long$species1 != dist_acoustic_long$species2, ]
# sympatry flag: 1 if same population, 0 otherwise
dist_acoustic_long$sympatry <- as.integer(
dist_acoustic_long$population1 == dist_acoustic_long$population2
)
# canonical species pair label (sorted alphabetically)
dist_acoustic_long$species_pair <- apply(
dist_acoustic_long[, c("species1", "species2")],
1,
function(x) paste(sort(x), collapse = "_")
)
dist_acoustic_long$individual1 <- factor(dist_acoustic_long$individual1)
dist_acoustic_long$individual2 <- factor(dist_acoustic_long$individual2)
dist_acoustic_long$population1 <- factor(dist_acoustic_long$population1)
dist_acoustic_long$population2 <- factor(dist_acoustic_long$population2)
dist_acoustic_long$species_pair <- factor(dist_acoustic_long$species_pair)1.1.5 Statistical analysis
To evaluate whether acoustic divergence between heterospecific songs differed between sympatric and allopatric population comparisons, we fitted a series of Bayesian mixed-effects models while accounting for the non-independence inherent to pairwise distance data.
The global model was specified as:
[ ]
where:
- acoustic distance(_{ij}) is the Euclidean distance between songs (i) and (j) in multivariate acoustic space.
- sympatry(_{ij}) is a binary predictor indicating whether the two songs originated from populations occurring in sympatry (1) or allopatry (0).
- geographic distance(_{ij}) is the geographic distance separating the two populations being compared.
- species pair is a random intercept accounting for baseline differences in acoustic divergence among heterospecific species combinations.
- mm(population(_1), population(_2)) is a multi-membership random effect accounting for repeated use of the same populations across pairwise comparisons.
- mm(individual(_1), individual(_2)) is a multi-membership random effect accounting for repeated use of the same individuals across pairwise comparisons.
Model specifications:
- Models were fitted in a Bayesian framework using the brms package.
- The response variable was pairwise Euclidean acoustic distance and was modeled using a Gaussian error distribution.
- Only heterospecific comparisons were included.
- Analyses were restricted to species pairs for which both sympatric and allopatric population comparisons were available. This restriction ensured that sympatry effects were estimated within the same species-pair contrasts rather than being confounded by species pairs occurring exclusively in sympatry or exclusively in allopatry.
- Species-pair identity was included as a random intercept to account for inherent differences in acoustic divergence among species combinations.
- Population identity was modeled using a multi-membership random effect because each pairwise comparison simultaneously involves two populations, and each population contributes to multiple pairwise distances.
- Individual identity was modeled using a multi-membership random effect because each pairwise comparison simultaneously involves two individuals, and each individual contributes to multiple pairwise distances.
- Only unique pairwise comparisons were retained and self-comparisons were excluded.
- Weakly informative priors were specified for all parameters.
- Models were fitted using Hamiltonian Monte Carlo as implemented in Stan through the cmdstanr backend.
To assess the relative importance of sympatry and geographic distance, three competing models were fitted: (i) a sympatry-only model, (ii) a geographic-distance-only model, and (iii) a model including both predictors. Model performance was compared using approximate leave-one-out cross-validation (LOO). The coefficient for sympatry in the global model therefore represents differences in acoustic divergence between sympatric and allopatric population comparisons after accounting for geographic distance and the hierarchical structure of the data.
1.1.5.1 Model fitting
Prepare data for modeling by restricting to species pairs with both sympatric and allopatric comparisons and creating appropriate random effect structures.
Code
# restric to species pairs that have both sympatric and allopatric populations
tab <- table(dist_acoustic_long$species_pair,
dist_acoustic_long$sympatry)
colnames(tab) <- c("allopatric", "sympatric")
print(tab)| allopatric | sympatric |
|---|---|
| 1350 | 2088 |
| 2250 | 3480 |
| 30216 | 4164 |
| 0 | 540 |
| 3240 | 0 |
| 5400 | 0 |
Code
keep_pairs <- rownames(tab)[
tab[, "allopatric"] > 0 &
tab[, "sympatric"] > 0
]Only three species pairs were kept: Hypoxantha_Iberaensis, Hypoxantha_Palustris, Hypoxantha_Ruficollis
Code
sympatric_pairs <- dist_acoustic_long[
dist_acoustic_long$species_pair %in%
keep_pairs,
]
# create species-specific population IDs
sympatric_pairs$pop1 <- interaction(
sympatric_pairs$species1,
sympatric_pairs$population1,
drop = TRUE
)
sympatric_pairs$pop2 <- interaction(
sympatric_pairs$species2,
sympatric_pairs$population2,
drop = TRUE
)
# make species_pair a factor (fixed effect)
sympatric_pairs$species_pair <- factor(sympatric_pairs$species_pair)Code
# weak priors
priors <- c(
# Intercept
prior(normal(0, 5), class = "Intercept"),
# Fixed effect of sympatry
prior(normal(0, 2), class = "b"),
# Random-effect SDs
prior(exponential(1), class = "sd"),
# Residual SD
prior(exponential(1), class = "sigma")
)
sympatry_geo_model <- brm(
acoustic_distance ~ sympatry +
scale(geo_distance) +
(1 | species_pair) +
(1 | mm(pop1, pop2)) +
(1 | mm(individual1, individual2)),
data = sympatric_pairs,
family = gaussian(),
cores = 4,
chains = 4,
prior = priors,
iter = 10000,
backend = "cmdstanr",
threads = threading(8),
control = list(adapt_delta = 0.95, max_treedepth = 15),
file = "./data/processed/acoustic_distance_sympatry_geographic_distance_fit",
file_refit = "on_change"
)
sympatry_geo_model <- add_criterion(
sympatry_geo_model,
criterion = "loo"
)
sympatry_model <- brm(
acoustic_distance ~ sympatry +
(1 | species_pair) +
(1 | mm(pop1, pop2)) +
(1 | mm(individual1, individual2)),
data = sympatric_pairs,
family = gaussian(),
cores = 4,
chains = 4,
prior = priors,
iter = 10000,
backend = "cmdstanr",
threads = threading(8),
control = list(adapt_delta = 0.95, max_treedepth = 15),
file = "./data/processed/acoustic_distance_sympatry_fit",
file_refit = "on_change"
)
sympatry_model <- add_criterion(
sympatry_model,
criterion = "loo"
)
geo_model <- brm(
acoustic_distance ~ scale(geo_distance) +
(1 | species_pair) +
(1 | mm(pop1, pop2)) +
(1 | mm(individual1, individual2)),
data = sympatric_pairs,
family = gaussian(),
cores = 4,
chains = 4,
prior = priors,
iter = 10000,
backend = "cmdstanr",
threads = threading(8),
control = list(adapt_delta = 0.95, max_treedepth = 15),
file = "./data/processed/acoustic_distance_geographic_distance_fit",
file_refit = "on_change"
)
geo_model <- add_criterion(
geo_model,
criterion = "loo"
)
# beepr::beep(3)1.1.5.2 Model comparison
Code
comparison <- loo_compare(
sympatry_geo_model,
sympatry_model,
geo_model
)
print(comparison[, c("elpd_diff", "se_diff")])| elpd_diff | se_diff |
|---|---|
| 0.0000 | 0.0000 |
| -2.7441 | 2.8513 |
| -18.6656 | 6.2121 |
1.1.6 Fit summary
Code
models <- list.files("./data/processed/", pattern = "fit", full.names = TRUE)
for (i in models){
extended_summary(
read.file = i,
highlight = TRUE,
trace.palette = viridis::mako,
remove.intercepts = TRUE,
print.name = FALSE
)
}| formula | family | priors | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | acoustic_distance ~ scale(geo_distance) + (1 | species_pair) + (1 | mm(pop1, pop2)) + (1 | mm(individual1, individual2)) | gaussian (identity) | b-normal(0, 2) Intercept-normal(0, 5) sd-exponential(1) sigma-exponential(1) | 10000 | 4 | 1 | 5000 | 346 (0.017%) | 0 | 30177.69 | 12252.41 | 1984467235 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_scalegeo_distance | -0.073 | -0.09 | -0.057 | 1 | 30177.69 | 12252.41 |
| formula | family | priors | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | acoustic_distance ~ sympatry + (1 | species_pair) + (1 | mm(pop1, pop2)) + (1 | mm(individual1, individual2)) | gaussian (identity) | b-normal(0, 2) Intercept-normal(0, 5) sd-exponential(1) sigma-exponential(1) | 10000 | 4 | 1 | 5000 | 196 (0.0098%) | 0 | 26697.54 | 12458.94 | 1786129554 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sympatry | 0.171 | 0.139 | 0.204 | 1 | 26697.54 | 12458.94 |
| formula | family | priors | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | acoustic_distance ~ sympatry + scale(geo_distance) + (1 | species_pair) + (1 | mm(pop1, pop2)) + (1 | mm(individual1, individual2)) | gaussian (identity) | b-normal(0, 2) Intercept-normal(0, 5) sd-exponential(1) sigma-exponential(1) | 10000 | 4 | 1 | 5000 | 291 (0.015%) | 0 | 24920 | 13041.9 | 1423353538 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sympatry | 0.297 | 0.205 | 0.387 | 1 | 25895.72 | 13041.90 |
| b_scalegeo_distance | 0.069 | 0.022 | 0.115 | 1 | 24920.00 | 13565.85 |
Summary
The model including both sympatry and geographic distance had the highest predictive performance and was therefore the best-supported model according to leave-one-out cross-validation (LOO).
The sympatry-only model performed nearly as well as the full model (ΔELPD = -2.7 ± 2.9), indicating that geographic distance contributed little additional predictive information once sympatry was accounted for.
The geographic-distance-only model performed substantially worse than the full model (ΔELPD = -18.7 ± 6.2), demonstrating that geographic distance alone was a poor predictor of acoustic divergence relative to models that included sympatry.
The large performance difference between the geographic-distance-only model and the models containing sympatry suggests that sympatry is a much stronger predictor of acoustic divergence than geographic distance.
Overall, these results indicate that variation in acoustic divergence among heterospecific populations is primarily associated with whether populations occur in sympatry rather than with the geographic distance separating them.
Estimates show that sympatric population pairs have, on average, acoustic distances greater than those of allopatric pairs
Estimates also show an increase in acoustic distance associated to increasing geographic distances
Although both predictors were positively associated with acoustic divergence, the effect of sympatry was approximately four times larger than the effect of geographic distance (0.30 vs. 0.07)
Session information
─ Session info ───────────────────────────────────────────────────────────────
setting value
version R version 4.5.2 (2025-10-31)
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-06-19
pandoc 3.2 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/x86_64/ (via rmarkdown)
quarto 1.7.31 @ /usr/local/bin/quarto
─ Packages ───────────────────────────────────────────────────────────────────
package * version date (UTC) lib source
abind 1.4-8 2024-09-12 [1] CRAN (R 4.5.0)
ape 5.8-1 2024-12-16 [1] CRAN (R 4.5.0)
arrayhelpers 1.1-0 2020-02-04 [1] CRAN (R 4.5.0)
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)
bridgesampling 1.2-1 2025-11-19 [1] CRAN (R 4.5.0)
brms * 2.23.0 2025-09-09 [1] CRAN (R 4.5.0)
brmsish * 1.0.0 2026-01-06 [1] CRANs (R 4.5.2)
Brobdingnag 1.2-9 2022-10-19 [1] CRAN (R 4.5.0)
cachem 1.1.0 2024-05-16 [1] CRAN (R 4.5.0)
checkmate 2.3.4 2026-02-03 [1] CRAN (R 4.5.2)
class 7.3-23 2025-01-01 [4] CRAN (R 4.4.2)
classInt 0.4-11 2025-01-08 [1] CRAN (R 4.5.0)
cli 3.6.6 2026-04-09 [1] CRAN (R 4.5.2)
cluster 2.1.8.2 2026-02-05 [4] CRAN (R 4.5.2)
cmdstanr 0.9.0 2025-08-21 [1] https://stan-dev.r-universe.dev (R 4.5.0)
coda 0.19-4.1 2024-01-31 [1] CRAN (R 4.5.0)
codetools 0.2-20 2024-03-31 [4] CRAN (R 4.5.0)
commonmark 1.9.5 2025-03-17 [3] CRAN (R 4.5.0)
cowplot 1.2.0 2025-07-07 [1] CRAN (R 4.5.0)
crayon 1.5.3 2024-06-20 [1] CRAN (R 4.5.0)
curl 7.1.0 2026-04-22 [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.0)
devtools 2.4.5 2022-10-11 [1] CRAN (R 4.5.0)
digest 0.6.39 2025-11-19 [1] CRAN (R 4.5.0)
distributional 0.5.0 2024-09-17 [1] CRAN (R 4.5.0)
dplyr * 1.1.4 2023-11-17 [1] CRAN (R 4.5.0)
e1071 1.7-17 2025-12-18 [1] CRAN (R 4.5.2)
ecodist * 2.1.3 2023-10-30 [1] CRAN (R 4.5.0)
ellipsis 0.3.2 2021-04-29 [3] CRAN (R 4.1.1)
emmeans 1.11.1 2025-05-04 [3] CRAN (R 4.5.0)
estimability 1.5.1 2024-05-12 [3] CRAN (R 4.5.0)
evaluate 1.0.5 2025-08-27 [1] CRAN (R 4.5.0)
farver 2.1.2 2024-05-13 [1] CRAN (R 4.5.0)
fastmap 1.2.0 2024-05-15 [1] CRAN (R 4.5.0)
forcats * 1.0.0 2023-01-29 [1] CRAN (R 4.5.0)
fs 2.1.0 2026-04-18 [1] CRAN (R 4.5.2)
generics 0.1.4 2025-05-09 [1] CRAN (R 4.5.0)
geosphere * 1.5-14 2021-10-13 [3] CRAN (R 4.1.1)
ggdist 3.3.3 2025-04-23 [1] CRAN (R 4.5.0)
ggplot2 * 4.0.3 2026-04-22 [1] CRAN (R 4.5.2)
ggtext * 0.1.2 2022-09-16 [1] RSPM (R 4.5.0)
glue 1.8.1 2026-04-17 [1] CRAN (R 4.5.2)
goftest 1.2-3 2021-10-07 [1] CRAN (R 4.5.2)
gridExtra 2.3 2017-09-09 [1] CRAN (R 4.5.0)
gridtext 0.1.5 2022-09-16 [1] RSPM (R 4.5.0)
gtable 0.3.6 2024-10-25 [1] CRAN (R 4.5.0)
hms 1.1.3 2023-03-21 [3] CRAN (R 4.5.0)
htmltools 0.5.9 2025-12-04 [1] CRAN (R 4.5.0)
htmlwidgets 1.6.4 2023-12-06 [1] RSPM (R 4.5.0)
httpuv 1.6.16 2025-04-16 [1] RSPM (R 4.5.0)
igraph 2.2.1 2025-10-27 [1] CRAN (R 4.5.0)
inline 0.3.21 2025-01-09 [1] CRAN (R 4.5.0)
jsonlite 2.0.0 2025-03-27 [1] CRAN (R 4.5.0)
kableExtra * 1.4.0 2024-01-24 [1] CRAN (R 4.5.0)
KernSmooth 2.23-26 2025-01-01 [4] CRAN (R 4.4.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.0)
later 1.4.2 2025-04-08 [1] RSPM (R 4.5.0)
lattice 0.22-9 2026-02-09 [4] CRAN (R 4.5.2)
lifecycle 1.0.5 2026-01-08 [1] CRAN (R 4.5.2)
litedown 0.7 2025-04-08 [1] CRAN (R 4.5.0)
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)
markdown 2.0 2025-03-23 [1] CRAN (R 4.5.0)
MASS 7.3-65 2025-02-28 [4] CRAN (R 4.4.3)
Matrix 1.7-4 2025-08-28 [4] CRAN (R 4.5.1)
matrixStats 1.5.0 2025-01-07 [1] CRAN (R 4.5.0)
memoise 2.0.1 2021-11-26 [1] CRAN (R 4.5.0)
mgcv 1.9-4 2025-11-07 [4] CRAN (R 4.5.2)
mime 0.13 2025-03-17 [1] CRAN (R 4.5.0)
miniUI 0.1.2 2025-04-17 [3] CRAN (R 4.5.0)
multcomp 1.4-28 2025-01-29 [3] CRAN (R 4.5.0)
mvtnorm 1.3-3 2025-01-10 [1] CRAN (R 4.5.0)
nicheROVER 1.1.2 2023-10-13 [1] CRAN (R 4.5.0)
nlme 3.1-168 2025-03-31 [4] CRAN (R 4.4.3)
otel 0.2.0 2025-08-29 [1] CRAN (R 4.5.2)
packrat 0.9.3 2025-06-16 [1] CRAN (R 4.5.2)
pbapply 1.7-4 2025-07-20 [1] CRAN (R 4.5.0)
permute * 0.9-10 2026-02-06 [1] CRAN (R 4.5.2)
PhenotypeSpace * 0.1.1 2026-02-25 [1] CRAN (R 4.5.2)
pillar 1.11.1 2025-09-17 [1] CRAN (R 4.5.0)
pkgbuild 1.4.8 2025-05-26 [1] CRAN (R 4.5.0)
pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.5.0)
pkgload 1.4.1 2025-09-23 [1] CRAN (R 4.5.0)
plyr 1.8.9 2023-10-02 [1] CRAN (R 4.5.0)
polyclip 1.10-7 2024-07-23 [1] CRAN (R 4.5.0)
posterior 1.6.1 2025-02-27 [1] CRAN (R 4.5.0)
processx 3.8.6 2025-02-21 [1] CRAN (R 4.5.0)
profvis 0.4.0 2024-09-20 [1] CRAN (R 4.5.0)
promises 1.3.3 2025-05-29 [1] RSPM (R 4.5.0)
proxy 0.4-29 2025-12-29 [1] CRAN (R 4.5.2)
ps 1.9.1 2025-04-12 [1] CRAN (R 4.5.0)
purrr * 1.2.0 2025-11-04 [1] CRAN (R 4.5.0)
QuickJSR 1.8.1 2025-09-20 [1] CRAN (R 4.5.0)
R6 2.6.1 2025-02-15 [1] CRAN (R 4.5.0)
raster 3.6-32 2025-03-28 [1] CRAN (R 4.5.0)
RColorBrewer 1.1-3 2022-04-03 [1] CRAN (R 4.5.0)
Rcpp * 1.1.1-1.1 2026-04-24 [1] CRAN (R 4.5.2)
RcppParallel 5.1.11-1 2025-08-27 [1] CRAN (R 4.5.0)
readr * 2.1.5 2024-01-10 [1] CRAN (R 4.5.0)
remotes 2.5.0 2024-03-17 [1] CRAN (R 4.5.0)
reshape2 1.4.5 2025-11-12 [1] CRAN (R 4.5.0)
rlang 1.2.0 2026-04-06 [1] CRAN (R 4.5.2)
rmarkdown 2.31 2026-03-26 [1] CRAN (R 4.5.2)
rstan 2.32.7 2025-03-10 [1] CRAN (R 4.5.0)
rstantools 2.5.0 2025-09-01 [1] CRAN (R 4.5.0)
rstudioapi 0.17.1 2024-10-22 [1] CRAN (R 4.5.0)
S7 0.2.2 2026-04-22 [1] CRAN (R 4.5.2)
sandwich 3.1-1 2024-09-15 [3] CRAN (R 4.5.0)
scales 1.4.0 2025-04-24 [1] CRAN (R 4.5.0)
sessioninfo 1.2.3 2025-02-05 [3] CRAN (R 4.5.0)
sf 1.1-0 2026-02-24 [1] CRAN (R 4.5.2)
shiny 1.10.0 2024-12-14 [1] CRAN (R 4.5.0)
sketchy 1.0.7 2026-01-28 [1] CRANs (R 4.5.2)
sp 2.2-1 2026-02-13 [1] CRAN (R 4.5.2)
spatstat.data 3.1-9 2025-10-18 [1] CRAN (R 4.5.2)
spatstat.explore 3.7-0 2026-01-22 [1] CRAN (R 4.5.2)
spatstat.geom 3.7-0 2026-01-20 [1] CRAN (R 4.5.2)
spatstat.random 3.4-4 2026-01-21 [1] CRAN (R 4.5.2)
spatstat.sparse 3.1-0 2024-06-21 [1] CRAN (R 4.5.0)
spatstat.univar 3.1-6 2026-01-17 [1] CRAN (R 4.5.2)
spatstat.utils 3.2-1 2026-01-10 [1] CRAN (R 4.5.2)
StanHeaders 2.32.10 2024-07-15 [1] CRAN (R 4.5.0)
stringi 1.8.7 2025-03-27 [1] CRAN (R 4.5.0)
stringr * 1.6.0 2025-11-04 [1] CRAN (R 4.5.0)
survival 3.8-6 2026-01-16 [4] CRAN (R 4.5.2)
svglite 2.2.2 2025-10-21 [1] CRAN (R 4.5.0)
svUnit 1.0.8 2025-08-26 [1] CRAN (R 4.5.0)
systemfonts 1.3.1 2025-10-01 [1] CRAN (R 4.5.0)
tensor 1.5.1 2025-06-17 [1] CRAN (R 4.5.0)
tensorA 0.36.2.1 2023-12-13 [1] CRAN (R 4.5.0)
terra 1.9-11 2026-03-26 [1] CRAN (R 4.5.2)
textshaping 1.0.4 2025-10-10 [1] CRAN (R 4.5.0)
TH.data 1.1-3 2025-01-17 [3] CRAN (R 4.5.0)
tibble * 3.3.0 2025-06-08 [1] RSPM (R 4.5.0)
tidybayes 3.0.7 2024-09-15 [1] CRAN (R 4.5.0)
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.0)
tidyverse * 2.0.0 2023-02-22 [1] RSPM (R 4.5.0)
timechange 0.4.0 2026-01-29 [1] CRAN (R 4.5.2)
tzdb 0.5.0 2025-03-15 [1] CRAN (R 4.5.0)
units 1.0-1 2026-03-11 [1] CRAN (R 4.5.2)
urlchecker 1.0.1 2021-11-30 [1] CRAN (R 4.5.0)
usethis 3.1.0 2024-11-26 [3] CRAN (R 4.5.0)
V8 6.0.4 2025-06-04 [1] RSPM (R 4.5.0)
vctrs 0.7.3 2026-04-11 [1] CRAN (R 4.5.2)
vegan * 2.7-2 2025-10-08 [1] CRAN (R 4.5.2)
viridis 0.6.5 2024-01-29 [1] CRAN (R 4.5.0)
viridisLite 0.4.3 2026-02-04 [1] CRAN (R 4.5.2)
withr 3.0.2 2024-10-28 [1] CRAN (R 4.5.0)
xaringanExtra 0.8.0 2024-05-19 [1] CRAN (R 4.5.0)
xfun 0.57 2026-03-20 [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-14 2025-04-10 [3] CRAN (R 4.5.0)
[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.
──────────────────────────────────────────────────────────────────────────────