LBH song type prevalence and degradation
Source code and data found at https://github.com/maRce10/lbh_songtype_prevalence_and_degradation
# 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", "rprojroot", "warbleR", "baRulho", "readxl", "Rraven", "dplyr", "dtw", github = "maRce10/PhenotypeSpace", "umap", "mclust", "viridis", "caret"))contours <- read.csv(file = "./data/processed/seltailor_output.csv", header = T)
contours$tailored <- NULL
rav_dat_raw <- imp_raven(path = "./data/processed", files = "FrameworkContourDegradationv1.4.txt", all.data = TRUE, warbler.format = TRUE)
rav_dat_raw$segment.id <- paste(rav_dat_raw$indxF, rav_dat_raw$selec, sep = "-")# check all have 1 sample
all(table(rav_dat_raw$segment.id) == 1)
# number of fundamental freq values
ff_n <- length(grep("ff", names(contours)))
segment_list <- pblapply_wrblr_int(X = unique(rav_dat_raw$segment.id), cl = 10, function(x){
segment_data <- rav_dat_raw[rav_dat_raw$segment.id == x, ]
segment_data[,c("start", "end")]
contour <- contours[contours$selec == segment_data$indxF[1], ]
times <- seq(contour$start, contour$end, length.out = ff_n)
freq_time <- data.frame(freq = t(contour[, grep("ff", names(contours))]), times)
names(freq_time)[1] <- "freq"
seg_freq_time <- freq_time$freq[freq_time$times >= segment_data$start & freq_time$times <= segment_data$end]
return(seg_freq_time)
})
names(segment_list) <- unique(rav_dat_raw$segment.id)
saveRDS(segment_list, "./data/processed/raw_segment_ff_list.RDS")All have 143 values (which is the maximum number of ff values in a segment)
segment_list <- readRDS("./data/processed/raw_segment_ff_list.RDS")
# temporarily remove "280" id segments
segment_list <- segment_list[!grepl("^280", names(segment_list))]
max_length <- max(sapply(segment_list, length))
pad_segment_list <- lapply(segment_list, function(x) {
pad_x <- c(x, rep(mean(x), max_length - length(x)))
return(pad_x)
})
names(pad_segment_list) <- names(segment_list)
range(sapply(pad_segment_list, length))
pad_segment_df <- do.call(rbind, pad_segment_list)
# no NAs
which(apply(pad_segment_df, 1, anyNA))
dtw_dist_pad_segments <- dtwDist(mx = pad_segment_df)
saveRDS(dtw_dist_pad_segments, "./data/processed/dtw_padded_segment_ff_distances.RDS")segment_list <- readRDS("./data/processed/raw_segment_ff_list.RDS")
# temporarily remove "280" id segments
segment_list <- segment_list[!grepl("^280", names(segment_list))]
max_length <- max(sapply(segment_list, length))
intrp_segment_list <- lapply(segment_list, function(x) {
if (length(x) > 1)
intrp_x <- stats::approx(x = x, n = max_length)$y
if (length(x) == 1)
intrp_x <- rep(x, max_length)
return(intrp_x)
})
names(intrp_segment_list) <- names(segment_list)
range(sapply(intrp_segment_list, length))
intrp_segment_df <- do.call(rbind, intrp_segment_list)
# no NAs
which(apply(intrp_segment_df, 1, anyNA))
dtw_dist_intrp_segments <- dtwDist(mx = intrp_segment_df)
saveRDS(dtw_dist_intrp_segments, "./data/processed/dtw_interpolated_segment_ff_distances.RDS")segment_list <- readRDS("./data/processed/raw_segment_ff_list.RDS")
# temporarily remove "280" id segments
segment_list <- segment_list[!grepl("^280", names(segment_list))]
max_length <- max(sapply(segment_list, length))
combs <- combn(x = names(segment_list), m = 2)
slided_segment_list <- pblapply_wrblr_int(seq_len(ncol(combs)), cl = 20, function(x) {
xs <- segment_list[combs[, x]]
min_x <- which.min(shrt_length <- sapply(xs, length))
min_length <- shrt_length[min_x]
lg_x <- xs[[-min_x]]
shrt_x <- xs[[min_x]]
stps <- abs(diff(shrt_length))
if (stps <= 1)
stps <- 1 else
stps <- 1:stps
dtw_dists <- sapply(stps, function(e, cor.method = cm) {
warbleR::try_na(dtw(lg_x[e:(e + min_length - 1)],
shrt_x)$distance)
})
return(min(dtw_dists, na.rm = TRUE))
})
slided_segment_df <- cbind(t(combs), unlist(slided_segment_list))
rownames(slided_segment_df) <- 1:nrow(slided_segment_df)
colnames(slided_segment_df) <- c("segment.id1", "segment.id2", "dtw.dist")
head(slided_segment_df)
dim(slided_segment_df)
slided_segment_df <- as.data.frame(slided_segment_df)
slided_segment_df$dtw.dist <- as.numeric(slided_segment_df$dtw.dist)
saveRDS(slided_segment_df, "./data/processed/dtw_slided_segment_ff_tabular.RDS")
slided_segment_dist <- rectangular_to_triangular(slided_segment_df)
# no NAs
which(apply(slided_segment_dist, 1, anyNA))
saveRDS(slided_segment_dist, "./data/processed/dtw_slided_segment_ff_distances.RDS")dtw_dist_pad_segments <- readRDS("./data/processed/dtw_padded_segment_ff_distances.RDS")
dtw_dist_intrp_segments <- readRDS("./data/processed/dtw_interpolated_segment_ff_distances.RDS")
dtw_dist_slided_segments <- readRDS("./data/processed/dtw_slided_segment_ff_distances.RDS")
mds_pad_segments <- cmdscale(slided_segment_dist, k = 10)
# mds_intrp_segments <- cmdscale(dtw_dist_intrp_segments, k = 10)
mds_slided_segments <- cmdscale(dtw_dist_slided_segments, k = 10)
mds_segments <- cbind(mds_pad_segments, mds_slided_segments, mds_intrp_segments)
mds_segments_df <- data.frame(segment.id = rownames(mds_pad_segments), mds_segments)
# head(mds_segments_df)
mds_segments_df$segment.cat <- sapply(mds_segments_df$segment.id, function(x) rav_dat_raw$Category[rav_dat_raw$segment.id == x])
# table(mds_segments_df$segment.cat)
# temporarily remove 199-681
mds_segments_df <- mds_segments_df[mds_segments_df$segment.id != "199-681",]
mds_segments_df$X1 <- scale(mds_segments_df$X1) * 1000
mds_segments_df$X11 <- scale(mds_segments_df$X11) * 1000
category_similarity <- space_similarity(
formula = segment.cat ~ X1 + X11,
data = mds_segments_df,
cores = 10,
method = "mean.density.overlap"
)
# set method parameters
ctrl <- caret::trainControl(method = "repeatedcv", repeats = 5)
seg.id <- mds_segments_df$segment.id
mds_segments_df$segment.id <- NULL
# get similarities ("overlap")
category_similarity <- space_similarity(
formula = segment.cat ~ .,
data = mds_segments_df,
cores = 4,
method = "rf",
trControl = ctrl,
tuneLength = 4,
seed = 123
)set.seed(123)
umap_result <- umap(mds_segments_df[, grep("^X", names(mds_segments_df))], n_components = 20)
mds_segments_df$UMAP1 <- umap_result$layout[, 1]
mds_segments_df$UMAP2 <- umap_result$layout[, 2]
mod_umap <- Mclust(umap_result$layout)
mds_segments_df$umap.segment.cat <- as.character(mod_umap$classification)
table(mds_segments_df$umap.segment.cat)
write.csv(mds_segments_df, "./data/processed/dtw_tabular_dimensions_umap_vectors_and_unsupervised_categories_segments.csv", row.names = FALSE)mds_segments_df <- read.csv("./data/processed/dtw_tabular_dimensions_umap_vectors_and_unsupervised_categories_segments.csv")
conf_mat <- confusionMatrix(data = factor(mds_segments_df$segment.cat),
reference = factor(as.factor(mds_segments_df$umap.segment.cat), levels = c(unique(mds_segments_df$umap.segment.cat), unique(mds_segments_df$segment.cat))))Warning in levels(reference) != levels(data): comprimento do objeto maior não é
múltiplo do comprimento do objeto menorWarning in confusionMatrix.default(data = factor(mds_segments_df$segment.cat),
: Levels are not in the same order for reference and data. Refactoring data to
match.conf_df <- as.data.frame(conf_mat$table)
conf_df$total <- sapply(conf_df$Reference, function(x) sum(mds_segments_df$umap.segment.cat ==
x))
conf_df$proportion <- conf_df$Freq/conf_df$total
conf_df <- conf_df[conf_df$Reference %in% unique(mds_segments_df$umap.segment.cat),
]
plot(mds_segments_df$UMAP1, mds_segments_df$UMAP2, col = viridis(length(unique(mds_segments_df$umap.segment.cat)))[as.numeric(mds_segments_df$umap.segment.cat)], cex = 4, pch = as.numeric(as.factor(mds_segments_df$umap.segment.cat)) + 10, main= paste("UMAP: 2 dimensions;", length(unique(mds_segments_df$umap.segment.cat)), "groups"), xlab = "UMAP1", ylab = "UMAP2")
conf_df <- conf_df[!conf_df$Prediction %in% 1:9, ]
ggplot(conf_df, aes(x = Reference, y = Prediction, fill = proportion)) +
geom_tile() + theme_bw() + coord_equal() + scale_fill_distiller(palette = "Greens",
direction = 1) + geom_text(aes(label = round(proportion, 2)),
color = "black") + labs(x = "Manually labeled category", y = "Unsupervised category") +
theme_classic() + theme(axis.text.x = element_text(color = "black",
size = 11, angle = 30, vjust = 0.8, hjust = 0.8))
Session information
R version 4.3.1 (2023-06-16)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.2 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3; LAPACK version 3.9.0
locale:
[1] LC_CTYPE=pt_BR.UTF-8 LC_NUMERIC=C
[3] LC_TIME=es_CR.UTF-8 LC_COLLATE=pt_BR.UTF-8
[5] LC_MONETARY=es_CR.UTF-8 LC_MESSAGES=pt_BR.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] caret_6.0-94 lattice_0.22-5 ggplot2_3.4.4
[4] viridis_0.6.5 viridisLite_0.4.2 mclust_6.0.1
[7] umap_0.2.10.0 PhenotypeSpace_0.1.0 dtw_1.23-1
[10] proxy_0.4-27 dplyr_1.1.4 Rraven_1.0.14
[13] readxl_1.4.3 baRulho_2.1.0 ohun_1.0.2
[16] warbleR_1.1.30 NatureSounds_1.0.4 seewave_2.2.3
[19] tuneR_1.4.6 rprojroot_2.0.4 formatR_1.14
[22] knitr_1.45
loaded via a namespace (and not attached):
[1] RColorBrewer_1.1-3 rstudioapi_0.15.0 jsonlite_1.8.8
[4] magrittr_2.0.3 sketchy_1.0.2 spatstat.utils_3.0-4
[7] farver_2.1.1 rmarkdown_2.25 vctrs_0.6.5
[10] RCurl_1.98-1.14 askpass_1.2.0 terra_1.7-71
[13] htmltools_0.5.7 raster_3.6-26 signal_1.8-0
[16] cellranger_1.1.0 pROC_1.18.5 parallelly_1.36.0
[19] KernSmooth_2.23-22 htmlwidgets_1.6.4 plyr_1.8.9
[22] testthat_3.2.1 lubridate_1.9.3 igraph_1.6.0
[25] lifecycle_1.0.4 iterators_1.0.14 pkgconfig_2.0.3
[28] Matrix_1.6-5 R6_2.5.1 fastmap_1.1.1
[31] future_1.33.1 digest_0.6.34 colorspace_2.1-0
[34] tensor_1.5 RSpectra_0.16-1 vegan_2.6-4
[37] labeling_0.4.3 fansi_1.0.6 spatstat.sparse_3.0-3
[40] timechange_0.3.0 polyclip_1.10-6 abind_1.4-5
[43] mgcv_1.9-1 compiler_4.3.1 remotes_2.4.2.1
[46] withr_3.0.0 backports_1.4.1 DBI_1.2.1
[49] MASS_7.3-60.0.1 lava_1.7.3 openssl_2.1.1
[52] rjson_0.2.21 classInt_0.4-10 permute_0.9-7
[55] ModelMetrics_1.2.2.2 tools_4.3.1 units_0.8-5
[58] fftw_1.0-8 future.apply_1.11.1 nnet_7.3-19
[61] goftest_1.2-3 glue_1.7.0 nlme_3.1-164
[64] grid_4.3.1 sf_1.0-15 checkmate_2.3.1
[67] cluster_2.1.6 reshape2_1.4.4 generics_0.1.3
[70] recipes_1.0.9 gtable_0.3.4 spatstat.data_3.0-4
[73] class_7.3-22 spatstat.core_2.4-4 data.table_1.15.0
[76] sp_2.1-3 Deriv_4.1.3 utf8_1.2.4
[79] spatstat.geom_3.2-8 stringr_1.5.1 foreach_1.5.2
[82] pillar_1.9.0 splines_4.3.1 survival_3.5-7
[85] deldir_2.0-2 tidyselect_1.2.0 pbapply_1.7-2
[88] gridExtra_2.3 stats4_4.3.1 xfun_0.41
[91] hardhat_1.3.0 timeDate_4032.109 brio_1.1.4
[94] stringi_1.8.3 yaml_2.3.8 xaringanExtra_0.7.0
[97] evaluate_0.23 codetools_0.2-19 tibble_3.2.1
[100] cli_3.6.2 rpart_4.1.23 reticulate_1.35.0
[103] Sim.DiffProc_4.8 munsell_0.5.0 Rcpp_1.0.12
[106] globals_0.16.2 spatstat.random_3.2-2 png_0.1-8
[109] parallel_4.3.1 gower_1.0.1 bitops_1.0-7
[112] listenv_0.9.1 ipred_0.9-14 scales_1.3.0
[115] prodlim_2023.08.28 e1071_1.7-14 packrat_0.9.2
[118] purrr_1.0.2 crayon_1.5.2 rlang_1.1.3
[121] rgeos_0.6-4