Measuring degradation
Simulated song degradation Manu 2022
Universidad de Costa Rica
Universidad de Costa Rica
Marcelo
Araya-Salas
“2022-10-11”
cutw_fun <- function(wave, from, to) {
wave1 <- wave@left[from:to]
wave1 <- Wave(left = wave1, samp.rate = wave@samp.rate, bit = wave@bit,
pcm = TRUE)
return(wave1)
}
filter_sound <- function(wave = x, start, end, bottom.freq, top.freq,
buffer = 0) {
start <- start - buffer
end <- end + buffer
bfr_wv <- cutw_fun(wave, from = 1, to = (start * wave@samp.rate) -
1)
bfr_wv <- normalize(object = bfr_wv)
aftr_wv <- cutw_fun(wave, from = (end * wave@samp.rate) + 1, to = length(wave@left))
aftr_wv <- normalize(object = aftr_wv)
wv <- cutw_fun(wave, from = start * wave@samp.rate, to = end *
wave@samp.rate)
wv <- normalize(object = wv)
flt_wv <- ffilter(wv, from = bottom.freq * 1000, to = top.freq *
1000, output = "Wave", bandpass = FALSE)
flt_wv <- normalize(object = flt_wv)
wv1 <- pastew(wave1 = flt_wv, wave2 = bfr_wv, output = "Wave")
out_wv <- pastew(wave1 = aftr_wv, wave2 = wv1, output = "Wave")
return(out_wv)
}
standard_error <- function(x) sd(x)/sqrt(length(x))
cols <- viridis(10, alpha = 0.7)
# print results brms models
# print results brms models
summary_brm_model <- function(x = NULL, gsub.pattern = NULL, gsub.replacement = NULL,
xlab = "Effect size", n.subposts = 2000, model_name = NULL, read.file = NULL,
plot.area.prop = 0.9, grep.filter.out = NULL) {
if (is.null(x) & !is.null(read.file))
x <- readRDS(read.file)
# extract info from model
summ <- summary(x)$fixed
fit <- x$fit
betas <- grep("^b_", names(fit@sim$samples[[1]]), value = TRUE)
if (!is.null(grep.filter.out)) {
betas <- grep(grep.filter.out, betas, invert = TRUE, value = TRUE)
summ <- summ[grep(grep.filter.out, rownames(summ), invert = TRUE),
]
}
# subsample posteriors
xdrws <- brms::as_draws(x)
# only apply thinning if length of posterior < n.subposts
if (length(xdrws[[1]][[1]]) < n.subposts)
xdrws <- posterior::thin_draws(xdrws, round(length(xdrws[[1]][[1]])/n.subposts,
0))
xdrws <- posterior::subset_draws(x = xdrws, variable = betas)
sub_posts_by_chain_list <- lapply(1:length(xdrws), function(x) {
X <- as.data.frame(xdrws[[x]])
X$chain <- paste("chain", x)
return(X)
})
sub_posts_by_chain <- do.call(rbind, sub_posts_by_chain_list)
merged_xdrws <- posterior::merge_chains(xdrws)
sub_posts <- as.data.frame(merged_xdrws)
names(sub_posts) <- betas
hdis <- t(sapply(betas, function(y) HDInterval::hdi(sub_posts[,
colnames(sub_posts) == y])))
coef_table <- data.frame(summ, hdis)
coef_table <- coef_table[rownames(coef_table) != "Intercept",
c("Estimate", "Rhat", "Bulk_ESS", "l.95..CI", "u.95..CI")]
# add priors to model table
pt <- prior_summary(x)
b_prior <- pt$prior[pt$class == "b" & pt$coef == ""]
b_prior <- if (b_prior == "")
"flat" else b_prior
sd_prior <- unique(pt$prior[pt$class == "sd" & pt$coef == ""])
sd_prior <- if (length(sd_prior) > 1)
sd_prior[sd_prior != ""] else "flat"
sd_prior <- if (sd_prior == "")
"flat" else sd_prior
model_table <- data.frame(b_prior, sd_prior, iterations = fit@stan_args[[1]]$iter,
chains = length(attr(fit, "stan_args")), thinning = fit@stan_args[[1]]$thin,
warmup = fit@stan_args[[1]]$warmup)
np <- brms::nuts_params(x)
model_table$diverg_transitions <- sum(subset(np, Parameter ==
"divergent__")$Value)
model_table$`rhats > 10.5` <- sum(na.omit(brms::rhat(x)) > 1.05)
coef_table <- as.data.frame(coef_table)
coef_table$Rhat <- round(coef_table$Rhat, digits = 3)
coef_table$CI_low <- round(unlist(coef_table$l.95..CI), digits = 3)
coef_table$CI_high <- round(unlist(coef_table$u.95..CI), digits = 3)
coef_table$l.95..CI <- coef_table$u.95..CI <- NULL
out <- lapply(betas, function(y) data.frame(variable = y, value = sort(sub_posts[,
colnames(sub_posts) == y], decreasing = FALSE)))
posteriors <- do.call(rbind, out)
posteriors <- posteriors[posteriors$variable != "b_Intercept",
]
posteriors$variable <- factor(posteriors$variable, levels = sort(unique(posteriors$variable)))
names(sub_posts_by_chain)[1:(ncol(sub_posts_by_chain) - 1)] <- betas
out2 <- lapply(betas, function(y) {
X <- data.frame(variable = y, chain = sub_posts_by_chain[,
"chain"], value = sub_posts_by_chain[, y])
X$iteration <- round(seq(1, fit@stan_args[[1]]$iter, length.out = nrow(X)/length(attr(fit,
"stan_args"))))
return(X)
})
posteriors_by_chain <- do.call(rbind, out2)
posteriors_by_chain <- posteriors_by_chain[posteriors_by_chain$variable !=
"b_Intercept", ]
if (!is.null(gsub.pattern) & !is.null(gsub.replacement))
posteriors$variable <- gsub(pattern = gsub.pattern, replacement = gsub.replacement,
posteriors$variable)
if (!is.null(gsub.pattern) & !is.null(gsub.replacement))
posteriors_by_chain$variable <- gsub(pattern = gsub.pattern,
replacement = gsub.replacement, posteriors$posteriors_by_chain)
posteriors_by_chain$variable <- factor(posteriors_by_chain$variable,
levels = sort(levels(posteriors$variable), decreasing = TRUE))
coef_table2 <- coef_table
coef_table2$variable <- factor(paste0("b_", rownames(coef_table2)))
coef_table2$value <- coef_table2$Estimate
coef_table2$significance <- ifelse(coef_table2$CI_low * coef_table2$CI_high >
0, "sig", "non-sig")
coef_table2$significance <- factor(coef_table2$significance, levels = c("non-sig",
"sig"))
col_pointrange <- if (all(coef_table2$significance == "non-sig"))
"gray" else if (all(coef_table2$significance == "sig"))
"black" else c("gray", "black")
gg_dists <- ggplot(data = posteriors, aes(y = variable, x = value,
color = significance)) + geom_vline(xintercept = 0, col = "black",
lty = 2) + ggdist::stat_halfeye(aes(x = value), .width = c(0.95),
fill = viridis::viridis(10, alpha = 0.5)[8], normalize = "panels",
color = "transparent") + geom_point(data = coef_table2) +
geom_errorbar(data = coef_table2, aes(xmin = CI_low, xmax = CI_high),
width = 0) + scale_color_manual(values = col_pointrange) +
theme(axis.ticks.length = unit(0, "pt"), plot.margin = margin(0,
0, 0, 0, "pt")) + labs(x = "Effect size", y = "Parameter") +
theme_classic() + facet_wrap(~variable, scales = "free_y",
ncol = 1, strip.position = "right") + theme(legend.position = "none",
strip.background = element_blank(), strip.text.y = element_blank()) +
xlim(range(c(posteriors_by_chain$value, 0)) * plot.area.prop)
gg_traces <- ggplot(data = posteriors_by_chain, aes(x = iteration,
y = value, color = chain)) + geom_line() + scale_color_viridis_d(alpha = 0.7,
begin = 0.2, end = 0.9) + facet_wrap(~variable, scales = "free_y",
ncol = 1, strip.position = "right") + theme_classic() + theme(legend.position = "none",
strip.background = element_blank(), strip.text.y = element_blank(),
axis.title.y = element_blank(), axis.text.y = element_blank())
gg <- cowplot::plot_grid(gg_dists, gg_traces, ncol = 2, rel_widths = c(2,
1))
if (!is.null(gsub.pattern) & !is.null(gsub.replacement))
rownames(coef_table) <- gsub(pattern = gsub.pattern, replacement = gsub.replacement,
rownames(coef_table))
coef_table$Rhat <- ifelse(coef_table$Rhat > 1.05, cell_spec(coef_table$Rhat,
"html", color = "white", background = "red", bold = TRUE,
font_size = 12), cell_spec(coef_table$Rhat, "html"))
signif <- coef_table[, "CI_low"] * coef_table[, "CI_high"] > 0
model_table$diverg_transitions <- ifelse(model_table$diverg_transitions >
0, cell_spec(model_table$diverg_transitions, "html", color = "white",
background = "red", bold = TRUE, font_size = 12), cell_spec(model_table$diverg_transitions,
"html"))
model_table$`rhats > 10.5` <- ifelse(model_table$`rhats > 10.5` >
0, cell_spec(model_table$`rhats > 10.5`, "html", color = "white",
background = "red", bold = TRUE, font_size = 12), cell_spec(model_table$`rhats > 10.5`,
"html"))
df1 <- kbl(model_table, row.names = TRUE, escape = FALSE, format = "html",
digits = 3)
df1 <- kable_styling(df1, bootstrap_options = c("striped", "hover",
"condensed", "responsive"), full_width = FALSE, font_size = 12)
df2 <- kbl(coef_table, row.names = TRUE, escape = FALSE, format = "html",
digits = 3)
df2 <- row_spec(kable_input = df2, row = which(coef_table$CI_low *
coef_table$CI_high > 0), background = adjustcolor(cols[9],
alpha.f = 0.3))
df2 <- kable_styling(df2, bootstrap_options = c("striped", "hover",
"condensed", "responsive"), full_width = FALSE, font_size = 12)
if (!is.null(model_name))
cat(paste("<font size=\"4\"><b>", model_name, "</b></font><br>"))
cat(paste("<font size=\"3\"><b>", x$formula[1], "</b></font>"))
print(df1)
print(df2)
print(gg)
}Purpose
Prepare degradation experiment data
Align re-recorded test sounds
Measure degradation
Do stats
Load packages:
Simulate sounds:
# simulation parameters
gap <- 0.05
set.seed(123)
frqs <- sample(rep(seq(0.5, 10, length.out = 20), 3))
steps <- 11
sig2 <- 0.3
# amplitude modulation with 2 peaks
mod.amps <- c(0.5, 1, 2, 3, 2, 1, 2, 3, 2, 1, 0.5)
# make all posible combinations
eg <- expand.grid(dur = c(0.1, 0.2), fm = c("pure.tone", "BB"), am = c("no.am",
"am"), stringsAsFactors = FALSE)
pboptions(type = "timer")
setwd("./data/raw/simulated_songs")
# simulate songs
sim.songs <- pblapply(1:nrow(eg), function(x) {
if (eg$am[x] == "no.am")
am.amps <- 1 else am.amps <- mod.amps
# if (eg$harm[x] == 'no.harm') harms <- 1 else harms <-
# nharms
sm.sng <- sim_songs(n = length(frqs), durs = eg$dur[x], freqs = frqs,
gaps = gap, am.amps = am.amps, harms = 1, diff.fun = eg$fm[x],
selec.table = TRUE, sig2 = sig2, steps = steps, file.name = paste(eg[x,
], collapse = "_"), bgn = 0, seed = seed)
# add freq room if pure tone
if (eg$fm[x] == "pure.tone") {
sm.sng$selec.table$bottom.freq <- sm.sng$selec.table$bottom.freq -
0.2
sm.sng$selec.table$top.freq <- sm.sng$selec.table$top.freq +
0.2
}
sm.sng$selec.table$bottom.freq[sm.sng$selec.table$bottom.freq <
0] <- 0.1
sm.sng$selec.table$sim.freq <- as.character(frqs)
return(sm.sng)
})
setwd("../../")
# name with parameters
names(sim.songs) <- sapply(1:nrow(eg), function(x) paste(eg[x, ],
collapse = "_"))
# shuffle sim.songs
sim.songs <- sim.songs[sample(1:length(sim.songs))]
# plot spectros
out <- pblapply(1:length(sim.songs), function(x) {
X <- sim.songs[[x]]
jpeg(filename = file.path("./data/raw/simulated_song_images/",
paste0(names(sim.songs)[x], ".jpeg")), width = 1100)
spectro(sim.songs[[x]]$wave, scale = FALSE, palette = reverse.gray.colors.1,
grid = FALSE, flim = c(0, 10.5), collevels = seq(-20, 0, 1),
main = names(sim.songs)[x], osc = TRUE, fastdisp = TRUE)
dev.off()
})
# extract select tables
sim.song.sts <- lapply(sim.songs, function(X) X$selec.table)
sim.song.st <- do.call(rbind, sim.song.sts)
cs <- check_sels(sim.song.st, path = "./data/raw/simulated_songs")
cs[cs$check.res != "OK", ]
# make a single extended selection table for simulation
sim.song.est <- selection_table(X = sim.song.st, extended = TRUE,
pb = FALSE, confirm.extended = FALSE, path = "./data/raw/simulated_songs")
# save
saveRDS(sim.song.est, "./data/raw/consolidated.est.RDS")Create master sound files:
cons.est <- readRDS("./data/raw/consolidated.est.RDS")
# create master sound file
master.sf <- master_sound_file(X = cons.est, file.name = "consolidated_master",
gap.duration = 0.05, dest.path = "./data/raw", overwrite = TRUE,
amp.marker = 3)
nrow(master.sf)
max(master.sf$end)/60
# lspec(flist = 'consolidated_master.wav')
lspec(X = master.sf, path = "./data/raw")
write.csv(master.sf, "./data/raw/consolidated.master.sf_selection_table.csv",
row.names = FALSE)rerecs <- list.files("./data/raw/re-recorded/", pattern = "WAV$")
rerec_df <- data.frame(sound.files = rerecs, transect = substr(0,
2, x = rerecs), habitat = substr(4, 5, x = rerecs), mic.height = as.numeric(substr(8,
8, x = rerecs)), speak.height = as.numeric(substr(11, 11, x = rerecs)),
distance = as.numeric(substr(14, 15, x = rerecs)))
write.csv(rerec_df, file = "./data/raw/rerecorded_files_metadata.csv",
row.names = FALSE)master.sf <- read.csv("./data/raw/consolidated.master.sf_selection_table.csv")
rerec_df <- read.csv(file = "./data/raw/rerecorded_files_metadata.csv")
setwd("/home/m/Dropbox/Projects/sound_degradation_manu/data/raw/re-recorded")
all(rerec_df$sound.files %in% rerecs)
warbleR_options(path = "/home/m/Dropbox/Projects/sound_degradation_manu/data/raw/re-recorded")
found.starts <- search_templates(X = master.sf, template.rows = which(master.sf$orig.sound.file ==
"start_marker"), test.files = rerec_df$sound.files, pb = TRUE,
path = "./data/raw/re-recorded", parallel = 3)
master.sf2 <- master.sf
master.sf2$end <- master.sf2$end - min(master.sf2$start)
master.sf2$start <- master.sf2$start - min(master.sf2$start)
master.sf2 <- master.sf2[!master.sf2$orig.sound.file %in% c("start_marker",
"end_marker"), ]
out <- lapply(1:nrow(found.starts), function(x) {
Y <- master.sf2
Y$sound.files <- found.starts$test.files[x]
Y$start <- Y$start + found.starts$start[x]
Y$end <- Y$end + found.starts$start[x]
return(Y)
})
align_seltab <- do.call(rbind, out)
exp_raven(X = align_seltab, path = "./data/raw/re-recorded", sound.file.path = "./data/raw/re-recorded")
write.csv(align_seltab, file = "./data/processed/aligned_test_files_selection_table.csv",
row.names = FALSE)
alg.tests <- align_test_files(X = master.sf, Y = found.starts, path = "./data/raw/re-recorded",
by.song = TRUE)
alg.tests$amplitude.modulation <- ifelse(grepl("no.am", alg.tests$template),
"no.am", "am")
alg.tests$frequency.modulation <- ifelse(grepl("BB", alg.tests$template),
"no.fm", "fm")
alg.tests$duration <- ifelse(alg.tests$end - alg.tests$start < 0.15,
"short", "long")
saveRDS(alg.tests, file = "./data/processed/algined_test_files_extended_selection_table.RDS")rerecs <- list.files("./data/raw/re-recorded", pattern = "WAV$", full.names = TRUE)
align_seltab <- read.csv("./data/processed/aligned_test_files_selection_table.csv")
out <- pblapply(rerecs, function(x) {
wv <- readWave(x)
st <- align_seltab[align_seltab$sound.files == basename(x), ]
for (i in 1:nrow(st)) wv <- filter_sound(wave = wv, start = st$start[i],
end = st$end[i], bottom.freq = st$bottom.freq[i], top.freq = st$top.freq[i],
buffer = 0.05)
flt_wv <- ffilter(wv, from = 2.5 * 1000, to = 6.7 * 1000, output = "Wave",
bandpass = TRUE)
flt_wv <- normalize(object = flt_wv, unit = "16")
env <- env(flt_wv, plot = FALSE, envt = "abs", ssmooth = 2000)
plot(env, type = "l")
wv <- readWave(x)
env2 <- env(wv, plot = FALSE, envt = "abs", ssmooth = 41000)
plot(env2, type = "l")
# writeWave(object = flt_wv, filename =
# file.path('./data/raw/filtered_re-recorded', basename(x)),
# extensible = FALSE)
})# spectro(wv, tlim = c(4.2, 5.45), flim = c(2.7, 6.2), scale =
# FALSE, osc = FALSE)
cicadas <- imp_raven(path = "./data/raw", files = "T0_HA_M3_P3_D01_cicadas.Table.1.selections.txt",
warbler.format = TRUE)
correlations <- template_correlator(templates = cicadas, files = "T0_HA_M3_P3_D01.WAV",
path = "./data/raw/re-recorded")
optimization <- optimize_template_detector(template.correlations = correlations,
reference = cicadas, threshold = seq(0.1, 0.9, 0.05))
optimization <- optimization[order(optimization$f1.score, decreasing = TRUE),
]
best_temp <- cicadas[cicadas$sound.files == "T0_HA_M3_P3_D01.WAV" &
cicadas$selec == 3, ]
align_seltab <- read.csv("./data/processed/aligned_test_files_selection_table.csv")
# get correlations
correlations <- template_correlator(templates = best_temp, files = unique(align_seltab$sound.files),
path = "./data/raw/re-recorded", pb = TRUE)
# run detection
detection <- template_detector(template.correlations = correlations,
threshold = 0.4)
tab <- table(detection$sound.files)
sort(tab)
high.tab <- tab[tab > 15]
# more than 15 detected cicada calls is high cicada noise, less
# is low cicada noise
alg.tests <- readRDS("./data/processed/algined_test_files_extended_selection_table.RDS")
alg.tests$cicadas <- ifelse(sapply(alg.tests$sound.files, function(x) strsplit(x,
"-song")[[1]][[1]]) %in% names(high.tab), "high", "low")
rerec_df <- read.csv(file = "./data/raw/rerecorded_files_metadata.csv")
alg.tests$habitat <- sapply(sapply(alg.tests$sound.files, function(x) strsplit(x,
"-song")[[1]][[1]]), function(x) rerec_df$habitat[rerec_df$sound.files ==
x])
alg.tests$mic.height <- sapply(sapply(alg.tests$sound.files, function(x) strsplit(x,
"-song")[[1]][[1]]), function(x) rerec_df$mic.height[rerec_df$sound.files ==
x])
alg.tests$transect <- sapply(sapply(alg.tests$sound.files, function(x) strsplit(x,
"-song")[[1]][[1]]), function(x) rerec_df$transect[rerec_df$sound.files ==
x])
alg.tests$speak.height <- sapply(sapply(alg.tests$sound.files, function(x) strsplit(x,
"-song")[[1]][[1]]), function(x) rerec_df$speak.height[rerec_df$sound.files ==
x])
alg.tests$distance <- sapply(sapply(alg.tests$sound.files, function(x) strsplit(x,
"-song")[[1]][[1]]), function(x) rerec_df$distance[rerec_df$sound.files ==
x])
alg.tests$site <- NA
alg.tests$site[alg.tests$transect %in% c("T1", "T2")] <- "S1"
alg.tests$site[alg.tests$transect %in% c("T3", "T4")] <- "S2"
alg.tests$site[alg.tests$transect %in% c("T5", "T6")] <- "S3"
alg.tests$site[alg.tests$transect %in% c("T7", "T8")] <- "S4"
alg.tests <- alg.tests[grep("marker", alg.tests$template, invert = TRUE),
]
set.seed(123)
frqs <- sample(rep(seq(0.5, 10, length.out = 20), 3))
alg.tests$target.frequency <- frqs
nrow(alg.tests)/480
alg.tests$treatment.replicates <- paste("frq:", alg.tests$target.frequency,
"dur:", alg.tests$duration, alg.tests$amplitude.modulation, alg.tests$frequency.modulation,
sep = "_")
saveRDS(alg.tests, file = "./data/processed/algined_test_files_extended_selection_table.RDS")Measuring degradation:
alg.tests <- readRDS("./data/processed/algined_test_files_extended_selection_table.RDS")
alg.tests <- signal_to_noise_ratio(alg.tests, mar = 0.02, pb = TRUE,
parallel = 10, type = 1, wl = 512, ovlp = 50)
alg.tests <- blur_ratio(alg.tests, pb = TRUE, parallel = 10, wl = 256,
ovlp = 50, ssmooth = NULL)
alg.tests <- tail_to_signal_ratio(alg.tests, pb = TRUE, parallel = 10,
mar = 0.04, type = 1, wl = 512, ovlp = 99)
alg.tests <- excess_attenuation(X = alg.tests, pb = FALSE, parallel = 10,
type = "Dabelsteen", ovlp = 50, wl = 1024)
alg.tests <- spcc(X = alg.tests, pb = TRUE, parallel = 10, ovlp = 50,
wl = 1024)
alg.tests <- envelope_correlation(X = alg.tests, pb = TRUE, parallel = 1,
ovlp = 50, wl = 1024)
alg.tests <- baRulho::spectrum_blur_ratio(X = alg.tests, pb = TRUE,
parallel = 10, ovlp = 50, wl = 1024)
alg.tests <- baRulho::spectrum_correlation(X = alg.tests, pb = TRUE,
parallel = 10, ovlp = 50, wl = 1024)
write.csv(alg.tests, file = "./data/processed/degradation_measures_test_files_extended_selection_table.csv")
saveRDS(alg.tests, file = "./data/processed/degradation_measures_test_files_extended_selection_table.RDS")Statistical analysis
degrad <- read.csv("./data/processed/degradation_measures_test_files_extended_selection_table.csv")
degrad$frequency <- (degrad$top.freq + degrad$bottom.freq)/2
degrad <- degrad[degrad$distance == 30, ]
pca <- prcomp(scale(degrad[, c("blur.ratio", "tail.to.signal.ratio",
"signal.to.noise.ratio", "excess.attenuation", "cross.correlation",
"spectral.blur.ratio", "spectrum.correlation", "envelope.correlation")]))
summary(pca)
pca$rotation
degrad$pc1 <- pca$x[, 1]
degrad$frequency.modulation <- factor(degrad$frequency.modulation,
levels = c("no.fm", "fm"))
degrad$amplitude.modulation <- factor(degrad$amplitude.modulation,
levels = c("no.am", "am"))
degrad$cicadas <- factor(degrad$cicadas, levels = c("low", "high"))
form <- pc1 ~ habitat + frequency + frequency.modulation + amplitude.modulation +
cicadas + speak.height * mic.height + duration + (1 | site) +
(1 + treatment.replicates)
iter <- 20000
chains <- 4
priors <- c(prior(normal(0, 10), class = "b"))
mod <- brm(formula = form, data = degrad, prior = priors, iter = iter,
chains = 4, cores = 4, control = list(adapt_delta = 0.99, max_treedepth = 15),
file = "./data/processed/regression_model.RDS", file_refit = "always")
form_int <- pc1 ~ habitat + frequency * habitat + frequency.modulation *
habitat + amplitude.modulation * habitat + cicadas * habitat +
speak.height * mic.height * habitat + duration * habitat + (1 |
site) + (1 + treatment.replicates)
mod_int <- brm(formula = form_int, data = degrad, prior = priors,
iter = iter, chains = 4, cores = 4, control = list(adapt_delta = 0.99,
max_treedepth = 15), file = "./data/processed/regression_model_habitat_interactions.RDS",
file_refit = "always")Global model
mod <- readRDS("./data/processed/regression_model.RDS")
summary_brm_model(mod, plot.area.prop = 0.8, grep.filter.out = "treatment.replicates")| b_prior | sd_prior | iterations | chains | thinning | warmup | diverg_transitions | rhats > 10.5 | |
|---|---|---|---|---|---|---|---|---|
| 1 | normal(0, 10) | student_t(3, 0, 2.5) | 20000 | 4 | 1 | 10000 | 0 | 0 |
| Estimate | Rhat | Bulk_ESS | CI_low | CI_high | |
|---|---|---|---|---|---|
| habitatHC | 0.635 | 1.001 | 5086.016 | 0.573 | 0.698 |
| frequency | 0.178 | 1.001 | 6341.744 | 0.170 | 0.186 |
| frequency.modulationfm | -1.628 | 1.005 | 910.539 | -4.382 | 1.037 |
| amplitude.modulationam | 0.405 | 1.01 | 655.352 | -2.239 | 2.967 |
| cicadashigh | 0.203 | 1.001 | 5795.884 | 0.122 | 0.284 |
| speak.height | 0.100 | 1.001 | 3479.690 | -0.149 | 0.356 |
| mic.height | -0.362 | 1.002 | 3459.015 | -0.611 | -0.104 |
| durationshort | -0.376 | 1.01 | 758.975 | -2.991 | 2.172 |
| speak.height:mic.height | 0.042 | 1.002 | 3436.050 | -0.058 | 0.140 |
# conditions <- data.frame(sc_temp = c(`Low temperature` = -1,
# `Mean temperature` = 0, `High temperature` = 1))
plot(conditional_effects(mod))
Habitat interaction model
mod_int <- readRDS("./data/processed/regression_model_habitat_interactions.RDS")
summary_brm_model(mod_int, plot.area.prop = 0.8, grep.filter.out = "treatment.replicates")| b_prior | sd_prior | iterations | chains | thinning | warmup | diverg_transitions | rhats > 10.5 | |
|---|---|---|---|---|---|---|---|---|
| 1 | normal(0, 10) | student_t(3, 0, 2.5) | 20000 | 4 | 1 | 10000 | 1 | 0 |
| Estimate | Rhat | Bulk_ESS | CI_low | CI_high | |
|---|---|---|---|---|---|
| habitatHC | -0.479 | 1 | 10255.808 | -1.787 | 0.805 |
| frequency | 0.063 | 1 | 20378.425 | 0.052 | 0.075 |
| frequency.modulationfm | -1.543 | 1.002 | 2826.835 | -4.235 | 1.120 |
| amplitude.modulationam | 0.232 | 1.003 | 2922.395 | -2.448 | 2.939 |
| cicadashigh | 0.376 | 1 | 19727.747 | 0.263 | 0.489 |
| speak.height | 0.074 | 1 | 10786.261 | -0.278 | 0.423 |
| mic.height | -0.461 | 1.001 | 10843.798 | -0.811 | -0.112 |
| durationshort | -0.345 | 1.001 | 2839.731 | -3.039 | 2.246 |
| habitatHC:frequency | 0.230 | 1 | 20416.332 | 0.214 | 0.247 |
| habitatHC:frequency.modulationfm | 0.102 | 1 | 24695.706 | 0.008 | 0.196 |
| habitatHC:amplitude.modulationam | 0.107 | 1 | 25169.311 | 0.011 | 0.202 |
| habitatHC:cicadashigh | -0.266 | 1 | 19772.312 | -0.390 | -0.143 |
| speak.height:mic.height | 0.093 | 1.001 | 10750.719 | -0.044 | 0.230 |
| habitatHC:speak.height | 0.030 | 1 | 10132.742 | -0.474 | 0.544 |
| habitatHC:mic.height | 0.174 | 1 | 10224.883 | -0.326 | 0.684 |
| habitatHC:durationshort | 0.030 | 1 | 24811.484 | -0.065 | 0.126 |
| habitatHC:speak.height:mic.height | -0.091 | 1.001 | 10067.591 | -0.292 | 0.105 |
Report overview
Takeaways
Session information
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
##
## locale:
## [1] LC_CTYPE=pt_BR.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=es_CR.UTF-8 LC_COLLATE=pt_BR.UTF-8
## [5] LC_MONETARY=es_CR.UTF-8 LC_MESSAGES=pt_BR.UTF-8
## [7] LC_PAPER=es_CR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] kableExtra_1.3.1 ggplot2_3.3.6 viridis_0.6.2
## [4] viridisLite_0.4.0 brms_2.16.3 Rcpp_1.0.9
## [7] Rraven_1.0.13 baRulho_1.0.7 xaringanExtra_0.7.0
## [10] rprojroot_2.0.2 formatR_1.11 warbleR_1.1.28
## [13] NatureSounds_1.0.4 knitr_1.39 seewave_2.2.0
## [16] tuneR_1.4.0
##
## loaded via a namespace (and not attached):
## [1] uuid_0.1-4 backports_1.3.0 plyr_1.8.6
## [4] igraph_1.2.6 splines_4.1.0 crosstalk_1.1.1
## [7] TH.data_1.1-0 rstantools_2.1.1 inline_0.3.19
## [10] digest_0.6.29 htmltools_0.5.3 rsconnect_0.8.27
## [13] fansi_1.0.3 magrittr_2.0.3 checkmate_2.0.0
## [16] RcppParallel_5.1.4 matrixStats_0.61.0 xts_0.12.1
## [19] sandwich_3.0-1 prettyunits_1.1.1 colorspace_2.0-3
## [22] rvest_1.0.0 signal_0.7-7 ggdist_3.0.0
## [25] xfun_0.31 dplyr_1.0.7 callr_3.7.0
## [28] crayon_1.5.1 RCurl_1.98-1.7 jsonlite_1.8.0
## [31] lme4_1.1-27.1 survival_3.2-11 zoo_1.8-9
## [34] glue_1.6.2 gtable_0.3.0 emmeans_1.8.1-1
## [37] webshot_0.5.2 V8_3.4.2 distributional_0.2.2
## [40] pkgbuild_1.2.0 rstan_2.21.2 abind_1.4-5
## [43] scales_1.2.0 mvtnorm_1.1-2 DBI_1.1.1
## [46] miniUI_0.1.1.1 dtw_1.22-3 xtable_1.8-4
## [49] HDInterval_0.2.2 diffobj_0.3.4 proxy_0.4-27
## [52] stats4_4.1.0 StanHeaders_2.21.0-7 DT_0.18
## [55] httr_1.4.2 htmlwidgets_1.5.3 threejs_0.3.3
## [58] posterior_1.1.0 ellipsis_0.3.2 pkgconfig_2.0.3
## [61] loo_2.4.1.9000 farver_2.1.1 sass_0.4.0
## [64] utf8_1.2.2 labeling_0.4.2 tidyselect_1.1.1
## [67] rlang_1.0.5 reshape2_1.4.4 later_1.2.0
## [70] munsell_0.5.0 tools_4.1.0 cli_3.3.0
## [73] generics_0.1.0 ggridges_0.5.3 evaluate_0.15
## [76] stringr_1.4.0 fastmap_1.1.0 yaml_2.3.5
## [79] processx_3.5.2 purrr_0.3.4 pbapply_1.5-0
## [82] nlme_3.1-152 mime_0.11 projpred_2.0.2
## [85] xml2_1.3.2 compiler_4.1.0 bayesplot_1.8.1
## [88] shinythemes_1.2.0 rstudioapi_0.13 curl_4.3.2
## [91] gamm4_0.2-6 tibble_3.1.8 bslib_0.2.5.1
## [94] stringi_1.7.8 highr_0.9 ps_1.6.0
## [97] Brobdingnag_1.2-6 lattice_0.20-44 Matrix_1.3-4
## [100] nloptr_1.2.2.2 markdown_1.1 shinyjs_2.0.0
## [103] fftw_1.0-7 tensorA_0.36.2 vctrs_0.4.1
## [106] pillar_1.8.0 lifecycle_1.0.1 jquerylib_0.1.4
## [109] bridgesampling_1.1-2 estimability_1.4.1 cowplot_1.1.1
## [112] bitops_1.0-7 httpuv_1.6.1 R6_2.5.1
## [115] promises_1.2.0.1 gridExtra_2.3 codetools_0.2-18
## [118] boot_1.3-28 colourpicker_1.1.0 MASS_7.3-54
## [121] gtools_3.9.2 assertthat_0.2.1 rjson_0.2.21
## [124] withr_2.5.0 shinystan_2.5.0 multcomp_1.4-17
## [127] mgcv_1.8-36 parallel_4.1.0 grid_4.1.0
## [130] coda_0.19-4 minqa_1.2.4 rmarkdown_2.16
## [133] shiny_1.6.0 base64enc_0.1-3 dygraphs_1.1.1.6