Variation in vocal behavior under varying stress levels V01

Budgie INBRE stress and contact calls

Load packages

## add 'developer' to packages to be installed from github
x <- c("remotes", "lubridate", "readxl", "pbapply", "viridis", "ggplot2", "kableExtra", "knitr", "formatR", "randomForest", "Rtsne", "MASS", "adehabitatHR", "sp", "GGally", "spatstat", "raster", "vegan", "brms", "lme4", "dplyr", "purrr", "forcats", "tidyr", "modelr", "tidybayes", "cowplot", "ggrepel", "posterior", "ggridges", "maRce10/PhenotypeSpace", "maRce10/brmsish")

aa <- lapply(x, function(y) {
  
  # get pakage name
  pkg <- strsplit(y, "/")[[1]]
  pkg <- pkg[length(pkg)]
  
  # check if installed, if not then install 
  if (!pkg %in% installed.packages()[,"Package"])  {

      if (grepl("/", y))  remotes::install_github(y, force = TRUE) else
    install.packages(y) 
    }

  # load package
  a <- try(require(pkg, character.only = T), silent = T)

  if (!a) remove.packages(pkg)
  })

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")


cols <- viridis(10, alpha = 0.7)

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

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

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

  }

# 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")

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"))

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, ]

Exploratory graph

Physiological parameters

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 <- factor(stress$Treatment, levels = c("Control", "Medium 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)

ggplot(data = stress, aes(x = day, y = weight, group = Bird.ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Weight (g)", x = "Day sample was taken") +
    theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
    labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))

ggplot(data = stress, aes(x = day, y = breath.count, group = Bird.ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Breath count", x = "Day sample was taken") +
    theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
    labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))

ggplot(data = stress, aes(x = day, y = cort.baseline, group = Bird.ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Baseline CORT", x = "Day sample was taken") +
    theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
    labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))

ggplot(data = stress, aes(x = day, y = cort.stress, group = Bird.ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Stress CORT", x = "Day sample was taken") +
    theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
    labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))

ggplot(data = stress, aes(x = day, y = stress.response, group = Bird.ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Stress response", x = "Day sample was taken") +
    theme_classic(base_size = 30) + scale_x_continuous(breaks = c(3, 7, 14, 21, 28),
    labels = c(3, 7, 14, 21, 28)) + theme(legend.position = c(0.8, 0.25))

# stress_l <- lapply(unique(stress$Bird.ID), function(x){ X <-
# stress[stress$Bird.ID == x, ] X$cort.baseline.change <- X$cort.baseline - X$
# })

standard_error <- function(x) sd(x)/sqrt(length(x))

agg_stress <- aggregate(cort.baseline ~ week + treatment, stress, mean)
agg_stress$se <- aggregate(cort.baseline ~ week + treatment, stress, standard_error)[,
    3]
agg_stress$Week <- 1:4

ggplot(data = agg_stress, aes(x = Week, y = cort.baseline, fill = treatment)) + geom_bar(stat = "identity") +
    geom_errorbar(aes(ymin = cort.baseline - se, ymax = cort.baseline + se), width = 0.1) +
    # geom_line(size = 1.2) +
scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~treatment, ncol = 3, scale = "fixed") +
    labs(y = "Baseline CORT", x = "Week") + theme_classic(base_size = 30) + theme(legend.position = "none")

FOXP2

foxp2 <- read.csv("./data/raw/stress_IHC_counts.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"))

ggplot(data = foxp2, aes(x = Treatment, y = FoxP2.Counts.MMSt.Striatum, fill = Treatment)) +
    geom_boxplot() + geom_point() + scale_fill_viridis_d(begin = 0.1, end = 0.9) +
    labs(y = "% FoxP2 expression", x = "Treatment") + theme_classic(base_size = 30) +
    # scale_x_discrete( labels = c(1:4)) +
theme(legend.position = "none")

Behavioral parameters

Call counts per treatment

call_count <- aggregate(selec ~ week + treatment + round, data = sels, length)
call_count$Week <- factor(call_count$week, levels = (unique(call_count$week))[5:1])

call_count$round <- paste("Round", call_count$round)

call_count$treatment <- factor(call_count$treatment, levels = c("Control", "Medium Stress",
    "High Stress"))

call_count$stress.response <- aggregate(selec ~ week + treatment + round, data = sels,
    mean)[, 4]

ggplot(data = call_count, aes(x = treatment, y = selec, fill = Week)) + geom_bar(stat = "identity") +
    scale_fill_viridis_d(begin = 0.1, end = 0.9) + coord_flip() + facet_wrap(~round,
    ncol = 2) + labs(y = "Number of calls", x = "Treatment") + theme_classic(base_size = 30) +
    theme(legend.position = c(0.8, 0.15))

call_count2 <- aggregate(selec ~ week + treatment + round + ID, data = sels, length)
agg_count <- aggregate(selec ~ week + treatment, call_count2, mean)
agg_count$se <- aggregate(selec ~ week + treatment, call_count2, standard_error)[,
    3]

agg_count$treatment <- factor(agg_count$treatment, levels = c("Control", "Medium Stress",
    "High Stress"))

ggplot(data = agg_count, aes(x = week, y = selec, fill = treatment)) + geom_bar(stat = "identity") +
    geom_errorbar(aes(ymin = selec - se, ymax = selec + se), width = 0.1) + scale_fill_viridis_d(begin = 0.1,
    end = 0.9) + facet_wrap(~treatment, ncol = 3, scale = "fixed") + labs(y = "Number of calls",
    x = "Week") + theme_classic(base_size = 30) + theme(legend.position = "none")

Call counts per bird and treatment

# Call counts per treatment by individual (in color)

call_count <- aggregate(selec ~ week + treatment + round + ID, data = sels, length)
call_count$Week <- factor(call_count$week, levels = (1:5))

call_count$round <- paste("Round", call_count$round)

call_count$treatment <- factor(call_count$treatment, levels = c("Control", "Medium Stress",
    "High Stress"))

ggplot(data = call_count, aes(x = week, y = selec, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Call count", x = "Week") + theme_classic(base_size = 30) +
    theme(legend.position = c(0.8, 0.25))

Acoustic parameters

call_week_params <- aggregate(cbind(duration, meanfreq, sd, freq.median, freq.IQR,
    time.IQR, skew, kurt, sp.ent, time.ent, entropy, meandom, mindom, maxdom, dfrange,
    modindx, meanpeakf) ~ week + treatment + round, data = sels, mean)

call_week_params_sd <- aggregate(cbind(duration, meanfreq, sd, freq.median, freq.IQR,
    time.IQR, skew, kurt, sp.ent, time.ent, entropy, meandom, mindom, maxdom, dfrange,
    modindx, meanpeakf) ~ week + treatment + round, data = sels, sd)

names(call_week_params_sd) <- names(call_week_params) <- c("week", "treatment", "round",
    "Duration", "Mean_frequency_(kHz)", "SD_frequency_(kHz)", "Median_frequency_(kHz)",
    "Interquantile_frequency_range_(kHz)", "Interquantile_time_range_(s)", "Skewness",
    "Kurtosis", "Spectral_entropy", "Time_entropy", "Overall entropy", "Mean_dominant_frequency_(kHz)",
    "Minimum_dominant_frequency_(kHz)", "Maximum_dominant_frequency_(kHz)", "Dominant_frequency_range_(kHz)",
    "Modulation_index", "Mean_peak_frequency_(kHz)")


call_week_params_sd <- call_week_params_sd[, c("Duration", "Mean_frequency_(kHz)",
    "SD_frequency_(kHz)", "Median_frequency_(kHz)", "Interquantile_frequency_range_(kHz)",
    "Interquantile_time_range_(s)", "Skewness", "Kurtosis", "Spectral_entropy", "Time_entropy",
    "Overall entropy", "Mean_dominant_frequency_(kHz)", "Minimum_dominant_frequency_(kHz)",
    "Maximum_dominant_frequency_(kHz)", "Dominant_frequency_range_(kHz)", "Modulation_index",
    "Mean_peak_frequency_(kHz)")]

names(call_week_params_sd) <- paste(names(call_week_params_sd), "sd", sep = ".")


call_week_params <- cbind(call_week_params, call_week_params_sd)
call_week_params$round <- paste("Round", call_week_params$round)
call_week_params$treatment <- factor(call_week_params$treatment, levels = c("Control",
    "Medium Stress", "High Stress"))


pd <- position_dodge(0.3)

ggs <- lapply(c("Duration", "Mean_frequency_(kHz)", "SD_frequency_(kHz)", "Median_frequency_(kHz)",
    "Interquantile_frequency_range_(kHz)", "Interquantile_time_range_(s)", "Skewness",
    "Kurtosis", "Spectral_entropy", "Time_entropy", "Overall entropy", "Mean_dominant_frequency_(kHz)",
    "Minimum_dominant_frequency_(kHz)", "Maximum_dominant_frequency_(kHz)", "Dominant_frequency_range_(kHz)",
    "Modulation_index", "Mean_peak_frequency_(kHz)"), function(i) {
    call_week_params$var <- call_week_params[, i]
    call_week_params$var.min <- call_week_params$var - call_week_params[, paste(i,
        "sd", sep = ".")]
    call_week_params$var.max <- call_week_params$var + call_week_params[, paste(i,
        "sd", sep = ".")]

    ggplot(call_week_params, aes(x = week, y = var, col = treatment)) + geom_point(size = 4,
        position = pd) + geom_errorbar(aes(ymin = var.min, ymax = var.max), width = 0.3,
        position = pd, size = 1.7) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
        facet_wrap(~round, ncol = 2) + labs(y = i, x = "Week") + theme_classic(base_size = 30) +
        theme(legend.position = c(0.8, 0.25))
})

ggs

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

## 
## [[12]]

## 
## [[13]]

## 
## [[14]]

## 
## [[15]]

## 
## [[16]]

## 
## [[17]]

Acoustic space projection

scale_acous_param <- scale(sels[, c("duration", "meanfreq", "sd", "freq.median",
    "freq.IQR", "time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom",
    "mindom", "maxdom", "dfrange", "modindx", "meanpeakf")])

urfmod <- randomForest(x = scale_acous_param, importance = TRUE, proximity = TRUE,
    ntree = 10000)

saveRDS(urfmod, "./data/processed/unsupervised_random_forest_model_budgie_calls_stress_jun_2021.RDS")

mds_urf_prox <- cmdscale(1 - urfmod$proximity, k = 2)

saveRDS(mds_urf_prox, "./data/processed/MDS_unsupervised_random_forest_model_budgie_calls_stress_jun_2021.RDS")

Random forest

urfmod <- readRDS("./data/processed/unsupervised_random_forest_model_budgie_calls_stress_jun_2021.RDS")

mds_urf_prox <- readRDS("./data/processed/MDS_unsupervised_random_forest_model_budgie_calls_stress_jun_2021.RDS")

sels$MDS1 <- mds_urf_prox[, 1]
sels$MDS2 <- mds_urf_prox[, 2]

print(1 - sum(urfmod$votes[, 1] > urfmod$votes[, 2])/nrow(urfmod$votes))

## [1] 0.005058366

Proportion of real data classified as synthetic data by the unsupervised random forest: 0.005058

ggplot(sels, aes(x = MDS1, y = MDS2, col = as.factor(round))) + geom_point() + labs(color = "Round") +
    scale_color_viridis_d(alpha = 0.4) + theme_classic(base_size = 30) + theme(legend.position = c(0.62,
    0.3)) + guides(colour = guide_legend(override.aes = list(size = 10)))

PCA

pca <- prcomp(x = sels[, c("duration", "meanfreq", "sd", "freq.median", "freq.IQR",
    "time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom", "mindom",
    "maxdom", "dfrange", "modindx", "meanpeakf")], scale. = TRUE)

Y <- as.data.frame(pca$x[, 1:2])
names(Y) <- c("PC1", "PC2")

sels <- data.frame(sels, Y)

ggplot(sels, aes(x = PC1, y = PC2, col = as.factor(round))) + geom_point() + labs(color = "Round") +
    scale_color_viridis_d(alpha = 0.4) + theme_classic(base_size = 25) + theme(legend.position = c(0.9,
    0.8)) + guides(colour = guide_legend(override.aes = list(size = 10)))

t-SNE

scale_param <- scale(sels[, c("duration", "meanfreq", "sd", "freq.median", "freq.IQR",
    "time.IQR", "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom", "mindom",
    "maxdom", "dfrange", "modindx", "meanpeakf")])

tsne <- Rtsne(scale_param, dims = 2, perplexity = 30, verbose = FALSE, max_iter = 5000)

saveRDS(tsne, "./data/processed/tsne_on_acoustic_parameters_jun_2021.RDS")

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)

ggplot(sels, aes(x = TSNE1, y = TSNE2, col = as.factor(round))) + geom_point() +
    labs(color = "Round") + scale_color_viridis_d(alpha = 0.4) + theme_classic(base_size = 25) +
    theme(legend.position = c(0.955, 0.25)) + guides(colour = guide_legend(override.aes = list(size = 10)))

sels$treatment <- factor(sels$treatment, levels = c("Control", "Medium Stress", "High Stress"))

ggplot(sels, aes(x = TSNE1, y = TSNE2, col = as.factor(treatment))) + geom_point() +
    labs(color = "Treatment") + scale_color_viridis_d(alpha = 0.4) + theme_classic(base_size = 25) +
    theme(legend.position = c(0.88, 0.2)) + guides(colour = guide_legend(override.aes = list(size = 10)))

Acoustic space by individual (each point) for the 3 dimension reduction methods

• Only individuals with at least 20 calls
• Using rarefaction with n = the minimum sample size across all individuals

tab <- table(sels$ID)

Y <- sels[sels$ID %in% names(tab[tab > 20]), ]

pca_areas <- rarefact_space_size(X = Y, dimensions = c("PC1", "PC2"), group = "ID",
    parallel = 4, pb = TRUE, type = "density")

rf_mds_areas <- rarefact_space_size(X = Y, dimensions = c("MDS1", "MDS2"), group = "ID",
    parallel = 4, pb = TRUE, type = "density")

tsne_areas <- rarefact_space_size(X = Y, dimensions = c("TSNE1", "TSNE2"), group = "ID",
    parallel = 4, pb = TRUE, type = "density")

areas <- data.frame(pca_areas[, -ncol(pca_areas)], pca.area = pca_areas$mean.size,
    rf.mds.area = rf_mds_areas$mean.size, tsne.area = tsne_areas$mean.size)

saveRDS(areas, "./data/processed/acoustic_space_area_by_individual_for_each_dimension_reduction_method.RDS")

areas <- readRDS("./data/processed/acoustic_space_area_by_individual_for_each_dimension_reduction_method.RDS")

ggpairs(areas, columns = which(names(areas) %in% c("pca.area", "rf.mds.area", "tsne.area")),
    columnLabels = c("PCA", "Rand.for. MDS", "t-SNE")) + theme_minimal(base_size = 30)

Individual acoustic space across weeks

• PCA and t-SNE used for dimensionality reduction
• Only weeks with at least 6 calls were included

PCA

sels$ID.week <- paste(sels$ID, sels$week, sep = "-")

tab <- table(sels$ID.week)

# sort(tab)

Y <- sels[sels$ID.week %in% names(tab)[tab >= 6], ]

# run with rarefaction
areas_by_week_raref <- rarefact_space_size(X = Y, dimensions = c("PC1", "PC2"), group = "ID.week",
    parallel = 1, pb = TRUE, type = "density")

# run without rarefaction
areas_by_week <- space_size(X = Y, dimensions = c("PC1", "PC2"), group = "ID.week",
    parallel = 1, pb = TRUE, type = "density", extent = 2)

areas_by_week$raref.area <- areas_by_week_raref$mean.size

# add metadata
areas_by_week$ID <- sapply(strsplit(areas_by_week$group, "-"), "[[", 1)

areas_by_week$week <- factor(sapply(strsplit(areas_by_week$group, "-"), "[[", 2),
    levels = 1:5)

areas_by_week$treatment <- sapply(1:nrow(areas_by_week), function(x) metadat$Treatment[metadat$Bird.ID ==
    areas_by_week$ID[x]][1])

areas_by_week$round <- sapply(1:nrow(areas_by_week), function(x) metadat$Round[metadat$Bird.ID ==
    areas_by_week$ID[x]][1])

areas_by_week$round <- paste("Round", areas_by_week$round)

call_count$treatment <- factor(call_count$treatment, levels = c("Control", "Medium Stress",
    "High Stress"))

areas_by_week$treatment <- factor(areas_by_week$treatment, levels = c("Control",
    "Medium Stress", "High Stress"))

saveRDS(areas_by_week, "./data/processed/acoustic_space_area_by_individual_and_week.RDS")

With rarefaction

areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week.RDS")

ggplot(data = areas_by_week, aes(x = week, y = raref.area, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
    theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

With rarefaction and compared to week 1

areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week.RDS")

areas_by_week <- do.call(rbind, lapply(unique(areas_by_week$ID), function(x) {

    X <- areas_by_week[areas_by_week$ID == x, ]
    X$raref.area <- X$raref.area - X$raref.area[X$week == min(as.numeric(as.character(X$week)))]
    return(X)
}))

ggplot(data = areas_by_week, aes(x = week, y = raref.area, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
    theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

Without rarefaction

ggplot(data = areas_by_week[areas_by_week$week != 5, ], aes(x = week, y = size, group = ID,
    col = treatment)) + geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1,
    end = 0.9) + facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
    theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

Without rarefaction and compared to week 1

areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week.RDS")

areas_by_week <- do.call(rbind, lapply(unique(areas_by_week$ID), function(x) {

    X <- areas_by_week[areas_by_week$ID == x, ]
    X$size <- X$size - X$size[X$week == min(as.numeric(as.character(X$week)))]
    return(X)
}))

ggplot(data = areas_by_week, aes(x = week, y = size, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
    theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

Acoustic space area by number of calls per individual and week

ggplot(data = areas_by_week, aes(x = n, y = size, col = treatment)) + geom_point() +
    # geom_smooth(size = 1.2, method = 'lm') +
scale_color_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~round, ncol = 2, scales = "free_x") +
    labs(y = "Acoustic space area", x = "Number of calls (n)") + theme_classic(base_size = 30) +
    theme(legend.position = c(0.8, 0.25))

• Acoustic space area by number of calls per individual and week by treatment
• Must have at least 2 calls

ggplot(data = areas_by_week, aes(x = n, y = log(size + 10), col = treatment)) + geom_point() +
    geom_smooth(size = 1.2, method = "lm") + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    # facet_wrap(~ round, ncol = 2, scales = 'free_x') +
labs(y = "log(acoustic space area)", x = "Number of calls (n)") + theme_classic(base_size = 30) +
    theme(legend.position = c(0.8, 0.7))

Save acoustic space plots

plot_space(X = sels, indices = which(sels$ID == "395WBHM"), dimensions = c("TSNE1",
    "TSNE2"), point.alpha = 0.2, pch = 1)
#
plot_space(X = sels, indices = which(sels$ID == "395WBHM"), dimensions = c("TSNE1",
    "TSNE2"), background.indices = which(sels$ID != "395WBHM" & sels$record.group ==
    sels$record.group[sels$ID == "395WBHM"][1]), point.alpha = 0.2, pch = 1)

# Y <- sels[sels$ID == '395WBHM',] Y <- sels[sels$ID %in% c('395WBHM',
# '394WBHM', '360YGHM', '385YBHM'), ]

# trait_overlap(Y, dimensions = c('TSNE1', 'TSNE2'), group = 'ID', type =
# 'density', overlap.type = 'assymetric')


# X = Y dimensions = c('TSNE1', 'TSNE2') level <- '395WBHM' basecex <- 1 group
# <- 'ID' indices <- which(X[, group] == level)

# with points
out <- pblapply(unique(sels$ID), function(x) {

    tiff(filename = file.path("./images", paste(x, sels$treatment[sels$ID == x][1],
        "v2.tiff", sep = "-")), res = 100, width = 1200, height = 1200)

    par(mfrow = c(2, 2))

    # W <- sels[sels$ID == x, ]

    for (i in 1:4) try(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]), title = paste(x,
        "week", i), basecex = 1, labels = c(x, "group")), silent = TRUE)

    dev.off()
})

TSNE

sels$ID.week <- paste(sels$ID, sels$week, sep = "-")

tab <- table(sels$ID.week)

# sort(tab)

Y <- sels[sels$ID.week %in% names(tab)[tab >= 6], ]

# run with rarefaction
areas_by_week_raref <- rarefact_space_size(X = Y, dimensions = c("TSNE1", "TSNE2"),
    group = "ID.week", parallel = 1, pb = TRUE, type = "density")

# run without rarefaction
areas_by_week <- space_size(X = Y, dimensions = c("TSNE1", "TSNE2"), group = "ID.week",
    parallel = 4, pb = TRUE, type = "density")

areas_by_week$raref.area <- areas_by_week_raref$mean.size

# add metadata
areas_by_week$ID <- sapply(strsplit(areas_by_week$group, "-"), "[[", 1)

areas_by_week$week <- factor(sapply(strsplit(areas_by_week$group, "-"), "[[", 2),
    levels = 1:5)

areas_by_week$treatment <- sapply(1:nrow(areas_by_week), function(x) metadat$Treatment[metadat$Bird.ID ==
    areas_by_week$ID[x]][1])

areas_by_week$round <- sapply(1:nrow(areas_by_week), function(x) metadat$Round[metadat$Bird.ID ==
    areas_by_week$ID[x]][1])

areas_by_week$round <- paste("Round", areas_by_week$round)

call_count$treatment <- factor(call_count$treatment, levels = c("Control", "Medium Stress",
    "High Stress"))

areas_by_week$treatment <- factor(areas_by_week$treatment, levels = c("Control",
    "Medium Stress", "High Stress"))

saveRDS(areas_by_week, "./data/processed/acoustic_space_area_by_individual_and_week_TSNE.RDS")

With rarefaction

areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week_TSNE.RDS")

ggplot(data = areas_by_week, aes(x = week, y = raref.area, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
    theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

With rarefaction and compared to week 1

areas_by_week <- readRDS("./data/processed/acoustic_space_area_by_individual_and_week_TSNE.RDS")

areas_by_week <- do.call(rbind, lapply(unique(areas_by_week$ID), function(x) {

    X <- areas_by_week[areas_by_week$ID == x, ]
    X$raref.area <- X$raref.area - X$raref.area[X$week == min(as.numeric(as.character(X$week)))]
    return(X)
}))

ggplot(data = areas_by_week, aes(x = week, y = raref.area, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
    theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

Without rarefaction

ggplot(data = areas_by_week[areas_by_week$week != 5, ], aes(x = week, y = size, group = ID,
    col = treatment)) + geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1,
    end = 0.9) + facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
    theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

Without rarefaction and compared to week 1

areas_by_week <- do.call(rbind, lapply(unique(areas_by_week$ID), function(x) {

    X <- areas_by_week[areas_by_week$ID == x, ]
    X$size <- X$size - X$size[X$week == min(as.numeric(as.character(X$week)))]
    return(X)
}))

ggplot(data = areas_by_week, aes(x = week, y = size, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Acoustic space area", x = "Week") +
    theme_classic(base_size = 30) + theme(legend.position = c(0.8, 0.25))

Acoustic space area by number of calls per individual and week

ggplot(data = areas_by_week, aes(x = n, y = size, col = treatment)) + geom_point() +
    # geom_smooth(size = 1.2, method = 'lm') +
scale_color_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~round, ncol = 2, scales = "free_x") +
    labs(y = "Acoustic space area", x = "Number of calls (n)") + theme_classic(base_size = 30) +
    theme(legend.position = c(0.8, 0.25))

• Acoustic space area by number of calls per individual and week by treatment
• Must have at least 2 calls

ggplot(data = areas_by_week, aes(x = n, y = log(size + 10), col = treatment)) + geom_point() +
    geom_smooth(size = 1.2, method = "lm") + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    # facet_wrap(~ round, ncol = 2, scales = 'free_x') +
labs(y = "log(acoustic space area)", x = "Number of calls (n)") + theme_classic(base_size = 30) +
    theme(legend.position = c(0.8, 0.7))

Acoustic space overlap

• 2 methods: pairwise distance between observations and kernel density overlap

Individual overlap to its own group

PCA

sels$ID.week <- paste(sels$ID, sels$week, sep = "-")
tab <- table(sels$ID)
Y <- sels[sels$ID %in% names(tab)[tab >= 2], ]

group_ovlp_l <- lapply(unique(Y$ID.week), function(x) {

    # group data
    X <- Y[Y$record.group == Y$record.group[Y$ID.week == x][1] & Y$week == Y$week[Y$ID.week ==
        x][1], ]

    # label all other individuals as group
    X$ID2 <- ifelse(X$ID.week == x, "focal", "group")

    # run with density
    ovlp <- if (min(table(X$ID2)) > 4)
        rarefact_space_similarity(X = X, dimensions = c("PC1", "PC2"), group = "ID2",
            parallel = 1, pb = FALSE, outliers = 0.95, type = "mean.density.overlap") else matrix(rep(NA, 3), nrow = 1)

    dsts <- if (min(table(X$ID2)) > 1)
        rarefact_space_similarity(X = X, dimensions = c("PC1", "PC2"), group = "ID2",
            parallel = 1, pb = FALSE, outliers = 0.95, type = "distance") else matrix(rep(NA, 3), nrow = 1)

    out <- data.frame(ID = Y$ID[Y$ID.week == x][1], group = Y$record.group[Y$ID.week ==
        x][1], week = Y$week[Y$ID.week == x][1], overlap.to.group = ovlp[1, 3], distance.to.group = dsts[1,
        3], treatment = Y$treatment[Y$ID.week == x][1], round = Y$round[Y$ID.week ==
        x][1])

    return(out)

})

group_ovlp <- do.call(rbind, group_ovlp_l)

saveRDS(group_ovlp, "./data/processed/acoustic_space_density_overlap_to_group_by_week.RDS")

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"))

ggplot(data = group_ovlp, aes(x = week, y = overlap.to.group, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Individual overlap to group acoustic space",
    x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
    0.25))

ggplot(data = group_ovlp, aes(x = week, y = overlap.to.group, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Individual distance to group acoustic space",
    x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
    0.25))

sels$ID.week <- paste(sels$ID, sels$week, sep = "-")

tab <- table(sels$ID.week)

Y <- sels[sels$ID.week %in% names(tab)[tab >= 6], ]

# run with rarefaction
overlap_by_week <- rarefact_space_similarity(X = Y, dimensions = c("PC1", "PC2"),
    group = "ID.week", parallel = 1, pb = TRUE, outliers = 0.95, type = "mean.density.overlap")

saveRDS(overlap_by_week, "./data/processed/acoustic_space_overlap_by_individual_and_week.RDS")

TSNE

sels$ID.week <- paste(sels$ID, sels$week, sep = "-")
tab <- table(sels$ID)
Y <- sels[sels$ID %in% names(tab)[tab >= 6], ]

group_ovlp_l <- lapply(unique(Y$ID.week), function(x) {

    # group data
    X <- Y[Y$record.group == Y$record.group[Y$ID.week == x][1] & Y$week == Y$week[Y$ID.week ==
        x][1], ]

    # label all other individuals as group
    X$ID2 <- ifelse(X$ID.week == x, "focal", "group")

    # run with density
    ovlp <- if (min(table(X$ID2)) >= 6)
        rarefact_space_similarity(X = X, dimensions = c("TSNE1", "TSNE2"), group = "ID2",
            parallel = 1, pb = FALSE, type = "mean.density.overlap")[, 1:3] else matrix(rep(NA, 3), nrow = 1)

    dsts <- if (min(table(X$ID2)) >= 2)
        rarefact_space_similarity(X = X, dimensions = c("TSNE1", "TSNE2"), group = "ID2",
            parallel = 1, pb = FALSE, type = "distance")[, 1:3] else matrix(rep(NA, 3), nrow = 1)

    out <- data.frame(ID = X$ID[X$ID.week == x][1], group = X$record.group[X$ID.week ==
        x][1], week = X$week[X$ID.week == x][1], overlap.to.group = ovlp[1, 3], distance.to.group = dsts[1,
        3], treatment = X$treatment[X$ID.week == x][1], round = X$round[X$ID.week ==
        x][1])

    return(out)

})

group_ovlp <- do.call(rbind, group_ovlp_l)


saveRDS(group_ovlp, "./data/processed/acoustic_space_density_overlap_to_group_by_week_TSNE.RDS")

group_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_group_by_week_TSNE.RDS")


# comparing to week 1 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"))

ggplot(data = group_ovlp, aes(x = week, y = overlap.to.group, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Individual overlap to group acoustic space",
    x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
    0.25))

Change in overlap with self over time

PCA

tab <- table(sels$ID)
Y <- sels[sels$ID %in% names(tab)[tab >= 6], ]

indiv_ovlp_l <- lapply(unique(Y$ID), function(x) {

    # group data
    X <- Y[Y$ID == x, ]
    tab_week <- table(X$week)
    X <- X[X$week %in% names(tab_week)[tab_week >= 6], ]

    if (any(X$week == 1)) {
        # run with density
        ovlp <- try(rarefact_space_similarity(X = X, dimensions = c("PC1", "PC2"),
            group = "week", parallel = 1, pb = FALSE, type = "mean.density.overlap"),
            silent = TRUE)

        dsts <- try(rarefact_space_similarity(X = X, dimensions = c("PC1", "PC2"),
            group = "week", parallel = 1, pb = FALSE, type = "distance"), silent = TRUE)

        out <- data.frame(ID = x, group = X$record.group[X$ID == x][1], week = if (is(ovlp,
            "try-error"))
            NA else ovlp$group.2[ovlp$group.1 == 1], overlap.to.first.week = if (is(ovlp,
            "try-error"))
            NA else ovlp$mean.overlap[ovlp$group.1 == 1], distance.to.first.week = if (is(dsts,
            "try-error"))
            NA else dsts$mean.distance[dsts$group.1 == 1], treatment = X$treatment[X$ID ==
            x][1], round = X$round[X$ID == x][1])
    } else out <- data.frame(ID = x, group = X$record.group[X$ID == x][1], week = NA,
        overlap.to.first.week = NA, distance.to.first.week = NA, treatment = X$treatment[X$ID ==
            x][1], round = X$round[X$ID == x][1])
    return(out)

})

indiv_ovlp <- do.call(rbind, indiv_ovlp_l)

saveRDS(indiv_ovlp, "./data/processed/acoustic_space_density_overlap_to_first_week_by_individual.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"))
indiv_ovlp <- indiv_ovlp[!is.na(indiv_ovlp$group), ]

ggplot(data = indiv_ovlp, aes(x = week, y = overlap.to.first.week, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Individual overlap to acoustic space from first week",
    x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
    0.25))

ggplot(data = indiv_ovlp, aes(x = week, y = distance.to.first.week, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Individual distance to acoustic space from first week",
    x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
    0.25))

agg_dist <- aggregate(distance.to.first.week ~ week + treatment, indiv_ovlp, mean)
agg_dist$se <- aggregate(distance.to.first.week ~ week + treatment, indiv_ovlp, standard_error)[,
    3]
agg_dist$Week <- 1:4

# ggplot(data = agg_dist, aes(x = Week, y = distance.to.first.week, fill =
# treatment)) + geom_bar(stat = 'identity') +
# geom_errorbar(aes(ymin=distance.to.first.week - se, ymax =
# distance.to.first.week + se), width=.1) + # geom_line(size = 1.2) +
# scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~ treatment, ncol =
# 3, scale = 'fixed') + labs(y = 'Distance to self first week', x = 'Week') +
# theme_classic(base_size = 30) + # scale_x_discrete( labels = c(1:4)) +
# theme(legend.position = 'none')


agg_ovlp <- aggregate(overlap.to.first.week ~ week + treatment, indiv_ovlp, mean)
agg_ovlp$se <- aggregate(overlap.to.first.week ~ week + treatment, indiv_ovlp, standard_error)[,
    3]
agg_ovlp$Week <- 1:4

# ggplot(data = agg_ovlp, aes(x = Week, y = overlap.to.first.week, fill =
# treatment)) + geom_bar(stat = 'identity') +
# geom_errorbar(aes(ymin=overlap.to.first.week - se, ymax =
# overlap.to.first.week + se), width=.1) + # geom_line(size = 1.2) +
# scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~ treatment, ncol =
# 3, scale = 'fixed') + labs(y = 'Overlap to \nself first week', x = 'Week') +
# theme_classic(base_size = 30) + # scale_x_discrete( labels = c(1:4)) +
# theme(legend.position = 'none')

TSNE

tab <- table(sels$ID)
Y <- sels[sels$ID %in% names(tab)[tab >= 2], ]

indiv_ovlp_l <- lapply(unique(Y$ID), function(x) {

    # group data
    X <- Y[Y$ID == x, ]
    tab_week <- table(X$week)
    X <- X[X$week %in% names(tab_week)[tab_week >= 6], ]

    if (any(X$week == 1)) {
        # run with density
        ovlp <- try(rarefact_space_similarity(X = X, dimensions = c("TSNE1", "TSNE2"),
            group = "week", parallel = 1, pb = FALSE, type = "mean.density.overlap"),
            silent = TRUE)

        dsts <- try(rarefact_space_similarity(X = X, dimensions = c("TSNE1", "TSNE2"),
            group = "week", parallel = 1, pb = FALSE, type = "distance"), silent = TRUE)

        out <- data.frame(ID = x, group = X$record.group[X$ID == x][1], week = if (is(ovlp,
            "try-error"))
            NA else ovlp$group.2[ovlp$group.1 == 1], overlap.to.first.week = if (is(ovlp,
            "try-error"))
            NA else ovlp$mean.overlap[ovlp$group.1 == 1], distance.to.first.week = if (is(dsts,
            "try-error"))
            NA else dsts$mean.distance[dsts$group.1 == 1], treatment = X$treatment[X$ID ==
            x][1], round = X$round[X$ID == x][1])
    } else out <- data.frame(ID = x, group = X$record.group[X$ID == x][1], week = NA,
        overlap.to.first.week = NA, distance.to.first.week = NA, treatment = X$treatment[X$ID ==
            x][1], round = X$round[X$ID == x][1])
    return(out)

})

indiv_ovlp <- do.call(rbind, indiv_ovlp_l)

saveRDS(indiv_ovlp, "./data/processed/acoustic_space_density_overlap_to_first_week_by_individual_TSNE.RDS")

indiv_ovlp <- readRDS("./data/processed/acoustic_space_density_overlap_to_first_week_by_individual_TSNE.RDS")

indiv_ovlp$treatment <- factor(indiv_ovlp$treatment, levels = c("Control", "Medium Stress",
    "High Stress"))
indiv_ovlp <- indiv_ovlp[!is.na(indiv_ovlp$group), ]

ggplot(data = indiv_ovlp, aes(x = week, y = overlap.to.first.week, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Individual overlap to acoustic space from first week",
    x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
    0.25))

ggplot(data = indiv_ovlp, aes(x = week, y = distance.to.first.week, group = ID, col = treatment)) +
    geom_point() + geom_line(size = 1.2) + scale_color_viridis_d(begin = 0.1, end = 0.9) +
    facet_wrap(~round, ncol = 2) + labs(y = "Individual distance to acoustic space from first week",
    x = "Week") + theme_classic(base_size = 30) + theme(legend.position = c(0.8,
    0.25))

agg_dist <- aggregate(distance.to.first.week ~ week + treatment, indiv_ovlp, mean)
agg_dist$se <- aggregate(distance.to.first.week ~ week + treatment, indiv_ovlp, standard_error)[,
    3]
# agg_stress$Week <- 1:4


# ggplot(data = agg_dist, aes(x = week, y = distance.to.first.week, fill =
# treatment)) + geom_bar(stat = 'identity') +
# geom_errorbar(aes(ymin=distance.to.first.week - se, ymax =
# distance.to.first.week + se), width=.1) + # geom_line(size = 1.2) +
# scale_fill_viridis_d(begin = 0.1, end = 0.9) + facet_wrap(~ treatment, ncol =
# 3, scale = 'fixed') + labs(y = 'Distance to first week', x = 'Week') +
# theme_classic(base_size = 30) + # scale_x_discrete( labels = c(1:4)) +
# theme(legend.position = 'none')

Statistical analyses

• on t-SNE dimensions
• predictors were mean-centered
• link functions: binomial distribution for overlap and gaussian for area and call count
• cort.stress and cort.baseline were dropped as they are highly correlated to stress response
• Not a lot of data makes having many predictors complicated
• No random factor (only 1 row per individual)
• No p values printed. Look at whether credibility intervals (l-95% CI & u-95% CI) include 0 to assess statistical significance

Overlap to self (week 1) as response

overlap.to.first.week ~ treatment + stress.response

sels$ID.week <- paste(sels$ID, sels$week, sep = "-")
indiv_ovlp$ID.week <- paste(indiv_ovlp$ID, indiv_ovlp$week, sep = "-")

# addd stress parameters
indiv_ovlp$cort.baseline <- sapply(indiv_ovlp$ID.week, function(x) sels$cort.baseline[sels$ID.week ==
    x][1])
indiv_ovlp$cort.stress <- sapply(indiv_ovlp$ID.week, function(x) sels$cort.stress[sels$ID.week ==
    x][1])
indiv_ovlp$breath.count <- sapply(indiv_ovlp$ID.week, function(x) sels$breath.count[sels$ID.week ==
    x][1])
indiv_ovlp$stress.response <- indiv_ovlp$cort.stress - indiv_ovlp$cort.baseline

# mean centering
indiv_ovlp$cort.baseline <- mean(indiv_ovlp$cort.baseline, na.rm = TRUE) - indiv_ovlp$cort.baseline
indiv_ovlp$cort.stress <- mean(indiv_ovlp$cort.stress, na.rm = TRUE) - indiv_ovlp$cort.stress
indiv_ovlp$breath.count <- mean(indiv_ovlp$breath.count, na.rm = TRUE) - indiv_ovlp$breath.count
indiv_ovlp$stress.response <- mean(indiv_ovlp$stress.response, na.rm = TRUE) - indiv_ovlp$stress.response


# only individuals present in week 1 indiv_ovlp <- indiv_ovlp[indiv_ovlp$ID
# %in% indiv_ovlp$ID[indiv_ovlp$week == 1], ]

# indiv_ovlp <- indiv_ovlp[complete.cases(indiv_ovlp), ]

indiv_ovlp$overlap.to.first.week[indiv_ovlp$overlap.to.first.week == 0] <- 1e-04
indiv_ovlp$overlap.to.first.week[indiv_ovlp$overlap.to.first.week == 1] <- 0.9999

week 2:

## week 2

overlap_week_2 <- brm(iter = 5000, silent = 2, overlap.to.first.week ~ treatment +
    stress.response, data = indiv_ovlp[indiv_ovlp$week == 2, ], family = Beta(),
    control = list(adapt_delta = 0.9), chains = 1)

saveRDS(overlap_week_2, "./data/processed/brms_model_overlap_week_2.RDS")

overlap_week_2 <- readRDS("./data/processed/brms_model_overlap_week_2.RDS")

summary_brm_model(overlap_week_2)

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
treatmentMediumStress	0.065	1	1527.558	5000	1	-0.723	0.798
treatmentHighStress	0.618	1	1453.076	5000	1	-0.131	1.334
stress.response	-0.006	1	2615.287	5000	1	-0.026	0.013

week 4:

## week 4
overlap_week_4 <- brm(iter = 5000, silent = 2, overlap.to.first.week ~ treatment +
    stress.response, data = indiv_ovlp[indiv_ovlp$week == 4, ], family = Beta(),
    control = list(adapt_delta = 0.9, max_treedepth = 12), chains = 1)

saveRDS(overlap_week_4, "./data/processed/brms_model_overlap_week_4.RDS")

overlap_week_4 <- readRDS("./data/processed/brms_model_overlap_week_4.RDS")

summary_brm_model(overlap_week_4)

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
treatmentMediumStress	-0.702	1	2257.376	5000	1	-1.971	0.613
treatmentHighStress	0.454	1	2095.748	5000	1	-0.639	1.556
stress.response	0.042	1	2357.324	5000	1	0.001	0.082

Acoustic space change (compare to week 1)

area ~ treatment + stress.response

areas_by_week$ID.week <- paste(areas_by_week$ID, areas_by_week$week, sep = "-")


# addd stress parameters
areas_by_week$cort.baseline <- sapply(areas_by_week$ID.week, function(x) sels$cort.baseline[sels$ID.week ==
    x][1])
areas_by_week$cort.stress <- sapply(areas_by_week$ID.week, function(x) sels$cort.stress[sels$ID.week ==
    x][1])
areas_by_week$breath.count <- sapply(areas_by_week$ID.week, function(x) sels$breath.count[sels$ID.week ==
    x][1])
areas_by_week$stress.response <- areas_by_week$cort.stress - areas_by_week$cort.baseline

# mean centering
areas_by_week$cort.baseline <- mean(areas_by_week$cort.baseline, na.rm = TRUE) -
    areas_by_week$cort.baseline
areas_by_week$cort.stress <- mean(areas_by_week$cort.stress, na.rm = TRUE) - areas_by_week$cort.stress
areas_by_week$breath.count <- mean(areas_by_week$breath.count, na.rm = TRUE) - areas_by_week$breath.count
areas_by_week$stress.response <- mean(areas_by_week$stress.response, na.rm = TRUE) -
    areas_by_week$stress.response


# only individuals present in week 1 areas_by_week <-
# areas_by_week[areas_by_week$ID %in% areas_by_week$ID[areas_by_week$week ==
# 1], ]


# areas_by_week <- areas_by_week[complete.cases(areas_by_week), ]

week 2:

## week 2
area_week_2 <- brm(iter = 5000, silent = 2, size ~ treatment + stress.response, data = areas_by_week[areas_by_week$week ==
    2, ], control = list(adapt_delta = 0.9), chains = 1)

saveRDS(area_week_2, "./data/processed/brms_model_area_week_2.RDS")

area_week_2 <- readRDS("./data/processed/brms_model_area_week_2.RDS")

summary_brm_model(area_week_2)

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
treatmentMediumStress	0.175	1	1787.606	5000	1	-0.314	0.672
treatmentHighStress	0.126	1	1701.261	5000	1	-0.384	0.614
stress.response	0.012	1	2390.428	5000	1	-0.002	0.026

week 4:

## week 4
area_week_4 <- brm(iter = 5000, silent = 2, size ~ treatment + stress.response, data = areas_by_week[areas_by_week$week ==
    4, ], control = list(adapt_delta = 0.9), chains = 1)

saveRDS(area_week_4, "./data/processed/brms_model_area_week_4.RDS")

area_week_4 <- readRDS("./data/processed/brms_model_area_week_4.RDS")

summary_brm_model(area_week_4)

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
treatmentMediumStress	0.461	1.004	1131.119	5000	1	-0.335	1.255
treatmentHighStress	0.669	1.006	1362.654	5000	1	-0.054	1.394
stress.response	0.017	1.001	1806.867	5000	1	-0.013	0.045

Call count

week 2:

## week 2
count_week_2 <- brm(iter = 5000, silent = 2, selec ~ treatment + stress.response,
    data = call_count[call_count$week == 2, ], control = list(adapt_delta = 0.9),
    chains = 1)

saveRDS(count_week_2, "./data/processed/brms_model_count_week_2.RDS")

count_week_2 <- readRDS("./data/processed/brms_model_count_week_2.RDS")

summary_brm_model(count_week_2)

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
treatmentMediumStress	355.784	1.003	1074.230	5000	1	-127.180	804.368
treatmentHighStress	613.795	1.001	1050.304	5000	1	131.824	1086.057
stress.response	0.971	1	1396.041	5000	1	-9.032	10.761

week 4:

## week 4
count_week_4 <- brm(iter = 5000, silent = 2, selec ~ treatment + stress.response,
    data = call_count[call_count$week == 4, ], control = list(adapt_delta = 0.9),
    chains = 1)

saveRDS(count_week_4, "./data/processed/brms_model_count_week_4.RDS")

count_week_4 <- readRDS("./data/processed/brms_model_count_week_4.RDS")

summary_brm_model(count_week_4)

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
treatmentMediumStress	-103.093	1	1001.866	5000	1	-335.787	123.953
treatmentHighStress	3.989	1.001	840.242	5000	1	-352.231	359.975
stress.response	2.278	1.001	897.594	5000	1	-1.867	5.900

 

Takeaways

• Stress has an effect on number of calls: high-stress individuals produce more calls than low-stress individuals in week 2
• Stress response has a positive effect on the overlap to week 1 (overlap to itself)

 

---

Proposed Models for Stress and Learning Pilot Data based on 9/28/21 discussion

Model set 1: effects of treatment

Predicted variables: - Physiology: Baseline CORT, stress response CORT (stress), difference baseline-stress response CORT (stress response), weight, breath rate - Learning: number of calls, total acoustic area, distance moved from week 1, overlap to self week 1, overlap self to group

General form: - Predicted ~ treatment + week + ind(rand) - Predicted ~ treatment + week + treat * week + ind(rand) - Predicted ~ round + week [to check for effect of round]

Note: Use week 1 as baseline, comparing weeks 2, 3, 4

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.area <- 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$acoustic.overlap <- 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]]

    if (length(overlap) < 1)
        overlap <- NA

    return(overlap)
})

agg_dat$acoustic.overlap.to.group <- 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]]))

# remove week 1
agg_dat <- agg_dat[agg_dat$week != 1, ]


responses <- setdiff(names(agg_dat), c("ID", "week", "treatment", "round"))

predictors <- c("~ treatment + week + (1|ID)", "~ treatment * week + (1|ID)", "~ round + week + (1|ID)")

formulas <- expand.grid(responses = responses, predictors = predictors, stringsAsFactors = FALSE)

vars_to_scale <- names(agg_dat)[sapply(agg_dat, is.numeric)]
vars_to_scale <- vars_to_scale[!vars_to_scale %in% c("week", "round")]

for (i in vars_to_scale) agg_dat[, vars_to_scale] <- scale(agg_dat[, vars_to_scale])

treatment_models <- lapply(1:nrow(formulas), function(x) {

    sub_dat <- agg_dat[!is.na(agg_dat[names(agg_dat) == formulas$responses[x]]),
        ]
    sub_dat

    mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
        iter = 5000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
        chains = 2), silent = TRUE)
    return(mod)
})


names(treatment_models) <- paste(formulas$responses, formulas$predictors)

saveRDS(treatment_models, "./data/processed/treatment_models.RDS")

agg_dat$week <- as.factor(agg_dat$week)

treatment_models_week_factor <- lapply(1:nrow(formulas), function(x) {

    sub_dat <- agg_dat[!is.na(agg_dat[names(agg_dat) == formulas$responses[x]]),
        ]

    mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
        iter = 5000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
        chains = 2), silent = TRUE)
    return(mod)
})

names(treatment_models_week_factor) <- paste(formulas$responses, formulas$predictors)

saveRDS(treatment_models_week_factor, "./data/processed/treatment_models_week_factor.RDS")

week as continuous

treatment_models <- readRDS("./data/processed/treatment_models.RDS")



for (x in treatment_models) {
    summary_brm_model(x)
}

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
baseline.CORT	-0.342	1.001	2032.331	5000	2	-0.623	-0.047

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
baseline.CORT	0.057	1	3386.12	5000	2	-0.198	0.305

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
baseline.CORT	-0.08	1.001	2172.071	5000	2	-0.368	0.209

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
baseline.CORT	0.018	1.001	3184.54	5000	2	-0.28	0.314

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
baseline.CORT	0.062	1	2669.238	5000	2	-0.227	0.366

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
stress.response	0.267	1.001	1967.507	5000	2	-0.04	0.569

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
stress.response	-0.14	1	3307.176	5000	2	-0.397	0.111

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
stress.response	0.067	1	2819.069	5000	2	-0.213	0.339

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
stress.response	-0.089	1	2859.569	5000	2	-0.384	0.212

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
stress.response	0.001	1	3164.738	5000	2	-0.309	0.302

 

Takeaways

• No significant nor marginal effects

 

Model set 2: individual stress and effects on learning

General form: - Learning variables ~ baseline CORT [Week 2 immediate stress and learning, week 3 chronic stress and learn]

agg_dat <- agg_dat[agg_dat$week != "4", ]

agg_dat$week <- as.factor(agg_dat$week)


responses <- c("call.count", "acoustic.area", "acoustic.distance", "acoustic.overlap",
    "acoustic.overlap.to.group")

predictors <- c("~ baseline.CORT + week + (1|ID)", "~ stress.response + week + (1|ID)")

formulas <- expand.grid(responses = responses, predictors = predictors, stringsAsFactors = FALSE)

treatment_models <- lapply(1:nrow(formulas), function(x) {

    sub_dat <- agg_dat[!is.na(agg_dat[names(agg_dat) == formulas$responses[x]]),
        ]

    mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
        iter = 5000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
        chains = 1), silent = TRUE)
    return(mod)
})


names(treatment_models) <- paste(formulas$responses, formulas$predictors)

saveRDS(treatment_models, "./data/processed/stress_on_learning_models.RDS")

treatment_models <- readRDS("./data/processed/stress_on_learning_models.RDS")



for (x in treatment_models) {
    summary_brm_model(x)
}

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
baseline.CORT	-0.192	1.002	997.111	5000	1	-0.546	0.177

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
baseline.CORT	-0.154	1	537.359	5000	1	-0.462	0.159

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
baseline.CORT	0.007	1.015	256.643	5000	1	-0.248	0.274

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
baseline.CORT	0.101	1	616.351	5000	1	-0.247	0.44

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
stress.response	0.299	1.001	1394.371	5000	1	-0.112	0.69

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
stress.response	0.013	1	556.089	5000	1	-0.25	0.287

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
stress.response	-0.058	1.001	677.164	5000	1	-0.391	0.286

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
stress.response	-0.138	1.002	869.593	5000	1	-0.49	0.209

 

Takeaways

• No significant effects

 

Model set 3: individual stress and effects on learning on week 2 and 3 independently

Week 2 immediate stress and learning, week 3 chronic stress and learn

General form: - Learning variables ~ baseline CORT - Learning variables ~ difference baseline-stress CORT

Week 2

agg_dat_w2 <- agg_dat[agg_dat$week == "2", ]

responses <- c("call.count", "acoustic.area", "acoustic.distance", "acoustic.overlap",
    "acoustic.overlap.to.group")

predictors <- c("~ baseline.CORT + (1|ID)", "~ stress.response + (1|ID)")

formulas <- expand.grid(responses = responses, predictors = predictors, stringsAsFactors = FALSE)

treatment_models <- lapply(1:nrow(formulas), function(x) {

    sub_dat <- agg_dat_w2[!is.na(agg_dat_w2[names(agg_dat_w2) == formulas$responses[x]]),
        ]

    mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
        iter = 5000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
        chains = 1), silent = TRUE)
    return(mod)
})


names(treatment_models) <- paste(formulas$responses, formulas$predictors)

saveRDS(treatment_models, "./data/processed/stress_on_learning_models_week_2.RDS")

treatment_models <- readRDS("./data/processed/stress_on_learning_models_week_2.RDS")

for (x in treatment_models) {
    summary_brm_model(x)
}

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
baseline.CORT	-0.177	1	1201.316	5000	1	-0.501	0.164

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
baseline.CORT	0.037	1.002	701.318	5000	1	-0.328	0.408

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
baseline.CORT	-0.138	1.012	135.412	5000	1	-0.426	0.133

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
baseline.CORT	-0.03	1.045	78.7	5000	1	-0.329	0.298

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
baseline.CORT	0.11	1.004	1054.756	5000	1	-0.257	0.458

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
stress.response	0.261	1.001	831.794	5000	1	-0.072	0.586

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
stress.response	-0.178	1	227.403	5000	1	-0.545	0.202

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
stress.response	0.013	1	1041.844	5000	1	-0.232	0.25

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
stress.response	-0.045	1.002	546.832	5000	1	-0.381	0.267

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
stress.response	-0.138	1.004	1319.857	5000	1	-0.529	0.261

Week 3

agg_dat_w3 <- agg_dat[agg_dat$week == "3", ]

responses <- c("call.count", "acoustic.area", "acoustic.distance", "acoustic.overlap",
    "acoustic.overlap.to.group")

predictors <- c("~ baseline.CORT + (1|ID)", "~ stress.response + (1|ID)")

formulas <- expand.grid(responses = responses, predictors = predictors, stringsAsFactors = FALSE)

treatment_models <- lapply(1:nrow(formulas), function(x) {

    sub_dat <- agg_dat_w3[!is.na(agg_dat_w3[names(agg_dat_w3) == formulas$responses[x]]),
        ]

    mod <- try(brm(formula = paste(formulas$responses[x], formulas$predictors[x]),
        iter = 5000, silent = 2, data = sub_dat, control = list(adapt_delta = 0.9),
        chains = 1), silent = TRUE)
    return(mod)
})


names(treatment_models) <- paste(formulas$responses, formulas$predictors)

saveRDS(treatment_models, "./data/processed/stress_on_learning_models_week_3.RDS")

treatment_models <- readRDS("./data/processed/stress_on_learning_models_week_3.RDS")

for (x in treatment_models) {
    summary_brm_model(x)
}

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
baseline.CORT	-0.192	1.002	997.111	5000	1	-0.546	0.177

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
baseline.CORT	-0.154	1	537.359	5000	1	-0.462	0.159

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
baseline.CORT	0.007	1.015	256.643	5000	1	-0.248	0.274

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
baseline.CORT	0.101	1	616.351	5000	1	-0.247	0.44

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
stress.response	0.299	1.001	1394.371	5000	1	-0.112	0.69

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
stress.response	0.013	1	556.089	5000	1	-0.25	0.287

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
stress.response	-0.058	1.001	677.164	5000	1	-0.391	0.286

	Estimate	Rhat	Bulk_ESS	iterations	chains	CI_low	CI_high
stress.response	-0.138	1.002	869.593	5000	1	-0.49	0.209

 

Takeaways

• No significant nor marginal effects

 

Individual models for each acoustic feature

# spectrographic parameters
responses <- c("duration", "meanfreq", "sd", "freq.median", "freq.IQR", "time.IQR",
    "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom", "mindom", "maxdom",
    "dfrange", "modindx", "meanpeakf")

frmla <- as.formula(paste("cbind(", paste(responses, collapse = ", "), ") ~ ID + week + round + treatment"))

agg_dat <- aggregate(x = frmla, data = sels, FUN = mean)

predictors <- c("~ treatment + week + (1|ID) + (1|round)")

formulas <- paste(responses, predictors)

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])

acous_models <- lapply(1:length(formulas), function(x) {

    sub_dat <- sub_agg_dat[!is.na(sub_agg_dat[names(sub_agg_dat) == responses[x]]),
        ]
    # sub_dat

    mod <- brm(formula = formulas[x], iter = 20000, silent = 2, data = sub_agg_dat,
        control = list(adapt_delta = 0.9), chains = 4, 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")), file = paste0("./data/processed/",
            responses[x], "_model"), file_refit = "on_change")

    return(mod)
})

# spectrographic parameters
responses <- c("duration", "meanfreq", "sd", "freq.median", "freq.IQR", "time.IQR",
    "skew", "kurt", "sp.ent", "time.ent", "entropy", "meandom", "mindom", "maxdom",
    "dfrange", "modindx", "meanpeakf")


for (x in responses) {
    extended_summary(read.file = paste0("./data/processed/", x, "_model.rds"))
}

duration_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10)	duration ~ treatment + week + (1 | ID) + (1 | round)	20000	4	1	10000	1276	0	329.889	86.68	1562788543

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	0.144	-2.197	1.103	1.019	329.889	86.68
b_treatmentMediumStress	-0.251	-0.941	0.439	1.001	13072.859	20001.96
b_treatmentHighStress	-0.363	-1.049	0.310	1.002	12586.678	18130.57
b_week	-0.121	-0.281	0.038	1.002	32974.887	15669.85

meanfreq_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10)	meanfreq ~ treatment + week + (1 | ID) + (1 | round)	20000	4	1	10000	1799	0	323.95	83.014	370950808

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	0.184	-1.284	1.106	1.013	323.950	83.014
b_treatmentMediumStress	-0.525	-1.255	0.196	1.002	5773.121	12671.379
b_treatmentHighStress	-0.216	-0.937	0.511	1.002	6624.105	13789.724
b_week	-0.036	-0.167	0.097	1.001	8588.690	11420.647

sd_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10)	sd ~ treatment + week + (1 | ID) + (1 | round)	20000	4	1	10000	433	0	8020.574	7841.084	1815845249

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	0.149	-0.856	1.166	1.001	8020.574	8032.576
b_treatmentMediumStress	-0.232	-0.954	0.482	1	11617.406	7841.084
b_treatmentHighStress	-0.134	-0.858	0.594	1	12050.008	9822.115
b_week	0.101	-0.055	0.256	1	35129.504	25461.028

freq.median_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10)	freq.median ~ treatment + week + (1 | ID) + (1 | round)	20000	4	1	10000	731	0	1955.866	520.189	1791416949

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	0.218	-0.606	1.125	1.004	1955.866	520.189
b_treatmentMediumStress	-0.494	-1.196	0.210	1	9637.794	14689.933
b_treatmentHighStress	-0.136	-0.835	0.563	1.001	10174.628	14158.669
b_week	-0.087	-0.234	0.064	1.001	10271.898	3480.762

freq.IQR_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10)	freq.IQR ~ treatment + week + (1 | ID) + (1 | round)	20000	4	1	10000	1594	0	293.519	73.825	1169392167

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	0.326	-1.341	1.403	1.017	293.519	73.825
b_treatmentMediumStress	-0.529	-1.233	0.186	1.002	11953.607	18575.794
b_treatmentHighStress	-0.491	-1.199	0.226	1.001	9847.216	13624.883
b_week	0.076	-0.078	0.228	1.001	23453.079	24347.250

time.IQR_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10)	time.IQR ~ treatment + week + (1 | ID) + (1 | round)	20000	4	1	10000	1048	0	815.082	206.517	685652361

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	0.131	-0.679	0.955	1.003	5530.109	3419.626
b_treatmentMediumStress	-0.055	-0.740	0.641	1.001	8811.131	14717.701
b_treatmentHighStress	-0.315	-1.006	0.390	1.002	3942.247	9511.849
b_week	-0.146	-0.334	0.012	1.005	815.082	206.517

skew_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10)	skew ~ treatment + week + (1 | ID) + (1 | round)	20000	4	1	10000	476	0	2897.495	1552.042	1565054318

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-0.020	-0.856	0.876	1.001	2897.495	1552.042
b_treatmentMediumStress	-0.015	-0.740	0.709	1	13629.868	19828.604
b_treatmentHighStress	0.112	-0.607	0.829	1	13127.974	20827.895
b_week	-0.087	-0.239	0.065	1	39466.763	27254.629

kurt_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10)	kurt ~ treatment + week + (1 | ID) + (1 | round)	20000	4	1	10000	1472	0	214.634	54.007	1042666801

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	0.095	-0.734	1.780	1.024	214.634	54.007
b_treatmentMediumStress	-0.083	-0.792	0.626	1.004	9361.053	9567.585
b_treatmentHighStress	0.047	-0.665	0.760	1.001	6902.980	9999.986
b_week	-0.081	-0.239	0.074	1.001	8216.088	18059.358

sp.ent_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10)	sp.ent ~ treatment + week + (1 | ID) + (1 | round)	20000	4	1	10000	553	0	4756.196	2661.672	1463071904

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	0.156	-0.828	1.072	1.001	4756.196	2661.672
b_treatmentMediumStress	-0.304	-1.082	0.470	1.001	7376.602	8769.636
b_treatmentHighStress	-0.226	-0.993	0.546	1	9758.257	15131.749
b_week	0.030	-0.107	0.167	1	36625.948	25446.555

time.ent_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10)	time.ent ~ treatment + week + (1 | ID) + (1 | round)	20000	4	1	10000	1432	0	528.473	214.778	1451168900

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-0.255	-1.337	0.604	1.007	528.473	214.778
b_treatmentMediumStress	0.275	-0.445	1.027	1	4754.476	8405.619
b_treatmentHighStress	0.367	-0.360	1.085	1.002	7887.367	16317.369
b_week	0.101	-0.031	0.236	1.002	2475.276	14457.123

entropy_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10)	entropy ~ treatment + week + (1 | ID) + (1 | round)	20000	4	1	10000	550	0	6105.411	3347.197	1366523000

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	0.084	-0.815	0.996	1.001	6105.411	3347.197
b_treatmentMediumStress	-0.207	-0.986	0.585	1.001	8284.506	14634.871
b_treatmentHighStress	-0.116	-0.908	0.687	1.001	7477.031	4962.430
b_week	0.061	-0.062	0.185	1.001	28764.137	24297.597

meandom_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10)	meandom ~ treatment + week + (1 | ID) + (1 | round)	20000	4	1	10000	837	0	2277.245	819.184	2111900018

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-0.039	-0.953	0.809	1.002	2277.245	819.184
b_treatmentMediumStress	-0.087	-0.850	0.689	1	5270.850	5265.297
b_treatmentHighStress	0.117	-0.653	0.879	1	7529.964	14256.442
b_week	-0.014	-0.139	0.113	1.001	15571.733	8330.339

mindom_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10)	mindom ~ treatment + week + (1 | ID) + (1 | round)	20000	4	1	10000	1438	0	650.372	406.045	2129214010

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-0.011	-1.160	1.191	1.009	650.372	406.045
b_treatmentMediumStress	-0.118	-0.911	0.654	1.001	5951.474	12470.824
b_treatmentHighStress	0.130	-0.667	0.905	1.001	5237.840	4040.680
b_week	-0.009	-0.133	0.114	1.002	2461.532	811.614

maxdom_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10)	maxdom ~ treatment + week + (1 | ID) + (1 | round)	20000	4	1	10000	913	0	751.599	221.718	1966969513

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	0.204	-0.900	1.617	1.008	751.599	221.718
b_treatmentMediumStress	-0.174	-0.893	0.556	1.002	2775.468	4470.150
b_treatmentHighStress	-0.325	-1.046	0.407	1.002	2137.051	1319.135
b_week	0.058	-0.099	0.215	1.004	1484.800	634.757

dfrange_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10)	dfrange ~ treatment + week + (1 | ID) + (1 | round)	20000	4	1	10000	686	0	3461.257	1845.981	1303591726

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	0.105	-0.764	0.936	1.001	3461.257	1845.981
b_treatmentMediumStress	-0.015	-0.823	0.784	1	8443.265	9377.307
b_treatmentHighStress	-0.273	-1.049	0.534	1	7634.496	9244.040
b_week	0.046	-0.089	0.183	1	22376.728	6598.531

modindx_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10)	modindx ~ treatment + week + (1 | ID) + (1 | round)	20000	4	1	10000	722	0	4334.738	3227.693	287973161

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	0.085	-0.910	1.113	1.001	4334.738	3227.693
b_treatmentMediumStress	-0.163	-0.907	0.590	1	11115.663	15768.107
b_treatmentHighStress	-0.153	-0.901	0.595	1	11515.166	14705.783
b_week	-0.112	-0.246	0.026	1	19928.754	8393.945

meanpeakf_model

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 5) Intercept-normal(0, 10) sd-student_t(3, 0, 10) sigma-student_t(3, 0, 10)	meanpeakf ~ treatment + week + (1 | ID) + (1 | round)	20000	4	1	10000	2437	0	66.176	20.647	2062562242

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_Intercept	-0.254	-1.689	0.524	1.043	66.176	20.647
b_treatmentMediumStress	0.067	-0.574	0.700	1.008	1615.619	7377.159
b_treatmentHighStress	0.423	-0.209	1.040	1.005	1989.396	7704.721
b_week	-0.059	-0.230	0.118	1.005	820.734	2145.340

---

Session information

## R version 4.3.2 (2023-10-31)
## 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;  LAPACK version 3.9.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=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       
## 
## time zone: America/Costa_Rica
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] brmsish_1.0.0          PhenotypeSpace_0.1.0   ggridges_0.5.6        
##  [4] posterior_1.5.0        ggrepel_0.9.5          cowplot_1.1.3         
##  [7] tidybayes_3.0.6        modelr_0.1.11          tidyr_1.3.1           
## [10] forcats_1.0.0          purrr_1.0.2            dplyr_1.1.4           
## [13] lme4_1.1-35.1          Matrix_1.6-5           brms_2.20.4           
## [16] Rcpp_1.0.12            vegan_2.6-4            lattice_0.22-5        
## [19] permute_0.9-7          raster_3.6-26          spatstat_3.0-7        
## [22] spatstat.linnet_3.1-4  spatstat.model_3.2-10  rpart_4.1.23          
## [25] spatstat.explore_3.2-6 nlme_3.1-164           spatstat.random_3.2-2 
## [28] spatstat.geom_3.2-8    spatstat.data_3.0-4    GGally_2.2.1          
## [31] adehabitatHR_0.4.21    adehabitatLT_0.3.27    CircStats_0.2-6       
## [34] boot_1.3-28.1          adehabitatMA_0.3.16    ade4_1.7-22           
## [37] sp_2.1-3               MASS_7.3-60.0.1        Rtsne_0.17            
## [40] randomForest_4.7-1.1   formatR_1.14           knitr_1.45            
## [43] kableExtra_1.4.0       ggplot2_3.5.0          viridis_0.6.5         
## [46] viridisLite_0.4.2      pbapply_1.7-2          readxl_1.4.3          
## [49] lubridate_1.9.3        remotes_2.4.2.1       
## 
## loaded via a namespace (and not attached):
##   [1] svUnit_1.0.6          shinythemes_1.2.0     splines_4.3.2        
##   [4] later_1.3.2           tibble_3.2.1          cellranger_1.1.0     
##   [7] polyclip_1.10-6       xts_0.13.2            lifecycle_1.0.4      
##  [10] StanHeaders_2.32.5    crosstalk_1.2.1       ggdist_3.3.1         
##  [13] backports_1.4.1       magrittr_2.0.3        sass_0.4.8           
##  [16] rmarkdown_2.25        jquerylib_0.1.4       yaml_2.3.8           
##  [19] httpuv_1.6.14         spatstat.core_2.4-4   pkgbuild_1.4.3       
##  [22] spatstat.sparse_3.0-3 rgeos_0.6-4           minqa_1.2.6          
##  [25] RColorBrewer_1.1-3    multcomp_1.4-25       abind_1.4-5          
##  [28] TH.data_1.1-2         tensorA_0.36.2.1      sandwich_3.1-0       
##  [31] inline_0.3.19         spatstat.utils_3.0-4  terra_1.7-71         
##  [34] goftest_1.2-3         bridgesampling_1.1-2  svglite_2.1.3        
##  [37] codetools_0.2-19      DT_0.31               xml2_1.3.6           
##  [40] tidyselect_1.2.0      bayesplot_1.11.0      farver_2.1.1         
##  [43] matrixStats_1.2.0     stats4_4.3.2          base64enc_0.1-3      
##  [46] jsonlite_1.8.8        ellipsis_0.3.2        survival_3.5-7       
##  [49] emmeans_1.10.0        systemfonts_1.0.5     tools_4.3.2          
##  [52] glue_1.7.0            gridExtra_2.3         xfun_0.42            
##  [55] mgcv_1.9-1            distributional_0.3.2  loo_2.6.0            
##  [58] withr_3.0.0           fastmap_1.1.1         fansi_1.0.6          
##  [61] shinyjs_2.1.0         digest_0.6.34         timechange_0.3.0     
##  [64] R6_2.5.1              mime_0.12             estimability_1.4.1   
##  [67] colorspace_2.1-0      gtools_3.9.5          tensor_1.5           
##  [70] markdown_1.12         threejs_0.3.3         utf8_1.2.4           
##  [73] generics_0.1.3        htmlwidgets_1.6.4     ggstats_0.5.1        
##  [76] pkgconfig_2.0.3       dygraphs_1.1.1.6      gtable_0.3.4         
##  [79] htmltools_0.5.7       scales_1.3.0          rstudioapi_0.15.0    
##  [82] reshape2_1.4.4        coda_0.19-4.1         checkmate_2.3.1      
##  [85] nloptr_2.0.3          proxy_0.4-27          cachem_1.0.8         
##  [88] zoo_1.8-12            stringr_1.5.1         parallel_4.3.2       
##  [91] miniUI_0.1.1.1        pillar_1.9.0          grid_4.3.2           
##  [94] vctrs_0.6.5           shinystan_2.6.0       promises_1.2.1       
##  [97] arrayhelpers_1.1-0    xtable_1.8-4          cluster_2.1.6        
## [100] evaluate_0.23         mvtnorm_1.2-4         cli_3.6.2            
## [103] compiler_4.3.2        rlang_1.1.3           rstantools_2.4.0     
## [106] labeling_0.4.3        plyr_1.8.9            stringi_1.8.3        
## [109] rstan_2.32.5          deldir_2.0-2          QuickJSR_1.1.3       
## [112] HDInterval_0.2.4      munsell_0.5.0         colourpicker_1.3.0   
## [115] Brobdingnag_1.2-9     shiny_1.8.0           highr_0.10           
## [118] broom_1.0.5           igraph_1.6.0          RcppParallel_5.1.7   
## [121] bslib_0.6.1           ape_5.7-1