Stats

Vocal activity and noise in birds of Santa Rosa NP

Author

Published

September 23, 2025

Code

# print link to github repo if any
if (file.exists("./.git/config")){
  config <- readLines("./.git/config")
  url <- grep("url",  config, value = TRUE)
  url <- gsub("\\turl = |.git$", "", url)
  cat("\nSource code and data found at [", url, "](", url, ")", sep = "")
  }

Source code and data found at

Code

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

1 Custom functions

Code

brms_summary <- function(...) {
    extended_summary(..., highlight = TRUE, remove.intercepts = TRUE,
        trace.palette = viridis::mako, fill = "#54C9ADFF")
}

Purpose

Fit regression models

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("knitr", "formatR", "brms", "ggplot2",
    "readxl", github = "maRce10/brmsish", "viridis"))

Code

dat <- read_excel("./data/raw/MS_Chaves-Acuñaetal_AvesStaRosa_julio2025_revWCA.xlsx",
    sheet = "datos_largos", na = "N/A")


dat <- dat[, 1:which(names(dat) == "numero_cantos")]


n_sp <- table(dat$`Especie Ave`)


n_sp <- sort(n_sp, decreasing = T)

# select those spp with more than 100 records
dat <- dat[dat$`Especie Ave` %in% names(n_sp)[n_sp > 100], ]


dat$traslape_frecuencia <- ifelse(dat$traslape_frecuencia == "Antrofonía",
    "antrofonia", "biofonia")

dat$traslape_frecuencia_bn <- ifelse(dat$traslape_frecuencia == dat$tipo_ruido,
    "si", "no")

dat$especie <- dat$`Especie Ave`

2 Regression models

2.1 Overlap

Code

# make dat == 0 a bit above 0
dat$prop_traslape[dat$prop_traslape == 0] <- 1e-04

# combine duration and frequency overlap
dat$duracion_traslape <- ifelse(dat$traslape_frecuencia_bn == "si",
    dat$duracion_tratamiento, 0)
dat$duracion_traslape <- scale(dat$duracion_traslape)

# scale duration of treatment and song duration
dat$duracion_tratamiento <- scale(dat$duracion_tratamiento)
# dat$duracion_canto <- scale(dat$duracion_canto)


cores <- chains <- 4
iter <- 2000


# proporcion de traslape
proporcion_simple <- brm(prop_traslape ~ traslape_frecuencia_bn +
    tipo_ruido + (1 | Sitio/Punto) + (1 | especie), data = dat, iter = iter,
    cores = cores, chains = chains, family = Beta(link = "logit",
        link_phi = "log"), backend = "cmdstanr", prior = c(prior(normal(0,
        5), "b"), prior(normal(0, 24), "Intercept")), control = list(adapt_delta = 0.99,
        max_treedepth = 15), file = "./data/processed/proporcion_simple.rds",
    file_refit = "on_change")


proporcion_interaccion_duracion_tratamiento <- brm(prop_traslape ~
    traslape_frecuencia_bn + tipo_ruido * duracion_tratamiento + (1 |
        Sitio/Punto) + (1 | especie), data = dat, iter = iter, cores = cores,
    chains = chains, family = Beta(link = "logit", link_phi = "log"),
    backend = "cmdstanr", prior = c(prior(normal(0, 5), "b"), prior(normal(0,
        24), "Intercept")), control = list(adapt_delta = 0.99, max_treedepth = 15),
    file = "./data/processed/proporcion_interaccion_duracion_tratamiento.rds",
    file_refit = "on_change")

proporcion_interaccion_traslape <- brm(prop_traslape ~ traslape_frecuencia_bn *
    tipo_ruido + duracion_tratamiento + (1 | Sitio/Punto) + (1 | especie),
    data = dat, iter = iter, cores = cores, chains = chains, family = Beta(link = "logit",
        link_phi = "log"), backend = "cmdstanr", prior = c(prior(normal(0,
        5), "b"), prior(normal(0, 24), "Intercept")), control = list(adapt_delta = 0.99,
        max_treedepth = 15), file = "./data/processed/proporcion_interaccion_traslape.rds",
    file_refit = "on_change")

proporcion_interaccion_duracion_traslape <- brm(prop_traslape ~ traslape_frecuencia_bn *
    tipo_ruido * duracion_tratamiento + (1 | Sitio/Punto) + (1 | especie),
    data = dat, iter = iter, cores = cores, chains = chains, family = Beta(link = "logit",
        link_phi = "log"), backend = "cmdstanr", prior = c(prior(normal(0,
        5), "b"), prior(normal(0, 24), "Intercept")), control = list(adapt_delta = 0.99,
        max_treedepth = 15), file = "./data/processed/proporcion_interaccion_duracion_traslape.rds",
    file_refit = "on_change")

proporcion_duracion_traslape <- brm(prop_traslape ~ tipo_ruido + duracion_traslape +
    (1 | Sitio/Punto) + (1 | especie), data = dat, iter = iter, cores = cores,
    chains = chains, family = Beta(link = "logit", link_phi = "log"),
    backend = "cmdstanr", prior = c(prior(normal(0, 5), "b"), prior(normal(0,
        24), "Intercept")), control = list(adapt_delta = 0.99, max_treedepth = 15),
    file = "./data/processed/proporcion_duracion_traslape.rds", file_refit = "on_change")

proporcion_duracion_traslape_interaccion <- brm(prop_traslape ~ tipo_ruido *
    duracion_traslape + (1 | Sitio/Punto) + (1 | especie), data = dat,
    iter = iter, cores = cores, chains = chains, family = Beta(link = "logit",
        link_phi = "log"), backend = "cmdstanr", prior = c(prior(normal(0,
        5), "b"), prior(normal(0, 24), "Intercept")), control = list(adapt_delta = 0.99,
        max_treedepth = 15), file = "./data/processed/proporcion_duracion_traslape_interaccion.rds",
    file_refit = "on_change")


proporcion_interaccion_nulo <- brm(prop_traslape ~ 1 + (1 | Sitio/Punto) +
    (1 | especie), data = dat, iter = iter, cores = cores, chains = chains,
    family = Beta(link = "logit", link_phi = "log"), backend = "cmdstanr",
    prior = c(prior(normal(0, 24), "Intercept")), control = list(adapt_delta = 0.99,
        max_treedepth = 15), file = "./data/processed/proporcion_interaccion_nulo.rds",
    file_refit = "on_change")

2.1.1 Compare models

Code

mod_list <- list(proporcion_simple = proporcion_simple, proporcion_interaccion_duracion_tratamiento = proporcion_interaccion_duracion_tratamiento,
    proporcion_interaccion_traslape = proporcion_interaccion_traslape,
    proporcion_interaccion_duracion_traslape = proporcion_interaccion_duracion_traslape,
    proporcion_interaccion_nulo = proporcion_interaccion_nulo, proporcion_duracion_traslape = proporcion_duracion_traslape,
    proporcion_duracion_traslape_interaccion = proporcion_duracion_traslape_interaccion)

mod_loo <- lapply(mod_list, loo)

loo_comparison <- loo::loo_compare(mod_loo)

loo_comparison

                                            elpd_diff se_diff
proporcion_interaccion_duracion_traslape      0.0       0.0  
proporcion_interaccion_duracion_tratamiento  -8.2       3.8  
proporcion_interaccion_traslape              -9.4       3.8  
proporcion_duracion_traslape_interaccion    -47.3      10.5  
proporcion_duracion_traslape                -48.5       9.9  
proporcion_simple                           -59.2      11.5  
proporcion_interaccion_nulo                 -78.1      12.7

2.1.1.1 Best model

Code

best_mod <- mod_list[[which(names(mod_list) == rownames(loo_comparison)[1])]]

brms_summary(best_mod, gsub.pattern = c("b_treatment|b_", "tipo_ruido"),
    gsub.replacement = c("", "antropof_"))

2.2 best_mod

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 5) Intercept-normal(0, 24) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 2.5)	prop_traslape ~ traslape_frecuencia_bn * tipo_ruido * duracion_tratamiento + (1 \| Sitio/Punto) + (1 \| especie)	2000	4	1	1000	32 (0.008%)	0	1321.085	1939.343	1592323471

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
traslape_frecuencia_bnsi	-0.045	-0.265	0.183	1.002	1440.718	2261.919
antropof_biofonia	-0.045	-0.262	0.169	1.002	1431.306	1939.343
antropof_ventana	-0.047	-0.175	0.086	1.001	1544.025	2529.985
duracion_tratamiento	0.117	0.025	0.213	1.003	1321.085	2278.416
traslape_frecuencia_bnsi:antropof_biofonia	-0.036	-0.409	0.338	1.001	1350.923	1949.386
traslape_frecuencia_bnsi:antropof_ventana	0.119	-9.620	9.805	1	4781.570	2823.767
traslape_frecuencia_bnsi:duracion_tratamiento	0.210	0.057	0.356	1.002	1681.703	2662.978
antropof_biofonia:duracion_tratamiento	0.342	0.187	0.487	1.003	1711.509	2374.647
antropof_ventana:duracion_tratamiento	0.128	0.002	0.249	1.003	1565.907	2538.278
traslape_frecuencia_bnsi:antropof_biofonia:duracion_tratamiento	-0.519	-0.738	-0.291	1.002	1946.716	2342.399
traslape_frecuencia_bnsi:antropof_ventana:duracion_tratamiento	-0.166	-9.815	9.256	1.001	4474.550	2990.662

2.2.0.0.1 Marginal effect plots

Code

conditions <- make_conditions(best_mod, "tipo_ruido")

conditional_effects(best_mod, "duracion_tratamiento", conditions = conditions)

Code

conditional_effects(best_mod, "traslape_frecuencia_bn", conditions = conditions)

Code

conditional_effects(best_mod, "traslape_frecuencia_bn:tipo_ruido",
    conditions = conditions)

Code

conditions <- make_conditions(best_mod, "duracion_tratamiento")

conditional_effects(best_mod, "traslape_frecuencia_bn", conditions = conditions)

Code

conditional_effects(best_mod, "tipo_ruido", conditions = conditions)

Code

conditional_effects(best_mod, "traslape_frecuencia_bn:tipo_ruido",
    conditions = conditions)

Code

conditions <- make_conditions(best_mod, "traslape_frecuencia_bn")

conditional_effects(best_mod, "duracion_tratamiento", conditions = conditions)

Code

conditional_effects(best_mod, "tipo_ruido", conditions = conditions)

Code

conditional_effects(best_mod, "duracion_tratamiento:tipo_ruido", conditions = conditions)

2.3 Duration

Code

# make duracion_canto == 0 a bit above 0
dat$duracion_canto[dat$duracion_canto == 0] <- 1e-04

# duracion del canto
duracion_canto_simple <- brm(duracion_canto ~ traslape_frecuencia_bn +
    tipo_ruido + (1 | Sitio/Punto) + (1 | especie), data = dat, iter = iter,
    cores = cores, chains = chains, family = Gamma(link = "log"),
    backend = "cmdstanr", prior = c(prior(normal(0, 5), "b"), prior(normal(0,
        24), "Intercept")), control = list(adapt_delta = 0.99, max_treedepth = 15),
    file = "./data/processed/duracion_canto_simple.rds", file_refit = "on_change")

duracion_canto_interaccion_duracion_tratamiento <- brm(duracion_canto ~
    traslape_frecuencia_bn + tipo_ruido * duracion_tratamiento + (1 |
        Sitio/Punto) + (1 | especie), data = dat, iter = iter, cores = cores,
    chains = chains, family = Gamma(link = "log"), backend = "cmdstanr",
    prior = c(prior(normal(0, 5), "b"), prior(normal(0, 24), "Intercept")),
    control = list(adapt_delta = 0.99, max_treedepth = 15), file = "./data/processed/duracion_canto_interaccion_duracion_tratamiento.rds",
    file_refit = "on_change")


duracion_canto_interaccion_traslape <- brm(duracion_canto ~ traslape_frecuencia_bn *
    tipo_ruido + duracion_tratamiento + (1 | Sitio/Punto) + (1 | especie),
    data = dat, iter = iter, cores = cores, chains = chains, family = Gamma(link = "log"),
    backend = "cmdstanr", prior = c(prior(normal(0, 5), "b"), prior(normal(0,
        24), "Intercept")), control = list(adapt_delta = 0.99, max_treedepth = 15),
    file = "./data/processed/duracion_canto_interaccion_traslape.rds",
    file_refit = "on_change")

duracion_canto_interaccion_duracion_traslape <- brm(duracion_canto ~
    traslape_frecuencia_bn * tipo_ruido * duracion_tratamiento + (1 |
        Sitio/Punto) + (1 | especie), data = dat, iter = iter, cores = cores,
    chains = chains, family = Gamma(link = "log"), backend = "cmdstanr",
    prior = c(prior(normal(0, 5), "b"), prior(normal(0, 24), "Intercept")),
    control = list(adapt_delta = 0.99, max_treedepth = 15), file = "./data/processed/duracion_canto_interaccion_duracion_traslape.rds",
    file_refit = "on_change")

duracion_canto_nulo <- brm(duracion_canto ~ 1 + (1 | Sitio/Punto) +
    (1 | especie), data = dat, iter = iter, cores = cores, chains = chains,
    family = Gamma(link = "log"), backend = "cmdstanr", prior = c(prior(normal(0,
        24), "Intercept")), control = list(adapt_delta = 0.99, max_treedepth = 15),
    file = "./data/processed/duracion_canto_nulo.rds", file_refit = "on_change")

duracion_canto_proporcion_duracion_traslape <- brm(duracion_canto ~
    tipo_ruido + duracion_traslape + (1 | Sitio/Punto) + (1 | especie),
    data = dat, iter = iter, cores = cores, chains = chains, family = Gamma(link = "log"),
    backend = "cmdstanr", prior = c(prior(normal(0, 5), "b"), prior(normal(0,
        24), "Intercept")), control = list(adapt_delta = 0.99, max_treedepth = 15),
    file = "./data/processed/duracion_canto_proporcion_duracion_traslape.rds",
    file_refit = "on_change")

duracion_canto_proporcion_duracion_traslape_interaccion <- brm(duracion_canto ~
    tipo_ruido * duracion_traslape + (1 | Sitio/Punto) + (1 | especie),
    data = dat, iter = iter, cores = cores, chains = chains, family = Gamma(link = "log"),
    backend = "cmdstanr", prior = c(prior(normal(0, 5), "b"), prior(normal(0,
        24), "Intercept")), control = list(adapt_delta = 0.99, max_treedepth = 15),
    file = "./data/processed/duracion_canto_proporcion_duracion_traslape_interaccion.rds",
    file_refit = "on_change")

2.3.1 Compare models

Code

mod_list <- list(duracion_canto_simple = duracion_canto_simple, duracion_canto_interaccion_duracion_tratamiento = duracion_canto_interaccion_duracion_tratamiento,
    duracion_canto_interaccion_traslape = duracion_canto_interaccion_traslape,
    duracion_canto_interaccion_duracion_traslape = duracion_canto_interaccion_duracion_traslape,
    duracion_canto_nulo = duracion_canto_nulo, duracion_canto_proporcion_duracion_traslape_interaccion = duracion_canto_proporcion_duracion_traslape_interaccion,
    duracion_canto_proporcion_duracion_traslape = duracion_canto_proporcion_duracion_traslape)

mod_loo <- lapply(mod_list, loo)

loo_comparison <- loo::loo_compare(mod_loo)

loo_comparison

                                                        elpd_diff se_diff
duracion_canto_interaccion_duracion_traslape               0.0       0.0 
duracion_canto_interaccion_duracion_tratamiento           -6.5       8.2 
duracion_canto_interaccion_traslape                     -103.4      33.1 
duracion_canto_proporcion_duracion_traslape_interaccion -664.5      55.5 
duracion_canto_simple                                   -721.6      59.8 
duracion_canto_proporcion_duracion_traslape             -736.3      59.7 
duracion_canto_nulo                                     -830.7      60.7

2.3.1.1 Best model

Code

best_mod <- mod_list[[which(names(mod_list) == rownames(loo_comparison)[1])]]

brms_summary(best_mod, gsub.pattern = c("b_treatment|b_", "tipo_ruido"),
    gsub.replacement = c("", "antropof_"))

2.4 best_mod

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 5) Intercept-normal(0, 24) sd-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01)	duracion_canto ~ traslape_frecuencia_bn * tipo_ruido * duracion_tratamiento + (1 \| Sitio/Punto) + (1 \| especie)	2000	4	1	1000	7 (0.0018%)	0	2411.682	2393.778	1044184679

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
traslape_frecuencia_bnsi	-0.919	-1.432	-0.408	1	2438.176	2597.105
antropof_biofonia	-1.210	-1.774	-0.610	1.001	2488.353	2393.778
antropof_ventana	-1.197	-1.506	-0.900	1	2989.377	3119.870
duracion_tratamiento	1.366	1.123	1.604	1.002	2837.441	3001.552
traslape_frecuencia_bnsi:antropof_biofonia	-0.623	-1.539	0.253	1.001	2411.682	2453.504
traslape_frecuencia_bnsi:antropof_ventana	0.032	-9.300	9.658	1.002	4996.488	2712.316
traslape_frecuencia_bnsi:duracion_tratamiento	0.564	0.170	0.953	1.001	3421.356	3158.577
antropof_biofonia:duracion_tratamiento	3.690	2.895	4.501	1.001	3165.664	2881.766
antropof_ventana:duracion_tratamiento	1.563	1.204	1.927	1.001	3319.948	2888.373
traslape_frecuencia_bnsi:antropof_biofonia:duracion_tratamiento	-2.286	-3.178	-1.394	1.001	2958.493	2882.381
traslape_frecuencia_bnsi:antropof_ventana:duracion_tratamiento	-0.090	-9.803	9.478	1	5172.947	2832.424

2.4.0.0.1 Marginal effect plots

Code

conditions <- make_conditions(best_mod, "tipo_ruido")

conditional_effects(best_mod, "duracion_tratamiento", conditions = conditions)

Code

conditional_effects(best_mod, "traslape_frecuencia_bn", conditions = conditions)

Code

conditional_effects(best_mod, "traslape_frecuencia_bn:tipo_ruido",
    conditions = conditions)

Code

conditions <- make_conditions(best_mod, "duracion_tratamiento")

conditional_effects(best_mod, "traslape_frecuencia_bn", conditions = conditions)

Code

conditional_effects(best_mod, "tipo_ruido", conditions = conditions)

Code

conditional_effects(best_mod, "traslape_frecuencia_bn:tipo_ruido",
    conditions = conditions)

Code

conditions <- make_conditions(best_mod, "traslape_frecuencia_bn")

conditional_effects(best_mod, "duracion_tratamiento", conditions = conditions)

Code

conditional_effects(best_mod, "tipo_ruido", conditions = conditions)

Code

conditional_effects(best_mod, "duracion_tratamiento:tipo_ruido", conditions = conditions)

3 Questions

Que quiere decir “ninguna” en la columna Especie Ave?

Takeaways

Session information

R version 4.5.0 (2025-04-11)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.5 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  LAPACK version 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] viridis_0.6.5     viridisLite_0.4.2 brmsish_1.0.0     readxl_1.4.5     
[5] ggplot2_3.5.2     brms_2.22.0       Rcpp_1.1.0        formatR_1.14     
[9] knitr_1.50       

loaded via a namespace (and not attached):
 [1] tidyselect_1.2.1     svUnit_1.0.6         dplyr_1.1.4         
 [4] farver_2.1.2         loo_2.8.0            tidybayes_3.0.7     
 [7] fastmap_1.2.0        TH.data_1.1-3        tensorA_0.36.2.1    
[10] digest_0.6.37        estimability_1.5.1   lifecycle_1.0.4     
[13] StanHeaders_2.32.10  processx_3.8.6       survival_3.8-3      
[16] magrittr_2.0.3       posterior_1.6.1      compiler_4.5.0      
[19] rlang_1.1.6          tools_4.5.0          yaml_2.3.10         
[22] labeling_0.4.3       bridgesampling_1.1-2 htmlwidgets_1.6.4   
[25] curl_7.0.0           pkgbuild_1.4.8       plyr_1.8.9          
[28] cmdstanr_0.9.0       xml2_1.3.8           RColorBrewer_1.1-3  
[31] multcomp_1.4-28      abind_1.4-8          withr_3.0.2         
[34] purrr_1.1.0          stats4_4.5.0         grid_4.5.0          
[37] inline_0.3.21        xtable_1.8-4         emmeans_1.11.1      
[40] scales_1.4.0         MASS_7.3-65          cli_3.6.5           
[43] mvtnorm_1.3-3        rmarkdown_2.29       crayon_1.5.3        
[46] generics_0.1.4       remotes_2.5.0        RcppParallel_5.1.10 
[49] rstudioapi_0.17.1    reshape2_1.4.4       pbapply_1.7-4       
[52] ape_5.8-1            rstan_2.32.7         stringr_1.5.1       
[55] splines_4.5.0        bayesplot_1.12.0     parallel_4.5.0      
[58] cellranger_1.1.0     matrixStats_1.5.0    vctrs_0.6.5         
[61] V8_6.0.4             Matrix_1.7-3         sandwich_3.1-1      
[64] jsonlite_2.0.0       arrayhelpers_1.1-0   systemfonts_1.2.3   
[67] ggdist_3.3.3         packrat_0.9.3        tidyr_1.3.1         
[70] glue_1.8.0           ps_1.9.1             codetools_0.2-20    
[73] cowplot_1.1.3        distributional_0.5.0 xaringanExtra_0.8.0 
[76] stringi_1.8.7        gtable_0.3.6         QuickJSR_1.7.0      
[79] tibble_3.3.0         pillar_1.11.0        htmltools_0.5.8.1   
[82] Brobdingnag_1.2-9    R6_2.6.1             evaluate_1.0.4      
[85] kableExtra_1.4.0     lattice_0.22-7       backports_1.5.0     
[88] sketchy_1.0.5        rstantools_2.4.0     gridExtra_2.3       
[91] svglite_2.1.3        coda_0.19-4.1        nlme_3.1-168        
[94] checkmate_2.3.3      xfun_0.53            zoo_1.8-14          
[97] pkgconfig_2.0.3