Statistical analysis

Bellbird mimicry

# options to customize chunk outputs
knitr::opts_chunk$set(
  message = FALSE
 )

Purpose

• Evaluate similarity of vocalizations of a mimic bellbird and its model species using birdNET embeddings
  • Are focal belbird vocalizations more similar to sympatric than to allopatric species?
  • Are mimic vocalizations (resembling sympatric species) more common in the focal bellbird than in other bellbirds?

Analysis flowchart

flowchart LR
  A[Read birdNET output] --> B(Format data) 
  B --> C(Graphs)
  C --> D(Statistical analysis)
  D --> E(Model summary) 

style A fill:#382A5433
style B fill:#395D9C33
style C fill:#3497A933
style D fill:#60CEAC33

Load packages and custom functions

# knitr is require for creating html/pdf/word reports
# formatR is used for soft-wrapping code

# install knitr package if not installed
if (!requireNamespace("sketchy", quietly = TRUE)) {
  install.packages("sketchy")
}


packages <- c(
  "knitr",
  "Rtsne",
  "caret",
  "ggplot2",
  "viridis",
  "doParallel",
  "brms",
  "brmsish",
  "cowplot",
  "Rraven",
  "kableExtra"
)

# install/ load packages
sketchy::load_packages(packages = packages)


options(knitr.kable.NA = '')

print <- function(x, row.names = FALSE) {
    kb <- kable(x, row.names = row.names, digits = 4, "html")
    kb <- kable_styling(kb,
                        bootstrap_options = c("striped", "hover", "condensed", "responsive"))
    scroll_box(kb, width = "100%")
}

common_names <- c(
  "Procnias tricarunculatus" = "Three-wattled Bellbird",
  "Psarocolius guatimozinus" = "Chestnut-headed Oropendola",
  "Psarocolius montezuma" = "Montezuma Oropendola",
  "Pteroglossus erythropygius" = "Pale-mandibled Aracari",
  "Pteroglossus torquatus" = "Collared Aracari",
  "Ramphastos brevis" = "Choco Toucan",
  "Ramphastos sulfuratus" = "Keel-billed Toucan",
  "nocall" = "Unclassified"
)

# to create several posterior predictive check plots out of a brms fit
custom_ppc <- function(fit, group = NULL, ndraws = 30) {
  by_group  <- if (!is.null(group)) {
    TRUE
  } else
    FALSE
  
  if (by_group)
    by_group  <-  if (any(names(fit$data) == group)) {
      TRUE
    } else
      FALSE
  
  if (by_group)
    by_group <-
      if (is.character(fit$data[, group]) |
          is.factor(fit$data[, group])) {
        TRUE
      } else
        FALSE
  
  
  if (by_group) {
    ppc_dens <- pp_check(fit,
                         ndraws = ndraws,
                         type = 'dens_overlay_grouped',
                         group = group)
    
    pp_mean <- pp_check(
      fit,
      type = "stat_grouped",
      stat = "mean",
      group = group,
      ndraws = ndraws
    )  + theme_classic()
    
  } else {
    ppc_dens <- pp_check(fit, ndraws = ndraws, type = 'dens_overlay')
    
    pp_mean <- pp_check(fit,
                        type = "stat",
                        stat = "mean",
                        ndraws = ndraws) + theme_classic()
    
  }
  
  pp_plot_list <-
    list(ppc_dens, pp_mean)
  
  pp_plot_list <-
    lapply(pp_plot_list, function(x)
      x   + scale_fill_viridis_d(
        begin = 0.3,
        end = 0.8,
        alpha = 0.5,
        option = "mako"
      ) + theme_classic())
  
  
  ppc_plot <- plot_grid(plotlist = pp_plot_list, ncol = 2)
 
   ppc_plot
}

1 Sympatric and allopatric species

1.1 BirdNET model

mod_eval_symp_allo <- read.csv("./data/processed/birdnet_models/sympatric_and_allopatric/sympatric_and_allopatric_evaluation.csv")

names(mod_eval_symp_allo) <- gsub("..", " ", names(mod_eval_symp_allo), fixed = TRUE)

print(mod_eval_symp_allo)

Class	Precision 0.5.	Recall 0.5.	F1.Score 0.5.	Precision opt.	Recall opt.	F1.Score opt.	AUPRC	AUROC	Optimal.Threshold	True.Positives	False.Positives	True.Negatives	False.Negatives	Samples	Percentage
OVERALL (Macro-avg)	0.8908	0.9627	0.9232	0.9180	0.9504	0.9312	0.9630	0.9968
Procnias tricarunculatus	0.9712	0.9854	0.9783	1.0000	0.9708	0.9852	0.9990	0.9998	0.80	133	0	884	4	137	13.42
Psarocolius guatimozinus	0.8909	0.9608	0.9245	0.9423	0.9608	0.9515	0.9867	0.9992	0.55	49	3	967	2	51	5.00
Psarocolius montezuma	0.9014	0.9846	0.9412	0.9538	0.9538	0.9538	0.9873	0.9991	0.75	62	3	953	3	65	6.37
Pteroglossus erythropygius	0.8943	0.9486	0.9206	0.9336	0.9206	0.9271	0.9844	0.9959	0.70	197	14	793	17	214	20.96
Pteroglossus torquatus	0.6216	0.8679	0.7244	0.6216	0.8679	0.7244	0.7844	0.9842	0.50	46	28	940	7	53	5.19
Ramphastos brevis	0.9840	0.9946	0.9892	0.9840	0.9946	0.9892	0.9995	0.9999	0.50	184	3	833	1	185	18.12
Ramphastos sulfuratus	0.9722	0.9968	0.9844	0.9904	0.9842	0.9873	0.9994	0.9997	0.80	311	3	702	5	316	30.95

Keeping only true positives (confidence >= optimal threshold) for the classification of vocalizations of sympatric and allopatric species

# read data
sympatric_allopatric_class <- imp_raven(
  path = "./data/processed/birdnet_classification_outputs/sympatrics_and_allopactrics/sympatrics_and_allopatrics_model/",
  all.data = TRUE,
  warbler.format = TRUE,
  pb = FALSE
)

# filter low confidence classifications
sympatric_allopatric_class$true.positive <- sapply(seq_len(nrow(sympatric_allopatric_class)), function(x){
    sympatric_allopatric_class$Confidence[x] >= mod_eval_symp_allo$Optimal.Threshold[mod_eval_symp_allo$Class == sympatric_allopatric_class$`Common Name`[x]]
    })

sympatric_allopatric_class <- sympatric_allopatric_class[sympatric_allopatric_class$true.positive, ]

# rename column
sympatric_allopatric_class$predicted_sp <- sympatric_allopatric_class$`Common Name`

sympatric_allopatric_class$`Common Name` <- NULL

# read metadata
est <- readRDS(
  "./data/processed/extended_selection_table_sympatric_allopatric_and_mimic_sounds.RDS"
)

# add metadata
sympatric_allopatric_class$sound.files <- basename(gsub("\\\\", "/", sympatric_allopatric_class$sound.files))


sympatric_allopatric_class$vocalizing_sp <- sapply(sympatric_allopatric_class$sound.files, function(x) {
  est$source[paste0(est$sound.files, ".wav") == x][1]
})


# save confusion matrix
conf_mat <- confusionMatrix(
  data = as.factor(sympatric_allopatric_class$predicted_sp),
  reference = as.factor(sympatric_allopatric_class$vocalizing_sp)
)

conf_df <- as.data.frame(conf_mat$table)

conf_df$total <- sapply(conf_df$Reference, function(x) {
  sum(sympatric_allopatric_class$vocalizing_sp == x)
})

conf_df$Proportion <- conf_df$Freq / conf_df$total

conf_df$pred_common_name <- common_names[as.character(conf_df$Prediction)]
conf_df$ref_common_name <- common_names[as.character(conf_df$Reference)]

# sort levels by genus
conf_df$pred_common_name <- factor(
  conf_df$pred_common_name,
  levels = c(
    "Three-wattled Bellbird",
    "Chestnut-headed Oropendola",
    "Montezuma Oropendola",
    "Collared Aracari",
    "Pale-mandibled Aracari",
    "Keel-billed Toucan",
    "Choco Toucan"
  )
)

conf_df$ref_common_name <- factor(
  conf_df$ref_common_name,
  levels = c(
    "Three-wattled Bellbird",
    "Chestnut-headed Oropendola",
    "Montezuma Oropendola",
    "Collared Aracari",
    "Pale-mandibled Aracari",
    "Keel-billed Toucan",
    "Choco Toucan"
  )
)

ggplot(conf_df, aes(x = ref_common_name, y = pred_common_name, fill = Proportion)) +
  geom_tile() +
  theme_bw() +
  coord_equal() +
  scale_fill_gradient2(low = mako(10)[3], high = mako(10)[7]) +
  labs(x = "Reference species", y = "Predicted species") +
  geom_text(aes(label = round(Proportion, 2)), color = "black") +
  theme_classic() +
  theme(
    axis.text.x = element_text(
      color = "black",
      size = 11,
      angle = 30,
      hjust = 1
    )
  )

1.1.1 Embeddings

embs <- read.csv("./data/processed/embeddings/embeddings_sympatric_allopatric_spp.csv")

# Split into separate columns
emb_matrix <- do.call(rbind, strsplit(embs$embedding, ","))
colnames(emb_matrix) <- paste0("emb", 1:ncol(emb_matrix))
emb_matrix <- apply(emb_matrix, 2, as.numeric)

# Combine with the first 3 columns
embs <- as.data.frame(cbind(embs[, 1:3], emb_matrix))

# add metadata
embs$file_path <- gsub("Nueva-grabaci\xf3n", "Nueva-grabación", embs$file_path, useBytes = TRUE)
embs$sound.files <- basename(gsub("\\\\", "/", embs$file_path, useBytes = TRUE))

embs$vocalizing_sp <- sapply(embs$sound.files, function(x) {
  est$source[paste0(est$sound.files, ".wav") == x][1]
})

set.seed(123)
tsne_out <- Rtsne(as.matrix(embs[, grep("emb", names(embs))]))

tsne_all <- data.frame(vocalizing_sp = embs$vocalizing_sp, tsne1 = tsne_out$Y[,1], tsne2 = tsne_out$Y[,2])

saveRDS(tsne_all, "./data/processed/tsne_dims_model_and_sister_species.rds")

# colors and shapes by species for scatterplots
cols <- c(adjustcolor("tomato", alpha.f = 0.4), mako(4, alpha = 0.4, begin = 0.2, end = 0.8))
cols2 <- c(adjustcolor("tomato", alpha.f = 0.9), mako(4, alpha = 0.9, begin = 0.2, end = 0.8))

cols <- c(cols[1:2], rep(cols[3:5], each = 2))
cols2 <- c(cols2[1:2], rep(cols2[3:5], each = 2))

shapes <- c(16, 16, 16, 1, 17, 2, 15, 0)

names(shapes) <- names(cols) <- names(cols2) <- c(
  "mimic_bellbird",
  "Three-wattled Bellbird",
  "Montezuma Oropendola",
  "Chestnut-headed Oropendola",
  "Collared Aracari",
  "Pale-mandibled Aracari",
  "Keel-billed Toucan",
  "Choco Toucan"
)

tsne_all <- readRDS("./data/processed/tsne_dims_model_and_sister_species.rds")

tsne_only_model_species <- tsne_all[tsne_all$vocalizing_sp != "mimic_bellbird", ]

tsne_only_model_species$common_name <- common_names[as.character(tsne_only_model_species$vocalizing_sp)] 

tsne_only_model_species$common_name <- factor(
  tsne_only_model_species$common_name,
  levels = c(
    "Three-wattled Bellbird",
    "Chestnut-headed Oropendola",
    "Montezuma Oropendola",
    "Collared Aracari",
    "Pale-mandibled Aracari",
    "Keel-billed Toucan",
    "Choco Toucan"
  )
)

ggplot(tsne_only_model_species, aes(x = tsne1, y = tsne2, color = common_name, shape = common_name)) +
  geom_point(size = 2) +
  theme_minimal() +
  labs(x = "t-SNE 1", y = "t-SNE 2", color = "", shape = "") +
    scale_color_manual(values = cols[names(cols) %in% unique(tsne_only_model_species$common_name)]) +
    scale_shape_manual(values = shapes[names(shapes) %in% unique(tsne_only_model_species$common_name)]) +
    theme_classic()

1.1.2 Focal bellbird classification

# read data
mimic_class <- imp_raven(
  path = "./data/processed/birdnet_classification_outputs/mimic_bellbird/sympatrics_and_allopatrics_model/",
  all.data = TRUE,
  warbler.format = TRUE,
  pb = FALSE
)

# filter low confidence classifications
mimic_class$true.positive <- sapply(seq_len(nrow(mimic_class)), function(x){
    mimic_class$Confidence[x] >= mod_eval_symp_allo$Optimal.Threshold[mod_eval_symp_allo$Class == mimic_class$`Common Name`[x]]
    })

mimic_class <- mimic_class[mimic_class$true.positive, ]

# add metadata
mimic_class$sound.files <- basename(gsub(
  "\\\\",
  "/",
  mimic_class$sound.files
))

# rename column
mimic_class$predicted_sp <- mimic_class$`Common Name`
mimic_class$vocalizing_sp <- "Procnias tricarunculatus"

mimic_class$`Common Name` <- NULL

prop_mimic_class <- table(mimic_class$predicted_sp)

prop_mimic_class <- data.frame(prop_mimic_class)
names(prop_mimic_class) <- c("Prediction", "Freq")
prop_mimic_class$total <- sum(prop_mimic_class$Freq)
prop_mimic_class$Proportion <- prop_mimic_class$Freq / prop_mimic_class$total
prop_mimic_class$Reference <- "Focal\nbellbird"

others <- conf_df[conf_df$Reference == "Procnias tricarunculatus", ]
others$Reference <- "Other\nbellbirds"

# proportions observed in adult bellbird
prop_class <- rbind(
  prop_mimic_class,
  others[, c("Prediction", "Freq", "total", "Proportion", "Reference")]
)

prop_class$genus <- sapply(
  strsplit(as.character(prop_class$Prediction), " "),
  "[[",
  1
)

prop_class$type <- ifelse(
  grepl("guati|eryth|brevis", prop_class$Prediction),
  "Allopatric",
  "Sympatric"
)

prop_class$type[prop_class$genus == "Procnias"] <- "Own\nspecies"

prop_class$type <- factor(
  prop_class$type,
  levels = c("Allopatric", "Sympatric", "Own\nspecies")
)

prop_class$genus <- factor(
  prop_class$genus,
  levels = c("Procnias", "Pteroglossus", "Psarocolius", "Ramphastos")
)

prop_class$Reference <- factor(
  prop_class$Reference,
  levels = c("Other\nbellbirds", "Focal\nbellbird")
)


genus_labels <- c(
  "Procnias" = "Three-wattled bellbird",
  "Psarocolius" = "Oropendolas",
  "Pteroglossus" = "Aracaris",
  "Ramphastos" = "Toucans"
)

prop_class$common_name <- common_names[as.character(prop_class$Prediction)]

prop_class$common_name <- gsub("-", "\n", prop_class$common_name)


prop_class$sp_label <- sapply(
  strsplit(as.character(prop_class$common_name), " "),
  "[[",
  1
)

prop_class$sp_label[prop_class$sp_label == "Three\nwattled"] <- ""


prop_class$sp_label_y <- ave(
  prop_class$Proportion,
  prop_class$genus,
  FUN = function(x) {
    x + log10(max(x) + 0.999) * 0.44
  }
)


bar_cols <- c("gray", mako(10, alpha = 0.6, begin = 0.2, end = 0.8)[c(2, 8)])

# dummy rows controlling y-range of each facet
dummy_df <- data.frame(
  genus = c("Procnias", "Psarocolius", "Pteroglossus", "Ramphastos"),
  Reference = "Focal\nbellbird",
  Proportion = c(1.18, 0.0024, 0.17, 0.50)
)

dummy_df$genus <- factor(
  dummy_df$genus,
  levels = c("Procnias", "Pteroglossus", "Psarocolius", "Ramphastos")
)


ggplot(prop_class, aes(x = Reference, y = Proportion)) +

  # invisible layer used only to expand scales
  geom_blank(data = dummy_df) +

  geom_bar(
    aes(group = interaction(Reference, type), fill = type),
    position = position_dodge(),
    stat = "identity"
  ) +

  geom_text(
    aes(
      group = interaction(Reference, type),
      label = paste0(round(Proportion * 100, 1), "%")
    ),
    position = position_dodge(width = 0.9),
    vjust = -0.5,
    size = 3
  ) +

  facet_wrap(
    . ~ genus,
    scales = "free_y",
    labeller = labeller(genus = genus_labels)
  ) +

  ylab("Percentage of vocalizations predicted as") +
  xlab("") +

  scale_y_continuous(
    labels = scales::percent,
    expand = expansion(mult = c(0, 0))
  ) +

  scale_fill_manual(
    name = "",
    values = c(
      "Own\nspecies" = bar_cols[1],
      "Allopatric" = bar_cols[2],
      "Sympatric" = bar_cols[3]
    ),
    breaks = c("Allopatric", "Sympatric")
  ) +

  theme_classic() +
  theme(
    plot.title = element_text(hjust = 0.5),
    legend.position = c(0.74, 0.78),
    legend.justification = c(1, 0),
    legend.background = element_rect(fill = "white", color = NA),
    panel.grid.major.x = element_blank()
  ) +

  geom_text(
    aes(
      group = interaction(Reference, type),
      label = sp_label,
      y = sp_label_y
    ),
    position = position_dodge(width = 0.9),
    vjust = 0.5,
    size = 2.5,
    angle = 90,
    hjust = 0
  )

1.1.3 Statistical analysis

Evaluate whether the proportion of vocalizations predicted as sympatric (model) species vs allopatric species is significantly higher in the focal bellbird than in other bellbirds

mimic_other_bellbirds <- sympatric_allopatric_class[sympatric_allopatric_class$predicted_sp != "Procnias tricarunculatus", ]

mimic_focal <- mimic_class[mimic_class$predicted_sp != "Procnias tricarunculatus", ]

mimic_focal$vocalizing_type <- "focal"
mimic_other_bellbirds$vocalizing_type <- "others"

mimic_bellbirds <- rbind(
  mimic_focal[, c("predicted_sp", "vocalizing_type")],
  mimic_other_bellbirds[, c("predicted_sp", "vocalizing_type")]
)

mimic_bellbirds$vocalizing_type <- factor(
  mimic_bellbirds$vocalizing_type,
  levels = c("others", "focal")
)

mimic_bellbirds$mimic_type <- 
ifelse(
  grepl("guati|eryth|brevis", mimic_bellbirds$predicted_sp),
  0,
  1
)

mimic_bellbirds$genus <- sapply(
  strsplit(as.character(mimic_bellbirds$predicted_sp), " "),
  "[[",
  1
)

cores <- 4
chains <- 4
iter <- 4000

priors <- c(
  prior(normal(0, 1.5), class = "Intercept"),
  prior(normal(0, 0.7), class = "b"),
  prior(exponential(1), class = "sd")
)

only_prior_fit  <- brm(
  mimic_type ~ vocalizing_type + (1 | genus),
  data = mimic_bellbirds,
  family = bernoulli(),
  prior = priors,
  sample_prior = "only",
  chains = chains,
  iter = iter,
  cores = cores,
          file = "./data/processed/mimicry_prob_by_simpatry_model_only_priors.RDS", 
  file_refit = "on_change"
)


mod1 <- brm(
  mimic_type ~ vocalizing_type  + (1 | genus),
  data = mimic_bellbirds,
  family = bernoulli(),
  cores = cores,
    prior = priors,
        chains = chains,
        iter = iter,
        file = "./data/processed/mimicry_prob_by_simpatry_model.RDS", 
  file_refit = "on_change"
)

1.1.3.1 Prior check

Priors

prior_summary(mod1)

prior	class	coef	group	resp	dpar	nlpar	lb	ub	tag	source
normal(0, 0.7)	b									user
	b	vocalizing_typefocal								default
normal(0, 1.5)	Intercept									user
exponential(1)	sd						0			user
	sd		genus							default
	sd	Intercept	genus							default

Prior predictive checks

custom_ppc(only_prior_fit, group = "vocalizing_type", ndraws = 100)

1.1.3.2 Fit summary

fill_color <- mako(10)[8]

extended_summary(
  mod1,
  n.posterior = 1000,
  fill = fill_color,
  trace.palette = mako,
  remove.intercepts = FALSE,
  highlight = TRUE,
  print.name = FALSE
)

	formula	family	priors	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	mimic_type ~ vocalizing_type + (1 | genus)	bernoulli (logit)	b-normal(0, 0.7) Intercept-normal(0, 1.5) sd-exponential(1)	4000	4	1	2000	12 (0.0015%)	0	1519.44	1939.21	1537639093

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-0.274	-1.405	0.846	1.003	1519.44	1939.21
b_vocalizing_typefocal	2.966	2.641	3.307	1.003	3701.09	3564.68

Posterior predictive check

custom_ppc(mod1, group = "vocalizing_type", ndraws = 100)

post <- as_draws_df(mod1)

# Odds ratio: focal vs other birds
focal_vs_others <- exp(post$b_vocalizing_typefocal)

# Predicted probability of sympatric mimicry in focal birds
p_focal <- posterior_epred(
  mod1,
  newdata = data.frame(
    vocalizing_type = "focal",
    genus = unique(mimic_bellbirds$genus)[1]
  ),
  re_formula = NA
)

p_focal <- as.numeric(p_focal)

# Predicted probability of sympatric mimicry in other birds
p_others <- posterior_epred(
  mod1,
  newdata = data.frame(
    vocalizing_type = "others",
    genus = unique(mimic_bellbirds$genus)[1]
  ),
  re_formula = NA
)

p_others <- as.numeric(p_others)

# Within-group sympatric vs allopatric odds
focal_ratio <- p_focal / (1 - p_focal)
other_ratio <- p_others / (1 - p_others)

1.1.3.3 Summary

• Focal bellbird vocalizations had 19.32 times greater odds [95% UI: 14.03–27.3] of being classified as the sympatric rather than the allopatric model species compared with vocalizations from other bellbirds (predicted probability of sympatric classification for focal birds = 0.94; for other bellbirds = 0.43).

• Within focal bellbird vocalizations, the odds of classification as the sympatric model species were 14.9 times greater [95% UI: 4.56–47.88] than the odds of classification as closely related allopatric species (predicted probability of sympatric classification for focal birds = 0.94; predicted probability of allopatric classification for focal birds = 0.06).

• For vocalizations from other bellbirds, the odds of classification as sympatric and allopatric species were similar (odds = 0.77 [95% UI: 0.25–2.33]; predicted probability of sympatric classification for other bellbirds = 0.43; predicted probability of allopatric classification for other bellbirds = 0.57).

2 Only sympatric species

2.1 BirdNET model

mod_eval_symp <- read.csv("./data/processed/birdnet_models/sympatric/sympatric_evaluation.csv")

names(mod_eval_symp) <- gsub("..", " ", names(mod_eval_symp), fixed = TRUE)

print(mod_eval_symp)

Class	Precision 0.5.	Recall 0.5.	F1.Score 0.5.	Precision opt.	Recall opt.	F1.Score opt.	AUPRC	AUROC	Optimal.Threshold	True.Positives	False.Positives	True.Negatives	False.Negatives	Samples	Percentage
OVERALL (Macro-avg)	0.9884	0.9908	0.9896	0.9939	0.9890	0.9914	0.9994	0.9998
Procnias tricarunculatus	0.9783	0.9854	0.9818	1.0000	0.9781	0.9889	0.9994	0.9998	0.60	134	0	434	3	137	23.99
Psarocolius montezuma	0.9848	1.0000	0.9924	0.9848	1.0000	0.9924	0.9985	0.9998	0.15	65	1	505	0	65	11.38
Pteroglossus torquatus	1.0000	0.9811	0.9905	1.0000	0.9811	0.9905	0.9997	1.0000	0.15	52	0	518	1	53	9.28
Ramphastos sulfuratus	0.9906	0.9968	0.9937	0.9906	0.9968	0.9937	0.9998	0.9998	0.50	315	3	252	1	316	55.34

# read data
sympatric_class <- imp_raven(
  path = "./data/processed/birdnet_classification_outputs/sympatrics/sympatrics_model/",
  all.data = TRUE,
  warbler.format = TRUE,
  pb = FALSE
)

sympatric_class$true.positive <- sapply(seq_len(nrow(sympatric_class)), function(x){
    sympatric_class$Confidence[x] >= mod_eval_symp$Optimal.Threshold[mod_eval_symp$Class == sympatric_class$`Common Name`[x]]
    })
    

sympatric_class <- sympatric_class[sympatric_class$true.positive,]

# rename column
sympatric_class$predicted_sp <- sympatric_class$`Common Name`

sympatric_class$`Common Name` <- NULL

# read metadata
est <- readRDS(
  "./data/processed/extended_selection_table_sympatric_allopatric_and_mimic_sounds.RDS"
)

# add metadata
sympatric_class$sound.files <- basename(gsub("\\\\", "/", sympatric_class$sound.files))


sympatric_class$vocalizing_sp <- sapply(sympatric_class$sound.files, function(x) {
  est$source[paste0(est$sound.files, ".wav") == x][1]
})


sympatric_class <- sympatric_class[sympatric_class$vocalizing_sp %in% unique(sympatric_class$predicted_sp),]

# save confusion matrix
conf_mat <- confusionMatrix(
  data = as.factor(sympatric_class$predicted_sp),
  reference = as.factor(sympatric_class$vocalizing_sp)
)

conf_df <- as.data.frame(conf_mat$table)

conf_df$total <- sapply(conf_df$Reference, function(x) {
  sum(sympatric_class$vocalizing_sp == x)
})

conf_df$Proportion <- conf_df$Freq / conf_df$total

conf_df$pred_common_name <- common_names[as.character(conf_df$Prediction)]
conf_df$ref_common_name <- common_names[as.character(conf_df$Reference)]

# sort levels by genus
conf_df$pred_common_name <- factor(
  conf_df$pred_common_name,
  levels = c(
    "Three-wattled Bellbird",
    "Montezuma Oropendola",
    "Collared Aracari",
    "Keel-billed Toucan"
  )
)

conf_df$ref_common_name <- factor(
  conf_df$ref_common_name,
  levels = c(
    "Three-wattled Bellbird",
    "Montezuma Oropendola",
    "Collared Aracari",
    "Keel-billed Toucan"
  )
)

ggplot(conf_df, aes(x = ref_common_name, y = pred_common_name, fill = Proportion)) +
  geom_tile() +
  theme_bw() +
  coord_equal() +
  scale_fill_gradient2(low = mako(10)[3], high = mako(10)[7]) +
  labs(x = "Reference species", y = "Predicted species") +
  geom_text(aes(label = round(Proportion, 2)), color = "black") +
  theme_classic() +
  theme(
    axis.text.x = element_text(
      color = "black",
      size = 11,
      angle = 30,
      hjust = 1
    )
  )

# read data
sympatric_focal <- imp_raven(
  path = "./data/processed/birdnet_classification_outputs/mimic_bellbird/sympatric_model/",
  all.data = TRUE,
  warbler.format = TRUE,
  pb = FALSE
)



# filter out low confidence calls
sympatric_focal$true.positive <- sapply(seq_len(nrow(sympatric_focal)), function(x){
    if (sympatric_focal$`Common Name`[x] == "nocall")
        return(NA)
        
    sympatric_focal$Confidence[x] >= mod_eval_symp$Optimal.Threshold[mod_eval_symp$Class == sympatric_focal$`Common Name`[x]]
    })
    


# remove nocalls
sympatric_focal_nocalls <- sympatric_focal

sympatric_focal <- sympatric_focal[sympatric_focal$true.positive | is.na(sympatric_focal$true.positive), ]


sympatric_focal$Reference <- "Focal\nbellbird"

# rename column
sympatric_focal$predicted_sp <- sympatric_focal$`Common Name`

sympatric_focal$`Common Name` <- NULL

# get other bellbirds from previous data
sympatric_others <- sympatric_class[sympatric_class$vocalizing_sp == "Procnias tricarunculatus", ]

sympatric_others$vocalizing_sp <- NULL

sympatric_others$Reference <- "Other\nbellbirds"

sympatric <- rbind(
  sympatric_focal,
  sympatric_others
)

# # add metadata
# sympatric$sound.files <- basename(gsub(
#   "\\\\",
#   "/",
#   sympatric$sound.files
# ))
# 
# # rename column
# sympatric$predicted_sp <- sympatric$`Common Name`
# sympatric$`Common Name` <- NULL

prop_sympatric <- table(sympatric$predicted_sp, sympatric$Reference)

prop_sympatric <- data.frame(prop_sympatric)
names(prop_sympatric) <- c("Prediction", "Vocalizing", "Freq")

# add Psarocolius
# prop_sympatric2 <- prop_sympatric[3:4, ]
# prop_sympatric2$Prediction <- "Psarocolius montezuma"
# prop_sympatric2$Freq <- 0

# prop_sympatric <- rbind(prop_sympatric, prop_sympatric2)

prop_sympatric$total <- sapply(prop_sympatric$Vocalizing, function(x) sum(prop_sympatric$Freq[prop_sympatric$Vocalizing == x]))

prop_sympatric$Proportion <- prop_sympatric$Freq / prop_sympatric$total

prop_sympatric$common_name <- common_names[as.character(
  prop_sympatric$Prediction
)]

# order levels
prop_sympatric$common_name <- factor(
  prop_sympatric$common_name,
  levels = c(
    "Three-wattled Bellbird",
    "Collared Aracari",
    "Keel-billed Toucan",
    "Montezuma Oropendola"
  )
)


# ggplot2 bar plot single panel colored by references
ggplot(prop_sympatric[prop_sympatric$Prediction != "nocall",], aes(x = common_name, y = Proportion, fill = Vocalizing)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9)) +

  geom_text(
    aes(label = paste0(round(Proportion * 100, 1), "%")),
    position = position_dodge(width = 0.9),
    vjust = -0.4,
    size = 4
  ) +

  labs(y = "Percentage of vocalizations predicted as", x = "", fill = "") +

  scale_y_continuous(
    labels = scales::percent,
    expand = expansion(mult = c(0, 0.16))
  ) +

  scale_fill_manual(values = mako(10, alpha = 0.6)[c(3, 8)]) +

  theme_classic() +
  theme(
    legend.position = c(0.1, 0.87),
    plot.title = element_text(hjust = 0.5),
    panel.grid.major.x = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank()
  ) +
  facet_wrap(. ~ common_name, nrow = 2, scales = "free")

2.1.1 Statistical analysis

sympatric$prediction_type <- ifelse(
  grepl("rocnias", sympatric$predicted_sp),
  0,
  1
)

sympatric$vocalizing_type <- ifelse(
  grepl("ocal", sympatric$Reference),
  "focal",
  "others"
)

# order levels
sympatric$vocalizing_type <- factor(
  sympatric$vocalizing_type,
  levels = c("others", "focal")
)

cores <- 4
chains <- 4
iter <- 4000

priors <- c(
  prior(normal(-3, 1), class = "Intercept"),
  prior(normal(0, 2), class = "b")
)

only_prior_fit2  <- brm(
  prediction_type ~ vocalizing_type + predicted_sp,
  data = sympatric[sympatric$predicted_sp != "nocall",],
  family = bernoulli(),
  prior = priors,
sample_prior = "only",
  chains = chains,
  iter = iter,
  cores = cores,
  file = "./data/processed/mimicry_simpatric_only_model_only_priors.RDS",
  file_refit = "on_change"
)


mod2 <- brm(
  prediction_type ~ vocalizing_type,
  data = sympatric[sympatric$predicted_sp != "nocall",],
  family = bernoulli(),
  prior = priors,
  chains = chains,
  iter = iter,
  cores = cores,
  control = list(
  adapt_delta = 0.999,
  max_treedepth = 15
),
  file = "./data/processed/mimicry_simpatric_only_model.RDS",
  file_refit = "on_change"
)

2.1.1.1 Prior check

Priors

prior_summary(mod2)

prior	class	coef	group	resp	dpar	nlpar	lb	ub	tag	source
normal(0, 2)	b									user
	b	vocalizing_typefocal								default
normal(-3, 1)	Intercept									user

Prior predictive checks

custom_ppc(only_prior_fit2, group = "vocalizing_type", ndraws = 100)

2.1.2 Fit summary

fill_color <- mako(10)[8]

extended_summary(
  mod2,
  n.posterior = 1000,
  fill = fill_color,
  trace.palette = mako,
  remove.intercepts = FALSE,
  highlight = TRUE,
  print.name = FALSE
)

	formula	family	priors	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	prediction_type ~ vocalizing_type	bernoulli (logit)	b-normal(0, 2) Intercept-normal(-3, 1)	4000	4	1	2000	0 (0%)	0	1390.19	1412.9	549624481

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-4.814	-5.603	-4.128	1.003	1390.19	1423.82
b_vocalizing_typefocal	4.988	4.287	5.786	1.003	1422.15	1412.90

Posterior predictive check

custom_ppc(mod1, group = "vocalizing_type", ndraws = 100)

# Odds ratio: focal vs other bellbirds
post2 <- as_draws_df(mod2)

# Odds ratio: focal vs other bellbirds
focal_vs_others <- exp(post2$b_vocalizing_typefocal)

# Predicted probability of sympatric classification for focal birds
p_focal <- posterior_epred(
  mod2,
  newdata = data.frame(vocalizing_type = "focal")
)

p_focal <- as.numeric(p_focal)

# Predicted probability of sympatric classification for other bellbirds
p_others <- posterior_epred(
  mod2,
  newdata = data.frame(vocalizing_type = "others")
)

p_others <- as.numeric(p_others)

2.1.2.1 Summary

• Focal bellbird vocalizations had 143.81 times greater odds [95% UI: 72.75–325.59] of being classified as a sympatric model species compared with vocalizations from other bellbirds (predicted probability of sympatric classification for focal birds = 0.54; for other bellbirds = 0.01).

3 Mimicry through time

mimic_data <- est[est$source == "mimetic_bellbird", ]

rec_info <- read.table("./data/raw/Recording_dates_.csv", head = TRUE, sep = ";")

rec_info$date <- gsub("\\/", "-", rec_info$date)
rec_info$date <- gsub("out", "ago", rec_info$date)
rec_info$date <- gsub("set", "sep", rec_info$date)

rec_info$org.date <- rec_info$date
rec_info$date <- as.Date(rec_info$date, format = "%d-%b-%y")


mimic_data$orig.sound.files <- attr(mimic_data, "check.res")$orig.sound.files


mimic_data$date <- sapply(mimic_data$orig.sound.files, function(x) as.character(rec_info$date[rec_info$sound.files ==
    x][1]))

mimic_data$date <- as.Date(mimic_data$date)

sympatric_focal$fixed.sound.files <- gsub("Nueva-grabaci\xf3n", "Nueva-grabación", sympatric_focal$sound.files, useBytes = TRUE)

sympatric_focal$fixed.sound.files <- basename(gsub("\\\\", "/", sympatric_focal$fixed.sound.files, useBytes = TRUE))

mimic_data$sound.files <- gsub("Nueva-grabaci\xf3n", "Nueva-grabación", mimic_data$sound.files, useBytes = TRUE)


mimic_data$sound.files <- paste0(mimic_data$sound.files, ".wav")

mimic_data$pred_mimic <- sapply(mimic_data$sound.files, function(x) {
  sympatric_focal$predicted_sp[sympatric_focal$fixed.sound.files == x][1]
})

mimic_data <- mimic_data[mimic_data$pred_mimic != "babbling",]

mimic_data$pred_mimic <- ifelse(mimic_data$pred_mimic == "babbling" | is.na(mimic_data$pred_mimic), "nocall", mimic_data$pred_mimic)

agg_pred <- aggregate(orig.sound.files ~ date + pred_mimic, mimic_data, length)


out <- lapply(unique(agg_pred$date), function(x) {

    X <- agg_pred[agg_pred$date == x, ]
    X$prop <- X$orig.sound.files/sum(X$orig.sound.files)

    if (length(unique(X$pred_mimic)) < 4) {

        new_df <- data.frame(date = X$date[1], pred_mimic = setdiff(c("Procnias tricarunculatus",
            "Pteroglossus torquatus", "Ramphastos sulfuratus", "Psarocolius montezuma", "nocall"), X$pred_mimic),
            orig.sound.files = 0, prop = 0)
  
        X <- rbind(new_df, X)
    }
    
      X$total <- sum(X$orig.sound.files)
    return(X)
})

prop_pred <- do.call(rbind, out)

prop_pred$date.fact <- as.factor(as.character(prop_pred$date))

count_date <- aggregate(orig.sound.files ~ date, prop_pred, sum)

prop_pred$pred_common_name <- common_names[as.character(prop_pred$pred_mimic)]


ggplot(prop_pred[prop_pred$date %in% count_date$date[count_date$orig.sound.files >=
    10], ], aes(x = date.fact, y = prop, fill = pred_common_name)) + geom_bar(stat = "identity") +
    scale_fill_viridis_d(alpha = 0.7, end = 0.8, option = "G") + theme_classic() + theme(axis.text.x = element_text(angle = 45,
    vjust = 0.5)) + labs(x = "Date", y = "Proportion of vocalizations\npredicted as")

3.1 Statistical analysis

We modeled the probability of producing original sounds as a function of Julian day while allowing both intercepts and temporal slopes to vary among mimic categories (pred_mimic). Rather than fitting a full fixed interaction between date and mimic category, we used a hierarchical varying-slope structure to account for heterogeneity among mimic groups while enabling partial pooling across categories. This approach reduces overparameterization and yields more stable estimates, particularly for mimic groups with limited sample sizes, while still allowing different groups to exhibit distinct temporal trajectories.

# weakly informative priors
priors <- c(
# Weakly informative logistic intercept
prior(normal(0, 2.5), class = "Intercept"),
# Temporal slope
prior(normal(0, 1), class = "b"),
# Random effect standard deviations
prior(exponential(1), class = "sd")
)

dat <- prop_pred[prop_pred$pred_mimic != "Psarocolius montezuma",]

dat$julian_day <- as.numeric(format(as.Date(dat$date), "%j"))


dat$pred_common_name <- common_names[as.character(dat$pred_mimic)]


dat$pred_common_name <- factor(dat$pred_common_name, levels = c( "Three-wattled Bellbird", "Collared Aracari",     "Keel-billed Toucan", "Unclassified"))


mod2 <- brm(
  orig.sound.files | trials(total) ~ 
    scale(julian_day) +
    (1 + scale(julian_day) | pred_common_name),
  family = binomial,
  data = dat,
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 4000,
  control = list(adapt_delta = 0.95), file_refit = "on_change",
  file = "data/processed/brms_model_prop_pred_mimicry.rds"
)

# mod2 <- add_criterion(
#       mod2,
#       criterion = "loo")


only_prior_mod2 <- brm(
  orig.sound.files | trials(total) ~ 
    scale(julian_day) +
    (1 + scale(julian_day) | pred_common_name),
  family = binomial,
  data = dat,
  prior = priors,
  sample_prior = "only",
  chains = 4,
  cores = 4,
  iter = 4000,
  control = list(adapt_delta = 0.95), file_refit = "on_change",
  file = "data/processed/brms_model_prop_pred_mimicry_only_prior.rds"
)

3.1.1 Prior check

Priors

prior_summary(mod2)

prior	class	coef	group	resp	dpar	nlpar	lb	ub	tag	source
normal(0, 1)	b									user
	b	scalejulian_day								default
normal(0, 2.5)	Intercept									user
lkj_corr_cholesky(1)	L									default
	L		pred_common_name							default
exponential(1)	sd						0			user
	sd		pred_common_name							default
	sd	Intercept	pred_common_name							default
	sd	scalejulian_day	pred_common_name							default

Prior predictive checks

custom_ppc(only_prior_mod2, group = "vocalizing_type", ndraws = 100)

3.1.2 Fit summary

fill_color <- mako(10)[8]

extended_summary(
  mod2,
  n.posterior = 1000,
  fill = fill_color,
  trace.palette = mako,
  remove.intercepts = FALSE,
  highlight = TRUE,
  print.name = FALSE
)

	formula	family	priors	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	orig.sound.files | trials(total) ~ scale(julian_day) + (1 + scale(julian_day) | pred_common_name)	binomial (logit)	b-normal(0, 1) Intercept-normal(0, 2.5) L-lkj_corr_cholesky(1) sd-exponential(1)	4000	4	1	2000	2 (0.00025%)	0	2549.73	2900.21	1896986874

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-1.204	-2.185	-0.153	1.001	2549.73	2900.21
b_scalejulian_day	-0.046	-0.967	0.909	1.001	3452.85	3753.83

Posterior predictive check

custom_ppc(mod2, group = "vocalizing_type", ndraws = 100)

# Extract posterior draws as base R data.frame
post <- as.data.frame(mod2)

# Get species levels
species_levels <- unique(mod2$data$pred_common_name)

# Global slope
beta_global <- post[, "b_scalejulian_day"]
# Find all random slope columns automatically
rand_cols <- grep("^r_pred_common_name\\[.*,scalejulian_day\\]$",
names(post),
value = TRUE)
# Extract species names from column names
species_levels <- sub("r_pred_common_name\\[(.*),scalejulian_day\\]", "\\1", rand_cols)


# Storage matrix
results <- matrix(NA, nrow = length(species_levels), ncol = 3)
colnames(results) <- c("mean_slope",
                       "l95",
                       "u95")

rownames(results) <- species_levels


# Loop through species
for (i in seq_along(rand_cols)) {

  slope_draws <- beta_global + post[, rand_cols[i]]

  results[i, "mean_slope"]    <- mean(slope_draws)
  results[i, "l95"]           <- quantile(slope_draws, 0.025)
  results[i, "u95"]           <- quantile(slope_draws, 0.975)
}

# Print results
print(round(results, 4), row.names = TRUE)

	mean_slope	l95	u95
Three-wattled.Bellbird	0.9856	0.8240	1.1513
Collared.Aracari	-0.0238	-0.1982	0.1574
Keel-billed.Toucan	-0.8677	-1.0282	-0.7091
Unclassified	-0.3365	-0.5525	-0.1239

# dat <- prop_pred[prop_pred$pred_mimic != "Psarocolius montezuma",]
# 
# dat$julian_day <- as.numeric(format(as.Date(dat$date), "%j"))


# Create prediction grid
newdata <- expand.grid(
  julian_day = seq(min(dat$julian_day),
                   max(dat$julian_day),
                   length.out = 100),
  pred_common_name = unique(dat$pred_common_name)
)

# IMPORTANT: match scaling used in model
newdata$scalejulian_day <- scale(dat$julian_day)[1]  # dummy to get attributes
newdata$scalejulian_day <- (newdata$julian_day - mean(dat$julian_day)) /
                           sd(dat$julian_day)

newdata$total <- round(mean(dat$total), 0)  # use mean total for predictions

# Get posterior expected probabilities
epred <- posterior_epred(
  mod2,
  newdata = newdata,
  re_formula = NULL
)

# Posterior mean prediction
newdata$fit <- colMeans(epred)

# Convert julian day back to date for plotting
origin_year <- format(as.Date(dat$date[1]), "%Y")
newdata$date <- as.Date(newdata$julian_day - 1,
                        origin = paste0(origin_year, "-01-01"))

newdata$fit_prop <- newdata$fit / 100

# include only those for which CI did not include 0
results_df <- as.data.frame(results)
sig_spp <- rownames(results_df)[results_df$l95 * results_df$u95 > 0]

newdata <- newdata[newdata$pred_common_name %in% c("Three-wattled Bellbird", "Collared Aracari",     "Keel-billed Toucan", "Unclassified"), ]


# Plot
ggplot(
  dat[dat$prop > 0, ],
  aes(
    x = as.Date(date),
    y = prop,
    fill = pred_common_name,
    col = pred_common_name,
    size = orig.sound.files
  )
) +
  geom_point(pch = 21) +
  geom_line(
    data = newdata,
    aes(x = date, y = fit_prop, color = pred_common_name),
    size = 1.2
  ) +
  scale_fill_viridis_d(alpha = 0.6, end = 0.8) +
  scale_color_viridis_d(alpha = 0.8, end = 0.8) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5)) +
  labs(
    x = "Date",
    y = "Proportion of vocalizations",
    fill = "Predicted species",
    color = "Predicted species",
    size = "Number of recorded\nvocalizations"
  )

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

3.1.3 Summary

• The proportion of vocalizations of its own species increased true time while the proportion for Killed-billed Toucan and unclassified calls decreased
• Proportions for Collared Araraci did not show a definitive pattern

Session information

─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.5.2 (2025-10-31)
 os       Ubuntu 22.04.5 LTS
 system   x86_64, linux-gnu
 ui       X11
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       America/Costa_Rica
 date     2026-05-21
 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)
 bitops           1.0-9      2024-10-03 [3] CRAN (R 4.5.0)
 bridgesampling   1.2-1      2025-11-19 [1] CRAN (R 4.5.0)
 brio             1.1.5      2024-04-24 [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)
 caret          * 7.0-1      2024-12-10 [1] CRAN (R 4.5.0)
 checkmate        2.3.4      2026-02-03 [1] CRAN (R 4.5.2)
 class            7.3-23     2025-01-01 [4] CRAN (R 4.4.2)
 cli              3.6.6      2026-04-09 [1] CRAN (R 4.5.2)
 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)
 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)
 data.table       1.18.2.1   2026-01-27 [1] CRAN (R 4.5.2)
 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)
 doParallel     * 1.0.17     2022-02-07 [1] CRAN (R 4.5.0)
 dplyr            1.1.4      2023-11-17 [1] CRAN (R 4.5.0)
 dtw              1.23-1     2022-09-19 [1] CRAN (R 4.5.0)
 e1071            1.7-17     2025-12-18 [1] CRAN (R 4.5.2)
 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)
 fftw             1.0-9      2024-09-20 [3] CRAN (R 4.5.0)
 foreach        * 1.5.2      2022-02-02 [3] CRAN (R 4.1.2)
 fs               2.1.0      2026-04-18 [1] CRAN (R 4.5.2)
 future           1.70.0     2026-03-14 [1] CRAN (R 4.5.2)
 future.apply     1.20.2     2026-02-20 [1] CRAN (R 4.5.2)
 generics         0.1.4      2025-05-09 [1] CRAN (R 4.5.0)
 ggdist           3.3.3      2025-04-23 [1] CRAN (R 4.5.0)
 ggplot2        * 4.0.3      2026-04-22 [1] CRAN (R 4.5.2)
 globals          0.19.1     2026-03-13 [1] CRAN (R 4.5.2)
 glue             1.8.1      2026-04-17 [1] CRAN (R 4.5.2)
 gower            1.0.2      2024-12-17 [3] CRAN (R 4.5.0)
 gridExtra        2.3        2017-09-09 [1] CRAN (R 4.5.0)
 gtable           0.3.6      2024-10-25 [1] CRAN (R 4.5.0)
 hardhat          1.4.1      2025-01-31 [1] 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)
 httr             1.4.7      2023-08-15 [1] CRAN (R 4.5.0)
 inline           0.3.21     2025-01-09 [1] CRAN (R 4.5.0)
 ipred            0.9-15     2024-07-18 [3] CRAN (R 4.5.0)
 iterators      * 1.0.14     2022-02-05 [3] CRAN (R 4.1.2)
 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)
 lava             1.9.0      2026-04-05 [1] CRAN (R 4.5.2)
 lifecycle        1.0.5      2026-01-08 [1] CRAN (R 4.5.2)
 listenv          0.10.1     2026-03-10 [1] CRAN (R 4.5.2)
 loo              2.9.0      2025-12-23 [1] CRAN (R 4.5.2)
 lubridate        1.9.5      2026-02-04 [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)
 ModelMetrics     1.2.2.2    2020-03-17 [3] CRAN (R 4.0.1)
 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)
 NatureSounds     1.0.5      2025-01-17 [1] CRAN (R 4.5.0)
 nlme             3.1-168    2025-03-31 [4] CRAN (R 4.4.3)
 nnet             7.3-20     2025-01-01 [4] CRAN (R 4.4.2)
 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)
 parallelly       1.47.0     2026-04-17 [1] CRAN (R 4.5.2)
 pbapply          1.7-4      2025-07-20 [1] CRAN (R 4.5.0)
 pillar           1.11.1     2025-09-17 [1] CRAN (R 4.5.0)
 pkgbuild         1.4.8      2025-05-26 [1] CRAN (R 4.5.0)
 pkgconfig        2.0.3      2019-09-22 [1] CRAN (R 4.5.0)
 pkgload          1.4.1      2025-09-23 [1] CRAN (R 4.5.0)
 plyr             1.8.9      2023-10-02 [1] CRAN (R 4.5.0)
 posterior        1.6.1      2025-02-27 [1] CRAN (R 4.5.0)
 pROC             1.18.0     2021-09-03 [3] CRAN (R 4.1.1)
 prodlim          2026.03.11 2026-03-11 [1] CRAN (R 4.5.2)
 profvis          0.4.0      2024-09-20 [1] CRAN (R 4.5.0)
 promises         1.3.3      2025-05-29 [1] RSPM (R 4.5.0)
 proxy            0.4-29     2025-12-29 [1] CRAN (R 4.5.2)
 purrr            1.2.0      2025-11-04 [1] CRAN (R 4.5.0)
 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)
 RCurl            1.98-1.17  2025-03-22 [1] CRAN (R 4.5.0)
 recipes          1.3.0      2025-04-17 [1] CRAN (R 4.5.0)
 remotes          2.5.0      2024-03-17 [1] CRAN (R 4.5.0)
 reshape2         1.4.5      2025-11-12 [1] CRAN (R 4.5.0)
 rjson            0.2.23     2024-09-16 [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)
 rpart            4.1.24     2025-01-07 [4] CRAN (R 4.4.2)
 Rraven         * 1.0.16     2017-11-17 [1] CRAN (R 4.5.0)
 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)
 Rtsne          * 0.17       2023-12-07 [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)
 seewave          2.2.4      2025-08-19 [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)
 signal           1.8-1      2024-06-26 [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)
 testthat         3.2.3      2025-01-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)
 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)
 timechange       0.4.0      2026-01-29 [1] CRAN (R 4.5.2)
 timeDate         4041.110   2024-09-22 [3] CRAN (R 4.5.0)
 tuneR            1.4.7      2024-04-17 [1] 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)
 warbleR          1.1.37     2025-10-22 [1] CRAN (R 4.5.0)
 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.

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