Degradation statistics

baRulho:quantifying habitat-induced degradation of (animal) acoustic signals

Author

Marcelo Araya-Salas

Published

October 10, 2024

Source code, data and annotation protocol found at https://github.com/maRce10/barulho_paper

Purposes

  • Summarize degradation metrics on re-recorded signals from playback experiment at Bosque de Tlalpan, Mexico City, 2019

  • Run statistical analyses

 

Load packages

Code 
# install/ load packages
sketchy::load_packages(packages = c("remotes", "kableExtra", "knitr",
    "formatR", "rprojroot", "maRce10/warbleR", "ggplot2", "tidyr",
    "viridis", "corrplot", "brms", "ggdist", "cowplot", "cmdstanr",
    "maRce10/brmsish", "emmeans", "ggsignif", "rstan", "loo"))

Custom functions

Code 
my.viridis <- function(...) mako(alpha = 0.5, begin = 0.3, end = 0.7,
    ...)

# to create several posterior predictive check plots out of a
# brms fit
custom_ppc1 <- function(fit, group = NULL, ndraws = 30) {

    ppc_dens <- pp_check(fit, ndraws = ndraws, type = "dens_overlay_grouped",
        group = group)  # shows dens_overlay plot by default
    pp_mean <- pp_check(fit, type = "stat_grouped", stat = "mean",
        group = group, ndraws = ndraws)
    pp_scat <- pp_check(fit, type = "scatter_avg_grouped", group = group,
        ndraws = ndraws)
    pp_stat2 <- pp_check(fit, type = "stat_2d", ndraws = ndraws)
    pp_loo <- pp_check(fit, type = "loo_pit_qq", ndraws = ndraws)
    pp_error <- pp_check(fit, type = "error_scatter_avg_grouped",
        ndraws = ndraws, group = group)
    plot_l <- list(ppc_dens, pp_mean, pp_scat, pp_error, pp_stat2,
        pp_loo)

    plot_l[c(1, 3:length(plot_l))] <- lapply(plot_l[c(1, 3:length(plot_l))],
        function(x) x + scale_color_viridis_d(begin = 0.1, end = 0.8,
            alpha = 0.5) + scale_fill_viridis_d(begin = 0.1, end = 0.8,
            alpha = 0.5) + theme_classic())

    print(plot_grid(plotlist = plot_l, ncol = 2))
}

custom_ppc <- function(...) suppressMessages(custom_ppc1(...))
Code 
degrad_df <- read.csv("./data/processed/tlalpan_degradation_metrics_v01.csv")

degrad_df <- degrad_df[grep("marker", degrad_df$sound.id, invert = TRUE),
    ]

degrad_params <- c("blur.ratio", "spectrum.blur.ratio", "envelope.correlation",
    "excess.attenuation", "signal.to.noise.ratio", "cross.correlation",
    "tail.to.signal.ratio", "spectrum.correlation")

comp.cases <- complete.cases(degrad_df[, names(degrad_df) %in% degrad_params])

pca <- prcomp(degrad_df[comp.cases, names(degrad_df) %in% degrad_params],
    scale. = TRUE)

# add to data
degrad_df$pc1 <- NA
degrad_df$pc1[comp.cases] <- pca$x[, 1]
degrad_df$pc1.1m.rate <- degrad_df$distance
# 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]), "%)")


ggplot(pca_rot_stck, aes(x = variable, y = rotation, fill = Sign)) +
    geom_col() + coord_flip() + scale_fill_viridis_d(alpha = 0.7,
    begin = 0.2, end = 0.8) + facet_wrap(~ind_var) + theme_classic()

1 Correlation among metrics

Raw metrics:

Code 
cormat <- cor(degrad_df[, degrad_params], use = "pairwise.complete.obs")

rownames(cormat) <- colnames(cormat) <- degrad_params

cols_corr <- colorRampPalette(c(mako(3, direction = 1, begin = 0.2,
    end = 0.5), "#BEBEBE1A", "white", "#BEBEBE1A", mako(3, direction = 1,
    begin = 0.7, end = 0.9)))(30)

cp <- corrplot.mixed(cormat, tl.cex = 0.7, upper.col = cols_corr,
    lower.col = cols_corr, order = "hclust", lower = "number", upper = "ellipse",
    tl.col = "black")

Code 
# sort parameters as in clusters for cross correlation
# degrad_params <- degrad_params[match(rownames(cp$corr),
# degrad_params)] degrad_df$treatment.replicates <-
# sapply(strsplit(degrad_df$sound.id, '.wav'), '[[', 1)
# degrad_df$frequency <-
# as.numeric(sapply(strsplit(degrad_df$sound.id, '_'), '[[', 1))
# degrad_df$frequency.modulation <-
# sapply(strsplit(degrad_df$sound.id, '_'), '[[', 2)
# degrad_df$amplitude.modulation <-
# sapply(strsplit(degrad_df$sound.id, '_'), '[[', 3)
degrad_df$habitat.structure <- degrad_df$habitat_structure
degrad_df$habitat_structure <- NULL
# degrad_df$duration <- ifelse(grepl('dur:0.1',
# degrad_df$sound.id), 'short', 'long')

2 Data description

  • 20 frequencies
  • 3 locations
  • 11520 test sounds
  • 160 treatment combinations

Sample sizes per location, transect and signal type:

Code 
agg <- aggregate(cbind(sound.id, treatment.replicates) ~ location +
    habitat.structure + distance, degrad_df, function(x) length(unique(x)))

agg$replicates <- agg$sound.id/agg$treatment.replicates

df1 <- knitr::kable(agg, row.names = FALSE, escape = FALSE, format = "html")

kable_styling(df1, bootstrap_options = c("striped", "hover", "condensed",
    "responsive"), full_width = FALSE, font_size = 15)
location habitat.structure distance sound.id treatment.replicates replicates
1 closed 10 480 160 3
2 closed 10 480 160 3
3 closed 10 480 160 3
1 open 10 480 160 3
2 open 10 480 160 3
3 open 10 480 160 3
1 closed 30 480 160 3
2 closed 30 480 160 3
3 closed 30 480 160 3
1 open 30 480 160 3
2 open 30 480 160 3
3 open 30 480 160 3
1 closed 65 480 160 3
2 closed 65 480 160 3
3 closed 65 480 160 3
1 open 65 480 160 3
2 open 65 480 160 3
3 open 65 480 160 3
1 closed 100 480 160 3
2 closed 100 480 160 3
3 closed 100 480 160 3
1 open 100 480 160 3
2 open 100 480 160 3
3 open 100 480 160 3

3 Statistical analysis

Code 
iter <- 20000
chains <- 4
priors <- c(prior(normal(0, 6), class = "b"), prior(cauchy(0, 4),
    class = "sd"))
null_priors <- c(prior(cauchy(0, 4), class = "sd"))

# set frequency to mean-centered and divide by twice the
# standard deviation
degrad_df$raw_frequency <- degrad_df$frequency
degrad_df$frequency <- (degrad_df$frequency - mean(degrad_df$frequency))/(2 *
    sd(degrad_df$frequency))
# set base level for factors
degrad_df$habitat.structure <- factor(degrad_df$habitat.structure,
    levels = c("open", "closed"))
degrad_df$frequency.modulation <- factor(degrad_df$frequency.modulation,
    levels = c("no_fm", "fm"))
degrad_df$amplitude.modulation <- factor(degrad_df$amplitude.modulation,
    levels = c("no_am", "am"))
degrad_df$duration <- factor(degrad_df$duration, levels = c("short",
    "long"))
degrad_df$location <- as.factor(degrad_df$location)
degrad_df$distance_f <- paste0(degrad_df$distance, "m")
degrad_df$distance_f <- factor(degrad_df$distance_f, levels = c("10m",
    "30m", "65m", "100m"), ordered = TRUE)
set.seed(123)

mod_pc1 <- brm(formula = pc1 ~ frequency + frequency.modulation +
    amplitude.modulation + duration + habitat.structure + mo(distance_f) +
    (1 | location) + (1 | treatment.replicates), data = degrad_df,
    prior = priors, iter = iter, chains = chains, cores = chains,
    backend = "cmdstanr", control = list(adapt_delta = 0.99, max_treedepth = 15),
    file = "./data/processed/regression_model_pc1.RDS", file_refit = "always")

mod_pc1 <- add_criterion(mod_pc1, criterion = c("loo"))


null_mod_pc1 <- brm(formula = pc1 ~ 1 + mo(distance_f) + (1 | location) +
    (1 | treatment.replicates), data = degrad_df, prior = null_priors,
    iter = iter, chains = chains, cores = chains, backend = "cmdstanr",
    control = list(adapt_delta = 0.99, max_treedepth = 15), file = "./data/processed/null_regression_model_pc1.RDS",
    file_refit = "always")

null_mod_pc1 <- add_criterion(null_mod_pc1, criterion = c("loo"))

mod_blurratio <- brm(formula = blur.ratio ~ frequency + frequency.modulation +
    amplitude.modulation + duration + habitat.structure + mo(distance_f) +
    (1 | location) + (1 | treatment.replicates), data = degrad_df,
    prior = priors, iter = iter, chains = chains, cores = chains,
    backend = "cmdstanr", control = list(adapt_delta = 0.99, max_treedepth = 15),
    file = "./data/processed/regression_model_blur_ratio.RDS", file_refit = "always")

mod_blurratio <- add_criterion(mod_blurratio, criterion = c("loo"))


null_mod_blurratio <- brm(formula = blur.ratio ~ 1 + mo(distance_f) +
    (1 | location) + (1 | treatment.replicates), data = degrad_df,
    prior = null_priors, iter = iter, chains = chains, cores = chains,
    backend = "cmdstanr", control = list(adapt_delta = 0.99, max_treedepth = 15),
    file = "./data/processed/null_regression_model_blur_ratio.RDS",
    file_refit = "always")

null_mod_blurratio <- add_criterion(null_mod_blurratio, criterion = c("loo"))


mod_ea <- brm(formula = excess.attenuation ~ frequency + frequency.modulation +
    amplitude.modulation + duration + habitat.structure + mo(distance_f) +
    (1 | location) + (1 | treatment.replicates), data = degrad_df,
    prior = priors, iter = iter, chains = chains, backend = "cmdstanr",
    cores = chains, control = list(adapt_delta = 0.99, max_treedepth = 15),
    file = "./data/processed/regression_model_excess_attenuation.RDS",
    file_refit = "always")

mod_ea <- add_criterion(mod_ea, criterion = c("loo"))


null_mod_ea <- brm(formula = excess.attenuation ~ 1 + mo(distance_f) +
    (1 | location) + (1 | treatment.replicates), data = degrad_df,
    prior = null_priors, iter = iter, chains = chains, backend = "cmdstanr",
    cores = chains, control = list(adapt_delta = 0.99, max_treedepth = 15),
    file = "./data/processed/null_regression_model_excess_attenuation.RDS",
    file_refit = "always")

null_mod_ea <- add_criterion(null_mod_ea, criterion = c("loo"))


mod_tsr <- brm(formula = tail.to.signal.ratio ~ frequency + frequency.modulation +
    amplitude.modulation + duration + habitat.structure + mo(distance_f) +
    (1 | location) + (1 | treatment.replicates), data = degrad_df,
    prior = priors, iter = iter, chains = chains, backend = "cmdstanr",
    cores = chains, control = list(adapt_delta = 0.99, max_treedepth = 15),
    file = "./data/processed/regression_model_tail_to_signal_ratio.RDS",
    file_refit = "always")

mod_tsr <- add_criterion(mod_tsr, criterion = c("loo"))

null_mod_tsr <- brm(formula = tail.to.signal.ratio ~ 1 + mo(distance_f) +
    (1 | location) + (1 | treatment.replicates), data = degrad_df,
    prior = null_priors, iter = iter, chains = chains, backend = "cmdstanr",
    cores = chains, control = list(adapt_delta = 0.99, max_treedepth = 15),
    file = "./data/processed/null_regression_model_tail_to_signal_ratio.RDS",
    file_refit = "always")

null_mod_tsr <- add_criterion(null_mod_tsr, criterion = c("loo"))

3.1 PC1 degradation

Model selection:

Code 
effects <- c("habitat structure", "frequency modulation", "amplitude modulation",
    "frequency", "duration")

mod <- readRDS("./data/processed/regression_model_pc1.RDS")

null_mod <- readRDS("./data/processed/null_regression_model_pc1.RDS")

loo_diff <- loo::loo_compare(mod, null_mod)

as.data.frame(loo_diff[, 1:2], row.names = c("model", "null model"))
elpd_diff se_diff
model 0.0000 0.00000
null model -927.2961 38.64656
Code 
extended_summary(mod, n.posterior = 1000, fill = mako(10)[7], trace.palette = my.viridis,
    remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_",
        "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm",
        "durationlong"), gsub.replacement = c("", "habitat structure",
        "amplitude modulation", "frequency modulation", "duration"),
    effects = effects, print.name = FALSE, trace = TRUE, return = FALSE)
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 6) Intercept-student_t(3, -0.5, 2.5) sd-cauchy(0, 4) sigma-student_t(3, 0, 2.5) simo-dirichlet(1) pc1 ~ frequency + frequency.modulation + amplitude.modulation + duration + habitat.structure + mo(distance_f) + (1 | location) + (1 | treatment.replicates) 20000 4 1 10000 45 (0.0011%) 0 3580.76 7958.255 401962879
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
habitat structure 1.105 1.056 1.154 1 50030.556 30071.858
frequency modulation 0.371 0.180 0.563 1.001 4010.818 9243.320
amplitude modulation 0.229 0.040 0.418 1.001 4241.791 8136.630
frequency 1.260 1.071 1.448 1.002 4177.398 8329.452
duration -0.154 -0.343 0.032 1.001 3580.760 7958.255

Posterior predictive checks:

Code 
custom_ppc(fit = mod, group = "habitat.structure")

3.2 Blur ratio

Model selection:

Code 
mod <- readRDS("./data/processed/regression_model_blur_ratio.RDS")

null_mod <- readRDS("./data/processed/null_regression_model_blur_ratio.RDS")

loo_diff <- loo::loo_compare(mod, null_mod)

as.data.frame(loo_diff[, 1:2], row.names = c("model", "null model"))
elpd_diff se_diff
model 0.0000 0.00000
null model -296.3981 23.65781
Code 
extended_summary(fit = mod, n.posterior = 1000, fill = mako(10)[7],
    trace.palette = my.viridis, remove.intercepts = TRUE, highlight = TRUE,
    gsub.pattern = c("b_", "habitat.structureclosed", "amplitude.modulationam",
        "frequency.modulationfm", "durationlong"), gsub.replacement = c("",
        "habitat structure", "amplitude modulation", "frequency modulation",
        "duration"), effects = effects, print.name = FALSE)
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 6) Intercept-student_t(3, 0.1, 2.5) sd-cauchy(0, 4) sigma-student_t(3, 0, 2.5) simo-dirichlet(1) blur.ratio ~ frequency + frequency.modulation + amplitude.modulation + duration + habitat.structure + mo(distance_f) + (1 | location) + (1 | treatment.replicates) 20000 4 1 10000 264 (0.0066%) 0 4774.373 9146.111 909976274
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
habitat structure 0.023 0.021 0.025 1 65083.502 27720.147
frequency modulation 0.055 0.048 0.062 1.001 4867.721 9146.111
amplitude modulation 0.024 0.018 0.031 1.001 4774.373 9811.482
frequency 0.015 0.008 0.021 1 5254.483 9784.451
duration 0.003 -0.004 0.010 1.001 4992.844 9899.127

Posterior predictive checks:

Code 
custom_ppc(fit = mod, group = "habitat.structure")

3.3 Excess attenuation

Model selection:

Code 
mod <- readRDS("./data/processed/regression_model_excess_attenuation.RDS")

null_mod <- readRDS("./data/processed/null_regression_model_excess_attenuation.RDS")

loo_diff <- loo::loo_compare(mod, null_mod)

as.data.frame(loo_diff[, 1:2], row.names = c("model", "null model"))
elpd_diff se_diff
model 0.0000 0.00000
null model -892.8759 41.65656
Code 
extended_summary(fit = mod, n.posterior = 1000, fill = mako(10)[7],
    trace.palette = my.viridis, remove.intercepts = TRUE, highlight = TRUE,
    gsub.pattern = c("b_", "habitat.structureclosed", "amplitude.modulationam",
        "frequency.modulationfm", "durationlong"), gsub.replacement = c("",
        "habitat structure", "amplitude modulation", "frequency modulation",
        "duration"), effects = effects, print.name = FALSE)
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 6) Intercept-student_t(3, 3.8, 8.4) sd-cauchy(0, 4) sigma-student_t(3, 0, 8.4) simo-dirichlet(1) excess.attenuation ~ frequency + frequency.modulation + amplitude.modulation + duration + habitat.structure + mo(distance_f) + (1 | location) + (1 | treatment.replicates) 20000 4 1 10000 19 (0.00048%) 0 5367.054 9998.283 292208958
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
habitat structure 4.150 3.964 4.338 1 60198.441 28202.820
frequency modulation -1.525 -2.177 -0.874 1.001 5367.054 9998.283
amplitude modulation -1.342 -1.981 -0.710 1.002 5512.652 10710.716
frequency 7.381 6.725 8.029 1.001 5581.344 10529.779
duration -1.134 -1.782 -0.487 1 5637.607 10639.855

Posterior predictive checks:

Code 
custom_ppc(fit = mod, group = "habitat.structure")

3.4 Tail-to-signal ratio

Model selection:

Code 
mod <- readRDS("./data/processed/regression_model_tail_to_signal_ratio.RDS")

null_mod <- readRDS("./data/processed/null_regression_model_tail_to_signal_ratio.RDS")

loo_diff <- loo::loo_compare(mod, null_mod)

as.data.frame(loo_diff[, 1:2], row.names = c("model", "null model"))
elpd_diff se_diff
model 0.0000 0.00000
null model -830.3539 38.27065
Code 
extended_summary(fit = mod, n.posterior = 1000, fill = mako(10)[7],
    trace.palette = my.viridis, remove.intercepts = TRUE, highlight = TRUE,
    gsub.pattern = c("b_", "habitat.structureclosed", "amplitude.modulationam",
        "frequency.modulationfm", "durationlong"), gsub.replacement = c("",
        "habitat structure", "amplitude modulation", "frequency modulation",
        "duration"), effects = effects, print.name = FALSE)
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 6) Intercept-student_t(3, -13.3, 11.8) sd-cauchy(0, 4) sigma-student_t(3, 0, 11.8) simo-dirichlet(1) tail.to.signal.ratio ~ frequency + frequency.modulation + amplitude.modulation + duration + habitat.structure + mo(distance_f) + (1 | location) + (1 | treatment.replicates) 20000 4 1 10000 8 (2e-04%) 0 3181.828 7022.082 1460541246
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
habitat structure 5.026 4.792 5.260 1 53141.037 29937.008
frequency modulation 1.110 0.193 2.027 1.001 3181.828 7022.082
amplitude modulation 0.748 -0.151 1.650 1.001 3923.471 8501.926
frequency 4.274 3.380 5.161 1.001 3746.194 8619.106
duration -0.719 -1.608 0.170 1.001 3798.823 7796.815

Code 
# mod <- readRDS('./data/processed/regression_model_int.RDS')

Posterior predictive checks:

Code 
custom_ppc(fit = mod, group = "habitat.structure")

3.5 Combined results

Code 
pc1 <- extended_summary(read.file = "./data/processed/regression_model_pc1.RDS",
    n.posterior = 1000, fill = mako(10)[7], trace.palette = my.viridis,
    remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_",
        "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm",
        "durationlong"), gsub.replacement = c("", "habitat structure",
        "amplitude modulation", "frequency modulation", "duration"),
    effects = effects, print.name = FALSE, trace = FALSE, return = TRUE)

# gg_pc1 <- pc1$plot + labs(x = 'Degradation (PC1)')

br <- extended_summary(read.file = "./data/processed/regression_model_blur_ratio.RDS",
    n.posterior = 1000, fill = mako(10)[7], trace.palette = my.viridis,
    remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_",
        "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm",
        "durationlong"), gsub.replacement = c("", "habitat structure",
        "amplitude modulation", "frequency modulation", "duration"),
    effects = effects, print.name = FALSE, trace = FALSE, return = TRUE)

# gg_br <- br$plot + labs(x = 'Blur ratio', y = '') +
# theme(axis.text.y = element_blank())

ea <- extended_summary(read.file = "./data/processed/regression_model_excess_attenuation.RDS",
    n.posterior = 1000, fill = mako(10)[7], trace.palette = my.viridis,
    remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_",
        "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm",
        "durationlong"), gsub.replacement = c("", "habitat structure",
        "amplitude modulation", "frequency modulation", "duration"),
    effects = effects, print.name = FALSE, trace = FALSE, return = TRUE)

# gg_ea <- ea$plot + labs(x = 'Excess attenuation', y = '') +
# theme(axis.text.y = element_blank())

tsr <- extended_summary(read.file = "./data/processed/regression_model_tail_to_signal_ratio.RDS",
    n.posterior = 1000, fill = mako(10)[7], trace.palette = my.viridis,
    remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_",
        "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm",
        "durationlong"), gsub.replacement = c("", "habitat structure",
        "amplitude modulation", "frequency modulation", "duration"),
    effects = effects, print.name = FALSE, trace = FALSE, return = TRUE)

# gg_tsr <- tsr$plot + labs(x = 'Tail-to-signal ratio', y = '')
# + theme(axis.text.y = element_blank()) plot_grid(gg_pc1,
# gg_br, gg_ea, gg_tsr, nrow = 1, rel_widths = c(6, 2, 2, 2))
Code 
estimates <- data.frame(pc1 = pc1$coef_table$Estimate, br = br$coef_table$Estimate,
    ea = ea$coef_table$Estimate, tsr = tsr$coef_table$Estimate)

names(estimates) <- c("Degradation (PC1)", "Blur ratio", "Excess attenuation",
    "Tail-to-signal ratio")

# make estimates relative to maximum estimate in data
rel_estimates <- data.frame(lapply(estimates, function(x) x/max(abs(x)) *
    0.8))

names(rel_estimates) <- c("Degradation (PC1)", "Blur ratio", "Excess attenuation",
    "Tail-to-signal ratio")

# corrplot(as.matrix(estimates), method = 'ellipse', )

signif <- data.frame(pc1 = pc1$coef_table$`l-95% CI` * pc1$coef_table$`u-95% CI` >
    0, br = br$coef_table$`l-95% CI` * br$coef_table$`u-95% CI` >
    0, ea = ea$coef_table$`l-95% CI` * ea$coef_table$`u-95% CI` >
    0, tsr = tsr$coef_table$`l-95% CI` * tsr$coef_table$`u-95% CI` >
    0)

names(signif) <- c("Degradation (PC1)", "Blur ratio", "Excess attenuation",
    "Tail-to-signal ratio")

estimates <- as.data.frame(cbind(rownames(pc1$coef_table), stack(estimates)[,
    1], stack(rel_estimates)[, 1], stack(signif)))

names(estimates) <- c("predictor", "est", "relavite.est", "sig", "response")
estimates$signif <- ifelse(estimates$sig, "p < 0.05", "N.S.")


# estimates$relavite.est <- ifelse(estimates$sig,
# estimates$relavite.est, 0)

estimates$response <- factor(estimates$response, levels = rev(c("Degradation (PC1)",
    "Blur ratio", "Excess attenuation", "Tail-to-signal ratio")))

firstup <- function(x) {
    substr(x, 1, 1) <- toupper(substr(x, 1, 1))
    x
}

estimates$predictor <- firstup(estimates$predictor)

estimates$predictor <- factor(estimates$predictor, levels = c("Habitat structure",
    "Frequency", "Duration", "Amplitude modulation", "Frequency modulation"))

gg_tiled <- ggplot(estimates, aes(x = predictor, y = response, fill = relavite.est)) +
    geom_tile() + coord_equal() + scale_fill_gradient2(low = mako(10)[3],
    high = mako(10)[7], guide = "none") + geom_text(aes(label = round(est,
    3), color = signif), size = 5) + scale_color_manual(values = c("gray",
    "black"), guide = "none") + labs(x = "", y = "", color = "Significance") +
    theme_classic(base_size = 15) + theme(axis.text.x = element_text(color = "black",
    angle = 30, vjust = 1, hjust = 1), axis.text.y = element_text(color = "black"))

gg_tiled

Code 
gg_tiled2 <- ggplot(estimates, aes(x = predictor, y = response, fill = relavite.est)) +
    geom_tile() + coord_equal() + scale_fill_gradient2(low = mako(10)[3],
    high = mako(10)[7], guide = "none") + geom_text(aes(label = round(est,
    3), color = signif), size = 7) + scale_color_manual(values = c("gray",
    "black"), guide = "none") + labs(x = "", y = "", color = "Significance") +
    theme_classic(base_size = 22) + theme(axis.text.x = element_text(color = "black",
    angle = 30, vjust = 1, hjust = 1), axis.text.y = element_text(color = "black"))


# # save for presentation ggsave(plot = gg_tiled2, filename =
# './output/presentation/images/pvalues_grid.jpeg', dpi = 300)

# save for paper ggsave(plot = gg_tiled2, filename =
# './output/pvalues_grid.jpeg', dpi = 300,width = 12, height =
# 8)

3.6 Combined model metadata

Code 
check_rds_fits(path = "./data/processed/", html = TRUE)
Some RDS files cannot be read (check '.Options$brmsish$not_read_rds')
priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 Intercept-student_t(3, 0.1, 2.5) sd-cauchy(0, 4) sigma-student_t(3, 0, 2.5) simo-dirichlet(1) blur.ratio ~ 1 + mo(distance_f) + (1 | location) + (1 | treatment.replicates) 20000 4 1 10000 267 (0.006675%) 0 3080.848 5752.072 1243849502
2 Intercept-student_t(3, 3.8, 8.4) sd-cauchy(0, 4) sigma-student_t(3, 0, 8.4) simo-dirichlet(1) excess.attenuation ~ 1 + mo(distance_f) + (1 | location) + (1 | treatment.replicates) 20000 4 1 10000 13 (0.000325%) 0 2091.222 4600.028 821488295
3 Intercept-student_t(3, -0.5, 2.5) sd-cauchy(0, 4) sigma-student_t(3, 0, 2.5) simo-dirichlet(1) pc1 ~ 1 + mo(distance_f) + (1 | location) + (1 | treatment.replicates) 20000 4 1 10000 51 (0.001275%) 0 3715.405 7221.713 401962879
4 Intercept-student_t(3, -13.3, 11.8) sd-cauchy(0, 4) sigma-student_t(3, 0, 11.8) simo-dirichlet(1) tail.to.signal.ratio ~ 1 + mo(distance_f) + (1 | location) + (1 | treatment.replicates) 20000 4 1 10000 21 (0.000525%) 0 5143.574 9772.178 776267316
5 b-normal(0, 6) Intercept-student_t(3, 0.1, 2.5) sd-cauchy(0, 4) sigma-student_t(3, 0, 2.5) simo-dirichlet(1) blur.ratio ~ frequency + frequency.modulation + amplitude.modulation + duration + habitat.structure + mo(distance_f) + (1 | location) + (1 | treatment.replicates) 20000 4 1 10000 264 (0.0066%) 0 4774.373 NA 909976274
6 b-normal(0, 6) Intercept-student_t(3, 0.1, 2.5) sd-cauchy(0, 4) sigma-student_t(3, 0, 2.5) simo-dirichlet(1) blur.ratio ~ frequency + frequency.modulation + amplitude.modulation + duration + habitat.structure + mo(distance_f) + (1 | location) + (1 | treatment.replicates) 20000 4 1 10000 264 (0.0066%) 0 4774.373 NA 909976274
7 b-normal(0, 6) Intercept-student_t(3, 3.8, 8.4) sd-cauchy(0, 4) sigma-student_t(3, 0, 8.4) simo-dirichlet(1) excess.attenuation ~ frequency + frequency.modulation + amplitude.modulation + duration + habitat.structure + mo(distance_f) + (1 | location) + (1 | treatment.replicates) 20000 4 1 10000 19 (0.000475%) 0 5367.054 9998.283 292208958
8 b-normal(0, 6) Intercept-student_t(3, -0.5, 2.5) sd-cauchy(0, 4) sigma-student_t(3, 0, 2.5) simo-dirichlet(1) pc1 ~ frequency + frequency.modulation + amplitude.modulation + duration + habitat.structure + mo(distance_f) + (1 | location) + (1 | treatment.replicates) 20000 4 1 10000 45 (0.001125%) 0 3580.760 7958.255 401962879
9 b-normal(0, 6) Intercept-student_t(3, -13.3, 11.8) sd-cauchy(0, 4) sigma-student_t(3, 0, 11.8) simo-dirichlet(1) tail.to.signal.ratio ~ frequency + frequency.modulation + amplitude.modulation + duration + habitat.structure + mo(distance_f) + (1 | location) + (1 | treatment.replicates) 20000 4 1 10000 8 (2e-04%) 0 3181.828 7022.082 1460541246

 

4 Takeaways

  • Habitat structure seems to drive most of the observed degradation
  • Frequency modulation and amplitude modulation were the signal structural features that more strongly affected transmission
  • Signal features interact with habitat structure in diverse ways, sometimes increasing and sometimes decreasing degradation
  • Most degradation metrics are robust to variation in signal duration

 

Session information

R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.4 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       

time zone: America/Costa_Rica
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] loo_2.8.0           rstan_2.32.6        StanHeaders_2.32.10
 [4] ggsignif_0.6.4      emmeans_1.10.3      brmsish_1.0.0      
 [7] cmdstanr_0.8.1.9000 cowplot_1.1.3       ggdist_3.3.2       
[10] brms_2.22.0         Rcpp_1.0.13         corrplot_0.92      
[13] viridis_0.6.5       viridisLite_0.4.2   tidyr_1.3.1        
[16] ggplot2_3.5.1       warbleR_1.1.32      NatureSounds_1.0.4 
[19] seewave_2.2.3       tuneR_1.4.7         rprojroot_2.0.4    
[22] formatR_1.14        knitr_1.48          kableExtra_1.4.0   
[25] remotes_2.5.0      

loaded via a namespace (and not attached):
 [1] bitops_1.0-9         pbapply_1.7-2        gridExtra_2.3       
 [4] inline_0.3.19        testthat_3.2.1.1     sandwich_3.1-0      
 [7] rlang_1.1.4          magrittr_2.0.3       multcomp_1.4-25     
[10] matrixStats_1.4.1    compiler_4.4.1       reshape2_1.4.4      
[13] systemfonts_1.1.0    vctrs_0.6.5          stringr_1.5.1       
[16] arrayhelpers_1.1-0   pkgconfig_2.0.3      crayon_1.5.3        
[19] fastmap_1.2.0        backports_1.5.0      labeling_0.4.3      
[22] utf8_1.2.4           rmarkdown_2.28       ps_1.8.0            
[25] purrr_1.0.2          xfun_0.48            jsonlite_1.8.9      
[28] tidybayes_3.0.7      parallel_4.4.1       R6_2.5.1            
[31] stringi_1.8.4        brio_1.1.5           estimability_1.5.1  
[34] zoo_1.8-12           bayesplot_1.11.1     Matrix_1.7-0        
[37] splines_4.4.1        tidyselect_1.2.1     rstudioapi_0.16.0   
[40] abind_1.4-8          yaml_2.3.10          dtw_1.23-1          
[43] codetools_0.2-20     curl_5.2.2           processx_3.8.4      
[46] pkgbuild_1.4.4       plyr_1.8.9           lattice_0.22-6      
[49] tibble_3.2.1         withr_3.0.1          bridgesampling_1.1-2
[52] posterior_1.6.0      coda_0.19-4.1        evaluate_1.0.0      
[55] signal_1.8-1         survival_3.7-0       sketchy_1.0.3       
[58] proxy_0.4-27         RcppParallel_5.1.9   xml2_1.3.6          
[61] pillar_1.9.0         tensorA_0.36.2.1     packrat_0.9.2       
[64] stats4_4.4.1         checkmate_2.3.2      distributional_0.5.0
[67] generics_0.1.3       RCurl_1.98-1.16      rstantools_2.4.0    
[70] munsell_0.5.1        scales_1.3.0         xtable_1.8-4        
[73] glue_1.8.0           tools_4.4.1          xaringanExtra_0.8.0 
[76] mvtnorm_1.3-1        grid_4.4.1           ape_5.8             
[79] QuickJSR_1.4.0       colorspace_2.1-1     nlme_3.1-165        
[82] cli_3.6.3            svUnit_1.0.6         fansi_1.0.6         
[85] svglite_2.1.3        Brobdingnag_1.2-9    dplyr_1.1.4         
[88] V8_4.4.2             gtable_0.3.5         fftw_1.0-9          
[91] digest_0.6.37        TH.data_1.1-2        farver_2.1.2        
[94] rjson_0.2.23         htmlwidgets_1.6.4    htmltools_0.5.8.1   
[97] lifecycle_1.0.4      MASS_7.3-61