Statistical analysis

Owl facial disc evolution

Author

Published

May 20, 2026

Source code and data found at https://github.com/maRce10/owl_vocalization_comparative_analysis

Code

# options to customize chunk outputs
knitr::opts_chunk$set(
  tidy.opts = list(width.cutoff = 65), 
  tidy = TRUE,
  message = FALSE
 )

Purpose

Infer the evolutionary scenario favoring facial disc in owls

Analysis flowchart

flowchart LR
  A[PCA on\nsong features] --> B(Define DAGs) 
  B --> C(Statistical analysis)
  C --> D(Compare model fit)
  D --> E(Model summary) 

style A fill:#382A5433
style B fill:#395D9C33
style C fill:#3497A933
style D fill:#60CEAC33

Load packages

Code

# knitr is require for creating html/pdf/word reports formatR is
# used for soft-wrapping code

# install/ load packages
sketchy::load_packages(packages = c("ggplot2", "viridis", "dagitty",
    "ggdag", "brms", "brmsish", "rstan", "ggraph", "ggparty", "rncl",
    "ape", "ggtext", "phangorn", "corrplot"))

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

1 DAGS

1.1 DAG 1 — Song influences facial disc

Code

dag_song_to_disc <- dagify(

    # background variables
    habitat_structure ~ latitude,
    niche_width ~ latitude,

    activity_rhythm ~ body_size,
    niche_width ~ body_size,

    # ecological relationships
    activity_rhythm ~ habitat_structure,
    niche_width ~ activity_rhythm,

    # ecology affects song
    song_structure ~ habitat_structure,
    song_structure ~ activity_rhythm,
    song_structure ~ niche_width,

    # ecology directly affects facial disc
    facial_disc ~ habitat_structure,
    facial_disc ~ activity_rhythm,
    facial_disc ~ niche_width,

    # communication hypothesis
    facial_disc ~ song_structure,

    labels = c(
        latitude = "Latitude",
        body_size = "Body\nweight/size",

        habitat_structure = "Habitat\nstructure\n(tree coverage)",
        activity_rhythm = "Activity rhythm",
        niche_width = "Niche width",

        song_structure = "Song\nstructure",

        facial_disc = "Facial disc"
    )
)

# tidy DAG
set.seed(5)

tidy_dag1 <- tidy_dagitty(dag_song_to_disc)

# node classes
tidy_dag1$data$type <- "predictor"

tidy_dag1$data$type[
    tidy_dag1$data$name %in% c("latitude", "body_size")
] <- "background"

tidy_dag1$data$type[
    tidy_dag1$data$name %in% c(
        "habitat_structure",
        "activity_rhythm",
        "niche_width"
    )
] <- "ecology"

tidy_dag1$data$type[
    tidy_dag1$data$name %in% c("song_structure")
] <- "communication"

tidy_dag1$data$type[
    tidy_dag1$data$name %in% c("facial_disc")
] <- "outcome"

# plot
ggplot(
    tidy_dag1,
    aes(
        x = x,
        y = y,
        xend = xend,
        yend = yend
    )
) +

    scale_color_viridis_d(
        option = "G",
        begin = 0.2,
        end = 0.8,
        alpha = 0.8
    ) +

    geom_dag_edges_fan(
        edge_color = mako(10, alpha = 0.5)[2],
        arrow = grid::arrow(
            length = grid::unit(10, "pt"),
            type = "closed"
        )
    ) +

    geom_dag_point(
        aes(color = type),
        size = 24
    ) +

    geom_dag_text(
        aes(label = label),
        color = "black",
        size = 4
    ) +

    theme_dag() +

    ggtitle("DAG 1: Song structure influence facial disc") +

    theme(
        legend.position = "bottom"
    )

1.2 DAG 2 — Facial disc influences song

Code

dag_disc_to_song <- dagify(

    # background variables
    habitat_structure ~ latitude,
    niche_width ~ latitude,

    activity_rhythm ~ body_size,
    niche_width ~ body_size,

    # ecological relationships
    activity_rhythm ~ habitat_structure,
    niche_width ~ activity_rhythm,

    # ecology affects facial disc
    facial_disc ~ habitat_structure,
    facial_disc ~ activity_rhythm,
    facial_disc ~ niche_width,

    # ecology affects songs
    song_structure ~ habitat_structure,
    song_structure ~ activity_rhythm,
    song_structure ~ niche_width,

    # alternative hypothesis
    song_structure ~ facial_disc,

    labels = c(
        latitude = "Latitude",
        body_size = "Body\nweight/size",

        habitat_structure = "Habitat\nstructure\n(tree coverage)",
        activity_rhythm = "Activity rhythm",
        niche_width = "Niche width",

        song_structure = "Song\nstructure",

        facial_disc = "Facial disc"
    )
)

# tidy DAG
tidy_dag2 <- tidy_dagitty(dag_disc_to_song)

# node classes
tidy_dag2$data$type <- "predictor"

tidy_dag2$data$type[
    tidy_dag2$data$name %in% c("latitude", "body_size")
] <- "background"

tidy_dag2$data$type[
    tidy_dag2$data$name %in% c(
        "habitat_structure",
        "activity_rhythm",
        "niche_width"
    )
] <- "ecology"

tidy_dag2$data$type[
    tidy_dag2$data$name %in% c("song_structure")
] <- "communication"

tidy_dag2$data$type[
    tidy_dag2$data$name %in% c("facial_disc")
] <- "outcome"

# plot
ggplot(
    tidy_dag2,
    aes(
        x = x,
        y = y,
        xend = xend,
        yend = yend
    )
) +

    scale_color_viridis_d(
        option = "G",
        begin = 0.2,
        end = 0.8,
        alpha = 0.8
    ) +

    geom_dag_edges_fan(
        edge_color = mako(10, alpha = 0.5)[2],
        arrow = grid::arrow(
            length = grid::unit(10, "pt"),
            type = "closed"
        )
    ) +

    geom_dag_point(
        aes(color = type),
        size = 24
    ) +

    geom_dag_text(
        aes(label = label),
        color = "black",
        size = 4
    ) +

    theme_dag() +

    ggtitle("DAG 2: Facial disc influences song structure") +

    theme(
        legend.position = "bottom"
    )

2 Reduce song dimensionality

PCA song features

Code

dat <- readxl::read_excel("./data/processed/final_database_disc_comparative.xlsx",
    na = "NA")

# remove rows with missing categorical variables (cannot
# currently be modeled with mi())
dat <- dat |>
    filter(!is.na(cover2), !is.na(activity), !is.na(disc_pres))


trees <- read_nexus_phylo("./data/raw/owls206.nex")

# setdiff(dat$Species, trees[[1]]$tip.label)
# setdiff(trees[[1]]$tip.label, dat$Species)

trees <- drop.tip.multiPhylo(trees, setdiff(trees[[1]]$tip.label,
    dat$Species))

# reorder data to match trees
dat <- dat[match(trees[[1]]$tip.label, dat$Species), ]

# verify matching all(trees[[1]]$tip.label == dat$Species)

## Data prep

# binary activity variable
dat$activity <- ifelse(dat$activity == 1, 1, 0)

# ordered habitat variable
dat$cover2 <- ordered(dat$cover2)

# binary facial disc response
dat$facial_disc <- ifelse(dat$disc_pres == 1, 1, 0)

# scale continuous predictors
dat$log_weight_sc <- scale(log(dat$weight))[, 1]
dat$latitude_sc <- scale(dat$latitude)[, 1]
dat$niche_width_sc <- scale(dat$niche_width)[, 1]

# sample 100 trees
set.seed(123)
subtree <- sample(trees, 30)

chains = 1
cores = 4
iters = 2000

Features that combined contributed to 50% or more of the variance in song structure are highlighted:

Code

## Song variables for PCA

song_feats <- c("peak_frequency", "frequency_range", "modulation",
    "song_duration", "num_elms", "elm_types", "umap_space_size")

complete_cases <- complete.cases(dat[, song_feats])

## Run PCA

pca <- prcomp(dat[complete_cases, song_feats], center = TRUE, scale. = TRUE)

## Run PCA Inspect variance explained summary(pca)

# plot rotation values by PC
pca_rot <- as.data.frame(pca$rotation[, 1:4])
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 = 4, 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, 0.6, 1)


# Build colored labels per row

# Colored labels
pca_rot_stck$label_col <- ifelse(pca_rot_stck$stable, 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") + theme_classic() +
    theme(axis.text.y = element_markdown())

Code

## PCA loadings

dat$song_pc1_freq <- NA
dat$song_pc2_diver <- NA
dat$song_pc4_length <- NA

dat$song_pc1_freq[complete_cases] <- pca$x[, 1]
dat$song_pc2_diver[complete_cases] <- pca$x[, 2]
dat$song_pc4_length[complete_cases] <- pca$x[, 4]

## Standardize PC1
dat$song_pc1_freq_sc <- scale(dat$song_pc1_freq)[, 1]
dat$song_pc2_diver_sc <- scale(dat$song_pc2_diver)[, 1]
dat$song_pc4_length_sc <- scale(dat$song_pc4_length)[, 1]

Takeaways

PC1 was mainly driven by frequency features, PC2 by acoustic diversity features, and PC4 by length features
These 3 PCs were used as song features in the subsequent analyses

3 Fit phylogenetic SEM models

Two alternative phylogenetic structural equation models (SEMs) were fitted to evaluate competing hypotheses regarding the causal relationship between song traits and facial disc morphology in owls. Both models shared the same ecological and allometric background structure, differing only in the directionality of the relationship between song features and facial disc traits.

Shared model structure:

[ \[\begin{split} \text{niche width} &\sim \text{latitude} + \text{body size} \\ \\ \text{activity rhythm} &\sim \text{monotonic(habitat structure)} + \text{body size} \end{split}\]

]

Alternative causal hypotheses:

Song → facial disc hypothesis

[ \[\begin{split} \text{song feature} &\sim \text{monotonic(habitat structure)} + \\ &\quad \text{monotonic(activity rhythm)} + \text{niche width} \\ \\ \text{facial disc} &\sim \text{monotonic(habitat structure)} + \\ &\quad \text{monotonic(activity rhythm)} + \\ &\quad \text{niche width} + \text{song feature} \end{split}\]

]

Facial disc → song hypothesis

[ \[\begin{split} \text{facial disc} &\sim \text{monotonic(habitat structure)} + \\ &\quad \text{monotonic(activity rhythm)} + \text{niche width} \\ \\ \text{song feature} &\sim \text{monotonic(habitat structure)} + \\ &\quad \text{monotonic(activity rhythm)} + \\ &\quad \text{niche width} + \text{facial disc} \end{split}\]

]

Specifications:

Models were implemented as Bayesian piecewise structural equation models (SEMs) using the brms package
Phylogenetic relationships among species were incorporated as a random effect using a maximum-clade credibility phylogeny derived from the 1000 phylogenies available
Habitat structure (cover2) and activity rhythm were modeled as ordered predictors using monotonic effects (mo())
Niche width was modeled as a Gaussian response variable
Activity rhythm was modeled using a cumulative ordinal distribution
Facial disc presence was modeled as a Bernoulli response variable
Different models were fitted for each PCA-derived song features (e.g., frequency, temporal features) as the outcome variable, modeled as Gaussian responses
Latitude and log-transformed body size were included as background predictors to account for large-scale ecological and allometric effects
Residual correlations among submodels were fixed to zero using set_rescor(FALSE) to preserve the directed acyclic graph structure
Weakly informative priors were used for all model parameters
Continuous predictors were standardized (z-transformed) prior to analysis
Model fit and relative support among alternative causal models were evaluated using leave-one-out cross-validation (LOO) criteria

3.1 Fit SEM models over a single consensus tree

3.1.1 On PC derived song features

Code

## Estimate MCC tree
# subtree should be a multiPhylo object
# containing posterior phylogenetic trees

mcc_tree <- maxCladeCred(subtree)

## Build phylogenetic covariance matrix

vcv.phylo <- ape::vcv.phylo(
  mcc_tree,
  corr = TRUE
)

## Song variables

features <- c(
  "song_pc1_freq_sc",
  "song_pc2_diver_sc",
  "song_pc4_length_sc"
)

chains = 1
cores = 4
iters = 2000

## Loop over song variables

for (song_var in features) {

  cat("\n========================\n")
  cat("Running:", song_var, "\n")
  cat("========================\n")

  ## Response names for priors

  resp_name <- gsub("_", "", song_var)

  ## Common models

  bf_niche <- bf(
    niche_width_sc | mi() ~
      mi(latitude_sc) +
      mi(log_weight_sc) +
      (1 | gr(Species, cov = vcv.phylo))
  )

  bf_activity <- bf(
    activity ~
      mo(cover2) +
      mi(log_weight_sc) +
      (1 | gr(Species, cov = vcv.phylo)),
    family = bernoulli()
  )

  bf_latitude <- bf(
    latitude_sc | mi() ~ 1
  )

  bf_weight <- bf(
    log_weight_sc | mi() ~ 1
  )

  ## Priors
  
  priors <- c(

    # niche width model
    prior(normal(0, 1), class = "b", resp = "nichewidthsc"),
    prior(normal(0, 1.5), class = "Intercept", resp = "nichewidthsc"),
    prior(exponential(1), class = "sigma", resp = "nichewidthsc"),

    # activity model
    prior(normal(0, 1), class = "b", resp = "activity"),
    prior(normal(0, 1.5), class = "Intercept", resp = "activity"),

    # song model
    do.call(
      brms::prior,
      list(
        "normal(0, 1)",
        class = "b",
        resp = resp_name
      )
    ),

    do.call(
      brms::prior,
      list(
        "normal(0, 1.5)",
        class = "Intercept",
        resp = resp_name
      )
    ),

    do.call(
      brms::prior,
      list(
        "exponential(1)",
        class = "sigma",
        resp = resp_name
      )
    ),

    # facial disc model
    prior(normal(0, 1), class = "b", resp = "facialdisc"),
    prior(normal(0, 1.5), class = "Intercept", resp = "facialdisc")
  )

  ## DAG 1
  ## Song -> facial disc
  
  bf_song <- brms::bf(
    stats::as.formula(
      paste0(
        song_var,
        " | mi() ~ ",
        "mo(cover2) + ",
        "activity + ",
        "mi(niche_width_sc) + ",
        "(1 | gr(Species, cov = vcv.phylo))"
      )
    )
  )

  bf_disc <- brms::bf(
    stats::as.formula(
      paste0(
        "facial_disc ~ ",
        "mo(cover2) + ",
        "activity + ",
        "mi(niche_width_sc) + ",
        "mi(", song_var, ") + ",
        "(1 | gr(Species, cov = vcv.phylo))"
      )
    ),
    family = bernoulli()
  )

  fit_name <- paste0(
    "sem_song_to_disc_",
    gsub("_sc", "", song_var),
    "_mcc"
  )

  if (file.exists(paste0("./data/processed/fits/", fit_name, ".rds"))) {

    cat("Model already exists. Skipping:", fit_name, "\n")

  } else {

    cat("Running model:", fit_name, "\n")

    sem_song_to_disc <- brm(

      bf_latitude +
        bf_weight +
        bf_niche +
        bf_activity +
        bf_song +
        bf_disc +
        set_rescor(FALSE),

      data = dat,

      data2 = list(
        vcv.phylo = vcv.phylo
      ),

      prior = priors,

      chains = chains,
      cores = cores,
      iter = iters,

      backend = "cmdstanr",

      control = list(
        adapt_delta = 0.99,
        max_treedepth = 15
      ),

      file = paste0(
        "./data/processed/fits/",
        fit_name
      ),
      file_refit = "on_change"
    )

    ## Complete cases for LOO

    newdata_loo <- dat[
      complete.cases(
        dat$niche_width_sc,
        dat[[song_var]],
        dat$facial_disc,
        dat$activity
      ),
    ]

    ## Add LOO

    sem_song_to_disc <- add_criterion(
      sem_song_to_disc,
      criterion = "loo",
      newdata = newdata_loo
    )

    ## Save model

    saveRDS(
      sem_song_to_disc,
      file = paste0(
        "./data/processed/fits/",
        fit_name,
        ".rds"
      )
    )

    rm(sem_song_to_disc)

    gc()
  }

  ## ======================================================
  ## DAG 2
  ## Facial disc -> song
  ## ======================================================

  bf_disc <- bf(
    facial_disc ~
      mo(cover2) +
      activity +
      mi(niche_width_sc) +
      (1 | gr(Species, cov = vcv.phylo)),
    family = bernoulli()
  )

  bf_song <- brms::bf(
    stats::as.formula(
      paste0(
        song_var,
        " | mi() ~ ",
        "mo(cover2) + ",
        "activity + ",
        "mi(niche_width_sc) + ",
        "facial_disc + ",
        "(1 | gr(Species, cov = vcv.phylo))"
      )
    )
  )

  fit_name <- paste0(
    "sem_disc_to_song_",
    gsub("_sc", "", song_var),
    "_mcc"
  )

  if (file.exists(paste0("./data/processed/fits/", fit_name, ".rds"))) {

    cat("Model already exists. Skipping:", fit_name, "\n")

  } else {

    cat("Running model:", fit_name, "\n")

    sem_disc_to_song <- brm(

      bf_latitude +
        bf_weight +
        bf_niche +
        bf_activity +
        bf_disc +
        bf_song +
        set_rescor(FALSE),

      data = dat,

      data2 = list(
        vcv.phylo = vcv.phylo
      ),

      prior = priors,

      chains = chains,
      cores = cores,
      iter = iters,

      backend = "cmdstanr",

      control = list(
        adapt_delta = 0.99,
        max_treedepth = 15
      ),

      file = paste0(
        "./data/processed/fits/",
        fit_name
      ),
      file_refit = "on_change"
    )

    ## Add LOO

    sem_disc_to_song <- add_criterion(
      sem_disc_to_song,
      criterion = "loo",
      newdata = newdata_loo
    )

    ## Save model

    saveRDS(
      sem_disc_to_song,
      file = paste0(
        "./data/processed/fits/",
        fit_name,
        ".rds"
      )
    )

    rm(sem_disc_to_song)

    gc()
  }
}

3.1.1.1 Results

3.1.1.1.1 Comparing DAG model fit

The following table summarizes the ELPD differences comparing the predictive performance of the two DAGs (song -> disc vs disc -> song) for each of the song feature models (higher ELP mean better performance):

Code

features <- c("song_pc1_freq_sc", "song_pc2_diver_sc", "song_pc4_length_sc")

## Empty dataframe to store results

agg_loo <- data.frame()

## Loop over song variables

for (song_var in features) {

    ## File names

    fit1_file <- paste0("./data/processed/fits/", "sem_song_to_disc_",
        gsub("_sc", "", song_var), "_mcc.rds")

    fit2_file <- paste0("./data/processed/fits/", "sem_disc_to_song_",
        gsub("_sc", "", song_var), "_mcc.rds")

    ## Read models

    sem_song_to_disc <- readRDS(fit1_file)

    sem_disc_to_song <- readRDS(fit2_file)

    ## Compare models

    loo_diff <- loo::loo_compare(sem_disc_to_song, sem_song_to_disc)

    ## Convert to dataframe

    rows <- rownames(loo_diff)

    loo_diff <- as.data.frame(loo_diff)

    loo_diff$model <- rows

    loo_diff$feature <- song_var

    ## Store results

    agg_loo <- rbind(agg_loo, loo_diff)

    ## Cleanup

    rm(sem_song_to_disc, sem_disc_to_song, loo_diff)

    gc()
}

## Final table
agg_loo <- agg_loo[, c("feature", "model", "elpd_diff", "se_diff",
    "looic")]

x_kbl <- kableExtra::kbl(agg_loo, row.names = FALSE, escape = FALSE,
    format = "html", digits = 3)

x_kbl <- kableExtra::row_spec(kable_input = x_kbl, row = which(agg_loo$elpd_diff ==
    0), background = "#A0DFB980")

x_kbl

feature	model	elpd_diff	se_diff	looic
song_pc1_freq_sc	sem_song_to_disc	0.000	0.000	2340.481
song_pc1_freq_sc	sem_disc_to_song	-0.609	1.502	2341.698
song_pc2_diver_sc	sem_song_to_disc	0.000	0.000	2436.150
song_pc2_diver_sc	sem_disc_to_song	-2.582	3.188	2441.314
song_pc4_length_sc	sem_song_to_disc	0.000	0.000	2436.424
song_pc4_length_sc	sem_disc_to_song	-0.522	1.697	2437.468

Takeaway

Differences in predictive fit between competing DAGs were small relative to their associated uncertainty (()ELPD < SE in all comparisons), indicating that both causal structures were similarly consistent with the observed comparative data

3.1.1.1.2 Model fits

The following table summarizes the fixed effects estimates for each of the song feature models (song -> disc and disc -> song). Estimates are on the log-odds scale, and “significant” effects (95% credible intervals not overlapping zero) are highlighted.
Column names:
- dag: indicates the direction of the DAG (song -> disc or disc -> song)
- response: the response variable in the model
- predictor: the predictor variable in the model
- song_feature_model: the specific feature that was used in that particular model
- Estimate: the estimated coefficient for the predictor
- Est.Error: the standard error of the estimate
- l-95% CI: the lower bound of the 95% credible interval
- u-95% CI: the upper bound of the 95% credible interval
- pd: the probability of direction, a Bayesian measure quantifying the proportion of the posterior distribution supporting an effect in the same direction (positive or negative) as the posterior estimate. Values of pd close to 1 indicate strong directional support, whereas values near 0.5 indicate little evidence for a consistent effect direction

Color code:

strong directional support (pd > 0.95) orange
moderate support (0.90 < pd ≤ 0.95) yellow
weak suggestive support (0.80 < pd ≤ 0.90) pale yellow.

3.1.1.1.2.1 Ecological predictors

Code

## Empty dataframe
effects_df <- data.frame()

## Loop over models
fls <- list.files("./data/processed/fits/", pattern = "mcc\\.rds$",
    full.names = TRUE)

# keep fls only matching features

features <- c("peak_frequency", "song_duration", "umap_space_size",
    "frequency_range")


fls <- fls[!grepl(paste(features, collapse = "|"), fls)]


for (fl in fls) {

    ## Read model
    fit <- readRDS(fl)

    ## Extract fixed effects

    fe <- brms::fixef(fit, summary = TRUE, probs = c(0.025, 0.975))

    ## Convert to dataframe

    fe <- as.data.frame(fe)

    fe$model <- rownames(fe)

    rownames(fe) <- NULL

    fe <- fe[grep("Intercept", fe$model, invert = TRUE), ]

    fe$response <- sub("^(.*?)_(.*)$", "\\1", fe$model)
    fe$predictor <- sub("^(.*?)_(.*)$", "\\2", fe$model)

    ## Add feature label

    fe$song_feature_model <- gsub(".rds", "", sapply(strsplit(basename(fl),
        "_"), function(x) paste(x[5:7], collapse = "_")))

    fe$dag <- ifelse(grepl("song_to_disc", fl), "song_to_disc", "disc_to_song")

    ## --------------------------------------------------
    ## Approximate posterior probability of direction
    ## --------------------------------------------------

    fe$pd <- sapply(seq_len(nrow(fe)), function(i) {

        mu <- fe[["Estimate"]][i]

        ## approximate posterior sd from 95% CI
        sd_est <- abs(fe[["Q2.5"]][i] - fe[["Q97.5"]][i])/(2 * 1.96)

        pdir <- 1 - pnorm(0, mean = mu, sd = sd_est)

        max(pdir, 1 - pdir)
    })

    ## Keep relevant columns
    fe <- fe[, c("dag", "response", "predictor", "song_feature_model",
        "Estimate", "Est.Error", "Q2.5", "Q97.5", "pd")]

    ## Store
    effects_df <- rbind(effects_df, fe)

    ## Cleanup

    rm(fit, fe)

    gc()
}

effects_df <- effects_df[order(effects_df$predictor, effects_df$response,
    effects_df$dag), ]

effects_df$`response ~ predictor` <- paste(effects_df$response, "~",
    effects_df$predictor)


# show only non-song related effects
highlight(x = effects_df[!(grepl("song", effects_df$predictor) | grepl("song",
    effects_df$response)), c("dag", "song_feature_model", "response ~ predictor",
    "Estimate", "Q2.5", "Q97.5", "pd")])

	dag	song_feature_model	response ~ predictor	Estimate	Q2.5	Q97.5	pd
7	disc_to_song	num_elms_mcc	facialdisc ~ activity	-0.414	-1.755	0.986	0.723
71	disc_to_song	song_pc1_freq	facialdisc ~ activity	-0.452	-1.774	0.872	0.748
72	disc_to_song	song_pc2_diver	facialdisc ~ activity	-0.403	-1.743	0.927	0.723
73	disc_to_song	song_pc4_length	facialdisc ~ activity	-0.415	-1.739	1.052	0.720
74	disc_to_song	song_rate_mcc	facialdisc ~ activity	-0.415	-1.794	0.920	0.726
85	song_to_disc	num_elms_mcc	facialdisc ~ activity	-0.444	-1.776	0.938	0.739
86	song_to_disc	song_pc1_freq	facialdisc ~ activity	-0.416	-1.782	0.972	0.723
87	song_to_disc	song_pc2_diver	facialdisc ~ activity	-0.377	-1.799	0.939	0.705
88	song_to_disc	song_pc4_length	facialdisc ~ activity	-0.344	-1.807	1.071	0.680
89	song_to_disc	song_rate_mcc	facialdisc ~ activity	-0.424	-1.909	1.047	0.713
8	disc_to_song	num_elms_mcc	numelmssc ~ activity	0.068	-0.227	0.365	0.674
75	song_to_disc	num_elms_mcc	numelmssc ~ activity	0.056	-0.228	0.331	0.654
9	disc_to_song	num_elms_mcc	numelmssc ~ facial_disc	0.280	-0.150	0.711	0.899
10	disc_to_song	num_elms_mcc	nichewidthsc ~ milatitude_sc	0.013	-0.123	0.146	0.576
101	disc_to_song	song_pc1_freq	nichewidthsc ~ milatitude_sc	0.013	-0.113	0.148	0.576
102	disc_to_song	song_pc2_diver	nichewidthsc ~ milatitude_sc	0.012	-0.121	0.143	0.572
103	disc_to_song	song_pc4_length	nichewidthsc ~ milatitude_sc	0.009	-0.130	0.155	0.549
104	disc_to_song	song_rate_mcc	nichewidthsc ~ milatitude_sc	0.010	-0.128	0.141	0.558
95	song_to_disc	num_elms_mcc	nichewidthsc ~ milatitude_sc	0.011	-0.124	0.138	0.567
96	song_to_disc	song_pc1_freq	nichewidthsc ~ milatitude_sc	0.011	-0.113	0.140	0.568
97	song_to_disc	song_pc2_diver	nichewidthsc ~ milatitude_sc	0.014	-0.110	0.140	0.587
98	song_to_disc	song_pc4_length	nichewidthsc ~ milatitude_sc	0.009	-0.127	0.150	0.549
99	song_to_disc	song_rate_mcc	nichewidthsc ~ milatitude_sc	0.011	-0.131	0.153	0.559
13	disc_to_song	num_elms_mcc	activity ~ milog_weight_sc	0.021	-0.615	0.668	0.525
131	disc_to_song	song_pc1_freq	activity ~ milog_weight_sc	0.016	-0.619	0.604	0.520
132	disc_to_song	song_pc2_diver	activity ~ milog_weight_sc	0.013	-0.641	0.730	0.515
133	disc_to_song	song_pc4_length	activity ~ milog_weight_sc	0.006	-0.662	0.642	0.508
134	disc_to_song	song_rate_mcc	activity ~ milog_weight_sc	0.013	-0.627	0.662	0.515
125	song_to_disc	num_elms_mcc	activity ~ milog_weight_sc	0.028	-0.641	0.688	0.533
126	song_to_disc	song_pc1_freq	activity ~ milog_weight_sc	0.023	-0.649	0.647	0.528
127	song_to_disc	song_pc2_diver	activity ~ milog_weight_sc	0.008	-0.688	0.633	0.509
128	song_to_disc	song_pc4_length	activity ~ milog_weight_sc	0.005	-0.649	0.639	0.506
129	song_to_disc	song_rate_mcc	activity ~ milog_weight_sc	0.021	-0.644	0.676	0.525
11	disc_to_song	num_elms_mcc	nichewidthsc ~ milog_weight_sc	0.313	0.121	0.517	0.999
111	disc_to_song	song_pc1_freq	nichewidthsc ~ milog_weight_sc	0.319	0.111	0.513	0.999
112	disc_to_song	song_pc2_diver	nichewidthsc ~ milog_weight_sc	0.309	0.117	0.499	0.999
113	disc_to_song	song_pc4_length	nichewidthsc ~ milog_weight_sc	0.314	0.119	0.523	0.999
114	disc_to_song	song_rate_mcc	nichewidthsc ~ milog_weight_sc	0.318	0.121	0.519	0.999
105	song_to_disc	num_elms_mcc	nichewidthsc ~ milog_weight_sc	0.313	0.118	0.514	0.999
106	song_to_disc	song_pc1_freq	nichewidthsc ~ milog_weight_sc	0.316	0.126	0.507	0.999
107	song_to_disc	song_pc2_diver	nichewidthsc ~ milog_weight_sc	0.309	0.112	0.522	0.998
108	song_to_disc	song_pc4_length	nichewidthsc ~ milog_weight_sc	0.314	0.118	0.529	0.999
109	song_to_disc	song_rate_mcc	nichewidthsc ~ milog_weight_sc	0.312	0.123	0.506	0.999
15	disc_to_song	num_elms_mcc	facialdisc ~ miniche_width_sc	0.359	-0.518	1.335	0.777
151	disc_to_song	song_pc1_freq	facialdisc ~ miniche_width_sc	0.349	-0.585	1.270	0.770
152	disc_to_song	song_pc2_diver	facialdisc ~ miniche_width_sc	0.347	-0.562	1.482	0.747
153	disc_to_song	song_pc4_length	facialdisc ~ miniche_width_sc	0.331	-0.580	1.302	0.755
154	disc_to_song	song_rate_mcc	facialdisc ~ miniche_width_sc	0.353	-0.533	1.330	0.771
165	song_to_disc	num_elms_mcc	facialdisc ~ miniche_width_sc	0.381	-0.587	1.443	0.769
166	song_to_disc	song_pc1_freq	facialdisc ~ miniche_width_sc	0.384	-0.525	1.375	0.786
167	song_to_disc	song_pc2_diver	facialdisc ~ miniche_width_sc	0.358	-0.562	1.355	0.768
168	song_to_disc	song_pc4_length	facialdisc ~ miniche_width_sc	0.366	-0.599	1.366	0.767
169	song_to_disc	song_rate_mcc	facialdisc ~ miniche_width_sc	0.407	-0.494	1.413	0.799
17	disc_to_song	num_elms_mcc	numelmssc ~ miniche_width_sc	0.023	-0.118	0.171	0.623
145	song_to_disc	num_elms_mcc	numelmssc ~ miniche_width_sc	0.027	-0.109	0.165	0.652
175	song_to_disc	num_elms_mcc	facialdisc ~ minum_elms_sc	0.560	-0.305	1.550	0.882
12	disc_to_song	num_elms_mcc	activity ~ mocover2	0.784	0.141	1.518	0.987
121	disc_to_song	song_pc1_freq	activity ~ mocover2	0.780	0.123	1.511	0.986
122	disc_to_song	song_pc2_diver	activity ~ mocover2	0.770	0.183	1.508	0.989
123	disc_to_song	song_pc4_length	activity ~ mocover2	0.782	0.185	1.569	0.987
124	disc_to_song	song_rate_mcc	activity ~ mocover2	0.793	0.155	1.585	0.985
115	song_to_disc	num_elms_mcc	activity ~ mocover2	0.779	0.167	1.558	0.986
116	song_to_disc	song_pc1_freq	activity ~ mocover2	0.794	0.178	1.569	0.987
117	song_to_disc	song_pc2_diver	activity ~ mocover2	0.780	0.183	1.491	0.990
118	song_to_disc	song_pc4_length	activity ~ mocover2	0.793	0.192	1.556	0.989
119	song_to_disc	song_rate_mcc	activity ~ mocover2	0.792	0.175	1.533	0.989
14	disc_to_song	num_elms_mcc	facialdisc ~ mocover2	1.241	0.271	2.370	0.990
141	disc_to_song	song_pc1_freq	facialdisc ~ mocover2	1.253	0.273	2.330	0.992
142	disc_to_song	song_pc2_diver	facialdisc ~ mocover2	1.241	0.335	2.287	0.994
143	disc_to_song	song_pc4_length	facialdisc ~ mocover2	1.270	0.185	2.560	0.982
144	disc_to_song	song_rate_mcc	facialdisc ~ mocover2	1.250	0.244	2.397	0.989
155	song_to_disc	num_elms_mcc	facialdisc ~ mocover2	1.296	0.221	2.492	0.987
156	song_to_disc	song_pc1_freq	facialdisc ~ mocover2	1.303	0.320	2.501	0.990
157	song_to_disc	song_pc2_diver	facialdisc ~ mocover2	1.251	0.276	2.399	0.990
158	song_to_disc	song_pc4_length	facialdisc ~ mocover2	1.211	0.140	2.336	0.985
159	song_to_disc	song_rate_mcc	facialdisc ~ mocover2	1.232	0.281	2.350	0.990
16	disc_to_song	num_elms_mcc	numelmssc ~ mocover2	-0.063	-0.292	0.152	0.710
135	song_to_disc	num_elms_mcc	numelmssc ~ mocover2	-0.028	-0.274	0.193	0.594

3.1.1.1.2.2 Song features

Code

highlight(x = effects_df[grepl("pc", effects_df$predictor) | grepl("pc",
    effects_df$response), c("dag", "song_feature_model", "response ~ predictor",
    "Estimate", "Q2.5", "Q97.5", "pd")])

	dag	song_feature_model	response ~ predictor	Estimate	Q2.5	Q97.5	pd
81	disc_to_song	song_pc1_freq	songpc1freqsc ~ activity	0.096	-0.159	0.353	0.769
76	song_to_disc	song_pc1_freq	songpc1freqsc ~ activity	0.073	-0.191	0.323	0.712
82	disc_to_song	song_pc2_diver	songpc2diversc ~ activity	-0.117	-0.437	0.194	0.766
77	song_to_disc	song_pc2_diver	songpc2diversc ~ activity	-0.174	-0.469	0.127	0.874
83	disc_to_song	song_pc4_length	songpc4lengthsc ~ activity	0.065	-0.242	0.359	0.663
78	song_to_disc	song_pc4_length	songpc4lengthsc ~ activity	0.129	-0.202	0.470	0.775
91	disc_to_song	song_pc1_freq	songpc1freqsc ~ facial_disc	0.236	-0.131	0.609	0.895
92	disc_to_song	song_pc2_diver	songpc2diversc ~ facial_disc	0.348	-0.068	0.764	0.950
93	disc_to_song	song_pc4_length	songpc4lengthsc ~ facial_disc	-0.309	-0.692	0.087	0.940
171	disc_to_song	song_pc1_freq	songpc1freqsc ~ miniche_width_sc	0.006	-0.129	0.139	0.536
146	song_to_disc	song_pc1_freq	songpc1freqsc ~ miniche_width_sc	0.005	-0.120	0.142	0.529
172	disc_to_song	song_pc2_diver	songpc2diversc ~ miniche_width_sc	0.023	-0.127	0.188	0.613
147	song_to_disc	song_pc2_diver	songpc2diversc ~ miniche_width_sc	0.036	-0.117	0.187	0.677
173	disc_to_song	song_pc4_length	songpc4lengthsc ~ miniche_width_sc	0.015	-0.139	0.171	0.575
148	song_to_disc	song_pc4_length	songpc4lengthsc ~ miniche_width_sc	0.019	-0.146	0.196	0.584
176	song_to_disc	song_pc1_freq	facialdisc ~ misong_pc1_freq_sc	0.417	-0.564	1.470	0.789
177	song_to_disc	song_pc2_diver	facialdisc ~ misong_pc2_diver_sc	0.889	-0.078	1.887	0.962
178	song_to_disc	song_pc4_length	facialdisc ~ misong_pc4_length_sc	-0.725	-1.854	0.183	0.918
161	disc_to_song	song_pc1_freq	songpc1freqsc ~ mocover2	0.032	-0.165	0.237	0.623
136	song_to_disc	song_pc1_freq	songpc1freqsc ~ mocover2	0.058	-0.116	0.244	0.735
162	disc_to_song	song_pc2_diver	songpc2diversc ~ mocover2	-0.078	-0.326	0.167	0.731
137	song_to_disc	song_pc2_diver	songpc2diversc ~ mocover2	-0.033	-0.265	0.181	0.613
163	disc_to_song	song_pc4_length	songpc4lengthsc ~ mocover2	-0.118	-0.359	0.116	0.834
138	song_to_disc	song_pc4_length	songpc4lengthsc ~ mocover2	-0.166	-0.382	0.044	0.936

Takeaways

Across both DAGs, habitat structure and activity rhythm showed consistent positive associations with facial disc morphology, while body size was positively associated with niche width
Associations between facial disc morphology and song structure were generally modest but directionally consistent across alternative causal models, with some relationships showing moderate directional support (e.g., between facial disc morphology and song diversity)
However, leave-one-out cross-validation provided little evidence favoring either the “song → facial disc” or “facial disc → song” evolutionary hypothesis, suggesting that the available comparative data do not clearly distinguish between these alternative causal scenarios

3.1.1.1.3 Individual model fits

3.1.1.1.3.1 Song → facial disc hypothesis

Code

## Find all MCC models

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

## Loop over models

for (i in seq_along(model_files)) {

    f <- model_files[i]

    ## Read model
    mod <- readRDS(f)

    cat(paste0("\n##### ", gsub(".rds", "", sapply(strsplit(basename(f),
        "_"), function(x) paste(x[5:7], collapse = "_"))), "\n\n"))

    ## Get summary
    extended_summary(mod, highlight = TRUE, remove.intercepts = TRUE,
        print.name = FALSE, beta.prefix = c("^b_", "^bsp_"))
}

3.1.1.1.4 frequency_range_mcc

	formula	family	priors	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	list(latitudesc = list(formula = latitude_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “latitudesc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), logweightsc = list(formula = log_weight_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “logweightsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), nichewidthsc = list(formula = niche_width_sc \| mi() ~ mi(latitude_sc) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “nichewidthsc”, family = list( family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), activity = list(formula = activity ~ mo(cover2) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “activity”, mecor = TRUE), frequencyrangesc = list(formula = frequency_range_sc \| mi() ~ mo(cover2) + activity + mi(niche_width_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “frequencyrangesc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), facialdisc = list(formula = facial_disc ~ mo(cover2) + activity + mi(niche_width_sc) + mi(frequency_range_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list( family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “facialdisc”, mecor = TRUE))	()	b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) Intercept-student_t(3, -0.2, 2.5) Intercept-student_t(3, -0.2, 2.5) Intercept-normal(0, 1.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sigma-exponential(1) sigma-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-exponential(1) simo-dirichlet(1) simo-dirichlet(1) simo-dirichlet(1)	4000	1	1	2000	0 (0%)	0	2231.872	1179.634	680754594

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_latitudesc_Intercept	-0.001	-0.144	0.145	1.001	4167.945	1374.829
b_logweightsc_Intercept	-0.001	-0.127	0.127	1.004	4707.042	1517.103
b_nichewidthsc_Intercept	-0.097	-0.707	0.450	1	2442.106	1415.618
b_activity_Intercept	-0.840	-2.796	0.952	1.002	3241.232	1609.371
b_frequencyrangesc_Intercept	0.699	-0.146	1.539	1	2231.872	1494.963
b_facialdisc_Intercept	1.028	-1.754	3.711	1.003	3145.118	1448.405
b_frequencyrangesc_activity	-0.065	-0.368	0.226	1	4683.566	1621.899
b_facialdisc_activity	-0.379	-1.759	1.019	1	3156.037	1426.992
bsp_nichewidthsc_milatitude_sc	0.013	-0.120	0.150	1.001	3894.530	1179.634
bsp_nichewidthsc_milog_weight_sc	0.313	0.126	0.514	1	2731.004	1390.159
bsp_activity_mocover2	0.793	0.162	1.519	1	3155.123	1561.107
bsp_activity_milog_weight_sc	0.017	-0.642	0.678	1	3726.135	1869.842
bsp_frequencyrangesc_mocover2	-0.023	-0.223	0.190	1	2985.806	1706.962
bsp_frequencyrangesc_miniche_width_sc	0.007	-0.142	0.156	1	3554.277	1293.184
bsp_facialdisc_mocover2	1.284	0.281	2.435	1	2505.874	1818.608
bsp_facialdisc_miniche_width_sc	0.363	-0.584	1.404	1	2916.467	1459.642
bsp_facialdisc_mifrequency_range_sc	0.370	-0.572	1.434	1.001	2383.793	1608.132

##### num_elms_mcc

	formula	family	priors	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	list(latitudesc = list(formula = latitude_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “latitudesc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), logweightsc = list(formula = log_weight_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “logweightsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), nichewidthsc = list(formula = niche_width_sc \| mi() ~ mi(latitude_sc) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “nichewidthsc”, family = list( family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), activity = list(formula = activity ~ mo(cover2) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “activity”, mecor = TRUE), numelmssc = list(formula = num_elms_sc \| mi() ~ mo(cover2) + activity + mi(niche_width_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “numelmssc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), facialdisc = list(formula = facial_disc ~ mo(cover2) + activity + mi(niche_width_sc) + mi(num_elms_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “facialdisc”, mecor = TRUE))	()	b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) Intercept-student_t(3, -0.2, 2.5) Intercept-student_t(3, -0.2, 2.5) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-exponential(1) sigma-exponential(1) simo-dirichlet(1) simo-dirichlet(1) simo-dirichlet(1)	4000	1	1	2000	0 (0%)	0	458.297	939.465	1183312160

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_latitudesc_Intercept	-0.002	-0.138	0.136	1	4918.989	1227.915
b_logweightsc_Intercept	0.000	-0.139	0.138	1	4796.173	1306.761
b_nichewidthsc_Intercept	-0.068	-0.659	0.523	1.001	1567.930	1054.309
b_activity_Intercept	-0.831	-2.624	0.881	1.003	1993.005	1567.437
b_numelmssc_Intercept	-0.080	-1.380	1.153	1.003	782.569	939.465
b_facialdisc_Intercept	1.124	-1.658	3.810	1	2149.980	1267.235
b_numelmssc_activity	0.056	-0.228	0.331	1	2081.116	1657.457
b_facialdisc_activity	-0.444	-1.776	0.938	1	2983.116	1574.322
bsp_nichewidthsc_milatitude_sc	0.011	-0.124	0.138	1	2831.825	1370.249
bsp_nichewidthsc_milog_weight_sc	0.313	0.118	0.514	1	3041.907	1517.847
bsp_activity_mocover2	0.779	0.167	1.558	1.001	2417.322	1639.373
bsp_activity_milog_weight_sc	0.028	-0.641	0.688	1	2168.640	1441.231
bsp_numelmssc_mocover2	-0.028	-0.274	0.193	1.002	458.297	1300.786
bsp_numelmssc_miniche_width_sc	0.027	-0.109	0.165	1.002	1570.626	1462.690
bsp_facialdisc_mocover2	1.296	0.221	2.492	1	1698.620	1591.066
bsp_facialdisc_miniche_width_sc	0.381	-0.587	1.443	1	1929.746	1519.758
bsp_facialdisc_minum_elms_sc	0.560	-0.305	1.550	1	2274.523	1265.840

##### peak_frequency_mcc

	formula	family	priors	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	list(latitudesc = list(formula = latitude_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “latitudesc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), logweightsc = list(formula = log_weight_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “logweightsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), nichewidthsc = list(formula = niche_width_sc \| mi() ~ mi(latitude_sc) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “nichewidthsc”, family = list( family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), activity = list(formula = activity ~ mo(cover2) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “activity”, mecor = TRUE), peakfrequencysc = list(formula = peak_frequency_sc \| mi() ~ mo(cover2) + activity + mi(niche_width_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “peakfrequencysc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), facialdisc = list(formula = facial_disc ~ mo(cover2) + activity + mi(niche_width_sc) + mi(peak_frequency_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list( family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “facialdisc”, mecor = TRUE))	()	b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) Intercept-student_t(3, -0.2, 2.5) Intercept-student_t(3, -0.2, 2.5) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-exponential(1) sigma-exponential(1) simo-dirichlet(1) simo-dirichlet(1) simo-dirichlet(1)	4000	1	1	2000	0 (0%)	0	762.976	1044.478	1834181007

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_latitudesc_Intercept	-0.001	-0.148	0.146	1	4370.793	1152.909
b_logweightsc_Intercept	0.001	-0.136	0.139	1	3381.870	1168.300
b_nichewidthsc_Intercept	-0.079	-0.668	0.483	1	1244.731	1291.823
b_activity_Intercept	-0.803	-2.683	0.998	1	1581.067	1515.897
b_peakfrequencysc_Intercept	0.703	-0.061	1.510	1	762.976	1044.478
b_facialdisc_Intercept	1.044	-1.636	3.846	1.001	2011.710	1496.895
b_peakfrequencysc_activity	-0.054	-0.334	0.228	1	2304.708	1462.674
b_facialdisc_activity	-0.394	-1.765	0.994	1	2452.438	1626.928
bsp_nichewidthsc_milatitude_sc	0.013	-0.120	0.154	1.005	2745.067	1385.627
bsp_nichewidthsc_milog_weight_sc	0.312	0.113	0.518	1	2065.652	1474.205
bsp_activity_mocover2	0.779	0.148	1.567	1.002	1678.053	1284.244
bsp_activity_milog_weight_sc	0.015	-0.598	0.645	1.002	1943.303	1486.396
bsp_peakfrequencysc_mocover2	0.000	-0.195	0.204	1.002	1828.947	1526.451
bsp_peakfrequencysc_miniche_width_sc	0.009	-0.128	0.137	1	2372.776	1650.427
bsp_facialdisc_mocover2	1.293	0.221	2.504	1	1143.134	1424.940
bsp_facialdisc_miniche_width_sc	0.365	-0.517	1.371	1	1993.072	1662.957
bsp_facialdisc_mipeak_frequency_sc	0.054	-0.878	1.083	1.001	1781.900	1308.604

##### song_duration_mcc

	formula	family	priors	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	list(latitudesc = list(formula = latitude_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “latitudesc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), logweightsc = list(formula = log_weight_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “logweightsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), nichewidthsc = list(formula = niche_width_sc \| mi() ~ mi(latitude_sc) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “nichewidthsc”, family = list( family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), activity = list(formula = activity ~ mo(cover2) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “activity”, mecor = TRUE), songdurationsc = list(formula = song_duration_sc \| mi() ~ mo(cover2) + activity + mi(niche_width_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “songdurationsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), facialdisc = list(formula = facial_disc ~ mo(cover2) + activity + mi(niche_width_sc) + mi(song_duration_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list( family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “facialdisc”, mecor = TRUE))	()	b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) Intercept-student_t(3, -0.2, 2.5) Intercept-student_t(3, -0.2, 2.5) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-exponential(1) sigma-exponential(1) simo-dirichlet(1) simo-dirichlet(1) simo-dirichlet(1)	4000	1	1	2000	0 (0%)	0	1083.251	800.878	223180684

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_latitudesc_Intercept	0.000	-0.141	0.147	1.005	3237.610	1227.281
b_logweightsc_Intercept	-0.002	-0.153	0.140	1.001	3289.139	1243.804
b_nichewidthsc_Intercept	-0.087	-0.698	0.503	1	1083.251	800.878
b_activity_Intercept	-0.853	-2.590	0.998	1	1398.798	1253.973
b_songdurationsc_Intercept	-0.044	-0.449	0.325	1	1757.019	1324.395
b_facialdisc_Intercept	1.199	-1.548	3.959	1	2536.668	1587.976
b_songdurationsc_activity	0.024	-0.293	0.342	1	2997.485	1742.219
b_facialdisc_activity	-0.451	-1.840	0.959	1	2511.911	1479.973
bsp_nichewidthsc_milatitude_sc	0.013	-0.116	0.149	1	2168.103	1395.604
bsp_nichewidthsc_milog_weight_sc	0.316	0.125	0.509	1	1657.317	1516.752
bsp_activity_mocover2	0.798	0.165	1.566	1.001	1810.354	1542.754
bsp_activity_milog_weight_sc	0.025	-0.617	0.678	1	1931.921	1457.503
bsp_songdurationsc_mocover2	0.005	-0.259	0.234	1	1446.308	1413.504
bsp_songdurationsc_miniche_width_sc	-0.029	-0.194	0.133	1	2385.773	1363.323
bsp_facialdisc_mocover2	1.215	0.274	2.283	1.001	1629.122	1502.989
bsp_facialdisc_miniche_width_sc	0.393	-0.450	1.335	1.001	2084.108	1496.586
bsp_facialdisc_misong_duration_sc	1.107	-0.076	2.523	1.002	2362.838	1513.892

##### song_pc1_freq

	formula	family	priors	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	list(latitudesc = list(formula = latitude_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “latitudesc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), logweightsc = list(formula = log_weight_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “logweightsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), nichewidthsc = list(formula = niche_width_sc \| mi() ~ mi(latitude_sc) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “nichewidthsc”, family = list( family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), activity = list(formula = activity ~ mo(cover2) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “activity”, mecor = TRUE), songpc1freqsc = list(formula = song_pc1_freq_sc \| mi() ~ mo(cover2) + activity + mi(niche_width_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “songpc1freqsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), facialdisc = list(formula = facial_disc ~ mo(cover2) + activity + mi(niche_width_sc) + mi(song_pc1_freq_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list( family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “facialdisc”, mecor = TRUE))	()	b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) Intercept-student_t(3, -0.2, 2.5) Intercept-student_t(3, -0.2, 2.5) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-exponential(1) sigma-exponential(1) simo-dirichlet(1) simo-dirichlet(1) simo-dirichlet(1)	2000	1	1	1000	0 (0%)	0	345.647	363.752	1235109842

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_latitudesc_Intercept	0.001	-0.143	0.144	0.999	1140.688	767.471
b_logweightsc_Intercept	0.000	-0.129	0.140	0.999	1739.710	789.445
b_nichewidthsc_Intercept	-0.095	-0.701	0.408	1.003	361.063	363.752
b_activity_Intercept	-0.737	-2.574	1.097	1.004	345.647	567.929
b_songpc1freqsc_Intercept	-0.922	-1.827	0.016	1.001	351.412	456.495
b_facialdisc_Intercept	1.167	-1.575	4.183	0.999	500.920	549.564
b_songpc1freqsc_activity	0.073	-0.191	0.323	1.001	577.442	658.549
b_facialdisc_activity	-0.416	-1.782	0.972	1.003	715.899	663.223
bsp_nichewidthsc_milatitude_sc	0.011	-0.113	0.140	0.999	1011.699	814.361
bsp_nichewidthsc_milog_weight_sc	0.316	0.126	0.507	1.01	529.133	761.165
bsp_activity_mocover2	0.794	0.178	1.569	1.002	634.769	530.470
bsp_activity_milog_weight_sc	0.023	-0.649	0.647	1.001	557.045	690.066
bsp_songpc1freqsc_mocover2	0.058	-0.116	0.244	1.005	533.342	655.714
bsp_songpc1freqsc_miniche_width_sc	0.005	-0.120	0.142	0.999	570.229	632.860
bsp_facialdisc_mocover2	1.303	0.320	2.501	1.001	468.710	542.948
bsp_facialdisc_miniche_width_sc	0.384	-0.525	1.375	1.003	479.737	531.116
bsp_facialdisc_misong_pc1_freq_sc	0.417	-0.564	1.470	1.001	601.166	617.039

##### song_pc2_diver

	formula	family	priors	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	list(latitudesc = list(formula = latitude_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “latitudesc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), logweightsc = list(formula = log_weight_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “logweightsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), nichewidthsc = list(formula = niche_width_sc \| mi() ~ mi(latitude_sc) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “nichewidthsc”, family = list( family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), activity = list(formula = activity ~ mo(cover2) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “activity”, mecor = TRUE), songpc2diversc = list(formula = song_pc2_diver_sc \| mi() ~ mo(cover2) + activity + mi(niche_width_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “songpc2diversc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), facialdisc = list(formula = facial_disc ~ mo(cover2) + activity + mi(niche_width_sc) + mi(song_pc2_diver_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list( family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “facialdisc”, mecor = TRUE))	()	b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) Intercept-student_t(3, -0.2, 2.5) Intercept-student_t(3, -0.2, 2.5) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-exponential(1) sigma-exponential(1) simo-dirichlet(1) simo-dirichlet(1) simo-dirichlet(1)	2000	1	1	1000	0 (0%)	0	829.254	371.223	739154812

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_latitudesc_Intercept	0.001	-0.136	0.144	0.999	1777.512	694.868
b_logweightsc_Intercept	-0.001	-0.138	0.129	1.003	1871.207	781.341
b_nichewidthsc_Intercept	-0.088	-0.719	0.482	0.999	1013.577	760.251
b_activity_Intercept	-0.815	-2.700	1.049	1.001	1140.865	682.637
b_songpc2diversc_Intercept	0.365	-0.319	1.060	1.013	923.418	371.223
b_facialdisc_Intercept	1.078	-1.916	3.683	1	1648.315	819.515
b_songpc2diversc_activity	-0.174	-0.469	0.127	1.002	1410.238	792.484
b_facialdisc_activity	-0.377	-1.799	0.939	0.999	1226.388	864.779
bsp_nichewidthsc_milatitude_sc	0.014	-0.110	0.140	1.002	1776.165	983.283
bsp_nichewidthsc_milog_weight_sc	0.309	0.112	0.522	1.005	1308.989	734.754
bsp_activity_mocover2	0.780	0.183	1.491	0.999	922.369	601.526
bsp_activity_milog_weight_sc	0.008	-0.688	0.633	0.999	1005.653	739.535
bsp_songpc2diversc_mocover2	-0.033	-0.265	0.181	1.004	1276.545	865.819
bsp_songpc2diversc_miniche_width_sc	0.036	-0.117	0.187	1.001	829.254	506.599
bsp_facialdisc_mocover2	1.251	0.276	2.399	0.999	1072.266	870.904
bsp_facialdisc_miniche_width_sc	0.358	-0.562	1.355	0.999	1021.361	584.738
bsp_facialdisc_misong_pc2_diver_sc	0.889	-0.078	1.887	1	1154.561	843.070

##### song_pc4_length

	formula	family	priors	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	list(latitudesc = list(formula = latitude_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “latitudesc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), logweightsc = list(formula = log_weight_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “logweightsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), nichewidthsc = list(formula = niche_width_sc \| mi() ~ mi(latitude_sc) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “nichewidthsc”, family = list( family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), activity = list(formula = activity ~ mo(cover2) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “activity”, mecor = TRUE), songpc4lengthsc = list(formula = song_pc4_length_sc \| mi() ~ mo(cover2) + activity + mi(niche_width_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “songpc4lengthsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), facialdisc = list(formula = facial_disc ~ mo(cover2) + activity + mi(niche_width_sc) + mi(song_pc4_length_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list( family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “facialdisc”, mecor = TRUE))	()	b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) Intercept-student_t(3, -0.2, 2.5) Intercept-student_t(3, -0.2, 2.5) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-exponential(1) sigma-exponential(1) simo-dirichlet(1) simo-dirichlet(1) simo-dirichlet(1)	2000	1	1	1000	0 (0%)	0	372.165	388.158	1511329040

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_latitudesc_Intercept	0.001	-0.132	0.137	1	1055.975	670.095
b_logweightsc_Intercept	-0.002	-0.147	0.124	1.002	1085.364	691.745
b_nichewidthsc_Intercept	-0.068	-0.731	0.524	1	372.165	413.707
b_activity_Intercept	-0.840	-2.670	0.756	1.005	449.588	625.115
b_songpc4lengthsc_Intercept	-0.043	-0.612	0.479	1.002	415.559	388.158
b_facialdisc_Intercept	1.044	-1.912	3.930	1.003	765.330	660.122
b_songpc4lengthsc_activity	0.129	-0.202	0.470	1	657.534	563.539
b_facialdisc_activity	-0.344	-1.807	1.071	0.999	900.308	781.060
bsp_nichewidthsc_milatitude_sc	0.009	-0.127	0.150	1.008	709.036	646.618
bsp_nichewidthsc_milog_weight_sc	0.314	0.118	0.529	1.003	575.835	518.145
bsp_activity_mocover2	0.793	0.192	1.556	1	574.269	591.348
bsp_activity_milog_weight_sc	0.005	-0.649	0.639	1.004	667.282	464.163
bsp_songpc4lengthsc_mocover2	-0.166	-0.382	0.044	1.002	616.826	550.720
bsp_songpc4lengthsc_miniche_width_sc	0.019	-0.146	0.196	1.006	655.297	601.685
bsp_facialdisc_mocover2	1.211	0.140	2.336	1	474.822	698.963
bsp_facialdisc_miniche_width_sc	0.366	-0.599	1.366	0.999	750.273	749.046
bsp_facialdisc_misong_pc4_length_sc	-0.725	-1.854	0.183	1.002	725.810	664.496

##### song_rate_mcc

	formula	family	priors	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	list(latitudesc = list(formula = latitude_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “latitudesc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), logweightsc = list(formula = log_weight_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “logweightsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), nichewidthsc = list(formula = niche_width_sc \| mi() ~ mi(latitude_sc) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “nichewidthsc”, family = list( family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), activity = list(formula = activity ~ mo(cover2) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “activity”, mecor = TRUE), songratesc = list(formula = song_rate_sc \| mi() ~ mo(cover2) + activity + mi(niche_width_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “songratesc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), facialdisc = list(formula = facial_disc ~ mo(cover2) + activity + mi(niche_width_sc) + mi(song_rate_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “facialdisc”, mecor = TRUE))	()	b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) Intercept-student_t(3, -0.2, 2.5) Intercept-student_t(3, -0.2, 2.5) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-exponential(1) sigma-exponential(1) simo-dirichlet(1) simo-dirichlet(1) simo-dirichlet(1)	4000	1	1	2000	0 (0%)	0	1589.897	1039.791	1306371089

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_latitudesc_Intercept	0.000	-0.134	0.134	1.001	3906.271	1520.078
b_logweightsc_Intercept	-0.001	-0.140	0.141	1	3157.170	1263.559
b_nichewidthsc_Intercept	-0.092	-0.719	0.440	1.001	1589.897	1326.483
b_activity_Intercept	-0.827	-2.640	0.946	1	2223.801	1265.333
b_songratesc_Intercept	-0.155	-0.792	0.515	1.002	2208.444	1512.444
b_facialdisc_Intercept	1.144	-1.450	3.746	1	2720.794	1371.040
b_songratesc_activity	-0.119	-0.436	0.198	1.004	3894.048	1449.124
b_facialdisc_activity	-0.424	-1.909	1.047	1.001	2694.916	1512.727
bsp_nichewidthsc_milatitude_sc	0.011	-0.131	0.153	1	3594.797	1039.791
bsp_nichewidthsc_milog_weight_sc	0.312	0.123	0.506	1	2572.175	1544.435
bsp_activity_mocover2	0.792	0.175	1.533	1.001	2211.314	1561.588
bsp_activity_milog_weight_sc	0.021	-0.644	0.676	1	2552.811	1269.559
bsp_songratesc_mocover2	-0.024	-0.247	0.182	1	2497.532	1513.840
bsp_songratesc_miniche_width_sc	-0.060	-0.213	0.090	1	2956.185	1401.253
bsp_facialdisc_mocover2	1.232	0.281	2.350	1.002	1983.591	1551.817
bsp_facialdisc_miniche_width_sc	0.407	-0.494	1.413	1.002	2512.305	1494.014
bsp_facialdisc_misong_rate_sc	0.499	-0.358	1.456	1	2043.143	1618.695

##### umap_space_size

	formula	family	priors	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	list(latitudesc = list(formula = latitude_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “latitudesc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), logweightsc = list(formula = log_weight_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “logweightsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), nichewidthsc = list(formula = niche_width_sc \| mi() ~ mi(latitude_sc) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “nichewidthsc”, family = list( family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), activity = list(formula = activity ~ mo(cover2) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “activity”, mecor = TRUE), umapspacesizesc = list(formula = umap_space_size_sc \| mi() ~ mo(cover2) + activity + mi(niche_width_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “umapspacesizesc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), facialdisc = list(formula = facial_disc ~ mo(cover2) + activity + mi(niche_width_sc) + mi(umap_space_size_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list( family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “facialdisc”, mecor = TRUE))	()	b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) Intercept-student_t(3, -0.2, 2.5) Intercept-student_t(3, -0.2, 2.5) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-exponential(1) sigma-exponential(1) simo-dirichlet(1) simo-dirichlet(1) simo-dirichlet(1)	4000	1	1	2000	0 (0%)	0	729.034	687.181	1204202469

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_latitudesc_Intercept	0.001	-0.136	0.137	1	2333.965	1547.484
b_logweightsc_Intercept	0.001	-0.139	0.141	1	1703.993	1250.598
b_nichewidthsc_Intercept	-0.091	-0.739	0.452	1.001	729.034	687.181
b_activity_Intercept	-0.829	-2.661	0.880	1.001	898.562	1135.474
b_umapspacesizesc_Intercept	0.065	-0.379	0.470	1.003	791.826	1211.632
b_facialdisc_Intercept	1.091	-1.667	3.935	1.002	1046.015	1363.899
b_umapspacesizesc_activity	-0.123	-0.451	0.199	1	1786.054	1443.755
b_facialdisc_activity	-0.374	-1.812	1.111	1.002	1379.600	1130.938
bsp_nichewidthsc_milatitude_sc	0.010	-0.126	0.149	1.001	1535.857	1306.784
bsp_nichewidthsc_milog_weight_sc	0.314	0.121	0.516	1	938.631	828.198
bsp_activity_mocover2	0.777	0.141	1.541	1.002	1067.793	854.111
bsp_activity_milog_weight_sc	0.008	-0.637	0.655	1	778.924	959.670
bsp_umapspacesizesc_mocover2	-0.107	-0.337	0.114	1	1097.516	1261.331
bsp_umapspacesizesc_miniche_width_sc	0.004	-0.152	0.168	1	1293.958	1496.773
bsp_facialdisc_mocover2	1.294	0.270	2.448	1	1132.801	1297.989
bsp_facialdisc_miniche_width_sc	0.352	-0.526	1.318	1	1236.842	1114.844
bsp_facialdisc_miumap_space_size_sc	0.412	-0.486	1.606	1	1236.755	1404.867

3.1.1.1.4.1 Facial disc → song hypothesis

Code

## Find all MCC models

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

## Loop over models

for (i in seq_along(model_files)) {

    f <- model_files[i]

    ## Read model
    mod <- readRDS(f)

    cat(paste0("\n##### ", gsub(".rds", "", sapply(strsplit(basename(f),
        "_"), function(x) paste(x[5:7], collapse = "_"))), "\n\n"))

    ## Get summary
    extended_summary(mod, highlight = TRUE, remove.intercepts = TRUE,
        print.name = FALSE, beta.prefix = c("^b_", "^bsp_"))
}

3.1.1.1.5 frequency_range_mcc

	formula	family	priors	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	list(latitudesc = list(formula = latitude_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “latitudesc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), logweightsc = list(formula = log_weight_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “logweightsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), nichewidthsc = list(formula = niche_width_sc \| mi() ~ mi(latitude_sc) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “nichewidthsc”, family = list( family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), activity = list(formula = activity ~ mo(cover2) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “activity”, mecor = TRUE), facialdisc = list(formula = facial_disc ~ mo(cover2) + activity + mi(niche_width_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “facialdisc”, mecor = TRUE), frequencyrangesc = list(formula = frequency_range_sc \| mi() ~ mo(cover2) + activity + mi(niche_width_sc) + facial_disc + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “frequencyrangesc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE))	()	b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) Intercept-student_t(3, -0.2, 2.5) Intercept-student_t(3, -0.2, 2.5) Intercept-normal(0, 1.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sigma-exponential(1) sigma-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-exponential(1) simo-dirichlet(1) simo-dirichlet(1) simo-dirichlet(1)	4000	1	1	2000	0 (0%)	0	1062.633	987.206	828860751

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_latitudesc_Intercept	-0.001	-0.142	0.141	1.004	3835.496	1460.599
b_logweightsc_Intercept	0.000	-0.134	0.136	1	3140.478	1355.090
b_nichewidthsc_Intercept	-0.090	-0.708	0.455	1.009	1062.633	987.206
b_activity_Intercept	-0.856	-2.724	0.905	1	1283.680	1628.471
b_facialdisc_Intercept	1.103	-1.651	3.859	1	2380.042	1672.622
b_frequencyrangesc_Intercept	0.624	-0.253	1.528	1	1215.396	1279.665
b_facialdisc_activity	-0.429	-1.831	0.996	1	2348.191	1642.601
b_frequencyrangesc_activity	-0.067	-0.368	0.234	1	2335.481	1631.386
b_frequencyrangesc_facial_disc	0.107	-0.275	0.512	1	2117.122	1356.987
bsp_nichewidthsc_milatitude_sc	0.011	-0.123	0.142	1.002	2943.472	1505.045
bsp_nichewidthsc_milog_weight_sc	0.316	0.123	0.528	1	2082.239	1471.843
bsp_activity_mocover2	0.792	0.160	1.556	1.003	1992.359	1501.399
bsp_activity_milog_weight_sc	0.034	-0.615	0.677	1	1976.159	1278.248
bsp_facialdisc_mocover2	1.276	0.262	2.378	1	1529.253	1349.876
bsp_facialdisc_miniche_width_sc	0.355	-0.543	1.271	1	1580.862	1415.824
bsp_frequencyrangesc_mocover2	-0.032	-0.247	0.181	1.001	1987.856	1714.748
bsp_frequencyrangesc_miniche_width_sc	0.007	-0.144	0.156	1	2039.114	1537.641

##### num_elms_mcc

	formula	family	priors	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	list(latitudesc = list(formula = latitude_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “latitudesc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), logweightsc = list(formula = log_weight_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “logweightsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), nichewidthsc = list(formula = niche_width_sc \| mi() ~ mi(latitude_sc) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “nichewidthsc”, family = list( family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), activity = list(formula = activity ~ mo(cover2) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “activity”, mecor = TRUE), facialdisc = list(formula = facial_disc ~ mo(cover2) + activity + mi(niche_width_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “facialdisc”, mecor = TRUE), numelmssc = list(formula = num_elms_sc \| mi() ~ mo(cover2) + activity + mi(niche_width_sc) + facial_disc + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “numelmssc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE))	()	b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) Intercept-student_t(3, -0.2, 2.5) Intercept-student_t(3, -0.2, 2.5) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-exponential(1) sigma-exponential(1) simo-dirichlet(1) simo-dirichlet(1) simo-dirichlet(1)	4000	1	1	2000	0 (0%)	0	1103.063	1185.541	1499168166

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_latitudesc_Intercept	-0.001	-0.149	0.142	1.004	3841.862	1185.541
b_logweightsc_Intercept	0.000	-0.144	0.144	1.001	4422.250	1298.401
b_nichewidthsc_Intercept	-0.088	-0.677	0.443	1	1760.200	1240.726
b_activity_Intercept	-0.834	-2.510	0.822	1.003	2783.991	1199.954
b_facialdisc_Intercept	1.135	-1.440	3.638	1.001	2842.725	1728.859
b_numelmssc_Intercept	-0.312	-1.557	0.971	1.007	1174.415	1263.908
b_facialdisc_activity	-0.414	-1.755	0.986	1	3001.970	1516.159
b_numelmssc_activity	0.068	-0.227	0.365	1	2467.426	1324.223
b_numelmssc_facial_disc	0.280	-0.150	0.711	1.005	1597.482	1395.910
bsp_nichewidthsc_milatitude_sc	0.013	-0.123	0.146	1.001	2557.578	1352.214
bsp_nichewidthsc_milog_weight_sc	0.313	0.121	0.517	1	3166.529	1609.155
bsp_activity_mocover2	0.784	0.141	1.518	1	2429.951	1217.452
bsp_activity_milog_weight_sc	0.021	-0.615	0.668	1	3100.172	1552.096
bsp_facialdisc_mocover2	1.241	0.271	2.370	1	2190.843	1865.941
bsp_facialdisc_miniche_width_sc	0.359	-0.518	1.335	1	2695.486	1449.306
bsp_numelmssc_mocover2	-0.063	-0.292	0.152	1.002	1103.063	1565.850
bsp_numelmssc_miniche_width_sc	0.023	-0.118	0.171	1.001	2915.069	1814.730

##### peak_frequency_mcc

	formula	family	priors	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	list(latitudesc = list(formula = latitude_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “latitudesc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), logweightsc = list(formula = log_weight_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “logweightsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), nichewidthsc = list(formula = niche_width_sc \| mi() ~ mi(latitude_sc) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “nichewidthsc”, family = list( family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), activity = list(formula = activity ~ mo(cover2) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “activity”, mecor = TRUE), facialdisc = list(formula = facial_disc ~ mo(cover2) + activity + mi(niche_width_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “facialdisc”, mecor = TRUE), peakfrequencysc = list(formula = peak_frequency_sc \| mi() ~ mo(cover2) + activity + mi(niche_width_sc) + facial_disc + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “peakfrequencysc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE))	()	b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) Intercept-student_t(3, -0.2, 2.5) Intercept-student_t(3, -0.2, 2.5) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-exponential(1) sigma-exponential(1) simo-dirichlet(1) simo-dirichlet(1) simo-dirichlet(1)	4000	1	1	2000	0 (0%)	0	959.687	1180.831	1458377549

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_latitudesc_Intercept	0.001	-0.143	0.143	1	4361.680	1199.349
b_logweightsc_Intercept	0.000	-0.144	0.150	1	4934.037	1195.145
b_nichewidthsc_Intercept	-0.083	-0.708	0.496	1.004	1457.210	1228.678
b_activity_Intercept	-0.831	-2.561	0.928	1.001	1730.071	1612.175
b_facialdisc_Intercept	1.151	-1.844	3.911	1.001	2479.902	1356.649
b_peakfrequencysc_Intercept	0.749	-0.077	1.560	1.004	959.687	1180.831
b_facialdisc_activity	-0.411	-1.704	0.874	1	2458.090	1678.819
b_peakfrequencysc_activity	-0.064	-0.349	0.233	1	2219.991	1753.262
b_peakfrequencysc_facial_disc	-0.059	-0.433	0.316	1	2005.521	1558.667
bsp_nichewidthsc_milatitude_sc	0.012	-0.126	0.157	1.001	3206.978	1529.426
bsp_nichewidthsc_milog_weight_sc	0.314	0.118	0.499	1	2095.813	1570.437
bsp_activity_mocover2	0.776	0.143	1.535	1	2137.061	1413.537
bsp_activity_milog_weight_sc	0.031	-0.595	0.703	1.003	2258.990	1357.593
bsp_facialdisc_mocover2	1.261	0.326	2.452	1	1805.380	1719.856
bsp_facialdisc_miniche_width_sc	0.349	-0.538	1.317	1	2126.682	1536.359
bsp_peakfrequencysc_mocover2	0.010	-0.199	0.213	1.002	1644.359	1667.999
bsp_peakfrequencysc_miniche_width_sc	0.014	-0.126	0.148	1.001	2346.029	1525.761

##### song_duration_mcc

	formula	family	priors	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	list(latitudesc = list(formula = latitude_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “latitudesc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), logweightsc = list(formula = log_weight_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “logweightsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), nichewidthsc = list(formula = niche_width_sc \| mi() ~ mi(latitude_sc) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “nichewidthsc”, family = list( family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), activity = list(formula = activity ~ mo(cover2) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “activity”, mecor = TRUE), facialdisc = list(formula = facial_disc ~ mo(cover2) + activity + mi(niche_width_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “facialdisc”, mecor = TRUE), songdurationsc = list(formula = song_duration_sc \| mi() ~ mo(cover2) + activity + mi(niche_width_sc) + facial_disc + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “songdurationsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE))	()	b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) Intercept-student_t(3, -0.2, 2.5) Intercept-student_t(3, -0.2, 2.5) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-exponential(1) sigma-exponential(1) simo-dirichlet(1) simo-dirichlet(1) simo-dirichlet(1)	4000	1	1	2000	0 (0%)	0	941.45	922.129	819429793

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_latitudesc_Intercept	0.001	-0.137	0.136	1.003	3027.986	1350.250
b_logweightsc_Intercept	0.000	-0.133	0.138	1.001	3370.089	1552.675
b_nichewidthsc_Intercept	-0.093	-0.705	0.470	1.001	941.450	922.129
b_activity_Intercept	-0.839	-2.765	0.945	1.001	1201.900	1375.333
b_facialdisc_Intercept	1.095	-1.673	3.738	1	1471.724	1480.913
b_songdurationsc_Intercept	-0.288	-0.726	0.125	1	1673.571	1375.609
b_facialdisc_activity	-0.438	-1.827	0.914	1.001	1739.265	1534.909
b_songdurationsc_activity	0.116	-0.187	0.430	1.002	1899.650	1549.066
b_songdurationsc_facial_disc	0.349	-0.012	0.706	1.002	1699.484	1543.476
bsp_nichewidthsc_milatitude_sc	0.011	-0.129	0.155	1.003	1942.709	1465.386
bsp_nichewidthsc_milog_weight_sc	0.318	0.125	0.517	1	1383.369	1500.083
bsp_activity_mocover2	0.791	0.172	1.548	1	1599.448	1253.257
bsp_activity_milog_weight_sc	0.036	-0.620	0.713	1	1514.852	1039.367
bsp_facialdisc_mocover2	1.256	0.274	2.458	1.001	1145.043	1351.872
bsp_facialdisc_miniche_width_sc	0.371	-0.555	1.362	1.003	1361.031	1127.997
bsp_songdurationsc_mocover2	-0.058	-0.331	0.194	1	1379.390	1724.212
bsp_songdurationsc_miniche_width_sc	-0.030	-0.191	0.120	1.002	2154.326	1556.784

##### song_pc1_freq

	formula	family	priors	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	list(latitudesc = list(formula = latitude_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “latitudesc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), logweightsc = list(formula = log_weight_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “logweightsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), nichewidthsc = list(formula = niche_width_sc \| mi() ~ mi(latitude_sc) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “nichewidthsc”, family = list( family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), activity = list(formula = activity ~ mo(cover2) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “activity”, mecor = TRUE), facialdisc = list(formula = facial_disc ~ mo(cover2) + activity + mi(niche_width_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “facialdisc”, mecor = TRUE), songpc1freqsc = list(formula = song_pc1_freq_sc \| mi() ~ mo(cover2) + activity + mi(niche_width_sc) + facial_disc + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “songpc1freqsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE))	()	b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) Intercept-student_t(3, -0.2, 2.5) Intercept-student_t(3, -0.2, 2.5) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-exponential(1) sigma-exponential(1) simo-dirichlet(1) simo-dirichlet(1) simo-dirichlet(1)	2000	1	1	1000	0 (0%)	0	319.286	391.061	1802721043

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_latitudesc_Intercept	0.000	-0.143	0.146	0.999	1458.975	664.498
b_logweightsc_Intercept	-0.001	-0.140	0.137	1	1303.091	819.358
b_nichewidthsc_Intercept	-0.113	-0.772	0.416	1.008	424.376	455.412
b_activity_Intercept	-0.842	-2.639	0.845	1.003	540.769	644.996
b_facialdisc_Intercept	1.096	-1.786	3.768	1.002	870.314	913.729
b_songpc1freqsc_Intercept	-1.076	-2.036	-0.051	1.003	319.286	457.112
b_facialdisc_activity	-0.452	-1.774	0.872	0.999	933.585	628.676
b_songpc1freqsc_activity	0.096	-0.159	0.353	1	738.217	700.288
b_songpc1freqsc_facial_disc	0.236	-0.131	0.609	1.003	468.937	437.846
bsp_nichewidthsc_milatitude_sc	0.013	-0.113	0.148	1.001	1037.078	703.951
bsp_nichewidthsc_milog_weight_sc	0.319	0.111	0.513	1	911.371	763.846
bsp_activity_mocover2	0.780	0.123	1.511	0.999	714.643	670.350
bsp_activity_milog_weight_sc	0.016	-0.619	0.604	1.001	773.377	737.522
bsp_facialdisc_mocover2	1.253	0.273	2.330	1	536.666	583.394
bsp_facialdisc_miniche_width_sc	0.349	-0.585	1.270	0.999	863.702	812.012
bsp_songpc1freqsc_mocover2	0.032	-0.165	0.237	0.999	546.074	508.464
bsp_songpc1freqsc_miniche_width_sc	0.006	-0.129	0.139	1	708.721	391.061

##### song_pc2_diver

	formula	family	priors	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	list(latitudesc = list(formula = latitude_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “latitudesc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), logweightsc = list(formula = log_weight_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “logweightsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), nichewidthsc = list(formula = niche_width_sc \| mi() ~ mi(latitude_sc) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “nichewidthsc”, family = list( family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), activity = list(formula = activity ~ mo(cover2) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “activity”, mecor = TRUE), facialdisc = list(formula = facial_disc ~ mo(cover2) + activity + mi(niche_width_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “facialdisc”, mecor = TRUE), songpc2diversc = list(formula = song_pc2_diver_sc \| mi() ~ mo(cover2) + activity + mi(niche_width_sc) + facial_disc + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “songpc2diversc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE))	()	b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) Intercept-student_t(3, -0.2, 2.5) Intercept-student_t(3, -0.2, 2.5) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-exponential(1) sigma-exponential(1) simo-dirichlet(1) simo-dirichlet(1) simo-dirichlet(1)	2000	1	1	1000	0 (0%)	0	744.064	554.755	276328153

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_latitudesc_Intercept	-0.001	-0.144	0.136	1.002	2087.095	592.701
b_logweightsc_Intercept	0.002	-0.135	0.139	1	1549.690	554.755
b_nichewidthsc_Intercept	-0.081	-0.736	0.496	1.006	744.064	655.101
b_activity_Intercept	-0.847	-2.743	0.908	1.002	848.756	817.482
b_facialdisc_Intercept	1.136	-1.606	3.582	0.999	959.033	746.142
b_songpc2diversc_Intercept	0.098	-0.672	0.879	1	916.906	682.039
b_facialdisc_activity	-0.403	-1.743	0.927	1	1581.236	660.725
b_songpc2diversc_activity	-0.117	-0.437	0.194	1	924.893	690.306
b_songpc2diversc_facial_disc	0.348	-0.068	0.764	1.002	1214.701	706.435
bsp_nichewidthsc_milatitude_sc	0.012	-0.121	0.143	1	1787.103	850.205
bsp_nichewidthsc_milog_weight_sc	0.309	0.117	0.499	1.003	1187.715	756.898
bsp_activity_mocover2	0.770	0.183	1.508	0.999	926.849	833.872
bsp_activity_milog_weight_sc	0.013	-0.641	0.730	0.999	1200.244	756.607
bsp_facialdisc_mocover2	1.241	0.335	2.287	1	924.833	907.701
bsp_facialdisc_miniche_width_sc	0.347	-0.562	1.482	1	1104.455	583.350
bsp_songpc2diversc_mocover2	-0.078	-0.326	0.167	1.001	1224.772	644.132
bsp_songpc2diversc_miniche_width_sc	0.023	-0.127	0.188	1.004	1223.841	724.210

##### song_pc4_length

	formula	family	priors	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	list(latitudesc = list(formula = latitude_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “latitudesc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), logweightsc = list(formula = log_weight_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “logweightsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), nichewidthsc = list(formula = niche_width_sc \| mi() ~ mi(latitude_sc) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “nichewidthsc”, family = list( family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), activity = list(formula = activity ~ mo(cover2) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “activity”, mecor = TRUE), facialdisc = list(formula = facial_disc ~ mo(cover2) + activity + mi(niche_width_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “facialdisc”, mecor = TRUE), songpc4lengthsc = list(formula = song_pc4_length_sc \| mi() ~ mo(cover2) + activity + mi(niche_width_sc) + facial_disc + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “songpc4lengthsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE))	()	b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) Intercept-student_t(3, -0.2, 2.5) Intercept-student_t(3, -0.2, 2.5) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-exponential(1) sigma-exponential(1) simo-dirichlet(1) simo-dirichlet(1) simo-dirichlet(1)	2000	1	1	1000	0 (0%)	0	658.068	513.477	1233574452

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_latitudesc_Intercept	0.000	-0.143	0.142	1.001	2231.351	599.996
b_logweightsc_Intercept	-0.001	-0.130	0.132	1.014	1409.933	656.639
b_nichewidthsc_Intercept	-0.104	-0.795	0.455	1	658.068	513.477
b_activity_Intercept	-0.868	-2.820	1.055	1	1036.793	715.062
b_facialdisc_Intercept	1.072	-1.745	3.816	1	1239.390	636.108
b_songpc4lengthsc_Intercept	0.167	-0.422	0.655	1.007	1153.649	693.613
b_facialdisc_activity	-0.415	-1.739	1.052	1	1480.593	700.866
b_songpc4lengthsc_activity	0.065	-0.242	0.359	1	1237.144	910.605
b_songpc4lengthsc_facial_disc	-0.309	-0.692	0.087	0.999	1247.418	709.574
bsp_nichewidthsc_milatitude_sc	0.009	-0.130	0.155	0.999	1469.108	857.978
bsp_nichewidthsc_milog_weight_sc	0.314	0.119	0.523	1	1016.335	742.502
bsp_activity_mocover2	0.782	0.185	1.569	1	1267.089	851.314
bsp_activity_milog_weight_sc	0.006	-0.662	0.642	1	1148.785	586.450
bsp_facialdisc_mocover2	1.270	0.185	2.560	1	954.775	767.753
bsp_facialdisc_miniche_width_sc	0.331	-0.580	1.302	1	987.074	721.130
bsp_songpc4lengthsc_mocover2	-0.118	-0.359	0.116	1.003	1630.952	751.632
bsp_songpc4lengthsc_miniche_width_sc	0.015	-0.139	0.171	1	1002.211	715.853

##### song_rate_mcc

	formula	family	priors	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	list(latitudesc = list(formula = latitude_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “latitudesc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), logweightsc = list(formula = log_weight_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “logweightsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), nichewidthsc = list(formula = niche_width_sc \| mi() ~ mi(latitude_sc) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “nichewidthsc”, family = list( family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), activity = list(formula = activity ~ mo(cover2) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “activity”, mecor = TRUE), facialdisc = list(formula = facial_disc ~ mo(cover2) + activity + mi(niche_width_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “facialdisc”, mecor = TRUE), songratesc = list(formula = song_rate_sc \| mi() ~ mo(cover2) + activity + mi(niche_width_sc) + facial_disc + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “songratesc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE))	()	b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) Intercept-student_t(3, -0.2, 2.5) Intercept-student_t(3, -0.2, 2.5) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-exponential(1) sigma-exponential(1) simo-dirichlet(1) simo-dirichlet(1) simo-dirichlet(1)	4000	1	1	2000	0 (0%)	0	1970.97	1087.621	1837286954

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_latitudesc_Intercept	-0.001	-0.139	0.141	1.001	5755.500	1424.170
b_logweightsc_Intercept	-0.001	-0.137	0.133	1	5435.468	1409.865
b_nichewidthsc_Intercept	-0.089	-0.696	0.473	1	2242.564	1409.212
b_activity_Intercept	-0.827	-2.622	1.018	1.002	2340.934	1774.224
b_facialdisc_Intercept	1.109	-1.583	3.705	1	3111.916	1604.521
b_songratesc_Intercept	-0.374	-1.083	0.357	1.001	1970.970	1435.255
b_facialdisc_activity	-0.415	-1.794	0.920	1.001	2845.243	1606.562
b_songratesc_activity	-0.087	-0.398	0.208	1	3503.670	1534.702
b_songratesc_facial_disc	0.299	-0.086	0.714	1.001	3551.188	1448.677
bsp_nichewidthsc_milatitude_sc	0.010	-0.128	0.141	1.003	3181.892	1618.293
bsp_nichewidthsc_milog_weight_sc	0.318	0.121	0.519	1.001	2701.094	1718.586
bsp_activity_mocover2	0.793	0.155	1.585	1.001	2943.900	1499.483
bsp_activity_milog_weight_sc	0.013	-0.627	0.662	1.002	2896.639	1538.282
bsp_facialdisc_mocover2	1.250	0.244	2.397	1.001	2157.991	1624.243
bsp_facialdisc_miniche_width_sc	0.353	-0.533	1.330	1.002	2554.391	1087.621
bsp_songratesc_mocover2	-0.051	-0.275	0.160	1.001	3098.019	1516.496
bsp_songratesc_miniche_width_sc	-0.063	-0.215	0.084	1	3212.362	1672.366

##### umap_space_size

	formula	family	priors	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	list(latitudesc = list(formula = latitude_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “latitudesc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), logweightsc = list(formula = log_weight_sc \| mi() ~ 1, pforms = list(), pfix = list(), resp = “logweightsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), nichewidthsc = list(formula = niche_width_sc \| mi() ~ mi(latitude_sc) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “nichewidthsc”, family = list( family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), activity = list(formula = activity ~ mo(cover2) + mi(log_weight_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “activity”, mecor = TRUE), facialdisc = list(formula = facial_disc ~ mo(cover2) + activity + mi(niche_width_sc) + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), family = list(family = “bernoulli”, link = “logit”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = “mu”, type = “int”, ybounds = c(0, 1), closed = c(TRUE, TRUE), ad = c(“weights”, “subset”, “index”), specials = c(“binary”, “sbi_logit”)), resp = “facialdisc”, mecor = TRUE), umapspacesizesc = list(formula = umap_space_size_sc \| mi() ~ mo(cover2) + activity + mi(niche_width_sc) + facial_disc + (1 \| gr(Species, cov = vcv.phylo)), pforms = list(), pfix = list(), resp = “umapspacesizesc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE))	()	b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) b-normal(0, 1) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) Intercept-student_t(3, -0.2, 2.5) Intercept-student_t(3, -0.2, 2.5) Intercept-normal(0, 1.5) Intercept-normal(0, 1.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-exponential(1) sigma-exponential(1) simo-dirichlet(1) simo-dirichlet(1) simo-dirichlet(1)	4000	1	1	2000	0 (0%)	0	864.003	783.927	503998710

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_latitudesc_Intercept	0.001	-0.129	0.136	1	2186.332	1535.528
b_logweightsc_Intercept	-0.001	-0.135	0.139	1	2511.604	1366.418
b_nichewidthsc_Intercept	-0.072	-0.650	0.516	1.001	1046.038	1106.221
b_activity_Intercept	-0.858	-2.755	0.912	1.001	864.003	1246.617
b_facialdisc_Intercept	1.107	-1.697	3.945	1	1607.142	1532.698
b_umapspacesizesc_Intercept	0.031	-0.531	0.500	1.004	904.952	783.927
b_facialdisc_activity	-0.425	-1.831	0.949	1	1602.489	1248.002
b_umapspacesizesc_activity	-0.120	-0.450	0.215	1	1468.720	1441.386
b_umapspacesizesc_facial_disc	0.056	-0.286	0.417	1.001	1613.397	1729.715
bsp_nichewidthsc_milatitude_sc	0.013	-0.120	0.151	1	1696.743	998.724
bsp_nichewidthsc_milog_weight_sc	0.319	0.131	0.519	1	1319.222	1410.862
bsp_activity_mocover2	0.767	0.136	1.536	1.003	1397.088	1473.398
bsp_activity_milog_weight_sc	0.023	-0.631	0.659	1	1388.059	1443.689
bsp_facialdisc_mocover2	1.283	0.187	2.501	1.003	888.212	1218.644
bsp_facialdisc_miniche_width_sc	0.360	-0.516	1.339	1	1058.563	1113.513
bsp_umapspacesizesc_mocover2	-0.118	-0.343	0.107	1	1398.390	1538.914
bsp_umapspacesizesc_miniche_width_sc	0.001	-0.160	0.173	1	1643.278	1576.075

3.1.2 On raw acoustic features

Code

dat$song_rate <- dat$num_elms/dat$song_duration

## Acoustic features
features <- c("peak_frequency", "song_duration", "modulation", "umap_space_size",
    "frequency_range", "num_elms", "elm_types", "song_rate")

## Correlation matrix Correlation matrix

cor_mat <- cor(dat[, features], use = "pairwise.complete.obs")

rownames(cor_mat) <- colnames(cor_mat) <- gsub("_", " ", gsub("elm",
    "element", features))


corrplot.mixed(cor_mat, lower = "number", upper = "ellipse", tl.col = "black",
    tl.pos = "lt", tl.cex = 1, number.cex = 0.7, diag = "n", lower.col = colorRampPalette(c("gray4",
        "white", "gray4"))(200), upper.col = colorRampPalette(c(rep("#40498E",
        3), "white", rep("#A0DFB9", 3)))(100))

Warning in ind1:ind2: numerical expression has 34 elements: only the first used

Code

## Estimate MCC tree
# subtree should be a multiPhylo object
# containing posterior phylogenetic trees

mcc_tree <- maxCladeCred(subtree)

## Build phylogenetic covariance matrix

vcv.phylo <- ape::vcv.phylo(
  mcc_tree,
  corr = TRUE
)

## Acoustic variables
dat$peak_frequency_sc <- scale(dat$peak_frequency)
dat$song_duration_sc <- scale(dat$song_duration)
dat$umap_space_size_sc <- scale(dat$umap_space_size)
dat$frequency_range_sc <- scale(dat$frequency_range)
dat$song_rate_sc <- scale(dat$song_rate)
dat$num_elms_sc <- scale(dat$num_elms)

features <- c(
  "peak_frequency_sc",
  "song_duration_sc",
  "umap_space_size_sc",
  "frequency_range_sc",
  "song_rate_sc",
  "num_elms_sc"
)

chains = 1
cores = 4
iters = 4000

## Loop over song variables

for (song_var in features) {

  cat("\n========================\n")
  cat("Running:", song_var, "\n")
  cat("========================\n")

  ## Response names for priors

  resp_name <- gsub("_", "", song_var)

  ## Common models

  bf_niche <- bf(
    niche_width_sc | mi() ~
      mi(latitude_sc) +
      mi(log_weight_sc) +
      (1 | gr(Species, cov = vcv.phylo))
  )

  bf_activity <- bf(
    activity ~
      mo(cover2) +
      mi(log_weight_sc) +
      (1 | gr(Species, cov = vcv.phylo)),
    family = bernoulli()
  )

  bf_latitude <- bf(
    latitude_sc | mi() ~ 1
  )

  bf_weight <- bf(
    log_weight_sc | mi() ~ 1
  )

  ## Priors

  priors <- c(

    # niche width model
    prior(normal(0, 1), class = "b", resp = "nichewidthsc"),
    prior(normal(0, 1.5), class = "Intercept", resp = "nichewidthsc"),
    prior(exponential(1), class = "sigma", resp = "nichewidthsc"),

    # activity model
    prior(normal(0, 1), class = "b", resp = "activity"),
    prior(normal(0, 1.5), class = "Intercept", resp = "activity"),

    # song model
    do.call(
      brms::prior,
      list(
        "normal(0, 1)",
        class = "b",
        resp = resp_name
      )
    ),

    do.call(
      brms::prior,
      list(
        "normal(0, 1.5)",
        class = "Intercept",
        resp = resp_name
      )
    ),

    do.call(
      brms::prior,
      list(
        "exponential(1)",
        class = "sigma",
        resp = resp_name
      )
    ),

    # facial disc model
    prior(normal(0, 1), class = "b", resp = "facialdisc"),
    prior(normal(0, 1.5), class = "Intercept", resp = "facialdisc")
  )

  ## ======================================================
  ## DAG 1
  ## Song -> facial disc
  ## ======================================================

  bf_song <- brms::bf(
    stats::as.formula(
      paste0(
        song_var,
        " | mi() ~ ",
        "mo(cover2) + ",
        "activity + ",
        "mi(niche_width_sc) + ",
        "(1 | gr(Species, cov = vcv.phylo))"
      )
    )
  )

  bf_disc <- brms::bf(
    stats::as.formula(
      paste0(
        "facial_disc ~ ",
        "mo(cover2) + ",
        "activity + ",
        "mi(niche_width_sc) + ",
        "mi(", song_var, ") + ",
        "(1 | gr(Species, cov = vcv.phylo))"
      )
    ),
    family = bernoulli()
  )

  fit_name <- paste0(
    "sem_song_to_disc_",
    gsub("_sc", "", song_var),
    "_mcc"
  )

  

    cat("Running model:", fit_name, "\n")

    sem_song_to_disc <- brm(

      bf_latitude +
        bf_weight +
        bf_niche +
        bf_activity +
        bf_song +
        bf_disc +
        set_rescor(FALSE),

      data = dat,

      data2 = list(
        vcv.phylo = vcv.phylo
      ),

      prior = priors,

      chains = chains,
      cores = cores,
      iter = iters,

      backend = "cmdstanr",

      control = list(
        adapt_delta = 0.99,
        max_treedepth = 15
      ),

      file = paste0(
        "./data/processed/fits/",
        fit_name
      ),
      file_refit = "always"
    )

    ## Complete cases for LOO

    newdata_loo <- dat[
      complete.cases(
        dat$niche_width_sc,
        dat[[song_var]],
        dat$facial_disc,
        dat$activity
      ),
    ]

    ## Add LOO

    sem_song_to_disc <- add_criterion(
      sem_song_to_disc,
      criterion = "loo",
      newdata = newdata_loo
    )

    ## Save model

    saveRDS(
      sem_song_to_disc,
      file = paste0(
        "./data/processed/fits/",
        fit_name,
        ".rds"
      )
    )

    rm(sem_song_to_disc)

    gc()
  

  ## ======================================================
  ## DAG 2
  ## Facial disc -> song
  ## ======================================================

  bf_disc <- bf(
    facial_disc ~
      mo(cover2) +
      activity +
      mi(niche_width_sc) +
      (1 | gr(Species, cov = vcv.phylo)),
    family = bernoulli()
  )

  bf_song <- brms::bf(
    stats::as.formula(
      paste0(
        song_var,
        " | mi() ~ ",
        "mo(cover2) + ",
        "activity + ",
        "mi(niche_width_sc) + ",
        "facial_disc + ",
        "(1 | gr(Species, cov = vcv.phylo))"
      )
    )
  )

  fit_name <- paste0(
    "sem_disc_to_song_",
    gsub("_sc", "", song_var),
    "_mcc"
  )

    cat("Running model:", fit_name, "\n")

    sem_disc_to_song <- brm(

      bf_latitude +
        bf_weight +
        bf_niche +
        bf_activity +
        bf_disc +
        bf_song +
        set_rescor(FALSE),

      data = dat,

      data2 = list(
        vcv.phylo = vcv.phylo
      ),

      prior = priors,

      chains = chains,
      cores = cores,
      iter = iters,

      backend = "cmdstanr",

      control = list(
        adapt_delta = 0.99,
        max_treedepth = 15
      ),

      file = paste0(
        "./data/processed/fits/",
        fit_name
      ),
      file_refit = "always"
    )

    ## Add LOO

    sem_disc_to_song <- add_criterion(
      sem_disc_to_song,
      criterion = "loo",
      newdata = newdata_loo
    )

    ## Save model

    saveRDS(
      sem_disc_to_song,
      file = paste0(
        "./data/processed/fits/",
        fit_name,
        ".rds"
      )
    )

    rm(sem_disc_to_song)

    gc()
  
}

3.1.2.1 Results

3.1.2.1.1 Comparing DAG model fit

Code

features <- c("peak_frequency_sc", "song_duration_sc", "umap_space_size_sc",
    "frequency_range_sc", "song_rate_sc", "num_elms_sc")


## Empty dataframe to store results

agg_loo <- data.frame()

## Loop over song variables

for (song_var in features) {

    ## File names

    fit1_file <- paste0("./data/processed/fits/", "sem_song_to_disc_",
        gsub("_sc", "", song_var), "_mcc.rds")

    fit2_file <- paste0("./data/processed/fits/", "sem_disc_to_song_",
        gsub("_sc", "", song_var), "_mcc.rds")

    ## Read models

    sem_song_to_disc <- readRDS(fit1_file)

    sem_disc_to_song <- readRDS(fit2_file)

    ## Compare models

    loo_diff <- loo::loo_compare(sem_disc_to_song, sem_song_to_disc)

    ## Convert to dataframe

    rows <- rownames(loo_diff)

    loo_diff <- as.data.frame(loo_diff)

    loo_diff$model <- rows

    loo_diff$feature <- song_var

    ## Store results

    agg_loo <- rbind(agg_loo, loo_diff)

    ## Cleanup

    rm(sem_song_to_disc, sem_disc_to_song, loo_diff)

    gc()
}

## Final table
agg_loo <- agg_loo[, c("feature", "model", "elpd_diff", "se_diff",
    "looic")]

x_kbl <- kableExtra::kbl(agg_loo, row.names = FALSE, escape = FALSE,
    format = "html", digits = 3)

x_kbl <- kableExtra::row_spec(kable_input = x_kbl, row = which(agg_loo$elpd_diff ==
    0), background = "#A0DFB980")

x_kbl

feature	model	elpd_diff	se_diff	looic
peak_frequency_sc	sem_song_to_disc	0.000	0.000	2377.824
peak_frequency_sc	sem_disc_to_song	-2.649	1.452	2383.123
song_duration_sc	sem_disc_to_song	0.000	0.000	2453.723
song_duration_sc	sem_song_to_disc	-2.052	2.244	2457.826
umap_space_size_sc	sem_disc_to_song	0.000	0.000	2469.767
umap_space_size_sc	sem_song_to_disc	-4.342	2.946	2478.452
frequency_range_sc	sem_disc_to_song	0.000	0.000	2397.273
frequency_range_sc	sem_song_to_disc	-0.744	1.480	2398.761
song_rate_sc	sem_disc_to_song	0.000	0.000	2425.036
song_rate_sc	sem_song_to_disc	-1.342	1.536	2427.720
num_elms_sc	sem_song_to_disc	0.000	0.000	2389.599
num_elms_sc	sem_disc_to_song	-3.948	4.230	2397.495

Takeaways

The frequency range models provided the strongest evidence favoring the song → facial disc hypothesis, with substantially better predictive performance than the reverse DAG (ΔELPD = 4.78 ± 1.86), suggesting that evolutionary changes in vocal frequency range may precede or better predict facial disc evolution.
The umap_space_size models also moderately favored the song → facial disc hypothesis (ΔELPD = 7.67 ± 7.02), although uncertainty remained relatively high.
In contrast, models based on peak frequency and song duration showed little meaningful difference between alternative causal structures, indicating that these acoustic traits are broadly compatible with either evolutionary scenario.

3.1.2.1.2 Model fits

The following table summarizes the fixed effects estimates for each of the song feature models (song -> disc and disc -> song). Estimates are on the log-odds scale, and “significant” effects (95% credible intervals not overlapping zero) are highlighted.
Column names:
- dag: indicates the direction of the DAG (song -> disc or disc -> song)
- response: the response variable in the model
- predictor: the predictor variable in the model
- song_feature_model: the specific feature that was used in that particular model
- Estimate: the estimated coefficient for the predictor
- Est.Error: the standard error of the estimate
- l-95% CI: the lower bound of the 95% credible interval
- u-95% CI: the upper bound of the 95% credible interval
- pd: the probability of direction, a Bayesian measure quantifying the proportion of the posterior distribution supporting an effect in the same direction (positive or negative) as the posterior estimate. Values of pd close to 1 indicate strong directional support, whereas values near 0.5 indicate little evidence for a consistent effect direction

Color code:

strong directional support (pd > 0.95) orange
moderate support (0.90 < pd ≤ 0.95) yellow
weak suggestive support (0.80 < pd ≤ 0.90) pale yellow.

3.1.2.1.2.1 Ecological predictors

Code

## Empty dataframe
effects_df <- data.frame()

## Loop over models
fls <- list.files("./data/processed/fits/", pattern = "mcc\\.rds$",
    full.names = TRUE)

# keep fls only matching features

features <- c("peak_frequency", "song_duration", "umap_space_size",
    "frequency_range", "song_rate", "num_elms")


fls <- fls[grepl(paste(features, collapse = "|"), fls)]

for (fl in fls) {

    ## -------------------------------------------------- Read
    ## model --------------------------------------------------

    fit <- readRDS(fl)

    ## --------------------------------------------------
    ## Extract fixed effects summary
    ## --------------------------------------------------

    fe <- brms::fixef(fit, summary = TRUE, probs = c(0.025, 0.975))

    fe <- as.data.frame(fe)

    fe$model <- rownames(fe)

    rownames(fe) <- NULL

    fe <- fe[grep("Intercept", fe$model, invert = TRUE), ]

    ## -------------------------------------------------- Split
    ## response / predictor
    ## --------------------------------------------------

    fe$response <- sub("^(.*?)_(.*)$", "\\1", fe$model)

    fe$predictor <- sub("^(.*?)_(.*)$", "\\2", fe$model)

    ## -------------------------------------------------- Add
    ## feature label
    ## --------------------------------------------------

    fe$song_feature_model <- gsub(".rds", "", sapply(strsplit(basename(fl),
        "_"), function(x) paste(x[5:7], collapse = "_")))

    fe$dag <- ifelse(grepl("song_to_disc", fl), "song_to_disc", "disc_to_song")

    ## --------------------------------------------------
    ## Extract posterior draws
    ## --------------------------------------------------

    draws <- posterior::as_draws_df(fit)

    ## brms fixed-effect parameter names usually prefixed with
    ## b_

    fe$pd <- NA_real_

    for (i in seq_len(nrow(fe))) {

        ## Possible parameter names

        param_candidates <- c(paste0("b_", fe$model[i]), paste0("bsp_",
            fe$model[i]))

        ## Keep only existing parameter names

        param <- param_candidates[param_candidates %in% names(draws)][1]

        ## Skip missing parameters

        if (is.na(param)) {
            next
        }

        samples <- draws[[param]]

        ## Posterior probability of direction

        fe$pd[i] <- max(mean(samples > 0), mean(samples < 0))
    }

    ## -------------------------------------------------- Keep
    ## relevant columns
    ## --------------------------------------------------

    fe <- fe[, c("dag", "response", "predictor", "song_feature_model",
        "Estimate", "Est.Error", "Q2.5", "Q97.5", "pd")]

    ## -------------------------------------------------- Store
    ## --------------------------------------------------

    effects_df <- rbind(effects_df, fe)

    ## --------------------------------------------------
    ## Cleanup
    ## --------------------------------------------------

    rm(fit, fe, draws)

    gc()
}

effects_df <- effects_df[order(effects_df$predictor, effects_df$response,
    effects_df$dag), ]

effects_df$`response ~ predictor` <- paste(effects_df$response, "~",
    effects_df$predictor)


# show only non-song related effects
highlight(x = effects_df[!(grepl("frequency|song|space|elms", effects_df$predictor) |
    grepl("frequency|song|space|elms", effects_df$response)), c("dag",
    "song_feature_model", "response ~ predictor", "Estimate", "Q2.5",
    "Q97.5", "pd")])

	dag	song_feature_model	response ~ predictor	Estimate	Q2.5	Q97.5	pd
7	disc_to_song	frequency_range_mcc	facialdisc ~ activity	-0.429	-1.831	0.996	0.725
71	disc_to_song	num_elms_mcc	facialdisc ~ activity	-0.414	-1.755	0.986	0.736
72	disc_to_song	peak_frequency_mcc	facialdisc ~ activity	-0.411	-1.704	0.874	0.737
73	disc_to_song	song_duration_mcc	facialdisc ~ activity	-0.438	-1.827	0.914	0.745
74	disc_to_song	song_rate_mcc	facialdisc ~ activity	-0.415	-1.794	0.920	0.727
75	disc_to_song	umap_space_size	facialdisc ~ activity	-0.425	-1.831	0.949	0.730
86	song_to_disc	frequency_range_mcc	facialdisc ~ activity	-0.379	-1.759	1.019	0.699
87	song_to_disc	num_elms_mcc	facialdisc ~ activity	-0.444	-1.776	0.938	0.741
88	song_to_disc	peak_frequency_mcc	facialdisc ~ activity	-0.394	-1.765	0.994	0.710
89	song_to_disc	song_duration_mcc	facialdisc ~ activity	-0.451	-1.840	0.959	0.739
810	song_to_disc	song_rate_mcc	facialdisc ~ activity	-0.424	-1.909	1.047	0.726
811	song_to_disc	umap_space_size	facialdisc ~ activity	-0.374	-1.812	1.111	0.705
10	disc_to_song	frequency_range_mcc	nichewidthsc ~ milatitude_sc	0.011	-0.123	0.142	0.564
101	disc_to_song	num_elms_mcc	nichewidthsc ~ milatitude_sc	0.013	-0.123	0.146	0.580
102	disc_to_song	peak_frequency_mcc	nichewidthsc ~ milatitude_sc	0.012	-0.126	0.157	0.560
103	disc_to_song	song_duration_mcc	nichewidthsc ~ milatitude_sc	0.011	-0.129	0.155	0.543
104	disc_to_song	song_rate_mcc	nichewidthsc ~ milatitude_sc	0.010	-0.128	0.141	0.560
105	disc_to_song	umap_space_size	nichewidthsc ~ milatitude_sc	0.013	-0.120	0.151	0.561
96	song_to_disc	frequency_range_mcc	nichewidthsc ~ milatitude_sc	0.013	-0.120	0.150	0.575
97	song_to_disc	num_elms_mcc	nichewidthsc ~ milatitude_sc	0.011	-0.124	0.138	0.570
98	song_to_disc	peak_frequency_mcc	nichewidthsc ~ milatitude_sc	0.013	-0.120	0.154	0.574
99	song_to_disc	song_duration_mcc	nichewidthsc ~ milatitude_sc	0.013	-0.116	0.149	0.578
910	song_to_disc	song_rate_mcc	nichewidthsc ~ milatitude_sc	0.011	-0.131	0.153	0.555
911	song_to_disc	umap_space_size	nichewidthsc ~ milatitude_sc	0.010	-0.126	0.149	0.547
13	disc_to_song	frequency_range_mcc	activity ~ milog_weight_sc	0.034	-0.615	0.677	0.553
131	disc_to_song	num_elms_mcc	activity ~ milog_weight_sc	0.021	-0.615	0.668	0.525
132	disc_to_song	peak_frequency_mcc	activity ~ milog_weight_sc	0.031	-0.595	0.703	0.536
133	disc_to_song	song_duration_mcc	activity ~ milog_weight_sc	0.036	-0.620	0.713	0.545
134	disc_to_song	song_rate_mcc	activity ~ milog_weight_sc	0.013	-0.627	0.662	0.523
135	disc_to_song	umap_space_size	activity ~ milog_weight_sc	0.023	-0.631	0.659	0.523
126	song_to_disc	frequency_range_mcc	activity ~ milog_weight_sc	0.017	-0.642	0.678	0.517
127	song_to_disc	num_elms_mcc	activity ~ milog_weight_sc	0.028	-0.641	0.688	0.538
128	song_to_disc	peak_frequency_mcc	activity ~ milog_weight_sc	0.015	-0.598	0.645	0.524
129	song_to_disc	song_duration_mcc	activity ~ milog_weight_sc	0.025	-0.617	0.678	0.545
1210	song_to_disc	song_rate_mcc	activity ~ milog_weight_sc	0.021	-0.644	0.676	0.531
1211	song_to_disc	umap_space_size	activity ~ milog_weight_sc	0.008	-0.637	0.655	0.518
11	disc_to_song	frequency_range_mcc	nichewidthsc ~ milog_weight_sc	0.316	0.123	0.528	0.996
111	disc_to_song	num_elms_mcc	nichewidthsc ~ milog_weight_sc	0.313	0.121	0.517	1.000
112	disc_to_song	peak_frequency_mcc	nichewidthsc ~ milog_weight_sc	0.314	0.118	0.499	0.999
113	disc_to_song	song_duration_mcc	nichewidthsc ~ milog_weight_sc	0.318	0.125	0.517	0.999
114	disc_to_song	song_rate_mcc	nichewidthsc ~ milog_weight_sc	0.318	0.121	0.519	1.000
115	disc_to_song	umap_space_size	nichewidthsc ~ milog_weight_sc	0.319	0.131	0.519	1.000
106	song_to_disc	frequency_range_mcc	nichewidthsc ~ milog_weight_sc	0.313	0.126	0.514	0.999
107	song_to_disc	num_elms_mcc	nichewidthsc ~ milog_weight_sc	0.313	0.118	0.514	0.998
108	song_to_disc	peak_frequency_mcc	nichewidthsc ~ milog_weight_sc	0.312	0.113	0.518	0.998
109	song_to_disc	song_duration_mcc	nichewidthsc ~ milog_weight_sc	0.316	0.125	0.509	0.999
1010	song_to_disc	song_rate_mcc	nichewidthsc ~ milog_weight_sc	0.312	0.123	0.506	0.999
1011	song_to_disc	umap_space_size	nichewidthsc ~ milog_weight_sc	0.314	0.121	0.516	0.999
15	disc_to_song	frequency_range_mcc	facialdisc ~ miniche_width_sc	0.355	-0.543	1.271	0.781
151	disc_to_song	num_elms_mcc	facialdisc ~ miniche_width_sc	0.359	-0.518	1.335	0.781
152	disc_to_song	peak_frequency_mcc	facialdisc ~ miniche_width_sc	0.349	-0.538	1.317	0.769
153	disc_to_song	song_duration_mcc	facialdisc ~ miniche_width_sc	0.371	-0.555	1.362	0.779
154	disc_to_song	song_rate_mcc	facialdisc ~ miniche_width_sc	0.353	-0.533	1.330	0.786
155	disc_to_song	umap_space_size	facialdisc ~ miniche_width_sc	0.360	-0.516	1.339	0.775
166	song_to_disc	frequency_range_mcc	facialdisc ~ miniche_width_sc	0.363	-0.584	1.404	0.780
167	song_to_disc	num_elms_mcc	facialdisc ~ miniche_width_sc	0.381	-0.587	1.443	0.781
168	song_to_disc	peak_frequency_mcc	facialdisc ~ miniche_width_sc	0.365	-0.517	1.371	0.766
169	song_to_disc	song_duration_mcc	facialdisc ~ miniche_width_sc	0.393	-0.450	1.335	0.807
1610	song_to_disc	song_rate_mcc	facialdisc ~ miniche_width_sc	0.407	-0.494	1.413	0.814
1611	song_to_disc	umap_space_size	facialdisc ~ miniche_width_sc	0.352	-0.526	1.318	0.774
12	disc_to_song	frequency_range_mcc	activity ~ mocover2	0.792	0.160	1.556	0.995
121	disc_to_song	num_elms_mcc	activity ~ mocover2	0.784	0.141	1.518	0.994
122	disc_to_song	peak_frequency_mcc	activity ~ mocover2	0.776	0.143	1.535	0.993
123	disc_to_song	song_duration_mcc	activity ~ mocover2	0.791	0.172	1.548	0.993
124	disc_to_song	song_rate_mcc	activity ~ mocover2	0.793	0.155	1.585	0.990
125	disc_to_song	umap_space_size	activity ~ mocover2	0.767	0.136	1.536	0.991
116	song_to_disc	frequency_range_mcc	activity ~ mocover2	0.793	0.162	1.519	0.993
117	song_to_disc	num_elms_mcc	activity ~ mocover2	0.779	0.167	1.558	0.997
118	song_to_disc	peak_frequency_mcc	activity ~ mocover2	0.779	0.148	1.567	0.992
119	song_to_disc	song_duration_mcc	activity ~ mocover2	0.798	0.165	1.566	0.993
1110	song_to_disc	song_rate_mcc	activity ~ mocover2	0.792	0.175	1.533	0.996
1111	song_to_disc	umap_space_size	activity ~ mocover2	0.777	0.141	1.541	0.992
14	disc_to_song	frequency_range_mcc	facialdisc ~ mocover2	1.276	0.262	2.378	0.995
141	disc_to_song	num_elms_mcc	facialdisc ~ mocover2	1.241	0.271	2.370	0.992
142	disc_to_song	peak_frequency_mcc	facialdisc ~ mocover2	1.261	0.326	2.452	0.995
143	disc_to_song	song_duration_mcc	facialdisc ~ mocover2	1.256	0.274	2.458	0.993
144	disc_to_song	song_rate_mcc	facialdisc ~ mocover2	1.250	0.244	2.397	0.992
145	disc_to_song	umap_space_size	facialdisc ~ mocover2	1.283	0.187	2.501	0.988
156	song_to_disc	frequency_range_mcc	facialdisc ~ mocover2	1.284	0.281	2.435	0.994
157	song_to_disc	num_elms_mcc	facialdisc ~ mocover2	1.296	0.221	2.492	0.991
158	song_to_disc	peak_frequency_mcc	facialdisc ~ mocover2	1.293	0.221	2.504	0.992
159	song_to_disc	song_duration_mcc	facialdisc ~ mocover2	1.215	0.274	2.283	0.993
1510	song_to_disc	song_rate_mcc	facialdisc ~ mocover2	1.232	0.281	2.350	0.995
1511	song_to_disc	umap_space_size	facialdisc ~ mocover2	1.294	0.270	2.448	0.991

3.1.2.1.2.2 Song features

Code

highlight(x = effects_df[grepl("frequency|song|space|elms", effects_df$predictor) |
    grepl("frequency|song|space|elms", effects_df$response), c("dag",
    "song_feature_model", "response ~ predictor", "Estimate", "Q2.5",
    "Q97.5", "pd")])

	dag	song_feature_model	response ~ predictor	Estimate	Q2.5	Q97.5	pd
8	disc_to_song	frequency_range_mcc	frequencyrangesc ~ activity	-0.067	-0.368	0.234	0.675
76	song_to_disc	frequency_range_mcc	frequencyrangesc ~ activity	-0.065	-0.368	0.226	0.669
81	disc_to_song	num_elms_mcc	numelmssc ~ activity	0.068	-0.227	0.365	0.677
77	song_to_disc	num_elms_mcc	numelmssc ~ activity	0.056	-0.228	0.331	0.655
82	disc_to_song	peak_frequency_mcc	peakfrequencysc ~ activity	-0.064	-0.349	0.233	0.666
78	song_to_disc	peak_frequency_mcc	peakfrequencysc ~ activity	-0.054	-0.334	0.228	0.634
83	disc_to_song	song_duration_mcc	songdurationsc ~ activity	0.116	-0.187	0.430	0.748
79	song_to_disc	song_duration_mcc	songdurationsc ~ activity	0.024	-0.293	0.342	0.554
84	disc_to_song	song_rate_mcc	songratesc ~ activity	-0.087	-0.398	0.208	0.710
710	song_to_disc	song_rate_mcc	songratesc ~ activity	-0.119	-0.436	0.198	0.779
85	disc_to_song	umap_space_size	umapspacesizesc ~ activity	-0.120	-0.450	0.215	0.763
711	song_to_disc	umap_space_size	umapspacesizesc ~ activity	-0.123	-0.451	0.199	0.776
9	disc_to_song	frequency_range_mcc	frequencyrangesc ~ facial_disc	0.107	-0.275	0.512	0.694
91	disc_to_song	num_elms_mcc	numelmssc ~ facial_disc	0.280	-0.150	0.711	0.904
92	disc_to_song	peak_frequency_mcc	peakfrequencysc ~ facial_disc	-0.059	-0.433	0.316	0.617
93	disc_to_song	song_duration_mcc	songdurationsc ~ facial_disc	0.349	-0.012	0.706	0.971
94	disc_to_song	song_rate_mcc	songratesc ~ facial_disc	0.299	-0.086	0.714	0.929
95	disc_to_song	umap_space_size	umapspacesizesc ~ facial_disc	0.056	-0.286	0.417	0.621
176	song_to_disc	frequency_range_mcc	facialdisc ~ mifrequency_range_sc	0.370	-0.572	1.434	0.780
17	disc_to_song	frequency_range_mcc	frequencyrangesc ~ miniche_width_sc	0.007	-0.144	0.156	0.549
146	song_to_disc	frequency_range_mcc	frequencyrangesc ~ miniche_width_sc	0.007	-0.142	0.156	0.534
171	disc_to_song	num_elms_mcc	numelmssc ~ miniche_width_sc	0.023	-0.118	0.171	0.616
147	song_to_disc	num_elms_mcc	numelmssc ~ miniche_width_sc	0.027	-0.109	0.165	0.642
172	disc_to_song	peak_frequency_mcc	peakfrequencysc ~ miniche_width_sc	0.014	-0.126	0.148	0.576
148	song_to_disc	peak_frequency_mcc	peakfrequencysc ~ miniche_width_sc	0.009	-0.128	0.137	0.563
173	disc_to_song	song_duration_mcc	songdurationsc ~ miniche_width_sc	-0.030	-0.191	0.120	0.644
149	song_to_disc	song_duration_mcc	songdurationsc ~ miniche_width_sc	-0.029	-0.194	0.133	0.650
174	disc_to_song	song_rate_mcc	songratesc ~ miniche_width_sc	-0.063	-0.215	0.084	0.788
1410	song_to_disc	song_rate_mcc	songratesc ~ miniche_width_sc	-0.060	-0.213	0.090	0.778
175	disc_to_song	umap_space_size	umapspacesizesc ~ miniche_width_sc	0.001	-0.160	0.173	0.508
1411	song_to_disc	umap_space_size	umapspacesizesc ~ miniche_width_sc	0.004	-0.152	0.168	0.512
177	song_to_disc	num_elms_mcc	facialdisc ~ minum_elms_sc	0.560	-0.305	1.550	0.903
178	song_to_disc	peak_frequency_mcc	facialdisc ~ mipeak_frequency_sc	0.054	-0.878	1.083	0.543
179	song_to_disc	song_duration_mcc	facialdisc ~ misong_duration_sc	1.107	-0.076	2.523	0.962
1710	song_to_disc	song_rate_mcc	facialdisc ~ misong_rate_sc	0.499	-0.358	1.456	0.861
1711	song_to_disc	umap_space_size	facialdisc ~ miumap_space_size_sc	0.412	-0.486	1.606	0.786
16	disc_to_song	frequency_range_mcc	frequencyrangesc ~ mocover2	-0.032	-0.247	0.181	0.616
136	song_to_disc	frequency_range_mcc	frequencyrangesc ~ mocover2	-0.023	-0.223	0.190	0.584
161	disc_to_song	num_elms_mcc	numelmssc ~ mocover2	-0.063	-0.292	0.152	0.702
137	song_to_disc	num_elms_mcc	numelmssc ~ mocover2	-0.028	-0.274	0.193	0.581
162	disc_to_song	peak_frequency_mcc	peakfrequencysc ~ mocover2	0.010	-0.199	0.213	0.556
138	song_to_disc	peak_frequency_mcc	peakfrequencysc ~ mocover2	0.000	-0.195	0.204	0.503
163	disc_to_song	song_duration_mcc	songdurationsc ~ mocover2	-0.058	-0.331	0.194	0.651
139	song_to_disc	song_duration_mcc	songdurationsc ~ mocover2	0.005	-0.259	0.234	0.549
164	disc_to_song	song_rate_mcc	songratesc ~ mocover2	-0.051	-0.275	0.160	0.681
1310	song_to_disc	song_rate_mcc	songratesc ~ mocover2	-0.024	-0.247	0.182	0.600
165	disc_to_song	umap_space_size	umapspacesizesc ~ mocover2	-0.118	-0.343	0.107	0.835
1311	song_to_disc	umap_space_size	umapspacesizesc ~ mocover2	-0.107	-0.337	0.114	0.833

3.1.2.1.2.3 Effect size heatmap

Cells show the average standardized effect size estimated across SEM formulations for each predictor–response relationship, with tile color representing the posterior probability of direction (pd), black text indicating effects whose 95% credible intervals did not overlap zero, and gray cells representing relationships that were not included in the evaluated causal models.

Code

plot_df <- effects_df

## Clean labels only

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

plot_df$response <- gsub("^mi", "", plot_df$response)
plot_df$response <- gsub("_sc$", "", plot_df$response)

## Significant rows
## TRUE if 95% CI does not overlap 0

plot_df$sig <- (
  plot_df$Q2.5 * plot_df$Q97.5
) > 0

## --------------------------------------------------
## Average repeated cells
## --------------------------------------------------

plot_df <- aggregate(
  cbind(
    Estimate,
    pd,
    sig
  ) ~ predictor + response,
  data = plot_df,
  FUN = mean,
  na.rm = TRUE
)

## Recompute significance after averaging

plot_df$sig <- plot_df$sig > 0.5

## --------------------------------------------------
## Add missing combinations
## --------------------------------------------------

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

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

## --------------------------------------------------
## Order levels
## --------------------------------------------------

predictor_order <- c(
  "facial_disc",
  "mocover2",
  "activity",
  "niche_width",
  "log_weight",
  "latitude",
  "num_elms",
  "peak_frequency",
  "song_duration",
  "song_rate",
  "umap_space_size",
  "frequency_range"
)

response_order <- c(
  "facialdisc",
  "activity",
  "nichewidthsc",
  "numelmssc",
  "peakfrequencysc",
  "songdurationsc",
  "songratesc",
  "umapspacesizesc",
  "frequencyrangesc"
)

plot_df$predictor <- factor(
  plot_df$predictor,
  levels = predictor_order
)

plot_df$response <- factor(
  plot_df$response,
  levels = response_order
)

## --------------------------------------------------
## Plot
## --------------------------------------------------

ggplot(
  plot_df,
  aes(
    x = predictor,
    y = response
  )
) +

  ## Gray background for missing combinations
  geom_tile(
    fill = "grey90",
    color = "white",
    linewidth = 0.7
  ) +

  ## Tiles for evaluated effects
  geom_tile(
    data = plot_df[!is.na(plot_df$pd), ],
    aes(fill = pd),
    color = "white",
    linewidth = 0.7
  ) +

  geom_text(
    data = plot_df[!is.na(plot_df$Estimate), ],
    aes(
      label = round(Estimate, 2),
      color = sig
    ),
    size = 3,
    fontface = "bold"
  ) +

  scale_fill_gradientn(
    colors = c(
      "#FCFFA4FF",
      "#FAC127FF",
      "#E8602DFF"
    ),
    limits = c(0.5, 1),
    name = "Posterior\nprobability\nof direction"
  ) + 
    
    scale_x_discrete(
    labels = c(
      mocover2 = "Habitat\nstructure",
      activity = "Activity\nrhythm",
      facial_disc = "Facial\ndisc",
      niche_width = "Niche\nwidth",
      log_weight = "Body\nsize",
      peak_frequency = "Peak\nfrequency",
      song_duration = "Song\nduration",
      umap_space_size = "Acoustic\nspace size",
      frequency_range = "Frequency\nrange",
       num_elms = "# of\nelements",
      song_rate = "Song rate",
      latitude = "Latitude"
    )
  ) +

  scale_y_discrete(
    labels = c(
      facialdisc = "Facial\ndisc",
      activity = "Activity\nrhythm",
      nichewidthsc = "Niche\nwidth",
      peakfrequencysc = "Peak\nfrequency",
      songdurationsc = "Song\nduration",
      umapspacesizesc = "Acoustic\nspace size",
      frequencyrangesc = "Frequency\nrange",
      numelmssc = "# of elements",
      songratesc = "Song rate"
    )
  ) +

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

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

  theme_classic(base_size = 12) +

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

    panel.grid = element_blank(),

    legend.position = "right"
  )

Takeaways

The frequency range models provided the strongest evidence favoring the song → facial disc hypothesis, with substantially better predictive performance than the reverse DAG (ΔELPD = 4.78 ± 1.86), suggesting that evolutionary changes in vocal frequency range may precede or better predict facial disc evolution.
The umap_space_size models also moderately favored the song → facial disc hypothesis (ΔELPD = 7.67 ± 7.02), although uncertainty remained relatively high.
In contrast, models based on peak frequency and song duration showed little meaningful difference between alternative causal structures, indicating that these acoustic traits are broadly compatible with either evolutionary scenario.

Overall takeaways

Comparative Bayesian phylogenetic SEMs identified consistent ecological and allometric relationships among habitat structure, activity rhythm, niche width, and facial disc morphology across owl species, particularly strong positive effects of habitat structure on both activity rhythm and facial disc traits.
Associations between facial disc morphology and song structure were generally weak to moderate and uncertain, and leave-one-out cross-validation provided little evidence favoring either the “song → facial disc” or “facial disc → song” causal hypothesis, suggesting that the available comparative data do not clearly distinguish between these alternative evolutionary scenarios.

Session information

─ 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-05-20
 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)
 boot             1.3-32    2025-08-29 [4] CRAN (R 4.5.1)
 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)
 cellranger       1.1.0     2016-07-27 [3] CRAN (R 4.0.1)
 checkmate        2.3.4     2026-02-03 [1] CRAN (R 4.5.2)
 cli              3.6.6     2026-04-09 [1] CRAN (R 4.5.2)
 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)
 corrplot       * 0.95      2024-10-14 [1] 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)
 dagitty        * 0.3-4     2023-12-07 [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)
 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)
 fastmatch        1.1-6     2024-12-23 [3] CRAN (R 4.5.0)
 formatR          1.14      2023-01-17 [1] CRAN (R 4.5.0)
 Formula          1.2-5     2023-02-24 [3] 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)
 ggdag          * 0.2.13    2024-07-22 [1] CRAN (R 4.5.0)
 ggdist           3.3.3     2025-04-23 [1] CRAN (R 4.5.0)
 ggforce          0.3.3     2021-03-05 [3] CRAN (R 4.1.1)
 ggparty        * 1.0.0.1   2025-07-10 [1] CRAN (R 4.5.2)
 ggplot2        * 4.0.3     2026-04-22 [1] CRAN (R 4.5.2)
 ggraph         * 2.0.5     2021-02-23 [3] CRAN (R 4.1.1)
 ggrepel          0.9.6     2024-09-07 [1] CRAN (R 4.5.0)
 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)
 graphlayouts     0.8.0     2022-01-03 [3] CRAN (R 4.1.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)
 inum             1.0-5     2023-03-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)
 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)
 libcoin        * 1.0-10    2023-09-27 [1] CRAN (R 4.5.0)
 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)
 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)
 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)
 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)
 partykit       * 1.2-24    2025-05-02 [1] CRAN (R 4.5.0)
 pbapply          1.7-4     2025-07-20 [1] CRAN (R 4.5.0)
 phangorn       * 2.12.1    2024-09-17 [1] CRAN (R 4.5.0)
 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)
 prettyunits      1.2.0     2023-09-24 [3] 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)
 progress         1.2.3     2023-12-06 [3] CRAN (R 4.5.0)
 promises         1.3.3     2025-05-29 [1] RSPM (R 4.5.0)
 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)
 quadprog         1.5-8     2019-11-20 [3] CRAN (R 4.0.1)
 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)
 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)
 readxl           1.4.5     2025-03-07 [3] 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)
 rncl           * 0.8.9     2026-01-21 [1] CRAN (R 4.5.2)
 rpart            4.1.24    2025-01-07 [4] CRAN (R 4.4.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)
 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)
 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)
 tensorA          0.36.2.1  2023-12-13 [1] CRAN (R 4.5.0)
 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)
 tidygraph        1.2.0     2020-05-12 [3] CRAN (R 4.0.1)
 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)
 tweenr           2.0.3     2024-02-26 [3] CRAN (R 4.5.0)
 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)
 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.

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