Statistical analysis v2

Owl facial disc evolution

Author

Published

June 17, 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$abs_latitude_sc <- scale(abs(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,
    facial_disc ~ body_size,
    facial_disc ~ latitude,
    
    # 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} \\ \\ \text{activity rhythm} &\sim \text{monotonic(habitat structure)} + \text{body size} \\ \text{facial disc} &\sim \text{monotonic(habitat structure)} + \\ &\quad \text{activity rhythm} + \text{niche width} + \text{body size} + \text{latitude} \\ \\ \text{song feature} &\sim \text{monotonic(habitat structure)} + \\ &\quad \text{activity rhythm} + \text{niche width} + \\ &\quad \text{body size} + \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(abs_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(
    abs_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 +
    abs_latitude_sc +
    mi(niche_width_sc) +
    mi(log_weight_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) + ",
      "mi(log_weight_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
9	frequency_range	facialdisc ~ abs_latitude_sc	0.502	-0.414	1.511	0.862
91	num_elms	facialdisc ~ abs_latitude_sc	0.496	-0.391	1.522	0.863
92	peak_frequency	facialdisc ~ abs_latitude_sc	0.511	-0.400	1.464	0.846
93	song_duration	facialdisc ~ abs_latitude_sc	0.534	-0.365	1.548	0.866
94	song_rate	facialdisc ~ abs_latitude_sc	0.501	-0.372	1.543	0.865
95	umap_space_size	facialdisc ~ abs_latitude_sc	0.506	-0.363	1.504	0.869
8	frequency_range	facialdisc ~ activity	-0.410	-1.797	1.088	0.727
81	num_elms	facialdisc ~ activity	-0.417	-1.837	1.087	0.711
82	peak_frequency	facialdisc ~ activity	-0.438	-1.910	1.016	0.715
83	song_duration	facialdisc ~ activity	-0.430	-1.952	0.952	0.727
84	song_rate	facialdisc ~ activity	-0.414	-1.789	1.075	0.725
85	umap_space_size	facialdisc ~ activity	-0.422	-1.916	1.062	0.708
7	frequency_range	nichewidthsc ~ activity	0.271	-0.027	0.569	0.963
71	num_elms	nichewidthsc ~ activity	0.272	-0.016	0.565	0.964
72	peak_frequency	nichewidthsc ~ activity	0.270	0.004	0.552	0.976
73	song_duration	nichewidthsc ~ activity	0.271	-0.012	0.551	0.968
74	song_rate	nichewidthsc ~ activity	0.268	-0.013	0.550	0.968
75	umap_space_size	nichewidthsc ~ activity	0.273	-0.022	0.575	0.964
12	frequency_range	nichewidthsc ~ miabs_latitude_sc	0.035	-0.108	0.175	0.702
121	num_elms	nichewidthsc ~ miabs_latitude_sc	0.037	-0.105	0.174	0.709
122	peak_frequency	nichewidthsc ~ miabs_latitude_sc	0.036	-0.100	0.174	0.692
123	song_duration	nichewidthsc ~ miabs_latitude_sc	0.039	-0.091	0.171	0.718
124	song_rate	nichewidthsc ~ miabs_latitude_sc	0.038	-0.099	0.180	0.707
125	umap_space_size	nichewidthsc ~ miabs_latitude_sc	0.037	-0.094	0.175	0.704
15	frequency_range	activity ~ milog_weight_sc	0.028	-0.657	0.667	0.549
151	num_elms	activity ~ milog_weight_sc	0.016	-0.595	0.660	0.507
152	peak_frequency	activity ~ milog_weight_sc	0.017	-0.609	0.641	0.519
153	song_duration	activity ~ milog_weight_sc	-0.003	-0.620	0.605	0.501
154	song_rate	activity ~ milog_weight_sc	0.019	-0.593	0.645	0.514
155	umap_space_size	activity ~ milog_weight_sc	0.019	-0.637	0.652	0.531
18	frequency_range	facialdisc ~ milog_weight_sc	-0.129	-1.473	1.214	0.574
181	num_elms	facialdisc ~ milog_weight_sc	-0.131	-1.514	1.187	0.589
182	peak_frequency	facialdisc ~ milog_weight_sc	-0.133	-1.502	1.135	0.562
183	song_duration	facialdisc ~ milog_weight_sc	-0.102	-1.450	1.271	0.556
184	song_rate	facialdisc ~ milog_weight_sc	-0.109	-1.499	1.222	0.551
185	umap_space_size	facialdisc ~ milog_weight_sc	-0.145	-1.520	1.194	0.569
13	frequency_range	nichewidthsc ~ milog_weight_sc	0.302	0.102	0.496	0.997
131	num_elms	nichewidthsc ~ milog_weight_sc	0.302	0.111	0.502	1.000
132	peak_frequency	nichewidthsc ~ milog_weight_sc	0.303	0.101	0.508	0.998
133	song_duration	nichewidthsc ~ milog_weight_sc	0.300	0.108	0.495	0.998
134	song_rate	nichewidthsc ~ milog_weight_sc	0.303	0.108	0.504	0.998
135	umap_space_size	nichewidthsc ~ milog_weight_sc	0.301	0.100	0.510	0.999
17	frequency_range	facialdisc ~ miniche_width_sc	0.453	-0.517	1.512	0.818
171	num_elms	facialdisc ~ miniche_width_sc	0.452	-0.507	1.479	0.816
172	peak_frequency	facialdisc ~ miniche_width_sc	0.450	-0.542	1.490	0.823
173	song_duration	facialdisc ~ miniche_width_sc	0.462	-0.473	1.468	0.830
174	song_rate	facialdisc ~ miniche_width_sc	0.447	-0.513	1.492	0.814
175	umap_space_size	facialdisc ~ miniche_width_sc	0.440	-0.499	1.510	0.810
14	frequency_range	activity ~ mocover2	0.783	0.174	1.583	0.997
141	num_elms	activity ~ mocover2	0.787	0.154	1.540	0.995
142	peak_frequency	activity ~ mocover2	0.793	0.154	1.575	0.993
143	song_duration	activity ~ mocover2	0.776	0.145	1.489	0.992
144	song_rate	activity ~ mocover2	0.779	0.179	1.561	0.997
145	umap_space_size	activity ~ mocover2	0.784	0.145	1.546	0.992
16	frequency_range	facialdisc ~ mocover2	1.247	0.244	2.492	0.990
161	num_elms	facialdisc ~ mocover2	1.251	0.173	2.523	0.986
162	peak_frequency	facialdisc ~ mocover2	1.228	0.114	2.406	0.986
163	song_duration	facialdisc ~ mocover2	1.193	0.090	2.376	0.984
164	song_rate	facialdisc ~ mocover2	1.233	0.142	2.450	0.985
165	umap_space_size	facialdisc ~ mocover2	1.234	0.152	2.388	0.988

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
11	frequency_range	frequencyrangesc ~ activity	-0.063	-0.361	0.245	0.665
111	num_elms	numelmssc ~ activity	0.070	-0.225	0.365	0.685
112	peak_frequency	peakfrequencysc ~ activity	-0.055	-0.333	0.224	0.658
113	song_duration	songdurationsc ~ activity	0.114	-0.195	0.430	0.755
114	song_rate	songratesc ~ activity	-0.088	-0.403	0.226	0.700
115	umap_space_size	umapspacesizesc ~ activity	-0.126	-0.453	0.208	0.775
10	frequency_range	frequencyrangesc ~ facial_disc	0.105	-0.286	0.502	0.694
101	num_elms	numelmssc ~ facial_disc	0.293	-0.118	0.702	0.913
102	peak_frequency	peakfrequencysc ~ facial_disc	-0.054	-0.444	0.324	0.629
103	song_duration	songdurationsc ~ facial_disc	0.360	0.030	0.707	0.986
104	song_rate	songratesc ~ facial_disc	0.301	-0.100	0.696	0.926
105	umap_space_size	umapspacesizesc ~ facial_disc	0.105	-0.259	0.460	0.716
21	frequency_range	frequencyrangesc ~ milog_weight_sc	0.084	-0.140	0.325	0.766
211	num_elms	numelmssc ~ milog_weight_sc	-0.050	-0.312	0.195	0.643
212	peak_frequency	peakfrequencysc ~ milog_weight_sc	-0.120	-0.337	0.094	0.849
213	song_duration	songdurationsc ~ milog_weight_sc	-0.042	-0.216	0.122	0.693
214	song_rate	songratesc ~ milog_weight_sc	-0.023	-0.237	0.192	0.582
215	umap_space_size	umapspacesizesc ~ milog_weight_sc	-0.171	-0.362	0.009	0.971
20	frequency_range	frequencyrangesc ~ miniche_width_sc	-0.006	-0.154	0.145	0.534
201	num_elms	numelmssc ~ miniche_width_sc	0.030	-0.114	0.181	0.660
202	peak_frequency	peakfrequencysc ~ miniche_width_sc	0.031	-0.112	0.169	0.661
203	song_duration	songdurationsc ~ miniche_width_sc	-0.021	-0.181	0.137	0.606
204	song_rate	songratesc ~ miniche_width_sc	-0.060	-0.224	0.095	0.768
205	umap_space_size	umapspacesizesc ~ miniche_width_sc	0.048	-0.117	0.223	0.705
19	frequency_range	frequencyrangesc ~ mocover2	-0.029	-0.247	0.182	0.606
191	num_elms	numelmssc ~ mocover2	-0.071	-0.310	0.137	0.714
192	peak_frequency	peakfrequencysc ~ mocover2	0.000	-0.222	0.200	0.506
193	song_duration	songdurationsc ~ mocover2	-0.058	-0.332	0.191	0.653
194	song_rate	songratesc ~ mocover2	-0.050	-0.263	0.173	0.667
195	umap_space_size	umapspacesizesc ~ mocover2	-0.130	-0.361	0.087	0.873

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",
  "abs_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",
      abs_latitude = "Absolute\nlatitude"
    )
  ) +

  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 was the strongest and most consistent predictor in the SEM, showing positive associations with both activity rhythm and facial disc morphology across all models, supporting a major role of habitat complexity in shaping owl sensory ecology.
Body size was consistently associated with broader niche width, but showed little evidence of direct effects on facial disc morphology or most acoustic traits after accounting for ecological variables, suggesting that allometric effects are primarily mediated through ecological pathways.
Facial disc morphology exhibited a clear positive association with song duration and weaker positive associations with song rate and number of song elements. These results suggest that species with more developed facial discs tend to produce longer and potentially more structurally elaborate songs, although support was not equally strong across all acoustic traits.
Ecological predictors explained relatively little variation in song structure overall. Most acoustic traits showed weak and inconsistent relationships with habitat structure, activity rhythm, niche width, and body size, indicating that the factors driving song evolution may differ from those shaping facial disc morphology.

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

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