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("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 <- 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
models <- list.files("./data/processed/", pattern = "fit", full.names = TRUE)
models <- grep("acoustic", models, value = TRUE)
models <- grep("all_data", models, invert = TRUE, value = TRUE)
for (i in models){
extended_summary(
read.file = i,
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 | 752 (0.038%) | 0 | 24515.86 | 13314.61 | 1871852522 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sympatry1 | 0.155 | 0.096 | 0.215 | 1 | 24623.43 | 13339.44 |
| b_scalegeo_distance | -0.001 | -0.031 | 0.029 | 1 | 24515.86 | 13314.61 |
Summary
Pending update: this summary previously reported results from comparing sympatry-only, geographic-distance-only, and combined models via leave-one-out cross-validation (LOO); since only the combined sympatry + geographic distance model is now fitted, that comparison no longer applies. Replace this box with an interpretation of the combined model’s coefficients (magnitude, sign, and credible interval for sympatry and geo_distance_sc) once the current fit finishes running.
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"))
# scale mst
distance_vars[distance_vars == "mst_distance"] <- "mst_distance_sc"
sympatric_pairs_simple$mst_distance_sc <- scale(sympatric_pairs_simple$mst_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")
)
for (resp in distance_vars) {
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)
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("_distance_sympatry_geographic_distance_fit", "", 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, sympatry had stronger effects than geographic distance for most acoustic traits, suggesting that coexistence explains acoustic divergence beyond simple geographic separation.
Element duration showed positive effects of both predictors. Sympatric species differed by an additional 170 ms in element duration. A one standard deviation increase in geographic distance was associated with an additional 80 ms difference in element duration.
Frequency range exhibited the strongest positive sympatry effect. Sympatric species differed by an additional 300 Hz in frequency range, whereas a one standard deviation increase in geographic distance was associated with an increase of 60 Hz.
Song duration also increased with both predictors. Sympatric species differed by an additional 140 ms in song duration, while a one standard deviation increase in geographic distance was associated with an increase of 50 ms.
Song rate showed modest positive effects of both predictors. Sympatric species differed by approximately 0.07 elements s⁻¹, while a one standard deviation increase in geographic distance was associated with an increase of 0.06 elements s⁻¹.
Number of elements was only weakly associated with the predictors. Sympatric species differed by approximately 0.05 elements, whereas geographic distance had essentially no effect.
Peak frequency showed negative effects of both predictors. Sympatric species were approximately 60 Hz more similar in peak frequency than allopatric species. Likewise, a one standard deviation increase in geographic distance was associated with a reduction of approximately 20 Hz in mean peak-frequency differences.
Gap duration also showed weak negative effects. Sympatric species were approximately 30 ms more similar in gap duration, while a one standard deviation increase in geographic distance was associated with a reduction of approximately 10 ms in gap-duration differences.
The heterogeneous pattern among song traits indicates that acoustic divergence is trait-specific. Some song features become more divergent in sympatry (particularly frequency range, element duration, and song duration), whereas others remain relatively conserved or even become slightly more similar (mean peak frequency and gap duration).
Collectively, these results suggest that selection associated with sympatry acts on a subset of song characteristics, rather than producing uniform divergence across the entire acoustic phenotype. This pattern is consistent with the hypothesis that only particular components of song contribute strongly to species recognition and reproductive isolation.
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
3.2.6 Interpretation
Both sympatry and geographic distance had negative effects on overall acoustic distance in the multivariate PCA space, indicating that complex songs tended to become more similar, rather than more divergent, under both conditions.
The effect of sympatry was stronger than that of geographic distance. After accounting for species pair, populations, and individuals, sympatric species exhibited an average reduction of approximately 0.12 standard deviations in multivariate acoustic distance relative to allopatric species.
A one standard deviation increase in geographic distance was associated with an average reduction of approximately 0.10 standard deviations in multivariate acoustic distance.
The 95% credible intervals for both predictors excluded zero, providing strong evidence that the observed effects are consistently negative.
Trace plots showed good mixing among chains and no obvious trends, indicating satisfactory MCMC convergence for both regression coefficients.
These multivariate results agree with the analyses of individual song traits, in which most components of complex songs showed either weak effects or evidence of convergence rather than divergence. Consequently, when considered jointly in multivariate acoustic space, complex songs exhibit an overall pattern of reduced acoustic differentiation among sympatric species.
This contrasts with the results for simple songs, where several acoustic traits exhibited positive effects of sympatry, suggesting that character displacement is expressed primarily in simple songs rather than in complex songs.
Biologically, these findings suggest that the acoustic features composing complex songs may be evolutionarily conserved or constrained, whereas simple songs appear to be the principal target of divergence associated with species coexistence and species recognition.
3.2.6.0.1 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[8]) {
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))"
)
)
)
sympatry_geo_models[[resp]] <- 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)
model_files <- grep("_distance_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("_distance_sympatry_geographic_distance_fit", "", 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, both sympatry and geographic distance were associated primarily with reduced acoustic differences among most complex-song traits, suggesting convergence or conservation rather than divergence.
Song rate exhibited the strongest negative effects of both predictors. Sympatric species differed by approximately 0.35 elements s⁻¹ less than allopatric species, while a one standard deviation increase in geographic distance was associated with a reduction of 0.19 elements s⁻¹ in song-rate differences.
Frequency range also showed consistent negative effects. Sympatric species were approximately 150 Hz more similar in frequency range than allopatric species, and a one standard deviation increase in geographic distance was associated with a reduction of 130 Hz in frequency-range differences.
Gap duration was the only trait exhibiting clear positive effects of both predictors. Sympatric species differed by approximately 160 ms more in gap duration than allopatric species, while a one standard deviation increase in geographic distance was associated with an increase of 70 ms.
Song duration showed only weak negative effects, with sympatric species differing by approximately 40 ms less than allopatric species and a one standard deviation increase in geographic distance associated with a reduction of 30 ms.
Number of elements exhibited negligible effects, with sympatric species differing by only 0.02 fewer elements and geographic distance associated with a reduction of approximately 0.01 elements.
Mean peak frequency was essentially unaffected by either predictor, with estimated effects close to zero (≤ 10 Hz).
Overall, these results suggest that most components of complex songs become slightly more similar rather than more divergent in sympatry, in marked contrast to the stronger divergence observed for several simple-song traits.
The only clear exception was gap duration, which became substantially more divergent in sympatry, indicating that this temporal feature may be particularly responsive to selection associated with species coexistence.
Taken together, these findings suggest that character displacement in complex songs is limited and trait-specific, with evidence for divergence confined largely to gap duration, whereas other acoustic characteristics appear to remain conserved or converge among sympatric species.
3.2.7 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)
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-03
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.
──────────────────────────────────────────────────────────────────────────────