update readme, small edits re-running code to check

Audiograms linked to anatomy.R

@@ -4,6 +4,7 @@ library(ggpubr)
 library(flextable)
 library(officer)
 library(dplyr)
+library(PerformanceAnalytics)
 
 ####create averaged values for instances where multiple species match a congener with audiogram####
 # add avg Corvus and Phalacrocorax values----------------------------------------
@@ -23,7 +24,6 @@ cong_avg<-cong_avg[-c(grep('Corvus_albus|Corvus_splendens', cong_avg$Binomial)),
 cong_avg<-cong_avg[-c(grep('Phalacrocorax_capensis|Phalacrocorax_lucidus|Phalacrocorax_neglectus', cong_avg$Binomial)), ]
 avgdf<-cong_avg
 
-avgdf$aud_spp<-distinctdforder$spp_audio
 
 # load audiograms ---------------------------------------------------------
 fig1<-read.csv(paste0(getwd(),"/audiograms.csv"), stringsAsFactors = FALSE)
@@ -154,7 +154,7 @@ limits$spp_aud[limits$binomial=="Phalacrocorax_carbo"]<-"Phalacrocorax_carbo"
 aud_data<- limits[,c("LowHzlimit","HighHzlimit","besthz","bestsensitivity")]
 audlog<-aud_data %>% mutate_at(vars(c("LowHzlimit","HighHzlimit","besthz")),log)
 
-library(PerformanceAnalytics)
+
 chart.Correlation(audlog, histogram = TRUE, method = "pearson")
 
 # p-values from correlation tests
@@ -257,7 +257,7 @@ categorylist_lf<-c("Stiffness",
                    "Columella size",
                    "Columella size")
 
-categorylist_bs<-categorylist_lf
+categorylist_bs<-categorylist_lf#category list is the same, regardless of which metric we are predicting
 categorylist_bh<-categorylist_lf
 categorylist_hf<-categorylist_lf
 

README.md

@@ -1,51 +1,50 @@
 # Workflow
-This work corresponds to in prep for Hearing Research  "Inter-specific scaling of ear morphology of birds and its implications for hearing performance".
+This work documents the statistics used for a manuscript in prep for submission to Hearing Research  "Scaling of ear morphology across 127 bird species and its implications for hearing performance".
 
-**To reproduce the the scaling pgls analyses:**
-
-Start with "Set up data_scl.R". This loads the data and phylogeny. This also runs the phylogenetic regressions (phylogenetic generalized least squares regressions, or PGLS) between ear measures and between each ear measure and head mass. This script also exports the summary statistics and tables to csv or word files.
-
-**To run the pgls between ear structures and audiometric measurements:**  
-
-First,  "Set up data_scl.R" also has to be run to get the data prepared (minus the pgls from those scripts). Next, the scripts in "Audiograms linked to anatomy.R" can be run. This creates a 'limits' dataframe with the audiogram metrics, which will be joined with the anatomy dataframes for analysis. Once the limits df is created the "pgls_resids re headmass.R" can be used to run the pgls analyses with head mass-corrected values.
-
-**Scripts for plots are in the plots folder.**  
-For scatterplots, the associated analysis files must be run before plotting.  
+## **To reproduce the the scaling pgls analyses:**
 
+Start with _**"Set up data_scl.R".**_ This loads the data and phylogeny. This also runs the phylogenetic regressions (phylogenetic generalized least squares regressions, or PGLS) between ear measures and between each ear measure and head mass. This script also exports the summary statistics and tables to csv or word files.
 
+## **To run the pgls between ear structures and audiometric measurements:**  
 
+First,  _**"Set up data_scl.R"**_ also has to be run to get the data prepared (minus the pgls from those scripts). Next, the scripts in _**"Audiograms linked to anatomy.R"**_ can be run. This creates a 'limits' dataframe with the audiogram metrics, which will be joined with the anatomy dataframes for analysis. Once the limits df is created the _**"pgls_resids re headmass.R"**_ can be used to run the pgls analyses with head mass-corrected values.  
 
+<br>
 
 
 ## Data files
 |File|Description|
 |-----|-----|
-|["databmadded.csv"](https://github.com/jzeyl/Scaling_2021/blob/main/databmadded.csv)|anatomical data|
-|["Column name descriptions.csv"](https://github.com/jzeyl/Scaling_2021/blob/main/Column%20name%20descriptions.csv)|a file describing the column names in data files|
-|["audiograms.csv"](https://github.com/jzeyl/Scaling_2021/blob/main/audiograms.csv)|threshold data from behavioural audiograms collected from literature|
-|["JZ tree Prum merged hackett.tree"](https://github.com/jzeyl/Scaling_2021/blob/main/JZ%20tree%20Prum%20merged%20hackett.tree)|phylogeny file|
+|["databmadded.csv"](https://github.com/jzeyl/Scaling_2021/blob/main/databmadded.csv)|Anatomical data|
+|["Column name descriptions.csv"](https://github.com/jzeyl/Scaling_2021/blob/main/Column%20name%20descriptions.csv)|File describing the column names in data files|
+|["audiograms.csv"](https://github.com/jzeyl/Scaling_2021/blob/main/audiograms.csv)|Audiogram threshold data collected from literature|
+|["JZ tree Prum merged hackett.tree"](https://github.com/jzeyl/Scaling_2021/blob/main/JZ%20tree%20Prum%20merged%20hackett.tree)|Phylogeny file|
 
+<br>
 
 ## Scaling analyses
 
 |File|Description|
 |-----|-----|
-|["Set up data_scl.R"](https://github.com/jzeyl/Scaling_2021/blob/main/Set%20up%20data_scl.R)|This is the main analysis script, which imports the data analysis and runs pgls regressions, calling formulas make in other R scripts. Outputs are tables with statistical results (hm, intra, bm) exported to csv or word.|
-|["pgls_intraear.R"](https://github.com/jzeyl/Scaling_2022/blob/main/pgls_intraear.R)|pgls models run between auditory structures|
-|["pgls_bm.R"](https://github.com/jzeyl/Scaling_2022/blob/main/pgls_bm.R)|pgls models run against body mass|
-|["pgls_HM.R"](https://github.com/jzeyl/Scaling_2022/blob/main/pgls_HM.R)|pgls models run against head mass|
+|["Set up data_scl.R"](https://github.com/jzeyl/Scaling_2021/blob/main/Set%20up%20data_scl.R)|This is the main analysis script, which imports the data analysis and runs pgls regressions, calling formulas from other R scripts (below). Outputs are tables with statistical results exported to csv or word.|
+|["pgls_intraear.R"](https://github.com/jzeyl/Scaling_2022/blob/main/pgls_intraear.R)|Runs pgls regressions between auditory structures|
+|["pgls_HM.R"](https://github.com/jzeyl/Scaling_2022/blob/main/pgls_HM.R)|Runs pgls regressions pgls models against head mass|
+|["pgls_bm.R"](https://github.com/jzeyl/Scaling_2022/blob/main/pgls_bm.R)|Runs pgls regressions pgls models against body mass|  
+<br>
 
 ## Audiometry analyses
 
 |File|Description|
 |-----|-----|
-|["Audiograms linked to anatomy.R"](https://github.com/jzeyl/Scaling_2022/blob/main/Audiograms%20linked%20to%20anatomy.R)|This is the main analysis script for analyses between anatomy and audiogram metrics. This script runs the pgls regressions, calling formulas make in other R scripts. Outputs are tables with statistical results exported to csv or word.|
-|["pgls_audiogram_bs.R"](https://github.com/jzeyl/Scaling_2021/blob/main/pgls_audiogram_bs.R)|pgls models run against best sensitivity|
-|["pgls_audiogram_hf.R"](https://github.com/jzeyl/Scaling_2021/blob/main/pgls_audiogram_hf.R)|pgls models run against high frequency limit|
-|["pgls_audiogram_lf.R"](https://github.com/jzeyl/Scaling_2021/blob/main/pgls_audiogram_lf.R)|pgls models run against low frequency limit|
-|["pgls_audiogram_bh.R"](https://github.com/jzeyl/Scaling_2021/blob/main/pgls_audiogram_hf.R)|pgls models run against best frequency| 
-
-
+|["Audiograms linked to anatomy.R"](https://github.com/jzeyl/Scaling_2022/blob/main/Audiograms%20linked%20to%20anatomy.R)|This is the main analysis script for pgls analyses between anatomy and audiogram metrics, which calls other scripts (below). Outputs are tables with statistical results exported to csv or word.|
+|["pgls_audiogram_bs.R"](https://github.com/jzeyl/Scaling_2021/blob/main/pgls_audiogram_bs.R)|Runs pgls models to predict best sensitivity|
+|["pgls_audiogram_hf.R"](https://github.com/jzeyl/Scaling_2021/blob/main/pgls_audiogram_hf.R)|Runs pgls models to predict high frequency limit|
+|["pgls_audiogram_lf.R"](https://github.com/jzeyl/Scaling_2021/blob/main/pgls_audiogram_lf.R)|Runs pgls models to predict low frequency limit|
+|["pgls_audiogram_bh.R"](https://github.com/jzeyl/Scaling_2021/blob/main/pgls_audiogram_hf.R)|Runs pgls models to predict best frequency| 
+|["pgls_resids re headmass.R"](https://github.com/jzeyl/Scaling_2021/blob/main/pgls_resids%20re%20headmass.R)|Runs the pgls analyses with head mass-corrected values|
+
+# Plots
+R scripts for reproducing plots are in the 'plots' folder. The associated analysis files (if applicable) must be run before the associated plotting scripts.    
 
 
 

Set up data_scl.R

@@ -27,8 +27,9 @@ source("load phylogeny and make CDO.R")
 
 #Some missing headmass values to be imputed using PGLS of skull width and head mass
 #Computed head mass from head mass~skullwidth pgls
+
 df$HM#without imputed values
-source("SW_HM_.R")#add phylogeny here
+source("SW_HM_.R")#
 df$HM#with imputed values
 
 
@@ -63,10 +64,11 @@ birdCDO<-comparative.data(phy = birdtreels,data = avgdf,#[avgdf$Category!="Terre
 #check any tips dropped between linking phylogeny and dataframe
 birdCDO$dropped
 
-######If doing audiogram analyses now proceed to
+######If doing audiogram analyses, you can now proceed to
 ######'Audiograms linked to anatomy.R'
 ######Otherwise, proceed to do scaling between structures and with head mass
-
+###To do scaling analyses corrected for head mass, proceed
+###to 'pgls_resids re headmass.R'
 
 #########scaling intraear##########
 #set up intra-ear analyses
@@ -105,6 +107,7 @@ categorylist_intra<-c(rep("Impedance match",5),
 #run the pgls models
 source("pgls_intraear.R")
 
+#Summarizing model output
 #remove intercept estimates, drop model column,
 #only keep significant relationships
 
@@ -129,6 +132,7 @@ options(scipen = 100, digits = 2)
 intra<-intra %>% select(category, ymodel_nolog,Coefficients,
                   geometric_exp, pglsslope,scalingtype,Adj_Rsquared,pval, Lambda) %>%
   filter(Coefficients!="(Intercept)")
+
 # remove the "log" from 'Coefficients'
 #intra$xmodel_nolog<-numeric()
 for(i in seq_along(intra$Coefficients)){
@@ -137,14 +141,14 @@ for(i in seq_along(intra$Coefficients)){
 
 #sort table by category and then adjusted R2
 intra$category<-as.factor(intra$category)
-levels(intra$category)<-c("Impedance match",
+
 intra<-arrange(intra,factor(intra$category, levels = c("Impedance match", "Stiffness", "Columella morphology")),
                desc(Adj_Rsquared))
 intra$pval<-format(round(intra$pval, 3), nsmall = 3)
 
 
 
-#visualize the table better using the flextable package
+#visualize the table using the flextable package
 flexall<-flextable(intra) %>%
   add_header_lines(values = "Table X. ") %>%
   #bold(i = ~ P.val < 0.05) %>% # select columns add: j = ~ Coefficients + P.val
@@ -232,6 +236,7 @@ options(scipen = 100, digits = 2)
 hm<-hm %>% select(category, ymodel_nolog,Coefficients,
                   geometric_exp, pglsslope,scalingtype,Adj_Rsquared,pval, Lambda) %>%
   filter(Coefficients!="(Intercept)")
+
 # remove the "log" from 'Coefficients'
 #hm$xmodel_nolog<-numeric()
 for(i in seq_along(hm$Coefficients)){
@@ -245,8 +250,8 @@ hm<-arrange(hm,factor(hm$category, levels = c(
   "Auditory endorgan length",
   "Input/output areas",
   "Stiffness",
-  "Impedance match")),desc(Adj_Rsquared)
-            )
+  "Impedance match")),desc(Adj_Rsquared))
+
 hm$pval<-format(round(hm$pval, 3), nsmall = 3)
 
 
@@ -264,7 +269,7 @@ par(mar=c(1,1,1,1))
 plots_hm<-lapply(pgls_models_list, plot)
 plots_hm
 
-#write table to word file
+
 #write table to word file
 toprint<-read_docx() #create word doc object
 body_add_flextable(toprint,flexall)#add pgls output table
@@ -273,8 +278,7 @@ body_end_section_landscape(toprint)
 print(toprint,target = paste0(choose.dir(),"/pgls_hm_scaling all_Apr4 2022.docx"))
 
 
-
-#body mass independent
+#body mass vs head mass
 bm_vs_hm<-pgls_models(log(HM)~log(BM_lit))
 summary(bm_vs_hm)
 

pgls_resids re headmass.R

@@ -1,19 +1,21 @@
 library(patchwork)
 library(ggrepel)
 
+###note the dataframe 'limits', created from 'Audiograms linked to anatomy.R'
+###is required to run this code
+
+birdCDO<-comparative.data(phy = birdtreels,data = avgdf,#[avgdf$Category!="Terrestrial",]
+                          names.col = Binomial,
+                          vcv = TRUE, na.omit = FALSE,
+                          warn.dropped = TRUE)
+
 
 #ensure phylogeny matches dataframe
 #remake comparative data frame object with averaged congeners
 #rename phylogeny tips to matching with the species for which audiogram is available
 birdtreels$tip.label[14]<-"Corvus_cornix" #renamed from Corvus_albus
 birdtreels$tip.label[51]<-"Phalacrocorax_carbo" #rename "phalacrocorax_lucidus"
 
-
-birdCDO<-comparative.data(phy = birdtreels,data = avgdf,#[avgdf$Category!="Terrestrial",]
-                          names.col = Binomial,
-                          vcv = TRUE, na.omit = FALSE,
-                          warn.dropped = TRUE)
-
 #check any tips dropped between linking phylogeny and dataframe
 birdCDO$dropped
 
@@ -89,8 +91,8 @@ for(i in seq_along(pgls_todo_hm)){
 #}
 
 #join residual data into single dataframe
-joined<-limits %>% full_join(.,resids_df_list[[1]],by = c("spp_aud" = "resid_bname"))%>% 
-  full_join(.,resids_df_list[[2]],by = c("spp_aud" = "resid_bname"))%>% 
+joined<-limits %>% full_join(.,resids_df_list[[1]],by = c("spp_aud" = "resid_bname"))%>%
+  full_join(.,resids_df_list[[2]],by = c("spp_aud" = "resid_bname"))%>%
   full_join(.,resids_df_list[[3]],by = c("spp_aud" = "resid_bname"))%>%
   full_join(.,resids_df_list[[4]],by = c("spp_aud" = "resid_bname"))%>%
   full_join(.,resids_df_list[[5]],by = c("spp_aud" = "resid_bname"))%>%
@@ -103,16 +105,17 @@ joined<-limits %>% full_join(.,resids_df_list[[1]],by = c("spp_aud" = "resid_bna
   full_join(.,resids_df_list[[12]],by = c("spp_aud" = "resid_bname"))
 
 #only keep audiogram species
-joined<-joined[which(!is.na(joined$aud_rel)),] 
+joined<-joined[which(!is.na(joined$aud_rel)),]
 
 names(joined)
 
 #remove tilda from names to not mess up pgls formulas based on names
-joined<-joined %>% rename_with(~ gsub("~", "vs", .x, fixed = TRUE)) 
+joined<-joined %>% rename_with(~ gsub("~", "vs", .x, fixed = TRUE))
 residlist<-gsub("[()]","_",names(joined)[25:36])
 joined<-setNames(joined,c(names(joined)[1:24],residlist))
 
 names(joined)
+
 birdCDO<-comparative.data(phy = birdtreels,data = joined,#[avgdf$Category!="Terrestrial",]
                           names.col =binomial,
                           vcv = TRUE, na.omit = F,