Code
# options to customize chunk outputs
knitr::opts_chunk$set(
tidy = TRUE, # soft-wrap code in output
tidy.opts = list(width.cutoff = 65) # width cutoff for soft-wrapping code
)Choosing the best distance metric for quantifying vocal convergence
Vocal convergence and social behavior in budgies
Purpose
Figure out the best way to quantify convergence
Quantify vocal convergence
Load packages
Code
# 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("PhenotypeSpace", "warbleR", "T4transport", "ggplot2",
"viridis", "pbapply", "energy", "beepr")
# install/ load packages
sketchy::load_packages(packages = packages)Loading required package: PhenotypeSpaceLoading required package: warbleRLoading required package: tuneRLoading required package: seewaveLoading required package: NatureSoundsLoading required package: knitrLoading required package: T4transport** ------------------------------------------------------- **** T4transport || Computational Optimal Transport in R **** Maintainer : Kisung You (kisung.you@outlook.com)** Version : 0.1.8 (2026)** Website : https://www.kisungyou.com/T4transport/**** Please share any bugs or suggestions to the maintainer.** ------------------------------------------------------- **Loading required package: ggplot2Loading required package: viridisLoading required package: viridisLiteLoading required package: pbapplyLoading required package: energyLoading required package: beepr
Attaching package: 'beepr'The following object is masked from 'package:seewave':
beep1 PCA dimensionality reduction
Code
est <- readRDS("./data/raw/extended_selection_table.RDS")
mfs <- mfcc_stats(X = est, parallel = 5, bp = c(0.5, 8))
saveRDS(mfs, "./data/processed/mfccs.RDS")Code
mfs <- readRDS("./data/processed/mfccs.RDS")
# run PCA on acoustic variables
pca <- prcomp(mfs[, c(-1, -2)], center = TRUE, scale. = TRUE)
saveRDS(pca, "./data/processed/pca.RDS")
# plot rotation values by PC
summ_pca <- summary(pca)
summ_pca$importance[2, ]
# select the PCs explaining at least 90%
pcs <- colnames(summ_pca$importance)[1:min(which(summ_pca$importance[3,
] > 0.9))]
pca_dat <- cbind(mfs[, c(1, 2)], as.data.frame(pca$x)[pcs])
pca_dat$record.block <- sapply(pca_dat$sound.files, function(x) unique(est$record.block[est$sound.files ==
x]))
pca_dat$record.group <- sapply(pca_dat$sound.files, function(x) unique(est$record.group[est$sound.files ==
x]))
pca_dat$date <- sapply(pca_dat$sound.files, function(x) unique(est$date[est$sound.files ==
x]))
pca_dat$treatment <- sapply(pca_dat$sound.files, function(x) unique(est$treatment[est$sound.files ==
x]))
pca_dat$round <- sapply(pca_dat$sound.files, function(x) unique(est$round[est$sound.files ==
x]))
pca_dat$ID <- sapply(pca_dat$sound.files, function(x) unique(est$ID[est$sound.files ==
x]))
saveRDS(pca_dat, "./data/processed/pca_dat.RDS")The first 69 PCs explain 90% of the variance in the data. The following graph shows the cumulative proportion of variance explained by each PC. The dashed red line indicates the 90% threshold.
Code
pca <- readRDS("./data/processed/pca.RDS")
# plot rotation values by PC
summ_pca <- summary(pca)
summ_pca_df <- as.data.frame(t(summ_pca$importance))
summ_pca_df$PC <- as.character(1:nrow(summ_pca_df))
summ_pca_df$PC <- factor(summ_pca_df$PC, levels = summ_pca_df$PC)
summ_pca_df <- summ_pca_df[summ_pca_df$PC %in% 1:80, ]
## ggplot barplot for cummulative proportion
ggplot(data = summ_pca_df, aes(x = PC, y = `Cumulative Proportion`)) +
geom_bar(stat = "identity", fill = "steelblue") + geom_hline(yintercept = 0.9,
linetype = "dashed", color = "red") + labs(x = "Principal Component",
y = "Cumulative Proportion of Variance Explained") + theme_classic()
2 Simulate low sample scenarios
Code
pca_dat <- readRDS("./data/processed/pca_dat.RDS")
tab <- table(pca_dat$ID[pca_dat$record.block != 1], pca_dat$record.block[pca_dat$record.block !=
1])
counts <- as.data.frame(tab)
names(counts) <- c("ID", "block", "count")
counts$group <- sapply(counts$ID, function(x) unique(pca_dat$record.group[pca_dat$ID ==
x]))
cutoff <- 200
subcounts <- counts[counts$count >= cutoff, ]
sum(c(table(subcounts$group, subcounts$block)) > 1)[1] 22Code
subcounts <- as.data.frame(table(subcounts$group, subcounts$block))
names(subcounts) <- c("group", "block", "Freq")
sel_groups <- subcounts[subcounts$Freq > 1, ]Code
source("~/Dropbox/R_package_testing/PhenotypeSpace/R/space_similarity.R")
source("~/Dropbox/R_package_testing/PhenotypeSpace/R/rarefact_space_similarity.R")
source("~/Dropbox/R_package_testing/PhenotypeSpace/R/internal_functions.R")
sample_cutoffs <- seq(5, 200, by = 1)
distance_methods <- c("Wasserstein", "Euclidean", "energy")
out <- lapply(seq_len(nrow(sel_groups)), function(x) {
subpca_dat <- pca_dat[pca_dat$record.group == sel_groups$group[x] &
pca_dat$record.block == sel_groups$block[x], ]
# only those above cutoff
id_count <- table(subpca_dat$ID)
subpca_dat <- subpca_dat[subpca_dat$ID %in% names(id_count)[id_count >
cutoff], ]
best_dists <- list()
for (dm in distance_methods) {
best_dists[[length(best_dists) + 1]] <- space_similarity(formula = as.formula(paste("ID ~",
paste(grep("PC", names(subpca_dat), value = TRUE), collapse = " + "))),
data = subpca_dat, method = "distance", distance.method = dm,
pb = FALSE)
}
names(best_dists) <- distance_methods
dists <- lapply(sample_cutoffs, function(y) {
sel_ID <- subpca_dat$ID[1]
X1 <- subpca_dat[subpca_dat$ID %in% sel_ID, ]
X2 <- subpca_dat[!subpca_dat$ID %in% sel_ID, ]
X1 <- X1[sample(seq_len(nrow(X1)), size = y, replace = FALSE),
]
W <- rbind(X1, X2)
out_dfs <- list()
for (dm in distance_methods) {
# default
dists <- space_similarity(formula = as.formula(paste("ID ~",
paste(grep("PC", names(subpca_dat), value = TRUE),
collapse = " + "))), data = W, method = "distance",
distance.method = dm, pb = FALSE)
dists$full.distance <- best_dists[[dm]]$distance
dists <- dists[dists$group.1 == sel_ID | dists$group.2 ==
sel_ID, ]
dists$method <- dm
# rarefacted
dists_rare <- rarefact_space_similarity(formula = as.formula(paste("ID ~",
paste(grep("PC", names(subpca_dat), value = TRUE),
collapse = " + "))), data = W, method = "distance",
distance.method = dm, pb = FALSE, n = 30)
dists_rare$full.distance <- best_dists[[dm]]$distance
dists_rare <- dists_rare[dists_rare$group.1 == sel_ID |
dists_rare$group.2 == sel_ID, ]
dists_rare$distance <- dists_rare$mean.distance
dists_rare$method <- paste0(dm, "_rare")
out_df <- rbind(dists, dists_rare[, names(dists)])
out_df$min.sample <- y
out_df$group <- sel_groups$group[x]
out_df$block <- sel_groups$block[x]
out_dfs[[length(out_dfs) + 1]] <- out_df
}
dists_df <- do.call(rbind, out_dfs)
return(dists_df)
})
dists_df <- do.call(rbind, dists)
return(dists_df)
})
results <- do.call(rbind, out)
saveRDS(results, "./data/processed/simulated_low_sample_scenarios.RDS")
beep(3)3 Graph results
Code
results <- readRDS("./data/processed/simulated_low_sample_scenarios.RDS")
results$prop.error <- abs(results$full.distance - results$distance)/results$full.distance
# ggplot with lines
agg <- aggregate(prop.error ~ min.sample + method, data = results,
FUN = mean)
agg$sd <- aggregate(prop.error ~ min.sample + method, data = results,
FUN = sd)$prop.error
# plot
ggplot(agg, aes(x = min.sample, y = prop.error, color = method)) +
geom_line() + geom_ribbon(aes(ymin = prop.error - sd, ymax = prop.error +
sd), alpha = 0.2) + labs(x = "Minimum sample size", y = "Mean proportional error") +
theme_minimal()
Code
# add column rarefaction
agg$rare <- grepl("rare", agg$method)
agg$method <- gsub("_rare", "", agg$method)
agg$method <- factor(agg$method, levels = c("Wasserstein", "energy",
"Euclidean"))
ggplot(agg, aes(x = min.sample, y = prop.error, color = rare)) + geom_line() +
geom_ribbon(aes(ymin = prop.error - sd, ymax = prop.error + sd),
alpha = 0.2) + labs(x = "Minimum sample size", y = "Mean proportional error") +
theme_minimal() + facet_wrap(~method)
Zoom in in the y axis:
Code
ggplot(agg, aes(x = min.sample, y = prop.error, color = rare)) + geom_line() +
geom_ribbon(aes(ymin = prop.error - sd, ymax = prop.error + sd),
alpha = 0.2) + labs(x = "Minimum sample size", y = "Mean proportional error") +
theme_minimal() + facet_wrap(~method) + ylim(c(0, 0.25))Warning: Removed 31 rows containing missing values or values outside the scale range
(`geom_line()`).Warning: Removed 159 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
4 Check correlation between full data set and down sampled data set distances
Similarity metrics can still be useful if the keep the relative difference between distances even if the absolute distance is not accurate. To check this we can calculate the correlation between the full data set distances and the down sampled distances.
The following graph shows the correlation between the full data set distances and the down sampled distances for each method (columns) and minimum sample size (rows). The dashed line represents the 1:1 line, where the inferred distance equals the true distance.
Code
# order levels
results$method <- factor(results$method, levels = c("Wasserstein",
"energy", "Euclidean", "Wasserstein_rare", "energy_rare", "Euclidean_rare"))
# do scatterplot distance vs full
ggplot(results[results$min.sample %in% c(5, 10, 20, 30), ], aes(y = full.distance,
x = distance, color = method)) + geom_point(alpha = 0.5) + scale_color_viridis_d(option = "G",
guide = "none", end = 0.8) + geom_abline(slope = 1, intercept = 0,
linetype = "dashed") + labs(y = "True distance", x = "Inferred distance") +
theme_minimal() + facet_wrap(~min.sample + method, scales = "free",
ncol = 6)
The following graph shows the pearson correlation (y axis) between the full data set distances and the subsampled data for each method (columns) and minimum sample size (x axis)
Code
# get pearson correlation for each min.sample
cors_list <- lapply(split(results, f = interaction(results$min.sample,
results$method)), function(x) {
cor(x$full.distance, x$distance, method = "pearson")
})
# put in data frame
cors_df <- as.data.frame(do.call(rbind, cors_list))
cors_df$min.sample <- sapply(strsplit(rownames(cors_df), "\\."), function(x) x[1])
cors_df$method <- sapply(strsplit(rownames(cors_df), "\\."), function(x) x[2])
cors_df$rare <- grepl("rare", cors_df$method)
cors_df$method <- gsub("_rare", "", cors_df$method)
cors_df$correlation <- as.numeric(cors_df$V1)
cors_df$min.sample <- as.numeric(cors_df$min.sample)
ggplot(cors_df, aes(x = min.sample, y = correlation, color = rare)) +
scale_color_viridis_d(option = "G", end = 0.8, alpha = 0.6) +
geom_line() + labs(x = "Minimum sample size", y = "Pearson correlation",
color = "Rarefaction") + theme_minimal() + facet_wrap(~method,
ncol = 6)
Takeaways
Rarefaction is not that great ???
A minimum sample size between 30-50 seems to be a good compromise between accuracy and sample size
Session information
─ Session info ───────────────────────────────────────────────────────────────
setting value
version R version 4.3.2 (2023-10-31)
os Ubuntu 22.04.2 LTS
system x86_64, linux-gnu
ui X11
language es_ES:en
collate en_US.UTF-8
ctype en_US.UTF-8
tz America/Costa_Rica
date 2026-07-09
pandoc 3.6.3 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/x86_64/ (via rmarkdown)
quarto 1.4.549 @ /usr/local/bin/quarto
─ Packages ───────────────────────────────────────────────────────────────────
package * version date (UTC) lib source
abind 1.4-8 2024-09-12 [1] CRAN (R 4.3.2)
audio 0.1-12 2025-12-15 [1] CRAN (R 4.3.2)
beepr * 2.0 2024-07-06 [1] CRAN (R 4.3.2)
bitops 1.0-9 2024-10-03 [1] CRAN (R 4.3.2)
boot 1.3-32 2025-08-29 [1] CRAN (R 4.3.2)
brio 1.1.5 2024-04-24 [1] CRAN (R 4.3.2)
cachem 1.1.0 2024-05-16 [1] CRAN (R 4.3.2)
class 7.3-23 2025-01-01 [1] CRAN (R 4.3.2)
classInt 0.4-11 2025-01-08 [1] CRAN (R 4.3.2)
cli 3.6.6 2026-04-09 [1] CRAN (R 4.3.2)
cluster 2.1.8.2 2026-02-05 [1] CRAN (R 4.3.2)
codetools 0.2-20 2024-03-31 [1] CRAN (R 4.3.2)
crayon 1.5.3 2024-06-20 [1] CRAN (R 4.3.2)
curl 7.1.0 2026-04-22 [1] CRAN (R 4.3.2)
DBI 1.3.0 2026-02-25 [1] CRAN (R 4.3.2)
deldir 2.0-4 2024-02-28 [1] CRAN (R 4.3.2)
devtools 2.4.6 2025-10-03 [1] CRAN (R 4.3.2)
digest 0.6.39 2025-11-19 [1] CRAN (R 4.3.2)
dplyr 1.2.1 2026-04-03 [1] CRAN (R 4.3.2)
dtw 1.23-3 2026-06-09 [1] CRAN (R 4.3.2)
e1071 1.7-17 2025-12-18 [1] CRAN (R 4.3.2)
ellipsis 0.3.3 2026-04-04 [1] CRAN (R 4.3.2)
energy * 1.7-12 2024-08-24 [1] CRAN (R 4.3.2)
evaluate 1.0.5 2025-08-27 [1] CRAN (R 4.3.2)
farver 2.1.2 2024-05-13 [1] CRAN (R 4.3.2)
fastmap 1.2.0 2024-05-15 [1] CRAN (R 4.3.2)
fftw 1.0-9 2024-09-20 [1] CRAN (R 4.3.2)
formatR 1.14 2023-01-17 [1] CRAN (R 4.3.2)
fs 2.0.1 2026-03-24 [1] CRAN (R 4.3.2)
generics 0.1.4 2025-05-09 [1] CRAN (R 4.3.2)
ggplot2 * 4.0.3 2026-04-22 [1] CRAN (R 4.3.2)
glue 1.8.1 2026-04-17 [1] CRAN (R 4.3.2)
goftest 1.2-3 2021-10-07 [3] CRAN (R 4.1.1)
gridExtra 2.3.1 2026-06-25 [1] CRAN (R 4.3.2)
gsl 2.1-8 2023-01-24 [1] CRAN (R 4.3.2)
gtable 0.3.6 2024-10-25 [1] CRAN (R 4.3.2)
htmltools 0.5.9 2025-12-04 [1] CRAN (R 4.3.2)
htmlwidgets 1.6.4 2023-12-06 [1] CRAN (R 4.3.2)
httr 1.4.8 2026-02-13 [1] CRAN (R 4.3.2)
jsonlite 2.0.0 2025-03-27 [1] CRAN (R 4.3.2)
KernSmooth 2.23-26 2025-01-01 [1] CRAN (R 4.3.2)
knitr * 1.51 2025-12-20 [1] CRAN (R 4.3.2)
labeling 0.4.3 2023-08-29 [1] CRAN (R 4.3.2)
lattice 0.22-9 2026-02-09 [1] CRAN (R 4.3.2)
lifecycle 1.0.5 2026-01-08 [1] CRAN (R 4.3.2)
magrittr 2.0.5 2026-04-04 [1] CRAN (R 4.3.2)
MASS 7.3-55 2022-01-13 [4] CRAN (R 4.1.2)
Matrix 1.6-5 2024-01-11 [1] CRAN (R 4.3.2)
memoise 2.0.1 2021-11-26 [3] CRAN (R 4.1.2)
mgcv 1.8-39 2022-02-24 [4] CRAN (R 4.1.2)
NatureSounds * 1.0.5 2025-01-17 [1] CRAN (R 4.3.2)
nlme 3.1-169 2026-03-27 [1] CRAN (R 4.3.2)
otel 0.2.0 2025-08-29 [1] CRAN (R 4.3.2)
packrat 0.9.3 2025-06-16 [1] CRAN (R 4.3.2)
pbapply * 1.7-4 2025-07-20 [1] CRAN (R 4.3.2)
permute 0.9-10 2026-02-06 [1] CRAN (R 4.3.2)
PhenotypeSpace * 0.1.0 2024-10-04 [1] CRAN (R 4.3.2)
pillar 1.11.1 2025-09-17 [1] CRAN (R 4.3.2)
pkgbuild 1.4.8 2025-05-26 [1] CRAN (R 4.3.2)
pkgconfig 2.0.3 2019-09-22 [3] CRAN (R 4.0.1)
pkgload 1.5.3 2026-06-15 [1] CRAN (R 4.3.2)
polyclip 1.10-7 2024-07-23 [1] CRAN (R 4.3.2)
proxy 0.4-29 2025-12-29 [1] CRAN (R 4.3.2)
purrr 1.2.2 2026-04-10 [1] CRAN (R 4.3.2)
R6 2.6.1 2025-02-15 [1] CRAN (R 4.3.2)
raster 3.6-32 2025-03-28 [1] CRAN (R 4.3.2)
rbibutils 2.4.1 2026-01-21 [1] CRAN (R 4.3.2)
RColorBrewer 1.1-3 2022-04-03 [1] CRAN (R 4.3.2)
Rcpp 1.1.1-1.1 2026-04-24 [1] CRAN (R 4.3.2)
RCurl 1.98-1.19 2026-06-03 [1] CRAN (R 4.3.2)
Rdpack 2.6.6 2026-02-08 [1] CRAN (R 4.3.2)
remotes 2.5.0 2024-03-17 [1] CRAN (R 4.3.2)
rjson 0.2.23 2024-09-16 [1] CRAN (R 4.3.2)
rlang 1.2.0 2026-04-06 [1] CRAN (R 4.3.2)
rmarkdown 2.31 2026-03-26 [1] CRAN (R 4.3.2)
rstudioapi 0.19.0 2026-06-11 [1] CRAN (R 4.3.2)
S7 0.2.2 2026-04-22 [1] CRAN (R 4.3.2)
scales 1.4.0 2025-04-24 [1] CRAN (R 4.3.2)
seewave * 2.2.4 2025-08-19 [1] CRAN (R 4.3.2)
sessioninfo 1.2.4 2026-06-04 [1] CRAN (R 4.3.2)
sf 1.0-22 2025-11-10 [1] CRAN (R 4.3.2)
signal 1.8-1 2024-06-26 [1] CRAN (R 4.3.2)
sketchy 1.0.6 2025-12-04 [1] local (/home/marce/Dropbox/R_package_testing/sketchy)
sp 2.2-1 2026-02-13 [1] CRAN (R 4.3.2)
spatstat.data 3.1-9 2025-10-18 [1] CRAN (R 4.3.2)
spatstat.explore 3.8-1 2026-05-24 [1] CRAN (R 4.3.2)
spatstat.geom 3.8-1 2026-05-23 [1] CRAN (R 4.3.2)
spatstat.random 3.5-0 2026-05-24 [1] CRAN (R 4.3.2)
spatstat.sparse 3.2-0 2026-05-21 [1] CRAN (R 4.3.2)
spatstat.univar 3.2-0 2026-05-18 [1] CRAN (R 4.3.2)
spatstat.utils 3.2-3 2026-05-10 [1] CRAN (R 4.3.2)
stringi 1.8.7 2025-03-27 [1] CRAN (R 4.3.2)
stringr 1.6.0 2025-11-04 [1] CRAN (R 4.3.2)
T4transport * 0.1.8 2026-01-11 [1] CRAN (R 4.3.2)
tensor 1.5.1 2025-06-17 [1] CRAN (R 4.3.2)
terra 1.8-80 2025-11-05 [1] CRAN (R 4.3.2)
testthat 3.3.2 2026-01-11 [1] CRAN (R 4.3.2)
tibble 3.3.1 2026-01-11 [1] CRAN (R 4.3.2)
tidyselect 1.2.1 2024-03-11 [1] CRAN (R 4.3.2)
tuneR * 1.4.7 2024-04-17 [1] CRAN (R 4.3.2)
units 1.0-1 2026-03-11 [1] CRAN (R 4.3.2)
usethis 3.2.1 2025-09-06 [1] CRAN (R 4.3.2)
vctrs 0.7.3 2026-04-11 [1] CRAN (R 4.3.2)
vegan 2.7-5 2026-05-25 [1] CRAN (R 4.3.2)
viridis * 0.6.5 2024-01-29 [1] CRAN (R 4.3.2)
viridisLite * 0.4.3 2026-02-04 [1] CRAN (R 4.3.2)
warbleR * 1.1.37 2025-10-22 [1] CRAN (R 4.3.2)
withr 3.0.3 2026-06-19 [1] CRAN (R 4.3.2)
xaringanExtra 0.8.0 2024-05-19 [1] CRAN (R 4.3.2)
xfun 0.59 2026-06-19 [1] CRAN (R 4.3.2)
yaml 2.3.12 2025-12-10 [1] CRAN (R 4.3.2)
[1] /home/marce/R/x86_64-pc-linux-gnu-library/4.3
[2] /usr/local/lib/R/site-library
[3] /usr/lib/R/site-library
[4] /usr/lib/R/library
* ── Packages attached to the search path.
──────────────────────────────────────────────────────────────────────────────