Statistical analysis

Vocal character displacement in southern capuchinos

Author

Marcelo Araya-Salas

Published

July 2, 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",
  "brmsish",
  "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
  )

  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_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

Code 
# select variables that are not highly correlated
simple_song_dat <- simple_songs_lat_long[, c("Individual", "Lat", "Lon", "species", "Population", "location",
                                  "meanpeakf", "num.elms", "elm.duration", 
                                  "song.duration", "song.rate", "gap.duration",
                                  "freq.range.Min5toMax95", "mst")]

# remove individuals from underrepresented populations
simple_song_dat <- filter(simple_song_dat, Population != "Pal_ER")
simple_song_dat <- filter(simple_song_dat, Population != "NA")

simple_song_dat$Individual  <- factor(simple_song_dat$Individual)
simple_song_dat$species     <- factor(simple_song_dat$species)
simple_song_dat$Population  <- factor(simple_song_dat$Population)
simple_song_dat$location    <- factor(simple_song_dat$location)

1.1.3 Principal Component Analysis

We used the first 5 principal components (PCs) for subsequent analyses, which together explained 93% of the variance in the data.

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


## Run PCA Inspect variance explained summary(pca)

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

pca_rot_stck <- stack(pca_rot)

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


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

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

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

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

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

        selected
    })


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

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

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

        df
    }))


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

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

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


# Build colored labels per row

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

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

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


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

Code 
# bind PCA scores with metadata
pca_scores <- cbind(
  pca$x,
  simple_song_dat[,c("Population", "species", "location", "Individual", "Lat", "Lon", "meanpeakf", "num.elms", "elm.duration", 
                                  "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 <- pca_scores |>
  mutate(
    across(all_of(pcs), as.numeric),
    Lat        = as.numeric(Lat),
    Lon        = as.numeric(Lon),
    species    = as.character(species),
    Individual = as.character(Individual)
  ) |>
  filter(complete.cases(across(all_of(c(pcs, "Lat", "Lon", "species", "Individual")))))

# table(df_clean$Population, df_clean$species)

# assign song IDs
df_clean$song_id <- 1

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

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

1.1.4.2 Convert to long table

Code 
# acoustic distance between songs (Euclidean multivariate, unscaled PCs)
dist_acoustic_mat <- as.matrix(dist(
  df_clean[, pcs],
  method = "euclidean"
))


# and for each feature separately (unscaled)
dist_meanpeakf_mat <- as.matrix(dist(
  df_clean[, "meanpeakf"],
  method = "euclidean"
))

dist_numelms_mat <- as.matrix(dist(
  df_clean[, "num.elms"],
  method = "euclidean"
))

dist_elmduration_mat <- as.matrix(dist(
  df_clean[, "elm.duration"],
  method = "euclidean"
))

dist_songduration_mat <- as.matrix(dist(
  df_clean[, "song.duration"],
  method = "euclidean"
))

dist_songrate_mat <- as.matrix(dist(
  df_clean[, "song.rate"],
  method = "euclidean"
))

dist_gapduration_mat <- as.matrix(dist(
  df_clean[, "gap.duration"],
  method = "euclidean"
))

dist_freqrange_mat <- as.matrix(dist(
  df_clean[, "freq.range.Min5toMax95"],
  method = "euclidean"
))

dist_mst_mat <- as.matrix(dist(
  df_clean[, "mst"],
  method = "euclidean"
))

# geographic distance (Haversine, km) between songs
coords   <- as.matrix(df_clean[, c("Lon", "Lat")])  # Lon first, Lat second
D_geo_km <- geosphere::distm(coords, fun = distHaversine) / 1000
dist_geo <- as.dist(D_geo_km)
dist_geo_mat <- as.matrix(dist_geo)

rownames(dist_acoustic_mat) <- colnames(dist_acoustic_mat) <- df_clean$song_id
rownames(dist_geo_mat) <- colnames(dist_geo_mat) <- df_clean$song_id

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

dist_acoustic_long <- data.frame(
  id1      = rownames(dist_acoustic_mat)[idx[, 1]],
  id2      = colnames(dist_acoustic_mat)[idx[, 2]],
  acoustic_distance = dist_acoustic_mat[idx],
  meanpeakf_distance = dist_meanpeakf_mat[idx],
  numelms_distance = dist_numelms_mat[idx],
  elmduration_distance = dist_elmduration_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 <- dist_acoustic_long[
    dist_acoustic_long$species_pair %in%
keep_pairs,
  ]

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

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

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

# scale acoustic_distance
sympatric_pairs$acoustic_distance_sc <- scale(sympatric_pairs$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,
  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,
  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,
  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
)

    }
formula family priors iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 acoustic_distance_sc ~ sympatry + scale(geo_distance) + (1 | species_pair) + (1 | mm(pop1, pop2)) + (1 | mm(individual1, individual2)) gaussian (identity) b-normal(0, 2) Intercept-normal(0, 5) sd-exponential(1) sigma-exponential(1) 10000 4 1 5000 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_all <- dist_acoustic_long

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

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

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

# scale acoustic_distance
sympatric_pairs_all$acoustic_distance_sc <- scale(sympatric_pairs_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_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,
  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,
  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

  • Estimates also show an increase in acoustic distance associated to increasing geographic distances

  • Although both predictors were positively associated with acoustic divergence, the effect of sympatry was approximately four times larger than the effect of geographic distance (0.30 vs. 0.07)

1.1.5.1.2 Feature-based acoustic distance
Code 
# identify all response variables
distance_vars <- grep(
  "_distance$",
  names(sympatric_pairs),
  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$mst_distance_sc <- scale(sympatric_pairs$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,
    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,
  #   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,
  #   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(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)

Code 
for (i in model_files){
  extended_summary(
    read.file = i,
    highlight = TRUE,
    trace.palette = viridis::mako,
    remove.intercepts = TRUE,
    print.name = TRUE
  )
    }

1.2 elmduration_distance_sympatry_geographic_distance_fit

formula family priors iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 elmduration_distance ~ sympatry + scale(geo_distance) + (1 | species_pair) + (1 | mm(pop1, pop2)) + (1 | mm(individual1, individual2)) gaussian (identity) b-normal(0, 2) Intercept-normal(0, 5) sd-exponential(1) sigma-exponential(1) 10000 4 1 5000 776 (0.039%) 0 30920.06 13745.04 893347603
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sympatry 0.170 0.130 0.210 1 31287.59 14207.05
b_scalegeo_distance 0.076 0.055 0.096 1 30920.06 13745.04

1.3 freqrange_distance_sympatry_geographic_distance_fit

formula family priors iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 freqrange_distance ~ sympatry + scale(geo_distance) + (1 | species_pair) + (1 | mm(pop1, pop2)) + (1 | mm(individual1, individual2)) gaussian (identity) b-normal(0, 2) Intercept-normal(0, 5) sd-exponential(1) sigma-exponential(1) 10000 4 1 5000 841 (0.042%) 0 28988.67 13981.01 1672401102
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sympatry1 0.299 0.242 0.356 1 28988.67 13981.01
b_scalegeo_distance -0.057 -0.085 -0.027 1 29232.14 14102.45

1.4 gapduration_distance_sympatry_geographic_distance_fit

formula family priors iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 gapduration_distance ~ sympatry + scale(geo_distance) + (1 | species_pair) + (1 | mm(pop1, pop2)) + (1 | mm(individual1, individual2)) gaussian (identity) b-normal(0, 2) Intercept-normal(0, 5) sd-exponential(1) sigma-exponential(1) 10000 4 1 5000 600 (0.03%) 0 25338.56 12964.53 1761610196
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sympatry1 -0.027 -0.092 0.040 1 25338.56 13794.21
b_scalegeo_distance -0.009 -0.043 0.025 1 25712.83 12964.53

1.5 meanpeakf_distance_sympatry_geographic_distance_fit

formula family priors iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 meanpeakf_distance ~ sympatry + scale(geo_distance) + (1 | species_pair) + (1 | mm(pop1, pop2)) + (1 | mm(individual1, individual2)) gaussian (identity) b-normal(0, 2) Intercept-normal(0, 5) sd-exponential(1) sigma-exponential(1) 10000 4 1 5000 2573 (0.13%) 0 307.037 11754.23 1871852522
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sympatry -0.058 -0.075 -0.041 1.007 754.049 11754.23
b_scalegeo_distance -0.021 -0.030 -0.012 1.011 307.037 12904.09

1.6 numelms_distance_sympatry_geographic_distance_fit

formula family priors iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 numelms_distance ~ sympatry + scale(geo_distance) + (1 | species_pair) + (1 | mm(pop1, pop2)) + (1 | mm(individual1, individual2)) gaussian (identity) b-normal(0, 2) Intercept-normal(0, 5) sd-exponential(1) sigma-exponential(1) 10000 4 1 5000 515 (0.026%) 0 20109.15 13287.54 1761610196
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sympatry 0.052 -0.011 0.116 1 20386.57 13608.95
b_scalegeo_distance 0.003 -0.030 0.036 1 20109.15 13287.54

1.7 songduration_distance_sympatry_geographic_distance_fit

formula family priors iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 songduration_distance ~ sympatry + scale(geo_distance) + (1 | species_pair) + (1 | mm(pop1, pop2)) + (1 | mm(individual1, individual2)) gaussian (identity) b-normal(0, 2) Intercept-normal(0, 5) sd-exponential(1) sigma-exponential(1) 10000 4 1 5000 885 (0.044%) 0 24130.38 13621.9 376777791
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sympatry1 0.140 0.079 0.201 1 24279.21 13626.4
b_scalegeo_distance 0.055 0.023 0.086 1 24130.38 13621.9

1.8 songrate_distance_sympatry_geographic_distance_fit

formula family priors iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 songrate_distance ~ sympatry + scale(geo_distance) + (1 | species_pair) + (1 | mm(pop1, pop2)) + (1 | mm(individual1, individual2)) gaussian (identity) b-normal(0, 2) Intercept-normal(0, 5) sd-exponential(1) sigma-exponential(1) 10000 4 1 5000 687 (0.034%) 0 27433.37 13498.67 1296828965
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
b_sympatry1 0.069 0.016 0.122 1 27595.15 14299.11
b_scalegeo_distance 0.060 0.032 0.087 1.001 27433.37 13498.67

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
#   )
# }
─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.5.2 (2025-10-31)
 os       Ubuntu 22.04.5 LTS
 system   x86_64, linux-gnu
 ui       X11
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       America/Costa_Rica
 date     2026-07-02
 pandoc   3.2 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/x86_64/ (via rmarkdown)
 quarto   1.7.31 @ /usr/local/bin/quarto

─ Packages ───────────────────────────────────────────────────────────────────
 package          * version   date (UTC) lib source
 abind              1.4-8     2024-09-12 [1] CRAN (R 4.5.0)
 ape                5.8-1     2024-12-16 [1] CRAN (R 4.5.0)
 arrayhelpers       1.1-0     2020-02-04 [1] CRAN (R 4.5.0)
 backports          1.5.1     2026-04-03 [1] CRAN (R 4.5.2)
 bayesplot          1.15.0    2025-12-12 [1] CRAN (R 4.5.2)
 bridgesampling     1.2-1     2025-11-19 [1] CRAN (R 4.5.0)
 brms             * 2.23.0    2025-09-09 [1] CRAN (R 4.5.0)
 brmsish          * 1.0.0     2026-01-06 [1] CRANs (R 4.5.2)
 Brobdingnag        1.2-9     2022-10-19 [1] CRAN (R 4.5.0)
 cachem             1.1.0     2024-05-16 [1] CRAN (R 4.5.0)
 checkmate          2.3.4     2026-02-03 [1] CRAN (R 4.5.2)
 class              7.3-23    2025-01-01 [4] CRAN (R 4.4.2)
 classInt           0.4-11    2025-01-08 [1] CRAN (R 4.5.0)
 cli                3.6.6     2026-04-09 [1] CRAN (R 4.5.2)
 cluster            2.1.8.2   2026-02-05 [4] CRAN (R 4.5.2)
 cmdstanr         * 0.9.0     2025-08-21 [1] https://stan-dev.r-universe.dev (R 4.5.0)
 coda               0.19-4.1  2024-01-31 [1] CRAN (R 4.5.0)
 codetools          0.2-20    2024-03-31 [4] CRAN (R 4.5.0)
 commonmark         1.9.5     2025-03-17 [3] CRAN (R 4.5.0)
 cowplot            1.2.0     2025-07-07 [1] CRAN (R 4.5.0)
 crayon             1.5.3     2024-06-20 [1] CRAN (R 4.5.0)
 curl               7.1.0     2026-04-22 [1] CRAN (R 4.5.2)
 DBI                1.3.0     2026-02-25 [1] CRAN (R 4.5.2)
 deldir             2.0-4     2024-02-28 [1] CRAN (R 4.5.0)
 devtools           2.4.5     2022-10-11 [1] CRAN (R 4.5.0)
 digest             0.6.39    2025-11-19 [1] CRAN (R 4.5.0)
 distributional     0.5.0     2024-09-17 [1] CRAN (R 4.5.0)
 dplyr            * 1.1.4     2023-11-17 [1] CRAN (R 4.5.0)
 e1071              1.7-17    2025-12-18 [1] CRAN (R 4.5.2)
 ecodist          * 2.1.3     2023-10-30 [1] CRAN (R 4.5.0)
 ellipsis           0.3.2     2021-04-29 [3] CRAN (R 4.1.1)
 emmeans            1.11.1    2025-05-04 [3] CRAN (R 4.5.0)
 estimability       1.5.1     2024-05-12 [3] CRAN (R 4.5.0)
 evaluate           1.0.5     2025-08-27 [1] CRAN (R 4.5.0)
 farver             2.1.2     2024-05-13 [1] CRAN (R 4.5.0)
 fastmap            1.2.0     2024-05-15 [1] CRAN (R 4.5.0)
 forcats          * 1.0.0     2023-01-29 [1] CRAN (R 4.5.0)
 fs                 2.1.0     2026-04-18 [1] CRAN (R 4.5.2)
 generics           0.1.4     2025-05-09 [1] CRAN (R 4.5.0)
 geosphere        * 1.5-14    2021-10-13 [3] CRAN (R 4.1.1)
 ggdist             3.3.3     2025-04-23 [1] CRAN (R 4.5.0)
 ggplot2          * 4.0.3     2026-04-22 [1] CRAN (R 4.5.2)
 ggtext           * 0.1.2     2022-09-16 [1] RSPM (R 4.5.0)
 glue               1.8.1     2026-04-17 [1] CRAN (R 4.5.2)
 goftest            1.2-3     2021-10-07 [1] CRAN (R 4.5.2)
 gridExtra          2.3       2017-09-09 [1] CRAN (R 4.5.0)
 gridtext           0.1.5     2022-09-16 [1] RSPM (R 4.5.0)
 gtable             0.3.6     2024-10-25 [1] CRAN (R 4.5.0)
 hms                1.1.3     2023-03-21 [3] CRAN (R 4.5.0)
 htmltools          0.5.9     2025-12-04 [1] CRAN (R 4.5.0)
 htmlwidgets        1.6.4     2023-12-06 [1] RSPM (R 4.5.0)
 httpuv             1.6.16    2025-04-16 [1] RSPM (R 4.5.0)
 igraph             2.2.1     2025-10-27 [1] CRAN (R 4.5.0)
 inline             0.3.21    2025-01-09 [1] CRAN (R 4.5.0)
 jsonlite           2.0.0     2025-03-27 [1] CRAN (R 4.5.0)
 kableExtra       * 1.4.0     2024-01-24 [1] CRAN (R 4.5.0)
 KernSmooth         2.23-26   2025-01-01 [4] CRAN (R 4.4.2)
 knitr            * 1.51      2025-12-20 [1] CRAN (R 4.5.2)
 labeling           0.4.3     2023-08-29 [1] CRAN (R 4.5.0)
 later              1.4.2     2025-04-08 [1] RSPM (R 4.5.0)
 lattice            0.22-9    2026-02-09 [4] CRAN (R 4.5.2)
 lifecycle          1.0.5     2026-01-08 [1] CRAN (R 4.5.2)
 litedown           0.7       2025-04-08 [1] CRAN (R 4.5.0)
 loo                2.9.0     2025-12-23 [1] CRAN (R 4.5.2)
 lubridate        * 1.9.5     2026-02-04 [1] CRAN (R 4.5.2)
 magrittr           2.0.5     2026-04-04 [1] CRAN (R 4.5.2)
 markdown           2.0       2025-03-23 [1] CRAN (R 4.5.0)
 MASS               7.3-65    2025-02-28 [4] CRAN (R 4.4.3)
 Matrix             1.7-4     2025-08-28 [4] CRAN (R 4.5.1)
 matrixStats        1.5.0     2025-01-07 [1] CRAN (R 4.5.0)
 memoise            2.0.1     2021-11-26 [1] CRAN (R 4.5.0)
 mgcv               1.9-4     2025-11-07 [4] CRAN (R 4.5.2)
 mime               0.13      2025-03-17 [1] CRAN (R 4.5.0)
 miniUI             0.1.2     2025-04-17 [3] CRAN (R 4.5.0)
 multcomp           1.4-28    2025-01-29 [3] CRAN (R 4.5.0)
 mvtnorm            1.3-3     2025-01-10 [1] CRAN (R 4.5.0)
 nicheROVER         1.1.2     2023-10-13 [1] CRAN (R 4.5.0)
 nlme               3.1-168   2025-03-31 [4] CRAN (R 4.4.3)
 otel               0.2.0     2025-08-29 [1] CRAN (R 4.5.2)
 packrat            0.9.3     2025-06-16 [1] CRAN (R 4.5.2)
 pbapply            1.7-4     2025-07-20 [1] CRAN (R 4.5.0)
 permute          * 0.9-10    2026-02-06 [1] CRAN (R 4.5.2)
 PhenotypeSpace   * 0.1.1     2026-02-25 [1] CRAN (R 4.5.2)
 pillar             1.11.1    2025-09-17 [1] CRAN (R 4.5.0)
 pkgbuild           1.4.8     2025-05-26 [1] CRAN (R 4.5.0)
 pkgconfig          2.0.3     2019-09-22 [1] CRAN (R 4.5.0)
 pkgload            1.4.1     2025-09-23 [1] CRAN (R 4.5.0)
 plyr               1.8.9     2023-10-02 [1] CRAN (R 4.5.0)
 polyclip           1.10-7    2024-07-23 [1] CRAN (R 4.5.0)
 posterior          1.6.1     2025-02-27 [1] CRAN (R 4.5.0)
 processx           3.8.6     2025-02-21 [1] CRAN (R 4.5.0)
 profvis            0.4.0     2024-09-20 [1] CRAN (R 4.5.0)
 promises           1.3.3     2025-05-29 [1] RSPM (R 4.5.0)
 proxy              0.4-29    2025-12-29 [1] CRAN (R 4.5.2)
 ps                 1.9.1     2025-04-12 [1] CRAN (R 4.5.0)
 purrr            * 1.2.0     2025-11-04 [1] CRAN (R 4.5.0)
 QuickJSR           1.8.1     2025-09-20 [1] CRAN (R 4.5.0)
 R6                 2.6.1     2025-02-15 [1] CRAN (R 4.5.0)
 raster             3.6-32    2025-03-28 [1] CRAN (R 4.5.0)
 RColorBrewer       1.1-3     2022-04-03 [1] CRAN (R 4.5.0)
 Rcpp             * 1.1.1-1.1 2026-04-24 [1] CRAN (R 4.5.2)
 RcppParallel       5.1.11-1  2025-08-27 [1] CRAN (R 4.5.0)
 readr            * 2.1.5     2024-01-10 [1] CRAN (R 4.5.0)
 remotes            2.5.0     2024-03-17 [1] CRAN (R 4.5.0)
 reshape2           1.4.5     2025-11-12 [1] CRAN (R 4.5.0)
 rlang              1.2.0     2026-04-06 [1] CRAN (R 4.5.2)
 rmarkdown          2.31      2026-03-26 [1] CRAN (R 4.5.2)
 rstan              2.32.7    2025-03-10 [1] CRAN (R 4.5.0)
 rstantools         2.5.0     2025-09-01 [1] CRAN (R 4.5.0)
 rstudioapi         0.17.1    2024-10-22 [1] CRAN (R 4.5.0)
 S7                 0.2.2     2026-04-22 [1] CRAN (R 4.5.2)
 sandwich           3.1-1     2024-09-15 [3] CRAN (R 4.5.0)
 scales             1.4.0     2025-04-24 [1] CRAN (R 4.5.0)
 sessioninfo        1.2.3     2025-02-05 [3] CRAN (R 4.5.0)
 sf                 1.1-0     2026-02-24 [1] CRAN (R 4.5.2)
 shiny              1.10.0    2024-12-14 [1] CRAN (R 4.5.0)
 sketchy            1.0.7     2026-01-28 [1] CRANs (R 4.5.2)
 sp                 2.2-1     2026-02-13 [1] CRAN (R 4.5.2)
 spatstat.data      3.1-9     2025-10-18 [1] CRAN (R 4.5.2)
 spatstat.explore   3.7-0     2026-01-22 [1] CRAN (R 4.5.2)
 spatstat.geom      3.7-0     2026-01-20 [1] CRAN (R 4.5.2)
 spatstat.random    3.4-4     2026-01-21 [1] CRAN (R 4.5.2)
 spatstat.sparse    3.1-0     2024-06-21 [1] CRAN (R 4.5.0)
 spatstat.univar    3.1-6     2026-01-17 [1] CRAN (R 4.5.2)
 spatstat.utils     3.2-1     2026-01-10 [1] CRAN (R 4.5.2)
 StanHeaders        2.32.10   2024-07-15 [1] CRAN (R 4.5.0)
 stringi            1.8.7     2025-03-27 [1] CRAN (R 4.5.0)
 stringr          * 1.6.0     2025-11-04 [1] CRAN (R 4.5.0)
 survival           3.8-6     2026-01-16 [4] CRAN (R 4.5.2)
 svglite            2.2.2     2025-10-21 [1] CRAN (R 4.5.0)
 svUnit             1.0.8     2025-08-26 [1] CRAN (R 4.5.0)
 systemfonts        1.3.1     2025-10-01 [1] CRAN (R 4.5.0)
 tensor             1.5.1     2025-06-17 [1] CRAN (R 4.5.0)
 tensorA            0.36.2.1  2023-12-13 [1] CRAN (R 4.5.0)
 terra              1.9-11    2026-03-26 [1] CRAN (R 4.5.2)
 textshaping        1.0.4     2025-10-10 [1] CRAN (R 4.5.0)
 TH.data            1.1-3     2025-01-17 [3] CRAN (R 4.5.0)
 tibble           * 3.3.0     2025-06-08 [1] RSPM (R 4.5.0)
 tidybayes          3.0.7     2024-09-15 [1] CRAN (R 4.5.0)
 tidyr            * 1.3.2     2025-12-19 [1] CRAN (R 4.5.2)
 tidyselect         1.2.1     2024-03-11 [1] CRAN (R 4.5.0)
 tidyverse        * 2.0.0     2023-02-22 [1] RSPM (R 4.5.0)
 timechange         0.4.0     2026-01-29 [1] CRAN (R 4.5.2)
 tzdb               0.5.0     2025-03-15 [1] CRAN (R 4.5.0)
 units              1.0-1     2026-03-11 [1] CRAN (R 4.5.2)
 urlchecker         1.0.1     2021-11-30 [1] CRAN (R 4.5.0)
 usethis            3.1.0     2024-11-26 [3] CRAN (R 4.5.0)
 V8                 6.0.4     2025-06-04 [1] RSPM (R 4.5.0)
 vctrs              0.7.3     2026-04-11 [1] CRAN (R 4.5.2)
 vegan            * 2.7-2     2025-10-08 [1] CRAN (R 4.5.2)
 viridis            0.6.5     2024-01-29 [1] CRAN (R 4.5.0)
 viridisLite        0.4.3     2026-02-04 [1] CRAN (R 4.5.2)
 withr              3.0.2     2024-10-28 [1] CRAN (R 4.5.0)
 xaringanExtra      0.8.0     2024-05-19 [1] CRAN (R 4.5.0)
 xfun               0.57      2026-03-20 [1] CRAN (R 4.5.2)
 xml2               1.5.2     2026-01-17 [1] CRAN (R 4.5.2)
 xtable             1.8-8     2026-02-22 [1] CRAN (R 4.5.2)
 yaml               2.3.12    2025-12-10 [1] CRAN (R 4.5.2)
 zoo                1.8-14    2025-04-10 [3] CRAN (R 4.5.0)

 [1] /home/m/R/x86_64-pc-linux-gnu-library/4.5
 [2] /usr/local/lib/R/site-library
 [3] /usr/lib/R/site-library
 [4] /usr/lib/R/library
 * ── Packages attached to the search path.

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