diff --git a/.RData b/.RData
index deef0ad..04b96c1 100644
Binary files a/.RData and b/.RData differ
diff --git a/.Rhistory b/.Rhistory
index e693854..ae9e94f 100644
--- a/.Rhistory
+++ b/.Rhistory
@@ -1,512 +1,512 @@
-xlab("Frequency(Hz)")+
-#theme(legend.position = "none")
-#ylim(c(-2,20))+
-geom_boxplot(data = limits,aes(x = LowHzlimit, y = 1), width = 0.8)+
-geom_boxplot(data = limits,aes(x =  besthz, y = 2), width = 0.8)+
-geom_boxplot(data = limits,aes(x = HighHzlimit, y = 3), width = 0.8)+
-geom_hline(yintercept = 3.5)+
-theme(axis.title.y = element_text(angle= 0, vjust = 0.5, hjust=1),
-axis.text.x = element_text(angle= 0, vjust = -2.5, hjust=0.5))
-bound2<-bound
-levels(bound2$Species)<-c(levels(bound2$Species),"High freq. limit","Best freq.","Low freq. limit")
-levels(bound2$Species)
-levels(bound2$Species)<-c(levels(bound2$Species),"High freq. limit","Best freq.","Low freq. limit")
-levels(bound2$Species)
-levels(bound$Species)
-bound$Hz<-bound$x
-bound$Species = with(bound, reorder(Species, besthz, median))
-bound2$Species = with(bound2, reorder(Species, besthz, median))
-bound2$`Threshold (dB)`<-bound2$y
-bound$Species<-as.factor(bound$Species)
-levels(bound$Species)
-bound2<-bound
-levels(bound2$Species)<-c(levels(bound2$Species),"High freq. limit","Best freq.","Low freq. limit")
-bound2[95001,c("Species","besthz")]<-c("Low freq. limit",10)
-bound2[95002,c("Species","besthz")]<-c("Best freq.", 10)
-bound2[95003,c("Species","besthz")]<-c("High freq. limit" ,10)
-bound$Hz<-bound$x
-bound$Species = with(bound, reorder(Species, besthz, median))
-bound2$Species = with(bound2, reorder(Species, besthz, median))
-bound2$`Threshold (dB)`<-bound2$y
-warnings
-warnings()
-ggplot(bound, aes_string(x = "x", y = "Species"))+   #, col= supraorder
-geom_path(data = bound2[c(seq(1, 95000, 10),95001,95002,95003), ],aes(col = `Threshold (dB)`), size = 2)+
-geom_point(data = limitslong, aes(x = Hz, y = Species),col = "black", size = 2)+
-scale_color_viridis()+
-scale_x_log10()+
-theme_classic()
-range<-ggplot(bound, aes_string(x = "x", y = "Species"))+   #, col= supraorder
-geom_path(data = bound2[c(seq(1, 95000, 10),95001,95002,95003), ],aes(col = `Threshold (dB)`), size = 2)+
-geom_point(data = limitslong, aes(x = Hz, y = Species),col = "black", size = 2)+
-scale_color_viridis()+
-scale_x_log10()+
-theme_classic()+
-coord_cartesian(clip = "off", ylim = c(1,22))+
-annotation_logticks(sides = "b", outside = TRUE, colour = "black")+
-geom_point(data = limits, aes(x = Hz, y = Species), shape = 21, size = 2, colour = "black", fill = "white")+
-#geom_vline(xintercept = min(limits$LowHzlimit))+
-#geom_vline(xintercept = max(limits$LowHzlimit))+
-#geom_hline(yintercept = 5)+
-ylab("Species")+
-xlab("Frequency(Hz)")+
-#theme(legend.position = "none")
-#ylim(c(-2,20))+
-geom_boxplot(data = limits,aes(x = LowHzlimit, y = 1), width = 0.8)+
-geom_boxplot(data = limits,aes(x =  besthz, y = 2), width = 0.8)+
-geom_boxplot(data = limits,aes(x = HighHzlimit, y = 3), width = 0.8)+
-geom_hline(yintercept = 3.5)+
-theme(axis.title.y = element_text(angle= 0, vjust = 0.5, hjust=1),
-axis.text.x = element_text(angle= 0, vjust = -2.5, hjust=0.5))
-range
-range<-ggplot(bound, aes_string(x = "x", y = "Species"))+   #, col= supraorder
-geom_path(data = bound2[c(seq(1, 95000, 10),95001,95002,95003), ],aes(col = `Threshold (dB)`), size = 2)+
-geom_point(data = limitslong, aes(x = Hz, y = Species),col = "black", size = 2)+
-scale_color_viridis()+
-scale_x_log10()+
-theme_classic()+
-coord_cartesian(clip = "off", ylim = c(1,22))+
-annotation_logticks(sides = "b", outside = TRUE, colour = "black")+
-#geom_point(data = limits, aes(x = Hz, y = Species), shape = 21, size = 2, colour = "black", fill = "white")+
-#geom_vline(xintercept = min(limits$LowHzlimit))+
-#geom_vline(xintercept = max(limits$LowHzlimit))+
-#geom_hline(yintercept = 5)+
-ylab("Species")+
-xlab("Frequency(Hz)")+
-#theme(legend.position = "none")
-#ylim(c(-2,20))+
-geom_boxplot(data = limits,aes(x = LowHzlimit, y = 1), width = 0.8)+
-geom_boxplot(data = limits,aes(x =  besthz, y = 2), width = 0.8)+
-geom_boxplot(data = limits,aes(x = HighHzlimit, y = 3), width = 0.8)+
-geom_hline(yintercept = 3.5)+
-theme(axis.title.y = element_text(angle= 0, vjust = 0.5, hjust=1),
-axis.text.x = element_text(angle= 0, vjust = -2.5, hjust=0.5))
-range
-range<-ggplot(bound, aes_string(x = "x", y = "Species"))+   #, col= supraorder
-geom_path(data = bound2[c(seq(1, 95000, 10),95001,95002,95003), ],aes(col = `Threshold (dB)`), size = 2)+
-geom_point(data = limitslong, aes(x = Hz, y = Species),col = "black", size = 2)+
-scale_color_viridis()+
-scale_x_log10()+
-theme_classic()+
-coord_cartesian(clip = "off", ylim = c(1,22))+
-annotation_logticks(sides = "b", outside = TRUE, colour = "black")+
-geom_point(data = limits, aes(x = besthz, y = Species), shape = 21, size = 2, colour = "black", fill = "white")+
-#geom_vline(xintercept = min(limits$LowHzlimit))+
-#geom_vline(xintercept = max(limits$LowHzlimit))+
-#geom_hline(yintercept = 5)+
-ylab("Species")+
-xlab("Frequency(Hz)")+
-#theme(legend.position = "none")
-#ylim(c(-2,20))+
-geom_boxplot(data = limits,aes(x = LowHzlimit, y = 1), width = 0.8)+
-geom_boxplot(data = limits,aes(x =  besthz, y = 2), width = 0.8)+
-geom_boxplot(data = limits,aes(x = HighHzlimit, y = 3), width = 0.8)+
-geom_hline(yintercept = 3.5)+
-theme(axis.title.y = element_text(angle= 0, vjust = 0.5, hjust=1),
-axis.text.x = element_text(angle= 0, vjust = -2.5, hjust=0.5))
-range
-range<-ggplot(bound, aes_string(x = "x", y = "Species"))+   #, col= supraorder
-geom_path(data = bound2[c(seq(1, 95000, 10),95001,95002,95003), ],aes(col = `Threshold (dB)`), size = 2)+
-geom_point(data = limitslong, aes(x = Hz, y = Species),col = "black", size = 2)+
-scale_color_viridis()+
-scale_x_log10()+
-theme_classic()+
-coord_cartesian(clip = "off", ylim = c(1,22))+
-annotation_logticks(sides = "b", outside = TRUE, colour = "black")+
-geom_point(data = limits, aes(x = besthz, y = Species), shape = 21, size = 2, colour = "black", fill = "white")+
-#geom_vline(xintercept = min(limits$LowHzlimit))+
-#geom_vline(xintercept = max(limits$LowHzlimit))+
-#geom_hline(yintercept = 5)+
-ylab("")+
-xlab("Frequency(Hz)")+
-#theme(legend.position = "none")
-#ylim(c(-2,20))+
-geom_boxplot(data = limits,aes(x = LowHzlimit, y = 1), width = 0.8)+
-geom_boxplot(data = limits,aes(x =  besthz, y = 2), width = 0.8)+
-geom_boxplot(data = limits,aes(x = HighHzlimit, y = 3), width = 0.8)+
-geom_hline(yintercept = 3.5)+
-theme(axis.title.y = element_text(angle= 0, vjust = 0.5, hjust=1),
-axis.text.x = element_text(angle= 0, vjust = -2.5, hjust=0.5))
-range
-range<-ggplot(bound, aes_string(x = "x", y = "Species"))+   #, col= supraorder
-geom_path(data = bound2[c(seq(1, 95000, 10),95001,95002,95003), ],aes(col = `Threshold (dB)`), size = 2)+
-geom_point(data = limitslong, aes(x = Hz, y = Species),col = "black", size = 2)+
-scale_color_viridis()+
-scale_x_log10()+
-theme_classic()+
-coord_cartesian(clip = "off", ylim = c(1,22))+
-annotation_logticks(sides = "b", outside = TRUE, colour = "black")+
-geom_point(data = limits, aes(x = besthz, y = Species), shape = 21, size = 2, colour = "black", fill = "white")+
-#geom_vline(xintercept = min(limits$LowHzlimit))+
-#geom_vline(xintercept = max(limits$LowHzlimit))+
-#geom_hline(yintercept = 5)+
-ylab("")+
-xlab("Frequency(Hz)")+
-#theme(legend.position = "none")
-#ylim(c(-2,20))+
-geom_boxplot(data = limits,aes(x = LowHzlimit, y = 1.5), width = 0.8)+
-geom_boxplot(data = limits,aes(x =  besthz, y = 2.5), width = 0.8)+
-geom_boxplot(data = limits,aes(x = HighHzlimit, y = 3.5), width = 0.8)+
-geom_hline(yintercept = 4)+
-theme(axis.title.y = element_text(angle= 0, vjust = 0.5, hjust=1),
-axis.text.x = element_text(angle= 0, vjust = -2.5, hjust=0.5))
-range
-range<-ggplot(bound, aes_string(x = "x", y = "Species"))+   #, col= supraorder
-geom_path(data = bound2[c(seq(1, 95000, 10),95001,95002,95003), ],aes(col = `Threshold (dB)`), size = 2)+
-geom_point(data = limitslong, aes(x = Hz, y = Species),col = "black", size = 2)+
-scale_color_viridis()+
-scale_x_log10()+
-theme_classic()+
-coord_cartesian(clip = "off", ylim = c(1,23))+
-annotation_logticks(sides = "b", outside = TRUE, colour = "black")+
-geom_point(data = limits, aes(x = besthz, y = Species), shape = 21, size = 2, colour = "black", fill = "white")+
-#geom_vline(xintercept = min(limits$LowHzlimit))+
-#geom_vline(xintercept = max(limits$LowHzlimit))+
-#geom_hline(yintercept = 5)+
-ylab("")+
-xlab("Frequency(Hz)")+
-#theme(legend.position = "none")
-#ylim(c(-2,20))+
-geom_boxplot(data = limits,aes(x = LowHzlimit, y = 1.5), width = 0.8)+
-geom_boxplot(data = limits,aes(x =  besthz, y = 2.5), width = 0.8)+
-geom_boxplot(data = limits,aes(x = HighHzlimit, y = 3.5), width = 0.8)+
-geom_hline(yintercept = 4)+
-theme(axis.title.y = element_text(angle= 0, vjust = 0.5, hjust=1),
-axis.text.x = element_text(angle= 0, vjust = -2.5, hjust=0.5))
-range
-range<-ggplot(bound, aes_string(x = "x", y = "Species"))+   #, col= supraorder
-geom_path(data = bound2[c(seq(1, 95000, 10),95001,95002,95003), ],aes(col = `Threshold (dB)`), size = 2)+
-geom_point(data = limitslong, aes(x = Hz, y = Species),col = "black", size = 2)+
-scale_color_viridis()+
-scale_x_log10()+
-theme_classic()+
-coord_cartesian(clip = "off", ylim = c(1,23))+
-annotation_logticks(sides = "b", outside = TRUE, colour = "black")+
-geom_point(data = limits, aes(x = besthz, y = Species), shape = 21, size = 2, colour = "black", fill = "white")+
-#geom_vline(xintercept = min(limits$LowHzlimit))+
-#geom_vline(xintercept = max(limits$LowHzlimit))+
-#geom_hline(yintercept = 5)+
-ylab("")+
-xlab("Frequency(Hz)")+
-#theme(legend.position = "none")
-#ylim(c(-2,20))+
-geom_boxplot(data = limits,aes(x = LowHzlimit, y = 0.5), width = 0.8)+
-geom_boxplot(data = limits,aes(x =  besthz, y = 1.5), width = 0.8)+
-geom_boxplot(data = limits,aes(x = HighHzlimit, y = 2.5), width = 0.8)+
-geom_hline(yintercept = 3)+
-theme(axis.title.y = element_text(angle= 0, vjust = 0.5, hjust=1),
-axis.text.x = element_text(angle= 0, vjust = -2.5, hjust=0.5))
-range
-range<-ggplot(bound, aes_string(x = "x", y = "Species"))+   #, col= supraorder
-geom_path(data = bound2[c(seq(1, 95000, 10),95001,95002,95003), ],aes(col = `Threshold (dB)`), size = 2)+
-geom_point(data = limitslong, aes(x = Hz, y = Species),col = "black", size = 2)+
-scale_color_viridis()+
-scale_x_log10()+
-theme_classic()+
-coord_cartesian(clip = "off", ylim = c(0,22))+
-annotation_logticks(sides = "b", outside = TRUE, colour = "black")+
-geom_point(data = limits, aes(x = besthz, y = Species), shape = 21, size = 2, colour = "black", fill = "white")+
-#geom_vline(xintercept = min(limits$LowHzlimit))+
-#geom_vline(xintercept = max(limits$LowHzlimit))+
-#geom_hline(yintercept = 5)+
-ylab("")+
-xlab("Frequency(Hz)")+
-#theme(legend.position = "none")
-#ylim(c(-2,20))+
-geom_boxplot(data = limits,aes(x = LowHzlimit, y = 0.5), width = 0.8)+
-geom_boxplot(data = limits,aes(x =  besthz, y = 1.5), width = 0.8)+
-geom_boxplot(data = limits,aes(x = HighHzlimit, y = 2.5), width = 0.8)+
-geom_hline(yintercept = 3)+
-theme(axis.title.y = element_text(angle= 0, vjust = 0.5, hjust=1),
-axis.text.x = element_text(angle= 0, vjust = -2.5, hjust=0.5))
-range
-#View(bound2[seq(1, nrow(bound2), 10), ])
-range<-ggplot(bound, aes_string(x = "x", y = "Species"))+   #, col= supraorder
-geom_path(data = bound2[c(seq(1, 95000, 10),95001,95002,95003), ],aes(col = `Threshold (dB)`), size = 2)+
-geom_point(data = limitslong, aes(x = Hz, y = Species),col = "black", size = 2)+
-scale_color_viridis()+
-scale_x_log10()+
+# filter(Coefficients!="(Intercept)" &
+#          P.val <0.05)
+#log transform anatomy data for the slope line
+#logged<-limits%>% mutate_at(vars(TM:UH),log)
+##########best Hz##############
+for(i in seq_along(anattraitx)){
+assign(paste0("slpline","_",as.character(i)),
+pgls_models_sig[i][[1]]$model$coef[1]+
+joined[,anattraitx[i]]*pgls_models_sig[i][[1]]$model$coef[2])
+}
+anattraitx2<-gsub("resid_log_", "Resid. ", anattraitx)
+anattraitx3<-gsub("_vslog_HM_", "",anattraitx2)
+runplot_audio<-function(e){
+p<-ggplot(joined,
+aes_string(x = spltmodel[[e]][2], y = spltmodel[[e]][1]))+
 theme_classic()+
-coord_cartesian(clip = "off", ylim = c(1,22))+
-annotation_logticks(sides = "b", outside = TRUE, colour = "black")+
-geom_point(data = limits, aes(x = besthz, y = Species), shape = 21, size = 2, colour = "black", fill = "white")+
-#geom_vline(xintercept = min(limits$LowHzlimit))+
-#geom_vline(xintercept = max(limits$LowHzlimit))+
-#geom_hline(yintercept = 5)+
-ylab("")+
-xlab("Frequency(Hz)")+
-#theme(legend.position = "none")
-#ylim(c(-2,20))+
-geom_boxplot(data = limits,aes(x = LowHzlimit, y = 1), width = 0.8)+
-geom_boxplot(data = limits,aes(x =  besthz, y = 2), width = 0.8)+
-geom_boxplot(data = limits,aes(x = HighHzlimit, y = 3), width = 0.8)+
-geom_hline(yintercept = 3.5)+
-theme(axis.title.y = element_text(angle= 0, vjust = 0.5, hjust=1),
-axis.text.x = element_text(angle= 0, vjust = -2.5, hjust=0.5))
-range
-range/bestsens+
-#plot_layout(widths = c(1,1))+
-plot_annotation(tag_levels = "A")
-range/bestsens+
-plot_layout(heights = c(2,1))+
-plot_annotation(tag_levels = "A")
-ggsave(file=paste0(choose.dir(),"/supplemental_bar.svg"), width=10, height=10)
-library(RColorBrewer)
-library(viridis)
-library(patchwork)
-library(tidyr)
+#theme(legend.position = "none")+
+#      axis.text.y = element_blank(),
+#      axis.title.y = element_blank())+
+geom_point(aes_string(shape="aud_rel"), size = 2)+
+geom_line(aes_string(x = anattraitx[e],
+y = paste0("slpline_",as.character(e))),
+col = "black", size = 2)+
+xlab(anattraitx3[e])+{
+if(e<11)
+ylab("Best Sensitivity\n (dB)")
+#else if(e >10 & e < 13)
+#  ylab("Best Frequency\n (Hz)")
+else if(e >10 & e < 21)
+ylab("High Frequency\n Limit (Hz)")
+else if(e > 20)
+ylab("Low Frequency\n Limit (Hz)")
+}
+p
+}
+runplot_audio(1)
+#load main dataframe
+df<-read.csv("databmadded.csv", stringsAsFactors = F, header = T) #, stringsAsFactors = FALSE
+#The pgls model function, which will be applied to list of formulas
+pgls_models<-function(i){
+pglsfit<-pgls(as.formula(i), data = birdCDO, #check comparative data object here<---
+lambda = 'ML', #find lambda using maximum likelihood
+bounds = list(lambda=c(0.0001,1)))#####
+}
+#Load phylogeny and correct names that were different between birdtree.org and the up-to-date species names
+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")#
+df$HM#with imputed values
+#Since PGLS uses one point per species,I make the dataframe to have average values for species with more than one specimen:
+#First I make a dataframe with only one species per line
+distinctdf<-distinct(df, Binomial, .keep_all = TRUE)
+distinctdforder<-arrange(distinctdf,Binomial)#sort by species name
+#Next get averages by species for columns with continuous data
+avgdf<-df %>% group_by(Binomial) %>%
+summarise_at(vars(Skull.width..mm.:TM_FP),mean, na.rm = TRUE)
+avgdf<-as.data.frame(avgdf)
+#Columns from the distinctdf dataframe which don't require averaging are added back
+avgdf$Species<-distinctdforder$Species
+avgdf$Low.Hz<-distinctdforder$Low.Hz
+avgdf$Order<-distinctdforder$Order
+avgdf$Family<-distinctdforder$Family
+avgdf$Category<-as.character(distinctdforder$Category)
+avgdf$birdtree<-gsub(" ","_",distinctdforder$Birdtree)
+avgdf$BM_lit<-distinctdforder$BM_lit
+avgdf$aud_spp<-distinctdforder$spp_audio
+avgdf$aud_rel<-distinctdforder$audio_relation
+avgdf$aud_spp<-distinctdforder$spp_audio
+#make comparative data frame object
+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
+######If doing audiogram analyses, you can now proceed to
+######'Audiograms linked to anatomy.R'
+######Ot
+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
+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
+birdCDO<-comparative.data(phy = birdtreels,data = avgdf,#[avgdf$Category!="Terrestrial",]
+names.col = Binomial,
+vcv = TRUE, na.omit = FALSE,
+warn.dropped = TRUE)
+#load main dataframe
+df<-read.csv("databmadded.csv", stringsAsFactors = F, header = T) #, stringsAsFactors = FALSE
+#The pgls model function, which will be applied to list of formulas
+pgls_models<-function(i){
+pglsfit<-pgls(as.formula(i), data = birdCDO, #check comparative data object here<---
+lambda = 'ML', #find lambda using maximum likelihood
+bounds = list(lambda=c(0.0001,1)))#####
+}
+#Load phylogeny and correct names that were different between birdtree.org and the up-to-date species names
+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")#
+df$HM#with imputed values
+#Since PGLS uses one point per species,I make the dataframe to have average values for species with more than one specimen:
+#First I make a dataframe with only one species per line
+distinctdf<-distinct(df, Binomial, .keep_all = TRUE)
+distinctdforder<-arrange(distinctdf,Binomial)#sort by species name
+#Next get averages by species for columns with continuous data
+avgdf<-df %>% group_by(Binomial) %>%
+summarise_at(vars(Skull.width..mm.:TM_FP),mean, na.rm = TRUE)
+avgdf<-as.data.frame(avgdf)
+#Columns from the distinctdf dataframe which don't require averaging are added back
+avgdf$Species<-distinctdforder$Species
+avgdf$Low.Hz<-distinctdforder$Low.Hz
+avgdf$Order<-distinctdforder$Order
+avgdf$Family<-distinctdforder$Family
+avgdf$Category<-as.character(distinctdforder$Category)
+avgdf$birdtree<-gsub(" ","_",distinctdforder$Birdtree)
+avgdf$BM_lit<-distinctdforder$BM_lit
+avgdf$aud_spp<-distinctdforder$spp_audio
+avgdf$aud_rel<-distinctdforder$audio_relation
+avgdf$aud_spp<-distinctdforder$spp_audio
+#make comparative data frame object
+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
+birdCDO<-comparative.data(phy = birdtreels,data = avgdf,#[avgdf$Category!="Terrestrial",]
+names.col = Binomial,
+vcv = TRUE, na.omit = FALSE,
+warn.dropped = TRUE)
+birdtreels$tip.label[14]<-"Corvus_cornix" #renamed from Corvus_albus
+birdtreels$tip.label[51]<-"Phalacrocorax_carbo" #rename "phalacrocorax_lucidus"
+birdCDO$dropped
+pgls_todo_hm<- c("log(TM_FP)~log(HM)",
+"log(COffset)~log(HM)",
+"log(UH)~log(HM)",
+"log(TMA)~log(HM)",
+"log(ECD)~log(HM)",
+"log(TM)~log(HM)",#
+"log(FP)~log(HM)",#
+"log(RW)~log(HM)",
+"log(ES)~log(HM)",
+"log(Air)~log(HM)",#
+"log(CL)~log(HM)",
+"log(CV)~log(HM)")
+geomcoefs<-c(0,#impedance-matching
+0.33,
+0.33,
+0,
+0.33,#auditory endorgan
+0.67,#input/output areas
+0.67,
+0.67,
+0.33,#stiffness
+1,
+0.33,#columella size
+1)
+#######functional category list#######
+categorylist<-c(rep("Impedance matching",4),
+"Auditory endorgan length",
+rep("Input/output areas",3),
+rep("Stiffness",2),
+rep("Columella size",2))
+#creates list of model outputs 'pgls_model_list'
+#dataframe with results 'hm'
+source("pgls_HM.R")
+pgls_models_list
+pgls_todo_hm
+getresids_as_df<-function(i){
+residtest<-as.data.frame(residuals(pgls_models_list[[i]]))
+residtest$resid_bname<-row.names(residtest)
+resid_measure<-function(){
+paste0("resid_",pgls_todo_hm[i])
+}
+residtest<-setNames(residtest,c(resid_measure(),"resid_bname"))
+}
+resids_df_list<-list()
+for(i in seq_along(pgls_todo_hm)){
+resids_df_list[[i]]<-assign("toadd",getresids_as_df(i))
+}
+#for(i in seq_along(resids_df_list)){
+#limit2<-limits
+#limit2<-full_join(limits,resids_df_list[[i]],by = c("spp_aud" = "resid_bname"))
+#
+#}
+#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"))%>%
+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"))%>%
+full_join(.,resids_df_list[[6]],by = c("spp_aud" = "resid_bname"))%>%
+full_join(.,resids_df_list[[7]],by = c("spp_aud" = "resid_bname"))%>%
+full_join(.,resids_df_list[[8]],by = c("spp_aud" = "resid_bname"))%>%
+full_join(.,resids_df_list[[9]],by = c("spp_aud" = "resid_bname"))%>%
+full_join(.,resids_df_list[[10]],by = c("spp_aud" = "resid_bname"))%>%
+full_join(.,resids_df_list[[11]],by = c("spp_aud" = "resid_bname"))%>%
+full_join(.,resids_df_list[[12]],by = c("spp_aud" = "resid_bname"))
+#only keep audiogram species
+joined<-joined[which(!is.na(joined$aud_rel)),]
+names(joined)
+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"))%>%
+full_join(.,resids_df_list[[6]],by = c("spp_aud" = "resid_bname"))%>%
+full_join(.,resids_df_list[[7]],by = c("spp_aud" = "resid_bname"))%>%
+full_join(.,resids_df_list[[8]],by = c("spp_aud" = "resid_bname"))%>%
+full_join(.,resids_df_list[[9]],by = c("spp_aud" = "resid_bname"))%>%
+full_join(.,resids_df_list[[10]],by = c("spp_aud" = "resid_bname"))%>%
+full_join(.,resids_df_list[[11]],by = c("spp_aud" = "resid_bname"))%>%
+full_join(.,resids_df_list[[12]],by = c("spp_aud" = "resid_bname"))
 library(ggrepel)
-##plot graph after creating the 'limits' dataframe (see 'Audiograms linked to anatomy' file)
-df_audiogrm_lst<-list()
+library(ggplot2)
+library(ggpubr)
+library(flextable)
+library(officer)
+library(dplyr)
+library(PerformanceAnalytics)
+phalacrocoraxavg<-avgdf[grepl('Phalacrocorax', avgdf$Binomial), ] %>%
+dplyr::select(where(is.numeric)) %>%
+summarise_all(mean, na.rm=T)
+corvusavg<-avgdf[grepl('Corvus', avgdf$Binomial), ] %>%
+dplyr::select(where(is.numeric)) %>%
+summarise_all(mean, na.rm=T)
+names(corvusavg)
+cong_avg<-dplyr::bind_rows(avgdf,corvusavg,phalacrocoraxavg)
+cong_avg$Binomial[128]<-"Corvus_cornix"
+cong_avg$Binomial[129]<-"Phalacrocorax_carbo"
+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
+fig1<-read.csv(paste0(getwd(),"/audiograms.csv"), stringsAsFactors = FALSE)
+#check how many reach cutoff
+#species not reaching lower Hz limit:
+minsubset<-fig1 %>% group_by(Species) %>% filter(Hz == min(Hz))
+minsubset$reachcutoff<-ifelse(minsubset$Threshold <35, "under35","over35")
+table(minsubset$reachcutoff)#count number
+minsubset$Species[minsubset$reachcutoff=="under35"]
+#species not reaching upper Hz limit
+maxsubset<-fig1 %>% group_by(Species) %>% filter(Hz == max(Hz))
+maxsubset$reachcutoff<-ifelse(maxsubset$Threshold <35, "under35","over35")
+table(maxsubset$reachcutoff)#count number
+maxsubset$Species[maxsubset$reachcutoff=="under35"]
+# get the high and low Hz limits from a cutoff level ----------------------
+splt<-split(fig1,fig1$Species)
+#set cutoff for the high and low Hz limits (35 dB)
+cutoff<-35
+#create new matrix to populate with data and convert to data audiogramram
+limits<-matrix(nrow=length(splt),ncol = 7)
 for(i in seq_along(splt)){
-df_audiogrm_lst[[i]]<-as.data.frame(approx(splt[[i]]$Hz,splt[[i]]$Threshold,n = 5000))
-df_audiogrm_lst[[i]]$Species<- rep(splt[[i]]$Species,length.out = 5000)
+audiogram<-data.frame()#
+df_audiogram<-as.data.frame(approx(splt[[i]]$Hz,splt[[i]]$Threshold,n = 5000))#approx function to interpolate audiogram
+df_audiogram
+#df_audiogram$y is sound level (dB)
+#df_audiogram$x is frequency (Hz)
+besthz<-df_audiogram$x[df_audiogram$y==min(df_audiogram$y)]
+bestsensitivity<-df_audiogram$y[df_audiogram$y==min(df_audiogram$y)]
+#calcualte low Hz limit
+if(nrow(df_audiogram[df_audiogram$y>cutoff & df_audiogram$x <besthz,])==0){#if the audiogram does not go above cutoff value, get minimum frequency tested
+lowlimit<-min(df_audiogram$x)
 }
-#single y example
-ggplot(df_audiogram, aes(x = x, y = y), col = y)+
-geom_point(aes(y = 60, col = y))+
-scale_color_viridis()
-#bind audiograms together
-bound<-do.call(rbind.data.frame,df_audiogrm_lst)
-ggplot(bound, aes(x = x, y = Species), col = y)+
-scale_x_log10()+
-geom_point()+
-geom_point(aes(col = y))+
-scale_color_viridis()
-#limits$Hz<-limits$besthz
-#dataframe holding the range
-tolong<-limits %>% select(Species,LowHzlimit,HighHzlimit,besthz)
-limitslong<-tolong %>% gather(key = "limit", value = "Hz", -Species)
-bound$LowHzlimit<-NA
-bound$HighHzlimit<-NA
-bound$besthz<-NA
-bound$bestsensitivity<-NA
-#add best hz to interpolated datset for sorting
-bound$besthz<-limits$besthz[match(bound$Species,limits$Species)]
-bound$bestsensitivity<-limits$bestsensitivity[match(bound$Species,limits$Species)]
-bound$HighHzlimit<-limits$HighHzlimit[match(bound$Species,limits$Species)]
-bound$LowHzlimit<-limits$LowHzlimit[match(bound$Species,limits$Species)]
-# Add new rows with audiogram metrics (best Hz, etc.)
-bound$Species<-as.factor(bound$Species)
-bound2<-bound
-levels(bound2$Species)<-c(levels(bound2$Species),"High freq. limit","Best freq.","Low freq. limit")
-bound2[95001,c("Species","besthz")]<-c("Low freq. limit",10)
-bound2[95002,c("Species","besthz")]<-c("Best freq.", 10)
-bound2[95003,c("Species","besthz")]<-c("High freq. limit" ,10)
-#bound$Species<-as.character(bound$Species)
-#bound$besthz<-as.numeric(bound$besthz)
-#reorder 'bound' df by besthz
-bound$Hz<-bound$x
-bound$Species = with(bound, reorder(Species, besthz, median))
-bound2$Species = with(bound2, reorder(Species, besthz, median))
-bound2$`Threshold (dB)`<-bound2$y
-#View(bound2[seq(1, nrow(bound2), 10), ])
-range<-ggplot(bound, aes_string(x = "x", y = "Species"))+   #, col= supraorder
-geom_path(data = bound2[c(seq(1, 95000, 10),95001,95002,95003), ],aes(col = `Threshold (dB)`), size = 2)+
-geom_point(data = limitslong, aes(x = Hz, y = Species),col = "black", size = 2)+
-scale_color_viridis()+
-scale_x_log10()+
-theme_classic()+
-coord_cartesian(clip = "off", ylim = c(1,22))+
-annotation_logticks(sides = "b", outside = TRUE, colour = "black")+
-geom_point(data = limits, aes(x = besthz, y = Species), shape = 21, size = 2, colour = "black", fill = "white")+
-#geom_vline(xintercept = min(limits$LowHzlimit))+
-#geom_vline(xintercept = max(limits$LowHzlimit))+
-#geom_hline(yintercept = 5)+
-ylab("")+
-xlab("Frequency(Hz)")+
-#theme(legend.position = "none")
-#ylim(c(-2,20))+
-geom_boxplot(data = limits,aes(x = LowHzlimit, y = 1), width = 0.8)+
-geom_boxplot(data = limits,aes(x =  besthz, y = 2), width = 0.8)+
-geom_boxplot(data = limits,aes(x = HighHzlimit, y = 3), width = 0.8)+
-geom_hline(yintercept = 3.5)+
-theme(axis.title.y = element_text(angle= 0, vjust = 0.5, hjust=1),
-axis.text.x = element_text(angle= 0, vjust = -2.5, hjust=0.5))
-range
-bound$Hz<-bound$x
-bound$Species = with(bound, reorder(Species, besthz, median))
-bound2$Species = with(bound2, reorder(Species, besthz, median))
-bound2$`Threshold (dB)`<-bound2$y
-range<-ggplot(bound, aes_string(x = "x", y = "Species"))+   #, col= supraorder
-geom_path(data = bound2[c(seq(1, 95000, 10),95001,95002,95003), ],aes(col = `Threshold (dB)`), size = 2)+
-geom_point(data = limitslong, aes(x = Hz, y = Species),col = "black", size = 2)+
-scale_color_viridis()+
-scale_x_log10()+
-theme_classic()+
-coord_cartesian(clip = "off", ylim = c(1,22))+
-annotation_logticks(sides = "b", outside = TRUE, colour = "black")+
-geom_point(data = limits, aes(x = besthz, y = Species), shape = 21, size = 2, colour = "black", fill = "white")+
-#geom_vline(xintercept = min(limits$LowHzlimit))+
-#geom_vline(xintercept = max(limits$LowHzlimit))+
-#geom_hline(yintercept = 5)+
-ylab("")+
-xlab("Frequency(Hz)")+
-#theme(legend.position = "none")
-#ylim(c(-2,20))+
-geom_boxplot(data = limits,aes(x = LowHzlimit, y = 1), width = 0.8)+
-geom_boxplot(data = limits,aes(x =  besthz, y = 2), width = 0.8)+
-geom_boxplot(data = limits,aes(x = HighHzlimit, y = 3), width = 0.8)+
-geom_hline(yintercept = 3.5)+
-theme(axis.title.y = element_text(angle= 0, vjust = 0.5, hjust=1),
-axis.text.x = element_text(angle= 0, vjust = -2.5, hjust=0.5))
-range
-#we will append data from the 'limits' df using species as a key
-bound$LowHzlimit<-NA
-bound$HighHzlimit<-NA
-bound$besthz<-NA
-bound$bestsensitivity<-NA
-#add best hz to interpolated datset for sorting
-bound$besthz<-limits$besthz[match(bound$Species,limits$Species)]
-bound$bestsensitivity<-limits$bestsensitivity[match(bound$Species,limits$Species)]
-bound$HighHzlimit<-limits$HighHzlimit[match(bound$Species,limits$Species)]
-bound$LowHzlimit<-limits$LowHzlimit[match(bound$Species,limits$Species)]
-bound$Species<-as.factor(bound$Species)
-bound2<-bound
-levels(bound2$Species)<-c(levels(bound2$Species),"High freq. limit","Best freq.","Low freq. limit")
-bound2[95001,c("Species","besthz")]<-c("Low freq. limit",10)
-bound2[95002,c("Species","besthz")]<-c("Best freq.", 10)
-bound2[95003,c("Species","besthz")]<-c("High freq. limit" ,10)
-bound$Hz<-bound$x
-bound$Species = with(bound, reorder(Species, besthz, median))
-bound2$Species = with(bound2, reorder(Species, besthz, median))
-bound2$`Threshold (dB)`<-bound2$y
-range<-ggplot(bound, aes_string(x = "x", y = "Species"))+   #, col= supraorder
-geom_path(data = bound2[c(seq(1, 95000, 10),95001,95002,95003), ],aes(col = `Threshold (dB)`), size = 2)+
-geom_point(data = limitslong, aes(x = Hz, y = Species),col = "black", size = 2)+
-scale_color_viridis()+
-scale_x_log10()+
-theme_classic()+
-coord_cartesian(clip = "off", ylim = c(1,22))+
-annotation_logticks(sides = "b", outside = TRUE, colour = "black")+
-geom_point(data = limits, aes(x = besthz, y = Species), shape = 21, size = 2, colour = "black", fill = "white")+
-#geom_vline(xintercept = min(limits$LowHzlimit))+
-#geom_vline(xintercept = max(limits$LowHzlimit))+
-#geom_hline(yintercept = 5)+
-ylab("")+
-xlab("Frequency(Hz)")+
-#theme(legend.position = "none")
-#ylim(c(-2,20))+
-geom_boxplot(data = limits,aes(x = LowHzlimit, y = 1), width = 0.8)+
-geom_boxplot(data = limits,aes(x =  besthz, y = 2), width = 0.8)+
-geom_boxplot(data = limits,aes(x = HighHzlimit, y = 3), width = 0.8)+
-geom_hline(yintercept = 3.5)+
-theme(axis.title.y = element_text(angle= 0, vjust = 0.5, hjust=1),
-axis.text.x = element_text(angle= 0, vjust = -2.5, hjust=0.5))
-range
-library(RColorBrewer)
-library(viridis)
-library(patchwork)
-library(tidyr)
-library(ggrepel)
-##plot graph after creating the 'limits' dataframe (see 'Audiograms linked to anatomy' file)
-df_audiogrm_lst<-list()
+else{
+lowflank<-df_audiogram[df_audiogram$y>cutoff & df_audiogram$x <besthz,]#get frequency where audiogram crosses cutoff value
+lowlimit<-max(lowflank$x)#lowhz limit
+}
+#calculate high Hz limit
+if(nrow(df_audiogram[df_audiogram$y>cutoff & df_audiogram$x >besthz,])==0){# #if the audiogram does not go above cutoff value, get max frequency tested
+highlimit<-max(df_audiogram$x)
+}
+else{
+highflank<-df_audiogram[df_audiogram$y>35 & df_audiogram$x >besthz,]#get frequency where audiogram crosses cutoff value
+highlimit<-min(highflank$x)#High hz limit
+}
+limits[i,1]<-lowlimit
+limits[i,2]<-highlimit
+limits[i,3]<-splt[[i]]$Species[1]
+limits[i,4]<-splt[[i]]$group[1]
+limits[i,5]<-splt[[i]]$Hz[1]
+limits[i,6]<-besthz
+limits[i,7]<-bestsensitivity
+}
+#View(limits)
+#convert to dataframe and give column names
+limits<-as.data.frame(limits)
+colnames(limits)<-c("LowHzlimit","HighHzlimit","Species","supraorder","Hz", "besthz","bestsensitivity")
+limits[,1]<-as.numeric(as.character(limits$LowHzlimit))
+limits[,2]<-as.numeric(as.character(limits$HighHzlimit))
+limits$Hz<-as.numeric(as.character(limits$Hz))
+limits$besthz<-as.numeric(as.character(limits$besthz))
+limits$bestsensitivity<-as.numeric(as.character(limits$bestsensitivity))
+###################add species from scan data that correspond with audiograms###############
+limits$binomial<-NA
+limits$binomial[limits$Species=="Barn owl"]<-"Tyto_alba"
+limits$binomial[limits$Species=="American kestrel"]<-"Falco_rupicolus" #
+limits$binomial[limits$Species=="Budgerigar"]<-"Melopsittacus_undulatus"
+limits$binomial[limits$Species=="Canary"]<-"Serinus_canaria"
+limits$binomial[limits$Species=="Chicken"]<-"Gallus_domesticus"
+limits$binomial[limits$Species=="Cockatiel"]<-"Nymphicus_hollandicus"
+limits$binomial[limits$Species=="Eurasian eagle owl"]<-"Bubo_africanus"
+limits$binomial[limits$Species=="Eurasian sparrowhawk"]<-"Accipiter_melanoleucus"
+limits$binomial[limits$Species=="Great cormorant"]<-"Phalacrocorax_carbo"
+limits$binomial[limits$Species=="Hooded crow"]<-"Corvus_cornix"
+limits$binomial[limits$Species=="Indian peafowl"]<-"Pavo_muticus"
+limits$binomial[limits$Species=="Mallard duck"]<-"Anas_georgica_georgica"
+limits$binomial[limits$Species=="Rock dove"]<-"Columba_livia"#
+limits$binomial[limits$Species=="Zebra finch"]<-"Taeniopygia_guttata"
+##################add anatomical data from anatomy df############
+limits$TM<-avgdf$TM[match(limits$binomial,avgdf$Binomial)]
+limits$RW<-avgdf$RW[match(limits$binomial,avgdf$Binomial)]
+limits$FP<-avgdf$FP[match(limits$binomial,avgdf$Binomial)]
+limits$Air<-avgdf$Air[match(limits$binomial,avgdf$Binomial)]
+limits$TM<-avgdf$TM[match(limits$binomial,avgdf$Binomial)]
+limits$HM<-avgdf$HM[match(limits$binomial,avgdf$Binomial)]
+limits$BM<-avgdf$BM_lit[match(limits$binomial,avgdf$Binomial)]
+limits$ES<-avgdf$ES[match(limits$binomial,avgdf$Binomial)]
+limits$TM_FP<-avgdf$TM_FP[match(limits$binomial,avgdf$Binomial)]
+limits$TMA<-avgdf$TMA[match(limits$binomial,avgdf$Binomial)]
+limits$COffset<-avgdf$COffset[match(limits$binomial,avgdf$Binomial)]
+limits$ECD<-avgdf$ECD[match(limits$binomial,avgdf$Binomial)]
+limits$CL<-avgdf$CL[match(limits$binomial,avgdf$Binomial)]
+limits$CV<-avgdf$CV[match(limits$binomial,avgdf$Binomial)]
+limits$UH<-avgdf$UH[match(limits$binomial,avgdf$Binomial)]
+limits$spp_aud<-avgdf$aud_spp[match(limits$binomial,avgdf$Binomial)]
+limits$aud_rel<-avgdf$aud_rel[match(limits$binomial,avgdf$Binomial)]
+#classification for the two species withaverage
+limits$aud_rel[limits$binomial=="Corvus_cornix"]<-"Congener"
+limits$spp_aud[limits$binomial=="Corvus_cornix"]<-"Corvus_cornix"
+limits$aud_rel[limits$binomial=="Phalacrocorax_carbo"]<-"Congener"
+limits$spp_aud[limits$binomial=="Phalacrocorax_carbo"]<-"Phalacrocorax_carbo"
+phalacrocoraxavg<-avgdf[grepl('Phalacrocorax', avgdf$Binomial), ] %>%
+dplyr::select(where(is.numeric)) %>%
+summarise_all(mean, na.rm=T)
+corvusavg<-avgdf[grepl('Corvus', avgdf$Binomial), ] %>%
+dplyr::select(where(is.numeric)) %>%
+summarise_all(mean, na.rm=T)
+names(corvusavg)
+cong_avg<-dplyr::bind_rows(avgdf,corvusavg,phalacrocoraxavg)
+cong_avg$Binomial[128]<-"Corvus_cornix"
+cong_avg$Binomial[129]<-"Phalacrocorax_carbo"
+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
+# load audiograms ---------------------------------------------------------
+fig1<-read.csv(paste0(getwd(),"/audiograms.csv"), stringsAsFactors = FALSE)
+#check how many reach cutoff
+#species not reaching lower Hz limit:
+minsubset<-fig1 %>% group_by(Species) %>% filter(Hz == min(Hz))
+minsubset$reachcutoff<-ifelse(minsubset$Threshold <35, "under35","over35")
+table(minsubset$reachcutoff)#count number
+minsubset$Species[minsubset$reachcutoff=="under35"]
+#species not reaching upper Hz limit
+maxsubset<-fig1 %>% group_by(Species) %>% filter(Hz == max(Hz))
+maxsubset$reachcutoff<-ifelse(maxsubset$Threshold <35, "under35","over35")
+table(maxsubset$reachcutoff)#count number
+maxsubset$Species[maxsubset$reachcutoff=="under35"]
+# get the high and low Hz limits from a cutoff level ----------------------
+splt<-split(fig1,fig1$Species)
+#set cutoff for the high and low Hz limits (35 dB)
+cutoff<-35
+#create new matrix to populate with data and convert to data audiogramram
+limits<-matrix(nrow=length(splt),ncol = 7)
 for(i in seq_along(splt)){
-df_audiogrm_lst[[i]]<-as.data.frame(approx(splt[[i]]$Hz,splt[[i]]$Threshold,n = 5000))
-df_audiogrm_lst[[i]]$Species<- rep(splt[[i]]$Species,length.out = 5000)
+audiogram<-data.frame()#
+df_audiogram<-as.data.frame(approx(splt[[i]]$Hz,splt[[i]]$Threshold,n = 5000))#approx function to interpolate audiogram
+df_audiogram
+#df_audiogram$y is sound level (dB)
+#df_audiogram$x is frequency (Hz)
+besthz<-df_audiogram$x[df_audiogram$y==min(df_audiogram$y)]
+bestsensitivity<-df_audiogram$y[df_audiogram$y==min(df_audiogram$y)]
+#calcualte low Hz limit
+if(nrow(df_audiogram[df_audiogram$y>cutoff & df_audiogram$x <besthz,])==0){#if the audiogram does not go above cutoff value, get minimum frequency tested
+lowlimit<-min(df_audiogram$x)
 }
-#single y audiogram example
-#ggplot(df_audiogram, aes(x = x, y = y), col = y)+
-#  geom_point(aes(y = 60, col = y))+
-#  scale_color_viridis()
-#bind audiograms together
-bound<-do.call(rbind.data.frame,df_audiogrm_lst)
-ggplot(bound, aes(x = x, y = Species), col = y)+
-scale_x_log10()+
-geom_point()+
-geom_point(aes(col = y))+
-scale_color_viridis()
-#convert "limits" dataframe to long format
-tolong<-limits %>% select(Species,LowHzlimit,HighHzlimit,besthz)
-limitslong<-tolong %>% gather(key = "limit", value = "Hz", -Species)
-#we will append data from the 'limits' df using species as a key
-bound$LowHzlimit<-NA
-bound$HighHzlimit<-NA
-bound$besthz<-NA
-bound$bestsensitivity<-NA
-#add best hz to interpolated datset for sorting
-bound$besthz<-limits$besthz[match(bound$Species,limits$Species)]
-bound$bestsensitivity<-limits$bestsensitivity[match(bound$Species,limits$Species)]
-bound$HighHzlimit<-limits$HighHzlimit[match(bound$Species,limits$Species)]
-bound$LowHzlimit<-limits$LowHzlimit[match(bound$Species,limits$Species)]
-# Add new rows with audiogram metrics (best Hz, etc.)
-bound$Species<-as.factor(bound$Species)
-bound2<-bound
-levels(bound2$Species)<-c(levels(bound2$Species),"High freq. limit","Best freq.","Low freq. limit")
-bound2[95001,c("Species","besthz")]<-c("Low freq. limit",10)
-bound2[95002,c("Species","besthz")]<-c("Best freq.", 10)
-bound2[95003,c("Species","besthz")]<-c("High freq. limit" ,10)
-#bound$Species<-as.character(bound$Species)
-#bound$besthz<-as.numeric(bound$besthz)
-#reorder 'bound' df by besthz
-bound$Hz<-bound$x #give appropriate naming for x
-bound$Species = with(bound, reorder(Species, besthz, median))
-bound2$Species = with(bound2, reorder(Species, besthz, median))
-bound2$`Threshold (dB)`<-bound2$y
-#plot audiograms as bars, per species
-range<-ggplot(bound, aes_string(x = "x", y = "Species"))+   #, col= supraorder
-geom_path(data = bound2[c(seq(1, 95000, 10),95001,95002,95003), ],aes(col = `Threshold (dB)`), size = 2)+
-geom_point(data = limitslong, aes(x = Hz, y = Species),col = "black", size = 2)+
-scale_color_viridis()+
-scale_x_log10()+
-theme_classic()+
-coord_cartesian(clip = "off", ylim = c(1,22))+
-annotation_logticks(sides = "b", outside = TRUE, colour = "black")+
-geom_point(data = limits, aes(x = besthz, y = Species), shape = 21, size = 2, colour = "black", fill = "white")+
-#geom_vline(xintercept = min(limits$LowHzlimit))+
-#geom_vline(xintercept = max(limits$LowHzlimit))+
-#geom_hline(yintercept = 5)+
-ylab("")+
-xlab("Frequency(Hz)")+
-#theme(legend.position = "none")
-#ylim(c(-2,20))+
-geom_boxplot(data = limits,aes(x = LowHzlimit, y = 1), width = 0.8)+
-geom_boxplot(data = limits,aes(x =  besthz, y = 2), width = 0.8)+
-geom_boxplot(data = limits,aes(x = HighHzlimit, y = 3), width = 0.8)+
-geom_hline(yintercept = 3.5)+
-theme(axis.title.y = element_text(angle= 0, vjust = 0.5, hjust=1),
-axis.text.x = element_text(angle= 0, vjust = -2.5, hjust=0.5))
-range
-bound$LowHzlimit<-NA
-bound$HighHzlimit<-NA
-bound$besthz<-NA
-bound$bestsensitivity<-NA
-bound$besthz<-limits$besthz[match(bound$Species,limits$Species)]
-bound$bestsensitivity<-limits$bestsensitivity[match(bound$Species,limits$Species)]
-bound$HighHzlimit<-limits$HighHzlimit[match(bound$Species,limits$Species)]
-bound$LowHzlimit<-limits$LowHzlimit[match(bound$Species,limits$Species)]
-bound$Species<-as.factor(bound$Species)
-bound2<-bound
-levels(bound2$Species)<-c(levels(bound2$Species),"High freq. limit","Best freq.","Low freq. limit")
-bound2[95001,c("Species","besthz")]<-c("Low freq. limit",10)
-bound2[95002,c("Species","besthz")]<-c("Best freq.", 10)
-bound2[95003,c("Species","besthz")]<-c("High freq. limit" ,10)
-ggplot(bound, aes(x = Hz, y = Species))
-#reorder 'bound' df by besthz
-bound$Hz<-bound$x #give appropriate naming for x
-bound$Species = with(bound, reorder(Species, besthz, median))
-bound2$Species = with(bound2, reorder(Species, besthz, median))
-bound2$`Threshold (dB)`<-bound2$y
-range<-ggplot(bound, aes(x = Hz, y = Species))+   #, col= supraorder
-geom_path(data = bound2[c(seq(1, 95000, 10),95001,95002,95003), ],aes(col = `Threshold (dB)`), size = 2)+
-geom_point(data = limitslong, aes(x = Hz, y = Species),col = "black", size = 2)+
-scale_color_viridis()+
-scale_x_log10()+
-theme_classic()+
-coord_cartesian(clip = "off", ylim = c(1,22))+
-annotation_logticks(sides = "b", outside = TRUE, colour = "black")+
-geom_point(data = limits, aes(x = besthz, y = Species), shape = 21, size = 2, colour = "black", fill = "white")+
-#geom_vline(xintercept = min(limits$LowHzlimit))+
-#geom_vline(xintercept = max(limits$LowHzlimit))+
-#geom_hline(yintercept = 5)+
-ylab("")+
-xlab("Frequency(Hz)")+
-#theme(legend.position = "none")
-#ylim(c(-2,20))+
-geom_boxplot(data = limits,aes(x = LowHzlimit, y = 1), width = 0.8)+
-geom_boxplot(data = limits,aes(x =  besthz, y = 2), width = 0.8)+
-geom_boxplot(data = limits,aes(x = HighHzlimit, y = 3), width = 0.8)+
-geom_hline(yintercept = 3.5)+
-theme(axis.title.y = element_text(angle= 0, vjust = 0.5, hjust=1),
-axis.text.x = element_text(angle= 0, vjust = -2.5, hjust=0.5))
-range
+else{
+lowflank<-df_audiogram[df_audiogram$y>cutoff & df_audiogram$x <besthz,]#get frequency where audiogram crosses cutoff value
+lowlimit<-max(lowflank$x)#lowhz limit
+}
+#calculate high Hz limit
+if(nrow(df_audiogram[df_audiogram$y>cutoff & df_audiogram$x >besthz,])==0){# #if the audiogram does not go above cutoff value, get max frequency tested
+highlimit<-max(df_audiogram$x)
+}
+else{
+highflank<-df_audiogram[df_audiogram$y>35 & df_audiogram$x >besthz,]#get frequency where audiogram crosses cutoff value
+highlimit<-min(highflank$x)#High hz limit
+}
+limits[i,1]<-lowlimit
+limits[i,2]<-highlimit
+limits[i,3]<-splt[[i]]$Species[1]
+limits[i,4]<-splt[[i]]$group[1]
+limits[i,5]<-splt[[i]]$Hz[1]
+limits[i,6]<-besthz
+limits[i,7]<-bestsensitivity
+}
+limits<-as.data.frame(limits)
+colnames(limits)<-c("LowHzlimit","HighHzlimit","Species","supraorder","Hz", "besthz","bestsensitivity")
+limits[,1]<-as.numeric(as.character(limits$LowHzlimit))
+limits[,2]<-as.numeric(as.character(limits$HighHzlimit))
+limits$Hz<-as.numeric(as.character(limits$Hz))
+limits$besthz<-as.numeric(as.character(limits$besthz))
+limits$bestsensitivity<-as.numeric(as.character(limits$bestsensitivity))
+###################add species from scan data that correspond with audiograms###############
+limits$binomial<-NA
+limits$binomial[limits$Species=="Barn owl"]<-"Tyto_alba"
+limits$binomial[limits$Species=="American kestrel"]<-"Falco_rupicolus" #
+limits$binomial[limits$Species=="Budgerigar"]<-"Melopsittacus_undulatus"
+limits$binomial[limits$Species=="Canary"]<-"Serinus_canaria"
+limits$binomial[limits$Species=="Chicken"]<-"Gallus_domesticus"
+limits$binomial[limits$Species=="Cockatiel"]<-"Nymphicus_hollandicus"
+limits$binomial[limits$Species=="Eurasian eagle owl"]<-"Bubo_africanus"
+limits$binomial[limits$Species=="Eurasian sparrowhawk"]<-"Accipiter_melanoleucus"
+limits$binomial[limits$Species=="Great cormorant"]<-"Phalacrocorax_carbo"
+limits$binomial[limits$Species=="Hooded crow"]<-"Corvus_cornix"
+limits$binomial[limits$Species=="Indian peafowl"]<-"Pavo_muticus"
+limits$binomial[limits$Species=="Mallard duck"]<-"Anas_georgica_georgica"
+limits$binomial[limits$Species=="Rock dove"]<-"Columba_livia"#
+limits$binomial[limits$Species=="Zebra finch"]<-"Taeniopygia_guttata"
+##################add anatomical data from anatomy df############
+limits$TM<-avgdf$TM[match(limits$binomial,avgdf$Binomial)]
+limits$RW<-avgdf$RW[match(limits$binomial,avgdf$Binomial)]
+limits$FP<-avgdf$FP[match(limits$binomial,avgdf$Binomial)]
+limits$Air<-avgdf$Air[match(limits$binomial,avgdf$Binomial)]
+limits$TM<-avgdf$TM[match(limits$binomial,avgdf$Binomial)]
+limits$HM<-avgdf$HM[match(limits$binomial,avgdf$Binomial)]
+limits$BM<-avgdf$BM_lit[match(limits$binomial,avgdf$Binomial)]
+limits$ES<-avgdf$ES[match(limits$binomial,avgdf$Binomial)]
+limits$TM_FP<-avgdf$TM_FP[match(limits$binomial,avgdf$Binomial)]
+limits$TMA<-avgdf$TMA[match(limits$binomial,avgdf$Binomial)]
+limits$COffset<-avgdf$COffset[match(limits$binomial,avgdf$Binomial)]
+limits$ECD<-avgdf$ECD[match(limits$binomial,avgdf$Binomial)]
+limits$CL<-avgdf$CL[match(limits$binomial,avgdf$Binomial)]
+limits$CV<-avgdf$CV[match(limits$binomial,avgdf$Binomial)]
+limits$UH<-avgdf$UH[match(limits$binomial,avgdf$Binomial)]
+limits$spp_aud<-avgdf$aud_spp[match(limits$binomial,avgdf$Binomial)]
+limits$aud_rel<-avgdf$aud_rel[match(limits$binomial,avgdf$Binomial)]
+#classification for the two species withaverage
+limits$aud_rel[limits$binomial=="Corvus_cornix"]<-"Congener"
+limits$spp_aud[limits$binomial=="Corvus_cornix"]<-"Corvus_cornix"
+limits$aud_rel[limits$binomial=="Phalacrocorax_carbo"]<-"Congener"
+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)
+chart.Correlation(audlog, histogram = TRUE, method = "pearson")
+# p-values from correlation tests
+cor.test(aud_data$LowHzlimit, aud_data$HighHzlimit)
+cor.test(aud_data$LowHzlimit, aud_data$besthz)
+cor.test(aud_data$LowHzlimit, aud_data$bestsensitivity)
+cor.test(aud_data$HighHzlimit, aud_data$bestsensitivity)
+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"))%>%
+full_join(.,resids_df_list[[6]],by = c("spp_aud" = "resid_bname"))%>%
+full_join(.,resids_df_list[[7]],by = c("spp_aud" = "resid_bname"))%>%
+full_join(.,resids_df_list[[8]],by = c("spp_aud" = "resid_bname"))%>%
+full_join(.,resids_df_list[[9]],by = c("spp_aud" = "resid_bname"))%>%
+full_join(.,resids_df_list[[10]],by = c("spp_aud" = "resid_bname"))%>%
+full_join(.,resids_df_list[[11]],by = c("spp_aud" = "resid_bname"))%>%
+full_join(.,resids_df_list[[12]],by = c("spp_aud" = "resid_bname"))
+#only keep audiogram species
+joined<-joined[which(!is.na(joined$aud_rel)),]
+names(joined)
+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,
+warn.dropped = TRUE)
+#check any tips dropped between linking phylogeny and dataframe
+birdCDO$dropped
+modellist_bs<-paste0("bestsensitivity~",residlist)
+modellist_lf<-paste0("log(LowHzlimit)~",residlist)
+modellist_hf<-paste0("log(HighHzlimit)~",residlist)
+modellist_bh<-paste0("log(besthz)~",residlist)
+residlist
+categorylist_lf<-c("Impedance match",
+"Impedance match",
+"Impedance match",
+"Impedance match",
+"Auditory endorgan length",
+"Input/output areas",
+"Input/output areas",
+"Input/output areas",
+"Stiffness",
+"Stiffness",
+"Columella size",
+"Columella size")
+categorylist_bs<-categorylist_lf
+categorylist_bh<-categorylist_lf
+categorylist_hf<-categorylist_lf
+pgls_models(modellist_bs[[1]])
+source("pgls_audiogram_lf.R")
+source("pgls_audiogram_lf.R")
diff --git a/Audiograms linked to anatomy.R b/Audiograms linked to anatomy.R
index e998757..24cd8c4 100644
--- a/Audiograms linked to anatomy.R	
+++ b/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
 
diff --git a/README.md b/README.md
index 89208e6..2370e37 100644
--- a/README.md
+++ b/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.    
                       
 
 
diff --git a/Scaling_2021.Rproj b/Scaling_2022.Rproj
similarity index 100%
rename from Scaling_2021.Rproj
rename to Scaling_2022.Rproj
diff --git a/Set up data_scl.R b/Set up data_scl.R
index ed56ed4..e2164db 100644
--- a/Set up data_scl.R	
+++ b/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)
 
diff --git a/pgls_resids re headmass.R b/pgls_resids re headmass.R
index b9e9fdd..d2e8765 100644
--- a/pgls_resids re headmass.R	
+++ b/pgls_resids re headmass.R	
@@ -1,6 +1,14 @@
 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
@@ -8,12 +16,6 @@ library(ggrepel)
 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,