Code
# options to customize chunk outputs
knitr::opts_chunk$set(
message = FALSE
)Statistical analysis
Vocal character displacement in southern capuchinos
Purpose
- Evaluate the role of sympatry in the vocal divergence 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",
github = "maRce10/brmsish",
"kableExtra",
"ggplot2",
"ggtext",
github = "stan-dev/cmdstanr"
)
# install/load packages
sketchy::load_packages(packages = packages)
options(knitr.kable.NA = '', brms.file_refit = "never")
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))
}
# number of PCs to retain using the broken-stick criterion (Frontier 1976; Jackson 1993):
# a PC is retained while its observed proportion of variance exceeds the proportion
# expected under a random division ("broken stick") of total variance among the same
# number of components
n_pcs_broken_stick <- function(pca) {
# observed proportion of variance explained by each PC
obs_var <- summary(pca)$importance[2, ]
p <- length(obs_var)
# expected proportion of variance under the broken-stick null model
bstick_expected <- sapply(seq_len(p), function(k) sum(1 / (k:p)) / p)
# first PC (if any) at which the observed variance no longer exceeds
# the broken-stick expectation
below <- which(obs_var <= bstick_expected)
n_keep <- if (length(below) == 0) p else below[1] - 1
# always retain at least one PC
max(n_keep, 1)
}
# 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)
}
plot_brms_heatmap <- function(
model_files,
remove_intercepts = TRUE
) {
# models may still be fitting in the background, so some (or all) of
# the expected files can be missing, and a file that is still being
# written by brms can be present but not yet readable - skip either
# case instead of erroring
if(length(model_files) == 0) {
message("plot_brms_heatmap: no model files found yet - skipping.")
return(invisible(NULL))
}
effects_df <- data.frame()
for(i in seq_along(model_files)) {
fit <- tryCatch(
readRDS(model_files[i]),
error = function(e) NULL
)
if(is.null(fit)) {
message(
"plot_brms_heatmap: could not read '", model_files[i],
"' (likely still being written) - skipping."
)
next
}
fe <- as.data.frame(fixef(fit))
fe$predictor <- rownames(fe)
if(remove_intercepts)
fe <- fe[fe$predictor != "Intercept", ]
response <- deparse(fit$formula$formula[[2]])
fe$response <- response
effects_df <- rbind(
effects_df,
fe
)
}
if(nrow(effects_df) == 0) {
message("plot_brms_heatmap: no readable model results yet - skipping.")
return(invisible(NULL))
}
rownames(effects_df) <- NULL
# significance
effects_df$sig <- with(
effects_df,
Q2.5 * Q97.5 > 0
)
# clean names
effects_df$predictor <- gsub("^scale\\(", "", effects_df$predictor)
effects_df$predictor <- gsub("\\)$", "", effects_df$predictor)
effects_df$predictor <- gsub("^mo", "", effects_df$predictor)
effects_df$predictor <- gsub("^mi", "", effects_df$predictor)
effects_df$predictor <- gsub("_sc$", "", effects_df$predictor)
effects_df$response <- gsub("^mi", "", effects_df$response)
effects_df$response <- gsub("_sc$", "", effects_df$response)
# effects_df$predictor <- ifelse(grepl("sympatry", effects_df$predictor), "sympatry", effects_df$predictor)
# average duplicated cells if present
effects_df <- aggregate(
cbind(
Estimate,
sig
) ~ predictor + response,
data = effects_df,
FUN = mean
)
effects_df$sig <- effects_df$sig > 0.5
# complete combinations
all_combos <- expand.grid(
predictor = unique(effects_df$predictor),
response = unique(effects_df$response)
)
plot_df <- merge(
all_combos,
effects_df,
by = c("predictor", "response"),
all.x = TRUE
)
plot_df$response <- gsub("_distance$", "", plot_df$response)
plot_df$predictor <- gsub("scalegeo_", "Geographic\n", plot_df$predictor)
plot_df$predictor <- gsub("geo_distance", "Geographic\ndistance", plot_df$predictor)
plot_df$predictor <- gsub("sympatry1", "Sympatry", plot_df$predictor)
lim <- max(abs(plot_df$Estimate), na.rm = TRUE)
ggplot(
plot_df,
aes(
predictor,
response
)
) +
geom_tile(
fill = "grey90",
colour = "white"
) +
geom_tile(
data = plot_df[!is.na(plot_df$Estimate), ],
aes(fill = Estimate),
colour = "white"
) +
geom_text(
data = plot_df[!is.na(plot_df$Estimate), ],
aes(
label = sprintf("%.2f", Estimate),
colour = sig
),
fontface = "bold",
size = 3
) +
# scale_fill_gradient2(
# low = rep("#403B78", 2),
# mid = "white",
# high = rep("#DEF5E5", rep = 2),
# midpoint = 0,
# name = "Estimate"
# ) +
#
scale_fill_gradient2(
low = "#403B78",
mid = "white",
high = "#A0DFB9CC",
midpoint = 0,
limits = c(-lim, lim),
oob = scales::squish,
name = "Estimate"
) +
scale_color_manual(
values = c(
"TRUE" = "black",
"FALSE" = "grey70"
),
guide = "none"
) +
labs(
x = "Predictor",
y = "Response"
) +
theme_classic() +
theme(
axis.text.x = element_text(
angle = 45,
hjust = 1
)
)
}1 Geographic distance transform
Geographic distance between recording locations is log-transformed before entering the models, since under isolation-by-distance expectations for populations distributed across a two-dimensional landscape (rather than along a linear transect), divergence is expected to scale with log(geographic distance) rather than with raw distance (Rousset 1997). The additive offset needed to log-transform pairs with zero distance (i.e. same-population comparisons) is estimated from the data as half the smallest non-zero pairwise distance actually observed between recording locations, pooled across both song types, rather than from a distribution-dependent quantile - this ties the offset to the finest spatial resolution present in the sampling instead of the shape of the bulk distribution, and keeps it identical regardless of which subset of pairs (simple or complex songs) it is applied to. The log-transformed distances are then centered and scaled using the pooled mean and SD (again computed once, across both song types) and divided by 2 SD (Gelman 2008), so that a unit change in geo_distance_sc means the same thing - the same number of km, on the same log scale - in every model, and is directly comparable to the binary sympatry/same_population predictors.
Code
coord_simple <- read.csv("./data/raw/Coordenadas_individuos_simple_songs.csv", stringsAsFactors = FALSE)
coord_complex <- read.csv("./data/raw/Coordenadas_individuos_complex_songs.csv", stringsAsFactors = FALSE)
pooled_coords <- unique(rbind(
coord_simple[, c("Lat", "Lon")],
coord_complex[, c("Lat", "Lon")]
))
pooled_coords <- pooled_coords[complete.cases(pooled_coords), ]
# all pairwise Haversine distances (km) among pooled recording locations
# (both song types combined)
D_pooled_km <- geosphere::distm(
as.matrix(pooled_coords[, c("Lon", "Lat")]),
fun = geosphere::distHaversine
) / 1000
# offset = half the smallest non-zero pairwise distance actually observed
# across the whole study, i.e. the finest spatial resolution in the
# sampling, rather than a quantile of the bulk distribution
nonzero_pooled_d <- D_pooled_km[upper.tri(D_pooled_km)]
nonzero_pooled_d <- nonzero_pooled_d[nonzero_pooled_d > 0]
geo_const <- min(nonzero_pooled_d) / 2
# pooled mean/SD of the log-transformed pooled distances, so "1 SD" of
# geo_distance_sc represents the same number of km in every model
log_pooled_d <- log(D_pooled_km[upper.tri(D_pooled_km)] + geo_const)
geo_log_mean <- mean(log_pooled_d)
geo_log_sd <- sd(log_pooled_d)
# shared transform applied to each song type's own geo_distance column:
# log distance, centered/scaled by the pooled mean & SD, divided by 2
# (Gelman 2008) so the coefficient is comparable to the binary
# sympatry / same_population predictors
transform_geo_distance <- function(geo_distance) {
((log(geo_distance + geo_const) - geo_log_mean) / geo_log_sd) / 2
}2 Simple songs
2.1 Prepare element level data
Code
simple_elm <- read.csv("./data/raw/Sporophila song data_elmt-level plus PCA_simple.csv")
# remove all columns that start with "PC"
simple_elm <- simple_elm[, !grepl("^PC", names(simple_elm))]
simple_elm$Population <- simple_elm$Lat <- simple_elm$Lon <- NULL
individual_coord <- read.csv("./data/raw/Coordenadas_individuos_simple_songs.csv")
# assign coordinates to each individual
simple_elm_lat_long <- simple_elm |>
left_join(individual_coord, by = "Individual")2.1.1 Select variables and filter populations
Code
vars <- c(
"Q1.Time..s.",
"Q3.Time..s.",
"Time.5...s.",
"Time.95...s.",
"Delta.Time..s.",
"Dur.90...s.",
"IQR.Dur..s.",
"Peak.Time..s.",
"Center.Time..s.",
"Q1.Freq..Hz.",
"Q3.Freq..Hz.",
"Center.Freq..Hz.",
"Freq.5...Hz.",
"Freq.95...Hz.",
"Delta.Freq..Hz.",
"IQR.BW..Hz.",
"BW.90...Hz.",
"Max.Freq..Hz.",
"Peak.Freq..Hz.",
"meanfreq",
"sd",
"freq.median",
"freq.Q25",
"freq.Q75",
"freq.IQR",
"time.median",
"time.Q25",
"time.Q75",
"time.IQR",
"skew",
"kurt",
"sp.ent",
"time.ent",
"entropy",
"sfm",
"meandom",
"mindom",
"maxdom",
"dfrange",
"modindx",
"startdom",
"enddom",
"dfslope",
"meanpeakf"
)
# select variables that are not highly correlated
simple_elm_dat <- simple_elm_lat_long[, c("Individual", "Lat", "Lon", "species", "Population", "location", "song", vars)]
# remove individuals from underrepresented populations
simple_elm_dat <- filter(simple_elm_dat, Population != "Pal_ER")
simple_elm_dat <- filter(simple_elm_dat, Population != "NA")
simple_elm_dat$Individual <- factor(simple_elm_dat$Individual)
simple_elm_dat$species <- factor(simple_elm_dat$species)
simple_elm_dat$Population <- factor(simple_elm_dat$Population)
simple_elm_dat$location <- factor(simple_elm_dat$location)
simple_elm_dat$song <- factor(simple_elm_dat$song)
simple_elm_dat <- simple_elm_dat[complete.cases(simple_elm_dat), ]2.1.2 Principal Component Analysis
Code
# run PCA on acoustic variables
pca <- prcomp(
simple_elm_dat[, vars],
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)We used the first 4 principal components (PCs), selected using the broken-stick criterion, for subsequent analyses, which together explained 78.8% of the variance in the data.
Code
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
elm_pca_scores <- cbind(
pca$x,
simple_elm_dat[, c(
"Population",
"species",
"location",
"Individual",
"song",
"Lat",
"Lon")]
)
# select the PCs retained by the broken-stick criterion
pcs_elm <- grep("^PC", names(elm_pca_scores), value = TRUE)[seq_len(n_pcs_broken_stick(pca))]
elm_pca_scores <- elm_pca_scores |>
mutate(
across(all_of(pcs_elm), ~ 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
)
)
df_clean_simple_elm <- elm_pca_scores |>
mutate(
across(all_of(pcs_elm), as.numeric),
Lat = as.numeric(Lat),
Lon = as.numeric(Lon),
species = as.character(species),
Individual = as.character(Individual),
song = as.character(song)
) |>
filter(complete.cases(across(all_of(c(pcs_elm, "Lat", "Lon", "species", "Individual", "song")))))2.2 Prepare Song level 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)2.2.1 Select variables and filter populations
The song level features used were: peak frequency, number of elements, song duration, song rate, gap duration, frequency range and element diversity (mst)
Element duration was excluded as it is an element level features
Code
# select variables that are not highly correlated
simple_song_dat <- simple_songs_lat_long[, c("Individual", "Lat", "Lon", "species", "Population", "location", "song",
"meanpeakf", "num.elms",
"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)
simple_song_dat$song <- factor(simple_song_dat$song)2.2.2 Principal Component Analysis
Code
# run PCA on acoustic variables
pca <- prcomp(
simple_song_dat[, c(
"meanpeakf",
"num.elms",
"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)We used the first 1 principal components (PCs), selected using the broken-stick criterion, for subsequent analyses, which together explained 47.6% of the variance in the data.
Code
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",
"song",
"Lat",
"Lon",
"meanpeakf",
"num.elms",
"song.duration",
"song.rate",
"gap.duration",
"freq.range.Min5toMax95",
"mst"
)]
)
# select the PCs retained by the broken-stick criterion
pcs <- grep("^PC", names(pca_scores), value = TRUE)[seq_len(n_pcs_broken_stick(pca))]
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
)
)2.2.3 Acoustic and geographic distances
2.2.3.1 Build pairwise distance dataset
Code
df_clean_simple <- 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),
population = as.character(location),
song = as.character(song)
) |>
filter(complete.cases(across(all_of(c(pcs, "Lat", "Lon", "species", "Individual", "song")))))
# table(df_clean_simple$Population, df_clean_simple$species)
# assign song IDs
# df_clean_simple$song_id <- 1
#
# for (i in 2:nrow(df_clean_simple)) {
# if (df_clean_simple$Individual[i] == df_clean_simple$Individual[i - 1] &&
# df_clean_simple$species[i] == df_clean_simple$species[i - 1]) {
# df_clean_simple$song_id[i] <- df_clean_simple$song_id[i - 1]
# } else {
# df_clean_simple$song_id[i] <- df_clean_simple$song_id[i - 1] + 1
# }
# }
#
# df_clean_simple$song_id <- paste(
# df_clean_simple$species,
# sapply(strsplit(as.character(df_clean_simple$Population), "_"), "[[", 2),
# df_clean_simple$Individual,
# df_clean_simple$song_id,
# sep = "-"
# )2.2.3.2 Convert to long table
Code
# acoustic distance between songs (Euclidean multivariate, unscaled PCs)
agg_df_clean_simple_elm <- aggregate(
. ~ song, df_clean_simple_elm[, c(pcs_elm, "song")],
FUN = mean
)
dist_acoustic_mat <- as.matrix(dist(
df_clean_simple[, pcs],
method = "euclidean"
))
# and for each feature separately (unscaled)
dist_meanpeakf_mat <- as.matrix(dist(
df_clean_simple[, "meanpeakf"],
method = "euclidean"
))
dist_numelms_mat <- as.matrix(dist(
df_clean_simple[, "num.elms"],
method = "euclidean"
))
dist_songduration_mat <- as.matrix(dist(
df_clean_simple[, "song.duration"],
method = "euclidean"
))
dist_songrate_mat <- as.matrix(dist(
df_clean_simple[, "song.rate"],
method = "euclidean"
))
dist_gapduration_mat <- as.matrix(dist(
df_clean_simple[, "gap.duration"],
method = "euclidean"
))
dist_freqrange_mat <- as.matrix(dist(
df_clean_simple[, "freq.range.Min5toMax95"],
method = "euclidean"
))
dist_mst_mat <- as.matrix(dist(
df_clean_simple[, "mst"],
method = "euclidean"
))
dist_elm_mat <- as.matrix(dist(
agg_df_clean_simple_elm[, pcs_elm],
method = "euclidean"
))
# geographic distance (Haversine, km) between songs
coords <- as.matrix(df_clean_simple[, 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_simple$song
rownames(dist_geo_mat) <- colnames(dist_geo_mat) <- df_clean_simple$song
rownames(dist_elm_mat) <- colnames(dist_elm_mat) <- agg_df_clean_simple_elm$song
# 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],
meanpeakf_distance = dist_meanpeakf_mat[idx],
numelms_distance = dist_numelms_mat[idx],
songduration_distance = dist_songduration_mat[idx],
songrate_distance = dist_songrate_mat[idx],
gapduration_distance = dist_gapduration_mat[idx],
freqrange_distance = dist_freqrange_mat[idx],
mst_distance = dist_mst_mat[idx],
geo_distance = dist_geo_mat[idx]
)
dist_acoustic_long$elm_pca_distance <- sapply(seq_len(nrow(dist_acoustic_long)), function(i) {
dist_elm_mat[rownames(dist_elm_mat) == dist_acoustic_long$id1[i], colnames(dist_elm_mat) == dist_acoustic_long$id2[i]]
})
# 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(dist_acoustic_long$id1, function(x)
df_clean_simple$species[df_clean_simple$song == x])
dist_acoustic_long$population1 <- sapply(dist_acoustic_long$id1, function(x)
df_clean_simple$population[df_clean_simple$song == x])
dist_acoustic_long$individual1 <- sapply(dist_acoustic_long$id1, function(x)
df_clean_simple$Individual[df_clean_simple$song == x])
dist_acoustic_long$species2 <- sapply(dist_acoustic_long$id2, function(x)
df_clean_simple$species[df_clean_simple$song == x])
dist_acoustic_long$population2 <- sapply(dist_acoustic_long$id2, function(x)
df_clean_simple$population[df_clean_simple$song == x])
dist_acoustic_long$individual2 <- sapply(dist_acoustic_long$id2, function(x)
df_clean_simple$Individual[df_clean_simple$song == x])
# 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.factor(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 = "_")
)
# offset, centering and scaling are estimated once from the pooled data
# across both song types (see "Geographic distance transform" section)
# so this transform is identical between the simple- and complex-song
# models
dist_acoustic_long$log_geo_distance <- log(dist_acoustic_long$geo_distance + geo_const)
dist_acoustic_long$geo_distance_sc <- transform_geo_distance(dist_acoustic_long$geo_distance)
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)2.2.4 Distance and sympatry
Code
ggplot(dist_acoustic_long, aes(x = log_geo_distance, fill = as.factor(sympatry))) + geom_histogram(bins = 50,
position = "identity", alpha = 0.6) + scale_x_log10() + scale_fill_viridis_d(begin = 0.2,
end = 0.75, labels = c("Allopatric", "Sympatric")) + labs(x = "Geographic distance (km, log scale)",
y = "Number of pairs", fill = "") + theme_classic()Warning in transformation$transform(x): NaNs producedWarning in scale_x_log10(): log-10 transformation introduced infinite values.Warning: Removed 1260 rows containing non-finite outside the scale range
(`stat_bin()`).
2.2.5 Collinearity
Code
vif_mod <- lm(acoustic_distance ~ sympatry + geo_distance_sc, data = dist_acoustic_long)
vif_vals <- car::vif(vif_mod)
vif_df <- data.frame(term = names(vif_vals), vif = as.numeric(vif_vals))
## conventional rule-of-thumb thresholds: <5 low concern, 5-10
## moderate, >10 high. Some fields use stricter cutoffs (2.5 /
## 4); doesn't matter for your numbers specifically, since
## ~1.7-2.2 clears either convention.
vif_df$severity <- cut(vif_df$vif, breaks = c(0, 5, 10, Inf), labels = c("low",
"moderate", "high"))
cols <- viridis::mako(10)
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)
2.2.6 Statistical analysis
To evaluate whether acoustic divergence between heterospecific songs differed between sympatric and allopatric population comparisons, we fitted a Bayesian mixed-effects model while accounting for the non-independence inherent to pairwise distance data.
The model was specified as:
\[ \begin{split} \text{acoustic distance} &\sim \text{sympatry} + \text{geographic distance} \\ &\quad + (1 \mid \text{species pair}) \\ &\quad + (1 \mid \text{mm(population}_1,\text{population}_2)) \\ &\quad + (1 \mid \text{mm(individual}_1,\text{individual}_2)) \end{split} \]
where:
(_{ij}) is the Euclidean distance between songs (i) and (j) in multivariate acoustic space.
(_{ij}) is a binary predictor indicating whether the populations from which songs (i) and (j) were recorded occur in sympatry (1) or allopatry (0).
(_{ij}) is the geographic distance between the populations from which songs (i) and (j) were recorded, log-transformed and scaled (see Geographic distance transform).
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:
- The model was fitted in a Bayesian framework using the brms package, with sympatry and geographic distance specified together as predictors a priori; no alternative (sympatry-only or geographic-distance-only) models were fitted or compared.
- 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, regularizing priors were specified for all parameters.
- The model was fitted using Hamiltonian Monte Carlo as implemented in Stan through the cmdstanr backend.
The coefficient for sympatry in this model represents differences in acoustic divergence between sympatric and allopatric population comparisons, after accounting for geographic distance and the hierarchical structure of the data. The same model structure was fitted separately for the overall PCA-based acoustic distance and for each individual acoustic feature, so that trait-specific effects of sympatry and geographic distance could be examined alongside the overall pattern.
2.2.6.1 Model fitting
2.2.6.1.1 PCA-based acoustic distance
Prepare data for modeling by restricting to species pairs with both sympatric and allopatric comparisons and creating appropriate random effect structures.
Species by location:
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")
keep_pairs <- rownames(tab)[
tab[, "allopatric"] > 0 &
tab[, "sympatric"] > 0
]
# Extract all species-population combinations
sp_pop <- unique(
rbind(
data.frame(
species = dist_acoustic_long$species1,
population = dist_acoustic_long$population1
),
data.frame(
species = dist_acoustic_long$species2,
population = dist_acoustic_long$population2
)
)
)
# Presence/absence table
tab <- with(
sp_pop,
table(species, population)
)
# Convert counts to X / blank
tab[] <- ifelse(tab > 0, "\u2713", "")
presence_table <- as.data.frame.matrix(tab)
presence_table$species <- rownames(presence_table)
presence_table <- presence_table[
, c("species", setdiff(names(presence_table), "species"))
]
print(presence_table)| species | E_Ibera | Entre_Rios | Esperanza | Mar_Chiquita | Salta |
|---|---|---|---|---|---|
| Hypoxantha | ✓ | ✓ | ✓ | ||
| Iberaensis | ✓ | ||||
| Palustris | ✓ | ||||
| Ruficollis | ✓ | ✓ | ✓ | ✓ |
Species pairs by sympatry:
Code
# Species x population occurrence matrix
occ <- as.matrix(tab)
species <- rownames(occ)
# All pairwise species combinations
pairs <- combn(species, 2, simplify = FALSE)
results <- data.frame(
Pair = character(),
Sympatric = character(),
Allopatric = character(),
stringsAsFactors = FALSE
)
for(p in pairs) {
sp1 <- p[1]
sp2 <- p[2]
pops1 <- colnames(occ)[occ[sp1, ] != ""]
pops2 <- colnames(occ)[occ[sp2, ] != ""]
sympatric <- intersect(pops1, pops2)
if(length(sympatric) == 0) {
sympatric_txt <- "none"
} else {
sympatric_txt <- paste(sympatric, collapse = ", ")
}
allopatric <- setdiff(union(pops1, pops2), sympatric)
if(length(allopatric) == 0) {
allopatric_txt <- "none"
} else {
allopatric_txt <- paste(allopatric, collapse = ", ")
}
results <- rbind(
results,
data.frame(
Pair = paste(sp1, sp2, sep = "–"),
Sympatric = sympatric_txt,
Allopatric = allopatric_txt,
stringsAsFactors = FALSE
)
)
}
print(results)| Pair | Sympatric | Allopatric |
|---|---|---|
| Hypoxantha–Iberaensis | E_Ibera | Entre_Rios, Mar_Chiquita |
| Hypoxantha–Palustris | E_Ibera | Entre_Rios, Mar_Chiquita |
| Hypoxantha–Ruficollis | Entre_Rios, Mar_Chiquita | E_Ibera, Esperanza, Salta |
| Iberaensis–Palustris | E_Ibera | none |
| Iberaensis–Ruficollis | none | E_Ibera, Entre_Rios, Esperanza, Mar_Chiquita, Salta |
| Palustris–Ruficollis | none | E_Ibera, Entre_Rios, Esperanza, Mar_Chiquita, Salta |
Only three species pairs were kept: Hypoxantha_Iberaensis, Hypoxantha_Palustris, Hypoxantha_Ruficollis
2.2.6.1.1.1 Species pairs with both sympatric and allopatric populations
Code
sympatric_pairs_simple <- dist_acoustic_long[
dist_acoustic_long$species_pair %in%
keep_pairs,
]
# create species-specific population IDs
sympatric_pairs_simple$pop1 <- interaction(
sympatric_pairs_simple$species1,
sympatric_pairs_simple$population1,
drop = TRUE
)
sympatric_pairs_simple$pop2 <- interaction(
sympatric_pairs_simple$species2,
sympatric_pairs_simple$population2,
drop = TRUE
)
# make species_pair a factor (fixed effect)
sympatric_pairs_simple$species_pair <- factor(sympatric_pairs_simple$species_pair)
# scale acoustic_distance
sympatric_pairs_simple$acoustic_distance_sc <- scale(sympatric_pairs_simple$acoustic_distance)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")
)
#fit model
sympatry_geo_model_simple <- brm(
acoustic_distance_sc ~ sympatry +
geo_distance_sc +
(1 | species_pair) +
(1 | mm(pop1, pop2)) +
(1 | mm(individual1, individual2)),
data = sympatric_pairs_simple,
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/fits/acoustic_distance_sympatry_geographic_distance_simple_fit"
)2.2.6.1.1.2 Fit summary
Code
extended_summary(
sympatry_geo_model_simple,
highlight = TRUE,
trace.palette = viridis::mako,
remove.intercepts = 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, 2) Intercept-normal(0, 5) sd-exponential(1) sigma-exponential(1) | acoustic_distance_sc ~ sympatry + scale(geo_distance) + (1 | species_pair) + (1 | mm(pop1, pop2)) + (1 | mm(individual1, individual2)) | 10000 | 4 | 1 | 5000 | 924 (0.046%) | 0 | 20124.31 | 10999.95 | 489747410 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sympatry1 | 0.462 | 0.396 | 0.527 | 1 | 21530.42 | 11474.75 |
| b_scalegeo_distance | 0.129 | 0.095 | 0.162 | 1 | 20124.31 | 10999.95 |
Summary
Both predictors had credible, positive effects on overall acoustic distance: sympatry (β = 0.462, 95% CI [0.396, 0.527]) and geographic distance (β = 0.129, 95% CI [0.095, 0.162]). Sympatric heterospecific pairs are substantially more acoustically divergent than allopatric pairs, and divergence also increases with geographic distance independent of sympatry — sympatry’s effect is roughly 3.5x the size of geographic distance’s. Trace plots showed good mixing among chains with no obvious trends, and Rhat = 1 with high effective sample sizes for both parameters, indicating satisfactory MCMC convergence. This is a clear multivariate signature of character displacement in simple songs, layered on top of an isolation-by-distance pattern.
2.2.6.1.2 Feature-based acoustic distance
The following plot summarizes the relationship between sympatry and acoustic divergence for each of the acoustic features used to calculate acoustic distance. The model was fit separately for each feature, and the estimated effect of sympatry on acoustic divergence is shown with 95% credible intervals. Rows correspond to individual acoustic traits and columns to the predictors included in the Bayesian mixed-effects models. Tile color represents the posterior mean regression coefficient (green = positive effect, purple = negative effect, white = no effect), while the numerical value within each tile indicates the estimated effect size. Models accounted for non-independence among pairwise comparisons by including multi-membership random effects for populations and individuals, as well as a random intercept for species pair.
Code
# identify all response variables
distance_vars <- grep(
"_distance$",
names(sympatric_pairs_simple),
value = TRUE
)
#remove geographic distance itself
distance_vars <- setdiff(distance_vars, c("geo_distance", "acoustic_distance", "log_geo_distance"))Code
# scale every per-feature response the same way (mirrors the complex-song
# loop), so all traits and both song types are on a comparable footing
distance_vars_sc <- paste0(distance_vars, "_sc")
for (v in distance_vars) {
sympatric_pairs_simple[[paste0(v, "_sc")]] <- as.numeric(scale(sympatric_pairs_simple[[v]]))
}
# tighter, genuinely regularizing priors now that the response is
# standardized (mean 0, SD 1) - a coefficient near |1| would already be
# an implausibly large effect
priors <- c(
prior(normal(0, 1), class = "Intercept"),
prior(normal(0, 0.5), class = "b"),
prior(exponential(2), class = "sd"),
prior(exponential(1), class = "sigma")
)
for (resp in distance_vars_sc) {
cat("\n=============================\n")
cat("Fitting:", resp, "\n")
cat("=============================\n")
# ----------------------------
# Sympatry + geographic distance
# ----------------------------
form_sympatry_geo <- bf(
as.formula(
paste0(
resp,
" ~ sympatry + geo_distance_sc + ",
"(1 | species_pair) + ",
"(1 | mm(pop1, pop2)) + ",
"(1 | mm(individual1, individual2))"
)
)
)
mod <- brm(
formula = form_sympatry_geo,
data = sympatric_pairs_simple,
family = gaussian(),
prior = priors,
cores = 4,
chains = 4,
iter = 10000,
backend = "cmdstanr",
threads = threading(8),
control = list(
adapt_delta = 0.95,
max_treedepth = 15
), file_refit = "never",
file = paste0(
"./data/processed/fits/",
resp,
"_sympatry_geographic_distance_simple_fit"
)
)
}Code
# Find and load all saved brms models
model_files <- list.files(
"./data/processed/fits",
pattern = "\\.rds$",
full.names = TRUE
)
model_files <- grep("_complex", model_files, value = TRUE, invert = TRUE)
model_files <- grep(paste(distance_vars, collapse = "|"), model_files, value = TRUE)
# not "_distance_sympatry_..." - a response scaled for its own fit (e.g.
# mst_distance_sc) has a suffix between "distance" and "sympatry", which
# would break that stricter match
model_files <- grep("_sympatry_geographic_distance", model_files, value = TRUE)
plot_brms_heatmap(model_files = model_files)
The expandable section below provides the complete Bayesian model summaries underlying those estimates, including posterior parameter estimates, credible intervals, and convergence diagnostics.
NoteModel results
Code
for (i in model_files) {
model_name <- tools::file_path_sans_ext(basename(i))
model_name <- gsub("_sympatry_geographic_distance.*", "", model_name)
cat("## ", model_name, "\n\n")
extended_summary(
read.file = i,
highlight = TRUE,
trace.palette = viridis::mako,
remove.intercepts = TRUE,
print.name = FALSE
)
cat("\n\n")
}Summary
- Overall acoustic distance (multivariate): both credible & positive — sympatry (β = 0.462, CI [0.396, 0.527]) and geographic distance (β = 0.129, CI [0.095, 0.162]); sympatry’s effect is roughly 3.5× the size of geographic distance’s.
- Frequency range: strongest sympatry divergence (β = 0.38, credible); also diverges with geo distance (β = -0.09, credible).
- Element-level acoustic distance (elm_pca): strong sympatry divergence (β = 0.39, credible); geo distance not credible.
- Number of elements: credible divergence with both sympatry (β = 0.13) and geo distance (β = 0.09).
- Song duration: credible divergence with both sympatry (β = 0.11) and geo distance (β = 0.08).
- Element diversity (mst): credible convergence with both sympatry (β = -0.10) and geo distance (β = -0.36, strongest geo effect overall).
- Peak frequency, gap duration, song rate: no credible effect of either predictor.
2.2.7 Combined results heatmap
The heatmap below combines every simple-song model fitted so far - the song-level PCA-distance model and every per-feature model, including elm_pca_distance - into a single view. Since models may still be running in the background, this only plots whichever fits have already been saved; it updates automatically as more finish, and produces no error if some (or all) are still missing.
Code
# combine every simple-song fit (main song-level model + all per-feature
# models) into a single heatmap; safe to run at any point while models
# are still fitting, since plot_brms_heatmap() skips missing/unreadable
# files rather than erroring
model_files_simple <- list.files(
"./data/processed/fits",
pattern = "\\.rds$",
full.names = TRUE
)
model_files_simple <- grep("simple", model_files_simple, value = TRUE)
plot_brms_heatmap(model_files = model_files_simple)
3 Complex songs
3.1 Prepare element level data
Code
complex_elm <- read.csv("./data/raw/Sporophila song data_elmt-level plus PCA_complex.csv")
# remove all columns that start with "PC"
complex_elm <- complex_elm[, !grepl("^PC", names(complex_elm))]
complex_elm$Population <- complex_elm$Lat <- complex_elm$Lon <- NULL
individual_coord <- read.csv("./data/raw/Coordenadas_individuos_complex_songs.csv")
# assign coordinates to each individual
complex_elm_lat_long <- complex_elm |>
left_join(individual_coord, by = "Individual")3.1.1 Select variables and filter populations
Code
vars <- c(
"Q1.Time..s.",
"Q3.Time..s.",
"Time.5...s.",
"Time.95...s.",
"Delta.Time..s.",
"Dur.90...s.",
"IQR.Dur..s.",
"Peak.Time..s.",
"Center.Time..s.",
"Q1.Freq..Hz.",
"Q3.Freq..Hz.",
"Center.Freq..Hz.",
"Freq.5...Hz.",
"Freq.95...Hz.",
"Delta.Freq..Hz.",
"IQR.BW..Hz.",
"BW.90...Hz.",
"Max.Freq..Hz.",
"Peak.Freq..Hz.",
"meanfreq",
"sd",
"freq.median",
"freq.Q25",
"freq.Q75",
"freq.IQR",
"time.median",
"time.Q25",
"time.Q75",
"time.IQR",
"skew",
"kurt",
"sp.ent",
"time.ent",
"entropy",
"sfm",
"meandom",
"mindom",
"maxdom",
"dfrange",
"modindx",
"startdom",
"enddom",
"dfslope",
"meanpeakf"
)
# select variables that are not highly correlated
complex_elm_dat <- complex_elm_lat_long[, c("Individual", "Lat", "Lon", "species", "Population", "location", "song", vars)]
# remove individuals from underrepresented populations
complex_elm_dat <- filter(complex_elm_dat, Population != "Pal_ER")
complex_elm_dat <- filter(complex_elm_dat, Population != "NA")
complex_elm_dat$Individual <- factor(complex_elm_dat$Individual)
complex_elm_dat$species <- factor(complex_elm_dat$species)
complex_elm_dat$Population <- factor(complex_elm_dat$Population)
complex_elm_dat$location <- factor(complex_elm_dat$location)
complex_elm_dat$song <- factor(complex_elm_dat$song)
complex_elm_dat <- complex_elm_dat[complete.cases(complex_elm_dat), ]3.1.2 Principal Component Analysis
Code
# run PCA on acoustic variables
pca <- prcomp(
complex_elm_dat[, vars],
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)We used the first 5 principal components (PCs), selected using the broken-stick criterion, for subsequent analyses, which together explained 81.5% of the variance in the data.
Code
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
elm_pca_scores <- cbind(
pca$x,
complex_elm_dat[, c(
"Population",
"species",
"location",
"Individual",
"song",
"Lat",
"Lon")]
)
# select the PCs retained by the broken-stick criterion
pcs_elm <- grep("^PC", names(elm_pca_scores), value = TRUE)[seq_len(n_pcs_broken_stick(pca))]
elm_pca_scores <- elm_pca_scores |>
mutate(
across(all_of(pcs_elm), ~ 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
)
)
df_clean_complex_elm <- elm_pca_scores |>
mutate(
across(all_of(pcs_elm), as.numeric),
Lat = as.numeric(Lat),
Lon = as.numeric(Lon),
species = as.character(species),
Individual = as.character(Individual),
song = as.character(song)
) |>
filter(complete.cases(across(all_of(c(pcs_elm, "Lat", "Lon", "species", "Individual", "song")))))3.2 Song level
3.2.1 Read and prepare data
Code
complex_songs <- read.csv(
"./data/raw/Sporophila song data_song-level plus MCP MST and PCA_complex songs.csv",
stringsAsFactors = FALSE
)
individual_coord <- read.csv(
"./data/raw/Coordenadas_individuos_complex_songs.csv",
stringsAsFactors = FALSE
)
# Asignar coordenadas a cada observación según Individual
complex_songs_lat_long <- complex_songs %>%
left_join(individual_coord, by = "Individual")3.2.2 Select variables and filter populations
The song level features used were: peak frequency, number of elements, song duration, song rate, gap duration, frequency range and element diversity (mst)
Element duration was excluded as it is an element level features
Code
# select variables that are not highly correlated
complex_song_dat <- complex_songs_lat_long[, c("Individual", "Lat", "Lon", "species", "Population", "location", "song",
"meanpeakf", "num.elms",
"song.duration", "song.rate", "gap.duration",
"freq.range.Min5toMax95", "mst")]
# remove individuals from underrepresented populations
complex_song_dat <- filter(complex_song_dat, Population != "Pal_ER")
complex_song_dat <- filter(complex_song_dat, Population != "NA")
complex_song_dat$Individual <- factor(complex_song_dat$Individual)
complex_song_dat$species <- factor(complex_song_dat$species)
complex_song_dat$Population <- factor(complex_song_dat$Population)
complex_song_dat$location <- factor(complex_song_dat$location)3.2.3 Principal Component Analysis
Code
# run PCA on acoustic variables
pca <- prcomp(
complex_song_dat[, c(
"meanpeakf",
"num.elms",
"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)We used the first 2 principal components (PCs), selected using the broken-stick criterion, for subsequent analyses, which together explained 73.7% of the variance in the data.
Code
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,
complex_song_dat[, c(
"Population",
"species",
"location",
"Individual",
"song",
"Lat",
"Lon",
"meanpeakf",
"num.elms",
"song.duration",
"song.rate",
"gap.duration",
"freq.range.Min5toMax95",
"mst"
)]
)
# select the PCs retained by the broken-stick criterion
pcs <- grep("^PC", names(pca_scores), value = TRUE)[seq_len(n_pcs_broken_stick(pca))]
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
)
)3.2.4 Acoustic and geographic distances
3.2.4.1 Build pairwise distance dataset
Code
df_clean_complex <- 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),
population = as.character(Population),
song = as.character(song)
) |>
filter(complete.cases(across(all_of(c(pcs, "Lat", "Lon", "species", "Individual", "song")))))
df_clean_complex$population <- sapply(strsplit(df_clean_complex$population, "_"), "[[", 2)
# table(df_clean_complex$Population, df_clean_complex$species)
# song identity uses the native "song" column carried through from the
# raw data (unique per recording: sound file + species + vocalization
# type + song number), rather than being re-derived - unlike the earlier
# row-order-based reconstruction, this can't silently collapse multiple
# distinct songs from the same individual into one song ID3.2.4.2 Convert to long table
Code
# acoustic distance between songs (Euclidean multivariate, unscaled PCs)
agg_df_clean_complex_elm <- aggregate(
. ~ song, df_clean_complex_elm[, c(pcs_elm, "song")],
FUN = mean
)
dist_acoustic_mat <- as.matrix(dist(
df_clean_complex[, pcs],
method = "euclidean"
))
# and for each feature separately (unscaled)
dist_meanpeakf_mat <- as.matrix(dist(
df_clean_complex[, "meanpeakf"],
method = "euclidean"
))
dist_numelms_mat <- as.matrix(dist(
df_clean_complex[, "num.elms"],
method = "euclidean"
))
dist_songduration_mat <- as.matrix(dist(
df_clean_complex[, "song.duration"],
method = "euclidean"
))
dist_songrate_mat <- as.matrix(dist(
df_clean_complex[, "song.rate"],
method = "euclidean"
))
dist_gapduration_mat <- as.matrix(dist(
df_clean_complex[, "gap.duration"],
method = "euclidean"
))
dist_freqrange_mat <- as.matrix(dist(
df_clean_complex[, "freq.range.Min5toMax95"],
method = "euclidean"
))
dist_mst_mat <- as.matrix(dist(
df_clean_complex[, "mst"],
method = "euclidean"
))
dist_elm_mat <- as.matrix(dist(
agg_df_clean_complex_elm[, pcs_elm],
method = "euclidean"
))
# geographic distance (Haversine, km) between songs
coords <- as.matrix(df_clean_complex[, 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_complex$song
rownames(dist_geo_mat) <- colnames(dist_geo_mat) <- df_clean_complex$song
rownames(dist_elm_mat) <- colnames(dist_elm_mat) <- agg_df_clean_complex_elm$song
# 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],
meanpeakf_distance = dist_meanpeakf_mat[idx],
numelms_distance = dist_numelms_mat[idx],
songduration_distance = dist_songduration_mat[idx],
songrate_distance = dist_songrate_mat[idx],
gapduration_distance = dist_gapduration_mat[idx],
freqrange_distance = dist_freqrange_mat[idx],
mst_distance = dist_mst_mat[idx],
geo_distance = dist_geo_mat[idx]
)
dist_acoustic_long$elm_pca_distance <- sapply(seq_len(nrow(dist_acoustic_long)), function(i) {
dist_elm_mat[rownames(dist_elm_mat) == dist_acoustic_long$id1[i], colnames(dist_elm_mat) == dist_acoustic_long$id2[i]]
})
# look up species, population, and individual for each song by matching
# against the native song identifier (mirrors the simple-song approach)
dist_acoustic_long$species1 <- sapply(dist_acoustic_long$id1, function(x)
df_clean_complex$species[df_clean_complex$song == x])
dist_acoustic_long$population1 <- sapply(dist_acoustic_long$id1, function(x)
df_clean_complex$population[df_clean_complex$song == x])
dist_acoustic_long$individual1 <- sapply(dist_acoustic_long$id1, function(x)
df_clean_complex$Individual[df_clean_complex$song == x])
dist_acoustic_long$species2 <- sapply(dist_acoustic_long$id2, function(x)
df_clean_complex$species[df_clean_complex$song == x])
dist_acoustic_long$population2 <- sapply(dist_acoustic_long$id2, function(x)
df_clean_complex$population[df_clean_complex$song == x])
dist_acoustic_long$individual2 <- sapply(dist_acoustic_long$id2, function(x)
df_clean_complex$Individual[df_clean_complex$song == x])
# 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.factor(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 = "_")
)
# offset, centering and scaling are estimated once from the pooled data
# across both song types (see "Geographic distance transform" section)
# so this transform is identical between the simple- and complex-song
# models
dist_acoustic_long$log_geo_distance <- log(dist_acoustic_long$geo_distance + geo_const)
dist_acoustic_long$geo_distance_sc <- transform_geo_distance(dist_acoustic_long$geo_distance)
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)3.2.5 Statistical analysis
To evaluate whether acoustic divergence between heterospecific songs differed between sympatric and allopatric population comparisons, we fitted a Bayesian mixed-effects model while accounting for the non-independence inherent to pairwise distance data.
The model was specified as:
\[ \begin{split} \text{acoustic distance} &\sim \text{sympatry} + \text{geographic distance} \\ &\quad + (1 \mid \text{species pair}) \\ &\quad + (1 \mid \text{mm(population}_1,\text{population}_2)) \\ &\quad + (1 \mid \text{mm(individual}_1,\text{individual}_2)) \end{split} \]
where:
(_{ij}) is the Euclidean distance between songs (i) and (j) in multivariate acoustic space.
(_{ij}) is a binary predictor indicating whether the populations from which songs (i) and (j) were recorded occur in sympatry (1) or allopatry (0).
(_{ij}) is the geographic distance between the populations from which songs (i) and (j) were recorded, log-transformed and scaled (see Geographic distance transform).
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:
- The model was fitted in a Bayesian framework using the brms package, with sympatry and geographic distance specified together as predictors a priori; no alternative (sympatry-only or geographic-distance-only) models were fitted or compared.
- 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, regularizing priors were specified for all parameters.
- The model was fitted using Hamiltonian Monte Carlo as implemented in Stan through the cmdstanr backend.
The coefficient for sympatry in this model represents differences in acoustic divergence between sympatric and allopatric population comparisons, after accounting for geographic distance and the hierarchical structure of the data. The same model structure was fitted separately for the overall PCA-based acoustic distance and for each individual acoustic feature, so that trait-specific effects of sympatry and geographic distance could be examined alongside the overall pattern.
3.2.5.1 Model fitting
3.2.5.1.1 PCA-based acoustic distance
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")
keep_pairs <- rownames(tab)[
tab[, "allopatric"] > 0 &
tab[, "sympatric"] > 0
]
# Extract all species-population combinations
sp_pop <- unique(
rbind(
data.frame(
species = dist_acoustic_long$species1,
population = dist_acoustic_long$population1
),
data.frame(
species = dist_acoustic_long$species2,
population = dist_acoustic_long$population2
)
)
)
# Presence/absence table
tab <- with(
sp_pop,
table(species, population)
)
# Convert counts to X / blank
tab[] <- ifelse(tab > 0, "\u2713", "")
presence_table <- as.data.frame.matrix(tab)
presence_table$species <- rownames(presence_table)
presence_table <- presence_table[
, c("species", setdiff(names(presence_table), "species"))
]
# print(presence_table)Code
# Species x population occurrence matrix
occ <- as.matrix(tab)
species <- rownames(occ)
# All pairwise species combinations
pairs <- combn(species, 2, simplify = FALSE)
results <- data.frame(
Pair = character(),
Sympatric = character(),
Allopatric = character(),
stringsAsFactors = FALSE
)
for(p in pairs) {
sp1 <- p[1]
sp2 <- p[2]
pops1 <- colnames(occ)[occ[sp1, ] != ""]
pops2 <- colnames(occ)[occ[sp2, ] != ""]
sympatric <- intersect(pops1, pops2)
if(length(sympatric) == 0) {
sympatric_txt <- "none"
} else {
sympatric_txt <- paste(sympatric, collapse = ", ")
}
allopatric <- setdiff(union(pops1, pops2), sympatric)
if(length(allopatric) == 0) {
allopatric_txt <- "none"
} else {
allopatric_txt <- paste(allopatric, collapse = ", ")
}
results <- rbind(
results,
data.frame(
Pair = paste(sp1, sp2, sep = "–"),
Sympatric = sympatric_txt,
Allopatric = allopatric_txt,
stringsAsFactors = FALSE
)
)
}
# print(results)Only three species pairs were kept: Hypoxantha_Iberaensis, Hypoxantha_Palustris, Hypoxantha_Ruficollis
3.2.5.1.1.1 Species pairs with both sympatric and allopatric populations
Code
sympatric_pairs_complex <- dist_acoustic_long[
dist_acoustic_long$species_pair %in%
keep_pairs,
]
# create species-specific population IDs
sympatric_pairs_complex$pop1 <- interaction(
sympatric_pairs_complex$species1,
sympatric_pairs_complex$population1,
drop = TRUE
)
sympatric_pairs_complex$pop2 <- interaction(
sympatric_pairs_complex$species2,
sympatric_pairs_complex$population2,
drop = TRUE
)
# make species_pair a factor (fixed effect)
sympatric_pairs_complex$species_pair <- factor(sympatric_pairs_complex$species_pair)
# scale acoustic_distance
sympatric_pairs_complex$acoustic_distance_sc <- scale(sympatric_pairs_complex$acoustic_distance)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")
)
#fit model
sympatry_geo_model_complex <- brm(
acoustic_distance_sc ~ sympatry +
geo_distance_sc +
(1 | species_pair) +
(1 | mm(pop1, pop2)) +
(1 | mm(individual1, individual2)),
data = sympatric_pairs_complex,
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/fits/acoustic_distance_sympatry_geographic_distance_complex_fit"
)
# beepr::beep(3)3.2.5.1.1.2 Fit summary
Code
extended_summary(
sympatry_geo_model_complex,
highlight = TRUE,
trace.palette = viridis::mako,
remove.intercepts = 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, 2) Intercept-normal(0, 5) sd-exponential(1) sigma-exponential(1) | acoustic_distance_sc ~ sympatry + scale(geo_distance) + (1 | species_pair) + (1 | mm(pop1, pop2)) + (1 | mm(individual1, individual2)) | 10000 | 4 | 1 | 5000 | 346 (0.017%) | 0 | 21025.34 | 13991.35 | 820698594 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sympatry1 | -0.062 | -0.144 | 0.019 | 1 | 21025.34 | 14481.96 |
| b_scalegeo_distance | -0.073 | -0.125 | -0.022 | 1 | 21063.64 | 13991.35 |
Summary
- Overall acoustic distance (multivariate): geo distance credible & negative (β = -0.073, CI [-0.125, -0.022]) → farther populations slightly more similar; sympatry not credible (β = -0.062, CI crosses zero).
- Element-level acoustic distance (elm_pca): by far the largest effect in the study — sympatry β = 1.08 (credible, strong divergence), geo distance β = -0.23 (credible).
- Frequency range: credible convergence with both sympatry (β = -0.17) and geo distance (β = -0.14).
- Peak frequency: credible convergence with both sympatry (β = -0.15) and geo distance (β = -0.10).
- Gap duration: credible sympatry divergence (β = 0.11); geo distance not credible.
- Song rate: credible sympatry convergence (β = -0.10); geo distance not credible.
- Song duration, number of elements, element diversity (mst): no credible effect of either predictor.
3.2.5.1.2 Feature-based acoustic distance
The following plot summarizes the relationship between sympatry and acoustic divergence for each of the acoustic features used to calculate acoustic distance. The model was fit separately for each feature, and the estimated effect of sympatry on acoustic divergence is shown with 95% credible intervals. Rows correspond to individual acoustic traits and columns to the predictors included in the Bayesian mixed-effects models. Tile color represents the posterior mean regression coefficient (green = positive effect, purple = negative effect, white = no effect), while the numerical value within each tile indicates the estimated effect size. Models accounted for non-independence among pairwise comparisons by including multi-membership random effects for populations and individuals, as well as a random intercept for species pair.
Code
# identify all response variables
distance_vars <- grep(
"_distance$",
names(sympatric_pairs_complex),
value = TRUE
)
#remove geographic distance itself
distance_vars <- setdiff(distance_vars, c("geo_distance", "acoustic_distance", "log_geo_distance"))
# scale mst
# distance_vars[distance_vars == "mst_distance"] <- "mst_distance_sc"
# sympatric_pairs_complex$mst_distance_sc <- scale(sympatric_pairs_complex$mst_distance)Code
distance_vars_sc <- paste0(distance_vars, "_sc")
for (v in distance_vars) {
sympatric_pairs_complex[[paste0(v, "_sc")]] <- as.numeric(scale(sympatric_pairs_complex[[v]]))
}
# tighter, genuinely regularizing priors now that the response is
# standardized (mean 0, SD 1) - a coefficient near |1| would already be
# an implausibly large effect
priors <- c(
prior(normal(0, 1), class = "Intercept"),
prior(normal(0, 0.5), class = "b"),
prior(exponential(2), class = "sd"),
prior(exponential(1), class = "sigma")
)
for (resp in distance_vars_sc) {
cat("\n=============================\n")
cat("Fitting:", resp, "\n")
cat("=============================\n")
form_sympatry_geo <- bf(
as.formula(
paste0(
resp,
" ~ sympatry + geo_distance_sc + ",
"(1 | species_pair) + ",
"(1 | mm(pop1, pop2)) + ",
"(1 | mm(individual1, individual2))"
)
)
)
mod <- brm(
formula = form_sympatry_geo,
data = sympatric_pairs_complex,
family = gaussian(),
prior = priors,
cores = 4, chains = 4, iter = 10000,
backend = "cmdstanr", threads = threading(8),
control = list(adapt_delta = 0.95, max_treedepth = 15),
file_refit = "always",
file = paste0("./data/processed/fits/", resp, "_sympatry_geographic_distance_complex_fit")
)
}Code
# Find and load all saved brms models
model_files <- list.files(
"./data/processed/fits",
pattern = "\\.rds$",
full.names = TRUE
)
model_files <- grep("_complex", model_files, value = TRUE)
model_files <- grep(paste(distance_vars, collapse = "|"), model_files, value = TRUE)
# not "_distance_sympatry_..." - every complex per-feature response is
# refit under a "_sc" suffix (see distance_vars_sc above), which sits
# between "distance" and "sympatry" and would break that stricter match
model_files <- grep("_sympatry_geographic_distance", model_files, value = TRUE)
plot_brms_heatmap(model_files = model_files)
The expandable section below provides the complete Bayesian model summaries underlying those estimates, including posterior parameter estimates, credible intervals, and convergence diagnostics.
NoteModel results
Code
for (i in model_files) {
model_name <- tools::file_path_sans_ext(basename(i))
model_name <- gsub("_sympatry_geographic_distance.*", "", model_name)
cat("## ", model_name, "\n\n")
extended_summary(
read.file = i,
highlight = TRUE,
trace.palette = viridis::mako,
remove.intercepts = TRUE,
print.name = FALSE
)
cat("\n\n")
}Summary
In contrast to the simple-song analyses, sympatry and geographic distance were associated primarily with reduced (convergent) acoustic differences among most complex-song traits — with one striking exception.
Element-level acoustic distance (elm_pca) showed by far the largest effect in the entire study: sympatry β = 1.08 (credible, strong divergence), geographic distance β = -0.23 (credible). Sympatric pairs are much more divergent in fine element-level structure, even as the broader song-level features converge or show no effect.
Frequency range and peak frequency both converged credibly with both predictors (frequency range: sympatry β = -0.17, geo β = -0.14; peak frequency: sympatry β = -0.15, geo β = -0.10).
Gap duration diverged credibly with sympatry (β = 0.11) but showed no credible geographic-distance effect (β = 0.05).
Song rate converged credibly with sympatry (β = -0.10) but showed no credible geographic-distance effect (β = -0.02).
Song duration, number of elements, and element diversity (mst) showed no credible effect of either predictor.
Overall, complex songs show a mixed pattern: general convergence or no effect across most individual song-level traits, but very strong, narrowly focused divergence in fine element-level structure — the opposite of the broad divergence seen across several simple-song traits.
Taken together, these findings suggest that character displacement in complex songs is limited and trait-specific, with clear evidence for divergence confined to element-level structure and, to a lesser extent, gap duration, whereas the broader suite of song-level acoustic characteristics appears conserved or convergent among sympatric species.
3.2.6 Combined results heatmap
The heatmap below combines every complex-song model fitted so far - the song-level PCA-distance model and every per-feature model, including elm_pca_distance - into a single view. Since models may still be running in the background, this only plots whichever fits have already been saved; it updates automatically as more finish, and produces no error if some (or all) are still missing.
Code
# combine every complex-song fit (main song-level model + all
# per-feature models) into a single heatmap; safe to run at any point
# while models are still fitting, since plot_brms_heatmap() skips
# missing/unreadable files rather than erroring
model_files_complex <- list.files(
"./data/processed/fits",
pattern = "\\.rds$",
full.names = TRUE
)
model_files_complex <- grep("complex", model_files_complex, value = TRUE)
plot_brms_heatmap(model_files = model_files_complex)
4 Overall summary
- Character displacement (sympatry-driven divergence) is clear and broad in simple songs, especially frequency range and element structure, and is layered on top of an isolation-by-distance effect at the overall multivariate level.
- In complex songs, the overall song-level signal leans toward convergence (mainly geography-driven), but element-level structure diverges very strongly with sympatry — the single biggest effect found.
NoteSession information
─ Session info ───────────────────────────────────────────────────────────────
setting value
version R version 4.5.2 (2025-10-31)
os Ubuntu 22.04.4 LTS
system x86_64, linux-gnu
ui X11
language (EN)
collate en_US.UTF-8
ctype en_US.UTF-8
tz America/Costa_Rica
date 2026-09-07
pandoc 3.6.3 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/x86_64/ (via rmarkdown)
quarto 1.8.25 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/quarto
─ Packages ───────────────────────────────────────────────────────────────────
package * version date (UTC) lib source
abind 1.4-8 2024-09-12 [1] CRAN (R 4.5.2)
ape 5.8-1 2024-12-16 [1] CRAN (R 4.5.2)
arrayhelpers 1.1-0 2020-02-04 [1] CRAN (R 4.5.2)
backports 1.5.1 2026-04-03 [1] CRAN (R 4.5.2)
bayesplot 1.15.0 2025-12-12 [1] CRAN (R 4.5.2)
bridgesampling 1.2-1 2025-11-19 [1] CRAN (R 4.5.2)
brms * 2.23.0 2025-09-09 [1] CRAN (R 4.5.2)
brmsish * 1.0.0 2026-03-03 [1] Github (maRce10/brmsish@81ab826)
Brobdingnag 1.2-9 2022-10-19 [1] CRAN (R 4.5.2)
cachem 1.1.0 2024-05-16 [1] CRAN (R 4.5.2)
car 3.1-5 2026-02-03 [1] CRAN (R 4.5.2)
carData 3.0-6 2026-01-30 [1] CRAN (R 4.5.2)
checkmate 2.3.4 2026-02-03 [1] CRAN (R 4.5.2)
cli 3.6.6 2026-04-09 [1] CRAN (R 4.5.2)
cluster 2.1.8.2 2026-02-05 [1] CRAN (R 4.5.2)
cmdstanr * 0.9.0 2025-03-30 [1] https://stan-dev.r-universe.dev (R 4.5.2)
coda 0.19-4.1 2024-01-31 [1] CRAN (R 4.5.2)
codetools 0.2-20 2024-03-31 [1] CRAN (R 4.5.2)
commonmark 2.0.0 2025-07-07 [1] CRAN (R 4.5.2)
cowplot 1.2.0 2025-07-07 [1] CRAN (R 4.5.2)
crayon 1.5.3 2024-06-20 [1] CRAN (R 4.5.2)
curl 7.1.0 2026-04-22 [1] CRAN (R 4.5.2)
devtools 2.4.6 2025-10-03 [1] CRAN (R 4.5.2)
digest 0.6.39 2025-11-19 [1] CRAN (R 4.5.2)
distributional 0.6.0 2026-01-14 [1] CRAN (R 4.5.2)
dplyr * 1.2.1 2026-04-03 [1] CRAN (R 4.5.2)
ecodist * 2.1.3 2023-10-30 [1] CRAN (R 4.5.2)
ellipsis 0.3.2 2021-04-29 [3] CRAN (R 4.1.1)
emmeans 2.0.4 2026-07-15 [1] CRAN (R 4.5.2)
estimability 2.0.0 2026-06-26 [1] CRAN (R 4.5.2)
evaluate 1.0.5 2025-08-27 [1] CRAN (R 4.5.2)
farver 2.1.2 2024-05-13 [1] CRAN (R 4.5.2)
fastmap 1.2.0 2024-05-15 [1] CRAN (R 4.5.2)
forcats * 1.0.1 2025-09-25 [1] CRAN (R 4.5.2)
Formula 1.2-5 2023-02-24 [1] CRAN (R 4.5.2)
fs 2.1.0 2026-04-18 [1] CRAN (R 4.5.2)
generics 0.1.4 2025-05-09 [1] CRAN (R 4.5.2)
geosphere * 1.6-5 2026-03-02 [1] CRAN (R 4.5.2)
ggdist 3.3.3 2025-04-23 [1] CRAN (R 4.5.2)
ggplot2 * 4.0.3 2026-04-22 [1] CRAN (R 4.5.2)
ggtext * 0.1.2 2022-09-16 [1] CRAN (R 4.5.2)
glue 1.8.1 2026-04-17 [1] CRAN (R 4.5.2)
gridExtra 2.3.1 2026-06-25 [1] CRAN (R 4.5.2)
gridtext 0.1.6 2026-02-19 [1] CRAN (R 4.5.2)
gtable 0.3.6 2024-10-25 [1] CRAN (R 4.5.2)
hms 1.1.4 2025-10-17 [1] CRAN (R 4.5.2)
htmltools 0.5.9 2025-12-04 [1] CRAN (R 4.5.2)
htmlwidgets 1.6.4 2023-12-06 [1] CRAN (R 4.5.2)
igraph 2.3.3 2026-06-26 [1] CRAN (R 4.5.2)
inline 0.3.21 2025-01-09 [1] CRAN (R 4.5.2)
jsonlite 2.0.0 2025-03-27 [1] CRAN (R 4.5.2)
kableExtra * 1.4.0 2024-01-24 [1] CRAN (R 4.5.2)
knitr * 1.51 2025-12-20 [1] CRAN (R 4.5.2)
labeling 0.4.3 2023-08-29 [1] CRAN (R 4.5.2)
lattice 0.22-9 2026-02-09 [1] CRAN (R 4.5.2)
lifecycle 1.0.5 2026-01-08 [1] CRAN (R 4.5.2)
litedown 0.9 2025-12-18 [1] CRAN (R 4.5.2)
loo 2.9.0 2025-12-23 [1] CRAN (R 4.5.2)
lubridate * 1.9.5 2026-02-04 [1] CRAN (R 4.5.2)
magrittr 2.0.5 2026-04-04 [1] CRAN (R 4.5.2)
markdown 2.0 2025-03-23 [1] CRAN (R 4.5.2)
MASS 7.3-65 2025-02-28 [1] CRAN (R 4.5.2)
Matrix 1.7-4 2025-08-28 [1] CRAN (R 4.5.2)
matrixStats 1.5.0 2025-01-07 [1] CRAN (R 4.5.2)
memoise 2.0.1 2021-11-26 [3] CRAN (R 4.1.2)
mgcv 1.9-4 2025-11-07 [1] CRAN (R 4.5.2)
multcomp 1.4-29 2025-10-20 [1] CRAN (R 4.5.2)
mvtnorm 1.3-3 2025-01-10 [1] CRAN (R 4.5.2)
nlme 3.1-168 2025-03-31 [1] CRAN (R 4.5.2)
otel 0.2.0 2025-08-29 [1] CRAN (R 4.5.2)
packrat 0.9.3 2025-06-16 [1] CRAN (R 4.5.2)
pbapply 1.7-4 2025-07-20 [1] CRAN (R 4.5.2)
permute * 0.9-10 2026-02-06 [1] CRAN (R 4.5.2)
pillar 1.11.1 2025-09-17 [1] CRAN (R 4.5.2)
pkgbuild 1.4.8 2025-05-26 [1] CRAN (R 4.5.2)
pkgconfig 2.0.3 2019-09-22 [3] CRAN (R 4.0.1)
pkgload 1.5.3 2026-06-15 [1] CRAN (R 4.5.2)
plyr 1.8.9 2023-10-02 [1] CRAN (R 4.5.2)
posterior 1.6.1 2025-02-27 [1] CRAN (R 4.5.2)
processx 3.9.0 2026-04-22 [1] CRAN (R 4.5.2)
ps 1.9.3 2026-04-20 [1] CRAN (R 4.5.2)
purrr * 1.2.2 2026-04-10 [1] CRAN (R 4.5.2)
QuickJSR 1.9.0 2026-01-25 [1] CRAN (R 4.5.2)
R6 2.6.1 2025-02-15 [1] CRAN (R 4.5.2)
RColorBrewer 1.1-3 2022-04-03 [1] CRAN (R 4.5.2)
Rcpp * 1.1.2 2026-07-05 [1] CRAN (R 4.5.2)
RcppParallel 5.1.11-2 2026-03-05 [1] CRAN (R 4.5.2)
readr * 2.2.0 2026-02-19 [1] CRAN (R 4.5.2)
remotes 2.5.0 2024-03-17 [1] CRAN (R 4.5.2)
reshape2 1.4.5 2025-11-12 [1] CRAN (R 4.5.2)
rlang 1.3.0 2026-07-05 [1] CRAN (R 4.5.2)
rmarkdown 2.31 2026-03-26 [1] CRAN (R 4.5.2)
rstan 2.32.7 2025-03-10 [1] CRAN (R 4.5.2)
rstantools 2.6.0 2026-01-10 [1] CRAN (R 4.5.2)
rstudioapi 0.18.0 2026-01-16 [1] CRAN (R 4.5.2)
S7 0.2.2 2026-04-22 [1] CRAN (R 4.5.2)
sandwich 3.1-1 2024-09-15 [1] CRAN (R 4.5.2)
scales 1.4.0 2025-04-24 [1] CRAN (R 4.5.2)
sessioninfo 1.2.3 2025-02-05 [1] CRAN (R 4.5.2)
sketchy 1.0.7 2026-03-03 [1] CRANs (R 4.5.2)
sp 2.2-3 2026-07-19 [1] CRAN (R 4.5.2)
StanHeaders 2.32.10 2024-07-15 [1] CRAN (R 4.5.2)
stringi 1.8.7 2025-03-27 [1] CRAN (R 4.5.2)
stringr * 1.6.0 2025-11-04 [1] CRAN (R 4.5.2)
survival 3.8-6 2026-01-16 [1] CRAN (R 4.5.2)
svglite 2.2.2 2025-10-21 [1] CRAN (R 4.5.2)
svUnit 1.0.8 2025-08-26 [1] CRAN (R 4.5.2)
systemfonts 1.3.1 2025-10-01 [1] CRAN (R 4.5.2)
tensorA 0.36.2.1 2023-12-13 [1] CRAN (R 4.5.2)
textshaping 1.0.4 2025-10-10 [1] CRAN (R 4.5.2)
TH.data 1.1-5 2025-11-17 [1] CRAN (R 4.5.2)
tibble * 3.3.1 2026-01-11 [1] CRAN (R 4.5.2)
tidybayes 3.0.7 2024-09-15 [1] CRAN (R 4.5.2)
tidyr * 1.3.2 2025-12-19 [1] CRAN (R 4.5.2)
tidyselect 1.2.1 2024-03-11 [1] CRAN (R 4.5.2)
tidyverse * 2.0.0 2023-02-22 [1] CRAN (R 4.5.2)
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.2)
usethis 3.2.1 2025-09-06 [1] CRAN (R 4.5.2)
V8 8.2.0 2026-04-21 [1] CRAN (R 4.5.2)
vctrs 0.7.3 2026-04-11 [1] CRAN (R 4.5.2)
vegan * 2.7-5 2026-05-25 [1] CRAN (R 4.5.2)
viridis 0.6.5 2024-01-29 [1] CRAN (R 4.5.2)
viridisLite 0.4.3 2026-02-04 [1] CRAN (R 4.5.2)
withr 3.0.3 2026-06-19 [1] CRAN (R 4.5.2)
xaringanExtra 0.8.0 2024-05-19 [1] CRAN (R 4.5.2)
xfun 0.60 2026-07-09 [1] CRAN (R 4.5.2)
xml2 1.5.2 2026-01-17 [1] CRAN (R 4.5.2)
xtable 1.8-8 2026-02-22 [1] CRAN (R 4.5.2)
yaml 2.3.12 2025-12-10 [1] CRAN (R 4.5.2)
zoo 1.8-15 2025-12-15 [1] CRAN (R 4.5.2)
[1] /home/m/R/x86_64-pc-linux-gnu-library/4.5
[2] /usr/local/lib/R/site-library
[3] /usr/lib/R/site-library
[4] /usr/lib/R/library
* ── Packages attached to the search path.
──────────────────────────────────────────────────────────────────────────────