.Rhistory

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
############RUN PGLS############
####scaling vs head mass########
####list of pgls models to run (only models with head mass are used)####
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)")
####list of expected geometric coefficients for___###
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")
# list of models predicted by hm are in an object 'pgls_model_list'
pgls_models_list
#list of the original regressions here:
pgls_todo_hm
#########get residuals of a pgls of measure~head mass as a dataframe#####
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"))
}
#############create list of dataframes containing residuals##############
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)
#remove tilda from names to not mess up pgls formulas based on names
joined<-joined %>% rename_with(~ gsub("~", "vs", .x, fixed = TRUE))
residlist<-gsub("[()]","_",names(joined)[27:38])
joined<-setNames(joined,c(names(joined)[1:26],residlist))
names(joined)
#make the comparative data object for pgls point to your new
#residual data
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
############get model list#########3
modellist_bs<-paste0("bestsensitivity~",residlist)
modellist_lf<-paste0("log(LowHzlimit)~",residlist)
modellist_hf<-paste0("log(HighHzlimit)~",residlist)
modellist_bh<-paste0("log(besthz)~",residlist)
residlist
#note same object name for category list but different order relative to non-resid pgls
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[[2]])
###########recreate the pgls forumula list to include the residuals#################
####now the same scripts can be run for analysis for each part of the audiogram####
source("pgls_audiogram_lf.R")
# high frequency limit ----------------------------------------------------
source("pgls_audiogram_hf.R")
source("pgls_audiogram_bh.R")
source("pgls_audiogram_bs.R")
### combine results from regressions for each measure into a single datagrame
audio_pgls_results<-bind_rows(audiogrampgls_bh,
audiogrampgls_bs,
audiogrampgls_lf,
audiogrampgls_hf)
audio_pgls_results$CI95_low<-audio_pgls_results$Estimate-audio_pgls_results$`Std. Error`*1.96
audio_pgls_results$CI95_high<-audio_pgls_results$Estimate+audio_pgls_results$`Std. Error`*1.96
#############formatting table###############
#combine estimate +/- 95 CI into one cell
audio_pgls_results$pglsslope<-paste0(audio_pgls_results$Estimate," (",
format(round(audio_pgls_results$CI95_low, 3), nsmall = 3),
",",
format(round(audio_pgls_results$CI95_high, 3), nsmall = 3),
")")
#split up model column
spltmodel<-strsplit(audio_pgls_results$Model,"~")
audio_pgls_results$`Audiogram metric`<-unlist(lapply(spltmodel, `[[`, 1))
audio_pgls_results$anattraitx<-unlist(lapply(spltmodel, `[[`, 2))
#only keep significant relationships
audio_pgls_results<-audio_pgls_results %>%
select(`Audiogram metric`,
category,
Coefficients,
pglsslope,
Adj_Rsquared,
P.val,
Lambda)%>%
filter(Coefficients!="(Intercept)" &
P.val <0.05)
#arrange by audiogram metric, category, and R2
audio_pgls_results <-audio_pgls_results %>%
arrange(factor(audio_pgls_results$`Audiogram metric`),
factor(category,levels = c("Impedance match",
"Stiffness",
"Input/output areas",
"Auditory endorgan length",
"Columella size")),
desc(Adj_Rsquared))
# remove the "log" from 'Coefficients'
#audio_pgls_results$xmodel_nolog<-numeric()
for(i in seq_along(audio_pgls_results$Coefficients)){
audio_pgls_results$Coefficients[i]<-gsub("[\\(\\)]", "", regmatches(audio_pgls_results$Coefficients, gregexpr("\\(.*?\\)", audio_pgls_results$Coefficients))[[i]])
}
#visualize the table better using the flextable package
flexall<-flextable(audio_pgls_results) %>% add_header_lines(
values = "Table X. ") %>%
#bold(i = ~ P.val < 0.05) %>% # select columns add: j = ~ Coefficients + P.val
autofit()
flexall
library(caper)
library(phytools)
library(ape)
library(dplyr)
library(ggpubr)
library(ggplot2)
library(tidyr)
library(purrr)
#########Line 73 - scaling between ear structures##########
#########Line 173 - scaling vs head mass########
################Set up data############
#set working directory and load data
#setwd(choose.dir())
#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.01,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'
######Other
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)")
####list of expected geometric coefficients for___###
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")
#################formatting table
#remove intercept estimates, drop model column,
#only keep significant relationships
#combine estimate +/- 95 CI into one cell
hm$pglsslope<-paste0(hm$Estimate," (",
format(round(hm$CI95_low, 3), nsmall = 3),
",",
format(round(hm$CI95_high, 3), nsmall = 3),
")")
#split up the model formula to get x and y components
splitmodel<-strsplit(hm$Model,"~")
hm$ymodel<-map(splitmodel,1)#left side of formula
# remove the "log" from each character string
hm$ymodel_nolog<-numeric()
for(i in seq_along(hm$ymodel)){
hm$ymodel_nolog[i]<-gsub("[\\(\\)]", "", regmatches(hm$ymodel, gregexpr("\\(.*?\\)", hm$ymodel))[[i]])
}
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)){
hm$Coefficients[i]<-gsub("[\\(\\)]", "", regmatches(hm$Coefficients, gregexpr("\\(.*?\\)", hm$Coefficients))[[i]])
}
#sort table by category and then adjusted R2
hm$category<-as.factor(hm$category)
hm<-arrange(hm,factor(hm$category, levels = c(
"Columella size",
"Auditory endorgan length",
"Input/output areas",
"Stiffness",
"Impedance match")),desc(Adj_Rsquared))
hm$pval<-format(round(hm$pval, 3), nsmall = 3)
####visualize the table better using the flextable package
flexall<-flextable(hm) %>%
add_header_lines(  values = "Table X. Models for selection") %>%
autofit()
flexall
####scaling vs head mass########
#list of pgls models to run (only models with head mass are used)
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)")
####list of expected geometric coefficients for___###
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")
#################formatting table
#remove intercept estimates, drop model column,
#only keep significant relationships
#combine estimate +/- 95 CI into one cell
hm$pglsslope<-paste0(hm$Estimate," (",
format(round(hm$CI95_low, 3), nsmall = 3),
",",
format(round(hm$CI95_high, 3), nsmall = 3),
")")
#split up the model formula to get x and y components
splitmodel<-strsplit(hm$Model,"~")
hm$ymodel<-map(splitmodel,1)#left side of formula
# remove the "log" from each character string
hm$ymodel_nolog<-numeric()
for(i in seq_along(hm$ymodel)){
hm$ymodel_nolog[i]<-gsub("[\\(\\)]", "", regmatches(hm$ymodel, gregexpr("\\(.*?\\)", hm$ymodel))[[i]])
}
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)){
hm$Coefficients[i]<-gsub("[\\(\\)]", "", regmatches(hm$Coefficients, gregexpr("\\(.*?\\)", hm$Coefficients))[[i]])
}
#sort table by category and then adjusted R2
hm$category<-as.factor(hm$category)
hm<-arrange(hm,factor(hm$category, levels = c(
"Columella size",
"Auditory endorgan length",
"Input/output areas",
"Stiffness",
"Impedance match")),desc(Adj_Rsquared))
hm$pval<-format(round(hm$pval, 3), nsmall = 3)
####visualize the table better using the flextable package
flexall<-flextable(hm) %>%
add_header_lines(  values = "Table X. Models for selection") %>%
autofit()
flexall
hm$tval
####scaling vs head mass########
#list of pgls models to run (only models with head mass are used)
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)")
####list of expected geometric coefficients for___###
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")
#################formatting table
#remove intercept estimates, drop model column,
#only keep significant relationships
#combine estimate +/- 95 CI into one cell
hm$pglsslope<-paste0(hm$Estimate," (",
format(round(hm$CI95_low, 3), nsmall = 3),
",",
format(round(hm$CI95_high, 3), nsmall = 3),
")")
#split up the model formula to get x and y components
splitmodel<-strsplit(hm$Model,"~")
hm$ymodel<-map(splitmodel,1)#left side of formula
# remove the "log" from each character string
hm$ymodel_nolog<-numeric()
for(i in seq_along(hm$ymodel)){
hm$ymodel_nolog[i]<-gsub("[\\(\\)]", "", regmatches(hm$ymodel, gregexpr("\\(.*?\\)", hm$ymodel))[[i]])
}
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)){
hm$Coefficients[i]<-gsub("[\\(\\)]", "", regmatches(hm$Coefficients, gregexpr("\\(.*?\\)", hm$Coefficients))[[i]])
}
#sort table by category and then adjusted R2
hm$category<-as.factor(hm$category)
hm<-arrange(hm,factor(hm$category, levels = c(
"Columella size",
"Auditory endorgan length",
"Input/output areas",
"Stiffness",
"Impedance match")),desc(Adj_Rsquared))
hm$pval<-format(round(hm$pval, 3), nsmall = 3)
####visualize the table better using the flextable package
flexall<-flextable(hm) %>%
add_header_lines(  values = "Table X. Models for selection") %>%
autofit()
flexall
#options(scipen = 100)
library(flextable)
library(officer)
#pgls_todo_hm<-pgls_todo_nogeomet[seq(2,length(pgls_todo_nogeomet),2)]
#Head mass only
modellist<-pgls_todo_hm
pgls_models_list<-lapply(modellist,pgls_models)#run pgls
#make list of dataframes with the PGLS outputs.
tbllist_HM<-list()
for (i in seq_along(pgls_models_list)){#change th 'Model' colume in this as appropriate
tbllist_HM[[i]]<-as.data.frame(summary(pgls_models_list[[i]])$'coefficients')
tbllist_HM[[i]]$Adj_Rsquared<-summary(pgls_models_list[[i]])$'adj.r.squared'[1]#rsquared
tbllist_HM[[i]]$Model<-modellist[i]#formula<_____________________CHECK FORMULA LIST HERE is correct
tbllist_HM[[i]]$Lambda<-summary(pgls_models_list[[i]])$'param'[[2]]#lambda
tbllist_HM[[i]]$Fstat<-summary(pgls_models_list[[i]])$fstatistic[1]
tbllist_HM[[i]]$Fstat_numdf<-summary(pgls_models_list[[i]])$fstatistic[2]
tbllist_HM[[i]]$Fstat_dendf<-summary(pgls_models_list[[i]])$fstatistic[3]
tbllist_HM[[i]]$AICc<-pgls_models_list[[i]]$aicc[1]
tbllist_HM[[i]]$T_1<-(coef(pgls_models_list[[i]])[2]-1)/pgls_models_list[[i]]$sterr[2]
tbllist_HM[[i]]$T_0<-(coef(pgls_models_list[[i]])[2]-0)/pgls_models_list[[i]]$sterr[2]
tbllist_HM[[i]]$p_slope_re_1<- 2*pt(abs(tbllist_HM[[i]]$T_1), pgls_models_list[[i]]$n-2, lower.tail = FALSE)
tbllist_HM[[i]]$p_slope_re_0<- 2*pt(abs(tbllist_HM[[i]]$T_0), pgls_models_list[[i]]$n-2, lower.tail = FALSE)
}
#organize the dataframe table (significant digist, remove redundant F stat & R squared)
for(i in seq_along(tbllist_HM)){
tbllist_HM[[i]]$Coefficients<-row.names(tbllist_HM[[i]])
tbllist_HM[[i]]$Coefficients<-gsub('[[:digit:]]+', '', tbllist_HM[[i]]$Coefficients)#regex to remove number automatically added during the loop
#identify numeric cols and character cols to apply the significant digits function
character_cols<-unlist(lapply(tbllist_HM[[i]], is.character))
numeric_cols <- unlist(lapply(tbllist_HM[[i]], is.numeric))# Identify numeric columns
tbllist_HM[[i]]<-cbind(tbllist_HM[[i]][,which(character_cols)],signif(tbllist_HM[[i]][,which(numeric_cols)], digits = 2))
#tbllist_HM[[i]] <- tbllist_HM[[i]][, c(6,11,8:10,7,5,1:4)]#change order of columns
#dplyr::select_if(tbllist_HM[[i]], is.numeric)#select only numeric data
colnames(tbllist_HM[[i]])[6]<-"P.val"#rename b/c flextable doesn't work will with the '>' sign
row.names(tbllist_HM[[i]])<-c()#remove row names
print(tbllist_HM[[i]])
}
hm<-do.call(rbind.data.frame,tbllist_HM)
hm$CI95_low<-hm$Estimate-hm$`Std. Error`*1.96
hm$CI95_high<-hm$Estimate+hm$`Std. Error`*1.96
hm$geometric_exp<-rep(geomcoefs,each =2)
hm$category<-rep(categorylist,each = 2)
hm$scalingtype<-ifelse(hm$CI95_high<hm$geometric_exp,"Hypo",
"other")
hm$scalingtype<-ifelse(hm$CI95_low>hm$geometric_exp,"Hyper",
hm$scalingtype)
hm$scalingtype<-ifelse(hm$CI95_high>hm$geometric_exp&hm$CI95_low<hm$geometric_exp,"Iso",
hm$scalingtype)
hm$tval<-(hm$Estimate-hm$geometric_exp)/hm$`Std. Error`#t-value of differnce between estimate and isometric slope
hm$pval<-2*pt(abs(hm$tval),df=hm$Fstat_dendf, lower.tail = FALSE)#two tailed p-val
#hm$pval<-round(hm$pval,digits = 5)
hm$tval
hm$Estimate