Effect of background noise on schroeder structure detection

Fine scale acoustic perception in zebra finches

Author

Marcelo Araya-Salas

Published

April 8, 2025

Source code and data found at https://github.com/maRce10/acoustic-fine-features-zebra-finch

 

1 Purpose

  • Evaluate the effect of degradation on the detection of acoustic fine structural variation in Schroeders

 

Code 
source("~/Dropbox/Projects/geographic_call_variation_yellow-naped_amazon/scripts/MRM2.R")

warbleR_options(parallel = 22)

2 Data analysis workflow

flowchart LR
  A[Raw Schroeder\nrecordings] --> B(Aritfically\nDecrease SNR) 
  B --> C(measure\nSchroeder\nsimilarity)
  C --> D(Matrix\nregression\nmodels)
  
style A fill:#44015466
style B fill:#3E4A894D
style C fill:#6DCD594D
style D fill:#FDE7254D

3 Degradation as a function of signal-to-noise ratio

  • Adding synthetic noise to the master sound file annotations to decrease the signal-to-noise
  • Evaluate at which signal-to-noise ratio Schroeder sign can still be distinguished

3.1 On 200ms schroeders

3.1.1 Add synthetic noise

3.1.1.1 Repeated Schroeders

Code 
master_annotations <- imp_raven(path = "./data/processed", files = "200ms_schroeders.txt",
    warbler.format = TRUE, all.data = TRUE)

master_annotations$sound.id <- paste0("f0:", master_annotations$f0,
    "_comp:", master_annotations$label)

master_annotations$sound.id[1] <- "start_marker"
master_annotations$sound.id[nrow(master_annotations)] <- "end_marker"

est_master <- selection_table(master_annotations[grep("marker", master_annotations$sound.id,
    invert = TRUE), ], extended = TRUE, mar = 0.1, path = "./data/processed")


# est_schr <-
# readRDS('./data/processed/extended_selection_table_schroeders.RDS')
# est_schr <- readRDS(file.path(path,
# 'extended_sel_table_degradation_exp.RDS')) est_schr$treatment
# <- ifelse(grepl('inside', est_schr$sound.files), 'inside',
# 'outside') est_schr$sound.id <- 1:nrow(est_schr) est_in <-
# est_schr[est_schr$treatment == 'inside', ]

est_master <- signal_to_noise_ratio(est_master, mar = 0.07, cores = 20)

hist(est_master$signal.to.noise.ratio)

adj_snr_schroeder_white <- lapply(30:1, function(i) {
    print(i)
    est_schr_adj <- baRulho::add_noise(X = est_master, target.snr = i,
        precision = 0.1, mar = 0.07, cores = 20, kind = "white", seed = NULL)
    est_schr_adj$target.snr <- i
    return(est_schr_adj)
})

beepr::beep(2)

saveRDS(adj_snr_schroeder_white, "./data/processed/200ms_schroeders_adjusted_snr_white_noise.RDS")

3.1.1.2 Single Schroders

Code 
est_schr <- readRDS("./data/processed/extended_selection_table_single_schroeders.RDS")

est_schr$sound.id <- "a"
est_schr <- signal_to_noise_ratio(X = est_schr, mar = 0.01, cores = 20,
    bp = NULL)

hist(est_schr$signal.to.noise.ratio)

adj_snr_schroeder_white <- lapply(30:1, function(i) {
    print(i)
    est_schr_adj <- add_noise(X = est_schr, target.snr = i, precision = 0.1,
        mar = 0.01, cores = 1, bp = NULL, kind = "white", seed = NULL)
    est_schr_adj$target.snr <- i
    return(est_schr_adj)
})

beepr::beep(2)

saveRDS(adj_snr_schroeder_white, "./data/processed/single_schroeder_adjusted_snr_white_noise.RDS")

3.1.2 Check SNR adjusted waves

3.1.2.1 Examples

3.1.2.1.1 Repeated Schroeders

200 Hz, 9 components, positive phase

Code 
adj_snr_schroeder <- readRDS("./data/processed/200ms_schroeders_adjusted_snr_white_noise.RDS")

# print(adj_snr_schroeder[[1]]$sound.id[10])

par(mar = c(0, 0, 0, 0), mfrow = c(5, 6))
for (i in 1:30) {
    warbleR:::spectro_wrblr_int(wave = read_sound_file(adj_snr_schroeder[[i]],
        10, from = 0, to = Inf), flim = c(0, 5), osc = FALSE, scale = FALSE,
        axisY = FALSE, axisX = FALSE, grid = FALSE, palette = viridis::mako,
        collevels = seq(-110, 0, 5))
    text(0.2, 4, paste("SNR:", adj_snr_schroeder[[i]]$target.snr[1]),
        col = "black", cex = 1.3)
}

500 Hz, 16 components, positive phase

Code 
# print(adj_snr_schroeder[[1]]$sound.id[150])

par(mar = c(0, 0, 0, 0), mfrow = c(5, 6))
for (i in 1:30) {
    warbleR:::spectro_wrblr_int(wave = read_sound_file(adj_snr_schroeder[[i]],
        150, from = 0, to = Inf), flim = c(0, 12), osc = FALSE, scale = FALSE,
        axisY = FALSE, axisX = FALSE, grid = FALSE, palette = viridis::mako,
        collevels = seq(-110, 0, 5))
    text(0.2, 11, paste("SNR:", adj_snr_schroeder[[i]]$target.snr[1]),
        col = "black", cex = 1.3)
}

700 Hz, 24 components, positive phase

Code 
# print(adj_snr_schroeder[[1]]$sound.id[250])

par(mar = c(0, 0, 0, 0), mfrow = c(5, 6))
for (i in 1:30) {
    warbleR:::spectro_wrblr_int(wave = read_sound_file(adj_snr_schroeder[[i]],
        250, from = 0, to = Inf), flim = c(0, 20), osc = FALSE, scale = FALSE,
        axisY = FALSE, axisX = FALSE, grid = FALSE, palette = viridis::mako,
        collevels = seq(-110, 0, 5))
    text(0.2, 18, paste("SNR:", adj_snr_schroeder[[i]]$target.snr[1]),
        col = "black", cex = 1.3)
}

3.1.2.1.2 Single Schroeders

200 Hz, 9 components, positive phase

Code 
adj_snr_schroeder <- readRDS("./data/processed/single_schroeder_adjusted_snr_white_noise.RDS")

# print(adj_snr_schroeder[[1]]$sound.id[10])

par(mar = c(0, 0, 0, 0), mfrow = c(5, 6))
for (i in 1:30) {

    wv <- read_sound_file(adj_snr_schroeder[[i]], 10)@left

    plot(wv, type = "l", lwd = 2, col = mako(10)[6], xlab = "", ylab = "",
        axes = FALSE, ylim = c(-2, 2))
    box()

    text(length(wv)/2, 1.7, paste("SNR:", adj_snr_schroeder[[i]]$target.snr[1]),
        col = "black", cex = 1.3)
}

200 Hz, 9 components, negative phase

Code 
# print(adj_snr_schroeder[[1]]$sound.id[9])

par(mar = c(0, 0, 0, 0), mfrow = c(5, 6))
for (i in 1:30) {

    wv <- read_sound_file(adj_snr_schroeder[[i]], 9)@left

    plot(wv, type = "l", lwd = 2, col = mako(10)[6], xlab = "", ylab = "",
        axes = FALSE, ylim = c(-2, 2))
    box()

    text(length(wv)/2, 1.7, paste("SNR:", adj_snr_schroeder[[i]]$target.snr[1]),
        col = "black", cex = 1.3)
}

500 Hz, 16 components, positive phase

Code 
# print(adj_snr_schroeder[[1]]$sound.files[150])

par(mar = c(0, 0, 0, 0), mfrow = c(5, 6))
for (i in 1:30) {

    wv <- read_sound_file(adj_snr_schroeder[[i]], 150)@left

    plot(wv, type = "l", lwd = 2, col = mako(10)[6], xlab = "", ylab = "",
        axes = FALSE, ylim = c(-2, 2))
    box()

    text(length(wv)/2, 1.7, paste("SNR:", adj_snr_schroeder[[i]]$target.snr[1]),
        col = "black", cex = 1.3)
}

500 Hz, 16 components, negative phase

Code 
# print(adj_snr_schroeder[[1]]$sound.files[149])

par(mar = c(0, 0, 0, 0), mfrow = c(5, 6))
for (i in 1:30) {

    wv <- read_sound_file(adj_snr_schroeder[[i]], 149)@left

    plot(wv, type = "l", lwd = 2, col = mako(10)[6], xlab = "", ylab = "",
        axes = FALSE, ylim = c(-2, 2))
    box()

    text(length(wv)/2, 1.7, paste("SNR:", adj_snr_schroeder[[i]]$target.snr[1]),
        col = "black", cex = 1.3)
}

700 Hz, 24 components, positive phase

Code 
# print(adj_snr_schroeder[[1]]$sound.files[250])

par(mar = c(0, 0, 0, 0), mfrow = c(5, 6))
for (i in 1:30) {

    wv <- read_sound_file(adj_snr_schroeder[[i]], 250)@left

    plot(wv, type = "l", lwd = 2, col = mako(10)[6], xlab = "", ylab = "",
        axes = FALSE, ylim = c(-2, 2))
    box()

    text(length(wv)/2, 1.7, paste("SNR:", adj_snr_schroeder[[i]]$target.snr[1]),
        col = "black", cex = 1.3)
}

700 Hz, 24 components, negative phase

Code 
# print(adj_snr_schroeder[[1]]$sound.files[249])

par(mar = c(0, 0, 0, 0), mfrow = c(5, 6))
for (i in 1:30) {

    wv <- read_sound_file(adj_snr_schroeder[[i]], 249)@left

    plot(wv, type = "l", lwd = 2, col = mako(10)[6], xlab = "", ylab = "",
        axes = FALSE, ylim = c(-2, 2))
    box()

    text(length(wv)/2, 1.7, paste("SNR:", adj_snr_schroeder[[i]]$target.snr[1]),
        col = "black", cex = 1.3)
}

3.1.3 Measure DTW distances

3.1.3.1 On repeated Schroeders

Code 
adj_snr_schroeder <- readRDS("./data/processed/200ms_schroeders_adjusted_snr_white_noise.RDS")

dtw_dist_white_noise <- lapply(adj_snr_schroeder, function(schroeders) {

    print(schroeders$target.snr[1])
    cmbs <- t(combn(schroeders$sound.files, 2))

    min_dist_l <- pbapply::pbsapply(cl = 22, 1:nrow(cmbs), function(x) {
        s1 <- read_sound_file(schroeders, index = which(schroeders$sound.files ==
            cmbs[x, 1]))@left
        s2 <- read_sound_file(schroeders, index = which(schroeders$sound.files ==
            cmbs[x, 2]))@left

        # normalize
        s1 <- s1/max(abs(s1))
        s2 <- s2/max(abs(s2))

        # make same length if (length(s1) != length(s2))
        s1 <- approx(s1, n = 100)$y
        s2 <- approx(s2, n = 100)$y

        # duplicate 1
        s1 <- rep(s1, 2)

        # run dtw over longer vector
        dists <- vapply(seq_len(length(s1) - length(s2)), function(x) {
            segment <- s1[x:min(c(x + length(s2) - 1), length(s1))]

            dtw_dist <- dtw::dtwDist(mx = rbind(s2, segment))
            return(dtw_dist[1, 2])
        }, FUN.VALUE = numeric(1))

        return(data.frame(schr1 = cmbs[x, 1], schr2 = cmbs[x, 2],
            min(dists)))
    })

    min_dists <- do.call(rbind, min_dist_l)

    min_dists <- as.data.frame(matrix(min_dists[, 1], ncol = 3, byrow = TRUE))

    names(min_dists) <- c("schr1", "schr2", "dtw.dist")

    min_dists$dtw.dist <- as.numeric(min_dists$dtw.dist)

    min_dists$target.snr <- schroeders$target.snr[1]

    return(min_dists)

})

beepr::beep(2)

saveRDS(dtw_dist_white_noise, "./data/processed/dtw_distance_adjusted_snr_200ms_schroeders_white_noise.RDS")

3.1.3.2 On single Schroeders

Code 
adj_snr_schroeder <- readRDS("./data/processed/single_schroeder_adjusted_snr_white_noise.RDS")

dtw_dist_white_noise <- lapply(adj_snr_schroeder, function(schroeders) {

    print(schroeders$target.snr[1])
    cmbs <- t(combn(schroeders$sound.files, 2))

    min_dist_l <- pbapply::pbsapply(cl = 22, 1:nrow(cmbs), function(x) {
        s1 <- read_sound_file(schroeders, index = which(schroeders$sound.files ==
            cmbs[x, 1]))@left
        s2 <- read_sound_file(schroeders, index = which(schroeders$sound.files ==
            cmbs[x, 2]))@left

        # normalize
        s1 <- s1/max(abs(s1))
        s2 <- s2/max(abs(s2))

        # make same length if (length(s1) != length(s2))
        s1 <- approx(s1, n = 100)$y
        s2 <- approx(s2, n = 100)$y

        # duplicate 1
        s1 <- rep(s1, 2)

        # run dtw over longer vector
        dists <- vapply(seq_len(length(s1) - length(s2)), function(x) {
            segment <- s1[x:min(c(x + length(s2) - 1), length(s1))]

            dtw_dist <- dtw::dtwDist(mx = rbind(s2, segment))
            return(dtw_dist[1, 2])
        }, FUN.VALUE = numeric(1))

        return(data.frame(schr1 = cmbs[x, 1], schr2 = cmbs[x, 2],
            min(dists)))
    })

    min_dists <- do.call(rbind, min_dist_l)

    min_dists <- as.data.frame(matrix(min_dists[, 1], ncol = 3, byrow = TRUE))

    names(min_dists) <- c("schr1", "schr2", "dtw.dist")

    min_dists$dtw.dist <- as.numeric(min_dists$dtw.dist)

    min_dists$target.snr <- schroeders$target.snr[1]

    return(min_dists)

})

beepr::beep(2)

saveRDS(dtw_dist_white_noise, "./data/processed/dtw_distance_adjusted_snr_single_schroeder_white_noise.RDS")

3.1.3.3 Statistical analysis

Modeling: - Multiple Regression on distance Matrices - Model:
\[\begin{align*} Dissimilarity &\sim frequency + components + phase \end{align*}\] - Response values scaled to make effect sizes comparable across models - Predictors were coded as pairwise binary matrices in which 0 means that calls in a dyad belong to the same level and 1 means calls belong to different levels - One model for each environment treatment as well as on the original model sounds

3.1.3.3.1 On repeated Schroeders
Code 
dtw_dists_snr_list <- readRDS("./data/processed/dtw_distance_adjusted_snr_200ms_schroeders_white_noise.RDS")

adj_snr_schroeder1 <- readRDS("./data/processed/200ms_schroeders_adjusted_snr_white_noise.RDS")[[1]]

adj_snr_schroeder1$sound.id <- gsub("f0:f0 ", "f0:", adj_snr_schroeder1$sound.id)


dtw_wv_mod_list <- warbleR:::pblapply_wrblr_int(dtw_dists_snr_list,
    cl = 20, function(x) {

        # x$dist <- scale(x$dist)
        target_snr <- x$target.snr
        x$target.snr <- NULL

        x$schr1 <- sapply(x$schr1, function(x) adj_snr_schroeder1$sound.id[adj_snr_schroeder1$sound.files ==
            x])
        x$schr2 <- sapply(x$schr2, function(x) adj_snr_schroeder1$sound.id[adj_snr_schroeder1$sound.files ==
            x])

        dist_tri <- PhenotypeSpace::rectangular_to_triangular(x)

        freq_bi_tri <- as.dist(binary_triangular_matrix(group = sapply(strsplit(rownames(dist_tri),
            "_"), "[[", 1)))

        comp_bi_tri <- as.dist(binary_triangular_matrix(group = gsub("_n|_p",
            "", sapply(strsplit(rownames(dist_tri), "_"), "[[", 2))))

        sign_bi_tri <- as.dist(binary_triangular_matrix(group = sapply(strsplit(rownames(dist_tri),
            "_"), "[[", 3)))

        rect_var <- cbind(as.dist(dist_tri), freq_bi_tri, comp_bi_tri,
            sign_bi_tri)

        colnames(rect_var) <- c("dtw_dist", "frequency", "components",
            "sign")

        dtw_wv_mod <- MRM2(formula = dtw_dist ~ frequency + components +
            sign, nperm = 10000, mat = rect_var)

        return(dtw_wv_mod)
    })

names(dtw_wv_mod_list) <- paste0("snr:", sapply(dtw_dists_snr_list,
    function(x) x$target.snr[1]))

saveRDS(dtw_wv_mod_list, "./data/processed/matrix_correlation_dtw_distance_varying_snr_white_noise.RDS")
Code 
mods <- readRDS("./data/processed/matrix_correlation_dtw_distance_varying_snr_white_noise.RDS")


estimates_dtw_wn <- do.call(rbind, lapply(seq_along(mods), function(x) {
    Y <- data.frame(mods[[x]]$coef[-1, ])
    # Y$rel_coef <- Y[, 1] / max(Y[, 1])
    Y$rel_coef <- Y[, 1]
    Y$predictor <- rownames(Y)
    Y$model <- as.numeric(gsub("snr:", "", names(mods)[x]))
    names(Y) <- c("coef", "p", "rel_coef", "predictor", "model")
    return(Y)
}))

estimates_dtw_wn$rel_coef <- ifelse(estimates_dtw_wn$p < 0.05, estimates_dtw_wn$rel_coef,
    0)
estimates_dtw_wn$signif <- ifelse(estimates_dtw_wn$p < 0.05, "p < 0.05",
    "p >= 0.05")
estimates_dtw_wn$signif2 <- ifelse(estimates_dtw_wn$p < 0.05, "*",
    "")

estimates_dtw_wn$model_f <- factor(estimates_dtw_wn$model, levels = 1:30)

ggplot(estimates_dtw_wn, aes(x = predictor, y = model_f, fill = rel_coef)) +
    geom_tile() + coord_equal() + scale_fill_gradient2(low = viridis(10)[3],
    high = viridis(10)[7], guide = "none") + geom_text(aes(label = signif2,
    color = signif), size = 3) + scale_color_manual(values = c("black",
    "gray")) + labs(y = "Signal-to-noise ratio", x = "Predictor",
    color = "P value") + theme_classic() + theme(axis.text.x = element_text(color = "black",
    size = 7, angle = 30, vjust = 0.8, hjust = 0.8)) + coord_flip()

Effect sizes by signal-to-noise ratio and predictor (color intensity shows effect size magnitude relative to the highest effect size within the model).
Code 
estimates_dtw_wn$model <- as.numeric(estimates_dtw_wn$model)
Code 
coef_by_signif <- aggregate(coef ~ signif, data = estimates_dtw_wn,
    FUN = range)

ggplot(estimates_dtw_wn, aes(y = coef, x = model, color = predictor,
    lty = predictor)) + geom_line(linewidth = 1.6) + labs(x = "Signal-to-noise ratio",
    y = "Effect size") + scale_color_viridis_d(alpha = 0.7, begin = 0.2,
    end = 0.8) + geom_hline(yintercept = mean(coef_by_signif[, 2][1,
    1], coef_by_signif[, 2][2, 2]), lty = 2, color = "gray") + theme_classic()

Effect sizes by signal-to-noise ratio and predictor.
3.1.3.3.1.1 Only phase
Code 
ggplot(estimates_dtw_wn[estimates_dtw_wn$predictor == "sign", ], aes(y = coef,
    x = model, color = predictor, lty = predictor)) + geom_line(linewidth = 1.6) +
    geom_smooth() + labs(x = "Signal-to-noise ratio", y = "Effect size") +
    scale_color_viridis_d(alpha = 0.7, begin = 0.2, end = 0.8) + geom_hline(yintercept = mean(coef_by_signif[,
    2][1, 1], coef_by_signif[, 2][2, 2]), lty = 2, color = "gray") +
    theme_classic() + geom_vline(xintercept = 16, lty = 2, color = "gray")

Effect sizes by signal-to-noise ratio and predictor.
3.1.3.3.1.2 By frequency
Code 
dtw_dists_snr_list <- readRDS("./data/processed/dtw_distance_adjusted_snr_200ms_schroeders_white_noise.RDS")

adj_snr_schroeder1 <- readRDS("./data/processed/200ms_schroeders_adjusted_snr_pink_noise.RDS")[[1]]

adj_snr_schroeder1$sound.id <- gsub("f0:f0 ", "f0:", adj_snr_schroeder1$sound.id)

dtw_wv_mod_list <- warbleR:::pblapply_wrblr_int(dtw_dists_snr_list,
    cl = 1, function(y) {

        # print(y$target.snr[1])
        freqs <- unique(substr(sapply(y[, 1], function(x) adj_snr_schroeder1$sound.id[adj_snr_schroeder1$sound.files ==
            x]), 0, 6))

        mods <- lapply(1:length(freqs), function(i) {

            y$schr1 <- sapply(y$schr1, function(x) adj_snr_schroeder1$sound.id[adj_snr_schroeder1$sound.files ==
                x])
            y$schr2 <- sapply(y$schr2, function(x) adj_snr_schroeder1$sound.id[adj_snr_schroeder1$sound.files ==
                x])

            # extract only those with the same freq
            x <- y[grepl(paste(freqs[i], collapse = "|"), y$schr1) &
                grepl(paste(freqs[i], collapse = "|"), y$schr2), ]

            target_snr <- x$target.snr
            x$target.snr <- NULL

            dist_tri <- PhenotypeSpace::rectangular_to_triangular(x)

            if (nrow(dist_tri) > 2) {
                comp_bi_tri <- as.dist(binary_triangular_matrix(group = gsub("_n|_p",
                  "", sapply(strsplit(rownames(dist_tri), "_"), "[[",
                    2))))

                sign_bi_tri <- as.dist(binary_triangular_matrix(group = sapply(strsplit(rownames(dist_tri),
                  "_"), "[[", 3)))

                rect_var <- cbind(as.dist(dist_tri), comp_bi_tri,
                  sign_bi_tri)

                colnames(rect_var) <- c("dtw_dist", "components",
                  "sign")

                dtw_wv_mod <- MRM2(formula = dtw_dist ~ components +
                  sign, nperm = 10000, mat = rect_var)
            } else dtw_wv_mod <- NULL
            return(dtw_wv_mod)
        })

        names(mods) <- paste0("snr:", y$target.snr[1], "_", freqs)

        mods <- mods[!sapply(mods, is.null)]

        return(mods)
    })

saveRDS(dtw_wv_mod_list, "./data/processed/matrix_correlation_dtw_distance_varying_snr_by_frequency_white_noise.RDS")

Discriminating Schroeder’s by phase grouped by frequency:

Code 
mods <- readRDS("./data/processed/matrix_correlation_dtw_distance_varying_snr_by_frequency_white_noise.RDS")

estimates_list <- lapply(mods, function(y) {
    estimates <- do.call(rbind, lapply(seq_along(y), function(x) {
        Y <- data.frame(y[[x]]$coef[-1, ])
        # Y$rel_coef <- Y[, 1] / max(Y[, 1])
        Y$rel_coef <- Y[, 1]
        Y$mod <- names(y)[x]
        Y$predictor <- rownames(Y)
        Y$freq <- strsplit(names(y)[x], "_")[[1]][2]
        Y$freq <- strsplit(Y$freq, ":")[[1]][2]
        names(Y) <- c("coef", "p", "rel_coef", "model", "predictor",
            "freq")
        Y$snr <- gsub("_", "", substr(names(y)[x], 0, 6))

        return(Y)
    }))

    estimates$rel_coef <- ifelse(estimates$p < 0.05, estimates$rel_coef,
        0)
    estimates$signif <- ifelse(estimates$p < 0.05, "p < 0.05", "p >= 0.05")

    return(estimates)
})
estimates <- do.call(rbind, estimates_list)

estimates$snr <- as.numeric(substr(estimates$snr, 5, nchar(estimates$snr)))

rownames(estimates) <- 1:nrow(estimates)

coef_by_signif <- aggregate(coef ~ signif, data = estimates[estimates$predictor ==
    "sign", ], FUN = range)

ggplot(estimates[estimates$predictor == "sign", ], aes(y = coef, x = snr,
    color = freq)) + geom_line(linewidth = 1.6) + ylim(c(-0.2, 0.5)) +
    labs(x = "Signal-to-noise ratio", y = "Effect size when discriminating phase",
        color = "F0\nfrequency") + scale_color_viridis_d(alpha = 0.7,
    begin = 0.2, end = 0.8) + theme_classic()
Warning: Removed 24 rows containing missing values or values outside the scale range
(`geom_line()`).

Effect sizes for phase by signal-to-noise ratio and frequency.

Discriminating Schroeder’s by number of components grouped by frequency:

Code 
coef_by_signif <- aggregate(coef ~ signif, data = estimates[estimates$predictor ==
    "components", ], FUN = range)


ggplot(estimates[estimates$predictor == "components", ], aes(y = coef,
    x = snr, color = freq)) + geom_line(linewidth = 1.6) + labs(x = "Signal-to-noise ratio",
    y = "Effect size when discriminating number of components", color = "F0\nfrequency") +
    scale_color_viridis_d(alpha = 0.7, begin = 0.1, end = 0.9) + geom_hline(yintercept = mean(coef_by_signif[,
    2][1, 1], coef_by_signif[, 2][2, 2]), lty = 2, color = "gray") +
    theme_classic()

Effect sizes for number of components by signal-to-noise ratio and frequency.
3.1.3.3.2 On single Schroeders
Code 
dtw_dists_snr_list <- readRDS("./data/processed/dtw_distance_adjusted_snr_single_schroeder_white_noise.RDS")

adj_snr_schroeder1 <- readRDS("./data/processed/200ms_schroeders_adjusted_snr_white_noise.RDS")[[1]]

adj_snr_schroeder1$sound.id <- gsub("f0:f0 ", "f0:", adj_snr_schroeder1$sound.id)


dtw_wv_mod_list <- warbleR:::pblapply_wrblr_int(dtw_dists_snr_list,
    cl = 20, function(x) {

        # x$dist <- scale(x$dist)
        target_snr <- x$target.snr
        x$target.snr <- NULL

        x$schr1 <- sapply(strsplit(x$schr1, "_"), function(x) paste(x[2],
            x[3], substr(x[4], 0, 1), sep = "_"))
        x$schr2 <- sapply(strsplit(x$schr2, "_"), function(x) paste(x[2],
            x[3], substr(x[4], 0, 1), sep = "_"))

        # x$schr1 <- sapply(x$schr1, function(x)
        # adj_snr_schroeder1$sound.id[adj_snr_schroeder1$sound.files
        # == x]) x$schr2 <- sapply(x$schr2, function(x)
        # adj_snr_schroeder1$sound.id[adj_snr_schroeder1$sound.files
        # == x])
        dist_tri <- PhenotypeSpace::rectangular_to_triangular(x)

        freq_bi_tri <- as.dist(binary_triangular_matrix(group = sapply(strsplit(rownames(dist_tri),
            "_"), "[[", 1)))

        comp_bi_tri <- as.dist(binary_triangular_matrix(group = gsub("_n|_p",
            "", sapply(strsplit(rownames(dist_tri), "_"), "[[", 2))))

        sign_bi_tri <- as.dist(binary_triangular_matrix(group = sapply(strsplit(rownames(dist_tri),
            "_"), "[[", 3)))

        rect_var <- cbind(as.dist(dist_tri), freq_bi_tri, comp_bi_tri,
            sign_bi_tri)

        colnames(rect_var) <- c("dtw_dist", "frequency", "components",
            "sign")

        dtw_wv_mod <- MRM2(formula = dtw_dist ~ frequency + components +
            sign, nperm = 10000, mat = rect_var)

        return(dtw_wv_mod)
    })

names(dtw_wv_mod_list) <- paste0("snr:", sapply(dtw_dists_snr_list,
    function(x) x$target.snr[1]))

saveRDS(dtw_wv_mod_list, "./data/processed/matrix_correlation_dtw_distance_single_schroeder_varying_snr_white_noise.RDS")
Code 
mods <- readRDS("./data/processed/matrix_correlation_dtw_distance_single_schroeder_varying_snr_white_noise.RDS")


estimates_dtw_wn <- do.call(rbind, lapply(seq_along(mods), function(x) {
    Y <- data.frame(mods[[x]]$coef[-1, ])
    # Y$rel_coef <- Y[, 1] / max(Y[, 1])
    Y$rel_coef <- Y[, 1]
    Y$predictor <- rownames(Y)
    Y$model <- as.numeric(gsub("snr:", "", names(mods)[x]))
    names(Y) <- c("coef", "p", "rel_coef", "predictor", "model")
    return(Y)
}))

estimates_dtw_wn$rel_coef <- ifelse(estimates_dtw_wn$p < 0.05, estimates_dtw_wn$rel_coef,
    0)
estimates_dtw_wn$signif <- ifelse(estimates_dtw_wn$p < 0.05, "p < 0.05",
    "p >= 0.05")
estimates_dtw_wn$signif2 <- ifelse(estimates_dtw_wn$p < 0.05, "*",
    "")

estimates_dtw_wn$model_f <- factor(estimates_dtw_wn$model, levels = 1:30)

ggplot(estimates_dtw_wn, aes(x = predictor, y = model_f, fill = rel_coef)) +
    geom_tile() + coord_equal() + scale_fill_gradient2(low = viridis(10)[3],
    high = viridis(10)[7], guide = "none") + geom_text(aes(label = signif2,
    color = signif), size = 3) + scale_color_manual(values = c("black",
    "gray")) + labs(y = "Signal-to-noise ratio", x = "Predictor",
    color = "P value") + theme_classic() + theme(axis.text.x = element_text(color = "black",
    size = 7, angle = 30, vjust = 0.8, hjust = 0.8)) + coord_flip()

Effect sizes by signal-to-noise ratio and predictor (color intensity shows effect size magnitude relative to the highest effect size within the model).
Code 
estimates_dtw_wn$model <- as.numeric(estimates_dtw_wn$model)
Code 
coef_by_signif <- aggregate(coef ~ signif, data = estimates_dtw_wn,
    FUN = range)

ggplot(estimates_dtw_wn, aes(y = coef, x = model, color = predictor,
    lty = predictor)) + geom_line(linewidth = 1.6) + labs(x = "Signal-to-noise ratio",
    y = "Effect size") + scale_color_viridis_d(alpha = 0.7, begin = 0.2,
    end = 0.8) + geom_hline(yintercept = mean(coef_by_signif[, 2][1,
    1], coef_by_signif[, 2][2, 2]), lty = 2, color = "gray") + theme_classic()

Effect sizes by signal-to-noise ratio and predictor.
3.1.3.3.2.1 Only phase
Code 
ggplot(estimates_dtw_wn[estimates_dtw_wn$predictor == "sign", ], aes(y = coef,
    x = model, color = predictor, lty = predictor)) + geom_line(linewidth = 1.6) +
    geom_smooth() + labs(x = "Signal-to-noise ratio", y = "Effect size") +
    scale_color_viridis_d(alpha = 0.7, begin = 0.2, end = 0.8) + geom_hline(yintercept = mean(coef_by_signif[,
    2][1, 1], coef_by_signif[, 2][2, 2]), lty = 2, color = "gray") +
    theme_classic() + geom_vline(xintercept = 16, lty = 2, color = "gray")

Effect sizes by signal-to-noise ratio and predictor.
3.1.3.3.2.2 By frequency
Code 
dtw_dists_snr_list <- readRDS("./data/processed/dtw_distance_adjusted_snr_200ms_schroeders_white_noise.RDS")

adj_snr_schroeder1 <- readRDS("./data/processed/200ms_schroeders_adjusted_snr_pink_noise.RDS")[[1]]

adj_snr_schroeder1$sound.id <- gsub("f0:f0 ", "f0:", adj_snr_schroeder1$sound.id)

dtw_wv_mod_list <- warbleR:::pblapply_wrblr_int(dtw_dists_snr_list,
    cl = 1, function(y) {

        # print(y$target.snr[1])
        freqs <- unique(substr(sapply(y[, 1], function(x) adj_snr_schroeder1$sound.id[adj_snr_schroeder1$sound.files ==
            x]), 0, 6))

        mods <- lapply(1:length(freqs), function(i) {

            y$schr1 <- sapply(y$schr1, function(x) adj_snr_schroeder1$sound.id[adj_snr_schroeder1$sound.files ==
                x])
            y$schr2 <- sapply(y$schr2, function(x) adj_snr_schroeder1$sound.id[adj_snr_schroeder1$sound.files ==
                x])

            # extract only those with the same freq
            x <- y[grepl(paste(freqs[i], collapse = "|"), y$schr1) &
                grepl(paste(freqs[i], collapse = "|"), y$schr2), ]

            target_snr <- x$target.snr
            x$target.snr <- NULL

            dist_tri <- PhenotypeSpace::rectangular_to_triangular(x)

            if (nrow(dist_tri) > 2) {
                comp_bi_tri <- as.dist(binary_triangular_matrix(group = gsub("_n|_p",
                  "", sapply(strsplit(rownames(dist_tri), "_"), "[[",
                    2))))

                sign_bi_tri <- as.dist(binary_triangular_matrix(group = sapply(strsplit(rownames(dist_tri),
                  "_"), "[[", 3)))

                rect_var <- cbind(as.dist(dist_tri), comp_bi_tri,
                  sign_bi_tri)

                colnames(rect_var) <- c("dtw_dist", "components",
                  "sign")

                dtw_wv_mod <- MRM2(formula = dtw_dist ~ components +
                  sign, nperm = 10000, mat = rect_var)
            } else dtw_wv_mod <- NULL
            return(dtw_wv_mod)
        })

        names(mods) <- paste0("snr:", y$target.snr[1], "_", freqs)

        mods <- mods[!sapply(mods, is.null)]

        return(mods)
    })

saveRDS(dtw_wv_mod_list, "./data/processed/matrix_correlation_dtw_distance_varying_snr_by_frequency_white_noise.RDS")

Discriminating Schroeder’s by phase grouped by frequency:

Code 
mods <- readRDS("./data/processed/matrix_correlation_dtw_distance_varying_snr_by_frequency_white_noise.RDS")

estimates_list <- lapply(mods, function(y) {
    estimates <- do.call(rbind, lapply(seq_along(y), function(x) {
        Y <- data.frame(y[[x]]$coef[-1, ])
        # Y$rel_coef <- Y[, 1] / max(Y[, 1])
        Y$rel_coef <- Y[, 1]
        Y$mod <- names(y)[x]
        Y$predictor <- rownames(Y)
        Y$freq <- strsplit(names(y)[x], "_")[[1]][2]
        Y$freq <- strsplit(Y$freq, ":")[[1]][2]
        names(Y) <- c("coef", "p", "rel_coef", "model", "predictor",
            "freq")
        Y$snr <- gsub("_", "", substr(names(y)[x], 0, 6))

        return(Y)
    }))

    estimates$rel_coef <- ifelse(estimates$p < 0.05, estimates$rel_coef,
        0)
    estimates$signif <- ifelse(estimates$p < 0.05, "p < 0.05", "p >= 0.05")

    return(estimates)
})
estimates <- do.call(rbind, estimates_list)

estimates$snr <- as.numeric(substr(estimates$snr, 5, nchar(estimates$snr)))

rownames(estimates) <- 1:nrow(estimates)

coef_by_signif <- aggregate(coef ~ signif, data = estimates[estimates$predictor ==
    "sign", ], FUN = range)

ggplot(estimates[estimates$predictor == "sign", ], aes(y = coef, x = snr,
    color = freq)) + geom_line(linewidth = 1.6) + ylim(c(-0.2, 0.5)) +
    labs(x = "Signal-to-noise ratio", y = "Effect size when discriminating phase",
        color = "F0\nfrequency") + scale_color_viridis_d(alpha = 0.7,
    begin = 0.2, end = 0.8) + theme_classic()
Warning: Removed 24 rows containing missing values or values outside the scale range
(`geom_line()`).

Effect sizes for phase by signal-to-noise ratio and frequency.

Discriminating Schroeder’s by number of components grouped by frequency:

Code 
coef_by_signif <- aggregate(coef ~ signif, data = estimates[estimates$predictor ==
    "components", ], FUN = range)


ggplot(estimates[estimates$predictor == "components", ], aes(y = coef,
    x = snr, color = freq)) + geom_line(linewidth = 1.6) + labs(x = "Signal-to-noise ratio",
    y = "Effect size when discriminating number of components", color = "F0\nfrequency") +
    scale_color_viridis_d(alpha = 0.7, begin = 0.1, end = 0.9) + geom_hline(yintercept = mean(coef_by_signif[,
    2][1, 1], coef_by_signif[, 2][2, 2]), lty = 2, color = "gray") +
    theme_classic()

Effect sizes for number of components by signal-to-noise ratio and frequency.

3.1.4 Measure waveform correlation

3.1.4.1 On repeated Schroeders

Code 
adj_snr_schroeder <- readRDS("./data/processed/200ms_schroeders_adjusted_snr_white_noise.RDS")

cor_noise <- lapply(adj_snr_schroeder, function(schroeders) {

    print(schroeders$target.snr[1])
    cmbs <- t(combn(schroeders$sound.files, 2))

    min_cor_l <- pbapply::pbsapply(cl = 20, 1:nrow(cmbs), function(x) {
        s1 <- read_sound_file(schroeders, index = which(schroeders$sound.files ==
            cmbs[x, 1]))@left
        s2 <- read_sound_file(schroeders, index = which(schroeders$sound.files ==
            cmbs[x, 2]))@left

        # normalize
        s1 <- s1/max(abs(s1))
        s2 <- s2/max(abs(s2))

        # make same length if (length(s1) != length(s2))
        s1 <- approx(s1, n = 100)$y
        s2 <- approx(s2, n = 100)$y

        # run correlation
        cor <- cor(s1, s2)

        return(data.frame(schr1 = cmbs[x, 1], schr2 = cmbs[x, 2],
            cor = cor))

    })

    cors <- do.call(rbind, min_cor_l)

    cors <- as.data.frame(matrix(cors[, 1], ncol = 3, byrow = TRUE))

    names(cors) <- c("schr1", "schr2", "cor")

    cors$cor <- as.numeric(cors$cor)

    cors$target.snr <- cors$target.snr[1]

    return(cors)

})

names(cor_noise) <- paste0("snr:", sapply(adj_snr_schroeder, function(x) x$target.snr[1]))


beepr::beep(2)

saveRDS(cor_noise, "./data/processed/waveform_correlation_adjusted_snr_200ms_schroeders_white_noise.RDS")

3.1.4.2 On single Schroeders

Code 
adj_snr_schroeder <- readRDS("./data/processed/single_schroeder_adjusted_snr_white_noise.RDS")

cor_noise <- lapply(adj_snr_schroeder, function(schroeders) {

    print(schroeders$target.snr[1])
    cmbs <- t(combn(schroeders$sound.files, 2))

    min_cor_l <- pbapply::pbsapply(cl = 20, 1:nrow(cmbs), function(x) {
        s1 <- read_sound_file(schroeders, index = which(schroeders$sound.files ==
            cmbs[x, 1]))@left
        s2 <- read_sound_file(schroeders, index = which(schroeders$sound.files ==
            cmbs[x, 2]))@left

        # normalize
        s1 <- s1/max(abs(s1))
        s2 <- s2/max(abs(s2))

        # make same length if (length(s1) != length(s2))
        s1 <- approx(s1, n = 100)$y
        s2 <- approx(s2, n = 100)$y

        # run correlation
        cor <- cor(s1, s2)

        return(data.frame(schr1 = cmbs[x, 1], schr2 = cmbs[x, 2],
            cor = cor))

    })

    cors <- do.call(rbind, min_cor_l)

    cors <- as.data.frame(matrix(cors[, 1], ncol = 3, byrow = TRUE))

    names(cors) <- c("schr1", "schr2", "cor")

    cors$cor <- as.numeric(cors$cor)

    cors$target.snr <- cors$target.snr[1]

    return(cors)

})

names(cor_noise) <- paste0("snr:", sapply(adj_snr_schroeder, function(x) x$target.snr[1]))


beepr::beep(2)

saveRDS(cor_noise, "./data/processed/waveform_correlation_adjusted_snr_single_schroeders_white_noise.RDS")

3.1.4.3 Statistical analysis

3.1.4.3.1 On repeated Schroeders
Code 
cor_snr_list <- readRDS("./data/processed/waveform_correlation_adjusted_snr_200ms_schroeders_white_noise.RDS")

adj_snr_schroeder1 <- readRDS("./data/processed/200ms_schroeders_adjusted_snr_white_noise.RDS")[[1]]

adj_snr_schroeder1$sound.id <- gsub("f0:f0 ", "f0:", adj_snr_schroeder1$sound.id)


cor_mod_list <- warbleR:::pblapply_wrblr_int(names(cor_snr_list),
    cl = 20, function(x) {

        x <- cor_snr_list[[x]]
        # x$dist <- scale(x$dist)
        target_snr <- "snr:30"

        x$schr1 <- sapply(x$schr1, function(x) adj_snr_schroeder1$sound.id[adj_snr_schroeder1$sound.files ==
            x])
        x$schr2 <- sapply(x$schr2, function(x) adj_snr_schroeder1$sound.id[adj_snr_schroeder1$sound.files ==
            x])
        x$cor <- x$cor * -1

        dist_tri <- PhenotypeSpace::rectangular_to_triangular(x)

        freq_bi_tri <- as.dist(binary_triangular_matrix(group = sapply(strsplit(rownames(dist_tri),
            "_"), "[[", 1)))

        comp_bi_tri <- as.dist(binary_triangular_matrix(group = gsub("_n|_p",
            "", sapply(strsplit(rownames(dist_tri), "_"), "[[", 2))))

        sign_bi_tri <- as.dist(binary_triangular_matrix(group = sapply(strsplit(rownames(dist_tri),
            "_"), "[[", 3)))

        rect_var <- cbind(as.dist(dist_tri), freq_bi_tri, comp_bi_tri,
            sign_bi_tri)

        colnames(rect_var) <- c("cor", "frequency", "components",
            "sign")

        cor_wv_mod <- MRM2(formula = cor ~ frequency + components +
            sign, nperm = 10000, mat = rect_var)

        return(cor_wv_mod)
    })

names(cor_mod_list) <- names(cor_snr_list)

saveRDS(cor_mod_list, "./data/processed/matrix_regression_waveform_correlation_varying_snr_white_noise.RDS")
Code 
mods <- readRDS("./data/processed/matrix_regression_waveform_correlation_varying_snr_white_noise.RDS")

estimates_wf_wn <- do.call(rbind, lapply(seq_along(mods), function(x) {
    Y <- data.frame(mods[[x]]$coef[-1, ])
    # Y$rel_coef <- Y[, 1] / max(Y[, 1])
    Y$rel_coef <- Y[, 1]
    Y$predictor <- rownames(Y)
    Y$model <- as.numeric(gsub("snr:", "", names(mods)[x]))
    names(Y) <- c("coef", "p", "rel_coef", "predictor", "model")
    return(Y)
}))

estimates_wf_wn$rel_coef <- ifelse(estimates_wf_wn$p < 0.05, estimates_wf_wn$rel_coef,
    0)
estimates_wf_wn$signif <- ifelse(estimates_wf_wn$p < 0.05, "p < 0.05",
    "p >= 0.05")
estimates_wf_wn$signif2 <- ifelse(estimates_wf_wn$p < 0.05, "*", "")

estimates_wf_wn$model_f <- factor(estimates_wf_wn$model, levels = 1:30)

ggplot(estimates_wf_wn, aes(x = predictor, y = model_f, fill = rel_coef)) +
    geom_tile() + coord_equal() + scale_fill_gradient2(low = viridis(10)[3],
    high = viridis(10)[7], guide = "none") + geom_text(aes(label = signif2,
    color = signif), size = 3) + scale_color_manual(values = c("black",
    "gray")) + labs(y = "Signal-to-noise ratio", x = "Predictor",
    color = "P value") + theme_classic() + theme(axis.text.x = element_text(color = "black",
    size = 7, angle = 30, vjust = 0.8, hjust = 0.8)) + coord_flip()

Effect sizes by signal-to-noise ratio and predictor (color intensity shows effect size magnitude relative to the highest effect size within the model).
Code 
estimates_wf_wn$model <- as.numeric(estimates_wf_wn$model)
Code 
coef_by_signif <- aggregate(coef ~ signif, data = estimates_wf_wn,
    FUN = range)

ggplot(estimates_wf_wn, aes(y = coef, x = model, color = predictor,
    lty = predictor)) + geom_line(linewidth = 1.6) + labs(x = "Signal-to-noise ratio",
    y = "Effect size") + scale_color_viridis_d(alpha = 0.7, begin = 0.2,
    end = 0.8) + geom_hline(yintercept = mean(coef_by_signif[, 2][1,
    1], coef_by_signif[, 2][2, 2]), lty = 2, color = "gray") + theme_classic()

Effect sizes by signal-to-noise ratio and predictor waveform correlation
3.1.4.3.1.1 Only phase
Code 
ggplot(estimates_wf_wn[estimates_wf_wn$predictor == "sign", ], aes(y = coef,
    x = model, color = predictor, lty = predictor)) + geom_line(linewidth = 1.6) +
    geom_smooth() + labs(x = "Signal-to-noise ratio", y = "Effect size") +
    scale_color_viridis_d(alpha = 0.7, begin = 0.2, end = 0.8) + geom_hline(yintercept = mean(coef_by_signif[,
    2][1, 1], coef_by_signif[, 2][2, 2]), lty = 2, color = "gray") +
    theme_classic() + geom_vline(xintercept = 2, lty = 2, color = "gray")

Effect sizes by signal-to-noise ratio and predictor waveform correlation.
3.1.4.3.1.2 By frequency
Code 
cor_snr_list <- readRDS("./data/processed/waveform_correlation_adjusted_snr_200ms_schroeders_white_noise.RDS")
adj_snr_schroeder1 <- readRDS("./data/processed/200ms_schroeders_adjusted_snr_pink_noise.RDS")[[1]]


cor_snr_list <- lapply(seq_along(cor_snr_list), function(x) {
    Y <- cor_snr_list[[x]]
    Y$target.snr <- gsub("snr:", "", names(cor_snr_list)[x])
    return(Y)
})

adj_snr_schroeder1$sound.id <- gsub("f0:f0 ", "f0:", adj_snr_schroeder1$sound.id)

cor_mod_list <- warbleR:::pblapply_wrblr_int(cor_snr_list, cl = 1,
    function(y) {

        freqs <- unique(substr(sapply(y[, 1], function(x) adj_snr_schroeder1$sound.id[adj_snr_schroeder1$sound.files ==
            x]), 0, 6))

        mods <- lapply(1:length(freqs), function(i) {

            y$schr1 <- sapply(y$schr1, function(x) adj_snr_schroeder1$sound.id[adj_snr_schroeder1$sound.files ==
                x])
            y$schr2 <- sapply(y$schr2, function(x) adj_snr_schroeder1$sound.id[adj_snr_schroeder1$sound.files ==
                x])

            # extract only those with the same freq
            x <- y[grepl(paste(freqs[i], collapse = "|"), y$schr1) &
                grepl(paste(freqs[i], collapse = "|"), y$schr2), ]

            target_snr <- x$target.snr
            x$target.snr <- NULL
            x$cor <- 1 - x$cor

            dist_tri <- PhenotypeSpace::rectangular_to_triangular(x)

            if (nrow(dist_tri) > 2) {
                comp_bi_tri <- as.dist(binary_triangular_matrix(group = gsub("_n|_p",
                  "", sapply(strsplit(rownames(dist_tri), "_"), "[[",
                    2))))

                sign_bi_tri <- as.dist(binary_triangular_matrix(group = sapply(strsplit(rownames(dist_tri),
                  "_"), "[[", 3)))

                rect_var <- cbind(as.dist(dist_tri), comp_bi_tri,
                  sign_bi_tri)

                colnames(rect_var) <- c("dtw_dist", "components",
                  "sign")

                cor_mod <- MRM2(formula = dtw_dist ~ components +
                  sign, nperm = 10000, mat = rect_var)
            } else cor_mod <- NULL
            return(cor_mod)
        })

        names(mods) <- paste0("snr:", y$target.snr[1], "_", freqs)

        mods <- mods[!sapply(mods, is.null)]

        return(mods)
    })

saveRDS(cor_mod_list, "./data/processed/matrix_correlation_waveforme_correlatione_varying_snr_by_frequency_white_noise.RDS")

Discriminating Schroeder’s by phase grouped by frequency:

Code 
mods <- readRDS("./data/processed/matrix_correlation_waveforme_correlatione_varying_snr_by_frequency_white_noise.RDS")

estimates_list <- lapply(mods, function(y) {
    estimates <- do.call(rbind, lapply(seq_along(y), function(x) {
        Y <- data.frame(y[[x]]$coef[-1, ])
        # Y$rel_coef <- Y[, 1] / max(Y[, 1])
        Y$rel_coef <- Y[, 1]
        Y$mod <- names(y)[x]
        Y$predictor <- rownames(Y)
        Y$freq <- strsplit(names(y)[x], "_")[[1]][2]
        Y$freq <- strsplit(Y$freq, ":")[[1]][2]
        names(Y) <- c("coef", "p", "rel_coef", "model", "predictor",
            "freq")
        Y$snr <- gsub("_", "", substr(names(y)[x], 0, 6))

        return(Y)
    }))

    estimates$rel_coef <- ifelse(estimates$p < 0.05, estimates$rel_coef,
        0)
    estimates$signif <- ifelse(estimates$p < 0.05, "p < 0.05", "p >= 0.05")

    return(estimates)
})
estimates <- do.call(rbind, estimates_list)

estimates$snr <- as.numeric(substr(estimates$snr, 5, nchar(estimates$snr)))

rownames(estimates) <- 1:nrow(estimates)

coef_by_signif <- aggregate(coef ~ signif, data = estimates[estimates$predictor ==
    "sign", ], FUN = range)

ggplot(estimates[estimates$predictor == "sign", ], aes(y = coef, x = snr,
    color = freq)) + geom_line(linewidth = 1.6) + labs(x = "Signal-to-noise ratio",
    y = "Effect size when discriminating phase", color = "F0\nfrequency") +
    scale_color_viridis_d(alpha = 0.7, begin = 0.2, end = 0.8) + theme_classic()

Effect sizes for phase by signal-to-noise ratio and frequency.

Discriminating Schroeder’s by number of components grouped by frequency:

Code 
coef_by_signif <- aggregate(coef ~ signif, data = estimates[estimates$predictor ==
    "components", ], FUN = range)


ggplot(estimates[estimates$predictor == "components", ], aes(y = coef,
    x = snr, color = freq)) + geom_line(linewidth = 1.6) + labs(x = "Signal-to-noise ratio",
    y = "Effect size when discriminating number of components", color = "F0\nfrequency") +
    scale_color_viridis_d(alpha = 0.7, begin = 0.2, end = 0.8) + geom_hline(yintercept = mean(coef_by_signif[,
    2][1, 1], coef_by_signif[, 2][2, 2]), lty = 2, color = "gray") +
    theme_classic()

Effect sizes for number of components by signal-to-noise ratio and frequency.
3.1.4.3.1.3 By component
Code 
dtw_dists_snr_list <- readRDS("./data/processed/dtw_distance_schroeders_degradation_experiment_varying_snr.RDS")


dtw_wv_mod_list <- warbleR:::pblapply_wrblr_int(dtw_dists_snr_list,
    cl = 1, function(y) {

        # print(y$target.snr[1])
        comps <- unique(sapply(strsplit(y[, 1], "_"), "[[", 2))

        mods <- lapply(1:length(comps), function(i) {

            # extract only those with the same freq
            x <- y[grepl(paste(comps[i], collapse = "|"), y$schr1) &
                grepl(paste(comps[i], collapse = "|"), y$schr2), ]

            target_snr <- x$target.snr
            x$target.snr <- NULL

            dist_tri <- PhenotypeSpace::rectangular_to_triangular(x)

            if (nrow(dist_tri) > 2) {

                freq_bi_tri <- as.dist(binary_triangular_matrix(group = sapply(strsplit(rownames(dist_tri),
                  "_"), "[[", 1)))

                sign_bi_tri <- as.dist(binary_triangular_matrix(group = sapply(strsplit(rownames(dist_tri),
                  "_"), "[[", 3)))

                rect_var <- cbind(as.dist(dist_tri), freq_bi_tri,
                  sign_bi_tri)

                colnames(rect_var) <- c("dtw_dist", "frequency", "sign")

                dtw_wv_mod <- MRM2(formula = dtw_dist ~ components +
                  sign, nperm = 10000, mat = rect_var)
            } else dtw_wv_mod <- NULL
            return(dtw_wv_mod)
        })

        names(mods) <- paste0("snr:", y$target.snr[1], "_", comps)

        mods <- mods[!sapply(mods, is.null)]

        return(mods)
    })

saveRDS(dtw_wv_mod_list, "./data/processed/matrix_correlation_dtw_distance_varying_snr_by_components.RDS")

Discriminating Schroeder’s by phase grouped by number of components:

Code 
mods <- readRDS("./data/processed/matrix_correlation_dtw_distance_varying_snr_by_components.RDS")

estimates_list <- lapply(mods, function(y) {
    estimates <- do.call(rbind, lapply(seq_along(y), function(x) {
        Y <- data.frame(y[[x]]$coef[-1, ])
        # Y$rel_coef <- Y[, 1] / max(Y[, 1])
        Y$rel_coef <- Y[, 1]
        Y$mod <- names(y)[x]
        Y$predictor <- rownames(Y)
        Y$comps <- strsplit(names(y)[x], "_")[[1]][2]
        Y$comps <- strsplit(Y$comps, ":")[[1]][2]
        names(Y) <- c("coef", "p", "rel_coef", "model", "predictor",
            "comps")
        Y$snr <- gsub("_", "", substr(names(y)[x], 0, 6))

        return(Y)
    }))

    estimates$rel_coef <- ifelse(estimates$p < 0.05, estimates$rel_coef,
        0)
    estimates$signif <- ifelse(estimates$p < 0.05, "p < 0.05", "p >= 0.05")

    return(estimates)
})
estimates <- do.call(rbind, estimates_list)

estimates$snr <- as.numeric(substr(estimates$snr, 5, nchar(estimates$snr)))

estimates$comps <- factor(estimates$comps, levels = unique(as.numeric(estimates$comps)))

rownames(estimates) <- 1:nrow(estimates)

coef_by_signif <- aggregate(coef ~ signif, data = estimates[estimates$predictor ==
    "sign", ], FUN = range)

ggplot(estimates[estimates$predictor == "sign", ], aes(y = coef, x = snr,
    color = comps)) + geom_line(linewidth = 1.6) + labs(x = "Signal-to-noise ratio",
    y = "Effect size when discriminating sign", color = "Components") +
    scale_color_viridis_d(alpha = 0.7, begin = 0.2, end = 0.8) + geom_hline(yintercept = mean(coef_by_signif[,
    2][1, 1], coef_by_signif[, 2][2, 2]), lty = 2, color = "gray") +
    theme_classic()

Effect sizes for phase by signal-to-noise ratio and number of components. Dotted line shows the cutoff for statistical significance.

Discriminating Schroeder’s by frequency grouped by number of components:

Code 
ggplot(estimates[estimates$predictor == "frequency", ], aes(y = coef,
    x = snr, color = comps)) + geom_line(linewidth = 1.6) + labs(x = "Signal-to-noise ratio",
    y = "Effect size when discriminating frequency", color = "Components") +
    scale_color_viridis_d(alpha = 0.7, begin = 0.2, end = 0.8) + theme_classic()

Effect sizes for frequency by signal-to-noise ratio and number of components.
3.1.4.3.2 On single Schroeders
Code 
cor_snr_list <- readRDS("./data/processed/waveform_correlation_adjusted_snr_single_schroeders_white_noise.RDS")

adj_snr_schroeder1 <- readRDS("./data/processed/200ms_schroeders_adjusted_snr_white_noise.RDS")[[1]]

adj_snr_schroeder1$sound.id <- gsub("f0:f0 ", "f0:", adj_snr_schroeder1$sound.id)


cor_mod_list <- warbleR:::pblapply_wrblr_int(names(cor_snr_list),
    cl = 20, function(y) {

        x <- cor_snr_list[[y]]
        target_snr <- "snr:30"

        x$schr1 <- sapply(strsplit(x$schr1, "_"), function(x) paste(x[2],
            x[3], substr(x[4], 0, 1), sep = "_"))
        x$schr2 <- sapply(strsplit(x$schr2, "_"), function(x) paste(x[2],
            x[3], substr(x[4], 0, 1), sep = "_"))
        x$cor <- x$cor * -1

        dist_tri <- PhenotypeSpace::rectangular_to_triangular(x)

        freq_bi_tri <- as.dist(binary_triangular_matrix(group = sapply(strsplit(rownames(dist_tri),
            "_"), "[[", 1)))

        comp_bi_tri <- as.dist(binary_triangular_matrix(group = gsub("_n|_p",
            "", sapply(strsplit(rownames(dist_tri), "_"), "[[", 2))))

        sign_bi_tri <- as.dist(binary_triangular_matrix(group = sapply(strsplit(rownames(dist_tri),
            "_"), "[[", 3)))

        rect_var <- cbind(as.dist(dist_tri), freq_bi_tri, comp_bi_tri,
            sign_bi_tri)

        colnames(rect_var) <- c("cor", "frequency", "components",
            "sign")

        cor_wv_mod <- MRM2(formula = cor ~ frequency + components +
            sign, nperm = 10000, mat = rect_var)

        return(cor_wv_mod)
    })

names(cor_mod_list) <- names(cor_snr_list)

saveRDS(cor_mod_list, "./data/processed/matrix_regression_waveform_correlation_varying_snr_single_schroeder_white_noise.RDS")
Code 
mods <- readRDS("./data/processed/matrix_regression_waveform_correlation_varying_snr_single_schroeder_white_noise.RDS")


estimates_wf_wn <- do.call(rbind, lapply(seq_along(mods), function(x) {
    Y <- data.frame(mods[[x]]$coef[-1, ])
    # Y$rel_coef <- Y[, 1] / max(Y[, 1])
    Y$rel_coef <- Y[, 1]
    Y$predictor <- rownames(Y)
    Y$model <- as.numeric(gsub("snr:", "", names(mods)[x]))
    names(Y) <- c("coef", "p", "rel_coef", "predictor", "model")
    return(Y)
}))

estimates_wf_wn$rel_coef <- ifelse(estimates_wf_wn$p < 0.05, estimates_wf_wn$rel_coef,
    0)
estimates_wf_wn$signif <- ifelse(estimates_wf_wn$p < 0.05, "p < 0.05",
    "p >= 0.05")
estimates_wf_wn$signif2 <- ifelse(estimates_wf_wn$p < 0.05, "*", "")

estimates_wf_wn$model_f <- factor(estimates_wf_wn$model, levels = 1:30)

ggplot(estimates_wf_wn, aes(x = predictor, y = model_f, fill = rel_coef)) +
    geom_tile() + coord_equal() + scale_fill_gradient2(low = viridis(10)[3],
    high = viridis(10)[7], guide = "none") + geom_text(aes(label = signif2,
    color = signif), size = 3) + scale_color_manual(values = c("black",
    "gray")) + labs(y = "Signal-to-noise ratio", x = "Predictor",
    color = "P value") + theme_classic() + theme(axis.text.x = element_text(color = "black",
    size = 7, angle = 30, vjust = 0.8, hjust = 0.8)) + coord_flip()

Effect sizes by signal-to-noise ratio and predictor (color intensity shows effect size magnitude relative to the highest effect size within the model).
Code 
estimates_wf_wn$model <- as.numeric(estimates_wf_wn$model)
Code 
coef_by_signif <- aggregate(coef ~ signif, data = estimates_wf_wn,
    FUN = range)

ggplot(estimates_wf_wn, aes(y = coef, x = model, color = predictor,
    lty = predictor)) + geom_line(linewidth = 1.6) + labs(x = "Signal-to-noise ratio",
    y = "Effect size") + scale_color_viridis_d(alpha = 0.7, begin = 0.2,
    end = 0.8) + geom_hline(yintercept = mean(coef_by_signif[, 2][1,
    1], coef_by_signif[, 2][2, 2]), lty = 2, color = "gray") + theme_classic()

Effect sizes by signal-to-noise ratio and predictor waveform correlation

3.1.5 Measure linear-frequency spectrogram cross-correlation

Code 
adj_snr_schroeder <- readRDS("./data/processed/200ms_schroeders_adjusted_snr_white_noise.RDS")

xc_list <- lapply(adj_snr_schroeder, function(schroeders) {

    print(schroeders$target.snr[1])

    xc <- cross_correlation(schroeders, ovlp = 90, parallel = 20,
        wn = "hamming")

})

names(xc_list) <- paste0("snr:", sapply(adj_snr_schroeder, function(x) x$target.snr[1]))

saveRDS(xc_list, "./data/processed/fourier_cross_correlation_varying_snr_200ms_schroeders_white_noise.RDS")


beepr::beep(2)
Code 
xc_list <- readRDS("./data/processed/fourier_cross_correlation_varying_snr_200ms_schroeders_white_noise.RDS")

adj_snr_schroeder1 <- readRDS("./data/processed/200ms_schroeders_adjusted_snr_pink_noise.RDS")[[1]]

adj_snr_schroeder1$sound.id <- gsub("f0:f0 ", "f0:", adj_snr_schroeder1$sound.id)

xc_mod_list <- warbleR:::pblapply_wrblr_int(seq_along(xc_list), cl = 20,
    function(x) {

        X <- xc_list[[x]]

        # x$dist <- scale(x$dist)
        target_snr <- gsub("snr:", "", names(xc_list)[x])

        colnames(X) <- sapply(gsub("-1$", "", colnames(X)), function(x) adj_snr_schroeder1$sound.id[adj_snr_schroeder1$sound.files ==
            x])
        rownames(X) <- sapply(gsub("-1$", "", rownames(X)), function(x) adj_snr_schroeder1$sound.id[adj_snr_schroeder1$sound.files ==
            x])
        # convert to distance
        X <- 1 - X

        # cor_tri <-
        # PhenotypeSpace::rectangular_to_triangular(x)

        freq_bi_tri <- as.dist(binary_triangular_matrix(group = sapply(strsplit(rownames(X),
            "_"), "[[", 1)))

        comp_bi_tri <- as.dist(binary_triangular_matrix(group = gsub("_n|_p",
            "", sapply(strsplit(rownames(X), "_"), "[[", 2))))

        sign_bi_tri <- as.dist(binary_triangular_matrix(group = sapply(strsplit(rownames(X),
            "_"), "[[", 3)))

        rect_var <- cbind(as.dist(X), freq_bi_tri, comp_bi_tri, sign_bi_tri)

        colnames(rect_var) <- c("xcor", "frequency", "components",
            "sign")

        dtw_wv_mod <- MRM2(formula = xcor ~ frequency + components +
            sign, nperm = 10000, mat = rect_var)

        return(dtw_wv_mod)
    })

names(xc_mod_list) <- names(xc_list)

saveRDS(xc_mod_list, "./data/processed/matrix_regression_cross_correlation_varying_snr_white_noise.RDS")

Discriminating Schroeder’s by phase grouped by number of components:

Code 
mods <- readRDS("./data/processed/matrix_regression_cross_correlation_varying_snr_white_noise.RDS")

estimates_xc_wn <- do.call(rbind, lapply(seq_along(mods), function(x) {
    Y <- data.frame(mods[[x]]$coef[-1, ])
    # Y$rel_coef <- Y[, 1] / max(Y[, 1])
    Y$rel_coef <- Y[, 1]
    Y$predictor <- rownames(Y)
    Y$model <- as.numeric(gsub("snr:", "", names(mods)[x]))
    names(Y) <- c("coef", "p", "rel_coef", "predictor", "model")
    return(Y)
}))

estimates_xc_wn$rel_coef <- ifelse(estimates_xc_wn$p < 0.05, estimates_xc_wn$rel_coef,
    0)
estimates_xc_wn$signif <- ifelse(estimates_xc_wn$p < 0.05, "p < 0.05",
    "p >= 0.05")
estimates_xc_wn$signif2 <- ifelse(estimates_xc_wn$p < 0.05, "*", "")

estimates_xc_wn$model_f <- factor(estimates_xc_wn$model, levels = 1:30)

ggplot(estimates_xc_wn, aes(x = predictor, y = model_f, fill = rel_coef)) +
    geom_tile() + coord_equal() + scale_fill_gradient2(low = viridis(10)[3],
    high = viridis(10)[7], guide = "none") + geom_text(aes(label = signif2,
    color = signif), size = 3) + scale_color_manual(values = c("black",
    "gray")) + labs(y = "Signal-to-noise ratio", x = "Predictor",
    color = "P value") + theme_classic() + theme(axis.text.x = element_text(color = "black",
    size = 7, angle = 30, vjust = 0.8, hjust = 0.8)) + coord_flip()

Effect sizes for phase by signal-to-noise ratio and number of components. Dotted line shows the cutoff for statistical significance.
Code 
estimates_xc_wn$model <- as.numeric(estimates_xc_wn$model)
Code 
coef_by_signif3 <- aggregate(coef ~ signif, data = estimates_xc_wn,
    FUN = range)

ggplot(estimates_xc_wn, aes(y = coef, x = model, color = predictor,
    lty = predictor)) + geom_line(linewidth = 1.6) + labs(x = "Signal-to-noise ratio",
    y = "Effect size") + scale_color_viridis_d(alpha = 0.7, begin = 0.2,
    end = 0.8) + theme_classic()

Effect sizes by signal-to-noise ratio and predictor waveform correlation

3.1.5.1 Only phase

Code 
ggplot(estimates_xc_wn[estimates_xc_wn$predictor == "sign", ], aes(y = abs(coef),
    x = model, color = predictor, lty = predictor)) + geom_line(linewidth = 1.6) +
    geom_smooth() + labs(x = "Signal-to-noise ratio", y = "Effect size") +
    scale_color_viridis_d(alpha = 0.7, begin = 0.2, end = 0.8) + geom_hline(yintercept = 0,
    lty = 2, color = "gray") + theme_classic()

Effect sizes by signal-to-noise ratio and predictor.

4 Combined results plot

Code 
estimates_wf_wn$method <- "waveform correlation"
estimates_dtw_wn$method <- "waveform DTW"
estimates_xc_wn$method <- "cross-correlation"

estims <- rbind(estimates_wf_wn, estimates_dtw_wn, estimates_xc_wn)

coef_by_signif2 <- aggregate(coef ~ signif + method, data = estims,
    FUN = range)


ggplot(estims[estims$predictor == "sign", ], aes(y = abs(coef), x = model,
    color = method)) + geom_line(linewidth = 1.6) + geom_smooth() +
    labs(x = "Signal-to-noise ratio", y = "Phase\neffect size") +
    scale_color_viridis_d(alpha = 0.7, begin = 0.2, end = 0.8) + xlim(c(1,
    30)) + geom_hline(yintercept = 0, lty = 2, color = "gray") + theme_classic() +
    scale_x_continuous(expand = c(0, 0))

5 Takeaways

 

Session information

R version 4.4.2 (2024-10-31)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.5 LTS

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

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       

time zone: America/Costa_Rica
tzcode source: system (glibc)

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

other attached packages:
 [1] numform_0.7.0        ecodist_2.1.3        PhenotypeSpace_0.1.0
 [4] ggplot2_3.5.1        baRulho_2.1.4        ohun_1.0.2          
 [7] Rraven_1.0.14        viridis_0.6.5        viridisLite_0.4.2   
[10] warbleR_1.1.34       NatureSounds_1.0.5   tuneR_1.4.7         
[13] seewave_2.2.3        formatR_1.14         knitr_1.50          
[16] kableExtra_1.4.0    

loaded via a namespace (and not attached):
 [1] DBI_1.2.3              bitops_1.0-9           deldir_2.0-4          
 [4] pbapply_1.7-2          gridExtra_2.3          remotes_2.5.0         
 [7] permute_0.9-7          testthat_3.2.3         rlang_1.1.5           
[10] magrittr_2.0.3         e1071_1.7-16           compiler_4.4.2        
[13] spatstat.geom_3.3-2    mgcv_1.9-1             png_0.1-8             
[16] systemfonts_1.1.0      vctrs_0.6.5            Sim.DiffProc_4.9      
[19] stringr_1.5.1          pkgconfig_2.0.3        crayon_1.5.3          
[22] fastmap_1.2.0          backports_1.5.0        labeling_0.4.3        
[25] rmarkdown_2.29         xfun_0.52              jsonlite_2.0.0        
[28] goftest_1.2-3          spatstat.utils_3.1-0   terra_1.7-78          
[31] Deriv_4.1.6            parallel_4.4.2         cluster_2.1.6         
[34] R6_2.6.1               stringi_1.8.7          spatstat.data_3.1-2   
[37] spatstat.univar_3.0-1  brio_1.1.5             Rcpp_1.0.14           
[40] tensor_1.5             splines_4.4.2          Matrix_1.7-0          
[43] igraph_2.1.4           tidyselect_1.2.1       rstudioapi_0.16.0     
[46] abind_1.4-8            yaml_2.3.10            vegan_2.6-8           
[49] codetools_0.2-20       dtw_1.23-1             spatstat.random_3.3-1 
[52] spatstat.explore_3.3-2 curl_6.2.2             lattice_0.22-6        
[55] tibble_3.2.1           withr_3.0.2            evaluate_1.0.3        
[58] signal_1.8-1           sf_1.0-20              sketchy_1.0.5         
[61] units_0.8-7            proxy_0.4-27           polyclip_1.10-7       
[64] xml2_1.3.6             pillar_1.10.1          packrat_0.9.2         
[67] KernSmooth_2.23-24     checkmate_2.3.2        generics_0.1.3        
[70] sp_2.1-4               RCurl_1.98-1.17        munsell_0.5.1         
[73] scales_1.3.0           class_7.3-22           glue_1.8.0            
[76] tools_4.4.2            xaringanExtra_0.8.0    grid_4.4.2            
[79] colorspace_2.1-1       nlme_3.1-165           raster_3.6-26         
[82] cli_3.6.4              spatstat.sparse_3.1-0  svglite_2.1.3         
[85] dplyr_1.1.4            gtable_0.3.6           fftw_1.0-9            
[88] digest_0.6.37          classInt_0.4-11        farver_2.1.2          
[91] rjson_0.2.23           htmlwidgets_1.6.4      htmltools_0.5.8.1     
[94] lifecycle_1.0.4        httr_1.4.7             MASS_7.3-61