Acoustic and statistical analysis

Phyllostomid call evolution

Author

Marcelo Araya-Salas

Published

October 21, 2023

Source code and data found at https://github.com/maRce10/Phyllostomid-call-evolution

1 Format acoustic data

Code

org_est <- readRDS("./data/processed/curated_extended_selection_table.RDS")

# feature_reference(org_est)

est <- resample_est(org_est, samp.rate = 420, parallel = 20)

# Rraven::exp_est(est[est$sound.files == "M00030.WAV-song_Arti.lit.GC.010.txt",], path = "./data/processed",file.name = "M00030.WAV-song_Arti.lit.GC.010", selection.table = TRUE)
# 
# Rraven::exp_est(est[est$sound.files == "R0029331__08-01-201300-42-12__M00017.WAV-song_Lich.obs.GC.006.txt",], path = "./data/processed",file.name = "R0029331__08-01-201300-42-12__M00017.WAV-song_Lich.obs.GC.006", selection.table = TRUE)
# 
fix_sfs <- imp_raven(path = "./data/processed", files = c("Lich.obs.GC.006_Fixed.txt", "Arti_lit.GC.010_Fixed.txt"), warbler.format = TRUE)
# 
# X<- fix_sfs
# 
# catalog(X = X, nrow = 9, ncol = 8, same.time.scale = T, mar = 0.001, res = 100, pb = FALSE, sub.legend = T, spec.mar = 1, ovlp = 90, wl = 100, width = 10 * 2.3, height = 5 * 2.5, hatching = 0, cex = 1.3, fast.spec = FALSE, pal = viridis, img.prefix = "delete", rm.axes = TRUE, path = "./data/processed", flim = c(min(X$bottom.freq), max(X$top.freq)) + c(-5, 5), highlight = TRUE, alpha = 0.25, collevels = seq(-120, 0, 5))


sub_est <- est[!est$sound.files %in% c("R0029331__08-01-201300-42-12__M00017.WAV-song_Lich.obs.GC.006.txt", "M00030.WAV-song_Arti.lit.GC.010.txt"),]

fix_sfs$selec.file[nrow(fix_sfs)] <- fix_sfs$selec.file[nrow(fix_sfs) - 1] 

fix_est <- selection_table(fix_sfs, extended = TRUE, confirm.extended = FALSE, by.song = "selec.file", path = "./data/processed", mar = 0.005)

setdiff(names(sub_est), names(fix_est))

fix_est$especie <- ifelse(fix_est$sound.files == "R0029331__08-01-201300-42-12__M00017.WAV-song_Lich.obs.GC.006.txt.wav-song_Lich.obs.GC.006.txt", "Lichonycteris obscura", "Artibeus lituratus")

fix_est$fuente <- "GChaverri"
fix_est$expansion <- NA

re_est <- rbind(sub_est, fix_est)


saveRDS(re_est, "./data/processed/curated_extended_selection_table_420_kHz.RDS")

Code

# leer el csv
selection_errors_list <- read.csv("/home/m/Downloads/SeleccionesAusentes.csv")

for (i in unique(X$archivo)
)
file.copy(from = file.path("smb://cinnas.local/neurobiología/marcelo_araya/recordings/Phyllostomid_recs/grabaciones/converted_sound_files", i), to = file.path("./data/raw/recordings", i))

2 check species

Code

est <- readRDS("./data/processed/curated_extended_selection_table_420_kHz.RDS")

phy <- ape::read.tree("./deprecated/Datos Rabosky/treePL_ML/chiroptera.no_outgroups.absolute.tre")

est$species <- gsub(" ", "_", est$especie)

dat_spp <- unique(est$species)

setdiff(dat_spp, phy$tip.label)


fix_names <- read_excel("./data/processed/Analysis_Species_Names_Errors_Database_Phylogeny.xlsx")

# #Check all Dermanura species available
# y <- "Microny"
# 
# # in tree
# unlist(lapply(phy$tip.label, function(x) grep(y, x, value = TRUE)))
# 
# # in  data
# unique(unlist(lapply(est$species, function(x) grep(y, x, value = TRUE))))
# 
# arts <- unlist(lapply(phy$tip.label, function(x) grep("Sturnira",x, value = TRUE)))


# fix name in phylogeny
phy$tip.label[phy$tip.label == "Micronycteris_nicefori"] <- "Trinycteris_nicefori"

# removing from data species missing in tree 
est <- est[!est$species %in% fix_names$Species_Data[is.na(fix_names$Species_Phylogenetic_Tree) ], ]

# keep species to rename in data
fix_names <- fix_names[!is.na(fix_names$Species_Phylogenetic_Tree), ]
fix_names <- fix_names[fix_names$wrong == "data set", ]

est$species[grep("nicefori", est$species)] <- "Trinycteris_nicefori"

# fix names in data set
for (i in 1:nrow(fix_names))
est$species[est$species == fix_names$Species_Data[i]] <- fix_names$Species_Phylogenetic_Tree[i]
  
length(unique(est$species))

est$species[est$species == "Micronycteris_brachyotis"] <- "Lampronycteris_brachyotis"
phy$tip.label[phy$tip.label == "Micronycteris_brachyotis"] <- "Lampronycteris_brachyotis"


saveRDS(est, "./data/processed/curated_extended_selection_table_72_species_sep_2023.RDS")



# trim phylo
setdiff(unique(est$species), phy$tip.label)

sub_phy <- drop.tip(phy, tip = setdiff(phy$tip.label, unique(est$species)))

write.tree(sub_phy, "./data/processed/trimmed_phylogeny_72_species_chiroptera.no_outgroups.absolute.tre")

3 Create spectrograms

Code

est <- readRDS("./data/processed/curated_extended_selection_table_72_species_sep_2023.RDS")


# loop to create catalogs
out <- pbapply::pblapply(unique(est$especie), function(i){
  
  X <- est[est$especie == i, ]

  X <- signal_2_noise(X, eq.dur = TRUE, bp = c(min(X$bottom.freq), max(X$top.freq)), pb = FALSE)

  X$xlabel <- paste(gsub(".txt$", "", X$selec.file), "snr:", round(X$SNR, 2))
 X$discrete.snr <- cut(X$SNR, breaks = 5)
  
  nrw <- ncl <- 2
  
  n <- nrow(X)
  if(n > 64) n <- 64
  
  if(n > 4) {
    nrw <- ncl <- ceiling(sqrt(n))
  if(((nrw-1) * ncl) >= n) nrw <- nrw - 1
  }

    
  X <- X[order(X$SNR), ]
  
  mako2 <- function(...) heat.colors(alpha = 0.5, ...)
  catalog(X = X, nrow = nrw, ncol = ncl, same.time.scale = T, mar = 0.001, res = 100, group.tag = "discrete.snr", pb = FALSE, sub.legend = T, spec.mar = 1, labels = c("xlabel", "selec"), max.group.cols = 5, title = i, ovlp = 90, wl = 100, width = 10 * 2.3, height = 5 * 2.5, tag.pal = list(viridis, mako, mako2), hatching = 0, cex = 1.3, fast.spec = FALSE, pal = viridis, img.prefix = i, rm.axes = TRUE, flim = c(min(X$bottom.freq), max(X$top.freq)) + c(-5, 5), highlight = TRUE, alpha = 0.25, collevels = seq(-120, 0, 5))  
  move_images(from = ".", to = "./output/catalogs", overwrite = TRUE, cut = TRUE, pb = FALSE)
  
 } 
)

# sp <- spectro_analysis(est[1:100, ], wl = 100, parallel = 20)

4 Check taxonomic coverage

Code

phy <- ape::read.tree("./deprecated/Datos Rabosky/treePL_ML/chiroptera.no_outgroups.absolute.tre")

phy$tip.label[phy$tip.label == "Micronycteris_nicefori"] <- "Trinycteris_nicefori"

sub_phy <- read.tree("./data/processed/trimmed_phylogeny_72_species_chiroptera.no_outgroups.absolute.tre")

all_phyllost <- read_excel("./data/processed/list_phyllostomidae_species_catalog_of_life.xlsx")[, 1:2]

all_phyllost$genus <- sapply(strsplit(all_phyllost$Species, " "), "[[", 1)

phyllost_in_tree <- grep(paste(unique(all_phyllost$genus), collapse = "|"), phy$tip.label, value = TRUE)

phyllost_mrca <- getMRCA(phy, tip = phyllost_in_tree)

phyllost_phy <- extract.clade(phy, node = phyllost_mrca)

tree2 <- groupOTU(phyllost_phy, .node = sub_phy$tip.label)

attr(tree2, "Condition")  <- ifelse(attributes(tree2)$group == 1, "In data", "Missing")

table(attr(tree2, "Condition"))

ggtree(tree2, aes(color=Condition), layout = "fan") + geom_tiplab(offset = 0) + scale_color_manual(values= c("black", "black")) + xlim(c(NA, 37))

# ggsave("file_")

ggtree(tree2, aes(color=Condition), layout = "fan") + geom_tiplab(offset = 0) + scale_color_manual(values= c(viridis(10)[7], "gray")) + xlim(c(NA, 37)) +
   theme(legend.position = c(.9, .9))

missing <- setdiff(phyllost_phy$tip.label, sub_phy$tip.label)

write.csv(data.frame(missing = missing), "./data/processed/missing_species.csv")

# search missing ones
l <- pbapply::pblapply(gsub("_", " ", missing), function(x) query_xenocanto(x, pb = FALSE, all_data = FALSE))

all(sapply(l, is.null))

5 Pantheria

Code

panth <- read.csv("./data/processed/PanTHERIA_1-0_WR05_Aug2008.csv", na.strings = "-999")

panth$species <- gsub(" ", "_", panth$MSW05_Binomial)

sub_panth <- panth[panth$species %in% (sub_phy$tip.label), ]
nrow(sub_panth)

naz_dat$species <- gsub(" ", "_", naz_dat$org.names)
naz_dat$alt_species <- gsub(" ", "_", naz_dat$alternative.name)

naz_dat <- read_excel("./deprecated/check_new_files/base de datos.xlsx")
sub_dat <- naz_dat[naz_dat$species %in% (sub_phy$tip.label) | naz_dat$alt_species %in% (sub_phy$tip.label), ]

write.csv(data.frame(sp = sub_phy$tip.label), "./data/processed/species_list_72_sp_sep_2023.csv", row.names = FALSE)

6 measure acoustic features

Code

est <- readRDS("./data/processed/curated_extended_selection_table_72_species_sep_2023.RDS")



acoust_data <- spectro_analysis(est, parallel = 2, ovlp = 70, wl = 100)

acoust_data$species <- sapply(seq_len(nrow(acoust_data)), function(x)
  unique(est$species[est$sound.files == acoust_data$sound.files[x]])
  )
acoust_data$expansion <- sapply(seq_len(nrow(acoust_data)), function(x)
  unique(est$expansion[est$sound.files == acoust_data$sound.files[x]])
  )

acoust_data$expansion[is.na(acoust_data$expansion)] <- 1
acoust_data$expansion <- as.factor(acoust_data$expansion)


write.csv(acoust_data, "./data/processed/acoustic_features_all_calls.csv", row.names = FALSE)

Code

acoust_data <- read.csv("./data/processed/acoustic_features_all_calls.csv")

# remove outliers
vars <-  c("duration", "meanfreq", "sd" , "freq.median", "freq.Q25", "freq.Q75", "freq.IQR", "time.median", "time.Q25", "time.Q75", "time.IQR", "skew", "kurt", "sfm", "meandom", "mindom",  "maxdom",  "dfrange", "modindx", "startdom", "enddom", "dfslope", "meanpeakf")
epsilon <- 1e-10

out <- warbleR:::pblapply_wrblr_int(unique(acoust_data$species), function(x){
  X <- acoust_data[acoust_data$species == x, ]
  
  Y <- scale(X[, vars])
  
  X <- X[complete.cases(Y), ]
  Y <- Y[complete.cases(Y), ]
  
  cov_mat <- cov(Y) + epsilon * diag(ncol(Y))

  # Compute the Mahalanobis distance
  mahalanobis_dist <- mahalanobis(x = Y, cov = cov_mat, center = TRUE)

  # Identify outliers
  outliers <- which(mahalanobis_dist > quantile(mahalanobis_dist, 0.90, na.rm = TRUE))
  
  X <- X[-outliers, ]

})

acoust_data <- do.call(rbind, out)

nrow(acoust_data)

[1] 8707

Code

acoust_data$file_number <- 1

for (i in 2:nrow(acoust_data)){
  if (acoust_data$sound.files[i - 1] == acoust_data$sound.files[i] & acoust_data$species[i - 1] == acoust_data$species[i]) acoust_data$file_number[i] <- acoust_data$file_number[i - 1]
  
   if (acoust_data$sound.files[i - 1] != acoust_data$sound.files[i] & acoust_data$species[i - 1] == acoust_data$species[i]) acoust_data$file_number[i] <- acoust_data$file_number[i - 1] + 1
  
     if (acoust_data$species[i - 1] != acoust_data$species[i]) acoust_data$file_number[i] <- 1
  
}

 theme_set(theme_classic(base_size = 20))


ggplot(data = acoust_data, aes(x = species, y = meanpeakf, shape = as.factor(expansion), color = as.factor(file_number))) +
  # geom_point(size = 4) +
  geom_jitter(size = 4) +
  geom_violin() +
  scale_color_viridis_d(alpha = 0.5, begin = 0.2, end = 0.8) +
     theme(axis.text.x = element_text(size = 14, angle=45, hjust = 1)) +
  geom_hline(yintercept = 20, lty = 4, col = "gray")

Code

ggplot(data = acoust_data, aes(x = species, y = duration, shape = as.factor(expansion), color = as.factor(file_number))) +
  # geom_point(size = 4) +
  geom_jitter(size = 4, width = 0.3) +
  geom_violin() +
  scale_color_viridis_d(alpha = 0.5, begin = 0.2, end = 0.8) +
   theme_classic() +
     theme(axis.text.x = element_text(size = 14, angle=45, hjust = 1))

Code

# raincloud plot:
fill_color <- adjustcolor(viridis(10)[6], 0.6)

ggplot(acoust_data, aes(x = as.factor(expansion), y = duration)) +
  # add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    fill = fill_color,
    alpha = 0.5,
    # custom bandwidth
    adjust = .5,
    # adjust height
    width = .6,
    .width = 0,
    # move geom to the cright
    justification = -.2,
    point_colour = NA
  ) +
  geom_boxplot(fill = fill_color,
    width = .15,
    # remove outliers
    outlier.shape = NA # `outlier.shape = NA` works as well
  ) +
  # add justified jitter from the {gghalves} package
  gghalves::geom_half_point(
    color = fill_color,
    # draw jitter on the left
    side = "l",
    # control range of jitter
    range_scale = .4,
    # add some transparency
    alpha = .5,
  )

Code

ggplot(acoust_data, aes(x = as.factor(expansion), y = meanpeakf)) +
  # add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    fill = fill_color,
    alpha = 0.5,
    # custom bandwidth
    adjust = .5,
    # adjust height
    width = .6,
    .width = 0,
    # move geom to the cright
    justification = -.2,
    point_colour = NA
  ) +
  geom_boxplot(fill = fill_color,
    width = .15,
    # remove outliers
    outlier.shape = NA # `outlier.shape = NA` works as well
  ) +
  # add justified jitter from the {gghalves} package
  gghalves::geom_half_point(
    color = fill_color,
    # draw jitter on the left
    side = "l",
    # control range of jitter
    range_scale = .4,
    # add some transparency
    alpha = .5,
  )

Code

ggplot(acoust_data, aes(x = as.factor(expansion), y = meandom)) +
  # add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    fill = fill_color,
    alpha = 0.5,
    # custom bandwidth
    adjust = .5,
    # adjust height
    width = .6,
    .width = 0,
    # move geom to the cright
    justification = -.2,
    point_colour = NA
  ) +
  geom_boxplot(fill = fill_color,
    width = .15,
    # remove outliers
    outlier.shape = NA # `outlier.shape = NA` works as well
  ) +
  # add justified jitter from the {gghalves} package
  gghalves::geom_half_point(
    color = fill_color,
    # draw jitter on the left
    side = "l",
    # control range of jitter
    range_scale = .4,
    # add some transparency
    alpha = .5,
  )

Code

ggplot(acoust_data, aes(x = as.factor(expansion), y = sp.ent)) +
  # add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    fill = fill_color,
    alpha = 0.5,
    # custom bandwidth
    adjust = .5,
    # adjust height
    width = .6,
    .width = 0,
    # move geom to the cright
    justification = -.2,
    point_colour = NA
  ) +
  geom_boxplot(fill = fill_color,
    width = .15,
    # remove outliers
    outlier.shape = NA # `outlier.shape = NA` works as well
  ) +
  # add justified jitter from the {gghalves} package
  gghalves::geom_half_point(
    color = fill_color,
    # draw jitter on the left
    side = "l",
    # control range of jitter
    range_scale = .4,
    # add some transparency
    alpha = .5,
  )

7 Compare GBIF and Pantheria data

7.1 Diet

Code

comp_data <- read.csv("./data/processed/Complete_Data_Sets/Pantheria_GBIF_Available_Data_Comparison.csv")

diet_comp_data <- comp_data[!is.na(comp_data$Diet_GBIF) & comp_data$Diet_Nazareth != "desconocido", ]


diet_comp_data$Diet_GBIF[diet_comp_data$Diet_GBIF == "posible frugivoro"] <- "frugivoro"
diet_comp_data$Diet_Nazareth[diet_comp_data$Diet_Nazareth == "nectarivoro/frugivoro"] <- "nectarivoro"
diet_comp_data$Diet_Nazareth[diet_comp_data$Diet_Nazareth == "insectivoro"] <- "carnivoro"

unique(diet_comp_data$Diet_Nazareth)

[1] "nectarivoro" "frugivoro"   "omnivoro"    "carnivoro"   "hematofago"

Code

unique(diet_comp_data$Diet_GBIF)

[1] "nectarivoro" "frugivoro"   "omnivoro"    "carnivoro"   "hematofago"

Code

diet_comp_data$Diet_Nazareth <- factor(diet_comp_data$Diet_Nazareth, levels = unique(c(diet_comp_data$Diet_Nazareth)))
diet_comp_data$Diet_GBIF <- factor(diet_comp_data$Diet_GBIF, levels = unique(c(diet_comp_data$Diet_Nazareth)))

conf_mat <- caret::confusionMatrix(diet_comp_data$Diet_GBIF, diet_comp_data$Diet_Nazareth)


conf_mat

Confusion Matrix and Statistics

             Reference
Prediction    nectarivoro frugivoro omnivoro carnivoro hematofago
  nectarivoro          13         0        2         1          0
  frugivoro             1        24        0         0          0
  omnivoro              0         1        2         2          0
  carnivoro             0         1        1        17          0
  hematofago            0         0        0         0          3

Overall Statistics
                                          
               Accuracy : 0.8676          
                 95% CI : (0.7636, 0.9377)
    No Information Rate : 0.3824          
    P-Value [Acc > NIR] : < 2.2e-16       
                                          
                  Kappa : 0.8165          
                                          
 Mcnemar's Test P-Value : NA              

Statistics by Class:

                     Class: nectarivoro Class: frugivoro Class: omnivoro
Sensitivity                      0.9286           0.9231         0.40000
Specificity                      0.9444           0.9762         0.95238
Pos Pred Value                   0.8125           0.9600         0.40000
Neg Pred Value                   0.9808           0.9535         0.95238
Prevalence                       0.2059           0.3824         0.07353
Detection Rate                   0.1912           0.3529         0.02941
Detection Prevalence             0.2353           0.3676         0.07353
Balanced Accuracy                0.9365           0.9496         0.67619
                     Class: carnivoro Class: hematofago
Sensitivity                    0.8500           1.00000
Specificity                    0.9583           1.00000
Pos Pred Value                 0.8947           1.00000
Neg Pred Value                 0.9388           1.00000
Prevalence                     0.2941           0.04412
Detection Rate                 0.2500           0.04412
Detection Prevalence           0.2794           0.04412
Balanced Accuracy              0.9042           1.00000

Code

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

conf_df$total <-
    sapply(conf_df$Reference, function(x)
        sum(diet_comp_data$Diet_Nazareth ==
                x))

conf_df$proportion <- conf_df$Freq / conf_df$total

ggplot(conf_df, aes(x = Reference, y = Prediction, fill = proportion)) +
    geom_tile() + theme_bw() + coord_equal() + labs(x = "External source", y = "GBIF") +
  scale_fill_distiller(palette = "Greens", direction = 1) + 
  geom_text(aes(label = round(proportion, 2)), color = "black", size = 3) + 
  theme_classic() + 
  theme(axis.text.x = element_text(color = "black", size = 11, angle = 30, vjust = 0.8, hjust = 0.8))

7.2 Adult body length

Code

cor.test(comp_data$AdultHeadBodyLen_mm_Pantheria, comp_data$AdultHeadBodyLen_mm_GBIF)


    Pearson's product-moment correlation

data:  comp_data$AdultHeadBodyLen_mm_Pantheria and comp_data$AdultHeadBodyLen_mm_GBIF
t = 4.8073, df = 12, p-value = 0.0004282
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.4929067 0.9380789
sample estimates:
      cor 
0.8113079

Code

ggplot(comp_data, aes(x = AdultHeadBodyLen_mm_Pantheria, y = AdultHeadBodyLen_mm_GBIF)) +
  geom_point() +
  geom_abline(slope = 1, intercept = 0, lty = 2, color = "gray")

7.3 Adult body mass

Code

cor.test(comp_data$AdultBodyMass_g_Pantheria, comp_data$AdultBodyMass_g_GBIF)


    Pearson's product-moment correlation

data:  comp_data$AdultBodyMass_g_Pantheria and comp_data$AdultBodyMass_g_GBIF
t = 19.973, df = 57, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.8932990 0.9612304
sample estimates:
      cor 
0.9354008

Code

ggplot(comp_data, aes(x = AdultBodyMass_g_Pantheria, y = AdultBodyMass_g_GBIF)) +
  geom_point() +
  geom_abline(slope = 1, intercept = 0, lty = 2, color = "gray")

7.4 Adult forearm length

Code

cor.test(comp_data$AdultForearmLen_mm_Pantheria, comp_data$AdultBodyMass_g_GBIF)


    Pearson's product-moment correlation

data:  comp_data$AdultForearmLen_mm_Pantheria and comp_data$AdultBodyMass_g_GBIF
t = 7.5656, df = 58, p-value = 3.283e-10
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.5491111 0.8131695
sample estimates:
     cor 
0.704768

Code

ggplot(comp_data, aes(x = AdultForearmLen_mm_Pantheria, y = AdultForearmLen_mm_GBIF)) +
  geom_point() +
  geom_abline(slope = 1, intercept = 0, lty = 2, color = "gray")

8 PCA on acoustic data

Code

# Taking only numeric data for the pca

acoust_data <- read.csv("./data/processed/acoustic_features_all_calls.csv")
acoust_data$species <- gsub(" ", "_", acoust_data$species)
acoust_data$freq_by_time <- acoust_data$freq.IQR / acoust_data$time.IQR

acoust_param <- c("duration", "meanfreq", "sd" , "freq.median", "freq.Q25", "freq.Q75", "freq.IQR", "time.median", "time.Q25", "time.Q75", "time.IQR", "skew", "kurt", "sp.ent",  "time.ent", "entropy", "sfm", "meandom", "mindom",  "maxdom",  "dfrange", "modindx", "startdom", "enddom", "dfslope", "meanpeakf", "freq_by_time")

pca_degrad <- prcomp(acoust_data[, acoust_param], scale = TRUE)

pca_var <- round(summary(pca_degrad)$importance[2, ] * 100)

acoust_data$pc1 <- pca_degrad$x[, 1]

# keep only columns needed

acoust_data <- acoust_data[, c("species", "sound.files", "pc1", acoust_param)]

# plot rotation values by PC
pca_rot <- as.data.frame(pca_degrad$rotation[, 1:4])

pca_rot_stck <- stack(pca_rot)

pca_rot_stck$variable <- rownames(pca_rot)
pca_rot_stck$values[pca_rot_stck$ind == "PC1"] <- pca_rot_stck$values[pca_rot_stck$ind ==
    "PC1"]
pca_rot_stck$Sign <- ifelse(pca_rot_stck$values > 0, "Positive", "Negative")
pca_rot_stck$rotation <- abs(pca_rot_stck$values)
pca_rot_stck$ind_var <- paste0(pca_rot_stck$ind, " (", sapply(pca_rot_stck$ind,
    function(x) pca_var[names(pca_var) == x]), "%)")

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

9 Estimate functional space

Code

trait_data <- as.data.frame(read_excel(path = "./data/processed/Complete_Data_Sets/Phyllostomidae_Analysis_Species_Temp_Base_13_10_23.xlsx", na = "NA"))

trait_data$species <- gsub(" ", "_", trait_data$MSW05_Binomial)


names(trait_data)

 [1] "MSW05_Genus"                         "MSW05_Species"                      
 [3] "MSW05_Binomial"                      "Phyllostomid_mean_mass_g"           
 [5] "Phyllostomid_mean_mass_male_g"       "Phyllostomid_mean_mass_female_g"    
 [7] "BodyMass_g_Phyllostomid_Book"        "AdultBodyMass_g"                    
 [9] "AdultBodyMass_g_EXT"                 "mean_forearm_female_mm"             
[11] "mean_forearm_male_mm"                "mean_forearm_mm"                    
[13] "AdultForearmLen_mm"                  "AdultHeadBodyLen_mm"                
[15] "AdultHeadBodyLen_mm_GBIF"            "average_ear_size_GBIF"              
[17] "Aggregation_size"                    "Roost_permanence"                   
[19] "BasalMetRate_mLO2hr"                 "BasalMetRateMass_g"                 
[21] "DietBreadth"                         "DispersalAge_d"                     
[23] "GestationLen_d"                      "HabitatBreadth"                     
[25] "HomeRange_km2"                       "HomeRange_Indiv_km2"                
[27] "InterbirthInterval_d"                "LitterSize"                         
[29] "LittersPerYear"                      "MaxLongevity_m"                     
[31] "NeonateBodyMass_g"                   "NeonateHeadBodyLen_mm"              
[33] "PopulationDensity_n.km2"             "PopulationGrpSize"                  
[35] "SexualMaturityAge_d"                 "SocialGrpSize"                      
[37] "TeatNumber"                          "Terrestriality"                     
[39] "TrophicLevel"                        "WeaningAge_d"                       
[41] "WeaningBodyMass_g"                   "WeaningHeadBodyLen_mm"              
[43] "References"                          "AdultBodyMass_g_EXT2"               
[45] "LittersPerYear_EXT"                  "NeonateBodyMass_g_EXT"              
[47] "WeaningBodyMass_g_EXT"               "GR_Area_km2"                        
[49] "GR_MaxLat_dd"                        "GR_MinLat_dd"                       
[51] "GR_MidRangeLat_dd"                   "GR_MaxLong_dd"                      
[53] "GR_MinLong_dd"                       "GR_MidRangeLong_dd"                 
[55] "HuPopDen_Min_n.km2"                  "HuPopDen_Mean_n.km2"                
[57] "HuPopDen_5p_n.km2"                   "HuPopDen_Change"                    
[59] "Precip_Mean_mm"                      "Temp_Mean_01degC"                   
[61] "AET_Mean_mm"                         "PET_Mean_mm"                        
[63] "Species_Phyllos_1_Name"              "Diet"                               
[65] "vertical_stratification"             "foraging_habit"                     
[67] "dietary_specialization_Farneda_2015" "source"                             
[69] "Roost permanence"                    "Foliage"                            
[71] "Tent"                                "Hollow tree"                        
[73] "Cave"                                "Aggregation size"                   
[75] "agg1-10"                             "agg11-25"                           
[77] "agg25-100"                           "agg100"                             
[79] "species"

Code

# remove variables with many NAs
# na_count <- sapply(trait_data, function(x) sum(is.na(x)))
# 
# na_count <- na_count[na_count < 20]
# 
# names(na_count)
# 
# trait_data <- trait_data[, names(na_count)]

funct_traits <-  c("BodyMass_g_Phyllostomid_Book", "mean_forearm_mm", "AdultHeadBodyLen_mm", "average_ear_size_GBIF", "HabitatBreadth", "Precip_Mean_mm", "Temp_Mean_01degC", "AET_Mean_mm", "Diet", "foraging_habit", "Foliage", "Tent", "Hollow tree", "Cave", "agg1-10", "agg11-25", "agg25-100", "agg100")   


num_traits <- names(trait_data)[sapply(trait_data, is.numeric)]
num_traits <- intersect(num_traits, funct_traits)

res.comp <- imputePCA(scale(trait_data[, num_traits]), ncp = 2)

imp_num_trait_data <- res.comp$completeObs

res.comp <- imputePCA(trait_data[, num_traits], ncp = 2)

imp_num_trait_data <- as.data.frame(res.comp$completeObs)

# Encode categorical data using one-hot encoding
dummy_diet <- model.matrix(~ Diet - 1, trait_data)
dummy_foraging_habit <- model.matrix(~ foraging_habit - 1, trait_data)


imp_num_trait_data$species <- trait_data$species
imp_num_trait_data <- imp_num_trait_data[, c("species", setdiff(names(imp_num_trait_data), "species"))]

imp_num_trait_data <- cbind(imp_num_trait_data, dummy_diet, dummy_foraging_habit)


head(imp_num_trait_data)

species	BodyMass_g_Phyllostomid_Book	mean_forearm_mm	AdultHeadBodyLen_mm	average_ear_size_GBIF	HabitatBreadth	Precip_Mean_mm	Temp_Mean_01degC	AET_Mean_mm	Foliage	Tent	Hollow tree	Cave	agg1-10	agg11-25	agg25-100	agg100	Dietfrugivoro	Dietnectarivoro	foraging_habitHighly cluttered space/gleaning nectarivore	foraging_habithighly cluttered/gleaning frugivore	foraging_habithighly cluttered/gleaning nectarivore	foraging_habithighly cluttered/gleaning nectivore
Ametrida_centurio	10.15	28.76	43.50	13.66667	1	159.0600	250.2900	1385.280	0.2595458	0.1855850	0.5686982	0.5910958	0.9161587	0.2396946	-0.0537009	-0.1064500	1	0	0	1	0	0
Anoura_caudifer	10.55	36.97	58.32	12.50000	1	150.8800	230.8500	1319.220	0.0000000	0.0000000	1.0000000	1.0000000	1.0000000	0.0000000	0.0000000	0.0000000	0	1	1	0	0	0
Anoura_cultrata	17.42	41.83	66.50	14.50000	1	164.2700	232.1100	1355.350	0.0000000	0.0000000	0.0000000	1.0000000	0.0000000	1.0000000	1.0000000	0.0000000	0	1	0	0	0	1
Anoura_geoffroyi	15.00	42.62	65.50	15.00000	1	128.9800	221.1400	1191.280	0.0000000	0.0000000	0.0000000	1.0000000	0.0000000	1.0000000	1.0000000	1.0000000	0	1	0	0	1	0
Ardops_nichollsi	17.63	46.56	63.00	14.75000	1	162.2361	242.3149	1375.200	1.0000000	0.0000000	0.0000000	0.0000000	1.0000000	0.0000000	0.0000000	0.0000000	1	0	0	1	0	0
Ariteus_flavescens	12.30	40.24	60.00	15.00000	1	104.7500	240.0000	1245.288	0.1174753	0.1316847	0.3554850	0.9504472	0.4908346	0.6828367	0.4702686	0.3611293	1	0	0	1	0	0

Code

# comps <- 10:40
# for(i in comps){
# Sys.sleep(1)
  
set.seed(123)
umap_result <- umap(imp_num_trait_data[, -1], n_components = 2)
pca_result <- prcomp(imp_num_trait_data[, -1], scale. = TRUE)


trait_data$UMAP1 <- umap_result$layout[, 1]
trait_data$UMAP2 <- umap_result$layout[, 2]
trait_data$PC1 <- pca_result$x[, 1]
trait_data$PC2 <- pca_result$x[, 2]


# plot(trait_data$UMAP1, trait_data$UMAP2)


phy <- read.tree("./data/processed/trimmed_phylogeny_72_species_chiroptera.no_outgroups.absolute.tre")

rownames(trait_data) <- trait_data$species

## set colors for mapped discrete character
# cols<-setNames(viridis(6),
    # unique(trait_data$Diet))

# diet <- trait_data$Diet
# names(diet) <- trait_data$species
# Perform stochastic character mapping
# simmap_phy <- make.simmap(phy, x = diet)
# plot(simmap_phy)

# phylomorphospace(tree = simmap_phy, trait_data[, c("UMAP1", "UMAP2")], colors = cols, bty="l",ftype="off",node.by.map=TRUE, node.size=c(0,1.2))
# 
cormat <- cor(trait_data[, c("BodyMass_g_Phyllostomid_Book", "mean_forearm_mm", "AdultHeadBodyLen_mm", "average_ear_size_GBIF", "HabitatBreadth", "Precip_Mean_mm", "Temp_Mean_01degC", "Foliage", "Tent", "Hollow tree", "Cave", "agg1-10", "agg11-25", "agg25-100", "agg100", "UMAP1", "UMAP2", "PC1", "PC2")], use = "pairwise.complete.obs")

rownames(cormat) <- colnames(cormat) <- c("BodyMass_g_Phyllostomid_Book", "mean_forearm_mm", "AdultHeadBodyLen_mm", "average_ear_size_GBIF", "HabitatBreadth", "Precip_Mean_mm", "Temp_Mean_01degC",  "Foliage", "Tent", "Hollow tree", "Cave", "agg1-10", "agg11-25", "agg25-100", "agg100", "UMAP1", "UMAP2", "PC1", "PC2")

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

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

Code

mod_umap <- Mclust(umap_result$layout)

guild_umap <- mod_umap$classification
names(guild_umap) <- trait_data$species


mod_pca <- Mclust(pca_result$x)

guild_pca <- mod_pca$classification
names(guild_pca) <- trait_data$species

# guild_phy <- make.simmap(phy, x = guild)
# plot(guild_phy)

# cols<-setNames(viridis(length(unique(guild))), unique(guild))

trait_data$guild_pca <- as.factor(sapply(trait_data$species, function(x) guild_pca[names(guild_pca) == x]))

trait_data$guild_umap <- as.factor(sapply(trait_data$species, function(x) guild_umap[names(guild_umap) == x]))

plot(trait_data$UMAP1, trait_data$UMAP2, col = viridis(length(unique(trait_data$guild_umap)))[trait_data$guild_umap],  cex = 4, pch = as.numeric(as.factor(trait_data$Diet)) + 10, main= paste("UMAP: 2 dimensions;", length(unique(trait_data$guild_umap)), "groups"))

Code

plot(trait_data$PC1, trait_data$PC2, col = viridis(length(unique(trait_data$guild_pca)))[trait_data$guild_pca],  cex = 4, pch = as.numeric(as.factor(trait_data$Diet)) + 10, main= paste("PCA;", length(unique(trait_data$guild_pca)), "groups"))

Code

# text(trait_data$UMAP1, trait_data$UMAP2, trait_data$species)
# 
# phylomorphospace(tree = guild_phy, trait_data[, c("UMAP1", "UMAP2")], colors = cols, bty="l",ftype="off",node.by.map=TRUE, node.size=c(0,1.2))
# 
# 
# summary(mod)
# }

write.csv(trait_data, "./data/processed/trait_and_guild_data.csv", row.names = FALSE)

10 Pool acoustic and ecological data together

Code

trait_data <- trait_data[, c("species", "BodyMass_g_Phyllostomid_Book", "mean_forearm_mm", "AdultHeadBodyLen_mm", "Diet", "Precip_Mean_mm", "Foliage", "Tent", "Hollow tree", "Cave", "agg1-10", "agg11-25", "agg25-100", "agg100", "UMAP1", "UMAP2", "PC1", "PC2", "guild_umap", "guild_pca")]


trait_data$log_mass <- log(trait_data$BodyMass_g_Phyllostomid_Book)
trait_data$log_mass_centered <- trait_data$log_mass - mean(trait_data$log_mass, na.rm = T)
trait_data$weight.forearm.ratio <- trait_data$BodyMass_g_Phyllostomid_Book/trait_data$mean_forearm_mm

phy <- read.tree("./data/processed/trimmed_phylogeny_72_species_chiroptera.no_outgroups.absolute.tre")

# 
# Ntip(phy)
# length(unique(trait_data$species))
# length(unique(acoust_data$species))
# intersect(trait_data$species, acoust_data$species)
# intersect(trait_data$species, phy$tip.label)
# intersect(acoust_data$species, phy$tip.label)
# 
# setdiff(trait_data$species, acoust_data$species)
# setdiff(acoust_data$species, trait_data$species)
# setdiff(acoust_data$species,  phy$tip.label)
# 

acoust_data$peakf_sc <- scale(acoust_data$meanpeakf)
# acoust_data$freq_by_time_sc <- scale(acoust_data$freq_by_time)
merged_data <- merge(trait_data, acoust_data)

head(merged_data)

species	BodyMass_g_Phyllostomid_Book	mean_forearm_mm	AdultHeadBodyLen_mm	Diet	Precip_Mean_mm	Foliage	Tent	Hollow tree	Cave	agg1-10	agg11-25	agg25-100	agg100	UMAP1	UMAP2	PC1	PC2	guild_umap	guild_pca	log_mass	log_mass_centered	weight.forearm.ratio	sound.files	pc1	duration	meanfreq	sd	freq.median	freq.Q25	freq.Q75	freq.IQR	time.median	time.Q25	time.Q75	time.IQR	skew	kurt	sp.ent	time.ent	entropy	sfm	meandom	mindom	maxdom	dfrange	modindx	startdom	enddom	dfslope	meanpeakf	freq_by_time	peakf_sc
Ametrida_centurio	10.15	28.76	43.5	frugivoro	159.06	NA	NA	NA	NA	NA	NA	NA	NA	1.183465	-2.894392	1.654208	2.46953	1	1	2.317474	-0.5855311	0.3529207	Phyllost_Stenodermatinae_Ametrida-centurio_2001_MD2-70_F-Leblanc.wav-song_Ame.cen.MB.001.txt	-1.6932120	0.0012589	90.57941	8.583930	90.19244	85.39639	94.98848	9.59209	0.0006298	0.0003599	0.0008997	0.0005398	0.7967709	2.304932	0.9404691	0.8049147	0.7569975	0.7680777	90.86000	77.7	107.1	29.4	1.142857	107.1	81.9	-20016.97	90.3	17769.84	0.6063032
Ametrida_centurio	10.15	28.76	43.5	frugivoro	159.06	NA	NA	NA	NA	NA	NA	NA	NA	1.183465	-2.894392	1.654208	2.46953	1	1	2.317474	-0.5855311	0.3529207	Phyllost_Stenodermatinae_Ametrida-centurio_2001_MD2-70_F-Leblanc.wav-song_Ame.cen.MB.001.txt	-1.8038959	0.0012808	91.47430	9.499180	90.83448	85.30384	97.15521	11.85137	0.0006405	0.0002745	0.0009150	0.0006405	0.8980049	2.699468	0.9450159	0.8081140	0.7636806	0.7809586	92.82000	77.7	111.3	33.6	1.000000	111.3	77.7	-26232.72	90.3	18503.99	0.6063032
Ametrida_centurio	10.15	28.76	43.5	frugivoro	159.06	NA	NA	NA	NA	NA	NA	NA	NA	1.183465	-2.894392	1.654208	2.46953	1	1	2.317474	-0.5855311	0.3529207	Phyllost_Stenodermatinae_Ametrida-centurio_2001_MD2-70_F-Leblanc.wav-song_Ame.cen.MB.001.txt	-0.6190427	0.0032170	102.85754	21.082184	99.82434	88.30942	114.14017	25.83075	0.0019614	0.0009415	0.0025890	0.0016476	0.8707555	3.010304	0.9595291	0.8852966	0.8494679	0.7644487	104.50000	69.3	153.3	84.0	6.550000	153.3	73.5	-24805.71	86.1	15678.17	0.4056895
Ametrida_centurio	10.15	28.76	43.5	frugivoro	159.06	NA	NA	NA	NA	NA	NA	NA	NA	1.183465	-2.894392	1.654208	2.46953	1	1	2.317474	-0.5855311	0.3529207	Phyllost_Stenodermatinae_Ametrida-centurio_2001_MD2-70_F-Leblanc.wav-song_Ame.cen.MB.001.txt	-1.9852120	0.0015489	93.32498	9.270191	94.11441	86.26659	100.00028	13.73369	0.0008611	0.0004306	0.0012056	0.0007750	0.3936817	2.111864	0.9531655	0.8186791	0.7803367	0.7827092	96.71053	81.9	111.3	29.4	1.000000	111.3	81.9	-18980.88	102.9	17720.89	1.2081441
Ametrida_centurio	10.15	28.76	43.5	frugivoro	159.06	NA	NA	NA	NA	NA	NA	NA	NA	1.183465	-2.894392	1.654208	2.46953	1	1	2.317474	-0.5855311	0.3529207	Phyllost_Stenodermatinae_Ametrida-centurio_2001_MD2-70_F-Leblanc.wav-song_Ame.cen.MB.001.txt	-0.0254849	0.0019362	82.50245	8.727464	81.85055	76.64307	89.14104	12.49797	0.0010099	0.0005891	0.0013466	0.0007575	0.8957588	3.250050	0.9472120	0.8264910	0.7828622	0.7594154	85.22500	65.1	102.9	37.8	1.111111	102.9	69.3	-17353.96	81.9	16499.98	0.2050759
Ametrida_centurio	10.15	28.76	43.5	frugivoro	159.06	NA	NA	NA	NA	NA	NA	NA	NA	1.183465	-2.894392	1.654208	2.46953	1	1	2.317474	-0.5855311	0.3529207	Phyllost_Stenodermatinae_Ametrida-centurio_2001_MD2-70_F-Leblanc.wav-song_Ame.cen.MB.001.txt	-0.1530456	0.0016681	80.64948	8.110954	81.59952	74.95068	86.43504	11.48436	0.0009167	0.0005000	0.0012500	0.0007500	0.5324313	2.156171	0.9503327	0.8199970	0.7792700	0.7760111	83.58000	69.3	94.5	25.2	1.000000	94.5	69.3	-15107.23	86.1	15312.48	0.4056895

Code

# nrow(trait_data)
# nrow(acoust_data)
# nrow(merged_data)

write.csv(merged_data, "./data/processed/acoustic_trait_guild_data.csv", row.names = FALSE)

11 Data exploration

The Bodyweight-forearm Ratio in bats can be calculated by dividing the body weight of the bat by its forearm length. This ratio is often used as a body condition index (BCI) to estimate the relative size of fat stores in bats [1] [2]. It is assumed that bats with longer forearms weigh more than bats with shorter forearms, and therefore a higher BCI indicates better body condition [3]. However, it has been found that the relationship between body mass and forearm length is weak in bats, and BCI is essentially equivalent to dividing body mass by a constant [4]. Despite this, BCI can still be an effective predictor of fat mass in bats, as fat mass is strongly correlated with body mass [5]. Therefore, researchers can calculate the Bodyweight-forearm Ratio as a simple and cost-effective way to estimate body condition in bats

Code

cols <- viridis(10)

agg_dat <- aggregate(log_mass ~ Diet, trait_data, mean)
agg_dat$n <- sapply(1:nrow(agg_dat), function(x) sum(trait_data$Diet == agg_dat$Diet[x]))
agg_dat$n.labels <- paste("n =", agg_dat$n)
agg_dat$Diet <- factor(agg_dat$Diet)

# raincoud plot:
fill_color <- adjustcolor("#e85307", 0.6)

ggplot(trait_data, aes(x = Diet, y = log_mass)) +
  # add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    fill = fill_color,
    alpha = 0.5,
    # custom bandwidth
    adjust = .5,
    # adjust height
    width = .6,
    .width = 0,
    # move geom to the cright
    justification = -.2,
    point_colour = NA
  ) +
  geom_boxplot(fill = fill_color,
    width = .15,
    # remove outliers
    outlier.shape = NA # `outlier.shape = NA` works as well
  ) +
  # add justified jitter from the {gghalves} package
  gghalves::geom_half_point(
    color = fill_color,
    # draw jitter on the left
    side = "l",
    # control range of jitter
    range_scale = .4,
    # add some transparency
    alpha = .5,
    transformation = ggplot2::position_jitter(height = 0)  
  ) +
   # ylim(c(-30, 310)) #+
  geom_text(data = agg_dat, aes(y = rep(1.6, 6), x = Diet, label = n.labels), nudge_x = -0.13, size = 6)

Code

agg_dat <- aggregate(weight.forearm.ratio ~ Diet, trait_data, mean)
agg_dat$n <- sapply(1:nrow(agg_dat), function(x) sum(trait_data$Diet == agg_dat$Diet[x]))
agg_dat$n.labels <- paste("n =", agg_dat$n)
agg_dat$Diet <- factor(agg_dat$Diet)

# raincoud plot:
fill_color <- adjustcolor("#e85307", 0.6)

ggplot(trait_data, aes(x = Diet, y = weight.forearm.ratio)) +
  # add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    fill = fill_color,
    alpha = 0.5,
    # custom bandwidth
    adjust = .5,
    # adjust height
    width = .6,
    .width = 0,
    # move geom to the cright
    justification = -.2,
    point_colour = NA
  ) +
  geom_boxplot(fill = fill_color,
    width = .15,
    # remove outliers
    outlier.shape = NA # `outlier.shape = NA` works as well
  ) +
  # add justified jitter from the {gghalves} package
  gghalves::geom_half_point(
    color = fill_color,
    # draw jitter on the left
    side = "l",
    # control range of jitter
    range_scale = .4,
    # add some transparency
    alpha = .5,
    transformation = ggplot2::position_jitter(height = 0)  
  ) +
   # ylim(c(-30, 310)) #+
  geom_text(data = agg_dat, aes(y = rep(0, 6), x = Diet, label = n.labels), nudge_x = -0.13, size = 6)

Code

agg_dat <- aggregate(mean_forearm_mm ~ Diet, trait_data, mean)
agg_dat$n <- sapply(1:nrow(agg_dat), function(x) sum(trait_data$Diet == agg_dat$Diet[x]))
agg_dat$n.labels <- paste("n =", agg_dat$n)
agg_dat$Diet <- factor(agg_dat$Diet)

ggplot(trait_data, aes(x = Diet, y = mean_forearm_mm)) +
  # add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    fill = fill_color,
    alpha = 0.5,
    # custom bandwidth
    adjust = .5,
    # adjust height
    width = .6,
    .width = 0,
    # move geom to the cright
    justification = -.2,
    point_colour = NA
  ) +
  geom_boxplot(fill = fill_color,
    width = .15,
    # remove outliers
    outlier.shape = NA # `outlier.shape = NA` works as well
  ) +
  # add justified jitter from the {gghalves} package
  gghalves::geom_half_point(
    color = fill_color,
    # draw jitter on the left
    side = "l",
    # control range of jitter
    range_scale = .4,
    # add some transparency
    alpha = .5,
    transformation = ggplot2::position_jitter(height = 0)  
  ) +
   # ylim(c(-30, 310)) #+
  geom_text(data = agg_dat, aes(y = rep(min(trait_data$mean_forearm_mm), 6), x = Diet, label = n.labels), nudge_x = -0.13, size = 6)

Code

agg_dat <- aggregate(pc1 ~ Diet, merged_data, mean)
agg_dat$n <- sapply(1:nrow(agg_dat), function(x) length(unique(merged_data$sound.files[merged_data$Diet == agg_dat$Diet[x]])))
agg_dat$n.labels <- paste("n =", agg_dat$n)
agg_dat$Diet <- factor(agg_dat$Diet)

# raincoud plot:
fill_color <- adjustcolor("#e85307", 0.6)

ggplot(merged_data, aes(x = Diet, y = pc1)) +
  # add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    fill = fill_color,
    alpha = 0.5,
    # custom bandwidth
    adjust = .5,
    # adjust height
    width = .6,
    .width = 0,
    # move geom to the cright
    justification = -.2,
    point_colour = NA
  ) +
  geom_boxplot(fill = fill_color,
    width = .15,
    # remove outliers
    outlier.shape = NA # `outlier.shape = NA` works as well
  ) +
  # add justified jitter from the {gghalves} package
  gghalves::geom_half_point(
    color = fill_color,
    # draw jitter on the left
    side = "l",
    # control range of jitter
    range_scale = .4,
    # add some transparency
    alpha = .5,
    transformation = ggplot2::position_jitter(height = 0)  
  ) +
   # ylim(c(-30, 310)) #+
  geom_text(data = agg_dat, aes(y = rep(-12, 6), x = Diet, label = n.labels), nudge_x = -0.13, size = 6) #+

Code

   # scale_x_discrete(labels=c("Control" = "Noise control", "Sound vision" = "Sound & vision", "Vision" = "Vision", "Lessen input" = "Lessen input")) +
  # labs(x = "Sensory input       ", y = "Time to enter roost (s)") + theme(axis.text.x = element_text(angle = 15, hjust = 1))

Code

cols <- viridis(10)
ggplot(data = merged_data, aes(x = as.factor(guild_umap), y = pc1, fill = as.factor(guild_umap))) +
  scale_fill_viridis_d(alpha = 0.5) +
  geom_violin() +
  geom_boxplot() +
  # scale_color_viridis_d(alpha = 0.5, begin = 0.2, end = 0.8) +
     theme(axis.text.x = element_text(size = 14, angle=45, hjust = 1)) +
  geom_hline(yintercept = 20, lty = 4, col = "gray")

Code

ggplot(data = merged_data, aes(x = UMAP1, y = pc1, fill = as.factor(guild_umap))) +
  scale_fill_viridis_d(alpha = 0.3) +
  geom_point(shape = 21, size = 3) +
  # scale_color_viridis_d(alpha = 0.5, begin = 0.2, end = 0.8) +
     theme(axis.text.x = element_text(size = 14, angle=45, hjust = 1)) +
  geom_hline(yintercept = 20, lty = 4, col = "gray")

Code

ggplot(data = merged_data, aes(x = UMAP2, y = pc1, fill = as.factor(guild_umap))) +
  scale_fill_viridis_d(alpha = 0.3) +
  geom_point(shape = 21, size = 3) +
  # scale_color_viridis_d(alpha = 0.5, begin = 0.2, end = 0.8) +
     theme(axis.text.x = element_text(size = 14, angle=45, hjust = 1)) +
  geom_hline(yintercept = 20, lty = 4, col = "gray")

Code

agg_pc1_umap <- aggregate(cbind(pc1, UMAP1, UMAP2, guild_umap) ~ species, merged_data, mean)

ggplot(data = agg_pc1_umap, aes(x = UMAP1, y = pc1, fill = as.factor(guild_umap))) +
  scale_fill_viridis_d(alpha = 0.3) +
  geom_point(shape = 21, size = 5) +
  # scale_color_viridis_d(alpha = 0.5, begin = 0.2, end = 0.8) +
     theme(axis.text.x = element_text(size = 14, angle=45, hjust = 1)) +
  geom_hline(yintercept = 20, lty = 4, col = "gray")

11.1 Phylogenetic signal

Code

trait_data$body.mass.log <- log(trait_data$AdultBodyMass_g)

trait_data$Diet[trait_data$species == "Hylonycteris_underwoodi"] <- "nectarivoro"
trait_data <- as.data.frame(trait_data)

rownames(trait_data) <- trait_data$species
traits <- names(trait_data)[-1]
traits <- c("AdultBodyMass_g", "mean_forearm_mm", "DietBreadth", "GR_MidRangeLat_dd", "TrophicLevel", "Precip_Mean_mm", "Temp_Mean_01degC")

A <- ape::vcv.phylo(phy)

lambdas <- sapply(traits, function(x){

  trait <- trait_data[,x]

names(trait) <- trait_data$species
ps <- phytools::phylosig(tree = phy, x = trait, method = "lambda")

form <- paste(x, "~ + (1|gr(species, cov = A))")
mod <- brm(
  formula = form,
  data = trait_data,
  chains = 1,
  family = gaussian(),
  data2 = list(A = A),
  prior = c(
    prior(normal(0, 50), "Intercept"),
    prior(student_t(3, 0, 20), "sd"),
    prior(student_t(3, 0, 20), "sigma")
  )
)

hyp <- "sd_species__Intercept^2 / (sd_species__Intercept^2 + sigma^2) = 0"

hyp <-hypothesis(mod, hyp, class = NULL)


return(data.frame(trait = x, freq.lambda = ps$lambda, bayes.lambda = hyp$hypothesis$Estimate))
}
)

res <- do.call(rbind, lambdas)


mat <- matrix(res, ncol = 3, byrow = T)

Code

A <- ape::vcv.phylo(phylo)

m1 <- brm(
  phen ~ cofactor + (1|gr(phylo, cov = A)),
  data = data_simple,
  family = gaussian(),
  data2 = list(A = A),
  prior = c(
    prior(normal(0, 10), "b"),
    prior(normal(0, 50), "Intercept"),
    prior(student_t(3, 0, 20), "sd"),
    prior(student_t(3, 0, 20), "sigma")
  )
)

A <- ape::vcv.phylo(phylo)
model_simple <- brm(
  phen ~ cofactor + (1|gr(phylo, cov = A)),
  data = data_simple,
  chains = 1,
  family = gaussian(),
  data2 = list(A = A),
  prior = c(
    prior(normal(0, 10), "b"),
    prior(normal(0, 50), "Intercept"),
    prior(student_t(3, 0, 20), "sd"),
    prior(student_t(3, 0, 20), "sigma")
  )
)


hyp <- "sd_phylo__Intercept^2 / (sd_phylo__Intercept^2 + sigma^2) = 0"
(hyp <- hypothesis(model_simple, hyp, class = NULL))

12 Stats

Code

phy <- read.tree("./data/processed/trimmed_phylogeny_72_species_chiroptera.no_outgroups.absolute.tre")

A <- ape::vcv.phylo(phy)
merged_data$phylo <- merged_data$species

iter <- 1000
chains <- 1

Code

# mod <- readRDS("./data/processed/regression_model_diet_vs_pc1.RDS")
# 
# cntrs <- as.data.frame(emmeans(mod, pairwise ~ Diet)$contrasts)
# names(cntrs)[3:4] <- c("l-95% CI", "u-95% CI")
# 
# coef_table <- html_format_coef_table(cntrs, highlight = TRUE, fill = viridis(10)[7])
# 
# 
# print(coef_table)

Code

models <- c(
  "400" = "pc1 ~ log_mass",
  #guild UMAP
  "400" = "pc1 ~ guild_umap + log_mass_centered",
"400" = "modindx ~ guild_umap + log_mass_centered",
"400" = "freq_by_time ~ guild_umap + log_mass_centered",
"1000" = "duration ~ guild_umap + log_mass_centered",
"400" = "meanpeakf ~ guild_umap + log_mass_centered",
"400" = "freq.IQR ~ guild_umap + log_mass_centered",
"1000" = "time.IQR ~ guild_umap + log_mass_centered",
"400" = "pc1 ~ umap + log_mass_centered",

# functional trait dimension UMAP1
"400" = "pc1 ~ UMAP1 + log_mass_centered",
"400" = "modindx ~ UMAP1 + log_mass_centered",
"400" = "freq_by_time ~ UMAP1 + log_mass_centered",
"1000" = "duration ~ UMAP1 + log_mass_centered",
"400" = "meanpeakf ~ UMAP1 + log_mass_centered",
"400" = "freq.IQR ~ UMAP1 + log_mass_centered",
"400" = "time.IQR ~ UMAP1 + log_mass_centered",

# functional trait dimension UMAP2
"400" = "pc1 ~ UMAP2 + log_mass_centered",
"400" = "modindx ~ UMAP2 + log_mass_centered",
"400" = "freq_by_time ~ UMAP2 + log_mass_centered",
"1000" = "duration ~ UMAP2 + log_mass_centered",
"400" = "meanpeakf ~ UMAP2 + log_mass_centered",
"400" = "freq.IQR ~ UMAP2 + log_mass_centered",
"1000" = "time.IQR ~ UMAP2 + log_mass_centered",

#guild PCA
#   "pc1 ~ guild_pca + log_mass_centered",
# "modindx ~ guild_pca + log_mass_centered",
# "freq_by_time ~ guild_pca + log_mass_centered",
# "duration ~ guild_pca + log_mass_centered",
# "meanpeakf ~ guild_pca + log_mass_centered",
# "freq.IQR ~ guild_pca + log_mass_centered",
# "time.IQR ~ guild_pca + log_mass_centered",
# "pc1 ~ umap + log_mass_centered",

# functional trait dimension UMAP1
"400" = "modindx ~ PC1 + log_mass_centered",
"400" = "freq_by_time ~ PC1 + log_mass_centered",
"1000" = "duration ~ PC1 + log_mass_centered",
"400" = "meanpeakf ~ PC1 + log_mass_centered",
"400" = "freq.IQR ~ PC1 + log_mass_centered",
"1000" = "time.IQR ~ PC1 + log_mass_centered",

# Diet
"400" = "pc1 ~ Diet + log_mass_centered",
"400" = "modindx ~ Diet + log_mass_centered",
"400" = "freq_by_time ~ Diet + log_mass_centered",
"400" = "duration ~ Diet + log_mass_centered",
"400" = "meanpeakf ~ Diet + log_mass_centered",
"400" = "freq.IQR ~ Diet + log_mass_centered",
"1000" = "time.IQR ~ Diet + log_mass_centered"
)

Code

A <- ape::vcv.phylo(phy)
merged_data$phylo <- merged_data$species

chains <- 1
iter <- 400
mods <- gsub("~", "by", models, fixed = TRUE)
mods <- gsub("+", "plus", mods, fixed = TRUE)
rds_names <- paste0("./data/processed/regression_models/", mods, ".RDS")
list_models <- list.files(path = "./data/processed/regression_models", pattern = paste(mods, collapse = "|"), full.names = TRUE)

submodels <- models[!rds_names %in% list_models]
submodels <- models[names(models) == "1000"]

out <- warbleR:::pblapply_wrblr_int(seq_along(submodels), cl = 5, function(y){
  
  x <- submodels[y]
  iter <- as.numeric(names(submodels)[y])
  
    # scale response
  resp <- strsplit(x, " ~")[[1]][1]
  
    if (resp != "pc1"){

  scale_resp <- scale(merged_data[, resp])
  
  merged_data$scale_resp <- scale_resp
  
  colnames(merged_data)[ncol(merged_data)] <- paste0(resp, "_sc")
  
  x <- gsub(resp, paste0(resp, "_sc"), x)
  }
  
  form <- paste(x, "+ (1|species / sound.files) + (1|gr(phylo, cov = A))")
  
  rds_name <- paste0("./data/processed/regression_models/", x, ".RDS")
   rds_name <- gsub("~", "by", rds_name, fixed = TRUE)
   rds_name <- gsub("+", "plus", rds_name, fixed = TRUE)
  
  try(brm(
  formula = form,
  data = merged_data,
  chains = chains,
  cores = chains,
  iter = iter,
  family = gaussian(),
  
  data2 = list(A = A),
  prior = c(
    prior(normal(0,10), "b"),
    prior(normal(0, 50), "Intercept"),
    prior(student_t(3, 0, 20), "sd"),
    prior(student_t(3, 0, 20), "sigma")
  ),
   control = list(adapt_delta = 0.99, max_treedepth = 15), 
  file = rds_name,
    file_refit = "always"
), silent = TRUE)
  
})

13 Regression model results

Code

rds_names <- paste0("./data/processed/regression_models/", models, ".RDS")
rds_names <- gsub("~", "by", rds_names, fixed = TRUE)
rds_names <- gsub("+", "plus", rds_names, fixed = TRUE)

models <- gsub("~", "by", models, fixed = TRUE)
models <- gsub("+", "plus", models, fixed = TRUE)

# list_models <- list.files(path = "./data/processed/regression_models", pattern = "RDS", full.names = TRUE)


list_models <- list.files(path = "./data/processed/regression_models", pattern = paste(models, collapse = "|"), full.names = TRUE)



for(i in list_models)
{
extended_summary(read.file = i,
    n.posterior = 1000, fill = viridis(10)[7], trace.palette = my.viridis,
    highlight = TRUE, remove.intercepts = TRUE)
}

13.1 duration by guild_umap plus log_mass_centered

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	duration ~ guild_umap + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	0	1	29.59	40.709	1012035855

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_guild_umap2	-0.001	-0.002	0.000	1.053	29.590	71.403
b_guild_umap3	-0.001	-0.002	0.000	1.021	54.801	95.185
b_guild_umap4	0.000	-0.001	0.001	1.012	50.979	72.103
b_guild_umap5	-0.001	-0.003	0.000	1	48.838	113.523
b_guild_umap6	0.000	0.000	0.001	0.995	65.715	78.955
b_log_mass_centered	0.001	0.000	0.001	1	32.059	40.709

13.2 duration by PC1 plus log_mass_centered

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	duration ~ PC1 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	0	1	12.388	25.898	501037413

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_PC1	0	0	0	1.044	16.760	25.898
b_log_mass_centered	0	0	0	1.063	12.388	50.958

13.3 duration by UMAP1 plus log_mass_centered

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	duration ~ UMAP1 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	0	0	12.892	40.945	211055999

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_UMAP1	0	0	0.000	1.039	12.892	40.945
b_log_mass_centered	0	0	0.001	1.03	34.253	78.067

13.4 duration by UMAP2 plus log_mass_centered

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	duration ~ UMAP2 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	0	1	9.003	11.643	736171389

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_UMAP2	0.000	0	0.000	0.997	9.003	12.759
b_log_mass_centered	0.001	0	0.001	1.181	12.338	11.643

13.5 freq_by_time by guild_umap plus log_mass_centered

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	freq_by_time ~ guild_umap + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	33	0	167.02	95.649	1430328635

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_guild_umap2	0.620	-15.242	19.556	1.006	227.849	116.691
b_guild_umap3	-0.411	-17.870	19.892	1.009	206.979	129.321
b_guild_umap4	0.897	-13.398	16.875	1.004	167.020	148.380
b_guild_umap5	0.673	-18.541	17.991	1.008	176.677	95.649
b_guild_umap6	0.254	-23.792	22.758	1.018	208.630	112.774
b_log_mass_centered	-0.623	-19.582	18.280	1.001	246.380	130.475

13.6 freq_by_time by PC1 plus log_mass_centered

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	freq_by_time ~ PC1 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	0	0	142.62	109.736	1603593775

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_PC1	0.344	-18.355	21.011	1.017	142.62	109.736
b_log_mass_centered	-0.460	-16.596	14.769	0.995	321.65	180.787

13.7 freq_by_time by UMAP1 plus log_mass_centered

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	freq_by_time ~ UMAP1 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	0	0	180.799	130.76	931547068

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_UMAP1	1.470	-20.306	21.792	0.995	184.647	130.760
b_log_mass_centered	0.428	-20.903	21.665	1.015	180.799	136.785

13.8 freq_by_time by UMAP2 plus log_mass_centered

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	freq_by_time ~ UMAP2 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	0	0	297.83	118.134	1513858676

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_UMAP2	-0.685	-18.566	19.509	1.02	297.830	118.134
b_log_mass_centered	0.265	-15.318	16.208	0.999	390.743	190.342

13.9 freq.IQR by guild_umap plus log_mass_centered

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	freq.IQR ~ guild_umap + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	0	0	146.023	86.9	1420864103

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_guild_umap2	-3.218	-7.463	1.019	1.003	283.167	217.056
b_guild_umap3	-3.430	-7.041	0.979	1.008	184.236	130.535
b_guild_umap4	-3.675	-7.501	0.720	1.003	231.903	118.977
b_guild_umap5	-1.366	-7.689	5.438	0.996	146.023	86.900
b_guild_umap6	-4.213	-8.975	0.070	0.996	204.411	147.713
b_log_mass_centered	-4.511	-8.709	-1.020	1.007	234.462	116.691

13.10 freq.IQR by PC1 plus log_mass_centered

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	freq.IQR ~ PC1 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	0	0	170.413	89.77	1696722880

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_PC1	0.147	-0.611	0.921	1.002	170.413	89.770
b_log_mass_centered	-4.805	-7.188	-2.825	1.003	246.851	154.496

13.11 freq.IQR by UMAP1 plus log_mass_centered

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	freq.IQR ~ UMAP1 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	0	0	72.823	54.185	582078136

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_UMAP1	0.382	-0.591	1.174	0.996	108.401	148.380
b_log_mass_centered	-4.481	-7.394	-2.026	1.025	72.823	54.185

13.12 freq.IQR by UMAP2 plus log_mass_centered

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	freq.IQR ~ UMAP2 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	0	0	48.316	104.961	89945138

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_UMAP2	-0.097	-0.549	0.470	1.002	48.316	104.961
b_log_mass_centered	-4.619	-6.973	-1.801	0.996	65.688	112.316

13.13 modindx by guild_umap plus log_mass_centered

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	modindx ~ guild_umap + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	0	0	88.126	103.271	2039654840

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_guild_umap2	-0.293	-1.138	0.607	1.014	129.963	103.271
b_guild_umap3	-0.555	-1.202	0.171	1.01	114.320	149.077
b_guild_umap4	-0.091	-0.899	0.663	1.02	121.873	150.415
b_guild_umap5	-1.108	-2.117	0.068	0.998	128.555	129.321
b_guild_umap6	0.029	-0.836	0.872	1.031	88.126	128.893
b_log_mass_centered	0.708	0.064	1.246	0.998	140.166	148.612

13.14 modindx by PC1 plus log_mass_centered

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	modindx ~ PC1 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	0	0	132.745	142.42	704416026

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_PC1	-0.054	-0.144	0.063	1.005	132.745	142.420
b_log_mass_centered	0.236	-0.114	0.595	0.995	172.002	160.145

13.15 modindx by UMAP1 plus log_mass_centered

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	modindx ~ UMAP1 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	0	0	144.327	105.272	400793765

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_UMAP1	-0.041	-0.200	0.105	1.016	144.327	117.800
b_log_mass_centered	0.242	-0.129	0.569	0.997	204.706	105.272

13.16 modindx by UMAP2 plus log_mass_centered

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	modindx ~ UMAP2 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	0	0	82.116	98.183	1789545808

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_UMAP2	-0.004	-0.112	0.074	1.016	153.185	127.637
b_log_mass_centered	0.269	-0.117	0.586	1.012	82.116	98.183

13.17 pc1 by Diet plus log_mass_centered

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	pc1 ~ Diet + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	0	0	125.713	115.819	1416271543

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Dietfrugivoro	1.861	-2.277	5.817	1.003	150.062	153.468
b_Diethematofago	2.821	-2.460	7.946	1.006	129.156	152.858
b_Dietinsectivoro	3.681	-0.142	7.660	1.006	125.713	115.819
b_Dietnectarivoro	3.864	-0.115	7.870	0.997	125.922	182.726
b_Dietomnivoro	2.582	-1.923	6.411	1.002	176.794	138.037
b_log_mass_centered	2.478	1.083	3.877	1.005	218.568	176.756

13.18 pc1 by guild_umap plus log_mass_centered

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	pc1 ~ guild_umap + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	4	1	200	0	0	530.273	308.424	1848056092

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_guild_umap2	-0.991	-3.495	1.664	1.01	602.451	308.424
b_guild_umap3	-1.668	-3.773	0.371	1.009	817.708	599.896
b_guild_umap4	-0.346	-2.586	1.935	1.004	581.214	460.788
b_guild_umap5	-1.883	-5.045	1.006	1.007	530.273	550.854
b_guild_umap6	1.112	-0.988	3.093	1.003	843.765	605.135
b_log_mass_centered	2.732	1.006	4.421	1.005	718.128	492.175

13.19 pc1 by log_mass

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	pc1 ~ log_mass + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	4	1	200	0	0	1517.812	622.247	1784415959

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_log_mass	1.832	0.649	3.023	1.003	1517.812	622.247

13.20 pc1 by UMAP1 plus log_mass_centered

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	pc1 ~ UMAP1 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	0	0	184.415	154.069	1568764005

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_UMAP1	0.244	-0.242	0.714	1.001	184.415	154.069
b_log_mass_centered	1.982	0.619	3.161	1.003	200.216	174.002

13.21 pc1 by UMAP2 plus log_mass_centered

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	pc1 ~ UMAP2 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	0	0	215.604	137.043	1102094814

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_UMAP2	-0.173	-0.463	0.100	1.009	215.604	178.629
b_log_mass_centered	2.109	0.965	3.305	0.996	220.660	137.043

13.22 time.IQR by guild_umap plus log_mass_centered

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	time.IQR ~ guild_umap + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	0	5	1.401	15.48	1945099016

	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_guild_umap2	-0.001	0.000	2.022	1.401	18.505
b_guild_umap3	-0.001	0.000	1.2	9.795	35.608
b_guild_umap4	0.000	0.000	1.491	2.012	25.169
b_guild_umap5	-0.001	0.000	1.026	8.097	33.546
b_guild_umap6	0.000	0.001	1.395	2.484	38.338
b_log_mass_centered	0.000	0.000	1.213	3.537	15.480

13.23 time.IQR by UMAP1 plus log_mass_centered

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	time.IQR ~ UMAP1 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	0	2	3.194	21.564	349114049

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_UMAP1	0	0	0	1.128	9.909	21.564
b_log_mass_centered	0	0	0	1.254	3.194	46.313

13.24 time.IQR by UMAP2 plus log_mass_centered

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	time.IQR ~ UMAP2 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	0	1	1.842	20.687	60477512

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_UMAP2	0	0	0	1.562	1.842	20.687
b_log_mass_centered	0	0	0	1.04	11.224	30.232

13.25 Models summary

Code

check_rds_fits(path = "./data/processed/regression_models")

fit	priors	formula	iterations	chains	thinning	warmup	n_parameters	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
duration by guild_umap plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	duration ~ guild_umap + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	566	0	171	1.414179	11.52363	1012035855
duration by PC1 plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	duration ~ PC1 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	562	0	189	1.452993	11.54226	501037413
duration by UMAP1 plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	duration ~ UMAP1 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	562	0	125	1.800062	12.24317	211055999
duration by UMAP2 plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	duration ~ UMAP2 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	562	0	429	1.379944	11.52363	736171389
duration_sc by Diet plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	duration_sc ~ Diet + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	566	0	7	16.238699	42.41654	231008385
duration_sc by guild_umap plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	duration_sc ~ guild_umap + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	1000	4	1	500	566	0	0	344.804294	620.67792	1803259205
duration_sc by PC1 plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	duration_sc ~ PC1 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	1000	1	1	500	562	0	2	61.698646	102.39185	753365269
duration_sc by UMAP1 plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	duration_sc ~ UMAP1 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	1000	1	1	500	562	0	0	87.188235	123.73231	925036368
duration_sc by UMAP2 plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	duration_sc ~ UMAP2 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	1000	1	1	500	562	0	0	96.471010	105.98151	1719960287
freq_by_time by guild_umap plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	freq_by_time ~ guild_umap + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	566	33	12	10.809682	11.52363	1430328635
freq_by_time by PC1 plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	freq_by_time ~ PC1 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	562	0	6	12.925361	32.44265	1603593775
freq_by_time by UMAP1 plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	freq_by_time ~ UMAP1 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	562	0	18	10.455200	12.94361	931547068
freq_by_time by UMAP2 plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	freq_by_time ~ UMAP2 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	562	0	5	22.702389	37.86705	1513858676
freq_by_time_sc by Diet plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	freq_by_time_sc ~ Diet + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	566	0	18	16.688169	25.22704	1538358049
freq_by_time_sc by guild_umap plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	freq_by_time_sc ~ guild_umap + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	4	1	200	566	0	0	67.925228	171.44743	223002341
freq.IQR by guild_umap plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	freq.IQR ~ guild_umap + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	566	0	40	10.920752	21.41411	1420864103
freq.IQR by PC1 plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	freq.IQR ~ PC1 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	562	0	48	6.420637	15.88104	1696722880
freq.IQR by UMAP1 plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	freq.IQR ~ UMAP1 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	562	0	33	6.935264	35.04760	582078136
freq.IQR by UMAP2 plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	freq.IQR ~ UMAP2 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	562	0	21	16.049081	21.27742	89945138
freq.IQR_sc by Diet plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	freq.IQR_sc ~ Diet + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	566	0	42	12.173199	17.97004	1663669132
meanpeakf_sc by Diet plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	meanpeakf_sc ~ Diet + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	566	0	11	10.155642	32.64134	990740690
meanpeakf_sc by guild_umap plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	meanpeakf_sc ~ guild_umap + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	566	0	13	23.201968	18.52600	1359788322
meanpeakf_sc by PC1 plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	meanpeakf_sc ~ PC1 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	562	0	7	44.796910	28.21815	1670378914
meanpeakf_sc by UMAP1 plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	meanpeakf_sc ~ UMAP1 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	562	0	3	35.171879	17.02180	1123658327
meanpeakf_sc by UMAP2 plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	meanpeakf_sc ~ UMAP2 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	562	0	21	17.018213	23.93021	1429849328
modindx by guild_umap plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	modindx ~ guild_umap + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	566	0	17	6.619650	25.41942	2039654840
modindx by PC1 plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	modindx ~ PC1 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	562	0	11	31.446701	26.07006	704416026
modindx by UMAP1 plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	modindx ~ UMAP1 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	562	0	46	7.099505	17.45561	400793765
modindx by UMAP2 plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	modindx ~ UMAP2 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	562	0	10	16.058580	20.02476	1789545808
modindx_sc by Diet plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	modindx_sc ~ Diet + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	566	0	8	22.764186	20.79081	1465706639
modindx_sc by guild_umap plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	modindx_sc ~ guild_umap + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	4	1	200	566	0	0	132.906590	173.90231	535014709
pc1 by Diet plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	pc1 ~ Diet + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	566	0	25	6.684182	26.61827	1416271543
pc1 by guild_umap plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	pc1 ~ guild_umap + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	4	1	200	566	0	0	138.315809	134.83295	1848056092
pc1 by log_mass.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	pc1 ~ log_mass + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	4	1	200	561	0	0	187.914234	308.42009	1784415959
pc1 by UMAP1 plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	pc1 ~ UMAP1 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	562	0	3	60.242819	17.02180	1568764005
pc1 by UMAP2 plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	pc1 ~ UMAP2 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	562	0	6	23.455441	14.26978	1102094814
time.IQR by guild_umap plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	time.IQR ~ guild_umap + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	566	0	462	1.362101	11.52363	1945099016
time.IQR by UMAP1 plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	time.IQR ~ UMAP1 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	562	0	433	1.381603	11.61195	349114049
time.IQR by UMAP2 plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	time.IQR ~ UMAP2 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	400	1	1	200	562	0	533	1.340174	11.52363	60477512
time.IQR_sc by Diet plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	time.IQR_sc ~ Diet + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	1000	1	1	500	566	0	0	76.429207	104.36029	1467582519
time.IQR_sc by guild_umap plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	time.IQR_sc ~ guild_umap + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	1000	1	1	500	566	0	0	64.936333	143.64890	658104680
time.IQR_sc by PC1 plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	time.IQR_sc ~ PC1 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	1000	1	1	500	562	0	1	42.851792	91.78522	247987981
time.IQR_sc by UMAP2 plus log_mass_centered.RDS	b-normal(0, 10) Intercept-normal(0, 50) sd-student_t(3, 0, 20) sigma-student_t(3, 0, 20)	time.IQR_sc ~ UMAP2 + log_mass_centered + (1 \| species/sound.files) + (1 \| gr(phylo, cov = A))	1000	1	1	500	562	0	0	130.887791	92.59945	549692700

Session information

R version 4.2.2 Patched (2022-11-10 r83330)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.5 LTS

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

locale:
 [1] LC_CTYPE=es_ES.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=es_ES.UTF-8    
 [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=es_ES.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       

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

other attached packages:
 [1] ggtree_3.6.2       cowplot_1.1.1      corrplot_0.92      mclust_6.0.0      
 [5] phytools_1.5-1     maps_3.4.1         umap_0.2.10.0      missMDA_1.18      
 [9] ggdist_3.2.1       emmeans_1.8.6      brms_2.18.0        Rcpp_1.0.11       
[13] sjPlot_2.8.14      lmerTest_3.1-3     lme4_1.1-33        Matrix_1.5-1      
[17] ggplot2_3.4.3      mice_3.16.0        brmsish_1.0.0      suwo_0.1.0        
[21] jsonlite_1.8.7     RCurl_1.98-1.12    ape_5.7-1          readxl_1.4.3      
[25] Rraven_1.0.14      ohun_1.0.0         viridis_0.6.4      viridisLite_0.4.2 
[29] warbleR_1.1.29     NatureSounds_1.0.4 seewave_2.2.2      tuneR_1.4.5       
[33] rprojroot_2.0.3    formatR_1.12       knitr_1.44         kableExtra_1.3.4  
[37] remotes_2.4.2.1   

loaded via a namespace (and not attached):
  [1] estimability_1.4.1      ModelMetrics_1.2.2.2    coda_0.19-4            
  [4] tidyr_1.3.0             clusterGeneration_1.3.7 dygraphs_1.1.1.6       
  [7] data.table_1.14.4       rpart_4.1.19            inline_0.3.19          
 [10] hardhat_1.3.0           doParallel_1.0.17       generics_0.1.3         
 [13] callr_3.7.3             combinat_0.0-8          future_1.28.0          
 [16] proxy_0.4-27            lubridate_1.9.2         webshot_0.5.4          
 [19] xml2_1.3.5              httpuv_1.6.6            StanHeaders_2.21.0-7   
 [22] gower_1.0.1             xfun_0.40               bayesplot_1.9.0        
 [25] evaluate_0.22           promises_1.2.0.1        fansi_1.0.4            
 [28] igraph_1.5.1            DBI_1.1.3               htmlwidgets_1.5.4      
 [31] tensorA_0.36.2          stats4_4.2.2            purrr_1.0.2            
 [34] ellipsis_0.3.2          RSpectra_0.16-1         crosstalk_1.2.0        
 [37] dplyr_1.1.0             backports_1.4.1         signal_0.7-7           
 [40] FactoMineR_2.9          insight_0.19.2          markdown_1.3           
 [43] RcppParallel_5.1.5      vctrs_0.6.3             sjlabelled_1.2.0       
 [46] abind_1.4-5             caret_6.0-94            withr_2.5.1            
 [49] checkmate_2.2.0         treeio_1.23.0           xts_0.12.2             
 [52] prettyunits_1.2.0       mnormt_2.1.1            svglite_2.1.0          
 [55] cluster_2.1.4           lazyeval_0.2.2          crayon_1.5.2           
 [58] recipes_1.0.8           labeling_0.4.3          glmnet_4.1-8           
 [61] pkgconfig_2.0.3         units_0.8-4             nlme_3.1-162           
 [64] nnet_7.3-18             globals_0.16.1          rlang_1.1.1            
 [67] lifecycle_1.0.3         mitml_0.4-5             miniUI_0.1.1.1         
 [70] colourpicker_1.2.0      modelr_0.1.11           cellranger_1.1.0       
 [73] distributional_0.3.1    matrixStats_0.62.0      aplot_0.1.9            
 [76] phangorn_2.11.1         loo_2.5.1               boot_1.3-28            
 [79] zoo_1.8-11              pan_1.9                 base64enc_0.1-3        
 [82] ggridges_0.5.4          processx_3.8.2          png_0.1-8              
 [85] rjson_0.2.21            bitops_1.0-7            pROC_1.18.0            
 [88] KernSmooth_2.23-20      shape_1.4.6             classInt_0.4-10        
 [91] stringr_1.5.0           multcompView_0.1-9      parallelly_1.32.1      
 [94] shinystan_2.6.0         gridGraphics_0.5-1      ggeffects_1.2.2        
 [97] scales_1.2.1            leaps_3.1               gghalves_0.1.4         
[100] magrittr_2.0.3          plyr_1.8.7              threejs_0.3.3          
[103] compiler_4.2.2          rstantools_2.2.0        RColorBrewer_1.1-3     
[106] plotrix_3.8-2           cli_3.6.1               dtw_1.23-1             
[109] listenv_0.8.0           patchwork_1.1.2         pbapply_1.7-2          
[112] ps_1.7.5                Brobdingnag_1.2-9       MASS_7.3-58.2          
[115] tidyselect_1.2.0        fftw_1.0-7              stringi_1.7.12         
[118] highr_0.10              yaml_2.3.7              askpass_1.1            
[121] ggrepel_0.9.3           bridgesampling_1.1-2    grid_4.2.2             
[124] timechange_0.2.0        fastmatch_1.1-3         tools_4.2.2            
[127] future.apply_1.11.0     parallel_4.2.2          rstudioapi_0.14        
[130] foreach_1.5.2           gridExtra_2.3           prodlim_2023.08.28     
[133] posterior_1.3.1         scatterplot3d_0.3-43    farver_2.1.1           
[136] digest_0.6.33           lava_1.7.2.1            shiny_1.7.3            
[139] quadprog_1.5-8          broom_1.0.4             performance_0.10.3     
[142] later_1.3.1             httr_1.4.6              sf_1.0-14              
[145] sjstats_0.18.2          colorspace_2.1-0        rvest_1.0.3            
[148] brio_1.1.3              reticulate_1.26         splines_4.2.2          
[151] yulab.utils_0.0.6       tidytree_0.4.2          expm_0.999-7           
[154] shinythemes_1.2.0       ggplotify_0.1.0         systemfonts_1.0.4      
[157] xtable_1.8-4            nloptr_2.0.3            timeDate_4021.106      
[160] rstan_2.21.7            flashClust_1.01-2       ipred_0.9-14           
[163] ggfun_0.0.9             testthat_3.1.10         R6_2.5.1               
[166] pillar_1.9.0            htmltools_0.5.6         mime_0.12              
[169] glue_1.6.2              fastmap_1.1.1           minqa_1.2.5            
[172] DT_0.26                 class_7.3-21            codetools_0.2-19       
[175] jomo_2.7-6              optimParallel_1.0-2     pkgbuild_1.4.2         
[178] mvtnorm_1.1-3           utf8_1.2.3              lattice_0.20-45        
[181] tibble_3.2.1            numDeriv_2016.8-1.1     gtools_3.9.3           
[184] shinyjs_2.1.0           openssl_2.0.6           survival_3.5-3         
[187] rmarkdown_2.21          munsell_0.5.0           e1071_1.7-13           
[190] iterators_1.0.14        sjmisc_2.8.9            reshape2_1.4.4         
[193] gtable_0.3.4            bayestestR_0.13.1