--- title: Statistical analysis subtitle: Chronic stress raises and vocal plasticity in budgerigars author: <a href="marce10.github.io">Marcelo Araya-Salas, PhD</a> & <a href="https://wrightbehaviorlab.org/">Timothy Wright, PhD</a> date: "`r Sys.Date()`" toc: true toc-depth: 2 toc-location: left number-sections: true highlight-style: pygments format: html: df-print: kable code-fold: true code-tools: true css: qmd.css --- < style > counter-reset : source-line 0 ; }.numberSource codecounter-reset : none ; }</ style > ```{r, echo = FALSE} # set working directory as project directory or one directory above, <- try (rprojroot:: find_rstudio_root_file (), silent = TRUE )if (is (rootdir, "try-error" )) rootdir <- ".." :: opts_knit$ set (root.dir = rootdir)``` <!-- skyblue box --> < div class = "alert alert-info" > # Data analysis for the paper {.unnumbered .unlisted} </ div > ```{r add link to github repo if any, echo = FALSE, results='asis'} # print link to github repo if any if (file.exists ("./.git/config" )){<- readLines ("./.git/config" )<- grep ("url" , config, value = TRUE )<- gsub (" \\ turl = |.git$" , "" , url)cat (" \n All data and source code is also found at: [" , url, "](" , url, ")" , sep = "" )``` ```{r setup style, echo = FALSE, message = FALSE, warning=FALSE} # options to customize chunk outputs :: opts_chunk$ set (class.source = "numberLines lineAnchors" , # for code line numbers tidy.opts = list (width.cutoff = 65 ), tidy = TRUE ,message = FALSE ``` ## Load packages ```{r packages, message = FALSE, warning = FALSE, echo = TRUE, eval = TRUE} # install "sketchy" if not install yet if (! requireNamespace ("sketchy" , quietly = TRUE )) {install.packages ("sketchy" )## add 'developer' to packages to be installed from github <- c ("remotes" ,"readxl" ,"pbapply" ,"viridis" ,"ggplot2" ,"kableExtra" ,"knitr" ,"formatR" ,"MASS" ,"brms" ,"lme4" ,"cowplot" ,"maRce10/PhenotypeSpace" ,"maRce10/brmsish" :: load_packages (x)# funnction to print tables in html format <- function (x) {<- kable (x, row.names = TRUE , digits = 4 , "html" )<- kable_styling (kb,bootstrap_options = c ("striped" , "hover" , "condensed" , "responsive" ))scroll_box (kb, width = "100%" )``` ## Functions and global parameters ```{r functions and parameters, eval = TRUE, echo = TRUE} $ set (root.dir = ".." )# set evaluation false $ set (fig.width = 10 ,fig.height = 6 ,warning = FALSE ,message = FALSE ,tidy = TRUE <- function (...)data.frame (read_excel (...))# for reading months in english format <- Sys.setlocale (locale = "en_US.UTF-8" )<- function (x)sd (x) / sqrt (length (x))<- viridis (10 , alpha = 0.7 )<- cols[7 ]``` ```{r paths, eval = TRUE, echo = TRUE} <- "./data/raw/" <- "./data/processed/" <- "./output" ``` ```{r read detections and prepare data, eval = TRUE} # read ext sel tab calls <- read.csv (file.path (raw_data_path, "acoustic_features_budgie_calls.csv" ))# keep only spectrographic parameters <- 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" $ ID <- sapply (sels$ sound.files, function (x)strsplit (x, "_" )[[1 ]][1 ])$ month <- sapply (sels$ sound.files, function (x)strsplit (x, "_" )[[1 ]][2 ])$ day <- sapply (sels$ sound.files, function (x)strsplit (x, "_" )[[1 ]][3 ])$ year <- sapply (sels$ sound.files, function (x)strsplit (x, "_" )[[1 ]][4 ])$ date <- paste (sels$ day, substr (sels$ month, 0 , 3 ), sels$ year, sep = "-" )$ date <- as.Date (sels$ date, format = "%d-%b-%Y" )# acoustic measurements <- read.csv (file.path (processed_data_path, "acoustic_space_area_by_id_and_week.csv" ))<- read.csv (file.path (processed_data_path,"acoustic_space_density_overlap_to_first_week_by_id.csv" ))$ treatment <- factor (indiv_ovlp$ treatment,levels = c ("Control" , "Medium Stress" , "High Stress" ))<- read.csv (file.path (processed_data_path, "acoustic_space_density_overlap_to_group_by_week_id.csv" ))$ treatment <- factor (group_ovlp$ treatment,levels = c ("Control" , "Medium Stress" , "High Stress" ))``` # Format data ```{r, eval = TRUE} <- read_excel_df (file.path (raw_data_path, "experiment_metadata.xlsx" ))$ ID[sels$ ID == "125YGMM" ] <- "125YGHM" $ ID[sels$ ID == "394YBHM" ] <- "394WBHM" $ treatment <- sapply (1 : nrow (sels), function (x) {$ Treatment[metadat$ Bird.ID == sels$ ID[x]][1 ]$ treatment.room <- sapply (1 : nrow (sels), function (x) {$ Treatment.Room[metadat$ Bird.ID == sels$ ID[x]][1 ]$ round <- sapply (1 : nrow (sels), function (x) {$ Round[metadat$ Bird.ID == sels$ ID[x]][1 ]$ source.room <- sapply (1 : nrow (sels), function (x) {$ Source.Room[metadat$ Bird.ID == sels$ ID[x]][1 ]$ record.group <- sapply (1 : nrow (sels), function (x) {$ Record.Group[metadat$ Bird.ID == sels$ ID[x]][1 ]# add week <- lapply (unique (sels$ round), function (x) {<- sels[sels$ round == x, ]<- min (Y$ date)<- min_date + seq (0 , 100 , by = 7 )$ week <- NA for (i in 2 : length (week_limits))$ week[Y$ date >= week_limits[i - 1 ] & Y$ date < <- i - 1 return (Y)<- do.call (rbind, out)$ cort.baseline <- sapply (1 : nrow (sels), function (x) {if (sels$ week[x] == 1 )<- metadat$ D3.CORT.Baseline[metadat$ Bird.ID == sels$ ID[x]][1 ]if (sels$ week[x] == 2 )<- metadat$ D7.CORT.Baseline[metadat$ Bird.ID == sels$ ID[x]][1 ]if (sels$ week[x] == 3 )<- metadat$ D14.CORT.Baseline[metadat$ Bird.ID == sels$ ID[x]][1 ]if (sels$ week[x] == 4 )<- metadat$ D21.CORT.Baseline[metadat$ Bird.ID == sels$ ID[x]][1 ]if (sels$ week[x] == 5 )<- metadat$ D28.CORT.Baseline[metadat$ Bird.ID == sels$ ID[x]][1 ]return (out)$ cort.stress <- sapply (1 : nrow (sels), function (x) {if (sels$ week[x] == 1 )<- metadat$ D3.CORT.Stress[metadat$ Bird.ID == sels$ ID[x]][1 ]if (sels$ week[x] == 2 )<- metadat$ D7.CORT.Stress[metadat$ Bird.ID == sels$ ID[x]][1 ]if (sels$ week[x] == 3 )<- metadat$ D14.CORT.Stress[metadat$ Bird.ID == sels$ ID[x]][1 ]if (sels$ week[x] == 4 )<- metadat$ D21.CORT.Stress[metadat$ Bird.ID == sels$ ID[x]][1 ]if (sels$ week[x] == 5 )<- metadat$ D28.CORT.Stress[metadat$ Bird.ID == sels$ ID[x]][1 ]return (out)$ stress.response <- sels$ cort.stress #- sels$cort.baseline $ weight <- sapply (1 : nrow (sels), function (x) {if (sels$ week[x] == 1 )<- metadat$ D3.Bird.Weight..g.[metadat$ Bird.ID == sels$ ID[x]][1 ]if (sels$ week[x] == 2 )<- metadat$ D7.Bird.Weight..g.[metadat$ Bird.ID == sels$ ID[x]][1 ]if (sels$ week[x] == 3 )<- metadat$ D14.Bird.Weight..g.[metadat$ Bird.ID == sels$ ID[x]][1 ]if (sels$ week[x] == 4 )<- metadat$ D21.Bird.Weight..g.[metadat$ Bird.ID == sels$ ID[x]][1 ]if (sels$ week[x] == 5 )<- metadat$ D28.Bird.Weight..g.[metadat$ Bird.ID == sels$ ID[x]][1 ]return (out)$ breath.count <- sapply (1 : nrow (sels), function (x) {if (sels$ week[x] == 1 )<- metadat$ D3.Breath.Count[metadat$ Bird.ID == sels$ ID[x]][1 ]if (sels$ week[x] == 2 )<- metadat$ D7.Breath.Count[metadat$ Bird.ID == sels$ ID[x]][1 ]if (sels$ week[x] == 3 )<- metadat$ D14.Bird.Weight..g.[metadat$ Bird.ID == sels$ ID[x]][1 ]if (sels$ week[x] == 4 )<- metadat$ D21.Bird.Weight..g.[metadat$ Bird.ID == sels$ ID[x]][1 ]if (sels$ week[x] == 5 )<- metadat$ D28.Bird.Weight..g.[metadat$ Bird.ID == sels$ ID[x]][1 ]return (out)# remove week 5 <- sels[sels$ week != 5 , ]<- read.csv (file.path (raw_data_path, "foxp2_levels.csv" ))$ Treatment[grep ("HIGH" , foxp2$ Treatment)] <- "High stress" $ Treatment[grep ("Control" , foxp2$ Treatment, ignore.case = T)] <- "Control" $ Treatment[grep ("MED" , foxp2$ Treatment)] <- "Medium stress" $ Treatment <- factor (foxp2$ Treatment,levels = c ("Control" , "Medium stress" , "High stress" ))``` ```{r prepare data for stats, eval = TRUE} <- aggregate (selec ~ ID + week, data = sels, length)# without comparing to week 1 $ call.count <- sapply (1 : nrow (agg_dat), function (x)$ selec[x])$ selec <- NULL $ baseline.CORT <- sapply (1 : nrow (agg_dat), function (x) {<- sels$ cort.baseline[sels$ week == 1 & $ ID == agg_dat$ ID[x]]<- sels$ cort.baseline[sels$ week == agg_dat$ week[x] & $ ID == agg_dat$ ID[x]]if (length (baseline) > 0 & length (current) > 0 )<- mean (current) - mean (baseline) else <- NA return (change)$ stress.response <- sapply (1 : nrow (agg_dat), function (x) {<- sels$ stress.response[sels$ week == 1 & $ ID == agg_dat$ ID[x]]<- sels$ stress.response[sels$ week == agg_dat$ week[x] & $ ID == agg_dat$ ID[x]]if (length (baseline) > 0 & length (current) > 0 )<- mean (current) - mean (baseline) else <- NA return (change)$ stress.CORT <- sapply (1 : nrow (agg_dat), function (x) {<- sels$ cort.stress[sels$ week == 1 & $ ID == agg_dat$ ID[x]]<- sels$ cort.stress[sels$ week == agg_dat$ week[x] & $ ID == agg_dat$ ID[x]]if (length (baseline) > 0 & length (current) > 0 )<- mean (current) - mean (baseline) else <- NA return (change)$ weight <- sapply (1 : nrow (agg_dat), function (x) {<- sels$ weight[sels$ week == 1 & sels$ ID == agg_dat$ ID[x]]<- sels$ weight[sels$ week == agg_dat$ week[x] & $ ID == agg_dat$ ID[x]]if (length (baseline) > 0 & length (current) > 0 )<- mean (current) - mean (baseline) else <- NA return (change)$ breath.rate <- sapply (1 : nrow (agg_dat), function (x) {<- sels$ breath.count[sels$ week == 1 & $ ID == agg_dat$ ID[x]]<- sels$ breath.count[sels$ week == agg_dat$ week[x] & $ ID == agg_dat$ ID[x]]if (length (baseline) > 0 & length (current) > 0 )<- mean (current) - mean (baseline) else <- NA return (change)$ acoustic.diversity <- sapply (1 : nrow (agg_dat), function (x) {<- areas_by_week$ raref.area[areas_by_week$ ID == agg_dat$ ID[x] & $ week == $ week[x]]if (length (area) < 1 )<- NA return (area)$ acoustic.distance <- sapply (1 : nrow (agg_dat), function (x) {<- indiv_ovlp$ distance.to.first.week[indiv_ovlp$ ID == agg_dat$ ID[x] & $ week == agg_dat$ week[x]]if (length (distance) < 1 )<- NA return (distance)$ acustic.plasticity <- sapply (1 : nrow (agg_dat), function (x) {<- indiv_ovlp$ overlap.to.first.week[indiv_ovlp$ ID == agg_dat$ ID[x] & $ week == agg_dat$ week[x]]<- 1 - overlapif (length (plasticity) < 1 )<- NA return (plasticity)$ acoustic.convergence <- sapply (1 : nrow (agg_dat), function (x) {<- group_ovlp$ overlap.to.group[group_ovlp$ ID == agg_dat$ ID[x] & $ week == $ week[x]]if (length (overlap) < 1 )<- NA return (overlap)$ treatment <- sapply (1 : nrow (agg_dat), function (x)unique (sels$ treatment[sels$ ID == $ ID[x]]))$ round <- sapply (1 : nrow (agg_dat), function (x)unique (sels$ round[sels$ ID == $ ID[x]]))$ foxp2 <- sapply (1 : nrow (agg_dat), function (x) {if (agg_dat$ week[x] == 4 )<- foxp2$ FoxP2.Counts.MMSt.Striatum[foxp2$ Bird.ID == agg_dat$ ID[x]][1 ] else <- NA # print(length(fp2)) if (length (fp2) == 0 )<- NA return (fp2)``` # Physiological parameters ## Stats ```{r, eval = FALSE} <- c ("baseline.CORT" ,"stress.response" ,"stress.CORT" ,"weight" ,"breath.rate" )<- c ("~ treatment + week + (1|ID) + (1|round)" )<- expand.grid (responses = responses,predictors = predictors,stringsAsFactors = FALSE <- c (responses, "week" )# remove week 1 <- agg_dat[agg_dat$ week != 1 , ]for (i in vars_to_scale)<- scale (sub_agg_dat[, vars_to_scale])<- lapply (1 : nrow (formulas), function (x) {<- sub_agg_dat[! is.na (sub_agg_dat[names (sub_agg_dat) == formulas$ responses[x]]), ]# replace ~ by "by" in the formula <- paste (formulas$ responses[x],gsub ("~" , "by" , formulas$ predictors[x]))<- brm (formula = paste (formulas$ responses[x], formulas$ predictors[x]),iter = 100000 ,silent = 2 ,family = student (),data = sub_dat,control = list (adapt_delta = 0.9 , max_treedepth = 15 ),chains = 4 ,cores = 4 ,file = paste0 (file.path (processed_data_path,"physio-" ),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)``` ```{r, results = 'asis'} <- list.files (file.path (processed_data_path, "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 ``` ```{r, eval = TRUE, out.width="120%", out.height="200%"} <- list.files (file.path (processed_data_path, "regressions" ),pattern = "physio-" ,full.names = TRUE # read all physio models <- lapply (physio_model_list, readRDS)names (physio_models) <- gsub (".rds" , "" , basename (physio_model_list))<- stack (metadat[, c ("D3.Breath.Count" ,"D7.Breath.Count" ,"D14.Breath.Count" ,"D21.Breath.Count" ,"D28.Breath.Count" <- stack (metadat[, c ("D3.Bird.Weight..g." ,"D7.Bird.Weight..g." ,"D14.Bird.Weight..g." ,"D21.Bird.Weight..g." ,"D28.Bird.Weight..g." <- stack (metadat[, c ("D3.CORT.Stress" ,"D7.CORT.Stress" ,"D14.CORT.Stress" ,"D21.CORT.Stress" ,"D28.CORT.Stress" <- stack (metadat[, c ("D3.CORT.Baseline" ,"D7.CORT.Baseline" ,"D14.CORT.Baseline" ,"D21.CORT.Baseline" ,"D28.CORT.Baseline" <- data.frame (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 - $ values# head(stress) $ week <- factor (sapply (strsplit (as.character (stress$ week), " \\ ." ), "[[" , 1 ),levels = c ("D3" , "D7" , "D14" , "D21" , "D28" ))$ day <- as.numeric (gsub ("D" , "" , as.character (stress$ week)))$ round <- paste ("Round" , stress$ Round)$ Treatment <- gsub ("Medium" , "Med." , stress$ Treatment)$ treatment <- factor (stress$ Treatment,levels = c ("Control" , "Med. Stress" , "High Stress" ))# remove 5th week <- stress[stress$ week != "D28" , ]<- lapply (stress$ Bird.ID, function (x) {<- stress[stress$ Bird.ID == x, ]$ breath.count <- X$ breath.count - X$ breath.count[X$ week == "D3" ]$ weight <- X$ weight - X$ weight[X$ week == "D3" ]$ cort.stress <- X$ cort.stress - X$ cort.stress[X$ week == "D3" ]$ cort.baseline <- X$ cort.baseline - X$ cort.baseline[X$ week == "D3" ]$ stress.response <- X$ stress.response - X$ stress.response[X$ week == "D3" ]return (X)<- do.call (rbind, stress_l)<- aggregate (cbind (breath.count, weight, stress.response, cort.baseline) ~ week + treatment,<- aggregate (cbind (breath.count, weight, stress.response, cort.baseline) ~ week + treatment,names (agg_stress_se) <- paste (names (agg_stress_se), ".se" , sep = "" )<- cbind (agg_stress, agg_stress_se[, 3 : 6 ])$ Week <- 1 : 4 <- if (isTRUE (getOption ('knitr.in.progress' )))10 else 20 <- 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.sewidth = 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 \n rate (breaths/min)" , x = "Week" ) + theme_classic (base_size = bs) + theme (legend.position = "none" )<- 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 \n weight (grams)" , x = "Week" ) + theme_classic (base_size = bs) + theme (legend.position = "none" )<- 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.sewidth = 0.1 + scale_fill_viridis_d (begin = 0.1 , end = 0.9 ) + facet_wrap ( ~ treatment, ncol = 3 , scale = "fixed" ) + labs (y = "Mean change in \n baseline CORT (ng/mL)" , x = "Week" ) + theme_classic (base_size = bs) + theme (legend.position = "none" )<- 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.sewidth = 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 \n response CORT (ng/mL)" , x = "Week" ) + theme_classic (base_size = bs) + theme (legend.position = "none" )<- lapply (physio_models, function (x) {<- grep ("b_" , posterior:: variables (x), value = TRUE )<- posterior:: as_draws_array (x, variable = vars)<- brmsish::: draw_summary (variables = vars,probs = c (0.025 , 0.975 ),robust = TRUE ,spread.type = "HPDI" $ predictor <- rownames (coef_table)$ predictor <- gsub ("b_treatment|b_" , "" , coef_table$ predictor)$ predictor <- gsub ("Stress" , " stress" , coef_table$ predictor)$ predictor <- gsub ("week" , "Week" , coef_table$ predictor)<- coef_table[coef_table$ predictor != "Intercept" , ]<- ggplot2:: ggplot (data = coef_table, aes (x = Estimate, y = predictor)) + geom_vline (xintercept = 0 , lty = 2 ) + :: geom_point (size = 3 , col = col_pointrange) + :: geom_errorbar (:: aes (xmin = ` l-95% CI ` , xmax = ` u-95% CI ` ),width = 0 ,col = col_pointrange+ :: theme_classic (base_size = bs) + :: 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 ()+ :: labs (x = "Effect size" , y = "" )return (gg_coef):: plot_grid (grep ("breath" , names (gg_coeffs_physio))]] + theme_classic (base_size = bs),grep ("weight" , names (gg_coeffs_physio))]] + theme_classic (base_size = bs),grep ("baseline" , names (gg_coeffs_physio))]] + theme_classic (base_size = bs),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' ))) {:: ggsave2 (filename = file.path (output_path, "bar_graphs_and_estimates_physiology_70dpi.jpeg" ),width = 23 ,height = 9 ,device = grDevices:: png:: ggsave2 (filename = file.path (output_path, "bar_graphs_and_estimates_physiology_300dpi.jpeg" ),dpi = 300 ,width = 23 ,height = 9 ,device = grDevices:: png``` < div class = "alert alert-info" > ## Takeaways - Breath rate decreases gradually across time after the first week- Stress response is higher in "high stress" birds compared to first week</ div > # Acoustic space projection ```{r tsne plot 2, eval = TRUE, fig.width= 8, fig.height= 8} <- read.csv (file.path (processed_data_path, "tsne_acoustic_space.csv" ))<- c ("371YYLM" , "395WBHM" )<- matrix (c (11 : 14 , 1 : 4 , 5 : 8 , 22 , 9 , 9 , 18 ),ncol = 4 ,byrow = TRUE )<- cbind (c (23 , 10 , 10 , 20 ), mat, c (21 , 15 , 16 , 19 ))<- rbind (rep (17 , 6 ), mat)if (! isTRUE (getOption ('knitr.in.progress' )))jpeg (file.path (output_path, "acoustic_space_plot.jpg" ),width = 1800 ,height = 2200 ,res = 300 layout (mat,widths = c (0.4 , 1 , 1 , 1 , 1 , 0.2 ),heights = c (4 , 0.2 , 1 , 1 , 0.4 ))par (mar = c (0 , 0 , 0 , 0 ))<- 1 <- viridis (3 , alpha = 0.8 )for (x in ids)for (i in 1 : 4 ) {plot_space (X = tsne_acou_space,dimensions = c ("TSNE1" , "TSNE2" ),indices = which (tsne_acou_space$ ID == x & $ week == i),basecex = 2 ,pch = 20 ,labels = c (x, "group" ),legend = NULL ,xaxt = if (panel_count %in% mat[4 , ])"s" else "n" ,yaxt = "n" ,colors = colors[1 : 2 ],point.cex = 1.7 ,point.alpha = 0.3 ,point.colors = c ("black" , "black" )# add y axis removing the upper tick label if (panel_count %in% mat[, 2 ])axis (2 ,at = seq (- 60 , 60 , 20 ),labels = c (seq (- 60 , 40 , 20 ), "" ),tck = - 0.02 ,lwd = 0.5 ,col = "black" ,col.axis = "black" <- 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 ,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 ($ TSNE1,$ 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 ($ TSNE1,$ TSNE2,col = viridis:: viridis (3 , alpha = 0.3 , direction = 1 )[as.numeric (as.factor (tsne_acou_space$ 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 ```{r} <- space_similarity (treatment ~ TSNE1 + TSNE2, data = tsne_acou_space, method = "density.overlap" )print_kable (ss)``` - Mean overlap: `r mean(ss$mean.overlap)` - Range of overlap: `r min(ss$mean.overlap)` to `r max(ss$mean.overlap)` ## Stats ```{r, eval = FALSE} <- c ("call.count" ,"acoustic.diversity" ,"acustic.plasticity" ,"acoustic.convergence" )<- c ("~ treatment + week + (1|ID) + (1|round)" )<- 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]) <- lapply (1 : nrow (formulas), function (x) {<- 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$ week != 1 , ]<- paste (formulas$ responses[x],gsub ("~" , "by" , formulas$ predictors[x]))<- 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 (processed_data_path, "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 )} )``` ``` {r, results = 'asis' }<- list.files (file.path (processed_data_path, "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 ``` Barplot and effect sizes graph ``` {r, eval = TRUE , out.width= "120%" , out.height= "200%" }<- list.files (file.path (processed_data_path, "regressions" ),pattern = "behav-" ,full.names = TRUE # read all physio models <- lapply (behav_model_list, readRDS)names (behav_models) <- gsub (".rds" , "" , basename (behav_model_list))<- aggregate (cbind (call.count, acoustic.convergence) ~ week + treatment,<- aggregate (cbind (acoustic.diversity, acustic.plasticity) ~ week + treatment,<- aggregate (cbind (call.count, acoustic.convergence) ~ week + treatment,<- aggregate (cbind (acoustic.diversity, acustic.plasticity) ~ week + treatment,<- merge (agg_call.count_se, agg_behav_se, all = TRUE )names (agg_behav_se) <- paste (names (agg_behav_se), ".se" , sep = "" )<- merge (agg_call.count, agg_behav, all = TRUE )<- cbind (agg_behav, agg_behav_se[, 3 : 6 ])<- if (isTRUE (getOption ('knitr.in.progress' )))10 else 20 $ treatment <- gsub ("Medium" , "Med." , agg_behav$ treatment)$ treatment <- factor (agg_behav$ treatment,levels = c ("Control" , "Med. Stress" , "High Stress" ))<- 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.sewidth = 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" )<- 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.sewidth = 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" )<- 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.sewidth = 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" )<- 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.sewidth = 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" )<- lapply (behav_models, function (x) {<- grep ("b_" , posterior:: variables (x), value = TRUE )<- posterior:: as_draws_array (x, variable = vars)<- brmsish::: draw_summary (variables = vars,probs = c (0.025 , 0.975 ),robust = TRUE ,spread.type = "HPDI" $ predictor <- rownames (coef_table)$ predictor <- gsub ("b_treatment|b_" , "" , coef_table$ predictor)$ predictor <- gsub ("Stress" , " stress" , coef_table$ predictor)$ predictor <- gsub ("week" , "Week" , coef_table$ predictor)<- coef_table[coef_table$ predictor != "Intercept" , ]<- ggplot2:: ggplot (data = coef_table, aes (x = Estimate, y = predictor)) + geom_vline (xintercept = 0 , lty = 2 ) + :: geom_point (size = 4 , col = col_pointrange) + :: geom_errorbar (:: aes (xmin = ` l-95% CI ` , xmax = ` u-95% CI ` ),width = 0 ,col = col_pointrange+ :: theme_classic (base_size = bs) + :: 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 ()+ :: labs (x = "Effect size" , y = "" )return (gg_coef):: plot_grid (grep ("count" , names (gg_coeffs_behav))]] + theme_classic (base_size = bs),grep ("diversity" , names (gg_coeffs_behav))]] + theme_classic (base_size = bs),grep ("plasticity" , names (gg_coeffs_behav))]] + theme_classic (base_size = bs),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' ))) {:: ggsave2 (filename = file.path (output_path, "bar_graphs_and_estimates_behavior_70dpi.jpeg" ),width = 23 ,height = 9 ,device = grDevices:: png# :: ggsave2 (filename = file.path (output_path,"bar_graphs_and_estimates_behavior_300dpi.jpeg" ),dpi = 300 ,width = 23 ,height = 9 ,device = grDevices:: png``` <div class="alert alert-info"> ## 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 </div> # FOXP2 ## FoxP2 positive cells in MMST:VSP ### Stats ``` {r, eval = FALSE }<- brm (formula = FoxP2.Counts.MMSt.Striatum ~ Treatment,iter = 100000 ,family = Gamma (link = "log" ),silent = 2 ,data = foxp2,control = list (adapt_delta = 0.9 , max_treedepth = 15 ),chains = 4 ,cores = 4 ,file = file.path (processed_data_path, "regressions/foxp2_by_treatment_model" ),file_refit = "on_change" ,prior = c (prior (normal (0 , 5 ), "b" ), prior (normal (0 , 10 ), "Intercept" ))``` ``` {r, eval = TRUE , results= 'asis' }extended_summary (read.file = file.path (processed_data_path, "regressions/foxp2_by_treatment_model.rds" ),gsub.pattern = "b_treatment|b_" ,gsub.replacement = "" ,highlight = TRUE ,remove.intercepts = TRUE ``` ## FoxP2 mRNA expression MMST:VSP ``` {r, eval = TRUE }<- read.csv (file.path (raw_data_path, "qpcr_foxp2_expression.csv" ))``` ### Stats ``` {r, eval = FALSE }<- brm (formula = FoxP2.mRNA.MMST.VSP ~ Treatment,iter = 100000 ,family = Gamma (link = "log" ),silent = 2 ,data = qpcr,control = list (adapt_delta = 0.9 , max_treedepth = 15 ),chains = 4 ,cores = 4 ,file = file.path (processed_data_path, "regressions/qpcr_by_treatment_model" ),file_refit = "on_change" ,prior = c (prior (normal (0 , 5 ), "b" ), prior (normal (0 , 10 ), "Intercept" ))``` ``` {r, eval = TRUE , results= 'asis' }extended_summary (read.file = file.path (processed_data_path, "regressions/qpcr_by_treatment_model.rds" ),gsub.pattern = "b_treatment|b_" ,gsub.replacement = "" ,highlight = TRUE ,remove.intercepts = TRUE ``` Barplot and effect sizes graph ``` {r, eval = TRUE , out.width= "120%" , out.height= "200%" }<- aggregate (FoxP2.Counts.MMSt.Striatum ~ Treatment, foxp2, mean)$ se <- aggregate (FoxP2.Counts.MMSt.Striatum ~ Treatment, foxp2, standard_error)$ FoxP2.Counts.MMSt.Striatum<- ggplot (data = agg_foxp2,aes (x = Treatment, y = FoxP2.Counts.MMSt.Striatum, fill = Treatment)) + geom_bar (stat = "identity" ) + geom_errorbar (aes (ymin = FoxP2.Counts.MMSt.Striatum - se, ymax = FoxP2.Counts.MMSt.Striatum + se),width = 0.1 + scale_fill_viridis_d (begin = 0.1 , end = 0.9 ) + scale_x_discrete (labels = c ("Control" = "Control" ,"Medium stress" = "Medium \n stress" ,"High stress" = "High \n stress" + scale_fill_viridis_d (begin = 0.1 , end = 0.9 ) + labs (y = "FoxP2 positive \n cells in MMST:VSP" , x = "Treatment" ) + theme_classic (base_size = bs) + theme (legend.position = "none" )<- aggregate (FoxP2.mRNA.MMST.VSP ~ Treatment, qpcr, mean)$ se <- aggregate (FoxP2.mRNA.MMST.VSP ~ Treatment, qpcr, standard_error)$ FoxP2.mRNA.MMST.VSP# order levels $ Treatment <- factor (agg_foxp2$ Treatment,levels = c ("Control" , "Medium stress" , "High stress" ))$ Treatment <- factor (agg_mrna$ Treatment,levels = c ("Control" , "Medium" , "High Stress" ))<- ggplot (data = agg_mrna,aes (x = Treatment, y = FoxP2.mRNA.MMST.VSP, fill = Treatment)) + geom_bar (stat = "identity" ) + geom_errorbar (aes (ymin = FoxP2.mRNA.MMST.VSP - se, ymax = FoxP2.mRNA.MMST.VSP + se),width = 0.1 ) + scale_x_discrete (labels = c ("Control" = "Control" ,"Medium" = "Medium \n stress" ,"High Stress" = "High \n stress" + scale_fill_viridis_d (begin = 0.1 , end = 0.9 ) + labs (y = "FoxP2 mRNA \n expression MMST:VSP" , x = "Treatment" ) + theme_classic (base_size = bs) + theme (legend.position = "none" )<- readRDS ( file.path (processed_data_path, "regressions/foxp2_by_treatment_model.rds" ))<- readRDS (file.path (processed_data_path, "regressions/qpcr_by_treatment_model.rds" ))<- lapply (list (foxp2_cells_mod, qpcr_mod), function (x) {<- grep ("b_" , posterior:: variables (x), value = TRUE )<- posterior:: as_draws_array (x, variable = vars)<- brmsish::: draw_summary (variables = vars,probs = c (0.025 , 0.975 ),robust = TRUE ,spread.type = "HPDI" $ predictor <- rownames (coef_table)$ predictor <- gsub ("b_treatment|b_" , "" , coef_table$ predictor)$ predictor <- gsub ("Stress" , " stress" , coef_table$ predictor)$ predictor <- gsub ("week" , "Week" , coef_table$ predictor)$ predictor <- gsub ("Treatment" , "" , coef_table$ predictor)$ predictor <- gsub ("Mediumstress|Medium" ,"Medium stress" ,$ predictor)$ predictor <- gsub ("Highstress" , "High stress" , coef_table$ predictor)<- coef_table[coef_table$ predictor != "Intercept" , ]<- ggplot2:: ggplot (data = coef_table, aes (x = Estimate, y = predictor)) + geom_vline (xintercept = 0 , lty = 2 ) + :: geom_point (size = 4 , col = col_pointrange) + :: geom_errorbar (:: aes (xmin = ` l-95% CI ` , xmax = ` u-95% CI ` ),width = 0 ,col = col_pointrange+ :: theme_classic (base_size = bs) + :: 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 ()+ :: labs (x = "Effect size" , y = "" )return (gg_coef):: plot_grid (1 ]] + theme_classic (base_size = bs),2 ]] + theme_classic (base_size = bs),nrow = 2 ,rel_heights = c (1.8 , 1 )# try bs = 20 for saving plots :: ggsave2 (filename = file.path (output_path,"bar_graphs_and_estimates_foxp2_70dpi.jpeg" ),width = 15 ,height = 10 ,device = grDevices:: png:: ggsave2 (filename = file.path (output_path,"bar_graphs_and_estimates_foxp2_300dpi.jpeg" ),dpi = 300 ,width = 15 ,height = 10 ,device = grDevices:: png``` # Combined model diagnostics ``` {r, eval = TRUE , results= 'asis' }check_rds_fits (path = file.path (processed_data_path, "regressions" ), html = TRUE )``` --- # Session information {.unnumbered .unlisted} ``` {r session info, echo= F, eval = TRUE }:: session_info ()``` 
