Variation in vocal behavior under varying stress levels
Budgie INBRE stress and contact calls
Marcelo Araya-Salas, PhD & Timothy Wright, PhD
2025-07-17
Source code and data found at https://github.com/maRce10/budgie_inbre_stress_vocal_output
Load packages
## add 'developer' to packages to be installed from github
x <- c("remotes", "lubridate", "readxl", "pbapply", "viridis", "ggplot2",
"kableExtra", "knitr", "formatR", "MASS", "sp", "GGally", "brms",
"lme4", "dplyr", "purrr", "forcats", "tidyr", "modelr", "tidybayes",
"cowplot", "ggrepel", "posterior", "ggridges", "maRce10/PhenotypeSpace",
"corrplot", "DT")
source("~/Dropbox/R_package_testing/sketchy/R/load_packages.R")
load_packages(x)
source("~/Dropbox/R_package_testing/brmsish/R/extended_summary.R")
source("~/Dropbox/R_package_testing/brmsish/R/helpers.R")
source("~/Dropbox/R_package_testing/brmsish/R/read_summary.R")
source("~/Dropbox/R_package_testing/brmsish/R/check_rds_fits.R")
source("~/Dropbox/R_package_testing/PhenotypeSpace/R/plot_space.R")
print_kable <- function(x) {
kb <- kable(x, row.names = TRUE, digits = 4, "html")
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover",
"condensed", "responsive"))
scroll_box(kb, width = "100%")
}
data_table <- function(x, digits = 3) {
if (is.matrix(x))
x <- as.data.frame(x)
datatable(x) |>
formatRound(columns = which(sapply(x, is.numeric)), digits = digits)
}Functions and global parameters
opts_knit$set(root.dir = "..")
# set evaluation false
opts_chunk$set(fig.width = 10, fig.height = 6, warning = FALSE, message = FALSE,
tidy = TRUE)
read_excel_df <- function(...) data.frame(read_excel(...))
# for reading months in english format
sl <- Sys.setlocale(locale = "en_US.UTF-8")
standard_error <- function(x) sd(x)/sqrt(length(x))
cols <- viridis(10, alpha = 0.7)
col_pointrange <- cols[7]# read ext sel tab calls
sels <- read.csv("./data/processed/tailored_budgie_calls_sel_tab.csv")
# keep only spectrographic parameters
sels <- sels[, c("sound.files", "selec", "duration", "meanfreq", "sd",
"freq.median", "freq.IQR", "time.IQR", "skew", "kurt", "sp.ent",
"time.ent", "entropy", "meandom", "mindom", "maxdom", "dfrange",
"modindx", "meanpeakf")]
sels$ID <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][1])
sels$month <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][2])
sels$day <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][3])
sels$year <- sapply(sels$sound.files, function(x) strsplit(x, "_")[[1]][4])
sels$date <- paste(sels$day, substr(sels$month, 0, 3), sels$year,
sep = "-")
sels$date <- as.Date(sels$date, format = "%d-%b-%Y")
# acoustic measurements
areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week.RDS")
indiv_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_first_week_by_individual.RDS")
indiv_ovlp$treatment <- factor(indiv_ovlp$treatment, levels = c("Control",
"Medium Stress", "High Stress"))
group_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_group_by_week.RDS")
# group_ovlp <- do.call(rbind, lapply(unique(group_ovlp$ID),
# function(x) { X <- group_ovlp[group_ovlp$ID == x, ]
# X$overlap.to.group <- X$overlap.to.group -
# X$overlap.to.group[X$week ==
# min(as.numeric(as.character(X$week)))] X$distance.to.group <-
# X$distance.to.group - X$distance.to.group[X$week ==
# min(as.numeric(as.character(X$week)))] return(X) }))
group_ovlp$treatment <- factor(group_ovlp$treatment, levels = c("Control",
"Medium Stress", "High Stress"))Counts per individual
df <- as.data.frame(table(sels$ID))
names(df) <- c("ID", "Sample_size")
df <- df[order(df$Sample_size, decreasing = FALSE), ]
kb <- kable(df, row.names = FALSE)
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover",
"condensed", "responsive"), position = "float_left")
print(kb)| ID | Sample_size |
|---|---|
| 132YGMM | 6 |
| 125YGMM | 12 |
| 160YGHM | 16 |
| 310YGHM | 21 |
| 393YGHM | 25 |
| 303YGHM | 37 |
| 398WBHM | 41 |
| 365YGHM | 44 |
| 399YGLM | 46 |
| 300YGHM | 47 |
| 270YGHM | 64 |
| 407YGHM | 97 |
| 386WBHM | 100 |
| 377WWLM | 108 |
| 367WWNM | 119 |
| 371YYLM | 149 |
| 397YGHM | 175 |
| 378YYLM | 195 |
| 404WBHM | 196 |
| 258YGHM | 223 |
| 389WWLM | 230 |
| 262YGHM | 306 |
| 400YGHM | 360 |
| 388YGLM | 377 |
| 327YYLM | 404 |
| 396YBHM | 455 |
| 403WBHM | 566 |
| 356WBHM | 761 |
| 361YGHM | 767 |
| 402YGHM | 849 |
| 395WBHM | 854 |
| 360YGHM | 900 |
| 390WBHM | 935 |
| 405YBHM | 975 |
| 385YBHM | 1256 |
| 394YBHM | 1693 |
Add metadata
metadat <- read_excel_df("./data/raw/INBREStress_MasterDataSheet_14Nov19.xlsx")
# head(metadat)
sels$ID[sels$ID == "125YGMM"] <- "125YGHM"
sels$ID[sels$ID == "394YBHM"] <- "394WBHM"
# setdiff(sels$ID, metadat$Bird.ID) setdiff(metadat$Bird.ID,
# sels$ID)
sels$treatment <- sapply(1:nrow(sels), function(x) {
metadat$Treatment[metadat$Bird.ID == sels$ID[x]][1]
})
sels$treatment.room <- sapply(1:nrow(sels), function(x) {
metadat$Treatment.Room[metadat$Bird.ID == sels$ID[x]][1]
})
sels$round <- sapply(1:nrow(sels), function(x) {
metadat$Round[metadat$Bird.ID == sels$ID[x]][1]
})
sels$source.room <- sapply(1:nrow(sels), function(x) {
metadat$Source.Room[metadat$Bird.ID == sels$ID[x]][1]
})
sels$record.group <- sapply(1:nrow(sels), function(x) {
metadat$Record.Group[metadat$Bird.ID == sels$ID[x]][1]
})
# add week
out <- lapply(unique(sels$round), function(x) {
Y <- sels[sels$round == x, ]
min_date <- min(Y$date)
week_limits <- min_date + seq(0, 100, by = 7)
Y$week <- NA
for (i in 2:length(week_limits)) Y$week[Y$date >= week_limits[i -
1] & Y$date < week_limits[i]] <- i - 1
return(Y)
})
sels <- do.call(rbind, out)
sels$cort.baseline <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.CORT.Baseline[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
sels$cort.stress <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.CORT.Stress[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
sels$stress.response <- sels$cort.stress #- sels$cort.baseline
sels$weight <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
sels$breath.count <- sapply(1:nrow(sels), function(x) {
if (sels$week[x] == 1)
out <- metadat$D3.Breath.Count[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 2)
out <- metadat$D7.Breath.Count[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 3)
out <- metadat$D14.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 4)
out <- metadat$D21.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
if (sels$week[x] == 5)
out <- metadat$D28.Bird.Weight..g.[metadat$Bird.ID == sels$ID[x]][1]
return(out)
})
# remove week 5
sels <- sels[sels$week != 5, ]
foxp2 <- read.csv("./data/raw/INBREmale_IHCcounts_9Jun25.csv")
foxp2$Treatment[grep("HIGH", foxp2$Treatment)] <- "High stress"
foxp2$Treatment[grep("Control", foxp2$Treatment, ignore.case = T)] <- "Control"
foxp2$Treatment[grep("MED", foxp2$Treatment)] <- "Medium stress"
foxp2$Treatment <- factor(foxp2$Treatment, levels = c("Control", "Medium stress",
"High stress"))Correlation matrices
agg_dat <- aggregate(selec ~ ID + week, data = sels, length)
# compare to week 1 agg_dat$call.count <-
# sapply(1:nrow(agg_dat), function(x) { baseline <-
# agg_dat$selec[agg_dat$week == 1 & agg_dat$ID == agg_dat$ID[x]]
# if (length(baseline) > 0) change <- agg_dat$selec[x] -
# baseline else change <- agg_dat$selec[x] return(change) } )
# without comparing to week 1
agg_dat$call.count <- sapply(1:nrow(agg_dat), function(x) agg_dat$selec[x])
agg_dat$selec <- NULL
agg_dat$baseline.CORT <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$cort.baseline[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$cort.baseline[sels$week == agg_dat$week[x] & sels$ID ==
agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$stress.response <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$stress.response[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$stress.response[sels$week == agg_dat$week[x] &
sels$ID == agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$stress.CORT <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$cort.stress[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$cort.stress[sels$week == agg_dat$week[x] & sels$ID ==
agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$weight <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$weight[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$weight[sels$week == agg_dat$week[x] & sels$ID ==
agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$breath.rate <- sapply(1:nrow(agg_dat), function(x) {
baseline <- sels$breath.count[sels$week == 1 & sels$ID == agg_dat$ID[x]]
current <- sels$breath.count[sels$week == agg_dat$week[x] & sels$ID ==
agg_dat$ID[x]]
if (length(baseline) > 0 & length(current) > 0)
change <- mean(current) - mean(baseline) else change <- NA
return(change)
})
agg_dat$acoustic.diversity <- sapply(1:nrow(agg_dat), function(x) {
area <- areas_by_week$raref.area[areas_by_week$ID == agg_dat$ID[x] &
areas_by_week$week == agg_dat$week[x]]
if (length(area) < 1)
area <- NA
return(area)
})
agg_dat$acoustic.distance <- sapply(1:nrow(agg_dat), function(x) {
distance <- indiv_ovlp$distance.to.first.week[indiv_ovlp$ID ==
agg_dat$ID[x] & indiv_ovlp$week == agg_dat$week[x]]
if (length(distance) < 1)
distance <- NA
return(distance)
})
agg_dat$acustic.plasticity <- sapply(1:nrow(agg_dat), function(x) {
overlap <- indiv_ovlp$overlap.to.first.week[indiv_ovlp$ID == agg_dat$ID[x] &
indiv_ovlp$week == agg_dat$week[x]]
plasticity <- 1 - overlap
if (length(plasticity) < 1)
plasticity <- NA
return(plasticity)
})
agg_dat$acoustic.convergence <- sapply(1:nrow(agg_dat), function(x) {
overlap <- group_ovlp$overlap.to.group[group_ovlp$ID == agg_dat$ID[x] &
group_ovlp$week == agg_dat$week[x]]
if (length(overlap) < 1)
overlap <- NA
return(overlap)
})
agg_dat$treatment <- sapply(1:nrow(agg_dat), function(x) unique(sels$treatment[sels$ID ==
agg_dat$ID[x]]))
agg_dat$round <- sapply(1:nrow(agg_dat), function(x) unique(sels$round[sels$ID ==
agg_dat$ID[x]]))
agg_dat$foxp2 <- sapply(1:nrow(agg_dat), function(x) {
if (agg_dat$week[x] == 4)
fp2 <- foxp2$FoxP2.Counts.MMSt.Striatum[foxp2$Bird.ID == agg_dat$ID[x]][1] else fp2 <- NA
# print(length(fp2))
if (length(fp2) == 0)
fp2 <- NA
return(fp2)
})Day 3
agg_dat$...foxp2 <- sapply(agg_dat$ID, function(x) {
fp2 <- foxp2$FoxP2.Counts.MMSt.Striatum[foxp2$Bird.ID == x][1]
if (length(fp2) == 0)
fp2 <- NA
return(fp2)
})
day3_agg <- agg_dat[agg_dat$week == 3, c("ID", "breath.rate", "weight",
"baseline.CORT", "stress.CORT", "call.count", "acustic.plasticity",
"acoustic.convergence", "acoustic.diversity", "...foxp2")]
print("Raw data:")## [1] "Raw data:"data_table(day3_agg, 4)day3_agg$ID <- NULL
corr_mat <- cor(day3_agg, use = "pairwise.complete.obs")
pvalue_mat <- corr_mat
pvalue_mat <- pvalue_mat * NA
for (i in colnames(pvalue_mat)) {
for (e in rownames(pvalue_mat)) {
if (sum(complete.cases(day3_agg[, i], day3_agg[, e])) > 2)
ct <- cor.test(x = day3_agg[, i], y = day3_agg[, e], method = "pearson",
use = "pairwise.complete.obs")$p.value else ct <- NA
pvalue_mat[i, e] <- ct
pvalue_mat[e, i] <- ct
}
}
n_mat <- corr_mat
n_mat <- n_mat * NA
for (i in colnames(pvalue_mat)) {
for (e in rownames(pvalue_mat)) {
n_mat[i, e] <- sum(complete.cases(day3_agg[, i], day3_agg[,
e]))
}
}
# write.csv(pm, './output/parameters_correlation_pvalues.csv',
# row.names = TRUE) write.csv(pm2, './output/sample_sizes.csv',
# row.names = TRUE)
# make color scale
colorrange <- colorRampPalette(c("skyblue", "white", "tomato"))(10)
colorrange2 <- colorRampPalette(c("skyblue", "gray", "tomato"))(10)
rownames(corr_mat) <- colnames(corr_mat) <- names(day3_agg) <- c("breath\nrate",
"weight", "baseline\nCORT", "stress\nCORT", "vocal\noutput", "vocal\nplasticity",
"vocal\nconver-\ngence", "vocal\ndiversity", "FOXP2")
# jpeg('./output/parameters_correlation_matrix.jpg', width =
# 1800, height = 1800, res = 300)
corrplot.mixed(corr_mat, p.mat = pvalue_mat, sig.level = 0.05, insig = "pch",
pch = "x", pch.col = "gray", pch.cex = 0.9, upper = "ellipse",
lower = "number", tl.col = "black", tl.srt = 45, number.cex = 0.7,
tl.cex = 0.7, order = "original", mar = c(1, 1, 1, 1), upper.col = colorrange,
lower.col = colorrange2)
# dev.off()
print("p values:")## [1] "p values:"data_table(pvalue_mat, 4)print("sample sizes:")## [1] "sample sizes:"data_table(n_mat, 0)Day 4
day4_agg <- agg_dat[agg_dat$week == 4, c("ID", "breath.rate", "weight",
"baseline.CORT", "stress.CORT", "call.count", "acustic.plasticity",
"acoustic.convergence", "acoustic.diversity", "foxp2")]
print("Raw data:")[1] “Raw data:”
data_table(day4_agg)day4_agg$ID <- NULL
corr_mat <- cor(day4_agg, use = "pairwise.complete.obs")
pvalue_mat <- corr_mat
pvalue_mat <- pvalue_mat * NA
for (i in colnames(pvalue_mat)) {
for (e in rownames(pvalue_mat)) {
ct <- cor.test(x = day4_agg[, i], y = day4_agg[, e], method = "pearson",
use = "pairwise.complete.obs")$p.value
pvalue_mat[i, e] <- ct
pvalue_mat[e, i] <- ct
}
}
n_mat <- corr_mat
n_mat <- n_mat * NA
for (i in colnames(pvalue_mat)) {
for (e in rownames(pvalue_mat)) {
n_mat[i, e] <- sum(complete.cases(day4_agg[, i], day4_agg[,
e]))
}
}
# write.csv(pm, './output/parameters_correlation_pvalues.csv',
# row.names = TRUE) write.csv(pm2, './output/sample_sizes.csv',
# row.names = TRUE)
# make color scale
colorrange <- colorRampPalette(c("skyblue", "white", "tomato"))(10)
colorrange2 <- colorRampPalette(c("skyblue", "gray", "tomato"))(10)
rownames(corr_mat) <- colnames(corr_mat) <- names(day4_agg) <- c("breath\nrate",
"weight", "baseline\nCORT", "stress\nCORT", "vocal\noutput", "vocal\nplasticity",
"vocal\nconver-\ngence", "vocal\ndiversity", "FOXP2")
# jpeg('./output/parameters_correlation_matrix.jpg', width =
# 1800, height = 1800, res = 300)
corrplot.mixed(corr_mat, p.mat = pvalue_mat, sig.level = 0.05, insig = "pch",
pch = "x", pch.col = "gray", pch.cex = 0.9, upper = "ellipse",
lower = "number", tl.col = "black", tl.srt = 45, number.cex = 0.7,
tl.cex = 0.7, order = "original", mar = c(1, 1, 1, 1), upper.col = colorrange,
lower.col = colorrange2)
# dev.off()
print("p values:")[1] “p values:”
data_table(pvalue_mat, 4)print("sample sizes:")[1] “sample sizes:”
data_table(n_mat, 0)Day 3 and 4
Averaged across the 2 weeks, excluding NAs
day34_agg <- agg_dat[agg_dat$week %in% 3:4, c("ID", "breath.rate",
"weight", "baseline.CORT", "stress.CORT", "call.count", "acustic.plasticity",
"acoustic.convergence", "acoustic.diversity", "...foxp2")]
print("Raw data:")[1] “Raw data:”
data_table(day34_agg)day34_agg <- aggregate(day34_agg[, -1], by = list(ID = day34_agg$ID),
FUN = mean, na.rm = TRUE)
print("Averaged data:")[1] “Averaged data:”
data_table(day34_agg)day34_agg$ID <- NULL
corr_mat <- cor(day34_agg, use = "pairwise.complete.obs")
pvalue_mat <- corr_mat
pvalue_mat <- pvalue_mat * NA
for (i in colnames(pvalue_mat)) {
for (e in rownames(pvalue_mat)) {
ct <- cor.test(x = day34_agg[, i], y = day34_agg[, e], method = "pearson",
use = "pairwise.complete.obs")$p.value
pvalue_mat[i, e] <- ct
pvalue_mat[e, i] <- ct
}
}
n_mat <- corr_mat
n_mat <- n_mat * NA
for (i in colnames(pvalue_mat)) {
for (e in rownames(pvalue_mat)) {
n_mat[i, e] <- sum(complete.cases(day34_agg[, i], day34_agg[,
e]))
}
}
# write.csv(pm, './output/parameters_correlation_pvalues.csv',
# row.names = TRUE) write.csv(pm2, './output/sample_sizes.csv',
# row.names = TRUE)
# make color scale
colorrange <- colorRampPalette(c("skyblue", "white", "tomato"))(10)
colorrange2 <- colorRampPalette(c("skyblue", "gray", "tomato"))(10)
rownames(corr_mat) <- colnames(corr_mat) <- names(day34_agg) <- c("breath\nrate",
"weight", "baseline\nCORT", "stress\nCORT", "vocal\noutput", "vocal\nplasticity",
"vocal\nconver-\ngence", "vocal\ndiversity", "FOXP2")
# jpeg('./output/parameters_correlation_matrix.jpg', width =
# 1800, height = 1800, res = 300)
corrplot.mixed(corr_mat, p.mat = pvalue_mat, sig.level = 0.05, insig = "pch",
pch = "x", pch.col = "gray", pch.cex = 0.9, upper = "ellipse",
lower = "number", tl.col = "black", tl.srt = 45, number.cex = 0.7,
tl.cex = 0.7, order = "original", mar = c(1, 1, 1, 1), upper.col = colorrange,
lower.col = colorrange2)
# dev.off()
print("p values:")[1] “p values:”
data_table(pvalue_mat, 4)print("sample sizes:")[1] “sample sizes:”
data_table(n_mat, 0)Physiological parameters
Barplot and effect sizes graph
physio_model_list <- list.files("./data/processed/updated_regressions",
pattern = "physio-", full.names = TRUE)
# read all physio models
physio_models <- lapply(physio_model_list, readRDS)
names(physio_models) <- gsub(".rds", "", basename(physio_model_list))
breath.count <- stack(metadat[, c("D3.Breath.Count", "D7.Breath.Count",
"D14.Breath.Count", "D21.Breath.Count", "D28.Breath.Count")])
weight <- stack(metadat[, c("D3.Bird.Weight..g.", "D7.Bird.Weight..g.",
"D14.Bird.Weight..g.", "D21.Bird.Weight..g.", "D28.Bird.Weight..g.")])
cort.stress <- stack(metadat[, c("D3.CORT.Stress", "D7.CORT.Stress",
"D14.CORT.Stress", "D21.CORT.Stress", "D28.CORT.Stress")])
cort.baseline <- stack(metadat[, c("D3.CORT.Baseline", "D7.CORT.Baseline",
"D14.CORT.Baseline", "D21.CORT.Baseline", "D28.CORT.Baseline")])
stress <- data.frame(metadat[, c("Bird.ID", "Treatment", "Round",
"Treatment.Room")], week = breath.count$ind, breath.count = breath.count$values,
weight = weight$values, cort.stress = cort.stress$values, cort.baseline = cort.baseline$values,
stress.response = cort.stress$values - cort.baseline$values)
# head(stress)
stress$week <- factor(sapply(strsplit(as.character(stress$week), "\\."),
"[[", 1), levels = c("D3", "D7", "D14", "D21", "D28"))
stress$day <- as.numeric(gsub("D", "", as.character(stress$week)))
stress$round <- paste("Round", stress$Round)
stress$Treatment <- gsub("Medium", "Med.", stress$Treatment)
stress$treatment <- factor(stress$Treatment, levels = c("Control",
"Med. Stress", "High Stress"))
# remove 5th week
stress <- stress[stress$week != "D28", ]
stress_l <- lapply(stress$Bird.ID, function(x) {
X <- stress[stress$Bird.ID == x, ]
X$breath.count <- X$breath.count - X$breath.count[X$week == "D3"]
X$weight <- X$weight - X$weight[X$week == "D3"]
X$cort.stress <- X$cort.stress - X$cort.stress[X$week == "D3"]
X$cort.baseline <- X$cort.baseline - X$cort.baseline[X$week ==
"D3"]
X$stress.response <- X$stress.response - X$stress.response[X$week ==
"D3"]
return(X)
})
stress <- do.call(rbind, stress_l)
agg_stress <- aggregate(cbind(breath.count, weight, stress.response,
cort.baseline) ~ week + treatment, stress, mean)
agg_stress_se <- aggregate(cbind(breath.count, weight, stress.response,
cort.baseline) ~ week + treatment, stress, standard_error)
names(agg_stress_se) <- paste(names(agg_stress_se), ".se", sep = "")
agg_stress <- cbind(agg_stress, agg_stress_se[, 3:6])
agg_stress$Week <- 1:4
bs <- if (isTRUE(getOption("knitr.in.progress"))) 10 else 20
gg_breath.count <- ggplot(data = agg_stress, aes(x = Week, y = breath.count,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = breath.count -
breath.count.se, ymax = breath.count + breath.count.se), width = 0.1) +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
ncol = 3, scale = "fixed") + labs(y = "Mean change in breath\nrate (breaths/min)",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_weight <- ggplot(data = agg_stress, aes(x = Week, y = weight, fill = treatment)) +
geom_bar(stat = "identity") + geom_errorbar(aes(ymin = weight -
weight.se, ymax = weight + weight.se), width = 0.1) + scale_fill_viridis_d(begin = 0.1,
end = 0.9) + facet_wrap(~treatment, ncol = 3, scale = "fixed") +
labs(y = "Mean change in \nweight (grams)", x = "Week") + theme_classic(base_size = bs) +
theme(legend.position = "none")
gg_cort.baseline <- ggplot(data = agg_stress, aes(x = Week, y = cort.baseline,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = cort.baseline -
cort.baseline.se, ymax = cort.baseline + cort.baseline.se), width = 0.1) +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
ncol = 3, scale = "fixed") + labs(y = "Mean change in\nbaseline CORT (ng/mL)",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_stress.response <- ggplot(data = agg_stress, aes(x = Week, y = stress.response,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = stress.response -
stress.response.se, ymax = stress.response + stress.response.se),
width = 0.1) + scale_fill_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Mean change in stress\nresponse CORT (ng/mL)",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_coeffs_physio <- lapply(physio_models, function(x) {
vars <- grep("b_", posterior::variables(x), value = TRUE)
draws <- posterior::as_draws_array(x, variable = vars)
coef_table <- draw_summary(draws, variables = vars, probs = c(0.025,
0.975), robust = TRUE, spread.type = "HPDI")
coef_table$predictor <- rownames(coef_table)
coef_table$predictor <- gsub("b_treatment|b_", "", coef_table$predictor)
coef_table$predictor <- gsub("Stress", " stress", coef_table$predictor)
coef_table$predictor <- gsub("week", "Week", coef_table$predictor)
coef_table <- coef_table[coef_table$predictor != "Intercept",
]
gg_coef <- ggplot2::ggplot(data = coef_table, aes(x = Estimate,
y = predictor)) + geom_vline(xintercept = 0, lty = 2) + ggplot2::geom_point(size = 3,
col = col_pointrange) + ggplot2::geom_errorbar(ggplot2::aes(xmin = `l-95% CI`,
xmax = `u-95% CI`), width = 0, col = col_pointrange) + ggplot2::theme_classic(base_size = bs) +
ggplot2::theme(axis.ticks.length = ggplot2::unit(0, "pt"),
plot.margin = ggplot2::margin(0, 0, 0, 0, "pt"), legend.position = "none",
strip.background = ggplot2::element_blank(), strip.text = ggplot2::element_blank()) +
ggplot2::labs(x = "Effect size", y = "")
return(gg_coef)
})
cowplot::plot_grid(gg_breath.count, gg_weight, gg_cort.baseline, gg_stress.response,
gg_coeffs_physio[[grep("breath", names(gg_coeffs_physio))]] +
theme_classic(base_size = bs), gg_coeffs_physio[[grep("weight",
names(gg_coeffs_physio))]] + theme_classic(base_size = bs),
gg_coeffs_physio[[grep("baseline", names(gg_coeffs_physio))]] +
theme_classic(base_size = bs), gg_coeffs_physio[[grep("response",
names(gg_coeffs_physio))]] + theme_classic(base_size = bs),
nrow = 2, rel_heights = c(1.8, 1))
# try bs = 20 for saving plots
if (!isTRUE(getOption("knitr.in.progress"))) {
cowplot::ggsave2(filename = "./output/bar_graphs_and_estimates_physiology_70dpi.jpeg",
width = 23, height = 9, device = grDevices::png)
cowplot::ggsave2(filename = "./output/bar_graphs_and_estimates_physiology_300dpi.jpeg",
dpi = 300, width = 23, height = 9, device = grDevices::png)
}Stats
Models: Predicted physio measure ~ treatment + week (continuous) + IndRandom
Variables (Difference from Week 1): weight, BR, baseline CORT, Stress CORT, Stress Response
responses <- c("baseline.CORT", "stress.response", "stress.CORT",
"weight", "breath.rate")
predictors <- c("~ treatment + week + (1|ID) + (1|round)")
formulas <- expand.grid(responses = responses, predictors = predictors,
stringsAsFactors = FALSE)
vars_to_scale <- c(responses, "week")
# remove week 1
sub_agg_dat <- agg_dat[agg_dat$week != 1, ]
for (i in vars_to_scale) sub_agg_dat[, vars_to_scale] <- scale(sub_agg_dat[,
vars_to_scale])
physio_models <- lapply(1:nrow(formulas), function(x) {
sub_dat <- sub_agg_dat[!is.na(sub_agg_dat[names(sub_agg_dat) ==
formulas$responses[x]]), ]
sub_dat
# replace ~ by 'by' in the formula
mod_name <- paste(formulas$responses[x], gsub("~", "by", formulas$predictors[x]))
mod <- brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
iter = 1e+05, silent = 2, family = student(), data = sub_dat,
control = list(adapt_delta = 0.9, max_treedepth = 15), chains = 4,
cores = 4, file = paste0("./data/processed/updated_regressions/physio-",
mod_name), file_refit = "always", prior = c(prior(normal(0,
5), "b"), prior(normal(0, 10), "Intercept"), prior(student_t(3,
0, 10), "sd"), prior(student_t(3, 0, 10), "sigma")))
return(mod)
})physio_models <- list.files("./data/processed/updated_regressions",
pattern = "physio-", full.names = TRUE)
for (x in 1:length(physio_models)) extended_summary(read.file = physio_models[x],
gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
remove.intercepts = TRUE)physio-acustic.plasticity by treatment + week + (1|ID) + (1|round)
| formula | family | priors | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | breath.rate ~ treatment + week + (1 | ID) + (1 | round) | student (identity) | b-normal(0, 5) Intercept-normal(0, 10) nu-gamma(2, 0.1) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) | 1e+05 | 4 | 1 | 50000 | 0 (0%) | 0 | 282877.3 | 144723.3 | 1702562168 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | 1.563 | -7.621 | 10.703 | 1 | 367754.7 | 147396.5 |
| MediumStress | -0.199 | -9.468 | 9.053 | 1 | 368223.6 | 144723.3 |
| week | -16.139 | -25.424 | -6.599 | 1 | 282877.3 | 162601.3 |
physio-baseline.CORT by treatment + week + (1|ID) + (1|round)
| formula | family | priors | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | baseline.CORT ~ treatment + week + (1 | ID) + (1 | round) | student (identity) | b-normal(0, 5) Intercept-normal(0, 10) nu-gamma(2, 0.1) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) | 1e+05 | 4 | 1 | 50000 | 6255 (0.031%) | 0 | 8017.135 | 31502.5 | 662141751 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | 0.874 | 0.110 | 1.637 | 1.001 | 21696.113 | 35836.92 |
| MediumStress | 0.194 | -0.620 | 0.995 | 1 | 25210.214 | 44566.31 |
| week | -0.010 | -0.139 | 0.119 | 1.001 | 8017.135 | 31502.50 |
physio-breath.rate by treatment + week + (1|ID) + (1|round)
| formula | family | priors | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | breath.rate ~ treatment + week + (1 | ID) + (1 | round) | student (identity) | b-normal(0, 5) Intercept-normal(0, 10) nu-gamma(2, 0.1) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) | 1e+05 | 4 | 1 | 50000 | 5808 (0.029%) | 0 | 4237.78 | 22921.74 | 1411987597 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | 0.228 | -0.167 | 0.620 | 1.001 | 6084.785 | 22921.74 |
| MediumStress | 0.115 | -0.310 | 0.544 | 1 | 15827.830 | 87358.76 |
| week | -0.768 | -0.909 | -0.625 | 1.001 | 4237.780 | 35802.47 |
physio-stress.CORT by treatment + week + (1|ID) + (1|round)
| formula | family | priors | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | stress.CORT ~ treatment + week + (1 | ID) + (1 | round) | student (identity) | b-normal(0, 5) Intercept-normal(0, 10) nu-gamma(2, 0.1) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) | 1e+05 | 4 | 1 | 50000 | 1796 (0.009%) | 0 | 22431.33 | 5699.755 | 1854922145 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | -0.259 | -1.044 | 0.529 | 1 | 22431.33 | 5699.755 |
| MediumStress | 0.037 | -0.852 | 0.918 | 1 | 37367.15 | 70541.436 |
| week | -0.065 | -0.170 | 0.039 | 1 | 44228.85 | 45826.122 |
physio-stress.response by treatment + week + (1|ID) + (1|round)
| formula | family | priors | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | stress.response ~ treatment + week + (1 | ID) + (1 | round) | student (identity) | b-normal(0, 5) Intercept-normal(0, 10) nu-gamma(2, 0.1) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) | 1e+05 | 4 | 1 | 50000 | 4766 (0.024%) | 0 | 4017.288 | 4369.013 | 1337746393 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | -0.244 | -1.037 | 0.553 | 1.001 | 8303.060 | 6123.450 |
| MediumStress | 0.048 | -0.840 | 0.917 | 1.002 | 4017.288 | 4369.013 |
| week | -0.064 | -0.168 | 0.038 | 1.001 | 4332.610 | 8721.908 |
physio-weight by treatment + week + (1|ID) + (1|round)
| formula | family | priors | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | weight ~ treatment + week + (1 | ID) + (1 | round) | student (identity) | b-normal(0, 5) Intercept-normal(0, 10) nu-gamma(2, 0.1) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) | 1e+05 | 4 | 1 | 50000 | 4876 (0.024%) | 0 | 1595.299 | 536.633 | 1479872948 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | -0.400 | -1.185 | 0.399 | 1.002 | 4731.667 | 36350.344 |
| MediumStress | -0.166 | -1.031 | 0.706 | 1.001 | 22880.693 | 69798.743 |
| week | -0.085 | -0.210 | 0.039 | 1.003 | 1595.299 | 536.633 |
Takeaways
Breath rate decreases gradually across time after the first week
Stress response is higher in “high stress” birds compared to first week
Acoustic space projection
t-SNE
tsne <- readRDS("./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")
Y <- as.data.frame(tsne$Y)
names(Y) <- c("TSNE1", "TSNE2")
sels <- data.frame(sels, Y)
ids <- c("371YYLM", "395WBHM")
mat <- matrix(c(11:14, 1:4, 5:8, 22, 9, 9, 18),
ncol = 4,
byrow = TRUE)
mat <- cbind(c(23, 10, 10, 20), mat, c(21, 15, 16, 19))
mat <- rbind(rep(17, 6), mat)
if (!isTRUE(getOption('knitr.in.progress')))
jpeg("./output/acoustic_space_plot.jpg", width = 1800, height = 1800, res = 300)
layout(mat,
widths = c(0.4, 1, 1, 1, 1, 0.2),
heights = c(2, 0.2, 1, 1, 0.4))
par(mar = c(0, 0, 0, 0))
panel_count <- 1
colors <- mako(2)
for (x in ids)
for (i in 1:4) {
plot_space(
X = sels,
dimensions = c("TSNE1", "TSNE2"),
indices = which(sels$ID == x &
sels$week == i),
# background.indices = which(
# sels$ID != x &
# sels$week == i &
# sels$record.group == sels$record.group[sels$ID == x][1]
# ),
basecex = 2,
pch = 20,
labels = c(x, "group"),
legend = NULL,
xaxt = if (panel_count %in% mat[4, ])
"s"
else
"n",
yaxt = if (panel_count %in% mat[, 2])
"s"
else
"n",
colors = colors[1:2], point.cex = 1.7, point.alpha = 0.3,
point.colors = c("black", "black")
)
panel_count <- panel_count + 1
}
par(mar = c(0, 0, 0, 0))
plot(
1,
frame.plot = FALSE,
type = "n",
yaxt = "n",
xaxt = "n"
)
text(x = 1, y = 0.75, "TSNE1", cex = 1.2)
par(mar = c(0, 0, 0, 0))
plot(
1,
frame.plot = FALSE,
type = "n",
yaxt = "n",
xaxt = "n"
)
text(
x = 0.8,
y = 1,
"TSNE2",
srt = 90,
cex = 1.2
)
for (u in paste("Week", 1:4)) {
par(mar = c(0, 0, 0, 0))
plot(
1,
frame.plot = FALSE,
type = "n",
yaxt = "n",
xaxt = "n"
)
if (u == "Week 1") {
mtext("b", side=2, line=2.2, at=1, las = 1, cex = 1.5)
}
rect(0, 0, 2, 2, col = colors[2], border = NA)
box()
text(x = 1, y = 1, u, cex = 1.2)
}
for (u in ids) {
par(mar = c(0, 0, 0, 0))
plot(
1,
frame.plot = FALSE,
type = "n",
yaxt = "n",
xaxt = "n"
)
rect(0, 0, 2, 2, col = colors[2], border = NA)
box()
text(
x = 1,
y = 1,
u,
srt = 270,
cex = 1.2
)
}
# set par to default
par(mar = c(2, 4, 1, 2) + 0.1)
# make scatter plot like gg_acou_space using base R plotting functions
plot(sels$TSNE1,
sels$TSNE2,
col = "white",
pch = 20,
cex = 1.2,
xlab = "",
ylab = "TSNE2", cex.lab = 1.2)
# add legend with no margin color
legend(
x = -75,
y = -25,
legend = c( "Control", "Medium Stress", "High Stress"),
col = viridis::viridis(3, alpha = 0.3, direction = 1),
pch = 20,
box.lwd = 0,
pt.cex = 1.8,
cex = 1.2
)
points(sels$TSNE1,
sels$TSNE2,
col = viridis::viridis(3, alpha = 0.3, direction = 1)[as.numeric(as.factor(sels$treatment))],
pch = 20,
cex = 1.2)
mtext("a", side=2, line=2.7, at=65, las = 1, cex = 1.5)
if (!isTRUE(getOption('knitr.in.progress')))
dev.off()
Overlap between treatments
ss <- space_similarity(treatment ~ TSNE1 + TSNE2, data = sels, method = "density.overlap")
print_kable(ss)| group.1 | group.2 | overlap.1in2 | overlap.2in1 | mean.overlap | |
|---|---|---|---|---|---|
| 1 | Control | High Stress | 1.0000 | 0.8080 | 0.9040 |
| 2 | Control | Medium Stress | 1.0000 | 0.8456 | 0.9228 |
| 3 | High Stress | Medium Stress | 0.9354 | 1.0000 | 0.9677 |
Mean overlap: 0.9315006
Range of overlap: 0.9039993 to 0.9677188
Barplot and effect sizes graph
behav_model_list <- list.files("./data/processed/updated_regressions",
pattern = "behav-", full.names = TRUE)
# read all physio models
behav_models <- lapply(behav_model_list, readRDS)
names(behav_models) <- gsub(".rds", "", basename(behav_model_list))
agg_call.count <- aggregate(cbind(call.count, acoustic.convergence) ~
week + treatment, agg_dat, mean)
agg_behav <- aggregate(cbind(acoustic.diversity, acustic.plasticity) ~
week + treatment, agg_dat, mean)
agg_call.count_se <- aggregate(cbind(call.count, acoustic.convergence) ~
week + treatment, agg_dat, standard_error)
agg_behav_se <- aggregate(cbind(acoustic.diversity, acustic.plasticity) ~
week + treatment, agg_dat, standard_error)
agg_behav_se <- merge(agg_call.count_se, agg_behav_se, all = TRUE)
names(agg_behav_se) <- paste(names(agg_behav_se), ".se", sep = "")
agg_behav <- merge(agg_call.count, agg_behav, all = TRUE)
agg_behav <- cbind(agg_behav, agg_behav_se[, 3:6])
bs <- if (isTRUE(getOption("knitr.in.progress"))) 10 else 20
agg_behav$treatment <- gsub("Medium", "Med.", agg_behav$treatment)
agg_behav$treatment <- factor(agg_behav$treatment, levels = c("Control",
"Med. Stress", "High Stress"))
gg_call.count <- ggplot(data = agg_behav, aes(x = week, y = call.count,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = call.count -
call.count.se, ymax = call.count + call.count.se), width = 0.1) +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
ncol = 3, scale = "fixed") + labs(y = "Vocal output", x = "Week") +
theme_classic(base_size = bs) + theme(legend.position = "none")
gg_acoustic.diversity <- ggplot(data = agg_behav, aes(x = week, y = acoustic.diversity,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = acoustic.diversity -
acoustic.diversity.se, ymax = acoustic.diversity + acoustic.diversity.se),
width = 0.1) + scale_fill_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Change in vocal diversity",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_acustic.plasticity <- ggplot(data = agg_behav, aes(x = week, y = acustic.plasticity,
fill = treatment)) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = acustic.plasticity -
acustic.plasticity.se, ymax = acustic.plasticity + acustic.plasticity.se),
width = 0.1) + scale_fill_viridis_d(begin = 0.1, end = 0.9) +
facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Vocal plasticity",
x = "Week") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_acoustic.convergence <- ggplot(data = agg_behav, aes(x = week,
y = acoustic.convergence, fill = treatment)) + geom_bar(stat = "identity") +
geom_errorbar(aes(ymin = acoustic.convergence - acoustic.convergence.se,
ymax = acoustic.convergence + acoustic.convergence.se), width = 0.1) +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment,
ncol = 3, scale = "fixed") + labs(y = "Vocal convergence", x = "Week") +
theme_classic(base_size = bs) + theme(legend.position = "none")
gg_coeffs_behav <- lapply(behav_models, function(x) {
vars <- grep("b_", posterior::variables(x), value = TRUE)
draws <- posterior::as_draws_array(x, variable = vars)
coef_table <- draw_summary(draws, variables = vars, probs = c(0.025,
0.975), robust = TRUE, spread.type = "HPDI")
coef_table$predictor <- rownames(coef_table)
coef_table$predictor <- gsub("b_treatment|b_", "", coef_table$predictor)
coef_table$predictor <- gsub("Stress", " stress", coef_table$predictor)
coef_table$predictor <- gsub("week", "Week", coef_table$predictor)
coef_table <- coef_table[coef_table$predictor != "Intercept",
]
gg_coef <- ggplot2::ggplot(data = coef_table, aes(x = Estimate,
y = predictor)) + geom_vline(xintercept = 0, lty = 2) + ggplot2::geom_point(size = 4,
col = col_pointrange) + ggplot2::geom_errorbar(ggplot2::aes(xmin = `l-95% CI`,
xmax = `u-95% CI`), width = 0, col = col_pointrange) + ggplot2::theme_classic(base_size = bs) +
ggplot2::theme(axis.ticks.length = ggplot2::unit(0, "pt"),
plot.margin = ggplot2::margin(0, 0, 0, 0, "pt"), legend.position = "none",
strip.background = ggplot2::element_blank(), strip.text = ggplot2::element_blank()) +
ggplot2::labs(x = "Effect size", y = "")
return(gg_coef)
})
cowplot::plot_grid(gg_call.count, gg_acoustic.diversity, gg_acustic.plasticity,
gg_acoustic.convergence, gg_coeffs_behav[[grep("count", names(gg_coeffs_behav))]] +
theme_classic(base_size = bs), gg_coeffs_behav[[grep("diversity",
names(gg_coeffs_behav))]] + theme_classic(base_size = bs),
gg_coeffs_behav[[grep("plasticity", names(gg_coeffs_behav))]] +
theme_classic(base_size = bs), gg_coeffs_behav[[grep("convergence",
names(gg_coeffs_behav))]] + theme_classic(base_size = bs),
nrow = 2, rel_heights = c(1.8, 1))
# try bs = 20 for saving plots
if (!isTRUE(getOption("knitr.in.progress"))) {
cowplot::ggsave2(filename = "./output/bar_graphs_and_estimates_behavior_70dpi.jpeg",
width = 23, height = 9, device = grDevices::png)
#
cowplot::ggsave2(filename = "./output/bar_graphs_and_estimates_behavior_300dpi.jpeg",
dpi = 300, width = 23, height = 9, device = grDevices::png)
}Stats
Model: Predicted behavior ~ treatment + week (continuous) + IndRandom
Variables: # calls, Distance moved from self in first week, Overlap to original acoustic space, Match to group repertoire, Maybe overall size of acoustic space
Do as comparison to week one using rarefacted calls and kernel density
responses <- c("call.count", "acoustic.diversity", "acustic.plasticity", "acoustic.convergence")
predictors <- c("~ treatment + week + (1|ID) + (1|round)")
formulas <- expand.grid(responses = responses, predictors = predictors, stringsAsFactors = FALSE)
# vars_to_scale <- c(responses, "week")
# for (i in vars_to_scale) agg_dat[, vars_to_scale] <- scale(agg_dat[, vars_to_scale])
behav_models <- lapply(1:nrow(formulas), function(x) {
sub_dat <- agg_dat[!is.na(agg_dat[names(agg_dat) == formulas$responses[x]]),]
# remove week 1
if (!grepl("count|group", formulas$responses[x]))
sub_dat <- sub_dat[sub_dat$week != 1, ]
mod_name <- paste(formulas$responses[x], gsub("~", "by", formulas$predictors[x]))
mod <- brm(
formula = paste(formulas$responses[x], formulas$predictors[x]),
iter = 100000,
family = if (formulas$responses[x] == "call.count") negbinomial() else Beta(link = "logit"),
silent = 2,
data = sub_dat,
control = list(adapt_delta = 0.9, max_treedepth = 15),
chains = 4,
cores = 4,
file = paste0("./data/processed/updated_regressions/behav-", mod_name),
file_refit = "always",
prior = c(
prior(normal(0, 5), "b"),
prior(normal(0, 10), "Intercept") #,
# prior(student_t(3, 0, 10), "sd"),
# prior(student_t(3, 0, 10), "sigma")
)
)
return(NULL)
})behav_models <- list.files("./data/processed/updated_regressions",
pattern = "behav-", full.names = TRUE)
for (x in 1:length(behav_models)) extended_summary(read.file = behav_models[x],
gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
remove.intercepts = TRUE)behav-acoustic.convergence by treatment + week + (1|ID) + (1|round)
| formula | family | priors | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | acoustic.convergence ~ treatment + week + (1 | ID) + (1 | round) | beta (logit) | b-normal(0, 5) Intercept-normal(0, 10) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 2.5) | 1e+05 | 4 | 1 | 50000 | 2682 (0.013%) | 0 | 11010.11 | 24870.45 | 1543370143 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | 0.284 | -0.380 | 0.970 | 1.001 | 11010.11 | 24870.45 |
| MediumStress | 0.328 | -0.334 | 1.016 | 1.001 | 22753.15 | 53261.45 |
| week | -0.316 | -0.564 | -0.063 | 1 | 24236.72 | 33083.95 |
behav-acoustic.diversity by treatment + week + (1|ID) + (1|round)
| formula | family | priors | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | acoustic.diversity ~ treatment + week + (1 | ID) + (1 | round) | beta (logit) | b-normal(0, 5) Intercept-normal(0, 10) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 2.5) | 1e+05 | 4 | 1 | 50000 | 6528 (0.033%) | 0 | 20296.19 | 34793.85 | 195705932 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | -0.132 | -0.730 | 0.461 | 1 | 20296.19 | 35122.53 |
| MediumStress | -0.558 | -1.160 | 0.049 | 1 | 22357.72 | 45182.51 |
| week | 0.043 | -0.107 | 0.188 | 1 | 35567.91 | 34793.85 |
behav-acustic.plasticity by treatment + week + (1|ID) + (1|round)
| formula | family | priors | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | acustic.plasticity ~ treatment + week + (1 | ID) + (1 | round) | beta (logit) | b-normal(0, 5) Intercept-normal(0, 10) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 2.5) | 1e+05 | 4 | 1 | 50000 | 4613 (0.023%) | 0 | 9587.206 | 25616.83 | 443607903 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | -1.068 | -2.065 | -0.065 | 1.001 | 9587.206 | 31322.73 |
| MediumStress | -0.245 | -1.315 | 0.840 | 1 | 24675.909 | 56131.47 |
| week | 0.202 | 0.045 | 0.361 | 1 | 18446.592 | 25616.83 |
behav-call.count by treatment + week + (1|ID) + (1|round)
| formula | family | priors | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | call.count ~ treatment + week + (1 | ID) + (1 | round) | negbinomial (log) | b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | 1e+05 | 4 | 1 | 50000 | 3227 (0.016%) | 0 | 4923.233 | 1834.92 | 498436169 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| HighStress | 1.341 | 0.325 | 2.351 | 1 | 31246.848 | 8083.761 |
| MediumStress | 0.872 | -0.144 | 1.893 | 1 | 55318.238 | 75546.858 |
| week | -0.220 | -0.440 | -0.004 | 1.002 | 4923.233 | 1834.920 |
Takeaways
Lower vocal output over time
Higher vocal output in “high stress” birds compared to control
Lower acoustic plasticity to themselves in week 1 for “high stress” birds compared to control
Increase in acoustic plasticity over time
FOXP2
FoxP2 positive cells in MMST:VSP
data_table(foxp2)Stats
foxp2_mod <- brm(formula = FoxP2.Counts.MMSt.Striatum ~ Treatment,
iter = 1e+05, family = Gamma(link = "log"), silent = 2, data = foxp2,
control = list(adapt_delta = 0.9, max_treedepth = 15), chains = 4,
cores = 4, file = "./data/processed/updated_regressions/foxp2_by_treatment_model",
file_refit = "on_change", prior = c(prior(normal(0, 5), "b"),
prior(normal(0, 10), "Intercept")))extended_summary(read.file = "./data/processed/updated_regressions/foxp2_by_treatment_model.rds",
gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
remove.intercepts = TRUE)foxp2_by_treatment_model
| formula | family | priors | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | FoxP2.Counts.MMSt.Striatum ~ Treatment | gamma (log) | b-normal(0, 5) Intercept-normal(0, 10) shape-gamma(0.01, 0.01) | 1e+05 | 4 | 1 | 50000 | 0 (0%) | 0 | 145217.7 | 129358 | 411186403 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| TreatmentMediumstress | -0.013 | -0.284 | 0.257 | 1 | 150039.8 | 129358.0 |
| TreatmentHighstress | 0.071 | -0.198 | 0.342 | 1 | 145217.7 | 132119.9 |
FoxP2 mRNA expression MMST:VSP
qpcr <- read.csv("./data/raw/INBREmale_qPCR_10Jun25.csv")Stats
qpcr <- brm(formula = FoxP2.mRNA.MMST.VSP ~ Treatment, iter = 1e+05,
family = Gamma(link = "log"), silent = 2, data = qpcr, control = list(adapt_delta = 0.9,
max_treedepth = 15), chains = 4, cores = 4, file = "./data/processed/updated_regressions/qpcr_by_treatment_model",
file_refit = "on_change", prior = c(prior(normal(0, 5), "b"),
prior(normal(0, 10), "Intercept")))extended_summary(read.file = "./data/processed/updated_regressions/qpcr_by_treatment_model.rds",
gsub.pattern = "b_treatment|b_", gsub.replacement = "", highlight = TRUE,
remove.intercepts = TRUE)qpcr_by_treatment_model
| formula | family | priors | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | FoxP2.mRNA.MMST.VSP ~ Treatment | gamma (log) | b-normal(0, 5) Intercept-normal(0, 10) shape-gamma(0.01, 0.01) | 1e+05 | 4 | 1 | 50000 | 0 (0%) | 0 | 107766.9 | 84215.43 | 519692150 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| TreatmentHighStress | 0.437 | -0.293 | 1.129 | 1 | 107766.9 | 84215.43 |
Barplot and effect sizes graph
gg_foxp2 <- ggplot(data = foxp2, aes(x = Treatment, y = FoxP2.Counts.MMSt.Striatum,
fill = Treatment)) + geom_boxplot() + geom_point() + scale_x_discrete(labels = c(Control = "Control",
`Medium stress` = "Medium\nstress", `High stress` = "High\nstress")) +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + labs(y = "FoxP2 positive\ncells in MMST:VSP",
x = "Treatment") + theme_classic(base_size = bs) + theme(legend.position = "none")
gg_qpcr <- ggplot(data = qpcr, aes(x = Treatment, y = FoxP2.mRNA.MMST.VSP,
fill = Treatment)) + geom_boxplot() + geom_point() + scale_x_discrete(labels = c(Control = "Control",
`High Stress` = "High\nstress")) + scale_fill_viridis_d(begin = 0.1,
end = 0.9) + labs(y = "FoxP2 mRNA\nexpression MMST:VSP", x = "Treatment") +
theme_classic(base_size = bs) + theme(legend.position = "none")
foxp2_cells_mod <- readRDS("./data/processed/updated_regressions/foxp2_by_treatment_model.rds")
qpcr_mod <- readRDS("./data/processed/updated_regressions/qpcr_by_treatment_model.rds")
gg_coeffs_foxp2 <- lapply(list(foxp2_cells_mod, qpcr_mod), function(x) {
vars <- grep("b_", posterior::variables(x), value = TRUE)
draws <- posterior::as_draws_array(x, variable = vars)
coef_table <- draw_summary(draws, variables = vars, probs = c(0.025,
0.975), robust = TRUE, spread.type = "HPDI")
coef_table$predictor <- rownames(coef_table)
coef_table$predictor <- gsub("b_treatment|b_", "", coef_table$predictor)
coef_table$predictor <- gsub("Stress", " stress", coef_table$predictor)
coef_table$predictor <- gsub("week", "Week", coef_table$predictor)
coef_table$predictor <- gsub("Treatment", "", coef_table$predictor)
coef_table$predictor <- gsub("Mediumstress", "Medium stress",
coef_table$predictor)
coef_table$predictor <- gsub("Highstress", "High stress", coef_table$predictor)
coef_table <- coef_table[coef_table$predictor != "Intercept",
]
gg_coef <- ggplot2::ggplot(data = coef_table, aes(x = Estimate,
y = predictor)) + geom_vline(xintercept = 0, lty = 2) + ggplot2::geom_point(size = 4,
col = col_pointrange) + ggplot2::geom_errorbar(ggplot2::aes(xmin = `l-95% CI`,
xmax = `u-95% CI`), width = 0, col = col_pointrange) + ggplot2::theme_classic(base_size = bs) +
ggplot2::theme(axis.ticks.length = ggplot2::unit(0, "pt"),
plot.margin = ggplot2::margin(0, 0, 0, 0, "pt"), legend.position = "none",
strip.background = ggplot2::element_blank(), strip.text = ggplot2::element_blank()) +
ggplot2::labs(x = "Effect size", y = "")
return(gg_coef)
})
cowplot::plot_grid(gg_foxp2, gg_qpcr, gg_coeffs_foxp2[[1]] + theme_classic(base_size = bs),
gg_coeffs_foxp2[[2]] + theme_classic(base_size = bs), nrow = 2,
rel_heights = c(1.8, 1))
# try bs = 20 for saving plots
cowplot::ggsave2(filename = "./output/bar_graphs_and_estimates_foxp2_70dpi.jpeg",
width = 15, height = 10)
cowplot::ggsave2(filename = "./output/bar_graphs_and_estimates_foxp2_300dpi.jpeg",
dpi = 300, width = 15, height = 10)Combined model diagnostics
check_rds_fits(path = "./data/processed/updated_regressions", html = TRUE)| fit | formula | family | priors | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | behav-acoustic.convergence by treatment + week + (1|ID) + (1|round).rds | acoustic.convergence ~ treatment + week + (1 | ID) + (1 | round) | beta (logit) | b-normal(0, 5) Intercept-normal(0, 10) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 2.5) | 1e+05 | 4 | 1 | 50000 | 2682 (0.01341%) | 0 | 1508 | 906 | 1543370143 |
| 2 | behav-acoustic.diversity by treatment + week + (1|ID) + (1|round).rds | acoustic.diversity ~ treatment + week + (1 | ID) + (1 | round) | beta (logit) | b-normal(0, 5) Intercept-normal(0, 10) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 2.5) | 1e+05 | 4 | 1 | 50000 | 6528 (0.03264%) | 0 | 1419 | 644 | 195705932 |
| 3 | behav-acustic.plasticity by treatment + week + (1|ID) + (1|round).rds | acustic.plasticity ~ treatment + week + (1 | ID) + (1 | round) | beta (logit) | b-normal(0, 5) Intercept-normal(0, 10) phi-gamma(0.01, 0.01) sd-student_t(3, 0, 2.5) | 1e+05 | 4 | 1 | 50000 | 4613 (0.023065%) | 0 | 973 | 311 | 443607903 |
| 4 | behav-call.count by treatment + week + (1|ID) + (1|round).rds | call.count ~ treatment + week + (1 | ID) + (1 | round) | negbinomial (log) | b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | 1e+05 | 4 | 1 | 50000 | 3227 (0.016135%) | 0 | 1720 | 764 | 498436169 |
| 5 | foxp2_by_treatment_model.rds | FoxP2.Counts.MMSt.Striatum ~ Treatment | gamma (log) | b-normal(0, 5) Intercept-normal(0, 10) shape-gamma(0.01, 0.01) | 1e+05 | 4 | 1 | 50000 | 0 (0%) | 0 | 77833 | 99167 | 411186403 |
| 6 | physio-acustic.plasticity by treatment + week + (1|ID) + (1|round).rds | breath.rate ~ treatment + week + (1 | ID) + (1 | round) | student (identity) | b-normal(0, 5) Intercept-normal(0, 10) nu-gamma(2, 0.1) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) | 1e+05 | 4 | 1 | 50000 | 0 (0%) | 0 | 67534 | 63482 | 1702562168 |
| 7 | physio-baseline.CORT by treatment + week + (1|ID) + (1|round).rds | baseline.CORT ~ treatment + week + (1 | ID) + (1 | round) | student (identity) | b-normal(0, 5) Intercept-normal(0, 10) nu-gamma(2, 0.1) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) | 1e+05 | 4 | 1 | 50000 | 6255 (0.031275%) | 0 | 1288 | 419 | 662141751 |
| 8 | physio-breath.rate by treatment + week + (1|ID) + (1|round).rds | breath.rate ~ treatment + week + (1 | ID) + (1 | round) | student (identity) | b-normal(0, 5) Intercept-normal(0, 10) nu-gamma(2, 0.1) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) | 1e+05 | 4 | 1 | 50000 | 5808 (0.02904%) | 0 | 555 | 156 | 1411987597 |
| 9 | physio-stress.CORT by treatment + week + (1|ID) + (1|round).rds | stress.CORT ~ treatment + week + (1 | ID) + (1 | round) | student (identity) | b-normal(0, 5) Intercept-normal(0, 10) nu-gamma(2, 0.1) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) | 1e+05 | 4 | 1 | 50000 | 1796 (0.00898%) | 0 | 6443 | 2528 | 1854922145 |
| 10 | physio-stress.response by treatment + week + (1|ID) + (1|round).rds | stress.response ~ treatment + week + (1 | ID) + (1 | round) | student (identity) | b-normal(0, 5) Intercept-normal(0, 10) nu-gamma(2, 0.1) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) | 1e+05 | 4 | 1 | 50000 | 4766 (0.02383%) | 0 | 733 | 109 | 1337746393 |
| 11 | physio-weight by treatment + week + (1|ID) + (1|round).rds | weight ~ treatment + week + (1 | ID) + (1 | round) | student (identity) | b-normal(0, 5) Intercept-normal(0, 10) nu-gamma(2, 0.1) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10) | 1e+05 | 4 | 1 | 50000 | 4876 (0.02438%) | 0 | 708 | 223 | 1479872948 |
| 12 | qpcr_by_treatment_model.rds | FoxP2.mRNA.MMST.VSP ~ Treatment | gamma (log) | b-normal(0, 5) Intercept-normal(0, 10) shape-gamma(0.01, 0.01) | 1e+05 | 4 | 1 | 50000 | 0 (0%) | 0 | 58270 | 63807 | 519692150 |
Session information
## R version 4.5.0 (2025-04-11)
## 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 LAPACK version 3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.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] DT_0.33 corrplot_0.95 PhenotypeSpace_0.1.0
## [4] ggridges_0.5.6 posterior_1.6.1 ggrepel_0.9.6
## [7] cowplot_1.1.3 tidybayes_3.0.7 modelr_0.1.11
## [10] tidyr_1.3.1 forcats_1.0.0 purrr_1.0.4
## [13] dplyr_1.1.4 lme4_1.1-37 Matrix_1.7-3
## [16] brms_2.22.0 Rcpp_1.0.14 GGally_2.2.1
## [19] sp_2.2-0 MASS_7.3-65 formatR_1.14
## [22] knitr_1.50 kableExtra_1.4.0 ggplot2_3.5.2
## [25] viridis_0.6.5 viridisLite_0.4.2 pbapply_1.7-2
## [28] readxl_1.4.5 lubridate_1.9.4 remotes_2.5.0
##
## loaded via a namespace (and not attached):
## [1] RColorBrewer_1.1-3 tensorA_0.36.2.1 rstudioapi_0.17.1
## [4] jsonlite_2.0.0 magrittr_2.0.3 TH.data_1.1-3
## [7] estimability_1.5.1 spatstat.utils_3.1-4 farver_2.1.2
## [10] nloptr_2.2.1 rmarkdown_2.29 ragg_1.4.0
## [13] vctrs_0.6.5 minqa_1.2.4 spatstat.explore_3.4-2
## [16] terra_1.8-42 htmltools_0.5.8.1 curl_6.4.0
## [19] distributional_0.5.0 broom_1.0.8 raster_3.6-32
## [22] cellranger_1.1.0 StanHeaders_2.32.10 sass_0.4.10
## [25] KernSmooth_2.23-26 bslib_0.9.0 htmlwidgets_1.6.4
## [28] plyr_1.8.9 sandwich_3.1-1 emmeans_1.11.1
## [31] zoo_1.8-14 cachem_1.1.0 lifecycle_1.0.4
## [34] pkgconfig_2.0.3 R6_2.6.1 fastmap_1.2.0
## [37] rbibutils_2.3 digest_0.6.37 rprojroot_2.0.4
## [40] tensor_1.5 textshaping_1.0.1 crosstalk_1.2.1
## [43] vegan_2.6-10 labeling_0.4.3 spatstat.sparse_3.1-0
## [46] timechange_0.3.0 mgcv_1.9-3 polyclip_1.10-7
## [49] abind_1.4-8 compiler_4.5.0 proxy_0.4-27
## [52] withr_3.0.2 inline_0.3.21 backports_1.5.0
## [55] DBI_1.2.3 ggstats_0.9.0 QuickJSR_1.7.0
## [58] pkgbuild_1.4.8 classInt_0.4-11 permute_0.9-7
## [61] loo_2.8.0 tools_4.5.0 units_0.8-7
## [64] goftest_1.2-3 glue_1.8.0 nlme_3.1-168
## [67] grid_4.5.0 sf_1.0-20 checkmate_2.3.2
## [70] reshape2_1.4.4 cluster_2.1.8.1 generics_0.1.4
## [73] gtable_0.3.6 spatstat.data_3.1-6 class_7.3-23
## [76] xml2_1.3.8 spatstat.geom_3.3-6 pillar_1.11.0
## [79] ggdist_3.3.3 stringr_1.5.1 splines_4.5.0
## [82] lattice_0.22-7 survival_3.8-3 deldir_2.0-4
## [85] tidyselect_1.2.1 reformulas_0.4.1 arrayhelpers_1.1-0
## [88] gridExtra_2.3 V8_6.0.4 svglite_2.1.3
## [91] stats4_4.5.0 xfun_0.52 bridgesampling_1.1-2
## [94] matrixStats_1.5.0 rstan_2.32.7 stringi_1.8.7
## [97] yaml_2.3.10 boot_1.3-31 evaluate_1.0.3
## [100] codetools_0.2-20 tibble_3.3.0 cli_3.6.5
## [103] RcppParallel_5.1.10 xtable_1.8-4 systemfonts_1.2.3
## [106] Rdpack_2.6.4 jquerylib_0.1.4 spatstat.random_3.4-1
## [109] coda_0.19-4.1 svUnit_1.0.6 spatstat.univar_3.1-2
## [112] parallel_4.5.0 rstantools_2.4.0 bayesplot_1.12.0
## [115] Brobdingnag_1.2-9 mvtnorm_1.3-3 scales_1.4.0
## [118] e1071_1.7-16 rlang_1.1.6 multcomp_1.4-28