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",
"ggh4x",
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("\\)$", "", 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)
# 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("acoustic_distance$", "acoustic\ndistance", plot_df$response)
plot_df$response <- gsub("_distance$", "", plot_df$response)
plot_df$predictor <- gsub("geo_", "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
)
)
}
# table of fixed-effect estimates (posterior mean, SE and 95% uncertainty
# interval) for a set of saved brms fits, one row per model x predictor;
# skips missing/unreadable files so it can be run while models are still
# fitting
fixef_table <- function(model_files, remove_intercepts = TRUE, digits = 3) {
rows <- lapply(model_files, function(f) {
fit <- tryCatch(readRDS(f), error = function(e) NULL)
if (is.null(fit)) return(NULL)
fe <- as.data.frame(fixef(fit))
fe$predictor <- rownames(fe)
if (remove_intercepts) fe <- fe[fe$predictor != "Intercept", ]
fe$response <- deparse(fit$formula$formula[[2]])
rm(fit)
invisible(gc(verbose = FALSE))
fe
})
out <- do.call(rbind, rows)
if (is.null(out) || nrow(out) == 0) {
message("fixef_table: no readable model results yet - skipping.")
return(invisible(NULL))
}
rownames(out) <- NULL
out$response <- gsub("_distance(_sc)?$", "", out$response)
out$predictor <- gsub("^sympatry1$", "Sympatry", out$predictor)
out$predictor <- gsub("^geo_distance_sc$|^scalegeo_distance$", "Geographic distance", out$predictor)
out$ui_excludes_zero <- ifelse(out$Q2.5 > 0 | out$Q97.5 < 0, "yes", "no")
out <- out[, c("response", "predictor", "Estimate", "Est.Error", "Q2.5", "Q97.5", "ui_excludes_zero")]
out[, 3:6] <- round(out[, 3:6], digits)
names(out) <- c("Response", "Predictor", "Estimate", "SE", "l-95% UI", "u-95% UI", "UI excludes 0")
out[order(out$Response, out$Predictor), ]
}
# ---- contrasts by species pair -------------------------------------------
# sympatry effect for each species pair from a model with species-pair-
# specific sympatry slopes, i.e. (1 + sympatry | species_pair): per-pair
# effect = fixed sympatry effect + the pair's deviation. Returns the per-pair
# effects (plus the pooled fixed effect), the pairwise differences between
# species pairs, and a plot.
species_pair_sympatry_effects <- function(fit, digits = 3) {
# summary per species pair (fixed + random slope)
cf <- coef(fit, summary = TRUE)$species_pair[, , "sympatry1", drop = TRUE]
cf <- as.data.frame(cf)
cf$species_pair <- rownames(cf)
rownames(cf) <- NULL
pooled <- fixef(fit)["sympatry1", ]
cf <- rbind(
cf,
data.frame(
Estimate = pooled[["Estimate"]], Est.Error = pooled[["Est.Error"]],
Q2.5 = pooled[["Q2.5"]], Q97.5 = pooled[["Q97.5"]],
species_pair = "Pooled (fixed effect)"
)
)
cf$ui_excludes_zero <- ifelse(cf$Q2.5 > 0 | cf$Q97.5 < 0, "yes", "no")
effects <- cf[, c("species_pair", "Estimate", "Est.Error", "Q2.5", "Q97.5", "ui_excludes_zero")]
effects[, 2:5] <- round(effects[, 2:5], digits)
names(effects) <- c("Species pair", "Sympatry effect", "SE", "l-95% UI", "u-95% UI", "UI excludes 0")
# pairwise differences between species pairs, from the posterior draws
draws <- coef(fit, summary = FALSE)$species_pair[, , "sympatry1"]
lev <- colnames(draws)
diffs <- list()
if (length(lev) > 1) {
cmb <- combn(lev, 2)
for (k in seq_len(ncol(cmb))) {
d <- draws[, cmb[1, k]] - draws[, cmb[2, k]]
diffs[[k]] <- data.frame(
contrast = paste(cmb[1, k], "-", cmb[2, k]),
Estimate = mean(d),
Q2.5 = unname(quantile(d, 0.025)),
Q97.5 = unname(quantile(d, 0.975))
)
}
}
diffs <- do.call(rbind, diffs)
if (!is.null(diffs)) {
diffs$ui_excludes_zero <- ifelse(diffs$Q2.5 > 0 | diffs$Q97.5 < 0, "yes", "no")
diffs[, 2:4] <- round(diffs[, 2:4], digits)
names(diffs) <- c("Contrast", "Difference", "l-95% UI", "u-95% UI", "UI excludes 0")
}
plot_df <- cf
plot_df$species_pair <- factor(plot_df$species_pair, levels = rev(cf$species_pair))
plot_df$type <- ifelse(plot_df$species_pair == "Pooled (fixed effect)", "Pooled", "Species pair")
p <- ggplot(plot_df, aes(x = Estimate, y = species_pair, color = type)) +
geom_vline(xintercept = 0, linetype = "dashed", color = "grey60") +
geom_errorbar(aes(xmin = Q2.5, xmax = Q97.5), width = 0, orientation = "y", linewidth = 1) +
geom_point(size = 3) +
scale_color_manual(values = c("Species pair" = "#403B78", "Pooled" = "#5FA98A"), guide = "none") +
labs(x = "Sympatry effect (SD of acoustic distance) with 95% UI", y = NULL) +
theme_classic()
list(effects = effects, differences = diffs, plot = p)
}
# ---- contrasts among population pairs -------------------------------------
# expected (standardized) acoustic distance for every observed population-
# pair combination within each species pair, from the fitted model
# including the species-pair and population (multi-membership) effects but
# not the individual effects; then sympatric vs allopatric population-pair
# contrasts within species pairs. geo = "observed" uses each population
# pair's mean geographic distance; geo = "mean" sets it to the pooled mean
# (geo_distance_sc = 0) so that contrasts reflect sympatry + population
# effects only.
population_pair_contrasts <- function(fit, dat, geo = c("observed", "mean"), ndraws = 2000, digits = 3) {
geo <- match.arg(geo)
# put the two populations of each comparison in a fixed (alphabetical)
# order, so that A-B and B-A comparisons are pooled into one population
# pair (the multiple-membership term is symmetric, so order does not
# affect the prediction)
dat <- as.data.frame(dat)
p1 <- as.character(dat$pop1)
p2 <- as.character(dat$pop2)
swap <- p1 > p2
dat$pop_a <- ifelse(swap, p2, p1)
dat$pop_b <- ifelse(swap, p1, p2)
dat$geo_distance_sc <- as.numeric(dat$geo_distance_sc)
nd <- aggregate(
geo_distance_sc ~ species_pair + pop_a + pop_b + sympatry,
data = dat, FUN = mean
)
n_obs <- aggregate(
geo_distance_sc ~ species_pair + pop_a + pop_b + sympatry,
data = dat, FUN = length
)
nd$n_comparisons <- n_obs$geo_distance_sc
names(nd)[names(nd) == "pop_a"] <- "pop1"
names(nd)[names(nd) == "pop_b"] <- "pop2"
if (geo == "mean") nd$geo_distance_sc <- 0
# individual columns are required by brms but excluded from the prediction
nd$individual1 <- dat$individual1[1]
nd$individual2 <- dat$individual2[1]
# keep the species-pair and population terms, drop the individual terms;
# if the model has species-pair-specific sympatry slopes, keep them too
has_sp_slopes <- "sympatry1" %in% dimnames(ranef(fit)$species_pair)[[3]]
re_form <- if (has_sp_slopes) {
~ (1 + sympatry | species_pair) + (1 | mm(pop1, pop2))
} else {
~ (1 | species_pair) + (1 | mm(pop1, pop2))
}
ep <- posterior_epred(
fit,
newdata = nd,
re_formula = re_form,
ndraws = ndraws
)
nd$Estimate <- colMeans(ep)
nd$Q2.5 <- apply(ep, 2, quantile, 0.025)
nd$Q97.5 <- apply(ep, 2, quantile, 0.975)
# pop1/pop2 are "species.population" (interaction(), needed to keep
# population codes that repeat across species unique for the mm() term);
# for display, drop the species prefix - the facet strip (species_pair)
# already identifies which two species are being compared
strip_species <- function(x) sub("^[^.]*\\.", "", as.character(x))
nd$population_pair <- paste(strip_species(nd$pop1), strip_species(nd$pop2), sep = " vs ")
nd$sympatry_label <- ifelse(nd$sympatry == "1", "Sympatric", "Allopatric")
# sympatric vs allopatric population pairs within each species pair
contrasts <- list()
for (sp in unique(as.character(nd$species_pair))) {
idx_s <- which(nd$species_pair == sp & nd$sympatry == "1")
idx_a <- which(nd$species_pair == sp & nd$sympatry == "0")
for (i in idx_s) for (j in idx_a) {
d <- ep[, i] - ep[, j]
contrasts[[length(contrasts) + 1]] <- data.frame(
species_pair = sp,
sympatric_pair = nd$population_pair[i],
allopatric_pair = nd$population_pair[j],
Estimate = mean(d),
Q2.5 = unname(quantile(d, 0.025)),
Q97.5 = unname(quantile(d, 0.975))
)
}
}
contrasts <- do.call(rbind, contrasts)
if (!is.null(contrasts)) {
contrasts$ui_excludes_zero <- ifelse(contrasts$Q2.5 > 0 | contrasts$Q97.5 < 0, "yes", "no")
contrasts[, 4:6] <- round(contrasts[, 4:6], digits)
names(contrasts) <- c("Species pair", "Sympatric population pair", "Allopatric population pair",
"Difference", "l-95% UI", "u-95% UI", "UI excludes 0")
}
expected <- nd[, c("species_pair", "population_pair", "sympatry_label", "n_comparisons",
"geo_distance_sc", "Estimate", "Q2.5", "Q97.5")]
expected[, 5:8] <- round(expected[, 5:8], digits)
names(expected) <- c("Species pair", "Population pair", "Sympatry", "N comparisons",
"Geographic distance (sc)", "Expected distance", "l-95% UI", "u-95% UI")
# one row of nd per (species_pair, population pair) combination, so the
# count per species_pair (in its factor-level order, matching the order
# facet_wrap2() draws panels in) is exactly the number of population-pair
# rows that panel needs
panel_counts <- as.numeric(table(nd$species_pair))
# ggh4x::facet_wrap2() keeps the facet_wrap()-style banner strip across
# the top of each panel (facet_grid()'s "space = free" doesn't), while
# ggh4x::force_panelsizes() still gives each panel a height proportional
# to how many population-pair rows it holds, instead of every panel
# getting the same height regardless of row count
p <- ggplot(nd, aes(x = Estimate, y = reorder(population_pair, Estimate), color = sympatry_label)) +
geom_errorbar(aes(xmin = Q2.5, xmax = Q97.5), width = 0, orientation = "y", linewidth = 1) +
geom_point(size = 3) +
ggh4x::facet_wrap2(~ species_pair, scales = "free_y", ncol = 1) +
ggh4x::force_panelsizes(rows = panel_counts) +
scale_color_manual(values = c("Allopatric" = "#403B78", "Sympatric" = "#5FA98A"), name = NULL) +
labs(x = "Expected acoustic distance (SD) with 95% UI", y = NULL) +
theme_classic() +
theme(legend.position = "top")
list(expected = expected, contrasts = contrasts, plot = p)
}
# ---- population-pair contrasts across every per-feature model -------------
# runs population_pair_contrasts() once per saved per-feature fit (the same
# model_files list assembled for plot_brms_heatmap()/fixef_table() next to
# the feature-based heatmap), tags each resulting contrast with the feature
# it came from, and returns the combined contrasts table plus a named list
# of plots - one per feature, not one plot faceted across features - each
# reusing population_pair_contrasts()'s own plot (which already facets by
# species pair) so every feature gets its own, fully legible figure. Skips
# any file that fails to load or fails to fit (e.g. a model still running,
# or one for which population_pair_contrasts() errors) rather than
# stopping the whole loop.
population_pair_contrasts_by_feature <- function(model_files, dat, geo = c("mean", "observed"), ndraws = 2000, digits = 3) {
geo <- match.arg(geo)
all_contrasts <- list()
plots <- list()
for (f in model_files) {
if (!file.exists(f)) next
feature_name <- tools::file_path_sans_ext(basename(f))
feature_name <- gsub("_distance(_sc)?_sympatry_geographic_distance.*", "", feature_name)
feature_name <- gsub("_", " ", feature_name)
fit <- tryCatch(readRDS(f), error = function(e) NULL)
if (is.null(fit)) next
res <- tryCatch(
population_pair_contrasts(fit = fit, dat = dat, geo = geo, ndraws = ndraws, digits = digits),
error = function(e) NULL
)
rm(fit)
if (is.null(res) || is.null(res$contrasts)) next
res$contrasts$Feature <- feature_name
all_contrasts[[feature_name]] <- res$contrasts
plots[[feature_name]] <- res$plot + labs(title = feature_name)
}
invisible(gc(verbose = FALSE))
if (length(all_contrasts) == 0) return(NULL)
contrasts <- do.call(rbind, all_contrasts)
rownames(contrasts) <- NULL
contrasts <- contrasts[, c("Feature", "Species pair", "Sympatric population pair",
"Allopatric population pair", "Difference", "l-95% UI",
"u-95% UI", "UI excludes 0")]
list(contrasts = contrasts, plots = plots)
}
# correlation check among the song-level features that feed the PCA and
# are used as responses in the per-feature models: pairwise Pearson
# correlations among songs. Prints the pairs above the correlation
# threshold, draws a correlation matrix plot, and returns the matrix and
# pair table invisibly.
check_feature_collinearity <- function(
dat,
features,
labels = NULL,
r_threshold = 0.7
) {
X <- dat[, features, drop = FALSE]
X <- X[complete.cases(X), ]
X[] <- lapply(X, as.numeric)
if (is.null(labels)) labels <- features
names(labels) <- features
# pairwise correlations
cor_mat <- cor(X, method = "pearson")
cor_long <- as.data.frame(as.table(cor_mat), stringsAsFactors = FALSE)
names(cor_long) <- c("feature_1", "feature_2", "r")
cor_long$feature_1 <- factor(labels[cor_long$feature_1], levels = labels)
cor_long$feature_2 <- factor(labels[cor_long$feature_2], levels = labels)
# unique pairs above threshold
pair_idx <- which(upper.tri(cor_mat), arr.ind = TRUE)
pairs_df <- data.frame(
feature_1 = labels[rownames(cor_mat)[pair_idx[, 1]]],
feature_2 = labels[colnames(cor_mat)[pair_idx[, 2]]],
r = round(cor_mat[pair_idx], 3),
row.names = NULL
)
pairs_df <- pairs_df[order(-abs(pairs_df$r)), ]
high_pairs <- pairs_df[abs(pairs_df$r) >= r_threshold, ]
cat("\nPairs with |r| >=", r_threshold, ":\n")
if (nrow(high_pairs) == 0) {
cat(" none\n")
} else {
base::print(high_pairs, row.names = FALSE)
}
cat("\nStrongest pairwise correlations:\n")
base::print(head(pairs_df, 5), row.names = FALSE)
# correlation matrix plot (lower triangle), same palette as the model heatmaps
cor_long$show <- as.integer(cor_long$feature_1) > as.integer(cor_long$feature_2)
cor_long$flag <- abs(cor_long$r) >= r_threshold
p <- ggplot(cor_long[cor_long$show, ], aes(x = feature_2, y = feature_1, fill = r)) +
geom_tile(color = "white") +
geom_text(
aes(
label = sprintf("%.2f", r),
fontface = ifelse(flag, "bold", "plain"),
color = flag
),
size = 3.2
) +
scale_color_manual(values = c("TRUE" = "black", "FALSE" = "grey40"), guide = "none") +
scale_fill_gradient2(
low = "#403B78",
mid = "white",
high = "#A0DFB9CC",
midpoint = 0,
limits = c(-1, 1),
name = "Pearson r"
) +
labs(x = NULL, y = NULL) +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
# base::print - the document redefines print() as a kable wrapper
base::print(p)
invisible(list(cor = cor_mat, pairs = pairs_df, high_pairs = high_pairs))
}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)
# variables missing from stable_df (e.g. never among the top loadings for
# that PC) come back as NA from the left join; treat those as "not stable"
# rather than letting NA propagate into the plot aesthetics below
pca_rot_stck$stable[is.na(pca_rot_stck$stable)] <- FALSE
# fully opaque bars for stable variables, more transparent for the rest
# (previously ifelse(stable, 1, 1) - always 1 regardless of stability)
pca_rot_stck$top_vars <- ifelse(pca_rot_stck$stable, 1, 0.4)
# Build colored labels per row
# escape characters that markdown/HTML treats specially (gridtext, used by
# element_markdown() below, parses these labels as markdown-with-HTML; an
# unescaped &, <, > or an odd number of _ / * in a variable name can break
# that parser and silently fail to render the whole plot)
escape_markdown <- function(x) {
x <- gsub("&", "&", x, fixed = TRUE)
x <- gsub("<", "<", x, fixed = TRUE)
x <- gsub(">", ">", x, fixed = TRUE)
x <- gsub("_", "_", x, fixed = TRUE)
x <- gsub("*", "*", x, fixed = TRUE)
x
}
# Colored labels: black for stable variables, gray for the rest
# (previously compared stable < 1.1, which is true whether stable is 0 or
# 1, so every label always took the black branch)
pca_rot_stck$label_col <- ifelse(
pca_rot_stck$stable,
paste0("<span style='color:black;'>", escape_markdown(pca_rot_stck$variable), "</span>"),
paste0("<span style='color:gray50;'>", escape_markdown(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
# alpha is used directly as a numeric value (scale_alpha_identity) rather
# than mapped through as.factor()/scale_alpha_manual(); the previous version
# passed the whole top_vars column as "values" to scale_alpha_manual, which
# expects one value per factor level, not per row - fragile even before the
# always-1 bug above was fixed
ggplot(pca_rot_stck, aes(x = var_facet, y = rotation, fill = Sign,
alpha = top_vars)) + geom_col() + coord_flip() + scale_alpha_identity(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 Collinearity among song-level features
Before running the analyses we checked whether any of the seven song-level features are redundant with each other. The features enter the analyses in two ways: jointly, in the PCA that defines the multivariate acoustic distance (where correlation among them is expected and is precisely what the PCA summarizes), and separately, as the response of one model each in the feature-based analysis. In the second case, two highly correlated features would yield two near-identical models rather than two independent pieces of evidence, so strongly correlated pairs (|r| ≥ 0.7) should be reduced to a single representative before fitting. The correlation matrix plot shows pairwise Pearson correlations among songs, with pairs above the threshold in bold; the table lists the strongest pairs.
Code
song_features <- c(
"meanpeakf", "num.elms", "song.duration", "song.rate",
"gap.duration", "freq.range.Min5toMax95", "mst"
)
song_feature_labels <- c(
"Peak frequency", "Number of elements", "Song duration", "Song rate",
"Gap duration", "Frequency range", "Element diversity (mst)"
)
collin_simple <- check_feature_collinearity(
dat = simple_song_dat,
features = song_features,
labels = song_feature_labels
)
Pairs with |r| >= 0.7 :
feature_1 feature_2 r
Song duration Song rate -0.800
Frequency range Element diversity (mst) 0.795
Number of elements Song duration 0.779
Strongest pairwise correlations:
feature_1 feature_2 r
Song duration Song rate -0.800
Frequency range Element diversity (mst) 0.795
Number of elements Song duration 0.779
Song rate Gap duration -0.643
Number of elements Element diversity (mst) 0.562
2.2.3 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)
# variables missing from stable_df (e.g. never among the top loadings for
# that PC) come back as NA from the left join; treat those as "not stable"
# rather than letting NA propagate into the plot aesthetics below
pca_rot_stck$stable[is.na(pca_rot_stck$stable)] <- FALSE
# fully opaque bars for stable variables, more transparent for the rest
# (previously ifelse(stable, 1, 1) - always 1 regardless of stability)
pca_rot_stck$top_vars <- ifelse(pca_rot_stck$stable, 1, 0.4)
# Build colored labels per row
# escape characters that markdown/HTML treats specially (gridtext, used by
# element_markdown() below, parses these labels as markdown-with-HTML; an
# unescaped &, <, > or an odd number of _ / * in a variable name can break
# that parser and silently fail to render the whole plot)
escape_markdown <- function(x) {
x <- gsub("&", "&", x, fixed = TRUE)
x <- gsub("<", "<", x, fixed = TRUE)
x <- gsub(">", ">", x, fixed = TRUE)
x <- gsub("_", "_", x, fixed = TRUE)
x <- gsub("*", "*", x, fixed = TRUE)
x
}
# Colored labels: black for stable variables, gray for the rest
# (previously compared stable < 1.1, which is true whether stable is 0 or
# 1, so every label always took the black branch)
pca_rot_stck$label_col <- ifelse(
pca_rot_stck$stable,
paste0("<span style='color:black;'>", escape_markdown(pca_rot_stck$variable), "</span>"),
paste0("<span style='color:gray50;'>", escape_markdown(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
# alpha is used directly as a numeric value (scale_alpha_identity) rather
# than mapped through as.factor()/scale_alpha_manual(); the previous version
# passed the whole top_vars column as "values" to scale_alpha_manual, which
# expects one value per factor level, not per row - fragile even before the
# always-1 bug above was fixed
ggplot(pca_rot_stck, aes(x = var_facet, y = rotation, fill = Sign,
alpha = top_vars)) + geom_col() + coord_flip() + scale_alpha_identity(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.4 Acoustic and geographic distances
2.2.4.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.4.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.5 Distance and sympatry
Code
ggplot(dist_acoustic_long, aes(x = geo_distance + geo_const, 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()
2.2.6 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.7 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.7.1 Model fitting
2.2.7.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.7.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
# options(brms.file_refit = "on_change")
# 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.7.1.1.2 Prior predictive check
Before looking at the posterior we check what the priors alone imply for the response. The model below is refitted with sample_prior = "only", so the likelihood is ignored and the draws come purely from the priors. For a standardized response (mean 0, SD 1) the prior predictive distribution should comfortably cover the observed range (roughly -3 to 3) without being absurdly wide; a prior that could only generate values within the observed range would be doing the data’s job, and one spanning hundreds of SDs would not be regularizing anything.
Code
prior_only_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(),
prior = priors,
sample_prior = "only",
chains = 2,
iter = 2000,
cores = 2,
seed = 123,
backend = "cmdstanr",
file = "./data/processed/fits/prior_only_acoustic_distance_simple"
)
pp_check(prior_only_simple, type = "dens_overlay", ndraws = 50) +
coord_cartesian(xlim = c(-20, 20)) +
labs(
title = "Prior predictive check (simple songs, multivariate model)",
subtitle = "y: observed standardized acoustic distance; y_rep: draws from the priors only"
)
2.2.7.1.1.3 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 + geo_distance_sc + (1 | species_pair) + (1 | mm(pop1, pop2)) + (1 | mm(individual1, individual2)) | 10000 | 4 | 1 | 5000 | 448 (0.022%) | 0 | 29004.5 | 12502.6 | 489747410 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sympatry1 | 0.421 | 0.364 | 0.477 | 1 | 29004.50 | 12502.60 |
| b_geo_distance_sc | 0.210 | 0.154 | 0.266 | 1 | 29650.36 | 13365.59 |
2.2.7.1.1.4 Posterior predictive check
Posterior predictive checks compare the observed data with data simulated from the fitted model (Gelman et al. 2014). Three views: the distribution of the response versus 50 replicated data sets (density overlay); the mean of the response within sympatric and allopatric pairs versus its posterior predictive distribution (the quantity the sympatry coefficient is about); and observed values against their posterior-averaged predictions. The observed data should fall within the range of the replicates.
Code
pp_check(sympatry_geo_model_simple, type = "dens_overlay", ndraws = 50) +
labs(title = "Posterior predictive check (simple songs): density overlay")
Code
pp_check(sympatry_geo_model_simple, type = "stat_grouped", stat = "mean", group = "sympatry", ndraws = 500) +
labs(title = "Posterior predictive check (simple songs): mean by sympatry (0 = allopatric, 1 = sympatric)")
Code
pp_check(sympatry_geo_model_simple, type = "scatter_avg", ndraws = 100) +
labs(title = "Posterior predictive check (simple songs): observed vs. posterior-averaged prediction")
Summary
Both predictors had credible, positive effects on overall acoustic distance: sympatry (β = 0.421, 95% CI [0.364, 0.477]) and geographic distance (β = 0.210, 95% CI [0.154, 0.266]). 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 2x 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.7.1.1.5 Contrasts by species pair
The model above estimates a single sympatry effect for all species pairs; species pair enters only as a varying intercept, which shifts sympatric and allopatric comparisons by the same amount and therefore cannot tell whether the effect is shared by all pairs or driven by one of them. To obtain a separate sympatry effect for each species pair, the model is refitted with species-pair-specific sympatry slopes (a varying slope of sympatry by species pair, correlated with the varying intercept, with an LKJ(2) prior on the correlation). The per-pair effect is the fixed sympatry effect plus the pair’s deviation; the pairwise differences between species pairs are computed from the posterior draws. With only three species pairs the between-pair standard deviation of the sympatry slope, and its correlation with the intercept, are weakly informed by the data and lean on the priors; the per-pair effects are partially pooled towards the overall effect, so they should be read as shrunken estimates rather than as three independent fits. The fitting chunk is not evaluated when rendering (it takes as long as the main model); run it once and the results are loaded from the saved fit.
Code
priors_sp <- c(
priors,
prior(lkj(2), class = "cor")
)
sympatry_by_species_pair_simple <- brm(
acoustic_distance_sc ~ sympatry +
geo_distance_sc +
(1 + sympatry | species_pair) +
(1 | mm(pop1, pop2)) +
(1 | mm(individual1, individual2)),
data = sympatric_pairs_simple,
family = gaussian(),
prior = priors_sp,
cores = 4,
chains = 4,
iter = 10000,
backend = "cmdstanr",
threads = threading(8),
control = list(adapt_delta = 0.95, max_treedepth = 15),
file = "./data/processed/fits/acoustic_distance_sympatry_by_species_pair_simple_fit"
)Code
sp_fit_file <- "./data/processed/fits/acoustic_distance_sympatry_by_species_pair_simple_fit.rds"
if (file.exists(sp_fit_file)) {
sympatry_by_species_pair_simple <- readRDS(sp_fit_file)
sp_contrasts_simple <- species_pair_sympatry_effects(sympatry_by_species_pair_simple)
base::print(sp_contrasts_simple$plot)
} else {
message("Species-pair slope model not fitted yet - run the chunk above first.")
}Code
if (exists("sp_contrasts_simple")) print(sp_contrasts_simple$effects)Code
if (exists("sp_contrasts_simple")) print(sp_contrasts_simple$differences)Code
if (exists("sympatry_by_species_pair_simple")) {
extended_summary(
sympatry_by_species_pair_simple,
highlight = TRUE,
trace.palette = viridis::mako,
remove.intercepts = TRUE,
print.name = FALSE
)
}Population-pair contrasts from the species-pair slope model (same procedure as in the next section, but each species pair now has its own sympatry effect):
Code
if (exists("sympatry_by_species_pair_simple")) {
pop_contrasts_sp_simple <- population_pair_contrasts(
fit = sympatry_by_species_pair_simple,
dat = sympatric_pairs_simple,
geo = "mean"
)
base::print(pop_contrasts_sp_simple$plot)
}Code
if (exists("pop_contrasts_sp_simple")) print(pop_contrasts_sp_simple$contrasts)2.2.7.1.1.6 Contrasts among population pairs
Using the fitted model as is, the expected acoustic distance is computed for every observed combination of populations within each species pair, including the species-pair and population (multi-membership) effects and setting the individual effects to zero (i.e. for an average individual). Sympatric population pairs (the two species recorded at the same locality) are then contrasted against each allopatric population pair of the same species pair; each contrast is computed draw by draw, so its uncertainty interval reflects the joint posterior. Contrasts that share a population pair (e.g. the same sympatric pair against several allopatric ones) are not independent of each other and should be read as a set rather than as separate tests. Geographic distance is held at the pooled mean (geo_distance_sc = 0) so that the contrasts reflect sympatry and population effects rather than the distance between localities; set geo = "observed" to use each population pair’s own mean distance instead.
Code
pop_contrasts_simple <- population_pair_contrasts(
fit = sympatry_geo_model_simple,
dat = sympatric_pairs_simple,
geo = "mean"
)
base::print(pop_contrasts_simple$plot)
Code
print(pop_contrasts_simple$expected)| Species pair | Population pair | Sympatry | N comparisons | Geographic distance (sc) | Expected distance | l-95% UI | u-95% UI |
|---|---|---|---|---|---|---|---|
| Hypoxantha_Iberaensis | Entre_Rios vs E_Ibera | Allopatric | 774 | 0 | -0.308 | -0.754 | 0.092 |
| Hypoxantha_Iberaensis | Mar_Chiquita vs E_Ibera | Allopatric | 576 | 0 | -0.142 | -0.596 | 0.326 |
| Hypoxantha_Palustris | Entre_Rios vs E_Ibera | Allopatric | 1290 | 0 | -0.179 | -0.509 | 0.143 |
| Hypoxantha_Palustris | Mar_Chiquita vs E_Ibera | Allopatric | 960 | 0 | -0.013 | -0.384 | 0.420 |
| Hypoxantha_Ruficollis | E_Ibera vs Entre_Rios | Allopatric | 8816 | 0 | -0.170 | -0.455 | 0.090 |
| Hypoxantha_Ruficollis | Mar_Chiquita vs Entre_Rios | Allopatric | 2432 | 0 | -0.067 | -0.415 | 0.299 |
| Hypoxantha_Ruficollis | E_Ibera vs Esperanza | Allopatric | 5336 | 0 | -0.141 | -0.473 | 0.180 |
| Hypoxantha_Ruficollis | Entre_Rios vs Esperanza | Allopatric | 1978 | 0 | -0.204 | -0.616 | 0.146 |
| Hypoxantha_Ruficollis | Mar_Chiquita vs Esperanza | Allopatric | 1472 | 0 | -0.038 | -0.425 | 0.388 |
| Hypoxantha_Ruficollis | E_Ibera vs Mar_Chiquita | Allopatric | 3248 | 0 | -0.181 | -0.507 | 0.126 |
| Hypoxantha_Ruficollis | Entre_Rios vs Mar_Chiquita | Allopatric | 1204 | 0 | -0.243 | -0.643 | 0.100 |
| Hypoxantha_Ruficollis | E_Ibera vs Salta | Allopatric | 3480 | 0 | 0.143 | -0.254 | 0.584 |
| Hypoxantha_Ruficollis | Entre_Rios vs Salta | Allopatric | 1290 | 0 | 0.081 | -0.294 | 0.536 |
| Hypoxantha_Ruficollis | Mar_Chiquita vs Salta | Allopatric | 960 | 0 | 0.246 | -0.235 | 0.844 |
| Hypoxantha_Iberaensis | E_Ibera vs E_Ibera | Sympatric | 2088 | 0 | 0.175 | -0.216 | 0.538 |
| Hypoxantha_Palustris | E_Ibera vs E_Ibera | Sympatric | 3480 | 0 | 0.304 | 0.015 | 0.583 |
| Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Sympatric | 3268 | 0 | 0.188 | -0.152 | 0.498 |
| Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Sympatric | 896 | 0 | 0.343 | -0.031 | 0.752 |
Code
print(pop_contrasts_simple$contrasts)| Species pair | Sympatric population pair | Allopatric population pair | Difference | l-95% UI | u-95% UI | UI excludes 0 |
|---|---|---|---|---|---|---|
| Hypoxantha_Iberaensis | E_Ibera vs E_Ibera | Entre_Rios vs E_Ibera | 0.483 | 0.215 | 0.788 | yes |
| Hypoxantha_Iberaensis | E_Ibera vs E_Ibera | Mar_Chiquita vs E_Ibera | 0.317 | -0.102 | 0.613 | no |
| Hypoxantha_Palustris | E_Ibera vs E_Ibera | Entre_Rios vs E_Ibera | 0.483 | 0.215 | 0.788 | yes |
| Hypoxantha_Palustris | E_Ibera vs E_Ibera | Mar_Chiquita vs E_Ibera | 0.317 | -0.102 | 0.613 | no |
| Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Entre_Rios | 0.358 | 0.042 | 0.644 | yes |
| Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Mar_Chiquita vs Entre_Rios | 0.255 | -0.178 | 0.568 | no |
| Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Esperanza | 0.329 | -0.123 | 0.744 | no |
| Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Entre_Rios vs Esperanza | 0.391 | 0.076 | 0.697 | yes |
| Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Mar_Chiquita vs Esperanza | 0.226 | -0.351 | 0.675 | no |
| Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Mar_Chiquita | 0.369 | -0.045 | 0.777 | no |
| Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Entre_Rios vs Mar_Chiquita | 0.431 | 0.137 | 0.728 | yes |
| Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Salta | 0.045 | -0.571 | 0.486 | no |
| Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Entre_Rios vs Salta | 0.107 | -0.370 | 0.462 | no |
| Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Mar_Chiquita vs Salta | -0.059 | -0.796 | 0.461 | no |
| Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Entre_Rios | 0.513 | 0.095 | 0.989 | yes |
| Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Mar_Chiquita vs Entre_Rios | 0.410 | 0.102 | 0.694 | yes |
| Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Esperanza | 0.484 | 0.022 | 0.997 | yes |
| Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Entre_Rios vs Esperanza | 0.546 | 0.090 | 1.092 | yes |
| Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Mar_Chiquita vs Esperanza | 0.381 | 0.020 | 0.716 | yes |
| Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Mar_Chiquita | 0.524 | 0.232 | 0.931 | yes |
| Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Entre_Rios vs Mar_Chiquita | 0.586 | 0.276 | 1.024 | yes |
| Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Salta | 0.200 | -0.364 | 0.634 | no |
| Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Entre_Rios vs Salta | 0.262 | -0.325 | 0.699 | no |
| Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Mar_Chiquita vs Salta | 0.096 | -0.433 | 0.467 | no |
2.2.7.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
# prior predictive check for the trait-level prior set (defined again here
# because the fitting chunk below is not evaluated when rendering), using
# one representative standardized response; the same priors and formula are
# used for every trait, so one check covers the set
priors_trait <- c(
prior(normal(0, 1), class = "Intercept"),
prior(normal(0, 0.5), class = "b"),
prior(exponential(2), class = "sd"),
prior(exponential(1), class = "sigma")
)
sympatric_pairs_simple$freqrange_distance_sc <- as.numeric(scale(sympatric_pairs_simple$freqrange_distance))
prior_only_trait_simple <- brm(
freqrange_distance_sc ~ sympatry +
geo_distance_sc +
(1 | species_pair) +
(1 | mm(pop1, pop2)) +
(1 | mm(individual1, individual2)),
data = sympatric_pairs_simple,
family = gaussian(),
prior = priors_trait,
sample_prior = "only",
chains = 2,
iter = 2000,
cores = 2,
seed = 123,
backend = "cmdstanr",
file = "./data/processed/fits/prior_only_trait_level_simple"
)
pp_check(prior_only_trait_simple, type = "dens_overlay", ndraws = 50) +
coord_cartesian(xlim = c(-10, 10)) +
labs(
title = "Prior predictive check (simple songs, trait-level priors)",
subtitle = "y: observed standardized frequency-range distance; y_rep: draws from the priors only"
)
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 rev(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 = 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 table below gives, for every trait-level model above, the posterior mean, standard error and 95% uncertainty interval of each predictor; the last column flags intervals that exclude zero. Like the heatmap, it only includes models whose fit file exists.
Code
fixef_tab_simple <- fixef_table(model_files)
if (!is.null(fixef_tab_simple)) print(fixef_tab_simple)| Response | Predictor | Estimate | SE | l-95% UI | u-95% UI | UI excludes 0 |
|---|---|---|---|---|---|---|
| elm_pca | Geographic distance | 0.002 | 0.025 | -0.048 | 0.051 | no |
| elm_pca | Sympatry | 0.394 | 0.026 | 0.344 | 0.444 | yes |
| freqrange | Geographic distance | -0.093 | 0.029 | -0.151 | -0.035 | yes |
| freqrange | Sympatry | 0.384 | 0.030 | 0.326 | 0.442 | yes |
| gapduration | Geographic distance | -0.031 | 0.028 | -0.087 | 0.025 | no |
| gapduration | Sympatry | -0.039 | 0.029 | -0.095 | 0.018 | no |
| meanpeakf | Geographic distance | 0.031 | 0.027 | -0.021 | 0.084 | no |
| meanpeakf | Sympatry | -0.038 | 0.027 | -0.091 | 0.014 | no |
| mst | Geographic distance | -0.018 | 0.029 | -0.074 | 0.038 | no |
| mst | Sympatry | 0.539 | 0.029 | 0.482 | 0.596 | yes |
| numelms | Geographic distance | 0.087 | 0.028 | 0.033 | 0.141 | yes |
| numelms | Sympatry | 0.125 | 0.028 | 0.070 | 0.180 | yes |
| songduration | Geographic distance | 0.075 | 0.029 | 0.019 | 0.132 | yes |
| songduration | Sympatry | 0.114 | 0.029 | 0.057 | 0.171 | yes |
| songrate | Geographic distance | 0.025 | 0.023 | -0.021 | 0.071 | no |
| songrate | Sympatry | -0.017 | 0.024 | -0.064 | 0.030 | no |
The expandable section below provides the complete Bayesian model summaries underlying those estimates, including posterior parameter estimates, uncertainty 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")
}The expandable section below shows a posterior predictive check (density overlay, 50 replicated data sets) for every trait-level model that has finished fitting; models whose fit file is missing are skipped.
NotePosterior predictive checks
Code
for (i in model_files) {
fit <- tryCatch(readRDS(i), error = function(e) NULL)
if (is.null(fit)) next
model_name <- gsub("_sympatry_geographic_distance.*", "", tools::file_path_sans_ext(basename(i)))
cat("\n\n**", model_name, "**\n\n", sep = "")
# base::print - the document redefines print() as a kable wrapper
base::print(
pp_check(fit, type = "dens_overlay", ndraws = 50) +
labs(title = model_name)
)
cat("\n\n")
rm(fit)
invisible(gc(verbose = FALSE))
}elm_pca_distance_sc
freqrange_distance_sc
gapduration_distance_sc
meanpeakf_distance_sc
mst_distance_sc
numelms_distance_sc
songduration_distance_sc
songrate_distance_sc
2.2.7.1.2.1 Contrasts among population pairs, per feature
The same population-pair contrasts computed for the overall PCA-based model above, now repeated for every individual acoustic feature that has finished fitting: expected acoustic distance for each observed population pair (species-pair and population multi-membership effects included, individual effects set to zero, geographic distance held at the pooled mean), then sympatric vs allopatric population-pair differences within each species pair. Each feature gets its own plot (faceted by species pair, as in the main-model version above) rather than being crowded into one combined figure. As with the main-model version, contrasts sharing a population pair are not independent of each other.
Code
pop_contrasts_by_feature_simple <- population_pair_contrasts_by_feature(
model_files = model_files,
dat = sympatric_pairs_simple,
geo = "mean"
)
if (!is.null(pop_contrasts_by_feature_simple)) {
for (feat in names(pop_contrasts_by_feature_simple$plots)) {
cat("\n\n**", feat, "**\n\n", sep = "")
base::print(pop_contrasts_by_feature_simple$plots[[feat]])
cat("\n\n")
}
}elm pca
freqrange
gapduration
meanpeakf
mst
numelms
songduration
songrate
Code
if (!is.null(pop_contrasts_by_feature_simple)) print(pop_contrasts_by_feature_simple$contrasts)| Feature | Species pair | Sympatric population pair | Allopatric population pair | Difference | l-95% UI | u-95% UI | UI excludes 0 |
|---|---|---|---|---|---|---|---|
| elm pca | Hypoxantha_Iberaensis | E_Ibera vs E_Ibera | Entre_Rios vs E_Ibera | 0.471 | 0.216 | 0.794 | yes |
| elm pca | Hypoxantha_Iberaensis | E_Ibera vs E_Ibera | Mar_Chiquita vs E_Ibera | 0.675 | 0.336 | 1.197 | yes |
| elm pca | Hypoxantha_Palustris | E_Ibera vs E_Ibera | Entre_Rios vs E_Ibera | 0.471 | 0.216 | 0.794 | yes |
| elm pca | Hypoxantha_Palustris | E_Ibera vs E_Ibera | Mar_Chiquita vs E_Ibera | 0.675 | 0.336 | 1.197 | yes |
| elm pca | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Entre_Rios | 0.317 | -0.010 | 0.584 | no |
| elm pca | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Mar_Chiquita vs Entre_Rios | 0.599 | 0.278 | 1.104 | yes |
| elm pca | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Esperanza | 0.395 | -0.039 | 0.834 | no |
| elm pca | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Entre_Rios vs Esperanza | 0.472 | 0.175 | 0.853 | yes |
| elm pca | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Mar_Chiquita vs Esperanza | 0.676 | 0.243 | 1.367 | yes |
| elm pca | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Mar_Chiquita | 0.221 | -0.276 | 0.604 | no |
| elm pca | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Entre_Rios vs Mar_Chiquita | 0.298 | -0.074 | 0.578 | no |
| elm pca | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Salta | 0.300 | -0.159 | 0.699 | no |
| elm pca | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Entre_Rios vs Salta | 0.377 | 0.047 | 0.683 | yes |
| elm pca | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Mar_Chiquita vs Salta | 0.582 | 0.161 | 1.209 | yes |
| elm pca | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Entre_Rios | 0.208 | -0.363 | 0.638 | no |
| elm pca | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Mar_Chiquita vs Entre_Rios | 0.490 | 0.212 | 0.854 | yes |
| elm pca | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Esperanza | 0.286 | -0.271 | 0.741 | no |
| elm pca | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Entre_Rios vs Esperanza | 0.363 | -0.159 | 0.825 | no |
| elm pca | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Mar_Chiquita vs Esperanza | 0.567 | 0.270 | 1.047 | yes |
| elm pca | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Mar_Chiquita | 0.112 | -0.409 | 0.450 | no |
| elm pca | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Entre_Rios vs Mar_Chiquita | 0.189 | -0.309 | 0.508 | no |
| elm pca | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Salta | 0.191 | -0.405 | 0.612 | no |
| elm pca | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Entre_Rios vs Salta | 0.268 | -0.302 | 0.688 | no |
| elm pca | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Mar_Chiquita vs Salta | 0.473 | 0.161 | 0.863 | yes |
| freqrange | Hypoxantha_Iberaensis | E_Ibera vs E_Ibera | Entre_Rios vs E_Ibera | 0.418 | 0.184 | 0.671 | yes |
| freqrange | Hypoxantha_Iberaensis | E_Ibera vs E_Ibera | Mar_Chiquita vs E_Ibera | -0.011 | -0.390 | 0.386 | no |
| freqrange | Hypoxantha_Palustris | E_Ibera vs E_Ibera | Entre_Rios vs E_Ibera | 0.418 | 0.184 | 0.671 | yes |
| freqrange | Hypoxantha_Palustris | E_Ibera vs E_Ibera | Mar_Chiquita vs E_Ibera | -0.011 | -0.390 | 0.386 | no |
| freqrange | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Entre_Rios | 0.350 | 0.082 | 0.578 | yes |
| freqrange | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Mar_Chiquita vs Entre_Rios | -0.045 | -0.483 | 0.386 | no |
| freqrange | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Esperanza | 0.318 | -0.050 | 0.638 | no |
| freqrange | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Entre_Rios vs Esperanza | 0.352 | 0.090 | 0.626 | yes |
| freqrange | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Mar_Chiquita vs Esperanza | -0.077 | -0.610 | 0.399 | no |
| freqrange | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Mar_Chiquita | 0.229 | -0.149 | 0.575 | no |
| freqrange | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Entre_Rios vs Mar_Chiquita | 0.263 | -0.026 | 0.517 | no |
| freqrange | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Salta | 0.149 | -0.262 | 0.493 | no |
| freqrange | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Entre_Rios vs Salta | 0.183 | -0.126 | 0.453 | no |
| freqrange | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Mar_Chiquita vs Salta | -0.245 | -0.792 | 0.378 | no |
| freqrange | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Entre_Rios | 0.899 | 0.384 | 1.416 | yes |
| freqrange | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Mar_Chiquita vs Entre_Rios | 0.505 | 0.242 | 0.794 | yes |
| freqrange | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Esperanza | 0.867 | 0.370 | 1.396 | yes |
| freqrange | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Entre_Rios vs Esperanza | 0.902 | 0.368 | 1.470 | yes |
| freqrange | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Mar_Chiquita vs Esperanza | 0.473 | 0.175 | 0.791 | yes |
| freqrange | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Mar_Chiquita | 0.778 | 0.383 | 1.152 | yes |
| freqrange | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Entre_Rios vs Mar_Chiquita | 0.813 | 0.387 | 1.247 | yes |
| freqrange | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Salta | 0.699 | 0.279 | 1.204 | yes |
| freqrange | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Entre_Rios vs Salta | 0.733 | 0.284 | 1.265 | yes |
| freqrange | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Mar_Chiquita vs Salta | 0.304 | -0.035 | 0.608 | no |
| gapduration | Hypoxantha_Iberaensis | E_Ibera vs E_Ibera | Entre_Rios vs E_Ibera | 0.015 | -0.197 | 0.340 | no |
| gapduration | Hypoxantha_Iberaensis | E_Ibera vs E_Ibera | Mar_Chiquita vs E_Ibera | -0.072 | -0.413 | 0.160 | no |
| gapduration | Hypoxantha_Palustris | E_Ibera vs E_Ibera | Entre_Rios vs E_Ibera | 0.015 | -0.197 | 0.340 | no |
| gapduration | Hypoxantha_Palustris | E_Ibera vs E_Ibera | Mar_Chiquita vs E_Ibera | -0.072 | -0.413 | 0.160 | no |
| gapduration | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Entre_Rios | -0.092 | -0.405 | 0.115 | no |
| gapduration | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Mar_Chiquita vs Entre_Rios | -0.126 | -0.608 | 0.100 | no |
| gapduration | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Esperanza | -0.080 | -0.487 | 0.239 | no |
| gapduration | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Entre_Rios vs Esperanza | -0.026 | -0.268 | 0.259 | no |
| gapduration | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Mar_Chiquita vs Esperanza | -0.113 | -0.621 | 0.186 | no |
| gapduration | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Mar_Chiquita | -0.088 | -0.512 | 0.212 | no |
| gapduration | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Entre_Rios vs Mar_Chiquita | -0.035 | -0.303 | 0.233 | no |
| gapduration | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Salta | -0.097 | -0.545 | 0.182 | no |
| gapduration | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Entre_Rios vs Salta | -0.044 | -0.311 | 0.200 | no |
| gapduration | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Mar_Chiquita vs Salta | -0.131 | -0.696 | 0.163 | no |
| gapduration | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Entre_Rios | -0.008 | -0.335 | 0.423 | no |
| gapduration | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Mar_Chiquita vs Entre_Rios | -0.042 | -0.313 | 0.232 | no |
| gapduration | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Esperanza | 0.004 | -0.316 | 0.474 | no |
| gapduration | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Entre_Rios vs Esperanza | 0.058 | -0.242 | 0.630 | no |
| gapduration | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Mar_Chiquita vs Esperanza | -0.030 | -0.286 | 0.288 | no |
| gapduration | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Mar_Chiquita | -0.004 | -0.245 | 0.344 | no |
| gapduration | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Entre_Rios vs Mar_Chiquita | 0.049 | -0.180 | 0.518 | no |
| gapduration | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Salta | -0.013 | -0.361 | 0.397 | no |
| gapduration | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Entre_Rios vs Salta | 0.040 | -0.273 | 0.575 | no |
| gapduration | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Mar_Chiquita vs Salta | -0.047 | -0.339 | 0.221 | no |
| meanpeakf | Hypoxantha_Iberaensis | E_Ibera vs E_Ibera | Entre_Rios vs E_Ibera | 0.065 | -0.172 | 0.441 | no |
| meanpeakf | Hypoxantha_Iberaensis | E_Ibera vs E_Ibera | Mar_Chiquita vs E_Ibera | 0.025 | -0.257 | 0.421 | no |
| meanpeakf | Hypoxantha_Palustris | E_Ibera vs E_Ibera | Entre_Rios vs E_Ibera | 0.065 | -0.172 | 0.441 | no |
| meanpeakf | Hypoxantha_Palustris | E_Ibera vs E_Ibera | Mar_Chiquita vs E_Ibera | 0.025 | -0.257 | 0.421 | no |
| meanpeakf | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Entre_Rios | -0.142 | -0.508 | 0.098 | no |
| meanpeakf | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Mar_Chiquita vs Entre_Rios | -0.079 | -0.475 | 0.232 | no |
| meanpeakf | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Esperanza | -0.014 | -0.403 | 0.459 | no |
| meanpeakf | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Entre_Rios vs Esperanza | 0.090 | -0.164 | 0.536 | no |
| meanpeakf | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Mar_Chiquita vs Esperanza | 0.050 | -0.335 | 0.614 | no |
| meanpeakf | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Mar_Chiquita | -0.147 | -0.659 | 0.210 | no |
| meanpeakf | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Entre_Rios vs Mar_Chiquita | -0.043 | -0.355 | 0.262 | no |
| meanpeakf | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Salta | -0.127 | -0.631 | 0.245 | no |
| meanpeakf | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Entre_Rios vs Salta | -0.024 | -0.323 | 0.292 | no |
| meanpeakf | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Mar_Chiquita vs Salta | -0.063 | -0.545 | 0.386 | no |
| meanpeakf | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Entre_Rios | -0.098 | -0.618 | 0.314 | no |
| meanpeakf | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Mar_Chiquita vs Entre_Rios | -0.034 | -0.330 | 0.282 | no |
| meanpeakf | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Esperanza | 0.030 | -0.414 | 0.598 | no |
| meanpeakf | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Entre_Rios vs Esperanza | 0.134 | -0.270 | 0.803 | no |
| meanpeakf | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Mar_Chiquita vs Esperanza | 0.094 | -0.174 | 0.583 | no |
| meanpeakf | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Mar_Chiquita | -0.103 | -0.490 | 0.184 | no |
| meanpeakf | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Entre_Rios vs Mar_Chiquita | 0.001 | -0.313 | 0.376 | no |
| meanpeakf | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Salta | -0.083 | -0.611 | 0.353 | no |
| meanpeakf | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Entre_Rios vs Salta | 0.021 | -0.411 | 0.559 | no |
| meanpeakf | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Mar_Chiquita vs Salta | -0.019 | -0.332 | 0.318 | no |
| mst | Hypoxantha_Iberaensis | E_Ibera vs E_Ibera | Entre_Rios vs E_Ibera | 0.417 | 0.122 | 0.637 | yes |
| mst | Hypoxantha_Iberaensis | E_Ibera vs E_Ibera | Mar_Chiquita vs E_Ibera | 0.270 | -0.183 | 0.582 | no |
| mst | Hypoxantha_Palustris | E_Ibera vs E_Ibera | Entre_Rios vs E_Ibera | 0.417 | 0.122 | 0.637 | yes |
| mst | Hypoxantha_Palustris | E_Ibera vs E_Ibera | Mar_Chiquita vs E_Ibera | 0.270 | -0.183 | 0.582 | no |
| mst | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Entre_Rios | 0.662 | 0.450 | 0.953 | yes |
| mst | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Mar_Chiquita vs Entre_Rios | 0.392 | -0.005 | 0.662 | no |
| mst | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Esperanza | 0.645 | 0.344 | 1.032 | yes |
| mst | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Entre_Rios vs Esperanza | 0.522 | 0.267 | 0.771 | yes |
| mst | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Mar_Chiquita vs Esperanza | 0.375 | -0.121 | 0.733 | no |
| mst | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Mar_Chiquita | 0.619 | 0.293 | 0.980 | yes |
| mst | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Entre_Rios vs Mar_Chiquita | 0.497 | 0.227 | 0.739 | yes |
| mst | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Salta | 0.638 | 0.328 | 1.031 | yes |
| mst | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Entre_Rios vs Salta | 0.515 | 0.247 | 0.773 | yes |
| mst | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Mar_Chiquita vs Salta | 0.368 | -0.142 | 0.709 | no |
| mst | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Entre_Rios | 0.852 | 0.484 | 1.399 | yes |
| mst | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Mar_Chiquita vs Entre_Rios | 0.582 | 0.352 | 0.845 | yes |
| mst | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Esperanza | 0.835 | 0.476 | 1.407 | yes |
| mst | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Entre_Rios vs Esperanza | 0.713 | 0.368 | 1.216 | yes |
| mst | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Mar_Chiquita vs Esperanza | 0.565 | 0.296 | 0.858 | yes |
| mst | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Mar_Chiquita | 0.809 | 0.499 | 1.268 | yes |
| mst | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Entre_Rios vs Mar_Chiquita | 0.687 | 0.427 | 1.091 | yes |
| mst | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Salta | 0.828 | 0.467 | 1.397 | yes |
| mst | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Entre_Rios vs Salta | 0.706 | 0.356 | 1.216 | yes |
| mst | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Mar_Chiquita vs Salta | 0.558 | 0.293 | 0.837 | yes |
| numelms | Hypoxantha_Iberaensis | E_Ibera vs E_Ibera | Entre_Rios vs E_Ibera | 0.074 | -0.197 | 0.247 | no |
| numelms | Hypoxantha_Iberaensis | E_Ibera vs E_Ibera | Mar_Chiquita vs E_Ibera | 0.103 | -0.138 | 0.295 | no |
| numelms | Hypoxantha_Palustris | E_Ibera vs E_Ibera | Entre_Rios vs E_Ibera | 0.074 | -0.197 | 0.247 | no |
| numelms | Hypoxantha_Palustris | E_Ibera vs E_Ibera | Mar_Chiquita vs E_Ibera | 0.103 | -0.138 | 0.295 | no |
| numelms | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Entre_Rios | 0.176 | 0.009 | 0.438 | yes |
| numelms | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Mar_Chiquita vs Entre_Rios | 0.154 | -0.056 | 0.423 | no |
| numelms | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Esperanza | 0.159 | -0.118 | 0.472 | no |
| numelms | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Entre_Rios vs Esperanza | 0.108 | -0.124 | 0.316 | no |
| numelms | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Mar_Chiquita vs Esperanza | 0.137 | -0.181 | 0.475 | no |
| numelms | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Mar_Chiquita | 0.214 | -0.016 | 0.607 | no |
| numelms | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Entre_Rios vs Mar_Chiquita | 0.164 | -0.024 | 0.437 | no |
| numelms | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Salta | 0.181 | -0.060 | 0.552 | no |
| numelms | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Entre_Rios vs Salta | 0.130 | -0.082 | 0.351 | no |
| numelms | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Mar_Chiquita vs Salta | 0.159 | -0.115 | 0.537 | no |
| numelms | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Entre_Rios | 0.108 | -0.239 | 0.393 | no |
| numelms | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Mar_Chiquita vs Entre_Rios | 0.087 | -0.198 | 0.273 | no |
| numelms | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Esperanza | 0.091 | -0.289 | 0.371 | no |
| numelms | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Entre_Rios vs Esperanza | 0.040 | -0.412 | 0.293 | no |
| numelms | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Mar_Chiquita vs Esperanza | 0.069 | -0.247 | 0.251 | no |
| numelms | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Mar_Chiquita | 0.147 | -0.051 | 0.376 | no |
| numelms | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Entre_Rios vs Mar_Chiquita | 0.096 | -0.167 | 0.291 | no |
| numelms | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Salta | 0.113 | -0.234 | 0.418 | no |
| numelms | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Entre_Rios vs Salta | 0.063 | -0.328 | 0.331 | no |
| numelms | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Mar_Chiquita vs Salta | 0.092 | -0.187 | 0.295 | no |
| songduration | Hypoxantha_Iberaensis | E_Ibera vs E_Ibera | Entre_Rios vs E_Ibera | 0.052 | -0.251 | 0.296 | no |
| songduration | Hypoxantha_Iberaensis | E_Ibera vs E_Ibera | Mar_Chiquita vs E_Ibera | -0.062 | -0.492 | 0.209 | no |
| songduration | Hypoxantha_Palustris | E_Ibera vs E_Ibera | Entre_Rios vs E_Ibera | 0.052 | -0.251 | 0.296 | no |
| songduration | Hypoxantha_Palustris | E_Ibera vs E_Ibera | Mar_Chiquita vs E_Ibera | -0.062 | -0.492 | 0.209 | no |
| songduration | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Entre_Rios | 0.175 | -0.075 | 0.465 | no |
| songduration | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Mar_Chiquita vs Entre_Rios | 0.000 | -0.403 | 0.290 | no |
| songduration | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Esperanza | 0.114 | -0.276 | 0.504 | no |
| songduration | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Entre_Rios vs Esperanza | 0.053 | -0.265 | 0.324 | no |
| songduration | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Mar_Chiquita vs Esperanza | -0.061 | -0.617 | 0.340 | no |
| songduration | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Mar_Chiquita | 0.169 | -0.193 | 0.549 | no |
| songduration | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Entre_Rios vs Mar_Chiquita | 0.107 | -0.165 | 0.384 | no |
| songduration | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Salta | -0.055 | -0.603 | 0.288 | no |
| songduration | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Entre_Rios vs Salta | -0.117 | -0.581 | 0.176 | no |
| songduration | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Mar_Chiquita vs Salta | -0.230 | -0.947 | 0.203 | no |
| songduration | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Entre_Rios | 0.295 | -0.063 | 0.808 | no |
| songduration | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Mar_Chiquita vs Entre_Rios | 0.120 | -0.149 | 0.416 | no |
| songduration | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Esperanza | 0.235 | -0.154 | 0.738 | no |
| songduration | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Entre_Rios vs Esperanza | 0.173 | -0.232 | 0.659 | no |
| songduration | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Mar_Chiquita vs Esperanza | 0.059 | -0.288 | 0.332 | no |
| songduration | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Mar_Chiquita | 0.289 | 0.020 | 0.724 | yes |
| songduration | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Entre_Rios vs Mar_Chiquita | 0.227 | -0.059 | 0.640 | no |
| songduration | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Salta | 0.065 | -0.405 | 0.473 | no |
| songduration | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Entre_Rios vs Salta | 0.004 | -0.504 | 0.389 | no |
| songduration | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Mar_Chiquita vs Salta | -0.110 | -0.626 | 0.195 | no |
| songrate | Hypoxantha_Iberaensis | E_Ibera vs E_Ibera | Entre_Rios vs E_Ibera | 0.072 | -0.355 | 0.508 | no |
| songrate | Hypoxantha_Iberaensis | E_Ibera vs E_Ibera | Mar_Chiquita vs E_Ibera | -0.233 | -0.778 | 0.315 | no |
| songrate | Hypoxantha_Palustris | E_Ibera vs E_Ibera | Entre_Rios vs E_Ibera | 0.072 | -0.355 | 0.508 | no |
| songrate | Hypoxantha_Palustris | E_Ibera vs E_Ibera | Mar_Chiquita vs E_Ibera | -0.233 | -0.778 | 0.315 | no |
| songrate | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Entre_Rios | -0.106 | -0.537 | 0.333 | no |
| songrate | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Mar_Chiquita vs Entre_Rios | -0.322 | -0.914 | 0.227 | no |
| songrate | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Esperanza | -0.291 | -0.925 | 0.355 | no |
| songrate | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Entre_Rios vs Esperanza | -0.203 | -0.685 | 0.263 | no |
| songrate | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Mar_Chiquita vs Esperanza | -0.507 | -1.334 | 0.216 | no |
| songrate | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Mar_Chiquita | -0.222 | -0.870 | 0.406 | no |
| songrate | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Entre_Rios vs Mar_Chiquita | -0.133 | -0.598 | 0.326 | no |
| songrate | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | E_Ibera vs Salta | -1.164 | -1.927 | -0.406 | yes |
| songrate | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Entre_Rios vs Salta | -1.075 | -1.656 | -0.414 | yes |
| songrate | Hypoxantha_Ruficollis | Entre_Rios vs Entre_Rios | Mar_Chiquita vs Salta | -1.380 | -2.292 | -0.512 | yes |
| songrate | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Entre_Rios | 0.314 | -0.390 | 1.056 | no |
| songrate | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Mar_Chiquita vs Entre_Rios | 0.099 | -0.343 | 0.568 | no |
| songrate | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Esperanza | 0.129 | -0.643 | 0.912 | no |
| songrate | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Entre_Rios vs Esperanza | 0.218 | -0.559 | 1.020 | no |
| songrate | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Mar_Chiquita vs Esperanza | -0.087 | -0.649 | 0.457 | no |
| songrate | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Mar_Chiquita | 0.199 | -0.360 | 0.744 | no |
| songrate | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Entre_Rios vs Mar_Chiquita | 0.288 | -0.251 | 0.879 | no |
| songrate | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | E_Ibera vs Salta | -0.744 | -1.583 | 0.034 | no |
| songrate | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Entre_Rios vs Salta | -0.655 | -1.477 | 0.148 | no |
| songrate | Hypoxantha_Ruficollis | Mar_Chiquita vs Mar_Chiquita | Mar_Chiquita vs Salta | -0.960 | -1.576 | -0.284 | yes |
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.
- Element diversity (mst): now the strongest sympatry divergence among individual traits (β = 0.54, credible); no credible geographic-distance effect (β = -0.02). This reverses the earlier result from a stale cached fit, which had shown mst converging with both predictors.
- Element-level acoustic distance (elm_pca): strong sympatry divergence (β = 0.39, credible); geo distance not credible.
- Frequency range: strong sympatry divergence (β = 0.38, credible); also diverges with geo distance (β = -0.09, 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).
- Peak frequency, gap duration, song rate: no credible effect of either predictor.
2.2.8 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)
# variables missing from stable_df (e.g. never among the top loadings for
# that PC) come back as NA from the left join; treat those as "not stable"
# rather than letting NA propagate into the plot aesthetics below
pca_rot_stck$stable[is.na(pca_rot_stck$stable)] <- FALSE
# fully opaque bars for stable variables, more transparent for the rest
# (previously ifelse(stable, 1, 1) - always 1 regardless of stability)
pca_rot_stck$top_vars <- ifelse(pca_rot_stck$stable, 1, 0.4)
# Build colored labels per row
# escape characters that markdown/HTML treats specially (gridtext, used by
# element_markdown() below, parses these labels as markdown-with-HTML; an
# unescaped &, <, > or an odd number of _ / * in a variable name can break
# that parser and silently fail to render the whole plot)
escape_markdown <- function(x) {
x <- gsub("&", "&", x, fixed = TRUE)
x <- gsub("<", "<", x, fixed = TRUE)
x <- gsub(">", ">", x, fixed = TRUE)
x <- gsub("_", "_", x, fixed = TRUE)
x <- gsub("*", "*", x, fixed = TRUE)
x
}
# Colored labels: black for stable variables, gray for the rest
# (previously compared stable < 1.1, which is true whether stable is 0 or
# 1, so every label always took the black branch)
pca_rot_stck$label_col <- ifelse(
pca_rot_stck$stable,
paste0("<span style='color:black;'>", escape_markdown(pca_rot_stck$variable), "</span>"),
paste0("<span style='color:gray50;'>", escape_markdown(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
# alpha is used directly as a numeric value (scale_alpha_identity) rather
# than mapped through as.factor()/scale_alpha_manual(); the previous version
# passed the whole top_vars column as "values" to scale_alpha_manual, which
# expects one value per factor level, not per row - fragile even before the
# always-1 bug above was fixed
ggplot(pca_rot_stck, aes(x = var_facet, y = rotation, fill = Sign,
alpha = top_vars)) + geom_col() + coord_flip() + scale_alpha_identity(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 Collinearity among song-level features
Same check as for simple songs (see Collinearity among song-level features above): pairwise Pearson correlations among the seven song-level features, computed on the complex-song data set. Pairs with |r| ≥ 0.7 would be reduced to a single representative before the feature-based models are fitted.
Code
collin_complex <- check_feature_collinearity(
dat = complex_song_dat,
features = song_features,
labels = song_feature_labels
)
Pairs with |r| >= 0.7 :
feature_1 feature_2 r
Number of elements Song duration 0.964
Number of elements Element diversity (mst) 0.827
Song duration Element diversity (mst) 0.817
Frequency range Element diversity (mst) 0.795
Song rate Gap duration -0.718
Strongest pairwise correlations:
feature_1 feature_2 r
Number of elements Song duration 0.964
Number of elements Element diversity (mst) 0.827
Song duration Element diversity (mst) 0.817
Frequency range Element diversity (mst) 0.795
Song rate Gap duration -0.718
3.2.4 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)
# variables missing from stable_df (e.g. never among the top loadings for
# that PC) come back as NA from the left join; treat those as "not stable"
# rather than letting NA propagate into the plot aesthetics below
pca_rot_stck$stable[is.na(pca_rot_stck$stable)] <- FALSE
# fully opaque bars for stable variables, more transparent for the rest
# (previously ifelse(stable, 1, 1) - always 1 regardless of stability)
pca_rot_stck$top_vars <- ifelse(pca_rot_stck$stable, 1, 0.4)
# Build colored labels per row
# escape characters that markdown/HTML treats specially (gridtext, used by
# element_markdown() below, parses these labels as markdown-with-HTML; an
# unescaped &, <, > or an odd number of _ / * in a variable name can break
# that parser and silently fail to render the whole plot)
escape_markdown <- function(x) {
x <- gsub("&", "&", x, fixed = TRUE)
x <- gsub("<", "<", x, fixed = TRUE)
x <- gsub(">", ">", x, fixed = TRUE)
x <- gsub("_", "_", x, fixed = TRUE)
x <- gsub("*", "*", x, fixed = TRUE)
x
}
# Colored labels: black for stable variables, gray for the rest
# (previously compared stable < 1.1, which is true whether stable is 0 or
# 1, so every label always took the black branch)
pca_rot_stck$label_col <- ifelse(
pca_rot_stck$stable,
paste0("<span style='color:black;'>", escape_markdown(pca_rot_stck$variable), "</span>"),
paste0("<span style='color:gray50;'>", escape_markdown(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
# alpha is used directly as a numeric value (scale_alpha_identity) rather
# than mapped through as.factor()/scale_alpha_manual(); the previous version
# passed the whole top_vars column as "values" to scale_alpha_manual, which
# expects one value per factor level, not per row - fragile even before the
# always-1 bug above was fixed
ggplot(pca_rot_stck, aes(x = var_facet, y = rotation, fill = Sign,
alpha = top_vars)) + geom_col() + coord_flip() + scale_alpha_identity(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.5 Acoustic and geographic distances
3.2.5.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.5.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.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.
3.2.6.1 Model fitting
3.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.
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.6.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.6.1.1.2 Prior predictive check
Before looking at the posterior we check what the priors alone imply for the response. The model below is refitted with sample_prior = "only", so the likelihood is ignored and the draws come purely from the priors. For a standardized response (mean 0, SD 1) the prior predictive distribution should comfortably cover the observed range (roughly -3 to 3) without being absurdly wide; a prior that could only generate values within the observed range would be doing the data’s job, and one spanning hundreds of SDs would not be regularizing anything.
Code
prior_only_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(),
prior = priors,
sample_prior = "only",
chains = 2,
iter = 2000,
cores = 2,
seed = 123,
backend = "cmdstanr",
file = "./data/processed/fits/prior_only_acoustic_distance_complex"
)
pp_check(prior_only_complex, type = "dens_overlay", ndraws = 50) +
coord_cartesian(xlim = c(-20, 20)) +
labs(
title = "Prior predictive check (complex songs, multivariate model)",
subtitle = "y: observed standardized acoustic distance; y_rep: draws from the priors only"
)
3.2.6.1.1.3 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 + geo_distance_sc + (1 | species_pair) + (1 | mm(pop1, pop2)) + (1 | mm(individual1, individual2)) | 10000 | 4 | 1 | 5000 | 301 (0.015%) | 0 | 25506.31 | 13722.64 | 820698594 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sympatry1 | 0.046 | -0.040 | 0.131 | 1.001 | 25577.98 | 13796.50 |
| b_geo_distance_sc | 0.004 | -0.088 | 0.095 | 1.001 | 25506.31 | 13722.64 |
3.2.6.1.1.4 Posterior predictive check
Posterior predictive checks compare the observed data with data simulated from the fitted model (Gelman et al. 2014). Three views: the distribution of the response versus 50 replicated data sets (density overlay); the mean of the response within sympatric and allopatric pairs versus its posterior predictive distribution (the quantity the sympatry coefficient is about); and observed values against their posterior-averaged predictions. The observed data should fall within the range of the replicates.
Code
pp_check(sympatry_geo_model_complex, type = "dens_overlay", ndraws = 50) +
labs(title = "Posterior predictive check (complex songs): density overlay")
Code
pp_check(sympatry_geo_model_complex, type = "stat_grouped", stat = "mean", group = "sympatry", ndraws = 500) +
labs(title = "Posterior predictive check (complex songs): mean by sympatry (0 = allopatric, 1 = sympatric)")
Code
pp_check(sympatry_geo_model_complex, type = "scatter_avg", ndraws = 100) +
labs(title = "Posterior predictive check (complex songs): observed vs. posterior-averaged prediction")
Summary
- Overall acoustic distance (multivariate): neither predictor credible — sympatry β = 0.046 (CI [-0.040, 0.131]) and geographic distance β = 0.004 (CI [-0.088, 0.095]), both intervals cross zero. This is a change from the earlier stale fit: no credible evidence of either character displacement or isolation-by-distance at the overall multivariate level in complex songs.
- 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.6.1.1.5 Contrasts by species pair
The model above estimates a single sympatry effect for all species pairs; species pair enters only as a varying intercept, which shifts sympatric and allopatric comparisons by the same amount and therefore cannot tell whether the effect is shared by all pairs or driven by one of them. To obtain a separate sympatry effect for each species pair, the model is refitted with species-pair-specific sympatry slopes (a varying slope of sympatry by species pair, correlated with the varying intercept, with an LKJ(2) prior on the correlation). The per-pair effect is the fixed sympatry effect plus the pair’s deviation; the pairwise differences between species pairs are computed from the posterior draws. With only three species pairs the between-pair standard deviation of the sympatry slope, and its correlation with the intercept, are weakly informed by the data and lean on the priors; the per-pair effects are partially pooled towards the overall effect, so they should be read as shrunken estimates rather than as three independent fits. The fitting chunk is not evaluated when rendering (it takes as long as the main model); run it once and the results are loaded from the saved fit.
Code
priors_sp <- c(
priors,
prior(lkj(2), class = "cor")
)
sympatry_by_species_pair_complex <- brm(
acoustic_distance_sc ~ sympatry +
geo_distance_sc +
(1 + sympatry | species_pair) +
(1 | mm(pop1, pop2)) +
(1 | mm(individual1, individual2)),
data = sympatric_pairs_complex,
family = gaussian(),
prior = priors_sp,
cores = 4,
chains = 4,
iter = 10000,
backend = "cmdstanr",
threads = threading(8),
control = list(adapt_delta = 0.95, max_treedepth = 15),
file = "./data/processed/fits/acoustic_distance_sympatry_by_species_pair_complex_fit"
)Code
sp_fit_file <- "./data/processed/fits/acoustic_distance_sympatry_by_species_pair_complex_fit.rds"
if (file.exists(sp_fit_file)) {
sympatry_by_species_pair_complex <- readRDS(sp_fit_file)
sp_contrasts_complex <- species_pair_sympatry_effects(sympatry_by_species_pair_complex)
base::print(sp_contrasts_complex$plot)
} else {
message("Species-pair slope model not fitted yet - run the chunk above first.")
}Code
if (exists("sp_contrasts_complex")) print(sp_contrasts_complex$effects)Code
if (exists("sp_contrasts_complex")) print(sp_contrasts_complex$differences)Code
if (exists("sympatry_by_species_pair_complex")) {
extended_summary(
sympatry_by_species_pair_complex,
highlight = TRUE,
trace.palette = viridis::mako,
remove.intercepts = TRUE,
print.name = FALSE
)
}Population-pair contrasts from the species-pair slope model (same procedure as in the next section, but each species pair now has its own sympatry effect):
Code
if (exists("sympatry_by_species_pair_complex")) {
pop_contrasts_sp_complex <- population_pair_contrasts(
fit = sympatry_by_species_pair_complex,
dat = sympatric_pairs_complex,
geo = "mean"
)
base::print(pop_contrasts_sp_complex$plot)
}Code
if (exists("pop_contrasts_sp_complex")) print(pop_contrasts_sp_complex$contrasts)3.2.6.1.1.6 Contrasts among population pairs
Using the fitted model as is, the expected acoustic distance is computed for every observed combination of populations within each species pair, including the species-pair and population (multi-membership) effects and setting the individual effects to zero (i.e. for an average individual). Sympatric population pairs (the two species recorded at the same locality) are then contrasted against each allopatric population pair of the same species pair; each contrast is computed draw by draw, so its uncertainty interval reflects the joint posterior. Contrasts that share a population pair (e.g. the same sympatric pair against several allopatric ones) are not independent of each other and should be read as a set rather than as separate tests. Geographic distance is held at the pooled mean (geo_distance_sc = 0) so that the contrasts reflect sympatry and population effects rather than the distance between localities; set geo = "observed" to use each population pair’s own mean distance instead.
Code
pop_contrasts_complex <- population_pair_contrasts(
fit = sympatry_geo_model_complex,
dat = sympatric_pairs_complex,
geo = "mean"
)
base::print(pop_contrasts_complex$plot)
Code
print(pop_contrasts_complex$expected)| Species pair | Population pair | Sympatry | N comparisons | Geographic distance (sc) | Expected distance | l-95% UI | u-95% UI |
|---|---|---|---|---|---|---|---|
| Hypoxantha_Iberaensis | ER vs EI | Allopatric | 1462 | 0 | 0.141 | -0.285 | 0.596 |
| Hypoxantha_Iberaensis | MC vs EI | Allopatric | 816 | 0 | 0.072 | -0.482 | 0.649 |
| Hypoxantha_Palustris | ER vs EI | Allopatric | 7310 | 0 | 0.201 | -0.103 | 0.529 |
| Hypoxantha_Palustris | MC vs EI | Allopatric | 4080 | 0 | 0.132 | -0.354 | 0.630 |
| Hypoxantha_Ruficollis | EI vs ER | Allopatric | 935 | 0 | -0.525 | -0.947 | -0.099 |
| Hypoxantha_Ruficollis | MC vs ER | Allopatric | 408 | 0 | -0.378 | -0.913 | 0.140 |
| Hypoxantha_Ruficollis | EI vs Sal | Allopatric | 2145 | 0 | -0.121 | -0.467 | 0.244 |
| Hypoxantha_Ruficollis | ER vs Sal | Allopatric | 1677 | 0 | 0.095 | -0.293 | 0.491 |
| Hypoxantha_Ruficollis | MC vs Sal | Allopatric | 936 | 0 | 0.026 | -0.477 | 0.575 |
| Hypoxantha_Iberaensis | EI vs EI | Sympatric | 1870 | 0 | -0.028 | -0.443 | 0.391 |
| Hypoxantha_Palustris | EI vs EI | Sympatric | 9350 | 0 | 0.032 | -0.263 | 0.331 |
| Hypoxantha_Ruficollis | ER vs ER | Sympatric | 731 | 0 | -0.262 | -0.645 | 0.105 |
Code
print(pop_contrasts_complex$contrasts)| Species pair | Sympatric population pair | Allopatric population pair | Difference | l-95% UI | u-95% UI | UI excludes 0 |
|---|---|---|---|---|---|---|
| Hypoxantha_Iberaensis | EI vs EI | ER vs EI | -0.169 | -0.545 | 0.150 | no |
| Hypoxantha_Iberaensis | EI vs EI | MC vs EI | -0.100 | -0.616 | 0.348 | no |
| Hypoxantha_Palustris | EI vs EI | ER vs EI | -0.169 | -0.545 | 0.150 | no |
| Hypoxantha_Palustris | EI vs EI | MC vs EI | -0.100 | -0.616 | 0.348 | no |
| Hypoxantha_Ruficollis | ER vs ER | EI vs ER | 0.264 | -0.052 | 0.641 | no |
| Hypoxantha_Ruficollis | ER vs ER | MC vs ER | 0.117 | -0.352 | 0.587 | no |
| Hypoxantha_Ruficollis | ER vs ER | EI vs Sal | -0.141 | -0.638 | 0.316 | no |
| Hypoxantha_Ruficollis | ER vs ER | ER vs Sal | -0.357 | -0.773 | 0.049 | no |
| Hypoxantha_Ruficollis | ER vs ER | MC vs Sal | -0.288 | -0.921 | 0.280 | no |
3.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_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
# prior predictive check for the trait-level prior set (defined again here
# because the fitting chunk below is not evaluated when rendering), using
# one representative standardized response; the same priors and formula are
# used for every trait, so one check covers the set
priors_trait <- c(
prior(normal(0, 1), class = "Intercept"),
prior(normal(0, 0.5), class = "b"),
prior(exponential(2), class = "sd"),
prior(exponential(1), class = "sigma")
)
sympatric_pairs_complex$freqrange_distance_sc <- as.numeric(scale(sympatric_pairs_complex$freqrange_distance))
prior_only_trait_complex <- brm(
freqrange_distance_sc ~ sympatry +
geo_distance_sc +
(1 | species_pair) +
(1 | mm(pop1, pop2)) +
(1 | mm(individual1, individual2)),
data = sympatric_pairs_complex,
family = gaussian(),
prior = priors_trait,
sample_prior = "only",
chains = 2,
iter = 2000,
cores = 2,
seed = 123,
backend = "cmdstanr",
file = "./data/processed/fits/prior_only_trait_level_complex"
)
pp_check(prior_only_trait_complex, type = "dens_overlay", ndraws = 50) +
coord_cartesian(xlim = c(-10, 10)) +
labs(
title = "Prior predictive check (complex songs, trait-level priors)",
subtitle = "y: observed standardized frequency-range distance; y_rep: draws from the priors only"
)
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 = 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 table below gives, for every trait-level model above, the posterior mean, standard error and 95% uncertainty interval of each predictor; the last column flags intervals that exclude zero. Like the heatmap, it only includes models whose fit file exists.
Code
fixef_tab_complex <- fixef_table(model_files)
if (!is.null(fixef_tab_complex)) print(fixef_tab_complex)| Response | Predictor | Estimate | SE | l-95% UI | u-95% UI | UI excludes 0 |
|---|---|---|---|---|---|---|
| elm_pca | Geographic distance | -0.226 | 0.037 | -0.299 | -0.152 | yes |
| elm_pca | Sympatry | 1.083 | 0.035 | 1.014 | 1.153 | yes |
| freqrange | Geographic distance | -0.016 | 0.050 | -0.115 | 0.081 | no |
| freqrange | Sympatry | 0.024 | 0.047 | -0.069 | 0.116 | no |
| gapduration | Geographic distance | 0.051 | 0.047 | -0.040 | 0.143 | no |
| gapduration | Sympatry | 0.108 | 0.044 | 0.022 | 0.193 | yes |
| meanpeakf | Geographic distance | -0.096 | 0.047 | -0.188 | -0.003 | yes |
| meanpeakf | Sympatry | -0.146 | 0.045 | -0.233 | -0.059 | yes |
| mst | Geographic distance | 0.046 | 0.048 | -0.051 | 0.130 | no |
| mst | Sympatry | 0.082 | 0.044 | -0.007 | 0.162 | no |
| numelms | Geographic distance | -0.004 | 0.045 | -0.092 | 0.083 | no |
| numelms | Sympatry | -0.003 | 0.042 | -0.085 | 0.080 | no |
| songduration | Geographic distance | -0.011 | 0.046 | -0.102 | 0.081 | no |
| songduration | Sympatry | -0.004 | 0.043 | -0.090 | 0.081 | no |
| songrate | Geographic distance | -0.020 | 0.046 | -0.109 | 0.070 | no |
| songrate | Sympatry | -0.100 | 0.043 | -0.184 | -0.015 | yes |
The expandable section below provides the complete Bayesian model summaries underlying those estimates, including posterior parameter estimates, uncertainty 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")
}The expandable section below shows a posterior predictive check (density overlay, 50 replicated data sets) for every trait-level model that has finished fitting; models whose fit file is missing are skipped.
NotePosterior predictive checks
Code
for (i in model_files) {
fit <- tryCatch(readRDS(i), error = function(e) NULL)
if (is.null(fit)) next
model_name <- gsub("_sympatry_geographic_distance.*", "", tools::file_path_sans_ext(basename(i)))
cat("\n\n**", model_name, "**\n\n", sep = "")
# base::print - the document redefines print() as a kable wrapper
base::print(
pp_check(fit, type = "dens_overlay", ndraws = 50) +
labs(title = model_name)
)
cat("\n\n")
rm(fit)
invisible(gc(verbose = FALSE))
}elm_pca_distance_sc
freqrange_distance_sc
gapduration_distance_sc
meanpeakf_distance_sc
mst_distance_sc
numelms_distance_sc
songduration_distance_sc
songrate_distance_sc
3.2.6.1.2.1 Contrasts among population pairs, per feature
The same population-pair contrasts computed for the overall PCA-based model above, now repeated for every individual acoustic feature that has finished fitting: expected acoustic distance for each observed population pair (species-pair and population multi-membership effects included, individual effects set to zero, geographic distance held at the pooled mean), then sympatric vs allopatric population-pair differences within each species pair. Each feature gets its own plot (faceted by species pair, as in the main-model version above) rather than being crowded into one combined figure. As with the main-model version, contrasts sharing a population pair are not independent of each other.
Code
pop_contrasts_by_feature_complex <- population_pair_contrasts_by_feature(
model_files = model_files,
dat = sympatric_pairs_complex,
geo = "mean"
)
if (!is.null(pop_contrasts_by_feature_complex)) {
for (feat in names(pop_contrasts_by_feature_complex$plots)) {
cat("\n\n**", feat, "**\n\n", sep = "")
base::print(pop_contrasts_by_feature_complex$plots[[feat]])
cat("\n\n")
}
}elm pca
freqrange
gapduration
meanpeakf
mst
numelms
songduration
songrate
Code
if (!is.null(pop_contrasts_by_feature_complex)) print(pop_contrasts_by_feature_complex$contrasts)| Feature | Species pair | Sympatric population pair | Allopatric population pair | Difference | l-95% UI | u-95% UI | UI excludes 0 |
|---|---|---|---|---|---|---|---|
| elm pca | Hypoxantha_Iberaensis | EI vs EI | ER vs EI | 1.045 | 0.820 | 1.252 | yes |
| elm pca | Hypoxantha_Iberaensis | EI vs EI | MC vs EI | 1.123 | 0.858 | 1.461 | yes |
| elm pca | Hypoxantha_Palustris | EI vs EI | ER vs EI | 1.045 | 0.820 | 1.252 | yes |
| elm pca | Hypoxantha_Palustris | EI vs EI | MC vs EI | 1.123 | 0.858 | 1.461 | yes |
| elm pca | Hypoxantha_Ruficollis | ER vs ER | EI vs ER | 1.122 | 0.915 | 1.347 | yes |
| elm pca | Hypoxantha_Ruficollis | ER vs ER | MC vs ER | 1.161 | 0.900 | 1.518 | yes |
| elm pca | Hypoxantha_Ruficollis | ER vs ER | EI vs Sal | 0.985 | 0.617 | 1.242 | yes |
| elm pca | Hypoxantha_Ruficollis | ER vs ER | ER vs Sal | 0.947 | 0.640 | 1.158 | yes |
| elm pca | Hypoxantha_Ruficollis | ER vs ER | MC vs Sal | 1.025 | 0.650 | 1.362 | yes |
| freqrange | Hypoxantha_Iberaensis | EI vs EI | ER vs EI | -0.407 | -0.784 | -0.016 | yes |
| freqrange | Hypoxantha_Iberaensis | EI vs EI | MC vs EI | -0.173 | -0.633 | 0.251 | no |
| freqrange | Hypoxantha_Palustris | EI vs EI | ER vs EI | -0.407 | -0.784 | -0.016 | yes |
| freqrange | Hypoxantha_Palustris | EI vs EI | MC vs EI | -0.173 | -0.633 | 0.251 | no |
| freqrange | Hypoxantha_Ruficollis | ER vs ER | EI vs ER | 0.459 | 0.047 | 0.852 | yes |
| freqrange | Hypoxantha_Ruficollis | ER vs ER | MC vs ER | 0.260 | -0.180 | 0.789 | no |
| freqrange | Hypoxantha_Ruficollis | ER vs ER | EI vs Sal | 0.169 | -0.275 | 0.651 | no |
| freqrange | Hypoxantha_Ruficollis | ER vs ER | ER vs Sal | -0.264 | -0.643 | 0.077 | no |
| freqrange | Hypoxantha_Ruficollis | ER vs ER | MC vs Sal | -0.030 | -0.598 | 0.561 | no |
| gapduration | Hypoxantha_Iberaensis | EI vs EI | ER vs EI | 0.014 | -0.283 | 0.220 | no |
| gapduration | Hypoxantha_Iberaensis | EI vs EI | MC vs EI | 0.054 | -0.280 | 0.301 | no |
| gapduration | Hypoxantha_Palustris | EI vs EI | ER vs EI | 0.014 | -0.283 | 0.220 | no |
| gapduration | Hypoxantha_Palustris | EI vs EI | MC vs EI | 0.054 | -0.280 | 0.301 | no |
| gapduration | Hypoxantha_Ruficollis | ER vs ER | EI vs ER | 0.204 | 0.002 | 0.515 | yes |
| gapduration | Hypoxantha_Ruficollis | ER vs ER | MC vs ER | 0.149 | -0.139 | 0.489 | no |
| gapduration | Hypoxantha_Ruficollis | ER vs ER | EI vs Sal | 0.118 | -0.218 | 0.453 | no |
| gapduration | Hypoxantha_Ruficollis | ER vs ER | ER vs Sal | 0.023 | -0.316 | 0.238 | no |
| gapduration | Hypoxantha_Ruficollis | ER vs ER | MC vs Sal | 0.063 | -0.375 | 0.420 | no |
| meanpeakf | Hypoxantha_Iberaensis | EI vs EI | ER vs EI | -0.151 | -0.372 | 0.076 | no |
| meanpeakf | Hypoxantha_Iberaensis | EI vs EI | MC vs EI | -0.130 | -0.366 | 0.172 | no |
| meanpeakf | Hypoxantha_Palustris | EI vs EI | ER vs EI | -0.151 | -0.372 | 0.076 | no |
| meanpeakf | Hypoxantha_Palustris | EI vs EI | MC vs EI | -0.130 | -0.366 | 0.172 | no |
| meanpeakf | Hypoxantha_Ruficollis | ER vs ER | EI vs ER | -0.142 | -0.362 | 0.073 | no |
| meanpeakf | Hypoxantha_Ruficollis | ER vs ER | MC vs ER | -0.126 | -0.368 | 0.181 | no |
| meanpeakf | Hypoxantha_Ruficollis | ER vs ER | EI vs Sal | -0.161 | -0.496 | 0.120 | no |
| meanpeakf | Hypoxantha_Ruficollis | ER vs ER | ER vs Sal | -0.165 | -0.414 | 0.042 | no |
| meanpeakf | Hypoxantha_Ruficollis | ER vs ER | MC vs Sal | -0.145 | -0.493 | 0.180 | no |
| mst | Hypoxantha_Iberaensis | EI vs EI | ER vs EI | -0.082 | -0.403 | 0.269 | no |
| mst | Hypoxantha_Iberaensis | EI vs EI | MC vs EI | -0.066 | -0.555 | 0.374 | no |
| mst | Hypoxantha_Palustris | EI vs EI | ER vs EI | -0.082 | -0.403 | 0.269 | no |
| mst | Hypoxantha_Palustris | EI vs EI | MC vs EI | -0.066 | -0.555 | 0.374 | no |
| mst | Hypoxantha_Ruficollis | ER vs ER | EI vs ER | 0.243 | -0.122 | 0.574 | no |
| mst | Hypoxantha_Ruficollis | ER vs ER | MC vs ER | 0.097 | -0.406 | 0.533 | no |
| mst | Hypoxantha_Ruficollis | ER vs ER | EI vs Sal | -0.201 | -0.758 | 0.217 | no |
| mst | Hypoxantha_Ruficollis | ER vs ER | ER vs Sal | -0.363 | -0.797 | 0.077 | no |
| mst | Hypoxantha_Ruficollis | ER vs ER | MC vs Sal | -0.347 | -1.063 | 0.179 | no |
| numelms | Hypoxantha_Iberaensis | EI vs EI | ER vs EI | -0.018 | -0.302 | 0.261 | no |
| numelms | Hypoxantha_Iberaensis | EI vs EI | MC vs EI | -0.009 | -0.338 | 0.341 | no |
| numelms | Hypoxantha_Palustris | EI vs EI | ER vs EI | -0.018 | -0.302 | 0.261 | no |
| numelms | Hypoxantha_Palustris | EI vs EI | MC vs EI | -0.009 | -0.338 | 0.341 | no |
| numelms | Hypoxantha_Ruficollis | ER vs ER | EI vs ER | 0.013 | -0.253 | 0.290 | no |
| numelms | Hypoxantha_Ruficollis | ER vs ER | MC vs ER | 0.007 | -0.312 | 0.381 | no |
| numelms | Hypoxantha_Ruficollis | ER vs ER | EI vs Sal | -0.145 | -0.657 | 0.176 | no |
| numelms | Hypoxantha_Ruficollis | ER vs ER | ER vs Sal | -0.160 | -0.590 | 0.104 | no |
| numelms | Hypoxantha_Ruficollis | ER vs ER | MC vs Sal | -0.151 | -0.712 | 0.213 | no |
| songduration | Hypoxantha_Iberaensis | EI vs EI | ER vs EI | -0.054 | -0.370 | 0.201 | no |
| songduration | Hypoxantha_Iberaensis | EI vs EI | MC vs EI | -0.021 | -0.409 | 0.354 | no |
| songduration | Hypoxantha_Palustris | EI vs EI | ER vs EI | -0.054 | -0.370 | 0.201 | no |
| songduration | Hypoxantha_Palustris | EI vs EI | MC vs EI | -0.021 | -0.409 | 0.354 | no |
| songduration | Hypoxantha_Ruficollis | ER vs ER | EI vs ER | 0.048 | -0.233 | 0.358 | no |
| songduration | Hypoxantha_Ruficollis | ER vs ER | MC vs ER | 0.030 | -0.324 | 0.442 | no |
| songduration | Hypoxantha_Ruficollis | ER vs ER | EI vs Sal | -0.157 | -0.648 | 0.223 | no |
| songduration | Hypoxantha_Ruficollis | ER vs ER | ER vs Sal | -0.209 | -0.662 | 0.084 | no |
| songduration | Hypoxantha_Ruficollis | ER vs ER | MC vs Sal | -0.175 | -0.775 | 0.262 | no |
| songrate | Hypoxantha_Iberaensis | EI vs EI | ER vs EI | -0.108 | -0.359 | 0.127 | no |
| songrate | Hypoxantha_Iberaensis | EI vs EI | MC vs EI | -0.148 | -0.544 | 0.099 | no |
| songrate | Hypoxantha_Palustris | EI vs EI | ER vs EI | -0.108 | -0.359 | 0.127 | no |
| songrate | Hypoxantha_Palustris | EI vs EI | MC vs EI | -0.148 | -0.544 | 0.099 | no |
| songrate | Hypoxantha_Ruficollis | ER vs ER | EI vs ER | -0.093 | -0.326 | 0.142 | no |
| songrate | Hypoxantha_Ruficollis | ER vs ER | MC vs ER | -0.141 | -0.508 | 0.113 | no |
| songrate | Hypoxantha_Ruficollis | ER vs ER | EI vs Sal | -0.124 | -0.465 | 0.188 | no |
| songrate | Hypoxantha_Ruficollis | ER vs ER | ER vs Sal | -0.132 | -0.395 | 0.083 | no |
| songrate | Hypoxantha_Ruficollis | ER vs ER | MC vs Sal | -0.172 | -0.649 | 0.144 | no |
Summary
At the aggregate (multivariate) level, neither predictor showed a credible effect: sympatry β = 0.046 (CI [-0.040, 0.131]) and geographic distance β = 0.004 (CI [-0.088, 0.095]), both crossing zero. Any signal in complex songs is trait-specific rather than reflecting a broad shift in overall acoustic distance.
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: no credible effect at the aggregate multivariate level, 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 — and the overall multivariate signal — appears conserved, convergent, or simply absent 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)
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-25
pandoc 3.8.3 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/x86_64/ (via rmarkdown)
quarto 1.9.38 @ /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)
ggh4x * 0.3.1 2025-05-30 [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.
──────────────────────────────────────────────────────────────────────────────