install/load packages

x <- c("warbleR", "baRulho", "readxl", "pbapply", "ape", "brms", "phytools", "ggtree", "viridis", "ggplot2", "corrplot", "cowplot", "kableExtra", "tidybayes", "cowplot", "ggrepel", "posterior", "ggridges")

aa <- lapply(x, function(y){
  if(!y %in% installed.packages()[,"Package"])  {if(y != "warbleR") install.packages(y) else devtools::install_github("maRce10/warbleR")
}
try(require(y, character.only = T), silent = T)
  })

knitr::opts_knit$set(root.dir = normalizePath(".."))

knitr::opts_chunk$set(dpi = 80, fig.width = 12, fig.height = 8) 

theme_set(theme_classic(base_size = 30))

tibble <- function(x, ...) { 
  x <- kbl(x, digits=4, align= 'c', row.names = FALSE) 
   x <- kable_styling(x, position ="center", full_width = FALSE,  bootstrap_options = c("striped", "hover", "condensed", "responsive")) 
   asis_output(x)
}

registerS3method("knit_print", "data.frame", tibble)
se <- function(x) sd(x) / sqrt(length(x))

summary_brm_model <- function(x, plot = TRUE){
  
   summ <- summary(x)$fixed
  fit <- x$fit  
  betas <- grep("^b_", names(fit@sim$samples[[1]]), value = TRUE)  
  hdis <- t(sapply(betas, function(y)   hdi(fit@sim$samples[[1]][[y]]))
)
  summ_table <- data.frame(summ, hdis, iterations = attr(fit, "stan_args")[[1]]$iter, chains = length(attr(fit, "stan_args")))
  summ_table <- summ_table[rownames(summ_table) != "Intercept", c("Estimate", "Rhat", "Bulk_ESS", "l.95..CI", "u.95..CI", "iterations", "chains")]
  
  summ_table <- as.data.frame(summ_table)  
  summ_table$Rhat <- round(summ_table$Rhat, digits = 3)  
  summ_table$CI_low <- round(unlist(summ_table$l.95..CI), digits = 3)  
  summ_table$CI_high <- round(unlist(summ_table$u.95..CI), digits = 3)  
  summ_table$l.95..CI <- summ_table$u.95..CI <- NULL
  
    out <- lapply(betas, function(y)  data.frame(variable = y, value = sort(fit@sim$samples[[1]][[y]], decreasing = FALSE)))
  
  posteriors <- do.call(rbind, out)
  posteriors <- posteriors[posteriors$variable != "b_Intercept", ]
  
   gg <- ggplot(data = posteriors, aes(y = variable, x = value, fill = stat(quantile))) + 
  geom_density_ridges_gradient(quantile_lines = TRUE, quantile_fun = HDInterval::hdi,  vline_linetype = 2) + 
  theme_classic() + 
     xlim(c(min(summ_table[ , "CI_low"]), max(summ_table[ , "CI_high"]))) +
  geom_vline(xintercept = 0, col = "red") + 
  scale_fill_manual(values = c("transparent", "lightblue", "transparent"), guide = "none") + ggtitle(x$formula)

  
  summ_table$Rhat <- ifelse(summ_table$Rhat > 1.05, cell_spec(summ_table$Rhat, "html", color ="white", background = "red", bold = TRUE,  font_size = 12),  cell_spec(summ_table$Rhat, "html"))
  
  signif <- summ_table[,"CI_low"] * summ_table[,"CI_high"] > 0
  
  df1 <- kbl(summ_table, row.names = TRUE, escape = FALSE, format = "html", digits = 3)
    
  df1 <- row_spec(kable_input = df1,row =  which(summ_table$CI_low * summ_table$CI_high > 0), background = adjustcolor(cols[9], alpha.f = 0.3))

  df1 <- kable_styling(df1, bootstrap_options = c("striped", "hover", "condensed", "responsive"), full_width = FALSE, font_size = 12)
  
  print(df1)
  
  if (plot)
  print(gg)

  }

Color palette

cols <- viridis(10, begin = 0.3, end = 0.8)

print("Closed color:")
## [1] "Closed color:"
cols[10]
## [1] "#7AD151FF"
print("Open color:")
## [1] "Open color:"
cols[1]
## [1] "#35608DFF"

Degradation analysis

Aliging re-recorded signals

master_sels <- read.csv("./data/processed/selection_table_open_and_closed_habitat_shuffled.csv")

wi <- wav_info(path = "./data/raw/recordings")
table(wi$sample.rate)

test_files <- list.files(path = "./data/raw/recordings", pattern = "wav$", ignore.case = TRUE)

test_files <- test_files[grep("TB22", test_files, invert = TRUE)]

# search template on incomplete re-recorded files (those not including start marker)
ends <- search_templates(X = master_sels, template.rows = which(master_sels$orig.sound.file =="end_marker"), test.files = test_files, path = "./data/raw/recordings", parallel = 1)

alg_complete <- align_test_files(X = master_sels, Y = ends, path = "./data/raw/recordings", by.song = TRUE, marker = "end")

alg_complete_df <- align_test_files(X = master_sels, Y = ends, path = "./data/raw/recordings", by.song = TRUE, output = "data.frame")

full_alg <- alg_complete

full_spectrograms(X = alg_complete_df, flim = c(1, 13), sxrow = 6, rows = 15, ovlp = 50, fast.spec = TRUE, path = "./data/raw/recordings", wl = 300, parallel = 3, horizontal = TRUE, suffix = "first_day")

metadat <- read_excel("./data/raw/re-recorded_files_metadata.xlsx")

full_alg$transect <- sapply(full_alg$sound.files, function(x) strsplit(x, "-")[[1]][1])

full_alg$distance <- sapply(full_alg$sound.files, function(x) strsplit(x, "-")[[1]][2])

full_alg$distance <- as.numeric(gsub("d", "", full_alg$distance))

full_alg$playback.habitat <- sapply(full_alg$sound.files, function(x) strsplit(x, "-")[[1]][3])

full_alg$playback.habitat <- sub(".WAV", "", full_alg$playback.habitat, ignore.case = TRUE)

full_alg$signal.type <- full_alg$template

saveRDS(full_alg, "./data/processed/extended_selection_table_re-recorded_aligned_signals.RDS")

Measuring degradation

full_alg <- readRDS("./data/processed/extended_selection_table_re-recorded_aligned_signals.RDS")

# remove duplicates to make it run faster
sub_alg <- full_alg[grep("-dup$", invert = TRUE, full_alg$signal.type), ]

prl <- 20

# test
full_alg <- full_alg[full_alg$transect == "t2" & full_alg$signal.type %in% unique(full_alg$signal.type)[1:3], ]

sub_alg_df <- data.frame(sub_alg)

sub_alg_df$blur.ratio <- blur_ratio(sub_alg, method = 1, pb = TRUE, parallel = prl, output = "data.frame")$blur.ratio

write.csv(sub_alg_df, "./data/processed/degradation_measurements.csv", row.names = FALSE)


sub_alg_df$spectral.blur.ratio  <- spectral_blur_ratio(sub_alg, method = 1, pb = TRUE, parallel = prl, output = "data.frame")$spectral.blur.ratio 

write.csv(sub_alg_df, "./data/processed/degradation_measurements.csv", row.names = FALSE)

sub_alg_df$envelope.correlation <- envelope_correlation(sub_alg, method = 1,  pb = TRUE, parallel = prl, output = "data.frame")$envelope.correlation 

write.csv(sub_alg_df, "./data/processed/degradation_measurements.csv", row.names = FALSE)

sub_alg_df$spectral_correlation <- spectral_correlation(sub_alg, method = 1,  pb = TRUE, parallel = prl, output = "data.frame")$spectral.correlation

write.csv(sub_alg_df, "./data/processed/degradation_measurements.csv", row.names = FALSE)

sub_alg_df$signal.to.noise.ratio  <- signal_to_noise_ratio(sub_alg, pb = TRUE, parallel = prl, noise.ref = "adjacent", mar = 0.01, output = "data.frame")$signal.to.noise.ratio 

write.csv(sub_alg_df, "./data/processed/degradation_measurements.csv", row.names = FALSE)

sub_alg_df$tail.to.signal.ratio <- tail_to_signal_ratio(sub_alg, pb = TRUE, parallel = prl, mar = 0.01, output = "data.frame")$tail.to.signal.ratio

write.csv(sub_alg_df, "./data/processed/degradation_measurements.csv", row.names = FALSE)

sub_alg_df$cross.correlation <- spcc(X = sub_alg, method = 1, pb = TRUE, parallel = prl, output = "data.frame")$cross.correlation 

write.csv(sub_alg_df, "./data/processed/degradation_measurements.csv", row.names = FALSE)

sub_alg_df$excess.attenuation <- excess_attenuation(X = sub_alg, method = 1, pb = TRUE, parallel = prl, output = "data.frame")$excess.attenuation


write.csv(sub_alg_df, "./data/processed/degradation_measurements.csv", row.names = FALSE)

Statistical analysis

degrad <- read.csv("./data/processed/degradation_measurements.csv")

# force them to be numeric
degrad$blur.ratio <- as.numeric(degrad$blur.ratio)
degrad$spectral.blur.ratio <- as.numeric(degrad$spectral.blur.ratio)
degrad$envelope.correlation <- as.numeric(degrad$envelope.correlation)
degrad$spectral_correlation <- as.numeric(degrad$spectral_correlation)

master_sels <- read.csv("./data/processed/selection_table_open_and_closed_habitat_shuffled.csv")

# remove markers
degrad <- degrad[grep("marker$", degrad$signal.type, invert = TRUE), ]

degrad$clade <- sapply(1:nrow(degrad), function(x){
    
    master_sels$clade[master_sels$orig.sound.file == degrad$template[x]]
    
})

degrad$sp.habitat <- sapply(1:nrow(degrad), function(x){
    
    master_sels$habitat[master_sels$orig.sound.file == degrad$template[x]]
    
})

degrad$species <- sapply(1:nrow(degrad), function(x){
    
    master_sels$species[master_sels$orig.sound.file == degrad$template[x]]
    
})

# remove species not in phylogeny
degrad <- degrad[grep("amirae$", degrad$species, invert = TRUE), ]


degrad$rerec.sound.file <-  sapply(degrad$sound.files, function(x) strsplit(x, "-song")[[1]][1])


transc_metadata <- read.csv("./data/raw/transect_metadata.csv")

degrad$temperature <- sapply(1:nrow(degrad), function(x){
    
    transc_metadata$TEMPERATURE[tolower(transc_metadata$SOUND_FILE_NAME) == tolower(degrad$rerec.sound.file[x])]
    
})

degrad$humidity <- sapply(1:nrow(degrad), function(x){
    
    transc_metadata$HUMIDITY[tolower(transc_metadata$SOUND_FILE_NAME) == tolower(degrad$rerec.sound.file[x])]
    
})

# remove noise transects
degrad <- degrad[!(degrad$transect %in% c("t10", "t11") & degrad$distance == 10), ]

humm_mass <- read_excel("./data/raw/hummingbird_species_body_mass.xlsx")

degrad$size <- sapply(degrad$species, function(x) {
    
    y <- humm_mass$`Body mass (g)`[humm_mass$sp == x]
    if (length(y) == 0) y <- NA
    
    return(y)
    })

write.csv(degrad, "./data/processed/degradation_metrics_and_metadata_hummingbirds_and_frogs.csv", row.names = FALSE)

Prunning trees

degrad <- read.csv("./data/processed/degradation_metrics_and_metadata_hummingbirds_and_frogs.csv")

amph_tree <- read.tree("./data/raw/amph_shl_new_Consensus_7238.tre")

full_hmm_tree <- read.tree("./data/raw/consensus_tree_swifts_hummingbirds_and_nighjars_max_cred.tree")

sp_list_humm <- unique(degrad$species[!grepl("marker", degrad$species) & degrad$clade == "hummingbirds" & !is.na(degrad$clade == "hummingbirds")])

sp_list_frogs <- unique(degrad$species[!grepl("marker", degrad$species) & degrad$clade == "frogs" & !is.na(degrad$clade == "frogs")])

hmm_tree <- drop.tip(full_hmm_tree, tip = setdiff(full_hmm_tree$tip.label, sp_list_humm))

all(sp_list_humm %in% hmm_tree$tip.label)
all(hmm_tree$tip.label %in% sp_list_humm)

write.tree(hmm_tree, "./data/processed/prunned_hmm_tree.tre")

frog_tree <- drop.tip(amph_tree, tip = setdiff(amph_tree$tip.label, sp_list_frogs))

all(sp_list_frogs %in% frog_tree$tip.label)
all(frog_tree$tip.label %in% sp_list_frogs)

write.tree(frog_tree, "./data/processed/prunned_frog_tree.tre")

Plot trees

master_sels <- read.csv("./data/processed/selection_table_open_and_closed_habitat_shuffled.csv")

master_sels <- master_sels[master_sels$in.phlogeny == "yes", ]

# remove markers
master_sels <- master_sels[grep("marker$", invert = TRUE, master_sels$sound.files),]

hab_humm <- as.numeric(as.factor(master_sels$habitat[master_sels$clade == "hummingbirds"]))

names(hab_humm) <- master_sels$species[master_sels$clade == "hummingbirds"]

hab_frogs <- as.numeric(as.factor(master_sels$habitat[master_sels$clade == "frogs"]))

names(hab_frogs) <- master_sels$species[master_sels$clade == "frogs"]

humm_tree <- read.tree("./data/processed/prunned_hmm_tree.tre")

ggtree(humm_tree, layout = "slanted") %<+% master_sels[master_sels$clade == "hummingbirds", c("species", "habitat")] + geom_tippoint(aes(color=habitat), size=5, shape = 20) +
 scale_colour_viridis(discrete = TRUE, begin = 0.3, end = 0.8, alpha = 0.7) + ggtitle("Hummingbirds") + theme_classic(base_size = 30) + labs(color = "Species\n habitat") + theme(axis.ticks = element_blank(), axis.text = element_blank(), axis.line = element_blank())

frog_tree <- read.tree("./data/processed/prunned_frog_tree.tre")

ggtree(frog_tree, layout = "slanted", yscale = "none") %<+% master_sels[master_sels$clade == "frogs", c("species", "habitat")] + geom_tippoint(aes(color=habitat), size=5, shape = 20) +
 scale_colour_viridis(discrete = TRUE, begin = 0.3, end = 0.8, alpha = 0.7) + ggtitle("Frogs") + labs(color = "Species\n habitat") + theme_classic(base_size = 30) + theme(axis.ticks = element_blank(), axis.text = element_blank(), axis.line = element_blank())

Explore data

Hummingbirds

degrad <- read.csv("./data/processed/degradation_metrics_and_metadata_hummingbirds_and_frogs.csv")

humm_degrad <- degrad[degrad$clade == "hummingbirds", ]

humm_degrad$phylo <- humm_degrad$species

preds <- c("blur.ratio", "spectral.blur.ratio", "envelope.correlation", "spectral_correlation", "signal.to.noise.ratio", "tail.to.signal.ratio", "excess.attenuation")

cormat <- cor(humm_degrad[, preds], use = "pairwise.complete.obs")

corrplot.mixed(cormat)

which.complete <- complete.cases(humm_degrad[, preds])

pca <- prcomp(humm_degrad[which.complete, preds],scale. = TRUE)

summary(pca)
## Importance of components:
##                           PC1    PC2    PC3     PC4    PC5     PC6     PC7
## Standard deviation     1.8376 1.0010 0.9922 0.78357 0.6551 0.58907 0.49689
## Proportion of Variance 0.4824 0.1431 0.1406 0.08771 0.0613 0.04957 0.03527
## Cumulative Proportion  0.4824 0.6255 0.7661 0.85386 0.9152 0.96473 1.00000
pca$rotation
##                               PC1         PC2         PC3         PC4
## blur.ratio            -0.42422661  0.07708185 -0.11033943  0.58391475
## spectral.blur.ratio   -0.43871973  0.01576826  0.01605090 -0.56325611
## envelope.correlation   0.44548111 -0.13492961  0.23426212 -0.38769472
## spectral_correlation   0.44614461 -0.03985244  0.09208775  0.40368913
## signal.to.noise.ratio  0.16539564  0.40041103 -0.84841333 -0.16228868
## tail.to.signal.ratio   0.44940815  0.05981115 -0.12607154  0.04528832
## excess.attenuation    -0.02711006 -0.90006065 -0.43417678  0.01119672
##                                PC5         PC6         PC7
## blur.ratio            -0.442567606 -0.39391984 -0.33154189
## spectral.blur.ratio   -0.142409894 -0.54896246  0.41002553
## envelope.correlation  -0.224770836 -0.36571142 -0.62760487
## spectral_correlation   0.370724734 -0.61467362  0.33565044
## signal.to.noise.ratio  0.149197488 -0.13645927 -0.15903293
## tail.to.signal.ratio  -0.757353751  0.11432303  0.43572402
## excess.attenuation    -0.007069975 -0.01442204  0.01621906
humm_degrad$pca1.degradation <- NA

humm_degrad$pca1.degradation[which.complete] <- pca$x[, 1] * -1


preds <- c("blur.ratio", "spectral.blur.ratio", "envelope.correlation", "spectral_correlation", "signal.to.noise.ratio", "tail.to.signal.ratio", "excess.attenuation", "pca1.degradation")

ggs <- lapply(preds, function(x){
    
# plot histogram
ggplot(data = humm_degrad, mapping = aes(x = get(x))) +
        labs(x = x) +
    geom_histogram(fill = cols[10])
    
})


cowplot::plot_grid(plotlist = ggs, nrow = 4)

print("After transforming")
## [1] "After transforming"
ggs <- lapply(preds, function(x){
    
    Y <- humm_degrad
    Y[!is.na(Y[, x]), x] <- log(Y[!is.na(Y[, x]), x] + abs(min(Y[!is.na(Y[, x]), x])) + 1)
    
# plot histogram
ggplot(data = Y, mapping = aes(x = get(x))) +
        labs(x = x) +
    geom_histogram(fill = cols[10])
    
})


cowplot::plot_grid(plotlist = ggs, nrow = 4)

Frogs

degrad <- read.csv("./data/processed/degradation_metrics_and_metadata_hummingbirds_and_frogs.csv")

frog_degrad <- degrad[degrad$clade == "frogs", ]

frog_degrad$phylo <- frog_degrad$species

preds <- c("blur.ratio", "spectral.blur.ratio", "envelope.correlation", "spectral_correlation", "signal.to.noise.ratio", "tail.to.signal.ratio", "excess.attenuation")

cormat <- cor(frog_degrad[, preds], use = "pairwise.complete.obs")

corrplot.mixed(cormat)

which.complete <- complete.cases(frog_degrad[, preds])

pca <- prcomp(frog_degrad[which.complete, preds],scale. = TRUE)

summary(pca)
## Importance of components:
##                           PC1    PC2    PC3    PC4     PC5     PC6     PC7
## Standard deviation     1.5773 1.1715 1.0006 0.8546 0.82804 0.66793 0.52558
## Proportion of Variance 0.3554 0.1961 0.1430 0.1043 0.09795 0.06373 0.03946
## Cumulative Proportion  0.3554 0.5515 0.6945 0.7988 0.89680 0.96054 1.00000
pca$rotation
##                               PC1         PC2          PC3         PC4
## blur.ratio            -0.45812142  0.38623868 -0.007587817 -0.25390189
## spectral.blur.ratio   -0.39324518 -0.29739172 -0.045422275  0.52861424
## envelope.correlation   0.35445934 -0.62272717  0.007517892  0.19404464
## spectral_correlation   0.42073609 -0.07482098  0.030458467 -0.54772403
## signal.to.noise.ratio  0.32265464  0.53032340 -0.022184720  0.53527684
## tail.to.signal.ratio   0.47810182  0.29558897 -0.037764220  0.16975946
## excess.attenuation    -0.01163411  0.01935948  0.997485191  0.05573429
##                              PC5         PC6         PC7
## blur.ratio            0.36220879  0.27672919 -0.60716262
## spectral.blur.ratio   0.61507278  0.09437001  0.29698479
## envelope.correlation  0.05936817  0.11374776 -0.65756862
## spectral_correlation  0.66422584 -0.22462553  0.15748347
## signal.to.noise.ratio 0.19003270 -0.48497383 -0.23733009
## tail.to.signal.ratio  0.09662412  0.78418133  0.17182791
## excess.attenuation    0.01791857  0.03130659  0.01027921
frog_degrad$pca1.degradation <- NA

frog_degrad$pca1.degradation[which.complete] <- pca$x[, 1] * -1


preds <- c("blur.ratio", "spectral.blur.ratio", "envelope.correlation", "spectral_correlation", "signal.to.noise.ratio", "tail.to.signal.ratio", "excess.attenuation", "pca1.degradation")

ggs <- lapply(preds, function(x){
    
# plot histogram
ggplot(data = frog_degrad, mapping = aes(x = get(x))) +
        labs(x = x) +
    geom_histogram(fill = cols[10])
    
})


cowplot::plot_grid(plotlist = ggs, nrow = 4)

print("After transforming")
## [1] "After transforming"
ggs <- lapply(preds, function(x){
    
    Y <- frog_degrad
    Y[!is.na(Y[, x]), x] <- log(Y[!is.na(Y[, x]), x] + abs(min(Y[!is.na(Y[, x]), x])) + 1)
    
# plot histogram
ggplot(data = Y, mapping = aes(x = get(x))) +
        labs(x = x) +
    geom_histogram(fill = cols[10])
    
})


cowplot::plot_grid(plotlist = ggs, nrow = 4)

Stats hummingbirds

humm_tree <- read.tree("./data/processed/prunned_hmm_tree.tre")

all(all(humm_tree$tip.label %in% humm_degrad$species),
all(humm_degrad$species %in% humm_tree$tip.label))

inv_phylo <- ape::vcv.phylo(humm_tree)

itrn <- 3000
chains <- 4

humm_degrad$phylo <- humm_degrad$species
humm_degrad$distance.cat <- as.character(humm_degrad$distance)

preds <- c("blur.ratio", "spectral.blur.ratio", "envelope.correlation", "spectral_correlation", "signal.to.noise.ratio", "tail.to.signal.ratio", "excess.attenuation", "pca1.degradation")


base_models <- c(habitat.only = "~ playback.habitat + (1|gr(phylo, cov = A)) + (1|species) + (1 | transect) + (1 | signal.type)",
            habitat.by.species = "~ playback.habitat * sp.habitat + (1|gr(phylo, cov = A)) + (1|species) + (1 | transect) + (1 | signal.type)", 
            null.model = 
            "~ 1 + (1|gr(phylo, cov = A)) + (1|species) + (1 | transect) + (1 | signal.type)")

# pca.1.models <- gsub("response", "pca1.degradation", models)
models <- paste(rep(preds, each = 3), base_models)

names(models) <- paste(rep(preds, each = 3), names(base_models))

mods <- lapply(models, function(x){

  md <- brm(x,
  data = humm_degrad[humm_degrad$distance == 20, ],
  family = gaussian(),
  data2 = list(A = inv_phylo),
  sample_prior = TRUE, chains = chains, cores = chains,
  iter = itrn
)
  md <- add_criterion(md, "loo")
  
  return(md)
})

names(mods) <- names(models)

saveRDS(mods, "./output/models_habitat_sound_degradation_hummingbirds.RDS")

Results

mods <- readRDS("./output/models_habitat_sound_degradation_hummingbirds.RDS")

for(i in 1:length(mods)){
    print(names(mods)[i])
    summary_brm_model(mods[[i]], plot = FALSE)
  }
[1] “blur.ratio habitat.only”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen -0.02 1.01 1194.041 3000 4 -0.062 0.02
[1] “blur.ratio habitat.by.species”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen -0.021 1.002 2968.173 3000 4 -0.064 0.020
sp.habitatopen 0.006 1.001 1696.861 3000 4 -0.023 0.035
playback.habitatopen:sp.habitatopen 0.001 1.001 7669.907 3000 4 -0.010 0.011
[1] “blur.ratio null.model”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
[1] “spectral.blur.ratio habitat.only”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen -0.073 1.003 1896.638 3000 4 -0.3 0.149
[1] “spectral.blur.ratio habitat.by.species”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen -0.077 1.003 1887.722 3000 4 -0.309 0.160
sp.habitatopen -0.062 1.001 2263.998 3000 4 -0.109 -0.015
playback.habitatopen:sp.habitatopen 0.000 1 10419.003 3000 4 -0.023 0.023
[1] “spectral.blur.ratio null.model”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
[1] “envelope.correlation habitat.only”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen 0.069 1.002 1972.58 3000 4 -0.075 0.221
[1] “envelope.correlation habitat.by.species”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen 0.081 1.001 2245.202 3000 4 -0.073 0.232
sp.habitatopen 0.092 1.001 2037.779 3000 4 0.017 0.165
playback.habitatopen:sp.habitatopen -0.015 1 9733.527 3000 4 -0.048 0.019
[1] “envelope.correlation null.model”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
[1] “spectral_correlation habitat.only”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen 0.097 1.004 1527.642 3000 4 -0.272 0.446
[1] “spectral_correlation habitat.by.species”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen 0.108 1.003 1549.442 3000 4 -0.248 0.458
sp.habitatopen 0.093 1.001 1695.011 3000 4 -0.031 0.218
playback.habitatopen:sp.habitatopen -0.006 1.001 7878.486 3000 4 -0.070 0.060
[1] “spectral_correlation null.model”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
[1] “signal.to.noise.ratio habitat.only”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen -0.138 1.01 1599.121 3000 4 -2.391 2.07
[1] “signal.to.noise.ratio habitat.by.species”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen -0.316 1.001 3371.965 3000 4 -2.579 2.084
sp.habitatopen -0.350 1.001 2568.952 3000 4 -1.733 0.984
playback.habitatopen:sp.habitatopen 0.282 1 11354.265 3000 4 -0.206 0.775
[1] “signal.to.noise.ratio null.model”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
[1] “tail.to.signal.ratio habitat.only”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen 2.296 1 1688.103 3000 4 -0.128 4.847
[1] “tail.to.signal.ratio habitat.by.species”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen 2.360 1.002 1744.868 3000 4 -0.210 5.046
sp.habitatopen 1.317 1.001 1469.278 3000 4 -0.154 2.770
playback.habitatopen:sp.habitatopen 0.021 1.001 9233.236 3000 4 -0.677 0.725
[1] “tail.to.signal.ratio null.model”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
[1] “excess.attenuation habitat.only”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen 1.308 1.002 2499.761 3000 4 -4.075 6.716
[1] “excess.attenuation habitat.by.species”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen 4.330 1 2953.673 3000 4 -2.255 10.964
sp.habitatopen 2.877 1 5558.623 3000 4 -0.386 6.197
playback.habitatopen:sp.habitatopen -4.084 1 4993.579 3000 4 -9.178 1.092
[1] “excess.attenuation null.model”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
[1] “pca1.degradation habitat.only”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen -0.751 1.001 2199.061 3000 4 -2.528 0.992
[1] “pca1.degradation habitat.by.species”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen -0.792 1.001 2271.619 3000 4 -2.560 0.952
sp.habitatopen -0.627 1.001 2379.023 3000 4 -1.162 -0.087
playback.habitatopen:sp.habitatopen 0.038 1.001 14944.611 3000 4 -0.210 0.291
[1] “pca1.degradation null.model”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high

Stats frogs

frog_tree <- read.tree("./data/processed/prunned_frog_tree.tre")

all(all(frog_tree$tip.label %in% frog_degrad$species),
all(frog_degrad$species %in% frog_tree$tip.label))

inv_phylo <- ape::vcv.phylo(frog_tree)

itrn <- 4000
chains <- 4

frog_degrad$phylo <- frog_degrad$species
frog_degrad$distance.cat <- as.character(frog_degrad$distance)

preds <- c("blur.ratio", "spectral.blur.ratio", "envelope.correlation", "spectral_correlation", "signal.to.noise.ratio", "tail.to.signal.ratio", "excess.attenuation", "pca1.degradation")

base_models <- c(habitat.only = "~ playback.habitat + (1|gr(phylo, cov = A)) + (1|species) + (1 | transect) + (1 | signal.type)",
            habitat.by.species = "~ playback.habitat * sp.habitat + (1|gr(phylo, cov = A)) + (1|species) + (1 | transect) + (1 | signal.type)", 
            null.model = 
            "~ 1 + (1|gr(phylo, cov = A)) + (1|species) + (1 | transect) + (1 | signal.type)")

models <- paste(rep(preds, each = 3), base_models)

names(models) <- paste(rep(preds, each = 3), names(base_models))

mods <- lapply(models, function(x){

  md <- brm(x,
  data = frog_degrad[frog_degrad$distance == 20, ],
  family = gaussian(),
  data2 = list(A = inv_phylo),
  sample_prior = TRUE, chains = chains, cores = chains,
  iter = itrn
)
  md <- add_criterion(md, "loo")
  
  return(md)
})

names(mods) <- names(models)

saveRDS(mods, "./output/models_habitat_sound_degradation_frogs.RDS")

Results

mods <- readRDS("./output/models_habitat_sound_degradation_hummingbirds.RDS")

for(i in 1:length(mods)){
    print(names(mods)[i])
    summary_brm_model(mods[[i]], plot = FALSE)
  }
[1] “blur.ratio habitat.only”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen -0.02 1.01 1194.041 3000 4 -0.062 0.02
[1] “blur.ratio habitat.by.species”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen -0.021 1.002 2968.173 3000 4 -0.064 0.020
sp.habitatopen 0.006 1.001 1696.861 3000 4 -0.023 0.035
playback.habitatopen:sp.habitatopen 0.001 1.001 7669.907 3000 4 -0.010 0.011
[1] “blur.ratio null.model”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
[1] “spectral.blur.ratio habitat.only”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen -0.073 1.003 1896.638 3000 4 -0.3 0.149
[1] “spectral.blur.ratio habitat.by.species”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen -0.077 1.003 1887.722 3000 4 -0.309 0.160
sp.habitatopen -0.062 1.001 2263.998 3000 4 -0.109 -0.015
playback.habitatopen:sp.habitatopen 0.000 1 10419.003 3000 4 -0.023 0.023
[1] “spectral.blur.ratio null.model”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
[1] “envelope.correlation habitat.only”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen 0.069 1.002 1972.58 3000 4 -0.075 0.221
[1] “envelope.correlation habitat.by.species”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen 0.081 1.001 2245.202 3000 4 -0.073 0.232
sp.habitatopen 0.092 1.001 2037.779 3000 4 0.017 0.165
playback.habitatopen:sp.habitatopen -0.015 1 9733.527 3000 4 -0.048 0.019
[1] “envelope.correlation null.model”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
[1] “spectral_correlation habitat.only”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen 0.097 1.004 1527.642 3000 4 -0.272 0.446
[1] “spectral_correlation habitat.by.species”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen 0.108 1.003 1549.442 3000 4 -0.248 0.458
sp.habitatopen 0.093 1.001 1695.011 3000 4 -0.031 0.218
playback.habitatopen:sp.habitatopen -0.006 1.001 7878.486 3000 4 -0.070 0.060
[1] “spectral_correlation null.model”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
[1] “signal.to.noise.ratio habitat.only”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen -0.138 1.01 1599.121 3000 4 -2.391 2.07
[1] “signal.to.noise.ratio habitat.by.species”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen -0.316 1.001 3371.965 3000 4 -2.579 2.084
sp.habitatopen -0.350 1.001 2568.952 3000 4 -1.733 0.984
playback.habitatopen:sp.habitatopen 0.282 1 11354.265 3000 4 -0.206 0.775
[1] “signal.to.noise.ratio null.model”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
[1] “tail.to.signal.ratio habitat.only”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen 2.296 1 1688.103 3000 4 -0.128 4.847
[1] “tail.to.signal.ratio habitat.by.species”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen 2.360 1.002 1744.868 3000 4 -0.210 5.046
sp.habitatopen 1.317 1.001 1469.278 3000 4 -0.154 2.770
playback.habitatopen:sp.habitatopen 0.021 1.001 9233.236 3000 4 -0.677 0.725
[1] “tail.to.signal.ratio null.model”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
[1] “excess.attenuation habitat.only”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen 1.308 1.002 2499.761 3000 4 -4.075 6.716
[1] “excess.attenuation habitat.by.species”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen 4.330 1 2953.673 3000 4 -2.255 10.964
sp.habitatopen 2.877 1 5558.623 3000 4 -0.386 6.197
playback.habitatopen:sp.habitatopen -4.084 1 4993.579 3000 4 -9.178 1.092
[1] “excess.attenuation null.model”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
[1] “pca1.degradation habitat.only”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen -0.751 1.001 2199.061 3000 4 -2.528 0.992
[1] “pca1.degradation habitat.by.species”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high
playback.habitatopen -0.792 1.001 2271.619 3000 4 -2.560 0.952
sp.habitatopen -0.627 1.001 2379.023 3000 4 -1.162 -0.087
playback.habitatopen:sp.habitatopen 0.038 1.001 14944.611 3000 4 -0.210 0.291
[1] “pca1.degradation null.model”
Estimate Rhat Bulk_ESS iterations chains CI_low CI_high 
mods_closed <- pblapply(preds, cl = 1, function(i){

  subdat <- degrad[!(is.na(degrad[, i]) | is.infinite(degrad[, i])), ]
  subdat <- subdat[subdat$habitat == "closed", ]
  
  mod <- try(MCMCglmm(as.formula(paste(i, " ~ distance * sp_habitat + size")), random = ~ equipment + transect + signal.type + phylo, data = subdat, verbose = FALSE, ginverse = list(phylo = inv.phylo$Ainv), nitt = itrn), silent = TRUE)

  if(!is(mod, "try-error"))
    return(mod) else
      return(NULL)
})
names(mods_closed) <- preds

saveRDS(mods_closed, "./output/mcmcglmm_models_habitat_by_distance_sound_degradation_in_closed_habitat.RDS")

mods_open <- pblapply(preds, cl = 1, function(i){

  subdat <- degrad[!(is.na(degrad[, i]) | is.infinite(degrad[, i])), ]
  subdat <- subdat[subdat$habitat == "open", ]
  
  mod <- try(MCMCglmm(as.formula(paste(i, " ~ distance * sp_habitat + size")), random = ~ equipment + transect + signal.type + phylo, data = subdat, verbose = FALSE, ginverse = list(phylo = inv.phylo$Ainv), nitt = itrn), silent = TRUE)

  if(!is(mod, "try-error"))
    return(mod) else
      return(NULL)
})
names(mods_open) <- preds

saveRDS(mods_open, "./output/mcmcglmm_models_habitat_by_distance_sound_degradation_in_open_habitat.RDS")
wv <- read_wave("./data/raw/recordings/ETC21_duplicated_master_44.1.wav")

library(dynaSpec)
scrolling_spectro(wave = wv, wl = 400,
    t.display = 1, ovlp = 95, pal = inferno,
    grid = FALSE, flim = c(1, 13), width = 900,
    height = 250, res = 100, collevels = seq(-40, 0, 5),
    file.name = "./output/playback_full_osc.mp4", colbg = "black", lcol = "#FFFFFFE6",
    speed = 0.7, fps = 200, buffer = 0, loop = 1, lty = 1,
    osc = TRUE, colwave = inferno(10, alpha = 0.9)[3])

scrolling_spectro(wave = wv, wl = 400, lower.spectro = TRUE,
    t.display = 1, ovlp = 95, pal = inferno,
    grid = FALSE, flim = c(1, 13), width = 900,
    height = 400, res = 100, collevels = seq(-40, 0, 5),
    file.name = "./output/playback_full.mp4", colbg = "black", lcol = "#FFFFFFE6",
    speed = 0.7, fps = 200, buffer = 0, loop = 1, lty = 1,
    osc = FALSE, colwave = inferno(10, alpha = 0.9)[3])

scrolling_spectro(wave = wv, wl = 400, lower.spectro = TRUE,axis.type = "minimal",
    t.display = 1, ovlp = 95, pal = inferno,
    grid = FALSE, flim = c(1, 13), width = 900,
    height = 400, res = 100, collevels = seq(-40, 0, 5),
    file.name = "./output/playback_full_minimal.mp4", colbg = "black", lcol = "#FFFFFFE6",
    speed = 0.7, fps = 200, buffer = 0, loop = 1, lty = 1,
    osc = FALSE, colwave = inferno(10, alpha = 0.9)[3])

scrolling_spectro(wave = wv, wl = 400, lower.spectro = TRUE,axis.type = "none",
    t.display = 1, ovlp = 95, pal = inferno,
    grid = FALSE, flim = c(1, 13), width = 900,
    height = 400, res = 100, collevels = seq(-40, 0, 5),
    file.name = "./output/playback_full_no_axis.mp4", colbg = "black", lcol = "#FFFFFFE6",
    speed = 0.7, fps = 200, buffer = 0, loop = 1, lty = 1,
    osc = FALSE, colwave = inferno(10, alpha = 0.9)[3])

Descriptive stats

# master_sels <- read.csv("./data/processed/selection_table_open_and_closed_habitat_shuffled.csv")

degrad <- read.csv("./data/processed/degradation_metrics_and_metadata_hummingbirds_and_frogs.csv")


# master_sels <- master_sels[master_sels$in.phlogeny == "yes", ]

# # remove markers
# master_sels <- master_sels[grep("marker$", invert = TRUE, master_sels$sound.files),]

# number of especies per clade in open habitat
print("number of especies per clade in open habitat")
## [1] "number of especies per clade in open habitat"
table(degrad$sp.habitat[!duplicated(degrad$species)], degrad$clade[!duplicated(degrad$species)])
##         
##          frogs hummingbirds
##   closed    20           14
##   open      20           35
# mean number of sounds per especies
print("mean number of sounds per especies")
## [1] "mean number of sounds per especies"
mean(table(degrad$species[degrad$transect == "t1" & degrad$distance == 1]))
## [1] 2.505618
print("range of number of sounds per especies")
## [1] "range of number of sounds per especies"
range(table(degrad$species[degrad$transect == "t1" & degrad$distance == 1]))
## [1] 1 5
# number of signals
print("number of signals (sounds)")
## [1] "number of signals (sounds)"
length(unique(degrad$signal.type))
## [1] 223
print("number of signals per clade")
## [1] "number of signals per clade"
aggregate(signal.type ~ clade, data = degrad, FUN = function(x) length(unique(x)))
clade signal.type
frogs 130
hummingbirds 93 
print("total number of signals")
## [1] "total number of signals"
nrow(degrad)
## [1] 10481
print("total number of signals per clade")
## [1] "total number of signals per clade"
tapply(1:nrow(degrad), degrad$clade, length)
##        frogs hummingbirds 
##         6110         4371
print("species per habitat")
## [1] "species per habitat"
sub_degrad <- degrad[!duplicated(degrad$species), c("clade", "species", "sp.habitat")]

sub_degrad[order(sub_degrad$clade, sub_degrad$sp.habitat), ]
clade species sp.habitat
frogs Craugastor_melanostictus closed
frogs Aplastodiscus_cochranae closed
frogs Proceratophrys_concavitympanum closed
frogs Dendropsophus_ebraccatus closed
frogs Diasporus_ventrimaculatus closed
frogs Isthmohyla_zeteki closed
frogs Craugastor_crassidigitus closed
frogs Odontophrynus_carvalhoi closed
frogs Hypsiboas_rufitelus closed
frogs Craugastor_talamancae closed
frogs Diasporus_hylaeformis closed
frogs Diasporus_diastema closed
frogs Craugastor_fitzingeri closed
frogs Anotheca_spinosa closed
frogs Pristimantis_cruentus closed
frogs Syncope_bassleri closed
frogs Oophaga_pumilio closed
frogs Pristimantis_altae closed
frogs Dendrobates_auratus closed
frogs Allobates_talamancae closed
frogs Leptodactylus_fragilis open
frogs Scinax_x-signatus open
frogs Scinax_boulengeri open
frogs Hypsiboas_raniceps open
frogs Leptodactylus_savagei open
frogs Phyllomedusa_nordestina open
frogs Leptodactylus_insularum open
frogs Hypsiboas_albopunctatus open
frogs Scinax_staufferi open
frogs Smilisca_puma open
frogs Smilisca_phaeota open
frogs Dendropsophus_phlebodes open
frogs Hypsiboas_rosenbergi open
frogs Hypsiboas_albomarginatus open
frogs Engystomops_pustulosus open
frogs Isthmohyla_pseudopuma open
frogs Scinax_elaeochroa open
frogs Incilius_coniferus open
frogs Leptodactylus_poecilochilus open
frogs Incilius_luetkenii open
hummingbirds Calypte_anna closed
hummingbirds Oreotrochilus_melanogaster closed
hummingbirds Amazilia_yucatanensis closed
hummingbirds Phaethornis_koepckeae closed
hummingbirds Phaethornis_malaris closed
hummingbirds Glaucis_hirsutus closed
hummingbirds Sappho_sparganurus closed
hummingbirds Polytmus_guainumbi closed
hummingbirds Metallura_phoebe closed
hummingbirds Eupetomena_macroura closed
hummingbirds Phaethornis_bourcieri closed
hummingbirds Glaucis_dohrnii closed
hummingbirds Threnetes_leucurus closed
hummingbirds Oxypogon_stubelii closed
hummingbirds Amazilia_cyanifrons open
hummingbirds Campylopterus_villaviscensio open
hummingbirds Lepidopyga_goudoti open
hummingbirds Oreotrochilus_estella open
hummingbirds Calypte_costae open
hummingbirds Lepidopyga_lilliae open
hummingbirds Eriocnemis_cupreoventris open
hummingbirds Anthracothorax_veraguensis open
hummingbirds Chlorostilbon_alice open
hummingbirds Chalcostigma_stanleyi open
hummingbirds Metallura_eupogon open
hummingbirds Heliodoxa_gularis open
hummingbirds Anthracothorax_nigricollis open
hummingbirds Myrtis_fanny open
hummingbirds Ramphodon_naevius open
hummingbirds Oxypogon_guerinii open
hummingbirds Phaethornis_squalidus open
hummingbirds Metallura_theresiae open
hummingbirds Doryfera_johannae open
hummingbirds Heliactin_bilophus open
hummingbirds Coeligena_wilsoni open
hummingbirds Phaethornis_aethopygus open
hummingbirds Selasphorus_flammula open
hummingbirds Oreotrochilus_adela open
hummingbirds Heliomaster_furcifer open
hummingbirds Phaethornis_hispidus open
hummingbirds Heliodoxa_schreibersii open
hummingbirds Anthocephala_floriceps open
hummingbirds Amazilia_versicolor open
hummingbirds Metallura_williami open
hummingbirds Phaethornis_longuemareus open
hummingbirds Chrysolampis_mosquitus open
hummingbirds Chlorestes_notata open
hummingbirds Androdon_aequatorialis open
hummingbirds Aglaeactis_cupripennis open

Graphs

Hummingbirds

theme_set(theme_classic(base_size = 24))

humm_degrad <- humm_degrad[humm_degrad$distance != 40, ]

agg_deg_humm <- aggregate(cbind(blur.ratio, spectral.blur.ratio, envelope.correlation, spectral_correlation, signal.to.noise.ratio, tail.to.signal.ratio, excess.attenuation, pca1.degradation) ~ playback.habitat + sp.habitat, humm_degrad, FUN = se)

se.humm <- aggregate(cbind(blur.ratio, spectral.blur.ratio, envelope.correlation, spectral_correlation, signal.to.noise.ratio, tail.to.signal.ratio, excess.attenuation, pca1.degradation) ~ playback.habitat + sp.habitat, humm_degrad, FUN = se)

names(se.humm) <- paste0(names(se.humm), ".se")

agg_deg_humm <- data.frame(agg_deg_humm, se.humm)

agg_deg_humm$se <- aggregate(pca1.degradation ~ playback.habitat + sp.habitat, humm_degrad, FUN = se)$pca1.degradation

preds <- c("blur.ratio", "spectral.blur.ratio", "envelope.correlation", "spectral_correlation", "signal.to.noise.ratio", "tail.to.signal.ratio", "excess.attenuation", "pca1.degradation")

ggs <- lapply(preds, function(x){
    
    x.se <- paste(x, "se", sep = ".")
    
    ggplot(agg_deg_humm, aes(x = playback.habitat, y = get(x), color = sp.habitat)) + geom_pointrange(aes(ymin = get(x) - get(x.se), ymax = get(x) + get(x.se)), position = position_dodge(width = 0.5), size = 2) + labs(x = "Playback habitat", y = x, color = "Species\n habitat") + scale_color_viridis_d(begin = 0.3, end = 0.8, alpha = 0.7)

})

cowplot::plot_grid(plotlist = ggs, ncol = 2)

agg_deg_humm <- aggregate(cbind(blur.ratio, spectral.blur.ratio, envelope.correlation, spectral_correlation, signal.to.noise.ratio, tail.to.signal.ratio, excess.attenuation, pca1.degradation) ~ playback.habitat + sp.habitat + distance, humm_degrad, FUN = se)

se.humm <- aggregate(cbind(blur.ratio, spectral.blur.ratio, envelope.correlation, spectral_correlation, signal.to.noise.ratio, tail.to.signal.ratio, excess.attenuation, pca1.degradation) ~ playback.habitat + sp.habitat + distance, humm_degrad, FUN = se)

names(se.humm) <- paste0(names(se.humm), ".se")

agg_deg_humm <- data.frame(agg_deg_humm, se.humm)

# agg_deg_humm$se <- aggregate(pca1.degradation ~ playback.habitat + sp.habitat, humm_degrad, FUN = se)$pca1.degradation

preds <- c("blur.ratio", "spectral.blur.ratio", "envelope.correlation", "spectral_correlation", "signal.to.noise.ratio", "tail.to.signal.ratio", "excess.attenuation", "pca1.degradation")

ggs <- lapply(preds, function(x){
    
    x.se <- paste(x, "se", sep = ".")
    
    ggplot(agg_deg_humm, aes(x = distance, y = get(x), color = sp.habitat, shape = playback.habitat)) + geom_pointrange(aes(ymin = get(x) - get(x.se), ymax = get(x) + get(x.se)), position = position_dodge(width = 1.5), size = 1.2) + geom_line(position = position_dodge(width = 1.5)) +
        labs(x = "Distance", y = x, color = "Species\n habitat", shape = "Playback\n habitat") + scale_color_viridis_d(begin = 0.3, end = 0.8, alpha = 0.7)

})

cowplot::plot_grid(plotlist = ggs, ncol = 1)

Frogs

frog_degrad <- frog_degrad[frog_degrad$distance != 40, ]

agg_deg_frog <- aggregate(cbind(blur.ratio, spectral.blur.ratio, envelope.correlation, spectral_correlation, signal.to.noise.ratio, tail.to.signal.ratio, excess.attenuation, pca1.degradation) ~ playback.habitat + sp.habitat, frog_degrad, FUN = se)

se.frog <- aggregate(cbind(blur.ratio, spectral.blur.ratio, envelope.correlation, spectral_correlation, signal.to.noise.ratio, tail.to.signal.ratio, excess.attenuation, pca1.degradation) ~ playback.habitat + sp.habitat, frog_degrad, FUN = se)

names(se.frog) <- paste0(names(se.frog), ".se")

agg_deg_frog <- data.frame(agg_deg_frog, se.frog)

agg_deg_frog$se <- aggregate(pca1.degradation ~ playback.habitat + sp.habitat, frog_degrad, FUN = se)$pca1.degradation

preds <- c("blur.ratio", "spectral.blur.ratio", "envelope.correlation", "spectral_correlation", "signal.to.noise.ratio", "tail.to.signal.ratio", "excess.attenuation", "pca1.degradation")

ggs <- lapply(preds, function(x){
    
    x.se <- paste(x, "se", sep = ".")
    
    ggplot(agg_deg_frog, aes(x = playback.habitat, y = get(x), color = sp.habitat)) + geom_pointrange(aes(ymin = get(x) - get(x.se), ymax = get(x) + get(x.se)), position = position_dodge(width = 0.5), size = 2) + labs(x = "Playback habitat", y = x, color = "Species\n habitat") + scale_color_viridis_d(begin = 0.3, end = 0.8, alpha = 0.7)

})

cowplot::plot_grid(plotlist = ggs, ncol = 2)

agg_deg_frog <- aggregate(cbind(blur.ratio, spectral.blur.ratio, envelope.correlation, spectral_correlation, signal.to.noise.ratio, tail.to.signal.ratio, excess.attenuation, pca1.degradation) ~ playback.habitat + sp.habitat + distance, frog_degrad, FUN = se)

se.frog <- aggregate(cbind(blur.ratio, spectral.blur.ratio, envelope.correlation, spectral_correlation, signal.to.noise.ratio, tail.to.signal.ratio, excess.attenuation, pca1.degradation) ~ playback.habitat + sp.habitat + distance, frog_degrad, FUN = se)

names(se.frog) <- paste0(names(se.frog), ".se")

agg_deg_frog <- data.frame(agg_deg_frog, se.frog)

# agg_deg_frog$se <- aggregate(pca1.degradation ~ playback.habitat + sp.habitat, frog_degrad, FUN = se)$pca1.degradation

preds <- c("blur.ratio", "spectral.blur.ratio", "envelope.correlation", "spectral_correlation", "signal.to.noise.ratio", "tail.to.signal.ratio", "excess.attenuation", "pca1.degradation")

ggs <- lapply(preds, function(x){
    
    x.se <- paste(x, "se", sep = ".")
    
    ggplot(agg_deg_frog, aes(x = distance, y = get(x), color = sp.habitat, shape = playback.habitat)) + geom_pointrange(aes(ymin = get(x) - get(x.se), ymax = get(x) + get(x.se)), position = position_dodge(width = 1.5), size = 1.2) + geom_line(position = position_dodge(width = 1.5)) +
        labs(x = "Distance", y = x, color = "Species\n habitat", shape = "Playback\n habitat") + scale_color_viridis_d(begin = 0.3, end = 0.8, alpha = 0.7)

})

cowplot::plot_grid(plotlist = ggs, ncol = 1)

Prediction graphs

  • No acoustic adaptation
theme_set(theme_classic(base_size = 30))

# set seed
set.seed(13)

# number of observations
n <- 300
b0 <- 4
b1 <- 0
b2 <- 0
error <- rnorm(n = n, mean = 0, sd = 1)

# random variables
x1_num <- sample(0:1, size = n, replace = TRUE)
x2_num <- sample(0:1, size = n, replace = TRUE)

y <- b0 + b1 * x1_num + b2* x2_num + error

sp_hab <- factor(x1_num, labels = c("Open", "Closed"))

playb_hab <- factor(x2_num, labels = c("Open", "Closed"))

# create data frame
xy_data_cat <- data.frame(sp_hab, x1_num, playb_hab, x2_num, y)

agg_xy <- aggregate(y ~ sp_hab + playb_hab, xy_data_cat, FUN = mean)

agg_xy$sd <- aggregate(y ~ sp_hab + playb_hab, xy_data_cat, FUN = sd)$y

ggplot(agg_xy, aes(x = playb_hab, y = y, color = sp_hab)) + geom_pointrange(aes(ymin = y - sd, ymax = y + sd), position = position_dodge(width = 0.5), size = 2) + labs(x = "Playback habitat", y = "-  Degradation  +", color = "Species\n habitat") + scale_color_viridis_d(begin = 0.3, end = 0.8, alpha = 0.7) + theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank())

  • Only habitat matters
# set seed
set.seed(13)

# number of observations
n <- 300
b0 <- 4
b1 <- 0
b2 <- 2
error <- rnorm(n = n, mean = 0, sd = 1)

# random variables
x1_num <- sample(0:1, size = n, replace = TRUE)
x2_num <- sample(0:1, size = n, replace = TRUE)

y <- b0 + b1 * x1_num + b2* x2_num + error

sp_hab <- factor(x1_num, labels = c("Open", "Closed"))

playb_hab <- factor(x2_num, labels = c("Open", "Closed"))

# create data frame
xy_data_cat <- data.frame(sp_hab, x1_num, playb_hab, x2_num, y)

agg_xy <- aggregate(y ~ sp_hab + playb_hab, xy_data_cat, FUN = mean)

agg_xy$sd <- aggregate(y ~ sp_hab + playb_hab, xy_data_cat, FUN = sd)$y

ggplot(agg_xy, aes(x = playb_hab, y = y, color = sp_hab)) + geom_pointrange(aes(ymin = y - sd, ymax = y + sd), position = position_dodge(width = 0.5), size = 2) + labs(x = "Playback habitat", y = "-  Degradation  +", color = "Species\n habitat") + scale_color_viridis_d(begin = 0.3, end = 0.8, alpha = 0.7) + theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank())

  • Adversity
# set seed
set.seed(13)

# number of observations
n <- 300
b0 <- 4
b1 <- -2
b2 <- 2
error <- rnorm(n = n, mean = 0, sd = 1)

# random variables
x1_num <- sample(0:1, size = n, replace = TRUE)
x2_num <- sample(0:1, size = n, replace = TRUE)

y <- b0 + b1 * x1_num + b2* x2_num + error

sp_hab <- factor(x1_num, labels = c("Open", "Closed"))

playb_hab <- factor(x2_num, labels = c("Open", "Closed"))

# create data frame
xy_data_cat <- data.frame(sp_hab, x1_num, playb_hab, x2_num, y)

agg_xy <- aggregate(y ~ sp_hab + playb_hab, xy_data_cat, FUN = mean)

agg_xy$sd <- aggregate(y ~ sp_hab + playb_hab, xy_data_cat, FUN = sd)$y

ggplot(agg_xy, aes(x = playb_hab, y = y, color = sp_hab)) + geom_pointrange(aes(ymin = y - sd, ymax = y + sd), position = position_dodge(width = 0.5), size = 2) + labs(x = "Playback habitat", y = "-  Degradation  +", color = "Species\n habitat") + scale_color_viridis_d(begin = 0.3, end = 0.8, alpha = 0.7) + theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank())

  • Home-team
# set seed
set.seed(13)

# number of observations
n <- 300
b0 <- 4
b1 <- 2
b2 <- 3
b3 <- -4
error <- rnorm(n = n, mean = 0, sd = 1)

# random variables
x1_num <- sample(0:1, size = n, replace = TRUE)
x2_num <- sample(0:1, size = n, replace = TRUE)

y <- b0 + b1 * x1_num + b2 * x2_num + b3 * x1_num * x2_num + error

sp_hab <- factor(x1_num, labels = c("Open", "Closed"))

playb_hab <- factor(x2_num, labels = c("Open", "Closed"))

# create data frame
xy_data_cat <- data.frame(sp_hab, x1_num, playb_hab, x2_num, y)

agg_xy <- aggregate(y ~ sp_hab + playb_hab, xy_data_cat, FUN = mean)

agg_xy$sd <- aggregate(y ~ sp_hab + playb_hab, xy_data_cat, FUN = sd)$y

ggplot(agg_xy, aes(x = playb_hab, y = y, color = sp_hab)) + geom_pointrange(aes(ymin = y - sd, ymax = y + sd), position = position_dodge(width = 0.5), size = 2) + labs(x = "Playback habitat", y = "-  Degradation  +", color = "Species\n habitat") + scale_color_viridis_d(begin = 0.3, end = 0.8, alpha = 0.7) + theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank())