--- title: Data analysis subtitle: Age, vocal learning and FoxP2 in budgerigards author: Marcelo Araya-Salas 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 editor_options: chunk_output_type: console --- <!-- skyblue box --> <div class = "alert alert-info" > </div> <!-- this code add line numbers to code blocks --> <!-- only works when code folding is not used in yaml (code_folding: show) --> ```{=html} <style> body { counter-reset: source-line 0; } pre.numberSource code { counter-reset: none; } </style> ``` ```{r set root directory, echo = FALSE} :: opts_knit$ set (root.dir = ".." )``` ```{r add link to github repo, 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 Source code and data 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 ,warning = FALSE ``` <!-- skyblue box --> <div class = "alert alert-info" > # Purpose - Testing for correlation between vocal measures and FoxP2 MMST/VSP Protein Expression and Expression Difference with Age</div> # Load packages ```{r load packages} # knitr is require for creating html/pdf/word reports # load function from sketchy source ("https://raw.githubusercontent.com/maRce10/sketchy/main/R/load_packages.R" )# install/ load packages load_packages (packages = c ("ggplot2" , "tidyverse" , "cowplot" , "pwr" , "brms" , "viridis" ))<- function (...) viridis (alpha = 0.5 , begin = 0.3 , end = 0.7 ,source ("https://raw.githubusercontent.com/maRce10/brmsish/master/R/extended_summary.R" )source ("https://raw.githubusercontent.com/maRce10/brmsish/master/R/helpers.R" )source ("https://raw.githubusercontent.com/maRce10/brmsish/master/R/check_rds_fits.R" )``` # Load data ```{r} <- read.csv ("./data/raw/Vocal_Brain_Corr_03_09_2023_Final.csv" , header = TRUE )head (voc_brain_corr)``` - vocal.plasticity = (1-acoustic.overlap) ```{r} $ vocal.plasticity <- (1 - (voc_brain_corr$ acoustic.overlap))``` # Exploratory graphs ```{r, fig.height=4} $ treatment <- factor (voc_brain_corr$ treatment, levels = c ("Young" , "Old" ))#Visualize data theme_set (theme_classic ())<- ggplot (voc_brain_corr,aes (x = mmst.vsp,y = acoustic.area, colour = treatment)) + geom_point (size = 6 ) + geom_smooth (method = "lm" , se = FALSE , size = 1.5 ) + labs (y = "Vocal diversity" , x = "FoxP2 MMSt/VSP expression" , color = "Adult age" ) + theme (text = element_text (size = 20 ),axis.line = element_line (colour = 'black' , size = 1.5 ),axis.ticks = element_line (colour = "black" , size = 1.5 ),axis.text = element_text (colour = "black" )+ scale_color_viridis_d (begin = 0.2 , end = 0.8 , alpha = 0.5 ) + coord_cartesian (xlim = c (0.25 , 1.25 ))ggsave ("./output/FOXP2_VocalDiversity_Correlation_Final_03_09_2023.png" ,width = 7 ,height = 6 ,units = "in" <- ggplot (voc_brain_corr,aes (x = mmst.vsp, y = vocal.plasticity, colour = treatment)) + geom_point (size = 6 ) + geom_smooth (method = "lm" , se = FALSE , size = 1.5 ) + labs (y = "Vocal plasticity" , x = "FoxP2 MMSt/VSP expression" , color = "Adult age" ) + theme (text = element_text (size = 20 ),axis.line = element_line (colour = 'black' , size = 1.5 ),axis.ticks = element_line (colour = "black" , size = 1.5 ),axis.text = element_text (colour = "black" )+ scale_color_viridis_d (begin = 0.2 , end = 0.8 , alpha = 0.5 ) + coord_cartesian (xlim = c (0.25 , 1.25 ), ylim = c (0 , 1 ))ggsave ("./output/FOXP2_VocalPlasticity_Correlation_Final__03_09_2023.png" ,width = 7 ,height = 6 ,units = "in" <- ggplot (voc_brain_corr,aes (x = mmst.vsp,y = acoustic.overlap.to.group,colour = treatment)) + geom_point (size = 6 ) + geom_smooth (method = "lm" , se = FALSE , size = 1.5 ) + labs (y = "Vocal convergence" , x = "FoxP2 MMSt/VSP expression" , color = "Adult age" ) + theme (text = element_text (size = 20 ),axis.line = element_line (colour = 'black' , size = 1.5 ),axis.ticks = element_line (colour = "black" , size = 1.5 ),axis.text = element_text (colour = "black" )+ scale_color_viridis_d (begin = 0.2 , end = 0.8 , alpha = 0.5 ) + coord_cartesian (xlim = c (0.25 , 1.25 ), ylim = c (0 , 1 ))ggsave ("./output/FOXP2_VocalConvergence_Correlation_Final__03_03_2023.png" ,width = 7 ,height = 6 ,units = "in" ``` # Statistical analysis ## Models without age as a predictor ### Fit models ```{r, eval = FALSE} <- 20000 <- 4 <- c (prior (normal (0 , 6 ), class = "b" ))# mean center expression $ mmst.vsp_cntr <- voc_brain_corr$ mmst.vsp - mean (voc_brain_corr$ mmst.vsp)# cmdstanr::set_cmdstan_path('~/Documentos/cmdstan/') # to run within-chain parallelization .01 <- brm (formula = acoustic.area ~ mmst.vsp_cntr, data = voc_brain_corr,prior = priors, iter = iter, chains = chains, cores = chains,control = list (adapt_delta = 0.99 , max_treedepth = 15 ), file = "./data/processed/regression_model_acoustic_area_by_expression.RDS" ,file_refit = "always" ).02 <- brm (formula = vocal.plasticity ~ mmst.vsp_cntr, data = voc_brain_corr,prior = priors, iter = iter, chains = chains, cores = chains, family = Beta (),control = list (adapt_delta = 0.99 , max_treedepth = 15 ), file = "./data/processed/regression_model_vocal_plasticity_by_expression.RDS" ,file_refit = "always" ).03 <- brm (formula = acoustic.overlap.to.group ~ mmst.vsp_cntr, data = voc_brain_corr,prior = priors, iter = iter, chains = chains, cores = chains, family = Beta (),control = list (adapt_delta = 0.99 , max_treedepth = 15 ), file = "./data/processed/regression_model_vocal_overlap_to_group_by_expression.RDS" ,file_refit = "always" )``` ### Change in acoustic area (diversity) ```{r, results='asis', fig.height=3} <- readRDS ("./data/processed/regression_model_acoustic_area_by_expression.RDS" )extended_summary (mod, n.posterior = 1000 , fill = viridis (10 )[7 ], trace.palette = my.viridis,remove.intercepts = TRUE , highlight = TRUE , print.name = FALSE )``` ### Overlap of acoustic area to self (plasticity) ```{r, results='asis', fig.height=3} <- readRDS ("./data/processed/regression_model_vocal_plasticity_by_expression.RDS" )extended_summary (mod, n.posterior = 1000 , fill = viridis (10 )[7 ], trace.palette = my.viridis,remove.intercepts = TRUE , highlight = TRUE , print.name = FALSE )``` ### Overlap to group acoustic area (convergence) ```{r, results='asis', fig.height=3} <- readRDS ("./data/processed/regression_model_vocal_overlap_to_group_by_expression.RDS" )extended_summary (mod, n.posterior = 1000 , fill = viridis (10 )[7 ], trace.palette = my.viridis,remove.intercepts = TRUE , highlight = TRUE , print.name = FALSE )``` ## Models with foxp2 expression & age as predictors ### Fit models ```{r, eval = FALSE} # to run within-chain parallelization .01 <- brm (formula = acoustic.area ~ mmst.vsp_cntr * treatment, data = voc_brain_corr,prior = priors, iter = iter, chains = chains, cores = chains,control = list (adapt_delta = 0.99 , max_treedepth = 15 ), file = "./data/processed/regression_model_acoustic_area_by_expression_and_treatment.RDS" ,file_refit = "always" ).02 <- brm (formula = vocal.plasticity ~ mmst.vsp_cntr * treatment, data = voc_brain_corr,prior = priors, iter = iter, chains = chains, cores = chains, family = Beta (),control = list (adapt_delta = 0.99 , max_treedepth = 15 ), file = "./data/processed/regression_model_vocal_plasticity_by_expression_and_treatment.RDS" ,file_refit = "always" ).03 <- brm (formula = acoustic.overlap.to.group ~ mmst.vsp_cntr * treatment, data = voc_brain_corr,prior = priors, iter = iter, chains = chains, cores = chains, family = Beta (),control = list (adapt_delta = 0.99 , max_treedepth = 15 ), file = "./data/processed/regression_model_vocal_overlap_to_group_by_expression_and_treatment.RDS" ,file_refit = "always" )``` ### Change in acoustic area (diversity) ```{r, results='asis', fig.height=3} <- readRDS ("./data/processed/regression_model_acoustic_area_by_expression_and_treatment.RDS" )extended_summary (mod, n.posterior = 1000 , fill = viridis (10 )[7 ], trace.palette = my.viridis,remove.intercepts = TRUE , highlight = TRUE , print.name = FALSE )``` ### Overlap of acoustic area to self (plasticity) ```{r, results='asis', fig.height=3} <- readRDS ("./data/processed/regression_model_vocal_plasticity_by_expression_and_treatment.RDS" )extended_summary (mod, n.posterior = 1000 , fill = viridis (10 )[7 ], trace.palette = my.viridis,remove.intercepts = TRUE , highlight = TRUE , print.name = FALSE )``` ### Overlap to group acoustic area (convergence) ```{r, results='asis', fig.height=3} <- readRDS ("./data/processed/regression_model_vocal_overlap_to_group_by_expression_and_treatment.RDS" )extended_summary (mod, n.posterior = 1000 , fill = viridis (10 )[7 ], trace.palette = my.viridis,remove.intercepts = TRUE , highlight = TRUE , print.name = FALSE )``` ## Combined model metadata ```{r, results='asis'} check_rds_fits (path = "./data/processed/" , html = TRUE )``` <div class = "alert alert-success" > # Takeaways {.unnumbered .unlisted} - Lower diversity in old birds</div> <!-- '---' adds a gray vertical line --> ## Posterior predictive checks ```{r, results='asis'} <- c (divergence = "./data/processed/regression_model_acoustic_area_by_expression.RDS" ,plasticity = "./data/processed/regression_model_vocal_plasticity_by_expression.RDS" ,convergence = "./data/processed/regression_model_vocal_overlap_to_group_by_expression.RDS" ,divergence_interaction = "./data/processed/regression_model_acoustic_area_by_expression_and_treatment.RDS" ,plasticity_interaction = "./data/processed/regression_model_vocal_plasticity_by_expression_and_treatment.RDS" ,convergence_interaction = "./data/processed/regression_model_vocal_overlap_to_group_by_expression_and_treatment.RDS" )<- 20 for (i in seq_len (length (model_list))) {cat (" \n\n ## " , names (model_list)[i], " \n\n " )<- readRDS (model_list[[i]])<- pp_check (fit, ndraws = ndraws, type = "dens_overlay" ) # shows dens_overlay plot by default <- pp_check (fit, type = "stat" , stat = "mean" , ndraws = ndraws)<- pp_check (fit, type = "scatter_avg" , ndraws = ndraws)<- pp_check (fit, type = "stat_2d" , ndraws = ndraws)<- pp_check (fit, type = "loo_pit_qq" , ndraws = ndraws)<- pp_check (fit, type = "error_scatter_avg" , ndraws = ndraws)<- list (ppc_dens, pp_mean, pp_scat, pp_error, pp_stat2,<- lapply (plot_l, function (x) x + scale_color_viridis_d (begin = 0.1 ,end = 0.8 , alpha = 0.5 ) + scale_fill_viridis_d (begin = 0.1 ,end = 0.8 , alpha = 0.5 ) + theme_classic ())print (plot_grid (plotlist = plot_l, ncol = 2 ))``` <!-- add packages used, system details and versions --> <font size = "4" > Session information</font> ```{r session info, echo=F} sessionInfo ()``` 
