Statistical analysis

Owl facial disc evolution

Author

Published

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

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

    X$`l-95% CI` <- X$Q2.5
    X$`u-95% CI` <- X$Q97.5
    X$Q2.5 <- X$Q97.5 <- NULL
    out <- brmsish:::html_format_coef_table(X, fill = "#A0DFB980",
        highlight = TRUE)
    print(out)
}

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]

Takeaway

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

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

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

  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 Results

3.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

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

    ## 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 = TRUE, 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
sem_song_to_disc	song_pc1_freq_sc	sem_song_to_disc	0.000	0.000	2340.481
sem_disc_to_song	song_pc1_freq_sc	sem_disc_to_song	-0.609	1.502	2341.698
sem_song_to_disc1	song_pc2_diver_sc	sem_song_to_disc	0.000	0.000	2436.150
sem_disc_to_song1	song_pc2_diver_sc	sem_disc_to_song	-2.582	3.188	2441.314
sem_song_to_disc2	song_pc4_length_sc	sem_song_to_disc	0.000	0.000	2436.424
sem_disc_to_song2	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.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

Code

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

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


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

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

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

    ## Cleanup

    rm(fit, fe)

    gc()
}

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


## Final table
highlight(effects_df)

	dag	response	predictor	song_feature_model	Estimate	Est.Error	l-95% CI	u-95% CI
13	disc_to_song	activity	milog_weight_sc	song_pc1_freq	0.016	0.318	-0.619	0.604
123	song_to_disc	activity	milog_weight_sc	song_pc1_freq	0.023	0.326	-0.649	0.647
131	disc_to_song	activity	milog_weight_sc	song_pc2_diver	0.013	0.330	-0.641	0.730
124	song_to_disc	activity	milog_weight_sc	song_pc2_diver	0.008	0.333	-0.688	0.633
132	disc_to_song	activity	milog_weight_sc	song_pc4_length	0.006	0.323	-0.662	0.642
125	song_to_disc	activity	milog_weight_sc	song_pc4_length	0.005	0.335	-0.649	0.639
12	disc_to_song	activity	mocover2	song_pc1_freq	0.780	0.361	0.123	1.511
113	song_to_disc	activity	mocover2	song_pc1_freq	0.794	0.360	0.178	1.569
121	disc_to_song	activity	mocover2	song_pc2_diver	0.770	0.353	0.183	1.508
114	song_to_disc	activity	mocover2	song_pc2_diver	0.780	0.349	0.183	1.491
122	disc_to_song	activity	mocover2	song_pc4_length	0.782	0.361	0.185	1.569
115	song_to_disc	activity	mocover2	song_pc4_length	0.793	0.362	0.192	1.556
7	disc_to_song	facialdisc	activity	song_pc1_freq	-0.452	0.668	-1.774	0.872
83	song_to_disc	facialdisc	activity	song_pc1_freq	-0.416	0.692	-1.782	0.972
71	disc_to_song	facialdisc	activity	song_pc2_diver	-0.403	0.673	-1.743	0.927
84	song_to_disc	facialdisc	activity	song_pc2_diver	-0.377	0.700	-1.799	0.939
72	disc_to_song	facialdisc	activity	song_pc4_length	-0.415	0.711	-1.739	1.052
85	song_to_disc	facialdisc	activity	song_pc4_length	-0.344	0.732	-1.807	1.071
15	disc_to_song	facialdisc	miniche_width_sc	song_pc1_freq	0.349	0.480	-0.585	1.270
163	song_to_disc	facialdisc	miniche_width_sc	song_pc1_freq	0.384	0.478	-0.525	1.375
151	disc_to_song	facialdisc	miniche_width_sc	song_pc2_diver	0.347	0.505	-0.562	1.482
164	song_to_disc	facialdisc	miniche_width_sc	song_pc2_diver	0.358	0.482	-0.562	1.355
152	disc_to_song	facialdisc	miniche_width_sc	song_pc4_length	0.331	0.474	-0.580	1.302
165	song_to_disc	facialdisc	miniche_width_sc	song_pc4_length	0.366	0.496	-0.599	1.366
173	song_to_disc	facialdisc	misong_pc1_freq_sc	song_pc1_freq	0.417	0.540	-0.564	1.470
174	song_to_disc	facialdisc	misong_pc2_diver_sc	song_pc2_diver	0.889	0.501	-0.078	1.887
175	song_to_disc	facialdisc	misong_pc4_length_sc	song_pc4_length	-0.725	0.504	-1.854	0.183
14	disc_to_song	facialdisc	mocover2	song_pc1_freq	1.253	0.540	0.273	2.330
153	song_to_disc	facialdisc	mocover2	song_pc1_freq	1.303	0.559	0.320	2.501
141	disc_to_song	facialdisc	mocover2	song_pc2_diver	1.241	0.504	0.335	2.287
154	song_to_disc	facialdisc	mocover2	song_pc2_diver	1.251	0.549	0.276	2.399
142	disc_to_song	facialdisc	mocover2	song_pc4_length	1.270	0.578	0.185	2.560
155	song_to_disc	facialdisc	mocover2	song_pc4_length	1.211	0.572	0.140	2.336
10	disc_to_song	nichewidthsc	milatitude_sc	song_pc1_freq	0.013	0.070	-0.113	0.148
93	song_to_disc	nichewidthsc	milatitude_sc	song_pc1_freq	0.011	0.067	-0.113	0.140
101	disc_to_song	nichewidthsc	milatitude_sc	song_pc2_diver	0.012	0.067	-0.121	0.143
94	song_to_disc	nichewidthsc	milatitude_sc	song_pc2_diver	0.014	0.066	-0.110	0.140
102	disc_to_song	nichewidthsc	milatitude_sc	song_pc4_length	0.009	0.068	-0.130	0.155
95	song_to_disc	nichewidthsc	milatitude_sc	song_pc4_length	0.009	0.071	-0.127	0.150
11	disc_to_song	nichewidthsc	milog_weight_sc	song_pc1_freq	0.319	0.105	0.111	0.513
103	song_to_disc	nichewidthsc	milog_weight_sc	song_pc1_freq	0.316	0.098	0.126	0.507
111	disc_to_song	nichewidthsc	milog_weight_sc	song_pc2_diver	0.309	0.102	0.117	0.499
104	song_to_disc	nichewidthsc	milog_weight_sc	song_pc2_diver	0.309	0.102	0.112	0.522
112	disc_to_song	nichewidthsc	milog_weight_sc	song_pc4_length	0.314	0.102	0.119	0.523
105	song_to_disc	nichewidthsc	milog_weight_sc	song_pc4_length	0.314	0.103	0.118	0.529
8	disc_to_song	songpc1freqsc	activity	song_pc1_freq	0.096	0.132	-0.159	0.353
73	song_to_disc	songpc1freqsc	activity	song_pc1_freq	0.073	0.132	-0.191	0.323
9	disc_to_song	songpc1freqsc	facial_disc	song_pc1_freq	0.236	0.192	-0.131	0.609
17	disc_to_song	songpc1freqsc	miniche_width_sc	song_pc1_freq	0.006	0.068	-0.129	0.139
143	song_to_disc	songpc1freqsc	miniche_width_sc	song_pc1_freq	0.005	0.068	-0.120	0.142
16	disc_to_song	songpc1freqsc	mocover2	song_pc1_freq	0.032	0.103	-0.165	0.237
133	song_to_disc	songpc1freqsc	mocover2	song_pc1_freq	0.058	0.094	-0.116	0.244
81	disc_to_song	songpc2diversc	activity	song_pc2_diver	-0.117	0.161	-0.437	0.194
74	song_to_disc	songpc2diversc	activity	song_pc2_diver	-0.174	0.153	-0.469	0.127
91	disc_to_song	songpc2diversc	facial_disc	song_pc2_diver	0.348	0.203	-0.068	0.764
171	disc_to_song	songpc2diversc	miniche_width_sc	song_pc2_diver	0.023	0.083	-0.127	0.188
144	song_to_disc	songpc2diversc	miniche_width_sc	song_pc2_diver	0.036	0.080	-0.117	0.187
161	disc_to_song	songpc2diversc	mocover2	song_pc2_diver	-0.078	0.122	-0.326	0.167
134	song_to_disc	songpc2diversc	mocover2	song_pc2_diver	-0.033	0.115	-0.265	0.181
82	disc_to_song	songpc4lengthsc	activity	song_pc4_length	0.065	0.160	-0.242	0.359
75	song_to_disc	songpc4lengthsc	activity	song_pc4_length	0.129	0.173	-0.202	0.470
92	disc_to_song	songpc4lengthsc	facial_disc	song_pc4_length	-0.309	0.198	-0.692	0.087
172	disc_to_song	songpc4lengthsc	miniche_width_sc	song_pc4_length	0.015	0.084	-0.139	0.171
145	song_to_disc	songpc4lengthsc	miniche_width_sc	song_pc4_length	0.019	0.086	-0.146	0.196
162	disc_to_song	songpc4lengthsc	mocover2	song_pc4_length	-0.118	0.117	-0.359	0.116
135	song_to_disc	songpc4lengthsc	mocover2	song_pc4_length	-0.166	0.108	-0.382	0.044

Takeaways

Across both DAG formulations, habitat structure showed consistent positive associations with activity rhythm and facial disc morphology, while body size was positively associated with niche width, suggesting robust ecological and allometric effects across owl species.
Effects linking facial disc morphology and song structure were generally moderate in magnitude but uncertain, and their direction and strength were highly consistent across alternative causal models, indicating that the comparative data provide limited support for distinguishing between competing evolutionary hypotheses.

3.1.1.3 Individual model fits

3.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.3.2 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

3.1.1.4 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.4.1 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

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-14
 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)
 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.

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