Statistical analysis

Vocal character displacement in southern capuchinos

Author

Published

July 15, 2026

Code

# options to customize chunk outputs
knitr::opts_chunk$set(
  message = FALSE
)

Purpose

Look at 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",
  github = "maRce10/PhenotypeSpace",
  "kableExtra",
  "ggplot2",
  "ggtext",
  github = "stan-dev/cmdstanr"
)

# install/load packages
sketchy::load_packages(packages = packages)

options(knitr.kable.NA = '', brms.file_refit = "never")



print <- function(x, row.names = FALSE) {
  kb <- kable(x, row.names = row.names, digits = 4, "html")
  kb <- kable_styling(kb,
                      bootstrap_options = c("striped", "hover", "condensed", "responsive"))
  scroll_box(kb, width = "100%")
}



# set theme globally
theme_set(theme_classic(base_size = 20))

get_stable_loadings <- function(pca, pcs = 4, B = 1000, cum_threshold = 0.5,
    freq_threshold = 0.75, seed = 123) {

    set.seed(seed)

    # Reconstruct centered data
    X <- pca$x %*% t(pca$rotation)

    p <- ncol(pca$rotation)
    orig_rot <- pca$rotation[, 1:pcs]

    boot_loadings <- array(NA, dim = c(p, pcs, B))

    # ----------------------------- Bootstrap PCA
    # -----------------------------

    for (b in 1:B) {

        idx <- sample(1:nrow(X), replace = TRUE)
        Xb <- X[idx, ]

        pca_b <- prcomp(Xb, scale. = TRUE)
        rot_b <- pca_b$rotation[, 1:pcs]

        # Align signs
        for (k in 1:pcs) {
            if (cor(rot_b[, k], orig_rot[, k]) < 0) {
                rot_b[, k] <- -rot_b[, k]
            }
        }

        boot_loadings[, , b] <- rot_b
    }

    # ----------------------------- Summary statistics
    # -----------------------------

    abs_boot <- abs(boot_loadings)

    mean_loading <- apply(abs_boot, c(1, 2), mean)
    ci_lower <- apply(abs_boot, c(1, 2), quantile, 0.025)
    ci_upper <- apply(abs_boot, c(1, 2), quantile, 0.975)

    # ----------------------------- Stability frequency
    # -----------------------------

    top_freq <- matrix(0, nrow = p, ncol = pcs)

    for (b in 1:B) {
        for (k in 1:pcs) {

            sq <- boot_loadings[, k, b]^2
            ord <- order(sq, decreasing = TRUE)
            cumprop <- cumsum(sq[ord])/sum(sq)

            selected <- ord[cumprop <= cum_threshold]
            selected <- c(selected, ord[min(which(cumprop >= cum_threshold))])

            top_freq[selected, k] <- top_freq[selected, k] + 1
        }
    }

    top_freq <- top_freq/B

    stable <- (top_freq >= freq_threshold) & (ci_lower > 0 | ci_upper <
        0)

    # ----------------------------- Return tidy dataframe
    # -----------------------------

    data.frame(variable = rep(rownames(pca$rotation), pcs), ind = rep(colnames(pca$rotation)[1:pcs],
        each = p), mean_loading = as.vector(mean_loading), ci_lower = as.vector(ci_lower),
        ci_upper = as.vector(ci_upper), freq = as.vector(top_freq),
        stable = as.vector(stable))
}

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

  effects_df <- data.frame()

  for(i in seq_along(model_files)) {

    fit <- readRDS(model_files[i])

    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
    )
  }

  rownames(effects_df) <- NULL


  # significance

  effects_df$sig <- with(
    effects_df,
    Q2.5 * Q97.5 > 0
  )

  # clean names

  effects_df$predictor <- gsub("^scale\\(", "", effects_df$predictor)
  effects_df$predictor <- gsub("\\)$", "", effects_df$predictor)

  effects_df$predictor <- gsub("^mo", "", effects_df$predictor)
  effects_df$predictor <- gsub("^mi", "", effects_df$predictor)
  effects_df$predictor <- gsub("_sc$", "", effects_df$predictor)

  effects_df$response <- gsub("^mi", "", effects_df$response)
  effects_df$response <- gsub("_sc$", "", effects_df$response)
# effects_df$predictor <- ifelse(grepl("sympatry", effects_df$predictor), "sympatry", effects_df$predictor)

  
  # average duplicated cells if present

  effects_df <- aggregate(
    cbind(
      Estimate,
      sig
    ) ~ predictor + response,
    data = effects_df,
    FUN = mean
  )

  effects_df$sig <- effects_df$sig > 0.5
  
  # complete combinations

  all_combos <- expand.grid(
    predictor = unique(effects_df$predictor),
    response = unique(effects_df$response)
  )

  plot_df <- merge(
    all_combos,
    effects_df,
    by = c("predictor", "response"),
    all.x = TRUE
  )

  plot_df$response <- gsub("_distance$", "", plot_df$response)

plot_df$predictor <- gsub("scalegeo_", "Geographic\n", plot_df$predictor)

plot_df$predictor <- gsub("sympatry1", "Sympatry", plot_df$predictor)
    
  lim <- max(abs(plot_df$Estimate), na.rm = TRUE)

  
  ggplot(
    plot_df,
    aes(
      predictor,
      response
    )
  ) +

    geom_tile(
      fill = "grey90",
      colour = "white"
    ) +

    geom_tile(
      data = plot_df[!is.na(plot_df$Estimate), ],
      aes(fill = Estimate),
      colour = "white"
    ) +

    geom_text(
      data = plot_df[!is.na(plot_df$Estimate), ],
      aes(
        label = sprintf("%.2f", Estimate),
        colour = sig
      ),
      fontface = "bold",
      size = 3
    ) +
      
    # scale_fill_gradient2(
    #   low = rep("#403B78", 2),
    #   mid = "white",
    #   high = rep("#DEF5E5", rep = 2),
    #   midpoint = 0,
    #   name = "Estimate"
    # ) +
    # 
      scale_fill_gradient2(
  low = "#403B78",
  mid = "white",
  high = "#A0DFB9CC",
  midpoint = 0,
  limits = c(-lim, lim),
  oob = scales::squish,
  name = "Estimate"
) +

    scale_color_manual(
      values = c(
        "TRUE" = "black",
        "FALSE" = "grey70"
      ),
      guide = "none"
    ) +

    labs(
      x = "Predictor",
      y = "Response"
    ) +

    theme_classic() +

    theme(
      axis.text.x = element_text(
        angle = 45,
        hjust = 1
      )
    )
}

1 Simple songs

1.1 Song level

1.1.1 Read and prepare data

Code

simple_songs <- read.csv("./data/raw/Sporophila song data_song-level plus MCP MST and PCA_simple songs.csv")
individual_coord <- read.csv("./data/raw/Coordenadas_individuos_simple_songs.csv")

# assign coordinates to each individual
simple_songs_lat_long <- simple_songs |>
  left_join(individual_coord, by = "Individual")

# names(simple_songs_lat_long)

1.1.2 Select variables and filter populations

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

1.1.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 4 principal components (PCs) for subsequent analyses, which together explained rround(cumsum(summary(pca)\(importance[2, ])[(which(cumsum(summary(pca)\)importance[2, ]) > 0.9))[1]] * 100, 1) `% of the variance in the data.

Code

pca_rot_stck <- stack(pca_rot)

pca_rot_stck$variable <- rownames(pca_rot)
pca_rot_stck$values[pca_rot_stck$ind == "PC1"] <- pca_rot_stck$values[pca_rot_stck$ind ==
    "PC1"]
pca_rot_stck$Sign <- ifelse(pca_rot_stck$values > 0, "Positive", "Negative")
pca_rot_stck$rotation <- abs(pca_rot_stck$values)
pca_rot_stck$ind_var <- paste0(pca_rot_stck$ind, " (", sapply(pca_rot_stck$ind,
    function(x) pca_var[names(pca_var) == x]), "%)")


pca_rot_stck$top_vars <- ave(abs(pca_rot_stck$values), pca_rot_stck$ind,
    FUN = function(x) {

        # Order decreasing
        ord <- order(x, decreasing = TRUE)
        x_sorted <- x[ord]

        # Cumulative proportion
        cumprop <- cumsum(x_sorted)/sum(x_sorted)

        selected_sorted <- cumprop <= 0.5  # variables with 50% of contribution 
        selected_sorted[which(cumprop >= 0.5)[1]] <- TRUE

        # Return logical vector in original order
        selected <- logical(length(x))
        selected[ord] <- selected_sorted

        selected
    })


# Create facet-specific variable
pca_rot_stck$var_facet <- paste(pca_rot_stck$ind, pca_rot_stck$variable,
    sep = "_")

# Reorder within each ind by rotation (largest at top after
# coord_flip)
pca_rot_stck <- do.call(rbind, lapply(split(pca_rot_stck, pca_rot_stck$ind),
    function(df) {

        df$var_facet <- factor(df$var_facet, levels = df$var_facet[order(df$rotation)])

        df
    }))


# add which variables are stable
stable_df <- get_stable_loadings(pca, pcs = 5, B = 1000, cum_threshold = 0.5,
    freq_threshold = 0.5, seed = 123)

pca_rot_stck <- merge(pca_rot_stck, stable_df, by = c("variable",
    "ind"), all.x = TRUE)

pca_rot_stck$top_vars <- ifelse(pca_rot_stck$stable, 1, 1)


# Build colored labels per row

# Colored labels
pca_rot_stck$label_col <- ifelse(pca_rot_stck$stable < 1.1, paste0("<span style='color:black;'>",
    pca_rot_stck$variable, "</span>"), paste0("<span style='color:gray50;'>",
    pca_rot_stck$variable, "</span>"))

# Named vector for labels
label_vec <- setNames(pca_rot_stck$label_col, pca_rot_stck$var_facet)

# absolute CI
pca_rot_stck$ci_low_plot <- abs(pca_rot_stck$ci_lower)
pca_rot_stck$ci_high_plot <- abs(pca_rot_stck$ci_upper)


# Plot
ggplot(pca_rot_stck, aes(x = var_facet, y = rotation, fill = Sign,
    alpha = as.factor(top_vars))) + geom_col() + coord_flip() + scale_alpha_manual(values = pca_rot_stck$top_vars,
    guide = NULL) + scale_x_discrete(labels = label_vec, name = "Variable") +
    labs(x = "Rotation") + scale_fill_viridis_d(alpha = 0.7, begin = 0.2,
    end = 0.8) + facet_wrap(~ind_var, scales = "free_y", nrow = 3) + theme_classic() +
    theme(axis.text.y = element_markdown())

Code

# bind PCA scores with metadata
pca_scores <- cbind(
  pca$x,
  simple_song_dat[, c(
    "Population",
    "species",
    "location",
    "Individual",
    "Lat",
    "Lon",
    "meanpeakf",
    "num.elms",
    "song.duration",
    "song.rate",
    "gap.duration",
    "freq.range.Min5toMax95",
    "mst"
  )]
)

summ_pca <- summary(pca)

# select the PCs explaining at least 90%
pcs <- grep("^PC", names(pca_scores), value = TRUE)[1:min(which(summ_pca$importance[3,] > 0.9))]

pca_scores <- pca_scores |>
  mutate(
    across(all_of(pcs), ~ as.numeric(.x)),
    Lat = as.numeric(Lat),
    Lon = as.numeric(Lon),
    species = as.character(species),
    location = tryCatch(
      iconv(location, from = "", to = "UTF-8"),
      error = function(e) location
    )
  )

1.1.4 Acoustic and geographic distances

1.1.4.1 Build pairwise distance dataset

Code

df_clean_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)
  ) |>
  filter(complete.cases(across(all_of(c(pcs, "Lat", "Lon", "species", "Individual")))))

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

1.1.4.2 Convert to long table

Code

# acoustic distance between songs (Euclidean multivariate, unscaled PCs)
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"
))

# 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_id
rownames(dist_geo_mat) <- colnames(dist_geo_mat) <- df_clean_simple$song_id

# use only upper triangle
idx <- which(upper.tri(dist_acoustic_mat), arr.ind = TRUE)

dist_acoustic_long <- data.frame(
  id1      = rownames(dist_acoustic_mat)[idx[, 1]],
  id2      = colnames(dist_acoustic_mat)[idx[, 2]],
  acoustic_distance = dist_acoustic_mat[idx],
  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]
)

# parse species, population, and individual from song IDs
parts1 <- strsplit(as.character(dist_acoustic_long$id1), "-")
parts2 <- strsplit(as.character(dist_acoustic_long$id2), "-")

dist_acoustic_long$species1    <- sapply(parts1, `[`, 1)
dist_acoustic_long$population1 <- sapply(parts1, `[`, 2)
dist_acoustic_long$individual1 <- sapply(parts1, `[`, 3)

dist_acoustic_long$species2    <- sapply(parts2, `[`, 1)
dist_acoustic_long$population2 <- sapply(parts2, `[`, 2)
dist_acoustic_long$individual2 <- sapply(parts2, `[`, 3)

# keep only between-species comparisons
dist_acoustic_long <- dist_acoustic_long[
  dist_acoustic_long$species1 != dist_acoustic_long$species2, ]

# sympatry flag: 1 if same population, 0 otherwise
dist_acoustic_long$sympatry <- as.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 = "_")
)

dist_acoustic_long$individual1  <- factor(dist_acoustic_long$individual1)
dist_acoustic_long$individual2  <- factor(dist_acoustic_long$individual2)
dist_acoustic_long$population1  <- factor(dist_acoustic_long$population1)
dist_acoustic_long$population2  <- factor(dist_acoustic_long$population2)
dist_acoustic_long$species_pair <- factor(dist_acoustic_long$species_pair)

1.1.5 Statistical analysis

To evaluate whether acoustic divergence between heterospecific songs differed between sympatric and allopatric population comparisons, we fitted a series of Bayesian mixed-effects models while accounting for the non-independence inherent to pairwise distance data.

The global model was specified as: \[ \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.
species pair is a random intercept accounting for baseline differences in acoustic divergence among heterospecific species combinations.
mm(population(_1), population(_2)) is a multi-membership random effect accounting for repeated use of the same populations across pairwise comparisons.
mm(individual(_1), individual(_2)) is a multi-membership random effect accounting for repeated use of the same individuals across pairwise comparisons.

Model specifications:

Models were fitted in a Bayesian framework using the brms package.
The response variable was pairwise Euclidean acoustic distance and was modeled using a Gaussian error distribution.
Only heterospecific comparisons were included.
Analyses were restricted to species pairs for which both sympatric and allopatric population comparisons were available. This restriction ensured that sympatry effects were estimated within the same species-pair contrasts rather than being confounded by species pairs occurring exclusively in sympatry or exclusively in allopatry.
Species-pair identity was included as a random intercept to account for inherent differences in acoustic divergence among species combinations.
Population identity was modeled using a multi-membership random effect because each pairwise comparison simultaneously involves two populations, and each population contributes to multiple pairwise distances.
Individual identity was modeled using a multi-membership random effect because each pairwise comparison simultaneously involves two individuals, and each individual contributes to multiple pairwise distances.
Only unique pairwise comparisons were retained and self-comparisons were excluded.
Weakly informative priors were specified for all parameters.
Models were fitted using Hamiltonian Monte Carlo as implemented in Stan through the cmdstanr backend.

To assess the relative importance of sympatry and geographic distance, three competing models were fitted: (i) a sympatry-only model, (ii) a geographic-distance-only model, and (iii) a model including both predictors. Model performance was compared using approximate leave-one-out cross-validation (LOO). The coefficient for sympatry in the global model therefore represents differences in acoustic divergence between sympatric and allopatric population comparisons after accounting for geographic distance and the hierarchical structure of the data.

The model was fitted with two datasets:

Using only the species pairs that had both sympatric and allopatric comparisons. This directly evaluates if heterospecific songs more divergent in sympatry than in allopatry.
Using the entire dataset, including species pairs that only had sympatric or allopatric comparisons. This is a complementary test of whether sympatric species are generally more divergent than allopatric species.

1.1.5.1 Model fitting

1.1.5.1.1 PCA-based acoustic distance

Prepare data for modeling by restricting to species pairs with both sympatric and allopatric comparisons and creating appropriate random effect structures.

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	EI	ER	MC	Sal
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	EI	ER, MC
Hypoxantha–Palustris	EI	ER, MC
Hypoxantha–Ruficollis	ER, MC	EI, E, Sal
Iberaensis–Palustris	EI	none
Iberaensis–Ruficollis	none	EI, E, ER, MC, Sal
Palustris–Ruficollis	none	EI, E, ER, MC, Sal

Only three species pairs were kept: Hypoxantha_Iberaensis, Hypoxantha_Palustris, Hypoxantha_Ruficollis

1.1.5.1.1.1 Species pairs with both sympatric and allopatric populations

Code

sympatric_pairs_simple <- dist_acoustic_long[
    dist_acoustic_long$species_pair %in%
keep_pairs,
  ]

# create species-specific population IDs
sympatric_pairs_simple$pop1 <- interaction(
  sympatric_pairs_simple$species1,
  sympatric_pairs_simple$population1,
  drop = TRUE
)

sympatric_pairs_simple$pop2 <- interaction(
  sympatric_pairs_simple$species2,
  sympatric_pairs_simple$population2,
  drop = TRUE
)

# make species_pair a factor (fixed effect)
sympatric_pairs_simple$species_pair <- factor(sympatric_pairs_simple$species_pair)

# scale acoustic_distance
sympatric_pairs_simple$acoustic_distance_sc <- scale(sympatric_pairs_simple$acoustic_distance)

Code

# weak priors
priors <- c(
  
  # Intercept
  prior(normal(0, 5), class = "Intercept"),
  
  # Fixed effect of sympatry
  prior(normal(0, 2), class = "b"),
  
  # Random-effect SDs
  prior(exponential(1), class = "sd"),
  
  # Residual SD
  prior(exponential(1), class = "sigma")
  
)

#fit model
sympatry_geo_model <- brm(
  acoustic_distance_sc ~ sympatry + 
    scale(geo_distance) +
    (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_fit"
)


sympatry_geo_model <- add_criterion(
    sympatry_geo_model,
    criterion = "loo"
  )


sympatry_model <- brm(
  acoustic_distance_sc ~ sympatry + 
    (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_fit"
)


sympatry_model <- add_criterion(
    sympatry_model,
    criterion = "loo"
  )

geo_model <- brm(
  acoustic_distance_sc ~ scale(geo_distance) +
    (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_geographic_distance_fit"
  )


geo_model <- add_criterion(
    geo_model,
    criterion = "loo"
  )

# beepr::beep(3)

Model comparison

Code

comparison <- as.data.frame(loo_compare(
  sympatry_geo_model,
  sympatry_model,
  geo_model
))

Warning: Not all models have the same y variable. ('yhash' attributes do not
match)

Code

comparison$model <- c(
    "sympatry + geographic distance",
    "sympatry only",
    "geographic distance only"
  )
  
print(comparison[, c("model", "elpd_diff", "se_diff")])

model	elpd_diff	se_diff
sympatry + geographic distance	0.00	0.0000
sympatry only	-19338.82	50.5643
geographic distance only	-19354.74	50.6316

1.1.5.1.1.2 Fit summary

Code

models <- list.files("./data/processed/", pattern = "fit", full.names = TRUE)

models <- grep("acoustic", models, value = TRUE)
models <- grep("all_data", models, invert = TRUE,  value = TRUE)


for (i in models){

extended_summary(
  read.file = i,
  highlight = TRUE,
  trace.palette = viridis::mako,
  remove.intercepts = TRUE,
  print.name = FALSE
)

    }

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-normal(0, 5) sd-exponential(1) sigma-exponential(1)	acoustic_distance_sc ~ sympatry + scale(geo_distance) + (1 \| species_pair) + (1 \| mm(pop1, pop2)) + (1 \| mm(individual1, individual2))	10000	4	1	5000	752 (0.038%)	0	24515.86	13314.61	1871852522

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_sympatry1	0.155	0.096	0.215	1	24623.43	13339.44
b_scalegeo_distance	-0.001	-0.031	0.029	1	24515.86	13314.61

Code

###### Entire dataset

sympatric_pairs_simple_all <- dist_acoustic_long

# create species-specific population IDs
sympatric_pairs_simple_all$pop1 <- interaction(
  sympatric_pairs_simple_all$species1,
  sympatric_pairs_simple_all$population1,
  drop = TRUE
)

sympatric_pairs_simple_all$pop2 <- interaction(
  sympatric_pairs_simple_all$species2,
  sympatric_pairs_simple_all$population2,
  drop = TRUE
)

# make species_pair a factor (fixed effect)
sympatric_pairs_simple_all$species_pair <- factor(sympatric_pairs_simple_all$species_pair)

# scale acoustic_distance
sympatric_pairs_simple_all$acoustic_distance_sc <- scale(sympatric_pairs_simple_all$acoustic_distance)

Code

# weak priors
priors <- c(
  
  # Intercept
  prior(normal(0, 5), class = "Intercept"),
  
  # Fixed effect of sympatry
  prior(normal(0, 2), class = "b"),
  
  # Random-effect SDs
  prior(exponential(1), class = "sd"),
  
  # Residual SD
  prior(exponential(1), class = "sigma")
  
)

#fit model
sympatry_geo_model <- brm(
  acoustic_distance_sc ~ sympatry + 
    scale(geo_distance) +
    (1 | species_pair) +
    (1 | mm(pop1, pop2)) +
    (1 | mm(individual1, individual2)),
  data = sympatric_pairs_simple_all,
  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_all_data_fit"
)


sympatry_geo_model <- add_criterion(
    sympatry_geo_model,
    criterion = "loo"
  )


sympatry_model <- brm(
  acoustic_distance_sc ~ sympatry + 
    (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_all_data_fit"
)


sympatry_model <- add_criterion(
    sympatry_model,
    criterion = "loo"
  )

geo_model <- brm(
  acoustic_distance_sc ~ scale(geo_distance) +
    (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_geographic_distance_all_data_fit"
  )


geo_model <- add_criterion(
    geo_model,
    criterion = "loo"
  )

# beepr::beep(3)

Code

####### Model comparison

comparison <- as.data.frame(loo_compare(
  sympatry_geo_model,
  sympatry_model,
  geo_model
))

comparison$model <- c(
    "sympatry + geographic distance",
    "sympatry only",
    "geographic distance only"
  )
  
print(comparison[, c("model", "elpd_diff", "se_diff")])

Code

models <- list.files("./data/processed/", pattern = "fit", full.names = TRUE)

models <- grep("all_data_", models, value = TRUE)

for (i in models){

extended_summary(
  read.file = i,
  highlight = TRUE,
  trace.palette = viridis::mako,
  remove.intercepts = TRUE,
  print.name = FALSE
)

    }

Summary

The model including both sympatry and geographic distance had the highest predictive performance and was therefore the best-supported model according to leave-one-out cross-validation (LOO).
The sympatry-only model performed nearly as well as the full model (ΔELPD = -2.7 ± 2.9), indicating that geographic distance contributed little additional predictive information once sympatry was accounted for.
The geographic-distance-only model performed substantially worse than the full model (ΔELPD = -18.7 ± 6.2), demonstrating that geographic distance alone was a poor predictor of acoustic divergence relative to models that included sympatry.
The large performance difference between the geographic-distance-only model and the models containing sympatry suggests that sympatry is a much stronger predictor of acoustic divergence than geographic distance.
Overall, these results indicate that variation in acoustic divergence among heterospecific populations is primarily associated with whether populations occur in sympatry rather than with the geographic distance separating them.
Estimates show that sympatric population pairs have, on average, acoustic distances greater than those of allopatric pairs

1.1.5.1.2 Feature-based acoustic distance

The following plot summarizes the relationship between sympatry and acoustic divergence for each of the acoustic features used to calculate acoustic distance. The model was fit separately for each feature, and the estimated effect of sympatry on acoustic divergence is shown with 95% credible intervals. Rows correspond to individual acoustic traits and columns to the predictors included in the Bayesian mixed-effects models. Tile color represents the posterior mean regression coefficient (green = positive effect, purple = negative effect, white = no effect), while the numerical value within each tile indicates the estimated effect size. Models accounted for non-independence among pairwise comparisons by including multi-membership random effects for populations and individuals, as well as a random intercept for species pair.

Code

# identify all response variables
distance_vars <- grep(
  "_distance$",
  names(sympatric_pairs_simple),
  value = TRUE
)

#remove geographic distance itself
distance_vars <- setdiff(distance_vars, c("geo_distance", "acoustic_distance"))

# scale mst
distance_vars[distance_vars == "mst_distance"] <- "mst_distance_sc"

sympatric_pairs_simple$mst_distance_sc <- scale(sympatric_pairs_simple$mst_distance)

Code

# weak priors
priors <- c(
  
  # Intercept
  prior(normal(0, 5), class = "Intercept"),
  
  # Fixed effect of sympatry
  prior(normal(0, 2), class = "b"),
  
  # Random-effect SDs
  prior(exponential(1), class = "sd"),
  
  # Residual SD
  prior(exponential(1), class = "sigma")
  
)

# store models
sympatry_geo_models <- list()
sympatry_models <- list()
geo_models <- list()

for (resp in distance_vars) {

  cat("\n=============================\n")
  cat("Fitting:", resp, "\n")
  cat("=============================\n")

  # ----------------------------
  # Sympatry + geographic distance
  # ----------------------------

  form_sympatry_geo <- bf(
    as.formula(
      paste0(
        resp,
        " ~ sympatry + scale(geo_distance) + ",
        "(1 | species_pair) + ",
        "(1 | mm(pop1, pop2)) + ",
        "(1 | mm(individual1, individual2))"
      )
    )
  )

  sympatry_geo_models[[resp]] <- brm(
    formula = form_sympatry_geo,
    data = sympatric_pairs_simple,
    family = gaussian(),
    prior = priors,
    cores = 4,
    chains = 4,
    iter = 10000,
    backend = "cmdstanr",
    threads = threading(8),
    control = list(
      adapt_delta = 0.95,
      max_treedepth = 15
    ), file_refit = "always",
    file = paste0(
      "./data/processed/fits/",
      resp,
      "_sympatry_geographic_distance_fit"
    )
  )

  sympatry_geo_models[[resp]] <- add_criterion(
    sympatry_geo_models[[resp]],
    criterion = "loo"
  )

  # ----------------------------
  # Sympatry only
  # ----------------------------

  # form_sympatry <- bf(
  #   as.formula(
  #     paste0(
  #       resp,
  #       " ~ sympatry + ",
  #       "(1 | species_pair) + ",
  #       "(1 | mm(pop1, pop2)) + ",
  #       "(1 | mm(individual1, individual2))"
  #     )
  #   )
  # )
  # 
  # sympatry_models[[resp]] <- brm(
  #   formula = form_sympatry,
  #   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/",
  #     resp,
  #     "_sympatry_fit"
  #   )
  # )
  # 
  # sympatry_models[[resp]] <- add_criterion(
  #   sympatry_models[[resp]],
  #   criterion = "loo"
  # )
  # 
  # # ----------------------------
  # # Geographic distance only
  # # ----------------------------
  # 
  # form_geo <- bf(
  #   as.formula(
  #     paste0(
  #       resp,
  #       " ~ scale(geo_distance) + ",
  #       "(1 | species_pair) + ",
  #       "(1 | mm(pop1, pop2)) + ",
  #       "(1 | mm(individual1, individual2))"
  #     )
  #   )
  # )
  # 
  # geo_models[[resp]] <- brm(
  #   formula = form_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/",
  #     resp,
  #     "_geographic_distance_fit"
  #   )
  # )
  # 
  # geo_models[[resp]] <- add_criterion(
  #   geo_models[[resp]],
  #   criterion = "loo"
  # )
}

Code

# Find and load all saved brms models

model_files <- list.files(
  "./data/processed/fits",
  pattern = "\\.rds$",
  full.names = TRUE
)


model_files <- grep("_complex", model_files, value = TRUE, invert = TRUE)


model_files <- grep(paste(distance_vars, collapse = "|"), model_files, value = TRUE)
model_files <- grep("_distance_sympatry_geographic_distance", model_files, value = TRUE)

plot_brms_heatmap(model_files = model_files)

The expandable section below provides the complete Bayesian model summaries underlying those estimates, including posterior parameter estimates, credible intervals, and convergence diagnostics.

Model results

#::: {.panel-tabset}

Code

for (i in model_files) {

  model_name <- tools::file_path_sans_ext(basename(i))

  model_name <- gsub("_distance_sympatry_geographic_distance_fit", "", model_name)
  
  cat("## ", model_name, "\n\n")

  extended_summary(
    read.file = i,
    highlight = TRUE,
    trace.palette = viridis::mako,
    remove.intercepts = TRUE,
    print.name = FALSE
  )

  cat("\n\n")
}

#:::

Code

# fits <- lapply(model_files, readRDS)
# 
# for (i in fits){
#   extended_summary(
#     i,
#     highlight = TRUE,
#     trace.palette = viridis::mako,
#     remove.intercepts = TRUE,
#     print.name = TRUE
#   )
#     }
# 
# 
# names(fits) <- sub(
#   "\\.rds$",
#   "",
#   basename(model_files)
# )
# 
# # Extract LOOIC values
# 
# # Parse response variable and model type
# 
# model_info <- data.frame(
#   model = names(fits),
#   stringsAsFactors = FALSE
# )
# 
# model_info$response <- NA_character_
# model_info$model_type <- NA_character_
# 
# for(i in seq_len(nrow(model_info))) {
# 
#   m <- model_info$model[i]
# 
#   if(grepl("_sympatry_geographic_distance_fit$", m)) {
# 
#     model_info$response[i] <-
#       sub("_sympatry_geographic_distance_fit$", "", m)
# 
#     model_info$model_type[i] <- "sympatry_geo"
# 
#   } else if(grepl("_geographic_distance_fit$", m)) {
# 
#     model_info$response[i] <-
#       sub("_geographic_distance_fit$", "", m)
# 
#     model_info$model_type[i] <- "geo"
# 
#   } else if(grepl("_sympatry_fit$", m)) {
# 
#     model_info$response[i] <-
#       sub("_sympatry_fit$", "", m)
# 
#     model_info$model_type[i] <- "sympatry"
#   }
# }
# 
# # Compare models using delta ELPD
# 
# responses <- sort(unique(model_info$response))
# 
# loo_comparison <- data.frame(
#   response = responses,
#   best_model = NA_character_,
# 
#   dELPD_sympatry = NA_real_,
#   dELPD_geo = NA_real_,
#   dELPD_sympatry_geo = NA_real_,
# 
#   SE_sympatry = NA_real_,
#   SE_geo = NA_real_,
#   SE_sympatry_geo = NA_real_,
# 
#   stringsAsFactors = FALSE
# )
# 
# responses <- sort(unique(model_info$response))
# 
# for(resp in responses) {
# 
#   cat(
#     paste0(
#       "\n\n###### Response: ",
#       gsub("_distance", "", resp),
#       "\n\n"
#     )
#   )
# 
# 
#   comparison <- loo_compare(
#     fits[[paste0(resp, "_sympatry_geographic_distance_fit")]],
#     fits[[paste0(resp, "_sympatry_fit")]],
#     fits[[paste0(resp, "_geographic_distance_fit")]], 
#     model_names = c("sympatry + geographic distance", "sympatry only", "geographic distance only") 
#   )
# 
#   comparison <- as.data.frame(
#     comparison[, c("elpd_diff", "se_diff")]
#   )
# 
#   comparison$model <- rownames(comparison)
# 
#   kb <- knitr::kable(
#     comparison[, c("model", "elpd_diff", "se_diff")],
#     format = "html",
#     row.names = FALSE,
#     digits = 4
#   )
# 
#   kb <- kableExtra::kable_styling(
#     kb,
#     bootstrap_options = c(
#       "striped",
#       "hover",
#       "condensed",
#       "responsive"
#     )
#   )
# 
#   cat(as.character(kb))
# 
#   # identify best model
#   best_model_name <- comparison$model[1]
# 
#   best_fit <- switch(
#     best_model_name,
#     "sympatry + geographic distance" = fits[[paste0(resp, "_sympatry_geographic_distance_fit")]],
#     "sympatry only" = fits[[paste0(resp, "_sympatry_fit")]],
#     "geographic distance only" = fits[[paste0(resp, "_geographic_distance_fit")]]
#   )
# 
#   cat("\n\n#### Best-supported model\n\n")
# 
#   extended_summary(
#     best_fit,
#     highlight = TRUE,
#     trace.palette = viridis::mako,
#     remove.intercepts = TRUE,
#     print.name = FALSE
#   )
# }

Summary

Overall, sympatry had stronger effects than geographic distance for most acoustic traits, suggesting that coexistence explains acoustic divergence beyond simple geographic separation.
Element duration showed positive effects of both predictors. Sympatric species differed by an additional 170 ms in element duration. A one standard deviation increase in geographic distance was associated with an additional 80 ms difference in element duration.
Frequency range exhibited the strongest positive sympatry effect. Sympatric species differed by an additional 300 Hz in frequency range, whereas a one standard deviation increase in geographic distance was associated with an increase of 60 Hz.
Song duration also increased with both predictors. Sympatric species differed by an additional 140 ms in song duration, while a one standard deviation increase in geographic distance was associated with an increase of 50 ms.
Song rate showed modest positive effects of both predictors. Sympatric species differed by approximately 0.07 elements s⁻¹, while a one standard deviation increase in geographic distance was associated with an increase of 0.06 elements s⁻¹.
Number of elements was only weakly associated with the predictors. Sympatric species differed by approximately 0.05 elements, whereas geographic distance had essentially no effect.
Peak frequency showed negative effects of both predictors. Sympatric species were approximately 60 Hz more similar in peak frequency than allopatric species. Likewise, a one standard deviation increase in geographic distance was associated with a reduction of approximately 20 Hz in mean peak-frequency differences.
Gap duration also showed weak negative effects. Sympatric species were approximately 30 ms more similar in gap duration, while a one standard deviation increase in geographic distance was associated with a reduction of approximately 10 ms in gap-duration differences.
The heterogeneous pattern among song traits indicates that acoustic divergence is trait-specific. Some song features become more divergent in sympatry (particularly frequency range, element duration, and song duration), whereas others remain relatively conserved or even become slightly more similar (mean peak frequency and gap duration).
Collectively, these results suggest that selection associated with sympatry acts on a subset of song characteristics, rather than producing uniform divergence across the entire acoustic phenotype. This pattern is consistent with the hypothesis that only particular components of song contribute strongly to species recognition and reproductive isolation.

2 Complex songs

2.1 Song level

2.1.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")

2.1.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",
                                  "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)

2.1.3 Principal Component Analysis

Code

# run PCA on acoustic variables
pca <- prcomp(
  complex_song_dat[, c(
    "meanpeakf",
    "num.elms",
    "song.duration",
    "song.rate",
    "gap.duration",
    "freq.range.Min5toMax95",
    "mst"
  )],
  center = TRUE,
  scale. = TRUE
)

## Run PCA Inspect variance explained summary(pca)

# plot rotation values by PC
pca_rot <- as.data.frame(pca$rotation[, 1:5])
pca_var <- round(summary(pca)$importance[2, ] * 100)

Code

pca_rot_stck <- stack(pca_rot)

pca_rot_stck$variable <- rownames(pca_rot)
pca_rot_stck$values[pca_rot_stck$ind == "PC1"] <- pca_rot_stck$values[pca_rot_stck$ind ==
    "PC1"]
pca_rot_stck$Sign <- ifelse(pca_rot_stck$values > 0, "Positive", "Negative")
pca_rot_stck$rotation <- abs(pca_rot_stck$values)
pca_rot_stck$ind_var <- paste0(pca_rot_stck$ind, " (", sapply(pca_rot_stck$ind,
    function(x) pca_var[names(pca_var) == x]), "%)")


pca_rot_stck$top_vars <- ave(abs(pca_rot_stck$values), pca_rot_stck$ind,
    FUN = function(x) {

        # Order decreasing
        ord <- order(x, decreasing = TRUE)
        x_sorted <- x[ord]

        # Cumulative proportion
        cumprop <- cumsum(x_sorted)/sum(x_sorted)

        selected_sorted <- cumprop <= 0.5  # variables with 50% of contribution 
        selected_sorted[which(cumprop >= 0.5)[1]] <- TRUE

        # Return logical vector in original order
        selected <- logical(length(x))
        selected[ord] <- selected_sorted

        selected
    })


# Create facet-specific variable
pca_rot_stck$var_facet <- paste(pca_rot_stck$ind, pca_rot_stck$variable,
    sep = "_")

# Reorder within each ind by rotation (largest at top after
# coord_flip)
pca_rot_stck <- do.call(rbind, lapply(split(pca_rot_stck, pca_rot_stck$ind),
    function(df) {

        df$var_facet <- factor(df$var_facet, levels = df$var_facet[order(df$rotation)])

        df
    }))


# add which variables are stable
stable_df <- get_stable_loadings(pca, pcs = 5, B = 1000, cum_threshold = 0.5,
    freq_threshold = 0.5, seed = 123)

pca_rot_stck <- merge(pca_rot_stck, stable_df, by = c("variable",
    "ind"), all.x = TRUE)

pca_rot_stck$top_vars <- ifelse(pca_rot_stck$stable, 1, 1)


# Build colored labels per row

# Colored labels
pca_rot_stck$label_col <- ifelse(pca_rot_stck$stable < 1.1, paste0("<span style='color:black;'>",
    pca_rot_stck$variable, "</span>"), paste0("<span style='color:gray50;'>",
    pca_rot_stck$variable, "</span>"))

# Named vector for labels
label_vec <- setNames(pca_rot_stck$label_col, pca_rot_stck$var_facet)

# absolute CI
pca_rot_stck$ci_low_plot <- abs(pca_rot_stck$ci_lower)
pca_rot_stck$ci_high_plot <- abs(pca_rot_stck$ci_upper)


# Plot
ggplot(pca_rot_stck, aes(x = var_facet, y = rotation, fill = Sign,
    alpha = as.factor(top_vars))) + geom_col() + coord_flip() + scale_alpha_manual(values = pca_rot_stck$top_vars,
    guide = NULL) + scale_x_discrete(labels = label_vec, name = "Variable") +
    labs(x = "Rotation") + scale_fill_viridis_d(alpha = 0.7, begin = 0.2,
    end = 0.8) + facet_wrap(~ind_var, scales = "free_y", nrow = 3) + theme_classic() +
    theme(axis.text.y = element_markdown())

Code

# bind PCA scores with metadata
pca_scores <- cbind(
  pca$x,
  complex_song_dat[, c(
    "Population",
    "species",
    "location",
    "Individual",
    "Lat",
    "Lon",
    "meanpeakf",
    "num.elms",
    "song.duration",
    "song.rate",
    "gap.duration",
    "freq.range.Min5toMax95",
    "mst"
  )]
)

summ_pca <- summary(pca)

# select the PCs explaining at least 90%
pcs <- grep("^PC", names(pca_scores), value = TRUE)[1:min(which(summ_pca$importance[3,] > 0.9))]

pca_scores <- pca_scores |>
  mutate(
    across(all_of(pcs), ~ as.numeric(.x)),
    Lat = as.numeric(Lat),
    Lon = as.numeric(Lon),
    species = as.character(species),
    location = tryCatch(
      iconv(location, from = "", to = "UTF-8"),
      error = function(e) location
    )
  )

2.1.4 Acoustic and geographic distances

2.1.4.1 Build pairwise distance dataset

Code

df_clean_complex <- pca_scores |>
  mutate(
    across(all_of(pcs), as.numeric),
    Lat        = as.numeric(Lat),
    Lon        = as.numeric(Lon),
    species    = as.character(species),
    Individual = as.character(Individual)
  ) |>
  filter(complete.cases(across(all_of(c(pcs, "Lat", "Lon", "species", "Individual")))))

# table(df_clean_complex$Population, df_clean_complex$species)

# assign song IDs
df_clean_complex$song_id <- 1

for (i in 2:nrow(df_clean_complex)) {
  if (df_clean_complex$Individual[i] == df_clean_complex$Individual[i - 1] &&
      df_clean_complex$species[i]    == df_clean_complex$species[i - 1]) {
    df_clean_complex$song_id[i] <- df_clean_complex$song_id[i - 1]
  } else {
    df_clean_complex$song_id[i] <- df_clean_complex$song_id[i - 1] + 1
  }
}

df_clean_complex$song_id <- paste(
  df_clean_complex$species,
  sapply(strsplit(as.character(df_clean_complex$Population), "_"), "[[", 2),
  df_clean_complex$Individual,
  df_clean_complex$song_id,
  sep = "-"
)

2.1.4.2 Convert to long table

Code

# acoustic distance between songs (Euclidean multivariate, unscaled PCs)
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"
))

# 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_id
rownames(dist_geo_mat) <- colnames(dist_geo_mat) <- df_clean_complex$song_id

# use only upper triangle
idx <- which(upper.tri(dist_acoustic_mat), arr.ind = TRUE)

dist_acoustic_long <- data.frame(
  id1      = rownames(dist_acoustic_mat)[idx[, 1]],
  id2      = colnames(dist_acoustic_mat)[idx[, 2]],
  acoustic_distance = dist_acoustic_mat[idx],
  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]
)

# parse species, population, and individual from song IDs
parts1 <- strsplit(as.character(dist_acoustic_long$id1), "-")
parts2 <- strsplit(as.character(dist_acoustic_long$id2), "-")

dist_acoustic_long$species1    <- sapply(parts1, `[`, 1)
dist_acoustic_long$population1 <- sapply(parts1, `[`, 2)
dist_acoustic_long$individual1 <- sapply(parts1, `[`, 3)

dist_acoustic_long$species2    <- sapply(parts2, `[`, 1)
dist_acoustic_long$population2 <- sapply(parts2, `[`, 2)
dist_acoustic_long$individual2 <- sapply(parts2, `[`, 3)

# keep only between-species comparisons
dist_acoustic_long <- dist_acoustic_long[
  dist_acoustic_long$species1 != dist_acoustic_long$species2, ]

# sympatry flag: 1 if same population, 0 otherwise
dist_acoustic_long$sympatry <- as.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 = "_")
)

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.1.5 Statistical analysis

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.
species pair is a random intercept accounting for baseline differences in acoustic divergence among heterospecific species combinations.
mm(population(_1), population(_2)) is a multi-membership random effect accounting for repeated use of the same populations across pairwise comparisons.
mm(individual(_1), individual(_2)) is a multi-membership random effect accounting for repeated use of the same individuals across pairwise comparisons.

Model specifications:

Models were fitted in a Bayesian framework using the brms package.
The response variable was pairwise Euclidean acoustic distance and was modeled using a Gaussian error distribution.
Only heterospecific comparisons were included.
Analyses were restricted to species pairs for which both sympatric and allopatric population comparisons were available. This restriction ensured that sympatry effects were estimated within the same species-pair contrasts rather than being confounded by species pairs occurring exclusively in sympatry or exclusively in allopatry.
Species-pair identity was included as a random intercept to account for inherent differences in acoustic divergence among species combinations.
Population identity was modeled using a multi-membership random effect because each pairwise comparison simultaneously involves two populations, and each population contributes to multiple pairwise distances.
Individual identity was modeled using a multi-membership random effect because each pairwise comparison simultaneously involves two individuals, and each individual contributes to multiple pairwise distances.
Only unique pairwise comparisons were retained and self-comparisons were excluded.
Weakly informative priors were specified for all parameters.
Models were fitted using Hamiltonian Monte Carlo as implemented in Stan through the cmdstanr backend.

The model was fitted with two datasets:

Using only the species pairs that had both sympatric and allopatric comparisons. This directly evaluates if heterospecific songs more divergent in sympatry than in allopatry.
Using the entire dataset, including species pairs that only had sympatric or allopatric comparisons. This is a complementary test of whether sympatric species are generally more divergent than allopatric species.

2.1.5.1 Model fitting

2.1.5.1.1 PCA-based acoustic distance

Prepare data for modeling by restricting to species pairs with both sympatric and allopatric comparisons and creating appropriate random effect structures.

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	EI	ER	MC	Sal
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	EI	ER, MC
Hypoxantha–Palustris	EI	ER, MC
Hypoxantha–Ruficollis	ER	EI, MC, Sal
Iberaensis–Palustris	EI	none
Iberaensis–Ruficollis	none	EI, ER, Sal
Palustris–Ruficollis	none	EI, ER, Sal

Only three species pairs were kept: Hypoxantha_Iberaensis, Hypoxantha_Palustris, Hypoxantha_Ruficollis

2.1.5.1.1.1 Species pairs with both sympatric and allopatric populations

Code

sympatric_pairs_complex <- dist_acoustic_long[
    dist_acoustic_long$species_pair %in%
keep_pairs,
  ]

# create species-specific population IDs
sympatric_pairs_complex$pop1 <- interaction(
  sympatric_pairs_complex$species1,
  sympatric_pairs_complex$population1,
  drop = TRUE
)

sympatric_pairs_complex$pop2 <- interaction(
  sympatric_pairs_complex$species2,
  sympatric_pairs_complex$population2,
  drop = TRUE
)

# make species_pair a factor (fixed effect)
sympatric_pairs_complex$species_pair <- factor(sympatric_pairs_complex$species_pair)

# scale acoustic_distance
sympatric_pairs_complex$acoustic_distance_sc <- scale(sympatric_pairs_complex$acoustic_distance)

Code

# weak priors
priors <- c(
  
  # Intercept
  prior(normal(0, 5), class = "Intercept"),
  
  # Fixed effect of sympatry
  prior(normal(0, 2), class = "b"),
  
  # Random-effect SDs
  prior(exponential(1), class = "sd"),
  
  # Residual SD
  prior(exponential(1), class = "sigma")
  
)

#fit model
sympatry_geo_model_complex <- brm(
  acoustic_distance_sc ~ sympatry + 
    scale(geo_distance) +
    (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"
)
# 
# 
# sympatry_geo_model <- add_criterion(
#     sympatry_geo_model,
#     criterion = "loo"
#   )
# 
# 
# sympatry_model <- brm(
#   acoustic_distance_sc ~ sympatry + 
#     (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_complex_fit"
# )
# 
# 
# sympatry_model <- add_criterion(
#     sympatry_model,
#     criterion = "loo"
#   )
# 
# geo_model <- brm(
#   acoustic_distance_sc ~ scale(geo_distance) +
#     (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_geographic_distance_fit"
#   )
# 
# 
# geo_model <- add_criterion(
#     geo_model,
#     criterion = "loo"
#   )

# beepr::beep(3)

2.1.5.1.1.2 Fit summary

Code

extended_summary(
  sympatry_geo_model_complex,
  highlight = TRUE,
  trace.palette = viridis::mako,
  remove.intercepts = TRUE,
  print.name = FALSE
)

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 2) Intercept-normal(0, 5) sd-exponential(1) sigma-exponential(1)	acoustic_distance_sc ~ sympatry + scale(geo_distance) + (1 \| species_pair) + (1 \| mm(pop1, pop2)) + (1 \| mm(individual1, individual2))	10000	4	1	5000	84 (0.0042%)	0	25366.13	13582.26	820698594

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_sympatry1	-0.122	-0.204	-0.041	1	26506.88	13582.26
b_scalegeo_distance	-0.091	-0.143	-0.039	1	25366.13	13667.09

Summary

2.1.6 Interpretation

Both sympatry and geographic distance had negative effects on overall acoustic distance in the multivariate PCA space, indicating that complex songs tended to become more similar, rather than more divergent, under both conditions.
The effect of sympatry was stronger than that of geographic distance. After accounting for species pair, populations, and individuals, sympatric species exhibited an average reduction of approximately 0.12 standard deviations in multivariate acoustic distance relative to allopatric species.
A one standard deviation increase in geographic distance was associated with an average reduction of approximately 0.10 standard deviations in multivariate acoustic distance.
The 95% credible intervals for both predictors excluded zero, providing strong evidence that the observed effects are consistently negative.
Trace plots showed good mixing among chains and no obvious trends, indicating satisfactory MCMC convergence for both regression coefficients.
These multivariate results agree with the analyses of individual song traits, in which most components of complex songs showed either weak effects or evidence of convergence rather than divergence. Consequently, when considered jointly in multivariate acoustic space, complex songs exhibit an overall pattern of reduced acoustic differentiation among sympatric species.
This contrasts with the results for simple songs, where several acoustic traits exhibited positive effects of sympatry, suggesting that character displacement is expressed primarily in simple songs rather than in complex songs.
Biologically, these findings suggest that the acoustic features composing complex songs may be evolutionarily conserved or constrained, whereas simple songs appear to be the principal target of divergence associated with species coexistence and species recognition.

2.1.6.0.1 Feature-based acoustic distance

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

# scale mst
distance_vars[distance_vars == "mst_distance"] <- "mst_distance_sc"

sympatric_pairs_complex$mst_distance_sc <- scale(sympatric_pairs_complex$mst_distance)

Code

# weak priors
priors <- c(
  
  # Intercept
  prior(normal(0, 5), class = "Intercept"),
  
  # Fixed effect of sympatry
  prior(normal(0, 2), class = "b"),
  
  # Random-effect SDs
  prior(exponential(1), class = "sd"),
  
  # Residual SD
  prior(exponential(1), class = "sigma")
  
)

# store models
sympatry_geo_models <- list()
sympatry_models <- list()
geo_models <- list()

for (resp in distance_vars) {

  cat("\n=============================\n")
  cat("Fitting:", resp, "\n")
  cat("=============================\n")

  # ----------------------------
  # Sympatry + geographic distance
  # ----------------------------

  form_sympatry_geo <- bf(
    as.formula(
      paste0(
        resp,
        " ~ sympatry + scale(geo_distance) + ",
        "(1 | species_pair) + ",
        "(1 | mm(pop1, pop2)) + ",
        "(1 | mm(individual1, individual2))"
      )
    )
  )

  sympatry_geo_models[[resp]] <- brm(
    formula = form_sympatry_geo,
    data = sympatric_pairs_complex,
    family = gaussian(),
    prior = priors,
    cores = 4,
    chains = 4,
    iter = 10000,
    backend = "cmdstanr",
    threads = threading(8),
    control = list(
      adapt_delta = 0.95,
      max_treedepth = 15
    ), file_refit = "always",
    file = paste0(
      "./data/processed/fits/",
      resp,
      "_sympatry_geographic_distance_complex_fit"
    )
  )
  # 
  # sympatry_geo_models[[resp]] <- add_criterion(
  #   sympatry_geo_models[[resp]],
  #   criterion = "loo"
  # )

  # ----------------------------
  # Sympatry only
  # ----------------------------

  # form_sympatry <- bf(
  #   as.formula(
  #     paste0(
  #       resp,
  #       " ~ sympatry + ",
  #       "(1 | species_pair) + ",
  #       "(1 | mm(pop1, pop2)) + ",
  #       "(1 | mm(individual1, individual2))"
  #     )
  #   )
  # )
  # 
  # sympatry_models[[resp]] <- brm(
  #   formula = form_sympatry,
  #   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/",
  #     resp,
  #     "_sympatry_fit"
  #   )
  # )
  # 
  # sympatry_models[[resp]] <- add_criterion(
  #   sympatry_models[[resp]],
  #   criterion = "loo"
  # )
  # 
  # # ----------------------------
  # # Geographic distance only
  # # ----------------------------
  # 
  # form_geo <- bf(
  #   as.formula(
  #     paste0(
  #       resp,
  #       " ~ scale(geo_distance) + ",
  #       "(1 | species_pair) + ",
  #       "(1 | mm(pop1, pop2)) + ",
  #       "(1 | mm(individual1, individual2))"
  #     )
  #   )
  # )
  # 
  # geo_models[[resp]] <- brm(
  #   formula = form_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/",
  #     resp,
  #     "_geographic_distance_fit"
  #   )
  # )
  # 
  # geo_models[[resp]] <- add_criterion(
  #   geo_models[[resp]],
  #   criterion = "loo"
  # )
}

Code

# Find and load all saved brms models

model_files <- list.files(
  "./data/processed/fits",
  pattern = "\\.rds$",
  full.names = TRUE
)

model_files <- grep("_complex", model_files, value = TRUE)

model_files <- grep(paste(distance_vars, collapse = "|"), model_files, value = TRUE)
model_files <- grep("_distance_sympatry_geographic_distance", model_files, value = TRUE)

plot_brms_heatmap(model_files = model_files)

The expandable section below provides the complete Bayesian model summaries underlying those estimates, including posterior parameter estimates, credible intervals, and convergence diagnostics.

Model results

#::: {.panel-tabset}

Code

for (i in model_files) {

  model_name <- tools::file_path_sans_ext(basename(i))

  model_name <- gsub("_distance_sympatry_geographic_distance_fit", "", model_name)
  
  cat("## ", model_name, "\n\n")

  extended_summary(
    read.file = i,
    highlight = TRUE,
    trace.palette = viridis::mako,
    remove.intercepts = TRUE,
    print.name = FALSE
  )

  cat("\n\n")
}

#:::

Code

# fits <- lapply(model_files, readRDS)
# 
# for (i in fits){
#   extended_summary(
#     i,
#     highlight = TRUE,
#     trace.palette = viridis::mako,
#     remove.intercepts = TRUE,
#     print.name = TRUE
#   )
#     }
# 
# 
# names(fits) <- sub(
#   "\\.rds$",
#   "",
#   basename(model_files)
# )
# 
# # Extract LOOIC values
# 
# # Parse response variable and model type
# 
# model_info <- data.frame(
#   model = names(fits),
#   stringsAsFactors = FALSE
# )
# 
# model_info$response <- NA_character_
# model_info$model_type <- NA_character_
# 
# for(i in seq_len(nrow(model_info))) {
# 
#   m <- model_info$model[i]
# 
#   if(grepl("_sympatry_geographic_distance_fit$", m)) {
# 
#     model_info$response[i] <-
#       sub("_sympatry_geographic_distance_fit$", "", m)
# 
#     model_info$model_type[i] <- "sympatry_geo"
# 
#   } else if(grepl("_geographic_distance_fit$", m)) {
# 
#     model_info$response[i] <-
#       sub("_geographic_distance_fit$", "", m)
# 
#     model_info$model_type[i] <- "geo"
# 
#   } else if(grepl("_sympatry_fit$", m)) {
# 
#     model_info$response[i] <-
#       sub("_sympatry_fit$", "", m)
# 
#     model_info$model_type[i] <- "sympatry"
#   }
# }
# 
# # Compare models using delta ELPD
# 
# responses <- sort(unique(model_info$response))
# 
# loo_comparison <- data.frame(
#   response = responses,
#   best_model = NA_character_,
# 
#   dELPD_sympatry = NA_real_,
#   dELPD_geo = NA_real_,
#   dELPD_sympatry_geo = NA_real_,
# 
#   SE_sympatry = NA_real_,
#   SE_geo = NA_real_,
#   SE_sympatry_geo = NA_real_,
# 
#   stringsAsFactors = FALSE
# )
# 
# responses <- sort(unique(model_info$response))
# 
# for(resp in responses) {
# 
#   cat(
#     paste0(
#       "\n\n###### Response: ",
#       gsub("_distance", "", resp),
#       "\n\n"
#     )
#   )
# 
# 
#   comparison <- loo_compare(
#     fits[[paste0(resp, "_sympatry_geographic_distance_fit")]],
#     fits[[paste0(resp, "_sympatry_fit")]],
#     fits[[paste0(resp, "_geographic_distance_fit")]], 
#     model_names = c("sympatry + geographic distance", "sympatry only", "geographic distance only") 
#   )
# 
#   comparison <- as.data.frame(
#     comparison[, c("elpd_diff", "se_diff")]
#   )
# 
#   comparison$model <- rownames(comparison)
# 
#   kb <- knitr::kable(
#     comparison[, c("model", "elpd_diff", "se_diff")],
#     format = "html",
#     row.names = FALSE,
#     digits = 4
#   )
# 
#   kb <- kableExtra::kable_styling(
#     kb,
#     bootstrap_options = c(
#       "striped",
#       "hover",
#       "condensed",
#       "responsive"
#     )
#   )
# 
#   cat(as.character(kb))
# 
#   # identify best model
#   best_model_name <- comparison$model[1]
# 
#   best_fit <- switch(
#     best_model_name,
#     "sympatry + geographic distance" = fits[[paste0(resp, "_sympatry_geographic_distance_fit")]],
#     "sympatry only" = fits[[paste0(resp, "_sympatry_fit")]],
#     "geographic distance only" = fits[[paste0(resp, "_geographic_distance_fit")]]
#   )
# 
#   cat("\n\n#### Best-supported model\n\n")
# 
#   extended_summary(
#     best_fit,
#     highlight = TRUE,
#     trace.palette = viridis::mako,
#     remove.intercepts = TRUE,
#     print.name = FALSE
#   )
# }

Summary

In contrast to the simple-song analyses, both sympatry and geographic distance were associated primarily with reduced acoustic differences among most complex-song traits, suggesting convergence or conservation rather than divergence.
Song rate exhibited the strongest negative effects of both predictors. Sympatric species differed by approximately 0.35 elements s⁻¹ less than allopatric species, while a one standard deviation increase in geographic distance was associated with a reduction of 0.19 elements s⁻¹ in song-rate differences.
Frequency range also showed consistent negative effects. Sympatric species were approximately 150 Hz more similar in frequency range than allopatric species, and a one standard deviation increase in geographic distance was associated with a reduction of 130 Hz in frequency-range differences.
Gap duration was the only trait exhibiting clear positive effects of both predictors. Sympatric species differed by approximately 160 ms more in gap duration than allopatric species, while a one standard deviation increase in geographic distance was associated with an increase of 70 ms.
Song duration showed only weak negative effects, with sympatric species differing by approximately 40 ms less than allopatric species and a one standard deviation increase in geographic distance associated with a reduction of 30 ms.
Number of elements exhibited negligible effects, with sympatric species differing by only 0.02 fewer elements and geographic distance associated with a reduction of approximately 0.01 elements.
Mean peak frequency was essentially unaffected by either predictor, with estimated effects close to zero (≤ 10 Hz).
Overall, these results suggest that most components of complex songs become slightly more similar rather than more divergent in sympatry, in marked contrast to the stronger divergence observed for several simple-song traits.
The only clear exception was gap duration, which became substantially more divergent in sympatry, indicating that this temporal feature may be particularly responsive to selection associated with species coexistence.
Taken together, these findings suggest that character displacement in complex songs is limited and trait-specific, with evidence for divergence confined largely to gap duration, whereas other acoustic characteristics appear to remain conserved or converge among sympatric species.

Session 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-07-15
 pandoc   3.6.3 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/x86_64/ (via rmarkdown)
 quarto   1.8.25 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/quarto

─ Packages ───────────────────────────────────────────────────────────────────
 package          * version  date (UTC) lib source
 abind              1.4-8    2024-09-12 [1] CRAN (R 4.5.2)
 ape                5.8-1    2024-12-16 [1] CRAN (R 4.5.2)
 arrayhelpers       1.1-0    2020-02-04 [1] CRAN (R 4.5.2)
 backports          1.5.1    2026-04-03 [1] CRAN (R 4.5.2)
 bayesplot          1.15.0   2025-12-12 [1] CRAN (R 4.5.2)
 bridgesampling     1.2-1    2025-11-19 [1] CRAN (R 4.5.2)
 brms             * 2.23.0   2025-09-09 [1] CRAN (R 4.5.2)
 brmsish          * 1.0.0    2026-03-03 [1] Github (maRce10/brmsish@81ab826)
 Brobdingnag        1.2-9    2022-10-19 [1] CRAN (R 4.5.2)
 cachem             1.1.0    2024-05-16 [1] CRAN (R 4.5.2)
 checkmate          2.3.4    2026-02-03 [1] CRAN (R 4.5.2)
 class              7.3-23   2025-01-01 [1] CRAN (R 4.5.2)
 classInt           0.4-11   2025-01-08 [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.0.0    2025-08-19 [1] CRAN (R 4.5.2)
 DBI                1.3.0    2026-02-25 [1] CRAN (R 4.5.2)
 deldir             2.0-4    2024-02-28 [1] CRAN (R 4.5.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)
 e1071              1.7-17   2025-12-18 [1] CRAN (R 4.5.2)
 ecodist          * 2.1.3    2023-10-30 [1] CRAN (R 4.5.2)
 ellipsis           0.3.2    2021-04-29 [3] CRAN (R 4.1.1)
 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)
 fs                 2.0.1    2026-03-24 [1] CRAN (R 4.5.2)
 generics           0.1.4    2025-05-09 [1] CRAN (R 4.5.2)
 geosphere        * 1.6-5    2026-03-02 [1] CRAN (R 4.5.2)
 ggdist             3.3.3    2025-04-23 [1] CRAN (R 4.5.2)
 ggplot2          * 4.0.2    2026-02-03 [1] CRAN (R 4.5.2)
 ggtext           * 0.1.2    2022-09-16 [1] CRAN (R 4.5.2)
 glue               1.8.0    2024-09-30 [1] CRAN (R 4.5.2)
 goftest            1.2-3    2021-10-07 [3] CRAN (R 4.1.1)
 gridExtra          2.3      2017-09-09 [3] CRAN (R 4.0.1)
 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.2.3    2026-04-07 [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)
 KernSmooth         2.23-26  2025-01-01 [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)
 mvtnorm            1.3-3    2025-01-10 [1] CRAN (R 4.5.2)
 nicheROVER         1.1.2    2023-10-13 [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)
 PhenotypeSpace   * 0.1.1    2026-03-03 [1] Github (maRce10/PhenotypeSpace@72f146e)
 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.1    2026-04-01 [1] CRAN (R 4.5.2)
 plyr               1.8.9    2023-10-02 [1] CRAN (R 4.5.2)
 polyclip           1.10-7   2024-07-23 [1] CRAN (R 4.5.2)
 posterior          1.6.1    2025-02-27 [1] CRAN (R 4.5.2)
 processx           3.8.7    2026-04-01 [1] CRAN (R 4.5.2)
 proxy              0.4-29   2025-12-29 [1] CRAN (R 4.5.2)
 ps                 1.9.2    2026-03-31 [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)
 raster             3.6-32   2025-03-28 [1] CRAN (R 4.5.2)
 RColorBrewer       1.1-3    2022-04-03 [1] CRAN (R 4.5.2)
 Rcpp             * 1.1.1    2026-01-10 [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.2.0    2026-04-06 [1] CRAN (R 4.5.2)
 rmarkdown          2.31     2026-03-26 [1] CRAN (R 4.5.2)
 rstan              2.32.7   2025-03-10 [1] CRAN (R 4.5.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.1    2025-11-14 [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)
 sf                 1.1-0    2026-02-24 [1] CRAN (R 4.5.2)
 sketchy            1.0.7    2026-03-03 [1] CRANs (R 4.5.2)
 sp                 2.2-1    2026-02-13 [1] CRAN (R 4.5.2)
 spatstat.data      3.1-9    2025-10-18 [1] CRAN (R 4.5.2)
 spatstat.explore   3.7-0    2026-01-22 [1] CRAN (R 4.5.2)
 spatstat.geom      3.7-0    2026-01-20 [1] CRAN (R 4.5.2)
 spatstat.random    3.4-4    2026-01-21 [1] CRAN (R 4.5.2)
 spatstat.sparse    3.1-0    2024-06-21 [1] CRAN (R 4.5.2)
 spatstat.univar    3.1-6    2026-01-17 [1] CRAN (R 4.5.2)
 spatstat.utils     3.2-1    2026-01-10 [1] CRAN (R 4.5.2)
 StanHeaders        2.32.10  2024-07-15 [1] CRAN (R 4.5.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)
 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)
 tensor             1.5.1    2025-06-17 [1] CRAN (R 4.5.2)
 tensorA            0.36.2.1 2023-12-13 [1] CRAN (R 4.5.2)
 terra              1.9-11   2026-03-26 [1] CRAN (R 4.5.2)
 textshaping        1.0.4    2025-10-10 [1] CRAN (R 4.5.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)
 units              1.0-1    2026-03-11 [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-3    2026-03-04 [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.2    2024-10-28 [1] CRAN (R 4.5.2)
 xaringanExtra      0.8.0    2024-05-19 [1] CRAN (R 4.5.2)
 xfun               0.57     2026-03-20 [1] CRAN (R 4.5.2)
 xml2               1.5.2    2026-01-17 [1] CRAN (R 4.5.2)
 yaml               2.3.12   2025-12-10 [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.

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