Statistical analysis v2

Owl facial disc evolution

Author

Published

June 12, 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 DAG) 
  B --> C("Fit phylogenetic\nSEM")
  C --> D(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 Read and format data

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

2 DAG

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

3 Fit univariate phylogenetic SEM models

A phylogenetic structural equation model (SEMs) was fitted to evaluate competing hypotheses regarding the causal relationship between song traits and facial disc morphology in owls

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}\]

]

[ \[\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

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

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 acoustic variables
## ======================================================

for (song_var in features) {

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

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

  ## --------------------------------------------------
  ## Shared submodels
  ## --------------------------------------------------

  bf_niche <- bf(
    niche_width_sc | mi() ~
      mi(latitude_sc) +
      mi(log_weight_sc) +
      activity +
      (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
    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
    prior(
      normal(0, 1),
      class = "b",
      resp = "activity"
    ),

    prior(
      normal(0, 1.5),
      class = "Intercept",
      resp = "activity"
    ),

    ## Song trait
    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
    prior(
      normal(0, 1),
      class = "b",
      resp = "facialdisc"
    ),

    prior(
      normal(0, 1.5),
      class = "Intercept",
      resp = "facialdisc"
    )
  )

  ## --------------------------------------------------
  ## Facial disc model
  ## --------------------------------------------------

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

  ## --------------------------------------------------
  ## Song model
  ## Facial disc -> song
  ## --------------------------------------------------

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

  ## --------------------------------------------------
  ## Model name
  ## --------------------------------------------------

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

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

  ## --------------------------------------------------
  ## Fit model
  ## --------------------------------------------------

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

  ## --------------------------------------------------
  ## Complete cases for LOO
  ## --------------------------------------------------

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

  ## --------------------------------------------------
  ## Add LOO
  ## --------------------------------------------------

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

  ## --------------------------------------------------
  ## Save model
  ## --------------------------------------------------

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

  rm(fit)

  gc()
}

3.1.0.1 Results

3.1.0.1.1 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.0.1.1.1 Ecological predictors

Code

## ====================================================== Empty
## dataframe
## ======================================================

effects_df <- data.frame()

## ====================================================== Models
## ======================================================

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

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

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

## ======================================================
## Extract effects
## ======================================================

for (fl in fls) {

    fit <- readRDS(fl)

    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)

    ## ------------------------------------------ Acoustic
    ## feature ------------------------------------------

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

    ## ------------------------------------------ Posterior
    ## probability of direction
    ## ------------------------------------------

    draws <- posterior::as_draws_df(fit)

    fe$pd <- NA_real_

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

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

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

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

        samples <- draws[[param]]

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

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

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

    effects_df <- rbind(effects_df, fe)

    rm(fit, fe, draws)

    gc()
}

## ======================================================
## Formatting
## ======================================================

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

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

## ======================================================
## Ecological effects
## ======================================================

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

	song_feature_model	response ~ predictor	Estimate	Q2.5	Q97.5	pd
8	frequency_range	facialdisc ~ activity	-0.422	-1.741	0.852	0.734
81	num_elms	facialdisc ~ activity	-0.414	-1.824	0.942	0.726
82	peak_frequency	facialdisc ~ activity	-0.386	-1.690	0.998	0.715
83	song_duration	facialdisc ~ activity	-0.396	-1.837	0.997	0.713
84	song_rate	facialdisc ~ activity	-0.436	-1.788	0.914	0.739
85	umap_space_size	facialdisc ~ activity	-0.415	-1.773	0.957	0.726
7	frequency_range	nichewidthsc ~ activity	0.276	-0.008	0.572	0.970
71	num_elms	nichewidthsc ~ activity	0.284	-0.015	0.584	0.971
72	peak_frequency	nichewidthsc ~ activity	0.280	-0.028	0.572	0.966
73	song_duration	nichewidthsc ~ activity	0.286	-0.010	0.572	0.972
74	song_rate	nichewidthsc ~ activity	0.283	-0.008	0.591	0.970
75	umap_space_size	nichewidthsc ~ activity	0.285	0.010	0.562	0.981
11	frequency_range	nichewidthsc ~ milatitude_sc	0.003	-0.136	0.142	0.507
111	num_elms	nichewidthsc ~ milatitude_sc	0.004	-0.137	0.142	0.518
112	peak_frequency	nichewidthsc ~ milatitude_sc	0.006	-0.122	0.132	0.534
113	song_duration	nichewidthsc ~ milatitude_sc	0.003	-0.137	0.137	0.514
114	song_rate	nichewidthsc ~ milatitude_sc	0.005	-0.132	0.142	0.521
115	umap_space_size	nichewidthsc ~ milatitude_sc	0.002	-0.130	0.131	0.521
14	frequency_range	activity ~ milog_weight_sc	0.010	-0.639	0.683	0.501
141	num_elms	activity ~ milog_weight_sc	0.015	-0.627	0.674	0.514
142	peak_frequency	activity ~ milog_weight_sc	0.015	-0.619	0.654	0.518
143	song_duration	activity ~ milog_weight_sc	0.036	-0.614	0.709	0.541
144	song_rate	activity ~ milog_weight_sc	0.023	-0.704	0.694	0.525
145	umap_space_size	activity ~ milog_weight_sc	0.033	-0.612	0.723	0.531
12	frequency_range	nichewidthsc ~ milog_weight_sc	0.309	0.111	0.509	1.000
121	num_elms	nichewidthsc ~ milog_weight_sc	0.309	0.118	0.505	0.999
122	peak_frequency	nichewidthsc ~ milog_weight_sc	0.308	0.121	0.501	1.000
123	song_duration	nichewidthsc ~ milog_weight_sc	0.309	0.114	0.505	1.000
124	song_rate	nichewidthsc ~ milog_weight_sc	0.309	0.111	0.494	0.998
125	umap_space_size	nichewidthsc ~ milog_weight_sc	0.313	0.133	0.505	1.000
16	frequency_range	facialdisc ~ miniche_width_sc	0.355	-0.541	1.372	0.769
161	num_elms	facialdisc ~ miniche_width_sc	0.361	-0.547	1.395	0.781
162	peak_frequency	facialdisc ~ miniche_width_sc	0.365	-0.539	1.343	0.782
163	song_duration	facialdisc ~ miniche_width_sc	0.354	-0.545	1.355	0.777
164	song_rate	facialdisc ~ miniche_width_sc	0.369	-0.546	1.355	0.771
165	umap_space_size	facialdisc ~ miniche_width_sc	0.362	-0.498	1.365	0.781
13	frequency_range	activity ~ mocover2	0.790	0.147	1.531	0.994
131	num_elms	activity ~ mocover2	0.801	0.131	1.589	0.994
132	peak_frequency	activity ~ mocover2	0.773	0.117	1.544	0.990
133	song_duration	activity ~ mocover2	0.789	0.169	1.560	0.996
134	song_rate	activity ~ mocover2	0.794	0.174	1.531	0.995
135	umap_space_size	activity ~ mocover2	0.773	0.126	1.518	0.992
15	frequency_range	facialdisc ~ mocover2	1.266	0.248	2.426	0.990
151	num_elms	facialdisc ~ mocover2	1.261	0.265	2.425	0.993
152	peak_frequency	facialdisc ~ mocover2	1.278	0.264	2.426	0.993
153	song_duration	facialdisc ~ mocover2	1.275	0.264	2.446	0.993
154	song_rate	facialdisc ~ mocover2	1.273	0.269	2.424	0.992
155	umap_space_size	facialdisc ~ mocover2	1.224	0.188	2.335	0.993

3.1.0.1.1.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("song_feature_model",
    "response ~ predictor", "Estimate", "Q2.5", "Q97.5", "pd")])

	song_feature_model	response ~ predictor	Estimate	Q2.5	Q97.5	pd
10	frequency_range	frequencyrangesc ~ activity	-0.058	-0.360	0.243	0.656
101	num_elms	numelmssc ~ activity	0.073	-0.229	0.362	0.685
102	peak_frequency	peakfrequencysc ~ activity	-0.058	-0.347	0.220	0.656
103	song_duration	songdurationsc ~ activity	0.109	-0.227	0.438	0.743
104	song_rate	songratesc ~ activity	-0.091	-0.432	0.218	0.718
105	umap_space_size	umapspacesizesc ~ activity	-0.120	-0.458	0.209	0.766
9	frequency_range	frequencyrangesc ~ facial_disc	0.115	-0.273	0.525	0.719
91	num_elms	numelmssc ~ facial_disc	0.274	-0.146	0.687	0.904
92	peak_frequency	peakfrequencysc ~ facial_disc	-0.057	-0.429	0.343	0.633
93	song_duration	songdurationsc ~ facial_disc	0.343	-0.004	0.704	0.974
94	song_rate	songratesc ~ facial_disc	0.287	-0.103	0.686	0.921
95	umap_space_size	umapspacesizesc ~ facial_disc	0.048	-0.317	0.412	0.600
18	frequency_range	frequencyrangesc ~ miniche_width_sc	0.007	-0.141	0.162	0.529
181	num_elms	numelmssc ~ miniche_width_sc	0.020	-0.128	0.161	0.614
182	peak_frequency	peakfrequencysc ~ miniche_width_sc	0.012	-0.126	0.153	0.562
183	song_duration	songdurationsc ~ miniche_width_sc	-0.029	-0.190	0.142	0.653
184	song_rate	songratesc ~ miniche_width_sc	-0.066	-0.223	0.084	0.798
185	umap_space_size	umapspacesizesc ~ miniche_width_sc	0.004	-0.161	0.163	0.512
17	frequency_range	frequencyrangesc ~ mocover2	-0.032	-0.241	0.169	0.620
171	num_elms	numelmssc ~ mocover2	-0.055	-0.292	0.154	0.681
172	peak_frequency	peakfrequencysc ~ mocover2	0.010	-0.203	0.219	0.549
173	song_duration	songdurationsc ~ mocover2	-0.054	-0.326	0.188	0.622
174	song_rate	songratesc ~ mocover2	-0.047	-0.267	0.166	0.667
175	umap_space_size	umapspacesizesc ~ mocover2	-0.119	-0.351	0.100	0.852

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

Habitat structure showed the strongest and most consistent associations across models, with positive effects on both activity rhythm and facial disc morphology, indicating that habitat complexity is a major predictor of sensory trait evolution in owls.
Body size was consistently and strongly associated with niche width, whereas latitude showed little evidence of direct effects on any focal trait.
Associations between facial disc morphology and acoustic traits were generally weak to moderate. The strongest positive associations were observed for song duration, song rate, and number of song elements, whereas frequency-based traits showed little evidence of association with facial disc morphology.
Activity rhythm and niche width exhibited only weak and inconsistent relationships with acoustic traits, suggesting that most song variation is not strongly explained by the ecological variables included in the SEM.

:::

Overall takeaways

Comparative Bayesian phylogenetic SEMs identified strong and highly consistent ecological drivers of facial disc evolution across owl species. Habitat structure was positively associated with both activity rhythm and facial disc morphology, while body size was strongly associated with niche width, indicating robust ecological and allometric influences on sensory trait evolution.
In contrast, associations between facial disc morphology and song structure were generally weaker and more uncertain. Although some acoustic traits (particularly song duration, song rate, and number of song elements) showed moderate positive associations with facial disc morphology, effect sizes were small relative to those associated with ecological predictors.
The overall pattern suggests that ecological conditions provide a substantially better explanation for variation in facial disc morphology than song traits themselves. Rather than supporting a strong causal role of vocal evolution in driving facial disc evolution, the results are more consistent with facial discs and song traits exhibiting modest correlated evolution within a broader ecological context.

:::

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-06-12
 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)
 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)
 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)
 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)
 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)
 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.

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