diff --git a/.RData b/.RData index e0e47a7..0048687 100644 Binary files a/.RData and b/.RData differ diff --git a/.RDataTmp b/.RDataTmp new file mode 100644 index 0000000..60686a3 Binary files /dev/null and b/.RDataTmp differ diff --git a/.Rhistory b/.Rhistory index 4c33caf..fec7190 100644 --- a/.Rhistory +++ b/.Rhistory @@ -1,512 +1,512 @@ -speciesPCAvalues$Binomial4<-ifelse(grepl("Vanellus", speciesPCAvalues$Binomial),as.character(speciesPCAvalues$Binomial),speciesPCAvalues$Binomial4) -speciesPCAvalues$Binomial4p5<-ifelse(grepl("Stercorarius_antarcticus", speciesPCAvalues$Binomial),as.character(speciesPCAvalues$Binomial),"") -speciesPCAvalues$Binomial4p5<-ifelse(grepl("Larus_dominicanus", speciesPCAvalues$Binomial),as.character(speciesPCAvalues$Binomial),speciesPCAvalues$Binomial4p5) -speciesPCAvalues$Binomial5<-ifelse(grepl("Alca", speciesPCAvalues$Binomial),as.character(speciesPCAvalues$Binomial),"") -speciesPCAvalues$Binomial5<-ifelse(grepl("Fratercula", speciesPCAvalues$Binomial),as.character(speciesPCAvalues$Binomial),speciesPCAvalues$Binomial5) -speciesPCAvalues$Binomial5<-ifelse(grepl("Cepphus", speciesPCAvalues$Binomial),as.character(speciesPCAvalues$Binomial),speciesPCAvalues$Binomial5) -speciesPCAvalues$Order3<-ifelse(grepl("Anseriformes", speciesPCAvalues$Order),as.character(speciesPCAvalues$Order),"") -speciesPCAvalues$divescore<-as.factor(speciesPCAvalues$divescore) -######################plot PC1 vs PC2 coloured by different factors -###############BIPLOT BASIC FUNCTION#################### -qualitative_hcl(5, palette = "Dark 3") -runPCAplot<-function(group, p1,p2,n1,n2){ -scattercat1<-ggplot(speciesPCAvalues, aes_string(x = p1, y = p2, label = "Binomial")) + +d$percentexplained<-d$`diag(pPCA$Eval)/sum(pPCA$Eval) * 100` +#eigenvalue plot +p<-ggplot(d, aes(x = reorder(PC1,-percentexplained), y = percentexplained))+ +geom_bar(stat = "identity", col = "white")+ +theme_minimal()+ +xlab("Principal component")+ +ylab("Percent variance \n explained") +p +#remove lavels for inset plot +p<-ggplot(d[1:6,], aes(x = reorder(PC1,-percentexplained), y = percentexplained))+ +geom_bar(stat = "identity")+ theme_classic()+ -xlab(label = paste0("p",as.character(p1),"(", -as.character(signif(d$percentexplained[n1], digits = 3)), -"%)"))+ -ylab(label = paste0("p",as.character(p2),"(", -as.character(signif(d$percentexplained[n2], digits = 3)), -"%)"))+ -#theme(legend.position = "none")+ -geom_encircle(data = speciesPCAvalues[which(is.na(speciesPCAvalues$Order2)),], -aes_string(),s_shape=1, expand=0, alpha = 0, col = "black")+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$Order=="Charadriiformes",], -aes_string(),s_shape=1, expand=0, fill = "#E16A86", alpha = 1, col = "black")+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$Order=="Procellariiformes",], -aes_string(),s_shape=1, expand=0, fill = "#00AA5A", alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$Order=="Suliformes",], -aes_string(),s_shape=1, expand=0, fill = "#B675E0", alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$Order=="Sphenisciformes",], -aes_string(),s_shape=1, expand=0, fill = "#AA9000", alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$Order=="Anseriformes",], -aes_string(),s_shape=1, expand=0, fill = "#00A6CA",alpha = 1)+ -geom_encircle(aes_string(fill = group),s_shape=1, expand=0, color = "black", alpha = 0.5)+#s_shape = 1 and expan = 0 are convex hull -geom_point(aes_string(color = group), shape = 21, size = 3, color = "black")+ -#scale_color_manual(values = alpha(c("black","black","black","black","black","black",0.2)) -geom_point(aes_string(fill = group), size = 3, shape = 21, col = "black")+ -scale_color_manual(values = c("#00AA5A","#E16A86","#AA9000", "#B675E0", "#00A6CA"), na.value = "white")+ -theme(legend.position = "none") +xlab("")+ +ylab("")+ +theme(axis.text.x = element_blank(), +axis.text.x = element_text(angle = 90), +panel.background = element_blank(), +panel.grid = element_blank(), +axis.ticks.y.left = element_line(color = "black")) +p +p<-ggplot(d[1:6,], aes(x = reorder(PC1,-percentexplained), y = percentexplained))+ +geom_bar(stat = "identity")+ +theme_classic()+ +xlab("")+ +ylab("")+ +theme(axis.text.x = element_blank(), +axis.text.x = element_text(angle = 90), +panel.background = element_blank(), +panel.grid = element_blank(), +axis.ticks.y.left = element_line(color = "black")) +p +p<-ggplot(d[1:6,], aes(x = reorder(PC1,-percentexplained), y = percentexplained))+ +geom_bar(stat = "identity")+ +theme_classic()+ +xlab("")+ +ylab("")+ +theme(axis.text.x = element_blank(), +axis.text.x = element_text(angle = 90), +panel.background = element_blank(), +panel.grid = element_blank()) +p +library(caper) +library(phytools) +library(ape) +library(dplyr) +library(ggpubr) +library(ggplot2) +library(tidyr) +#set working directory and load data +setwd("E:/ATclone/A_T-stats") +#load main dataframe +df<-read.csv("Eardata.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))) } -runPCAplot("Order2","PC1","PC2",1,2) -#set up color palette +#note some missing headmass values now imputed +df$Head.mass..g. +#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") +#add computed head mass from head mass~skullwidth pgls +source("SW_HM_.R")#add phylogeney here +df$Head.mass..g. +library(caper) +library(phytools) +library(ape) +library(dplyr) +library(ggpubr) +library(ggplot2) +library(tidyr) +#set working directory +setwd("E:/ATclone/A_T-stats") +#load main dataframe +df<-read.csv("Eardata.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))) +} +#note some missing headmass values +df$Head.mass..g. +#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") +#add computed head mass from head mass~skullwidth pgls +source("SW_HM_.R")#add phylogeney here +df$Head.mass..g. +#since PGLS uses one point per species, +#we create a dataframeaverage values for species with more than one specimen: +#first 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 for columns with continuous data +avgdf<-df %>% group_by(Binomial) %>% summarise_at(vars(Skull.width..mm.:area_ratio), +mean, na.rm = TRUE) +avgdf<-as.data.frame(avgdf) +#Nnow we add back columns from distinctdf which don't require averaging +avgdf$Species<-distinctdforder$Species +avgdf$Low.Hz<-distinctdforder$Low.Hz +avgdf$Order<-distinctdforder$Order +avgdf$Family<-distinctdforder$Family +avgdf$maxdivedepth<-distinctdforder$max +avgdf$Category<-as.character(distinctdforder$Category) +avgdf$birdtree<-gsub(" ","_",distinctdforder$Birdtree) +avgdf$IAC_detail<-distinctdforder$IAC_detail +avgdf$IBP_detail<-distinctdforder$IBP_detail +avgdf$Behind.TM<-distinctdforder$Behind.TM +avgdf$`fluid.filled.`<-distinctdforder$`fluid.filled.` +#add dive depth data to the main dataframe +#this script groups the 'surface foraging' based on more detailed ecologies +source("add_dive_depth_data.R") +#set 'Terrestrial' as reference level for the ecological grouping variable +avgdf$plungedistinct<-as.character(avgdf$catfeeding2) +avgdf$plungedistinct[which(is.na(avgdf$plungedistinct))]<-"Terrestrial" +avgdf$plungedistinct<-relevel(as.factor(avgdf$plungedistinct), ref = "Terrestrial") +avgdf$plungedistinct[avgdf$Binomial=="Ardea_melanocephala"]<-"Surface" +#make comparative data frame object for the pgls function +birdCDO<-comparative.data(phy = birdtreels,data = avgdf, +names.col = Binomial, +vcv = TRUE, na.omit = FALSE, +warn.dropped = TRUE) +#check any tips dropped between linking phylogeny and dataframe +birdCDO$dropped +#create list of pgls models to run (only models with head mass are used) +pgls_todo_nogeomet <- c("log(TMtotalarea)~log(Skull.width..mm.)", +"log(TMtotalarea)~log(Head.mass..g.)",# +"log(FPtotalarea)~log(Skull.width..mm.)", +"log(FPtotalarea)~log(Head.mass..g.)",# +"log(area_ratio)~log(Skull.width..mm.)", +"log(area_ratio)~log(Head.mass..g.)", +"log(dis_coltip_TMcentroid)~log(Skull.width..mm.)", +"log(dis_coltip_TMcentroid)~log(Head.mass..g.)", +"log(Umbo_distancetoTMplane)~log(Skull.width..mm.)", +"log(Umbo_distancetoTMplane)~log(Head.mass..g.)", +"log(meanTMangle)~log(Skull.width..mm.)", +"log(meanTMangle)~log(Head.mass..g.)", +"log(totalEClength)~log(Skull.width..mm.)", +"log(totalEClength)~log(Head.mass..g.)", +"log(RWtotalarea)~log(Skull.width..mm.)", +"log(RWtotalarea)~log(Head.mass..g.)", +"log(CAtotalarea)~log(Skull.width..mm.)", +"log(CAtotalarea)~log(Head.mass..g.)", +"log(Behind.TM)~log(Skull.width..mm.)", +"log(Behind.TM)~log(Head.mass..g.)",# +"log(Columella.length.mm)~log(Skull.width..mm.)", +"log(Columella.length.mm)~log(Head.mass..g.)", +"log(Columella.volume.mm3)~log(Skull.width..mm.)", +"log(Columella.volume.mm3)~log(Head.mass..g.)") +#select models with head mass +pgls_todo_hm<-pgls_todo_nogeomet[seq(2,length(pgls_todo_nogeomet),2)] +library(ggplot2) +library(RColorBrewer) +library(ggtree) +library(colorspace) +divecol<-c(sequential_hcl(6, palette = "Purple-Blue",rev = T)[2:6]) +#count # in each order +orderdf<-avgdf %>% count(Order) +#distinct df +g2<-avgdf[!duplicated(avgdf$Order),]#orders only +g2order<-arrange(g2,Order)#sort to match order df +g2order$n<-orderdf$n#attach number of species in order +str(g2order$Binomial) +#drop timps for order +g2order$full<-paste0(g2order$Order," (",as.character(g2order$n),")") +g2order$Order<-as.character(g2order$Order) +#fig sampling and grouping mypal <- colorRampPalette(brewer.pal(6, "Blues")) -qualitative_hcl(5, palette = "Dark 3") -#plot by orders -speciesPCAvalues$Order2 <- fct_explicit_na(as.factor(speciesPCAvalues$Order2)) -speciesPCAvalues$Order2 <-relevel(speciesPCAvalues$Order2,ref = "(Missing)") -speciesPCAvalues$Order2 <-factor(speciesPCAvalues$Order2,levels = c("(Missing)" , "Procellariiformes", "Charadriiformes" , -"Sphenisciformes" , "Suliformes", "Anseriformes" )) -order<-runPCAplot("Order2","PC1","PC2",1,2) + -geom_point(aes(fill = Order2), size = 3, shape = 21, col = "black")+ -scale_fill_manual(values = c("white","#00AA5A","#E16A86","#AA9000", "#B675E0", "#00A6CA"))+ -theme(legend.position = "none")+ -theme(legend.position = "right")+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$Order=="Anseriformes",], -aes_string(),s_shape=1, expand=0, fill = "#00AA5A",alpha = 0.5)+ -geom_text_repel(aes(label = Binomial5), -nudge_y = 30 - speciesPCAvalues$PC2)+ -geom_text_repel(aes(label = Binomial4), -nudge_y = 30 - speciesPCAvalues$PC2) -order -cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7") +mypal2 <- colorRampPalette(brewer.pal(6, "YlOrRd")) +#clade labels +A<-findMRCA(birdtreels,avgdf$Binomial[avgdf$Order=="Passeriformes"]) +B<-findMRCA(birdtreels,avgdf$Binomial[avgdf$Order=="Charadriiformes"]) +C<-findMRCA(birdtreels,avgdf$Binomial[avgdf$Order=="Falconiformes"]) +D<-findMRCA(birdtreels,avgdf$Binomial[avgdf$Order=="Accipitriformes"]) +E<-findMRCA(birdtreels,avgdf$Binomial[avgdf$Order=="Bucerotiformes"]) +F<-findMRCA(birdtreels,avgdf$Binomial[avgdf$Order=="Gruiformes"]) +G<--findMRCA(birdtreels,avgdf$Binomial[avgdf$Order=="Anseriformes"]) +H<-findMRCA(birdtreels,avgdf$Binomial[avgdf$Order=="Galliformes"]) +I<-findMRCA(birdtreels,avgdf$Binomial[avgdf$Order=="Pelecaniformes"]) +J<-findMRCA(birdtreels,avgdf$Binomial[avgdf$Order=="Procellariiformes"]) +K<-findMRCA(birdtreels,avgdf$Binomial[avgdf$Order=="Sphenisciformes"]) +L<-findMRCA(birdtreels,avgdf$Binomial[avgdf$Order=="Suliformes"]) +M<-findMRCA(birdtreels,avgdf$Binomial[avgdf$Order=="Psittaciformes"]) +N<-findMRCA(birdtreels,avgdf$Binomial[avgdf$Order=="Columbiformes"]) +O<-findMRCA(birdtreels,avgdf$Binomial[avgdf$Order=="Strigiformes"]) +P<-findMRCA(birdtreels,avgdf$Binomial[avgdf$Order=="Anseriformes"]) +p<-ggtree(birdtreels, layout = "circular", open.angle = 150) %<+% avgdf + ###########, layout = "circular" +#geom_tiplab(aes(label = Binomial2), linesize = 0.1, offset = 4) + #circular +#geom_text(aes(label = node))+ +#xlim(NA, 150) + +#ylim(NA,140) + +#geom_tippoint(aes(color = Category), size = 5, shape = 15) + +#geom_tippoint(aes(x = x+2.5,color = plungedistinct), size = 5, shape = 15) + +#geom_tippoint(aes(x = x+5, color = as.factor(divescore)), size = 5, shape = 15) + +scale_color_manual(values = c(mypal(5),"black","grey","green","blue"))+ +geom_cladelabel(A, "Passeriformes", offset=25, barsize=2, align = T, angle=0,offset.text=15, fontsize=3)+ +geom_cladelabel(B, "Charadriiformes", offset=25, barsize=2, align = T, angle=0, offset.text=15, fontsize=3)+ +geom_cladelabel(C, "Falconiformes", offset=25, barsize=2, angle=0, offset.text=15, fontsize=3)+ +geom_cladelabel(D, "Accipitriformes", offset=25, barsize=2, angle=0, offset.text=15, fontsize=3)+ +geom_cladelabel(E, "Bucerotiformes", offset=25, barsize=2, angle=0, offset.text=15, fontsize=3)+ +geom_cladelabel(F, "Gruiformes", offset=25, barsize=2, angle=0, offset.text=15, fontsize=3)+ +geom_cladelabel(G, "Anseriformes", offset=25, barsize=2, angle=0, offset.text=15, fontsize=3)+ +geom_cladelabel(H, "Galliformes", offset=25, barsize=2, angle=0, offset.text=15, fontsize=3)+ +geom_cladelabel(I, "Pelecaniformes", offset=25, barsize=2, angle=0, offset.text=15, fontsize=3)+ +geom_cladelabel(J, "Procellariiformes", offset=25, barsize=2, angle=0, offset.text=15, fontsize=3)+ +geom_cladelabel(K, "Sphenisciformes", offset=25, barsize=2, angle=0, offset.text=15, fontsize=3)+ +geom_cladelabel(L, "Suliformes", offset=25, barsize=2, angle=0, offset.text=15, fontsize=3)+ +geom_cladelabel(M, "Psittaciformes", offset=25, barsize=2, angle=0, offset.text=15, fontsize=3)+ +#geom_cladelabel(119, "Struthioniformes", offset=0, barsize=2, angle=0, offset.text=15, fontsize=3)+ +geom_cladelabel(N, "Columbiformes", offset=25, barsize=2, angle=0, offset.text=15, fontsize=3)+ +geom_cladelabel(O, "Strigiformes", offset=25, barsize=2, angle=0, offset.text=15, fontsize=3)+ +geom_cladelabel(P, "Anseriformes", offset=25, barsize=2, angle=0, offset.text=15, fontsize=3)+ +#geom_tiplab2(aes(size = 1, label = Binomialsinge), align = TRUE, geom = "text", angle = 0, offset=15, linetype = "dotted")+ +geom_strip("Phaethon_rubricauda", "Phaethon_rubricauda", offset=16, offset.text=15, hjust=0, fontsize=3, +label="Phaethontifomres")+ +geom_strip(g2order$Binomial[9], g2order$Binomial[9], offset=16, offset.text=15, hjust=0, fontsize=3, +label=g2order$full[9])+ +geom_strip(g2order$Binomial[20], g2order$Binomial[20], offset=16, offset.text=15, hjust=0, fontsize=3, +label=g2order$full[20])+ +geom_strip(g2order$Binomial[7], g2order$Binomial[7], offset=16, offset.text=15, hjust=0, fontsize=3, +label=g2order$full[7])+ +geom_strip(g2order$Binomial[3], g2order$Binomial[3], offset=16, offset.text=15, hjust=0, fontsize=3, +label=g2order$full[3])+ +geom_strip(g2order$Binomial[15], g2order$Binomial[15], offset=16, offset.text=15, hjust=0, fontsize=3, +label=g2order$full[15])+ +geom_strip(g2order$Binomial[10], g2order$Binomial[10], offset=16, offset.text=15, hjust=0, fontsize=3, +label=g2order$full[10])+ +geom_strip(g2order$Binomial[19], g2order$Binomial[19], offset=16, offset.text=15, hjust=0, fontsize=3, +label=g2order$full[19])+ +geom_strip(g2order$Binomial[25], g2order$Binomial[25], offset=16, offset.text=15, hjust=0, fontsize=3, +label=g2order$full[25])+ +geom_strip(g2order$Binomial[13], g2order$Binomial[13], offset=16, offset.text=15, hjust=0, fontsize=3, +label=g2order$full[13])+ +geom_strip(g2order$Binomial[5], g2order$Binomial[5], offset=16, offset.text=15, hjust=0, fontsize=3, +label=g2order$full[5]) +p +rownames(avgdf)<-avgdf$Binomial # The palette with black: cbbPalette <- c( "#FFFFFF","#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7") -runPCAplotPlunge<-function(group, p1,p2,n1,n2){ -scattercat1<-ggplot(speciesPCAvalues, aes_string(x = p1, y = p2, label = "Binomial")) + -theme_classic()+ -xlab(label = paste0("p",as.character(p1),"(", -as.character(signif(d$percentexplained[n1], digits = 3)), -"%)"))+ -ylab(label = paste0("p",as.character(p2),"(", -as.character(signif(d$percentexplained[n2], digits = 3)), -"%)"))+ -#theme(legend.position = "none")+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$plungedistinct =="Plunging",], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$plungedistinct =="Terrestrial",], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$plungedistinct =="Surface",], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$plungedistinct =="Underwater pursuit",], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(aes_string(fill = group),s_shape=1, expand=0, color = "black", alpha = 0.7)+#s_shape = 1 and expan = 0 are convex hull -geom_point(aes_string(color = group), shape = 21, size = 3, color = "black") -#scale_color_manual(values = alpha(c("black","black","black","black","black","black",0.2)) -scattercat1 -} -runPCAplotPlunge("plungedistinct","PC1","PC2",1,2) -#plot by aquatic groupings -plunge<-runPCAplotPlunge("plungedistinct","PC1","PC2",1,2)+ -geom_point(aes(fill = plungedistinct), size = 3, shape = 21, col = "black")+ -scale_fill_manual(values = cbbPalette)+ -theme(legend.position = "right")+ -geom_text(aes(label = Binomial2))+ -geom_text(aes(label = Binomial0)) -plunge -divecol<-c("white",sequential_hcl(6, palette = "Purple-Blue",rev = T)[2:6]) -speciesPCAvalues$divescore <- fct_explicit_na(speciesPCAvalues$divescore) -runPCAplotdive<-function(group, p1,p2,n1,n2){ -scattercat1<-ggplot(speciesPCAvalues, aes_string(x = p1, y = p2, label = "Binomial")) + -theme_classic()+ -xlab(label = paste0("p",as.character(p1),"(", -as.character(signif(d$percentexplained[n1], digits = 3)), -"%)"))+ -ylab(label = paste0("p",as.character(p2),"(", -as.character(signif(d$percentexplained[n2], digits = 3)), -"%)"))+ -#theme(legend.position = "none")+ -geom_encircle(data = speciesPCAvalues[which(is.na(speciesPCAvalues$divescore)),], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==0,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==1,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==2,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==3,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==4,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(aes_string(fill = group),s_shape=1, expand=0, color = "black", alpha = 0.7)+#s_shape = 1 and expan = 0 are convex hull -geom_point(aes_string(color = group), shape = 21, size = 3, color = "black") -#scale_color_manual(values = alpha(c("black","black","black","black","black","black",0.2)) -scattercat1 -} -runPCAplotdive("divescore","PC1","PC2",1,2) -divescore<-runPCAplotdive("divescore","PC1","PC2",1,2)+ -geom_point(aes(fill = divescore), size = 3, shape = 21,col = "black", alpha = 1)+ -scale_fill_manual(values = divecol)+ -#scale_alpha_manual(values = c(-.5,1,0,0,56)) -theme(legend.position = "right") -divescore -#mypal <- colorRampPalette(rev(brewer.pal(6, "Blues"))) -sequential_hcl(5, palette = "Purple-Blue", rev = T) -BLUE<-c("white","#d0d1e6", +BLUE<-c("#d0d1e6", "#a6bddb", "#74a9cf", "#2b8cbe", "#045a8d") -divecol<-c("white",sequential_hcl(6, palette = "Purple-Blue",rev = T)[2:6]) -divecol -runPCAplotdive<-function(group, p1,p2,n1,n2){ -scattercat1<-ggplot(speciesPCAvalues, aes_string(x = p1, y = p2, label = "Binomial")) + -theme_classic()+ -xlab(label = paste0("p",as.character(p1),"(", -as.character(signif(d$percentexplained[n1], digits = 3)), -"%)"))+ -ylab(label = paste0("p",as.character(p2),"(", -as.character(signif(d$percentexplained[n2], digits = 3)), -"%)"))+ -#theme(legend.position = "none")+ -geom_encircle(data = speciesPCAvalues[which(is.na(speciesPCAvalues$divescore)),], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==0,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==1,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==2,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==3,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==4,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(aes_string(fill = group),s_shape=1, expand=0, color = "black", alpha = 0.7)+#s_shape = 1 and expan = 0 are convex hull -geom_point(aes_string(color = group), shape = 21, size = 3, color = "black") -#scale_color_manual(values = alpha(c("black","black","black","black","black","black",0.2)) -scattercat1 -} -runPCAplotdive("divescore","PC1","PC2",1,2) -divescore<-runPCAplotdive("divescore","PC1","PC2",1,2)+ -geom_point(aes(fill = divescore), size = 3, shape = 21,col = "black", alpha = 1)+ -scale_fill_manual(values = divecol)+ -#scale_alpha_manual(values = c(-.5,1,0,0,56)) -theme(legend.position = "right") -divescore -runPCAplotdive("divescore","PC1","PC2",1,2) -divecol<-c(sequential_hcl(6, palette = "Purple-Blue",rev = T)[2:6],"white") -divescore<-runPCAplotdive("divescore","PC1","PC2",1,2)+ -geom_point(aes(fill = divescore), size = 3, shape = 21,col = "black", alpha = 1)+ -scale_fill_manual(values = divecol)+ -#scale_alpha_manual(values = c(-.5,1,0,0,56)) -theme(legend.position = "right") -divescore -runPCAplotdive<-function(group, p1,p2,n1,n2){ -scattercat1<-ggplot(speciesPCAvalues, aes_string(x = p1, y = p2, label = "Binomial")) + -theme_classic()+ -xlab(label = paste0("p",as.character(p1),"(", -as.character(signif(d$percentexplained[n1], digits = 3)), -"%)"))+ -ylab(label = paste0("p",as.character(p2),"(", -as.character(signif(d$percentexplained[n2], digits = 3)), -"%)"))+ -#theme(legend.position = "none")+ -geom_encircle(data = speciesPCAvalues[which(is.na(speciesPCAvalues$divescore)),], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==0,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==1,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==2,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==3,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==4,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(aes_string(fill = group),s_shape=1, expand=0, color = "black", alpha = 0.7)+#s_shape = 1 and expan = 0 are convex hull -geom_point(aes_string(color = group), shape = 21, size = 3, color = "black") -#scale_color_manual(values = alpha(c("black","black","black","black","black","black",0.2)) -scattercat1 -} -runPCAplotdive("divescore","PC1","PC2",1,2) -divescore<-runPCAplotdive("divescore","PC1","PC2",1,2)+ -geom_point(aes(fill = divescore), size = 3, shape = 21,col = "black", alpha = 1)+ -scale_fill_manual(values = divecol)+ -#scale_alpha_manual(values = c(-.5,1,0,0,56)) -theme(legend.position = "right") -divescore -runPCAplotdive("divescore","PC1","PC2",1,2) -runPCAplotdive<-function(group, p1,p2,n1,n2){ -scattercat1<-ggplot(speciesPCAvalues, aes_string(x = p1, y = p2, label = "Binomial")) + -theme_classic()+ -xlab(label = paste0("p",as.character(p1),"(", -as.character(signif(d$percentexplained[n1], digits = 3)), -"%)"))+ -ylab(label = paste0("p",as.character(p2),"(", -as.character(signif(d$percentexplained[n2], digits = 3)), -"%)"))+ -#theme(legend.position = "none")+ -geom_encircle(data = speciesPCAvalues[which(is.na(speciesPCAvalues$divescore)),], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==0,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==1,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==2,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==3,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==4,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(aes_string(fill = group),s_shape=1, expand=0, color = "black", alpha = 0.7)+#s_shape = 1 and expan = 0 are convex hull -#geom_point(aes_string(color = group), shape = 21, size = 3, color = "black") -#scale_color_manual(values = alpha(c("black","black","black","black","black","black",0.2)) -scattercat1 -} -runPCAplotdive("divescore","PC1","PC2",1,2) -runPCAplotdive<-function(group, p1,p2,n1,n2){ -scattercat1<-ggplot(speciesPCAvalues, aes_string(x = p1, y = p2, label = "Binomial")) + -theme_classic()+ -xlab(label = paste0("p",as.character(p1),"(", -as.character(signif(d$percentexplained[n1], digits = 3)), -"%)"))+ -ylab(label = paste0("p",as.character(p2),"(", -as.character(signif(d$percentexplained[n2], digits = 3)), -"%)"))+ -#theme(legend.position = "none")+ -geom_encircle(data = speciesPCAvalues[which(is.na(speciesPCAvalues$divescore)),], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==0,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==1,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==2,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==3,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==4,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(aes_string(fill = group),s_shape=1, expand=0, color = "black", alpha = 0.7)#s_shape = 1 and expan = 0 are convex hull -#geom_point(aes_string(color = group), shape = 21, size = 3, color = "black") -#scale_color_manual(values = alpha(c("black","black","black","black","black","black",0.2)) -scattercat1 -} -runPCAplotdive("divescore","PC1","PC2",1,2) -speciesPCAvalues$divescore <- fct_explicit_na(speciesPCAvalues$divescore) -speciesPCAvalues$divescore<-as.factor(dfwithresids[match(row.names(speciesPCAvalues), -dfwithresids$Binomial),"divescore"])#<-------------<------- -speciesPCAvalues$plungedistinct<-as.factor(dfwithresids[match(row.names(speciesPCAvalues), -dfwithresids$Binomial),"plungedistinct"])#<-------------<------- -divecol<-c(sequential_hcl(6, palette = "Purple-Blue",rev = T)[2:6],"white") -runPCAplotdive<-function(group, p1,p2,n1,n2){ -scattercat1<-ggplot(speciesPCAvalues, aes_string(x = p1, y = p2, label = "Binomial")) + -theme_classic()+ -xlab(label = paste0("p",as.character(p1),"(", -as.character(signif(d$percentexplained[n1], digits = 3)), -"%)"))+ -ylab(label = paste0("p",as.character(p2),"(", -as.character(signif(d$percentexplained[n2], digits = 3)), -"%)"))+ -#theme(legend.position = "none")+ -geom_encircle(data = speciesPCAvalues[which(is.na(speciesPCAvalues$divescore)),], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==0,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==1,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==2,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==3,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==4,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(aes_string(fill = group),s_shape=1, expand=0, color = "black", alpha = 0.7)#s_shape = 1 and expan = 0 are convex hull -#geom_point(aes_string(color = group), shape = 21, size = 3, color = "black") -#scale_color_manual(values = alpha(c("black","black","black","black","black","black",0.2)) -scattercat1 -} -runPCAplotdive("divescore","PC1","PC2",1,2) -runPCAplotdive<-function(group, p1,p2,n1,n2){ -scattercat1<-ggplot(speciesPCAvalues, aes_string(x = p1, y = p2, label = "Binomial")) + -theme_classic()+ -xlab(label = paste0("p",as.character(p1),"(", -as.character(signif(d$percentexplained[n1], digits = 3)), -"%)"))+ -ylab(label = paste0("p",as.character(p2),"(", -as.character(signif(d$percentexplained[n2], digits = 3)), -"%)"))+ -#theme(legend.position = "none")+ -geom_encircle(data = speciesPCAvalues[which(is.na(speciesPCAvalues$divescore)),], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==0,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==1,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==2,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==3,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==4,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(aes_string(fill = group),s_shape=1, expand=0, color = "black", alpha = 0.7)+#s_shape = 1 and expan = 0 are convex hull -geom_point(aes_string(color = group), shape = 21, size = 3, color = "black") -#scale_color_manual(values = alpha(c("black","black","black","black","black","black",0.2)) -scattercat1 -} -runPCAplotdive("divescore","PC1","PC2",1,2) -runPCAplotdive("divescore","PC1","PC2",1,2) -which(is.na(speciesPCAvalues$divescore)) -speciesPCAvalues[which(is.na(speciesPCAvalues$divescore)),] -str(speciesPCAvalues$divescore) -levels(speciesPCAvalues) -levels(speciesPCAvalues$divescore) -divecol<-c(sequential_hcl(6, palette = "Purple-Blue",rev = T)[2:6],"white") -speciesPCAvalues$divescore <- fct_explicit_na(speciesPCAvalues$divescore) -runPCAplotdive<-function(group, p1,p2,n1,n2){ -scattercat1<-ggplot(speciesPCAvalues, aes_string(x = p1, y = p2, label = "Binomial")) + -theme_classic()+ -xlab(label = paste0("p",as.character(p1),"(", -as.character(signif(d$percentexplained[n1], digits = 3)), -"%)"))+ -ylab(label = paste0("p",as.character(p2),"(", -as.character(signif(d$percentexplained[n2], digits = 3)), -"%)"))+ -#theme(legend.position = "none")+ -geom_encircle(data = speciesPCAvalues[which(is.na(speciesPCAvalues$divescore)),], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==0,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==1,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==2,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==3,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==4,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(aes_string(fill = group),s_shape=1, expand=0, color = "black", alpha = 0.7)+#s_shape = 1 and expan = 0 are convex hull -geom_point(aes_string(color = group), shape = 21, size = 3, color = "black") -#scale_color_manual(values = alpha(c("black","black","black","black","black","black",0.2)) -scattercat1 -} -runPCAplotdive("divescore","PC1","PC2",1,2) -divescore<-runPCAplotdive("divescore","PC1","PC2",1,2)+ -geom_point(aes(fill = divescore), size = 3, shape = 21,col = "black", alpha = 1)+ -scale_fill_manual(values = divecol)+ -#scale_alpha_manual(values = c(-.5,1,0,0,56)) -theme(legend.position = "right") -divescore -runPCAplotdive<-function(group, p1,p2,n1,n2){ -scattercat1<-ggplot(speciesPCAvalues, aes_string(x = p1, y = p2, label = "Binomial")) + -theme_classic()+ -xlab(label = paste0("p",as.character(p1),"(", -as.character(signif(d$percentexplained[n1], digits = 3)), -"%)"))+ -ylab(label = paste0("p",as.character(p2),"(", -as.character(signif(d$percentexplained[n2], digits = 3)), -"%)"))+ -#theme(legend.position = "none")+ -geom_encircle(data = speciesPCAvalues[which(is.na(speciesPCAvalues$divescore)),], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==0,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==1,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==2,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==3,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==4,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(aes_string(fill = group),s_shape=1, expand=0, color = "black", alpha = 0.7)+#s_shape = 1 and expan = 0 are convex hull -geom_point(aes_string(color = group), shape = 21, size = 3, color = "black")+ -scale_color_manual(values = alpha(c("black","black","black","black","black","black",0.2)) -scattercat1 -} -runPCAplotdive("divescore","PC1","PC2",1,2) -divescore<-runPCAplotdive("divescore","PC1","PC2",1,2)+ -geom_point(aes(fill = divescore), size = 3, shape = 21,col = "black", alpha = 1)+ -scale_fill_manual(values = divecol)+ -#scale_alpha_manual(values = c(-.5,1,0,0,56)) -theme(legend.position = "right") -divescore -pPCAloadings$factor<-row.names(pPCAloadings) -pPCAloadings$factor<-gsub("RES_log"," ",pPCAloadings$factor) -pPCAloadings$factor<-gsub("logHeadmassg","",pPCAloadings$factor) -loadings1<-ggplot(pPCAloadings, aes(x = PC1, y = PC2, label =factor)) + -xlab(label = "pPC1")+ -ylab(label = "pPC2")+ -#geom_point(aes(), size = 0.005) + -#geom_point(data = pPCAloadings)+ -geom_text_repel(aes()) + -theme_classic()+ -xlim(-1.5,+1)+ -ylim(-1,1)+ -geom_segment(aes(x = 0, y = 0,xend = PC1, yend = PC2), arrow = arrow(type = "closed", length = unit(0.10,"inches")))+ -theme(legend.position = "bottom") -loadings1 -ggarrange(loadings1,plunge,divescore, order,labels = c("A","B","C","D")) -levels(speciesPCAvalues$divescore) -speciesPCAvalues$divescore <- fct_explicit_na(speciesPCAvalues$divescore) -relevel(speciesPCAvalues$divescore,c( "(Missing)","0", "1", "2" , "3", "4" )) -?relevel -<-reorder(speciesPCAvalues$divescore,c( "(Missing)","0", "1", "2" , "3", "4" )) -reorder(speciesPCAvalues$divescore,c( "(Missing)","0", "1", "2" , "3", "4" )) -c( "(Missing)","0", "1", "2" , "3", "4" ) -reorder(speciesPCAvalues$divescore,c( "(Missing)","0", "1", "2" , "3", "4" )) -levels(speciesPCAvalues$divescore) -reorder(speciesPCAvalues$divescore, -c( "(Missing)","0", "1","2" , "3", "4" )) -relv<-c( "(Missing)","0", "1","2" , "3", "4" ) -speciesPCAvalues$divescore<-reorder(speciesPCAvalues$divescore,relv) -divescore -speciesPCAvalues$divescore<-relevel(speciesPCAvalues$divescore,ref = "(Missing)") -runPCAplotdive<-function(group, p1,p2,n1,n2){ -scattercat1<-ggplot(speciesPCAvalues, aes_string(x = p1, y = p2, label = "Binomial")) + -theme_classic()+ -xlab(label = paste0("p",as.character(p1),"(", -as.character(signif(d$percentexplained[n1], digits = 3)), -"%)"))+ -ylab(label = paste0("p",as.character(p2),"(", -as.character(signif(d$percentexplained[n2], digits = 3)), -"%)"))+ -#theme(legend.position = "none")+ -geom_encircle(data = speciesPCAvalues[which(is.na(speciesPCAvalues$divescore)),], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==0,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==1,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==2,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==3,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==4,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(aes_string(fill = group),s_shape=1, expand=0, color = "black", alpha = 0.7)+#s_shape = 1 and expan = 0 are convex hull -geom_point(aes_string(color = group), shape = 21, size = 3, color = "black")+ -scale_color_manual(values = alpha(c("black","black","black","black","black","black",0.2)) -scattercat1 -} -runPCAplotdive("divescore","PC1","PC2",1,2) -divescore<-runPCAplotdive("divescore","PC1","PC2",1,2)+ -geom_point(aes(fill = divescore), size = 3, shape = 21,col = "black", alpha = 1)+ -scale_fill_manual(values = divecol)+ -#scale_alpha_manual(values = c(-.5,1,0,0,56)) -theme(legend.position = "right") -divescore -divecol<-c("white",sequential_hcl(6, palette = "Purple-Blue",rev = T)[2:6]) -runPCAplotdive<-function(group, p1,p2,n1,n2){ -scattercat1<-ggplot(speciesPCAvalues, aes_string(x = p1, y = p2, label = "Binomial")) + -theme_classic()+ -xlab(label = paste0("p",as.character(p1),"(", -as.character(signif(d$percentexplained[n1], digits = 3)), -"%)"))+ -ylab(label = paste0("p",as.character(p2),"(", -as.character(signif(d$percentexplained[n2], digits = 3)), -"%)"))+ -#theme(legend.position = "none")+ -geom_encircle(data = speciesPCAvalues[which(is.na(speciesPCAvalues$divescore)),], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==0,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==1,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==2,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==3,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(data = speciesPCAvalues[speciesPCAvalues$divescore==4,], -aes_string(),s_shape=1, expand=0, alpha = 1)+ -geom_encircle(aes_string(fill = group),s_shape=1, expand=0, color = "black", alpha = 0.7)+#s_shape = 1 and expan = 0 are convex hull -geom_point(aes_string(color = group), shape = 21, size = 3, color = "black")+ -scale_color_manual(values = alpha(c("black","black","black","black","black","black",0.2)) -scattercat1 -} -runPCAplotdive("divescore","PC1","PC2",1,2) -divescore<-runPCAplotdive("divescore","PC1","PC2",1,2)+ -geom_point(aes(fill = divescore), size = 3, shape = 21,col = "black", alpha = 1)+ -scale_fill_manual(values = divecol)+ -#scale_alpha_manual(values = c(-.5,1,0,0,56)) -theme(legend.position = "right") -divescore +#add heatmap +k<-gheatmap(p,avgdf[,c("plungedistinct","divescore")], #"Category", +width = 0.25, offset = 0, +color = "black", +#colnames = FALSE, +colnames_position = "top", +colnames_angle = 0, +colnames_offset_x = 0, +colnames_offset_y = 0)+ +scale_fill_manual(values = c(divecol,"#000000", "#E69F00","#FFFFFF", "#56B4E9"), na.value = "#FFFFFF") +k +#open up the circle +l<-open_tree(k,30) +l +avgdf$Binomial2<-avgdf$Binomial +a<-ggtree(birdtreels) %<+% avgdf + ###########, layout = "circular" +geom_tiplab(aes(label = Binomial2), linesize = 0.1, offset = 30) + #circular +#geom_text(aes(label = node))+ +xlim(NA, 300) + +ylim(NA,140) +a +avgdf$IACN<-ifelse(avgdf$IAC_detail=="",0,1) +avgdf$IBPN<-ifelse(avgdf$IBP_detail=="",0,1) +names(avgdf) +sampled<-avgdf[,c("Head.mass..g.", +"TMtotalarea" , +"FPtotalarea" , +"dis_coltip_TMcentroid", +"meanTMangle" , +"Umbo_distancetoTMplane", +"CAtotalarea" , +"RWtotalarea" , +"totalEClength", +"Columella.length.mm" , +"Columella.volume.mm3" , +"Behind.TM" , +"IACN","IBPN" )]/avgdf[,c("Head.mass..g.", +"TMtotalarea" , +"FPtotalarea" , +"dis_coltip_TMcentroid", +"meanTMangle" , +"Umbo_distancetoTMplane", +"CAtotalarea" , +"RWtotalarea" , +"totalEClength", +"Columella.length.mm" , +"Columella.volume.mm3" , +"Behind.TM" , +"IACN","IBPN" )] +colnames(sampled) +colnames(sampled)<-seq(1:length(colnames(sampled))) +#current sampling +k2<-gheatmap(a,sampled, +width = 0.4, offset = 0, +color = "black", +low = "white", high = "blue", +colnames = T, +colnames_position = "top", +colnames_angle = 0, +colnames_offset_x = 0, +colnames_offset_y = 5) +k2 +ggtree() +?ggheatmap +?gheatmap/ +?gheatmap +a<-ggtree(birdtreels) %<+% avgdf + ###########, layout = "circular" +geom_tiplab(aes(label = Binomial2), linesize = 0.1, offset = 30) + #circular +#geom_text(aes(label = node))+ +xlim(NA, 300) + +ylim(NA,140) +a +a<-ggtree(birdtreels) %<+% avgdf + ###########, layout = "circular" +geom_tiplab(aes(label = Binomial2), linesize = 0.1, offset = 30) + #circular +#geom_text(aes(label = node))+ +xlim(NA, 200) + +ylim(NA,140) +a +a<-ggtree(birdtreels) %<+% avgdf + ###########, layout = "circular" +geom_tiplab(aes(label = Binomial2), linesize = 0.1, offset = 30) + #circular +#geom_text(aes(label = node))+ +xlim(NA, 100) + +ylim(NA,140) +a +a<-ggtree(birdtreels) %<+% avgdf + ###########, layout = "circular" +geom_tiplab(aes(label = Binomial2), linesize = 0.1, offset = 30) + #circular +#geom_text(aes(label = node))+ +xlim(NA, 500) + +ylim(NA,140) +a +a<-ggtree(birdtreels) %<+% avgdf + ###########, layout = "circular" +geom_tiplab(aes(label = Binomial2), linesize = 0.1, offset = 30) + #circular +#geom_text(aes(label = node))+ +xlim(NA, 150) + +ylim(NA,140) +a +a<-ggtree(birdtreels) %<+% avgdf + ###########, layout = "circular" +geom_tiplab(aes(label = Binomial2), linesize = 0.1, offset = 30) + #circular +#geom_text(aes(label = node))+ +xlim(NA, 160) + +ylim(NA,140) +a +avgdf$IACN<-ifelse(avgdf$IAC_detail=="",0,1) +avgdf$IBPN<-ifelse(avgdf$IBP_detail=="",0,1) +names(avgdf) +sampled<-avgdf[,c("Head.mass..g.", +"TMtotalarea" , +"FPtotalarea" , +"dis_coltip_TMcentroid", +"meanTMangle" , +"Umbo_distancetoTMplane", +"CAtotalarea" , +"RWtotalarea" , +"totalEClength", +"Columella.length.mm" , +"Columella.volume.mm3" , +"Behind.TM" , +"IACN","IBPN" )]/avgdf[,c("Head.mass..g.", +"TMtotalarea" , +"FPtotalarea" , +"dis_coltip_TMcentroid", +"meanTMangle" , +"Umbo_distancetoTMplane", +"CAtotalarea" , +"RWtotalarea" , +"totalEClength", +"Columella.length.mm" , +"Columella.volume.mm3" , +"Behind.TM" , +"IACN","IBPN" )] +colnames(sampled) +colnames(sampled)<-seq(1:length(colnames(sampled))) +#current sampling +k2<-gheatmap(a,sampled, +width = 0.4, offset = 0, +color = "black", +low = "white", high = "blue", +colnames = T, +colnames_position = "top", +colnames_angle = 0, +colnames_offset_x = 0, +colnames_offset_y = 5) +k2 +a<-ggtree(birdtreels) %<+% avgdf + ###########, layout = "circular" +geom_tiplab(aes(label = Binomial2), linesize = 0.1, offset = 30) + #circular +#geom_text(aes(label = node))+ +xlim(NA, 160) + +ylim(NA,140) +a +a<-ggtree(birdtreels) %<+% avgdf + ###########, layout = "circular" +geom_tiplab(aes(label = Binomial2), linesize = 0.1, offset = 30) + #circular +#geom_text(aes(label = node))+ +xlim(NA, 160) + +ylim(NA,120) +a +a<-ggtree(birdtreels) %<+% avgdf + ###########, layout = "circular" +geom_tiplab(aes(label = Binomial2), linesize = 0.1, offset = 30) + #circular +#geom_text(aes(label = node))+ +xlim(NA, 160) + +ylim(NA,130) +a +avgdf$IACN<-ifelse(avgdf$IAC_detail=="",0,1) +avgdf$IBPN<-ifelse(avgdf$IBP_detail=="",0,1) +names(avgdf) +sampled<-avgdf[,c("Head.mass..g.", +"TMtotalarea" , +"FPtotalarea" , +"dis_coltip_TMcentroid", +"meanTMangle" , +"Umbo_distancetoTMplane", +"CAtotalarea" , +"RWtotalarea" , +"totalEClength", +"Columella.length.mm" , +"Columella.volume.mm3" , +"Behind.TM" , +"IACN","IBPN" )]/avgdf[,c("Head.mass..g.", +"TMtotalarea" , +"FPtotalarea" , +"dis_coltip_TMcentroid", +"meanTMangle" , +"Umbo_distancetoTMplane", +"CAtotalarea" , +"RWtotalarea" , +"totalEClength", +"Columella.length.mm" , +"Columella.volume.mm3" , +"Behind.TM" , +"IACN","IBPN" )] +colnames(sampled) +colnames(sampled)<-seq(1:length(colnames(sampled))) +#current sampling +k2<-gheatmap(a,sampled, +width = 0.4, offset = 0, +color = "black", +low = "white", high = "blue", +colnames = T, +colnames_position = "top", +colnames_angle = 0, +colnames_offset_x = 0, +colnames_offset_y = 5) +k2 +a<-ggtree(birdtreels) %<+% avgdf + ###########, layout = "circular" +geom_tiplab(aes(label = Binomial2), linesize = 0.1, offset = 30) + #circular +#geom_text(aes(label = node))+ +xlim(NA, 160) + +ylim(NA,135) +a +avgdf$IACN<-ifelse(avgdf$IAC_detail=="",0,1) +avgdf$IBPN<-ifelse(avgdf$IBP_detail=="",0,1) +names(avgdf) +sampled<-avgdf[,c("Head.mass..g.", +"TMtotalarea" , +"FPtotalarea" , +"dis_coltip_TMcentroid", +"meanTMangle" , +"Umbo_distancetoTMplane", +"CAtotalarea" , +"RWtotalarea" , +"totalEClength", +"Columella.length.mm" , +"Columella.volume.mm3" , +"Behind.TM" , +"IACN","IBPN" )]/avgdf[,c("Head.mass..g.", +"TMtotalarea" , +"FPtotalarea" , +"dis_coltip_TMcentroid", +"meanTMangle" , +"Umbo_distancetoTMplane", +"CAtotalarea" , +"RWtotalarea" , +"totalEClength", +"Columella.length.mm" , +"Columella.volume.mm3" , +"Behind.TM" , +"IACN","IBPN" )] +colnames(sampled) +colnames(sampled)<-seq(1:length(colnames(sampled))) +#current sampling +k2<-gheatmap(a,sampled, +width = 0.4, offset = 0, +color = "black", +low = "white", high = "blue", +colnames = T, +colnames_position = "top", +colnames_angle = 0, +colnames_offset_x = 0, +colnames_offset_y = 5) +k2 +library(ggplot2) +theme_set(theme_classic()) +arearatio<-read.csv(file.choose()) +names(arearatio) +ddprocell<-subset(divedepth, order =='Procellariformes') +petrels<-subset(ddprocell, family=='Procellariidae') +arearatio$fullname<- paste(arearatio$Bird.group, arearatio$Source, sep = " - ") +# Plot +g<- ggplot(arearatio, aes(x = Bird.group, y = Area.ratio, fill = Bird.group))+ +geom_boxplot(aes(Bird.group, Area.ratio)) + +theme(axis.text.x=element_text(angle = 30, hjust = 1), +legend.position = "none")+ +# axis.title=element_text(size=14,face="bold"))+ +scale_fill_manual(values = c("lightblue","lightblue","lightgreen","lightgreen","lightgreen")) + +geom_point(aes(Bird.group)) +g diff --git a/data Nov 17_.csv b/Eardata.csv similarity index 100% rename from data Nov 17_.csv rename to Eardata.csv diff --git a/Load phylogeny and make CDO.R b/Load phylogeny and make CDO.R index 1895a7f..de97d1e 100644 --- a/Load phylogeny and make CDO.R +++ b/Load phylogeny and make CDO.R @@ -8,7 +8,7 @@ birdtreels<-birdtree distinctdf<-distinct(df, Binomial, .keep_all = TRUE) distinctdforder<-arrange(distinctdf,Binomial) -#put current correct binomials on the tree, updated from t he birdtree +#put current correct binomials on the tree, updated from the birdtree new<-cbind.data.frame(birdtreels$tip.label,distinctdf$Binomial,gsub(" ","_",distinctdf$Birdtree)) colnames(new)<-c("tiplabel","binomial","birdtree") #View(new) @@ -20,7 +20,7 @@ str(birdtreels$tip.label) birdtreels$tip.label<-new$binomialordered birdtreels$tip.label<-as.character(birdtreels$tip.label) -#make comparative data object for caper +#make comparative data object for caper's pgls function birdCDO<-comparative.data(phy = birdtreels,data = distinctdf, #avgdf[avgdf$Category!="Terrestrial",] names.col = Binomial, vcv = TRUE, na.omit = FALSE, diff --git a/Plot PCA.R b/Plot PCA.R index 1a5dc32..30aaf9f 100644 --- a/Plot PCA.R +++ b/Plot PCA.R @@ -10,13 +10,42 @@ d<-as.data.frame(diag(pPCA$Eval)/sum(pPCA$Eval)*100) d$PC1<-row.names(d) d$percentexplained<-d$`diag(pPCA$Eval)/sum(pPCA$Eval) * 100` + +#eigenvalue plot p<-ggplot(d, aes(x = reorder(PC1,-percentexplained), y = percentexplained))+ - geom_bar(stat = "identity")+ - theme_classic()+ + geom_bar(stat = "identity", col = "white")+ + theme_minimal()+ xlab("Principal component")+ ylab("Percent variance \n explained") p +#remove lavels for inset plot +p<-ggplot(d[1:6,], aes(x = reorder(PC1,-percentexplained), y = percentexplained))+ + geom_bar(stat = "identity")+ + theme_classic()+ + xlab("")+ + ylab("")+ + theme(axis.text.x = element_blank(), + axis.text.x = element_text(angle = 90), + panel.background = element_blank(), + panel.grid = element_blank()) +p + + +themejz<-function (base_size = 11, base_family = "", base_line_size = base_size/22, + base_rect_size = base_size/22) +{ + theme_grey(base_size = base_size, base_family = base_family, + base_line_size = base_line_size, base_rect_size = base_rect_size) %+replace% + theme(panel.background = element_blank(), panel.border = element_rect(), panel.grid = element_line(), + panel.grid.minor = element_blank(), + strip.background = element_rect(fill = "grey85", + colour = "grey20"), legend.key = element_rect(fill = "white", + colour = NA), complete = TRUE) +} + + + ################Make a grouping of orders that includes species with one underwater-pursuit species##################### speciesPCAvalues$Order<-as.character(speciesPCAvalues$Order) speciesPCAvalues$Order2<-NA @@ -238,8 +267,8 @@ pPCAloadings$factor<-gsub("logHeadmassg","",pPCAloadings$factor) loadings1<-ggplot(pPCAloadings, aes(x = PC1, y = PC2, label =factor)) + - xlab(label = "pPC1")+ - ylab(label = "pPC2")+ + xlab(label = "pPC1(43.4%)")+ + ylab(label = "pPC2(14%)")+ #geom_point(aes(), size = 0.005) + #geom_point(data = pPCAloadings)+ @@ -251,8 +280,13 @@ loadings1<-ggplot(pPCAloadings, aes(x = PC1, y = PC2, label =factor)) + theme(legend.position = "bottom") loadings1 +loadings2<-loadings1 + annotation_custom(ggplotGrob(p), + xmin = 0.25, xmax = 1, + ymin = 0, ymax = 1) +loadings2 -ggarrange(loadings1,plunge,divescore, order,labels = c("A","B","C","D")) +ggarrange(loadings2,plunge,divescore, order, + labels = c("A","B","C","D")) ggsave("D:/00_Manuscripts/0Avian aquatic hearing project/___Oct 1 version/PCAOct 4_noair.pdf",width = 10, height = 10) ggsave("D:/00_Manuscripts/0Avian aquatic hearing project/___Oct 1 version/PCAOct 4_withair.pdf",width = 10, height = 10) diff --git a/README.md b/README.md index 10a2910..240c7b2 100644 --- a/README.md +++ b/README.md @@ -1,9 +1,9 @@ # A_T-stats - +[![DOI](https://zenodo.org/badge/313591727.svg)](https://zenodo.org/badge/latestdoi/313591727) ## 1. Data files -There are 3 data files: ->Main dataframe: "data Nov 17_.csv" +There are three data files: +>Main dataframe: "Eardata.csv" Phylogeny: "JZ tree Prum merged hackett.tree" Dive depth data: "Depth list.csv" @@ -50,4 +50,6 @@ This runs the analysis comparing. There is a separate script for the interaural ## 7. Other plotting scripts > "plot ecology on circular phylogeny.R" - plot ecological groupings on circular phylogeny (Fig 2), and sampling by phylogeny (Fig S2) -"Residualphylogplots.R - plotting residual plot (Fig 3) +"Residualphylogplots_a panel top.R - plotting residual plot (Fig 3) + + diff --git a/Residualphylogplots.R b/Residualphylogplots.R deleted file mode 100644 index 08af812..0000000 --- a/Residualphylogplots.R +++ /dev/null @@ -1,308 +0,0 @@ -#Plotting adapted from: -#https://thackl.github.io/ggtree-composite-plots -library(patchwork) -library(viridis) - -#color palettes -# The palette with grey: -cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7") -# The palette with black: -cbbPalette <- c( "#FFFFFF","#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7") -barplot(1:5, col=cbbPalette) - -#subset the df with resids dataframe to include only the residuals -#run the PCA scripts to create the dfwith resids data frame -names(dfwithresids) - -#input column numbers that have residuals -residseach<-dfwithresids[,c(37:48)] -names(residseach) - -#add -residseach$order<-dfwithresids$Order -residseach$Binomial<-dfwithresids$Binomial -residseach$Category<-dfwithresids$Category -residseach$divescore<-dfwithresids$divescore -residseach$plungedistinct<-dfwithresids$plungedistinct -residseach$waterbirds<-dfwithresids$waterbirds -residseach$superorder<-dfwithresids$superorder -residseach$IAC<-dfwithresids$IAC_detail -residseach$IBP<-dfwithresids$IBP_detail -residseach$PC1<-dfwithresids$PC1 - - -residseach$plungedistinct<-relevel(as.factor(residseach$plungedistinct),"Terrestrial") -#residseach$Category<-relevel(as.factor(residseach$Category),"Terrestrial") - -#put in long format -longdfplotting<-gather(residseach,key = "earmeasuresresid", value = "earmeasureval", - -c(order,Binomial,Category,divescore,plungedistinct,IAC,IBP))#waterbirds,superorder, -longdfplotting$earmeasuresresid<-as.factor(longdfplotting$earmeasuresresid) -longdfplotting$earmeasureval<-as.numeric(longdfplotting$earmeasureval) -#longdfplotting$label<-longdfplotting$order -#longdfplotting$label<-longdfplotting$waterbirds -#longdfplotting$label<-longdfplotting$superorder - -longdfplotting$label<-longdfplotting$order -longdfplotting$earmeasuresresid<-as.character(longdfplotting$earmeasuresresid) - -longdfplotting$Terr<-longdfplotting$Category -longdfplotting$Terr<-ifelse(longdfplotting$Category=="Terrestrial","Terrestrial",NA) - -#make list of ear measures to plot -yvarnames<-c( - "RES_logTMtotalarealogHeadmassg" , - "RES_logFPtotalarealogHeadmassg" , - "RES_logarearatiologHeadmassg" , - "RES_logdiscoltipTMcentroidlogHeadmassg" , - "RES_logUmbodistancetoTMplanelogHeadmassg", - "RES_logmeanTManglelogHeadmassg" , - "RES_logtotalEClengthlogHeadmassg" , - "RES_logRWtotalarealogHeadmassg" , - "RES_logCAtotalarealogHeadmassg" , - "RES_logBehindTMlogHeadmassg" , - "RES_logColumellalengthmmlogHeadmassg" , - "RES_logColumellavolumemm3logHeadmassg" ,"PC1") - -#create functions for plotting -source("Tblog.R") -#add medians -source("add median.R") - -#############make the order-level cladogram for aligning the residuals########### -source("Order_level_cladogram.R") -### -ggtree(orderPhy)+ - geom_text(aes(label = node)) - -gg_tr <- ggtree(orderPhy, branch.length = "none") + - geom_tiplab(align=TRUE) + - #scale_x_continuous(expand=expand_scale(0.2)) + # make more room for the labels - scale_y_tree()+ - xlim(0,40)+ - ylim(-3.5,26.5) - #geom_text(aes(x = 13,y = -0,label = "Terrestrial"))+ - #geom_text(aes(x = 13,y = -1,label = "Surface-foraging"))+ - #geom_text(aes(x = 13,y = -2,label = "Plunge-diving"))+ - #geom_text(aes(x = 13,y = -3,label = "Underwater-pursuit")) -gg_tr - -#reversed phylogenetic tree -gg_tr_rev <- ggtree(orderPhy, branch.length = "none", col = "white") + - geom_tiplab(align=TRUE) + - scale_x_continuous(expand=expand_scale(0.2)) + # make more room for the labels - scale_y_tree()+ - xlim(10,20)+ - ylim(-3.5,26.5)+ - scale_x_reverse() - #geom_text(aes(x = 11,y = -0,label = "Terrestrial"))+ - #geom_text(aes(x = 11,y = -1,label = "Surface-foraging"))+ - #geom_text(aes(x = 11,y = -2,label = "Plunge-diving"))+ - #geom_text(aes(x = 11,y = -3,label = "Underwater-pursuit")) -gg_tr_rev - -###############scatterplot function########### -gg_plungedistinct<-function(index2, letter, box = "yes"){ - ggtreeplot(gg_tr, subset(longdfplotting, - longdfplotting$earmeasuresresid==yvarnames[index2]), aes(y=earmeasureval), flip=TRUE) +{ - - if(box == "yes") geom_rect(aes(xmin = -3.5, xmax = 26.5, ymin = Inf, ymax = -Inf), fill = "lightgrey", alpha = 0.5) else geom_rect(aes(xmin = 0.5, xmax = -4.5, ymin = Inf, ymax = -Inf), fill = "white", alpha = 0.001) - } + - #geom_vline(xintercept = -4:30, col = "black")+ - #geom_rect(aes(xmin = -3.5, xmax = 26.5, ymin = -Inf, ymax = -Inf), col = "black", fill = "white", alpha = 0.1)+ - #geom_rect(aes(xmin = -3.5, xmax = 26.5, ymin = Inf, ymax = Inf), col = "black", fill = "white", alpha = 0.1)+ - #geom_line(aes(xmin = -3.5, xmax = 26.5, ymin = 0, ymax = 0), col = "black", size = 1)+ - geom_segment(aes(x = -3.5,xend = 26.5, y = 0, yend = 0), col = "black")+ - geom_point(aes(fill = plungedistinct), size = 2,shape = 21, col = "black")+ - scale_fill_manual(values = cbbPalette)+ - geom_vline(xintercept = 0.5,col = "black")+ - coord_flip() + no_y_axis()+ - ylab("")+ - xlim(-3.5,26.5)+ - theme_classic() + - theme(axis.line.y = element_blank(), - axis.title.y = element_blank(), - axis.text.y = element_blank(), - axis.ticks.y = element_blank(), - axis.text.x = element_text(angle = 90, colour = "black"), - legend.position = "none", - plot.margin = margin(t = 0, r = 0, b = 0, l = 0, unit = "pt"))+ - geom_boxplot(data = subset(longdfplotting, - longdfplotting$earmeasuresresid==yvarnames[index2]& - longdfplotting$plungedistinct=="Terrestrial"), - aes(x = 0, y = earmeasureval),fill= "white", col = "black", outlier.size = 2, - outlier.colour = "black", outlier.fill = "white", - outlier.shape = 21)+ - geom_text(aes(x=Inf,y=-Inf,vjust = 1, - hjust = -1,label=letter)) - } - -gg_plungedistinct(1, letter = "K",box = "no") -addbxplt(1,1,"a") - -gg_plungedistinct(1, letter = "K", box = "yes") - -############function to add boxplots at bottom################ -addbxplt<-function(j,index2,letter,box = "yes"){ - d<-gg_plungedistinct(index2,letter,box)+ - geom_boxplot(data = subset(longdfplotting, - longdfplotting$earmeasuresresid==yvarnames[j]& - longdfplotting$plungedistinct=="Plunging"), - aes(x = -2, y = earmeasureval),fill= "black",col = "black", outlier.size = 2, - outlier.fill = "#000000", outlier.shape = 21, outlier.color = "black")+ - geom_boxplot(data = subset(longdfplotting, - longdfplotting$earmeasuresresid==yvarnames[j]& - longdfplotting$plungedistinct=="Underwater pursuit"), - aes(x = -3, y = earmeasureval),fill ="#56B4E9",col = "black", outlier.size = 2, - outlier.fill = "#56B4E9", outlier.shape = 21, outlier.color = "black")+ - - geom_boxplot(data = subset(longdfplotting, - longdfplotting$earmeasuresresid==yvarnames[j]& - longdfplotting$plungedistinct=="Surface"), - mapping = aes(x = -1, y = earmeasureval), fill = "#E69F00",col = "black", - outlier.size = 2, outlier.fill = "#E69F00", outlier.shape = 21, outlier.color = "black") - - - d -} - - - -#####################plots for interaural canal an interbullar passage#######333 -#data for plotting by order -summ2<-avgdf %>% group_by(Order,IBP_detail) %>% count(na.omit = T) -summ2$IBP_detail<-ifelse(summ2$IBP_detail=="Pneumaticity present"| - summ2$IBP_detail=="Pneumaticity absent"| - summ2$IBP_detail=="Y",summ2$IBP_detail,NA) - -summ2$label<-summ2$Order -summ2$number<-summ2$n -summ2<-as.data.frame(summ2) -names(summ2) -summb<-summ2[which(!is.na(summ2$IBP_detail)),] - -#data for plotting by ecology -summpl<-avgdf %>% group_by(IBP_detail,plungedistinct) %>% count(na.omit = T) -summpl$IBP_detail<-ifelse(summpl$IBP_detail=="Pneumaticity present"| - summpl$IBP_detail=="Pneumaticity absent"| - summpl$IBP_detail=="Y",summpl$IBP_detail,NA) - -summpl$number<-summpl$n -summpl<-as.data.frame(summpl) -names(summpl) -summpl_<-summpl[which(!is.na(summpl$IBP_detail)),] - - -IBP<-function(letter){ - d<-ggtreeplot(gg_tr, summb, aes(y = number), flip=TRUE) + - geom_bar(aes(fill = as.factor(IBP_detail)), position = "fill", - color = "black",stat = "identity")+ - scale_fill_manual(values = c("black","grey","white"))+ - #no_legend()+ - xlim(-3.5,26.5)+ - coord_flip() + no_y_axis()+ - theme(axis.text.x = element_text(angle = 90))+ - ylab("")+ - theme_classic() + - theme(axis.line.y = element_blank(), - axis.title.y = element_blank(), - axis.text.y = element_blank(), - axis.ticks.y = element_blank(), - axis.text.x = element_text(angle = 90, colour = "black"), - legend.position = "none", - plot.margin = margin(t = 0, r = 0, b = 0, l = 0, unit = "pt"))+ - geom_text(aes(x=Inf,y=-Inf,vjust = 1, - hjust = -1,label=letter)) - - d -} -IBP("d") - - -IBPfull<-function(letter){ - IBP(letter)+ - geom_col(data = summpl_[summpl_$plungedistinct=="Underwater pursuit",], - aes(x = -3, fill = IBP_detail),col = "black", position = "fill")+ - geom_col(data = summpl_[summpl_$plungedistinct=="Terrestrial",], - aes(x = 0, fill = IBP_detail),col = "black", position = "fill")+ - geom_col(data = summpl_[summpl_$plungedistinct=="Plunging",], - aes(x = -2, fill = IBP_detail),col = "black", position = "fill")+ - geom_col(data = summpl_[summpl_$plungedistinct=="Surface",], - aes(x = -1, fill = IBP_detail),col = "black", position = "fill") -} -IBPfull("l") - -### -summ<-avgdf %>% group_by(Order,IAC_detail) %>% count(na.omit = T) -summ$IAC_detail<-ifelse(summ$IAC_detail=="Pneumaticity present"| - summ$IAC_detail=="Pneumaticity absent"| - summ$IAC_detail=="Y",summ$IAC_detail,NA) - - -summ$label<-summ$Order -summ$number<-summ$n -summ<-as.data.frame(summ) -summa<-summ[which(!is.na(summ$IAC_detail)),] - - -summpl2<-avgdf %>% group_by(IAC_detail,plungedistinct) %>% count(na.omit = T) -summpl2$IAC_detail<-ifelse(summpl2$IAC_detail=="Pneumaticity present"| - summpl2$IAC_detail=="Pneumaticity absent"| - summpl2$IAC_detail=="Y",summpl2$IAC_detail,NA) - -summpl2$number<-summpl2$n -summpl2<-as.data.frame(summpl2) -names(summpl2) -summpl2_<-summpl2[which(!is.na(summpl2$IAC_detail)),] - - -#IAC -IAC<-function(letter){ - d<-ggtreeplot(gg_tr, summa, aes(y = number), flip=TRUE) + - geom_bar(aes(fill = as.factor(IAC_detail)), position = "fill", - color = "black",stat = "identity")+ - scale_fill_manual(values = c("black","grey","white"))+ - #no_legend()+ - xlim(-3.5,26.5)+ - coord_flip() + no_y_axis()+ - theme(axis.text.x = element_text(angle = 90))+ - ylab("")+ - theme_classic() + - theme(axis.line.y = element_blank(), - axis.title.y = element_blank(), - axis.text.y = element_blank(), - axis.ticks.y = element_blank(), - axis.text.x = element_text(angle = 90, colour = "black"), - legend.position = "none", - plot.margin = margin(t = 0, r = 0, b = 0, l = 0, unit = "pt"))+ - geom_text(aes(x=Inf,y=-Inf,vjust = 1, - hjust = -1,label=letter)) - - d -} -IAC("d") - - -IACfull<-function(letter){ - - IAC(letter)+ - geom_col(data = summpl2_[summpl2_$plungedistinct=="Underwater pursuit",], - aes(x = -3, fill = IAC_detail),col = "black", position = "fill")+ - geom_col(data = summpl2_[summpl2_$plungedistinct=="Terrestrial",], - aes(x = 0, fill = IAC_detail),col = "black", position = "fill")+ - geom_col(data = summpl2_[summpl2_$plungedistinct=="Plunging",], - aes(x = -2, fill = IAC_detail),col = "black", position = "fill")+ - geom_col(data = summpl2_[summpl2_$plungedistinct=="Surface",], - aes(x = -1, fill = IAC_detail),col = "black", position = "fill") -} -IACfull("m") - - -###################Plot all together####################### -gg_tr|addbxplt(1,1,"a")+addbxplt(2,2,"b")| - addbxplt(3,3,"c")+addbxplt(4,4,"d")|# - addbxplt(5,5,"e")+addbxplt(6,6,"f")|#umbo height and TM angle - addbxplt(8,8,"g")+addbxplt(9,9,"h", box = "no")|#ESlength & RW - addbxplt(7,7,"i")+addbxplt(11,11,"j")|#CA and collenght - addbxplt(12,12,"k", box = "no")+addbxplt(10,10,"l")|#colvol and air - IACfull("m")+IBPfull("n")|gg_tr_rev diff --git a/Residualphylogplots_a panel top.R b/Residualphylogplots_a panel top.R index 3ab8b65..8ec6508 100644 --- a/Residualphylogplots_a panel top.R +++ b/Residualphylogplots_a panel top.R @@ -67,7 +67,7 @@ yvarnames<-c( #create functions for plotting source("Tblog.R") #add medians -source("add median.R") +#source("add median.R") #############make the order-level cladogram for aligning the residuals########### source("Order_level_cladogram.R") @@ -108,7 +108,7 @@ gg_plungedistinct<-function(index2, letter, box = "yes"){ ggtreeplot(gg_tr, subset(longdfplotting, longdfplotting$earmeasuresresid==yvarnames[index2]), aes(y=earmeasureval), flip=TRUE) +{ - if(box == "yes") annotate("text",x = Inf, y = Inf, label = "*", hjust = 2, vjust = 1.25) + if(box == "yes") annotate("text",x = Inf, y = Inf, label = "*", hjust = 2, vjust = 1.25, size = 5) else annotate("text",x = Inf, y = Inf, label = "", hjust = 2, vjust = 1.25) } + #geom_label(aes(x=Inf,y=-Inf,vjust = 1, @@ -120,7 +120,7 @@ gg_plungedistinct<-function(index2, letter, box = "yes"){ #geom_line(aes(xmin = -3.5, xmax = 26.5, ymin = 0, ymax = 0), col = "black", size = 1)+ geom_segment(aes(x = 0.5,xend = 38, y = 0, yend = 0), col = "black")+ #geom_segment(aes(x = 28,xend = 36, y = 0, yend = 0), col = "black")+ - geom_point(aes(fill = plungedistinct), size = 2,shape = 21, col = "black")+ + geom_point(aes(fill = plungedistinct), size = 3,shape = 21, col = "black")+ scale_fill_manual(values = cbbPalette)+ #geom_vline(xintercept = 27,col = "black")+ @@ -139,10 +139,10 @@ gg_plungedistinct<-function(index2, letter, box = "yes"){ geom_boxplot(data = subset(longdfplotting, longdfplotting$earmeasuresresid==yvarnames[index2]& longdfplotting$plungedistinct=="Terrestrial"), - aes(x = 36, y = earmeasureval),fill= "white", col = "black", outlier.size = 2, + aes(x = 36, y = earmeasureval),fill= "white", col = "black", outlier.size = 3, outlier.colour = "black", width = 2.25, outlier.fill = "white", outlier.shape = 21)+ - annotate("text",x = Inf, y = -Inf, label = letter,hjust = -1,vjust = 1.25) + annotate("text",x = Inf, y = -Inf, label = letter,hjust = -1,vjust = 1.25, size = 5) } @@ -167,20 +167,20 @@ addbxplt<-function(j,index2,letter,box = "yes"){ geom_boxplot(data = subset(longdfplotting, longdfplotting$earmeasuresresid==yvarnames[j]& longdfplotting$plungedistinct=="Plunging"), - aes(x = 31, y = earmeasureval),fill= "black",col = "black", outlier.size = 2, + aes(x = 31, y = earmeasureval),fill= "black",col = "black", outlier.size = 3, outlier.fill = "#000000", outlier.shape = 21, width = 2.25, outlier.color = "black")+ geom_boxplot(data = subset(longdfplotting, longdfplotting$earmeasuresresid==yvarnames[j]& longdfplotting$plungedistinct=="Underwater pursuit"), - aes(x = 28.5, y = earmeasureval),fill ="#56B4E9",col = "black", outlier.size = 2, + aes(x = 28.5, y = earmeasureval),fill ="#56B4E9",col = "black", outlier.size = 3, outlier.fill = "#56B4E9", outlier.shape = 21, width = 2.25, outlier.color = "black")+ geom_boxplot(data = subset(longdfplotting, longdfplotting$earmeasuresresid==yvarnames[j]& longdfplotting$plungedistinct=="Surface"), mapping = aes(x = 33.5, y = earmeasureval), fill = "#E69F00",col = "black", - outlier.size = 2, outlier.fill = "#E69F00", width = 2.25, outlier.shape = 21, outlier.color = "black") + outlier.size = 3, outlier.fill = "#E69F00", width = 2.25, outlier.shape = 21, outlier.color = "black") d @@ -188,7 +188,7 @@ addbxplt<-function(j,index2,letter,box = "yes"){ -#####################plots for interaural canal an interbullar passage#######333 +#####################plots for interaural canal an interbullar passage####### #data for plotting by order summ2<-avgdf %>% group_by(Order,IBP_detail) %>% count(na.omit = T) summ2$IBP_detail<-ifelse(summ2$IBP_detail=="Pneumaticity present"| @@ -296,7 +296,7 @@ IAC<-function(letter){ axis.text.x = element_text(angle = 90, colour = "black"), legend.position = "none", plot.margin = margin(t = 0, r = 0, b = 0, l = 0, unit = "pt")) + - annotate("text",x = Inf, y = -Inf, label = letter,hjust = -1,vjust = 1.25) + annotate("text",x = Inf, y = -Inf, label = letter,hjust = -1,vjust = 1.25, size = 5) # geom_text(aes(x=Inf,y=-Inf,vjust = 1, # hjust = -1,label=letter)) @@ -316,7 +316,7 @@ IACfull<-function(letter){ aes(x = 31, fill = IAC_detail),col = "black", position = "fill",width = 2.25)+ geom_col(data = summpl2_[summpl2_$plungedistinct=="Surface",], aes(x = 33.5, fill = IAC_detail),col = "black", position = "fill",width = 2.25)+ - annotate("text",x = Inf, y = Inf, label = "*", hjust = 2, vjust = 1.25) + annotate("text",x = Inf, y = Inf, label = "*", hjust = 2, vjust = 1.25, size = 5) } IACfull("m") @@ -325,10 +325,10 @@ IACfull("m") ###################Plot all together####################### (gg_tr|addbxplt(1,1,"i")+addbxplt(2,2,"ii")| addbxplt(3,3,"iii")+addbxplt(4,4,"iv")|# - addbxplt(5,5,"v")+addbxplt(6,6,"vi")|gg_tr_rev) -/ -/#umbo height and TM angle + addbxplt(5,5,"v")+addbxplt(6,6,"vi")) + +#umbo height and TM angle (gg_tr|addbxplt(8,8,"vii")+addbxplt(9,9,"viii", box = "no")|#ESlength & RW addbxplt(7,7,"ix")+addbxplt(11,11,"x")|#CA and collenght addbxplt(12,12,"xi", box = "no")+addbxplt(10,10,"xii")|#colvol and air - IACfull("xiii")+IBPfull("xiv")|gg_tr_rev) + IACfull("xiii")+IBPfull("xiv")) diff --git a/Run PGLS and AIC.R b/Run PGLS and AIC.R index be013a2..f1b8f9d 100644 --- a/Run PGLS and AIC.R +++ b/Run PGLS and AIC.R @@ -1,7 +1,7 @@ library(flextable) library(officer) -options(scipen=999, digits = 3)# take out of scientific notation and allow 3 significant digist +options(scipen=999, digits = 3)# take out of scientific notation and allow 3 significant digits ########################################pgls function#################################### pgls_models<-function(i){ @@ -13,15 +13,15 @@ pgls_models<-function(i){ ###########run different pgls models and save as a list of tables######## ############Each table being the stats results for a different ear measure -source("pgls_HM.R") # -source("pgls_HM_plus_ecol.R") # -source("pgls_HM_times_ecol.R") # +source("pgls_HM.R") # model with only head mass +source("pgls_HM_plus_ecol.R") #model with ecological groupings +source("pgls_HM_times_ecol.R") #model with ecological groupings and groupin x head mass interaction term #for aquatic-only analysis, go back and create a new data object lacking Terrestrial species, and also include models with dive score: -source("pgls_HM_plus_divescore.R") # -source("pgls_HM_times_divescore.R") # +source("pgls_HM_plus_divescore.R") #model with dive score +source("pgls_HM_times_divescore.R") #model with divescore and head mass x dive score interaction -####for each ear measure, combine together the different models single dataframe +####for each ear measure, combine together the different model outputs into single dataframe grouped_eachmeasure<-list() for (i in 1:length(tbllist_HM)){ grouped_eachmeasure[[i]]<-rbind(tbllist_HM[[i]],tbllist_HM_plus_ecol[[i]], @@ -38,7 +38,7 @@ for (i in 1:length(tbllist_HM)){ ##########The following remaining code does two things:############### #(1) get the AIC values to select the best model -#(2)report the details of the best models (models with change AIC) +#(2) report the details of the best models (models with change AIC) ############(1)Get the summary stats for each model############ #Compute the AIC, relative likelihood and AIC weight @@ -71,7 +71,6 @@ all<-all[,c(ncol(all),(1:(ncol(all)-1)))] all$Measurementgroup[which(duplicated(all$Measurementgroup))]<-"" all$F_<-paste0(all$Fstat,"(",all$Fstat_numdf,",",all$Fstat_dendf,")") - #visualize the table better using the flextable package flexall<-flextable(all) %>% add_header_lines( values = "Table X. Models for selection") %>% diff --git a/Run phyPCA.R b/Run phyPCA.R index 3ac1da9..f3ce6a5 100644 --- a/Run phyPCA.R +++ b/Run phyPCA.R @@ -11,7 +11,6 @@ names(dfwithresids) #isolate columns with residuals to be used in PCA PCAset<-dfwithresids[,37:48]# data to be used in PCA row.names(PCAset)<-dfwithresids$Binomial#species names need to be in row names for use of phy.PCA function - names(PCAset) #remove air volume from PCA if desired (air volume has more missing values than other measures) diff --git a/SW_HM_.R b/SW_HM_.R index e1a2029..ba7edc2 100644 --- a/SW_HM_.R +++ b/SW_HM_.R @@ -4,14 +4,12 @@ pglsfit<-pgls(log(distinctdf$Head.mass..g.)~log(distinctdf$Skull.width..mm.), da bounds = list(lambda=c(0.00001,1))) print(summary(pglsfit)) - - -#ones without headmass +#Show species with no head massdata available length(df$Binomial[which(is.na(df$Head.mass..g.))]) #nrow(df[which(is.na(df$Head.mass..g.)),]) #df$Binomial[which(is.na(df$Skull.width..mm.))] -#add imputed head mass that are NA from the specimen's skull width +#add imputed head mass based on the specimen's skull width and the pgls regression df$Head.mass..g.[which(is.na(df$Head.mass..g.))]<- exp(log(df$Skull.width..mm.[which(is.na(df$Head.mass..g.))])*pglsfit$model$coef[2]+ pglsfit$model$coef[1])#WCP diff --git a/Set up data.R b/Set up data.R index bdb5b15..947c2ab 100644 --- a/Set up data.R +++ b/Set up data.R @@ -6,40 +6,43 @@ library(ggpubr) library(ggplot2) library(tidyr) -#set working directory and load data +#set working directory setwd("E:/ATclone/A_T-stats") #load main dataframe -df<-read.csv("data Nov 17_.csv", stringsAsFactors = F, header = T) #, stringsAsFactors = FALSE +df<-read.csv("Eardata.csv", stringsAsFactors = F, header = T) #, stringsAsFactors = FALSE -#####The pgls model function, which will be applied to list of formulas### +#####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))) } -#note some missing headmass values now imputed +#note some missing headmass values df$Head.mass..g. -#load phylogeny and correct names that were different between birdtree.org and the up-to-date species names +#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") #add computed head mass from head mass~skullwidth pgls source("SW_HM_.R")#add phylogeney here df$Head.mass..g. -#Since PGLS uses one point per species,I make the dataframe to have average values for species with more than one specimen: -#First make a dataframe with only one species per line +#since PGLS uses one point per species, +#we create a dataframeaverage values for species with more than one specimen: + +#first make a dataframe with only one species per line distinctdf<-distinct(df, Binomial, .keep_all = TRUE) distinctdforder<-arrange(distinctdf,Binomial)#sort by species name -#Then get averages for columns with continuous data +#next, get averages for columns with continuous data avgdf<-df %>% group_by(Binomial) %>% summarise_at(vars(Skull.width..mm.:area_ratio), mean, na.rm = TRUE) avgdf<-as.data.frame(avgdf) -#And add back columns from distinct df which don't require averaging +#Nnow we add back columns from distinctdf which don't require averaging avgdf$Species<-distinctdforder$Species avgdf$Low.Hz<-distinctdforder$Low.Hz avgdf$Order<-distinctdforder$Order @@ -52,18 +55,18 @@ avgdf$IBP_detail<-distinctdforder$IBP_detail avgdf$Behind.TM<-distinctdforder$Behind.TM avgdf$`fluid.filled.`<-distinctdforder$`fluid.filled.` -#add dive depth dataframe here. This script groups the 'surface foraging' based on more detailed ecologies +#add dive depth data to the main dataframe +#this script groups the 'surface foraging' based on more detailed ecologies source("add_dive_depth_data.R") -#'plungedistin -#set 'terrestrial' as reference level +#set 'Terrestrial' as reference level for the ecological grouping variable avgdf$plungedistinct<-as.character(avgdf$catfeeding2) avgdf$plungedistinct[which(is.na(avgdf$plungedistinct))]<-"Terrestrial" avgdf$plungedistinct<-relevel(as.factor(avgdf$plungedistinct), ref = "Terrestrial") avgdf$plungedistinct[avgdf$Binomial=="Ardea_melanocephala"]<-"Surface" -#made data frame object -birdCDO<-comparative.data(phy = birdtreels,data = avgdf,#[avgdf$Category!="Terrestrial",] +#make comparative data frame object for the pgls function +birdCDO<-comparative.data(phy = birdtreels,data = avgdf, names.col = Binomial, vcv = TRUE, na.omit = FALSE, warn.dropped = TRUE) @@ -112,7 +115,7 @@ pgls_todo_nogeomet <- c("log(TMtotalarea)~log(Skull.width..mm.)", pgls_todo_hm<-pgls_todo_nogeomet[seq(2,length(pgls_todo_nogeomet),2)] -################### Data object for aquatic-only analyses ----------------------------- +################### Data object for aquatic-only analyses ###################----------------------------- birdCDO<-comparative.data(phy = birdtreels,data = avgdf[avgdf$Category!="Terrestrial",], names.col = Binomial, vcv = TRUE, na.omit = FALSE, @@ -121,7 +124,7 @@ birdCDO<-comparative.data(phy = birdtreels,data = avgdf[avgdf$Category!="Terrest #check any tips dropped between linking phylogeny and dataframe birdCDO$dropped -#Remove 'terrestrial' as a level +#remove 'terrestrial' as a level birdCDO$data$plungedistinct<-droplevels(birdCDO$data$plungedistinct, exclude = "Terrestrial") levels(birdCDO$data$plungedistinct) birdCDO$data$plungedistinct<-relevel(birdCDO$data$plungedistinct, ref = "Surface") diff --git a/Tblog.R b/Tblog.R index 8b1177d..5752cfc 100644 --- a/Tblog.R +++ b/Tblog.R @@ -1,4 +1,4 @@ -#code used from +#code based on functions from: https://thackl.github.io/ggtree-composite-plots library(ggtree) library(patchwork) diff --git a/add_dive_depth_data.R b/add_dive_depth_data.R index 9a9cadf..85b24a5 100644 --- a/add_dive_depth_data.R +++ b/add_dive_depth_data.R @@ -1,7 +1,8 @@ +#load dive data file divedf<-read.csv("Depth list.csv") divedf<-divedf[1:57,] -#SURFACE +#Classify SURFACE category from more detailed foraging modes divedf$feedmode2<-NA #reduce ecological groups divedf$feedmode2<-ifelse(divedf$Categories_feeding=="Dipping", "Surface", divedf$feedmode2) divedf$feedmode2<-ifelse(divedf$Categories_feeding=="Surface seizing - plunge diving?", "Surface", divedf$feedmode2) @@ -10,20 +11,17 @@ divedf$feedmode2<-ifelse(divedf$Categories_feeding=="Pattering", "Surface", dive divedf$feedmode2<-ifelse(divedf$Categories_feeding=="Surface seizing", "Surface", divedf$feedmode2) divedf$feedmode2<-ifelse(divedf$Categories_feeding=="Surface diving", "Surface", divedf$feedmode2) -#UNDERWATER PURSUIT (non-plunging) +#Classify UNDERWATER PURSUIT category from more detailed foraging modes divedf$feedmode2<-ifelse(divedf$Categories_feeding=="Pursuit diving", "Underwater pursuit", divedf$feedmode2) divedf$feedmode2<-ifelse(divedf$Categories_feeding=="Shallow swimming", "Underwater pursuit", divedf$feedmode2) -#PLUNGING (surface & pursuit plunging) +#Classify PLUNGING (surface & pursuit plunging) category from more detailed foraging modes divedf$feedmode2<-ifelse(divedf$Categories_feeding=="Surface plunging", "Plunging", divedf$feedmode2) divedf$feedmode2<-ifelse(divedf$Categories_feeding=="Pursuit plunging", "Plunging", divedf$feedmode2) divedf$feedmode2<-ifelse(divedf$Categories_feeding=="Surface plunging", "Plunging", divedf$feedmode2) divedf$feedmode2<-ifelse(divedf$Categories_feeding=="Surface plunging", "Plunging", divedf$feedmode2) -#species that didn't match -#ind<-which(divedf$IOC_2020_Binomial %in% avgdf$Binomial == T)#corresponding row number in main dataframe - -#check of all from divedf match to main df +#Match the dive dataframe to the main dataframe by species name match(divedf$IOC_2020_Binomial,avgdf$Binomial) #a list of the rows in the main df. #add dive data to main df diff --git a/pgls_HM.R b/pgls_HM.R index eff9307..737e87d 100644 --- a/pgls_HM.R +++ b/pgls_HM.R @@ -1,6 +1,6 @@ pgls_todo_hm<-pgls_todo_nogeomet[seq(2,length(pgls_todo_nogeomet),2)] -#Model list with head mass as only independent variable +#model list with head mass as only independent variable modellist<-pgls_todo_hm pgls_models_list<-lapply(pgls_todo_hm,pgls_models)#run pgls @@ -25,16 +25,7 @@ for(i in seq_along(tbllist_HM)){ 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 - #tbllist_HM[[i]]$Fstat[2:nrow(tbllist_HM[[i]])]<-"" - #tbllist_HM[[i]]$Fstat_numdf[2:nrow(tbllist_HM[[i]])]<-"" - #tbllist_HM[[i]]$Fstat_dendf[2:nrow(tbllist_HM[[i]])]<-" " - ##tbllist_HM[[i]]$Model[2:nrow(tbllist_HM[[i]])]<-"" - #tbllist_HM[[i]]$Lambda[2:nrow(tbllist_HM[[i]])]<-"" - #tbllist_HM[[i]]$Adj_Rsquared[2:nrow(tbllist_HM[[i]])]<-"" - #tbllist_HM[[i]]$AICc[2:nrow(tbllist_HM[[i]])]<-"" - row.names(tbllist_HM[[i]])<-c()#remove row names + 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]]) } \ No newline at end of file diff --git a/pgls_HM_plus_divescore.R b/pgls_HM_plus_divescore.R index b520b31..ed46b74 100644 --- a/pgls_HM_plus_divescore.R +++ b/pgls_HM_plus_divescore.R @@ -22,16 +22,7 @@ for(i in seq_along(tbllist_HM_plus_divescore)){ character_cols<-unlist(lapply(tbllist_HM_plus_divescore[[i]], is.character)) numeric_cols <- unlist(lapply(tbllist_HM_plus_divescore[[i]], is.numeric))# Identify numeric columns tbllist_HM_plus_divescore[[i]]<-cbind(tbllist_HM_plus_divescore[[i]][,which(character_cols)],signif(tbllist_HM_plus_divescore[[i]][,which(numeric_cols)], digits = 2)) - #tbllist_HM_plus_divescore[[i]] <- tbllist_HM_plus_divescore[[i]][, c(6,11,8:10,7,5,1:4)]#change order of columns - #dplyr::select_if(tbllist_HM_plus_divescore[[i]], is.numeric)#select only numeric data - colnames(tbllist_HM_plus_divescore[[i]])[6]<-"P.val"#rename b/c flextable doesn't work will with the '>' sign - #tbllist_HM_plus_divescore[[i]]$Fstat[2:nrow(tbllist_HM_plus_divescore[[i]])]<-"" - #tbllist_HM_plus_divescore[[i]]$Fstat_numdf[2:nrow(tbllist_HM_plus_divescore[[i]])]<-"" - #tbllist_HM_plus_divescore[[i]]$Fstat_dendf[2:nrow(tbllist_HM_plus_divescore[[i]])]<-" " - ##tbllist_HM_plus_divescore[[i]]$Model[2:nrow(tbllist_HM_plus_divescore[[i]])]<-"" - #tbllist_HM_plus_divescore[[i]]$Lambda[2:nrow(tbllist_HM_plus_divescore[[i]])]<-"" - #tbllist_HM_plus_divescore[[i]]$Adj_Rsquared[2:nrow(tbllist_HM_plus_divescore[[i]])]<-"" - #tbllist_HM_plus_divescore[[i]]$AICc[2:nrow(tbllist_HM_plus_divescore[[i]])]<-"" - row.names(tbllist_HM_plus_divescore[[i]])<-c()#remove row names + colnames(tbllist_HM_plus_divescore[[i]])[6]<-"P.val"#rename b/c flextable doesn't work will with the '>' sign + row.names(tbllist_HM_plus_divescore[[i]])<-c()#remove row names print(tbllist_HM_plus_divescore[[i]]) } \ No newline at end of file diff --git a/pgls_HM_times_divescore.R b/pgls_HM_times_divescore.R index 4dc93b7..fb493d1 100644 --- a/pgls_HM_times_divescore.R +++ b/pgls_HM_times_divescore.R @@ -22,16 +22,7 @@ for(i in seq_along(tbllist_HM_times_divescore)){ character_cols<-unlist(lapply(tbllist_HM_times_divescore[[i]], is.character)) numeric_cols <- unlist(lapply(tbllist_HM_times_divescore[[i]], is.numeric))# Identify numeric columns tbllist_HM_times_divescore[[i]]<-cbind(tbllist_HM_times_divescore[[i]][,which(character_cols)],signif(tbllist_HM_times_divescore[[i]][,which(numeric_cols)], digits = 2)) - #tbllist_HM_times_divescore[[i]] <- tbllist_HM_times_divescore[[i]][, c(6,11,8:10,7,5,1:4)]#change order of columns - #dplyr::select_if(tbllist_HM_times_divescore[[i]], is.numeric)#select only numeric data - colnames(tbllist_HM_times_divescore[[i]])[6]<-"P.val"#rename b/c flextable doesn't work will with the '>' sign - #tbllist_HM_times_divescore[[i]]$Fstat[2:nrow(tbllist_HM_times_divescore[[i]])]<-"" - #tbllist_HM_times_divescore[[i]]$Fstat_numdf[2:nrow(tbllist_HM_times_divescore[[i]])]<-"" - #tbllist_HM_times_divescore[[i]]$Fstat_dendf[2:nrow(tbllist_HM_times_divescore[[i]])]<-" " - ##tbllist_HM_times_divescore[[i]]$Model[2:nrow(tbllist_HM_times_divescore[[i]])]<-"" - #tbllist_HM_times_divescore[[i]]$Lambda[2:nrow(tbllist_HM_times_divescore[[i]])]<-"" - #tbllist_HM_times_divescore[[i]]$Adj_Rsquared[2:nrow(tbllist_HM_times_divescore[[i]])]<-"" - #tbllist_HM_times_divescore[[i]]$AICc[2:nrow(tbllist_HM_times_divescore[[i]])]<-"" - row.names(tbllist_HM_times_divescore[[i]])<-c()#remove row names + colnames(tbllist_HM_times_divescore[[i]])[6]<-"P.val"#rename b/c flextable doesn't work will with the '>' sign + row.names(tbllist_HM_times_divescore[[i]])<-c()#remove row names print(tbllist_HM_times_divescore[[i]]) } \ No newline at end of file diff --git a/plot ecology on circular phylogeny.R b/plot ecology on circular phylogeny.R index 1239a0d..312070e 100644 --- a/plot ecology on circular phylogeny.R +++ b/plot ecology on circular phylogeny.R @@ -5,16 +5,14 @@ library(colorspace) divecol<-c(sequential_hcl(6, palette = "Purple-Blue",rev = T)[2:6]) -#count # in each order +#count number of orders present in the full phylogenetic tree orderdf<-avgdf %>% count(Order) -#distinct df +#make a phylogeny with each unique order present g2<-avgdf[!duplicated(avgdf$Order),]#orders only g2order<-arrange(g2,Order)#sort to match order df g2order$n<-orderdf$n#attach number of species in order str(g2order$Binomial) -#drop timps for order - g2order$full<-paste0(g2order$Order," (",as.character(g2order$n),")") g2order$Order<-as.character(g2order$Order) @@ -23,7 +21,7 @@ mypal <- colorRampPalette(brewer.pal(6, "Blues")) mypal2 <- colorRampPalette(brewer.pal(6, "YlOrRd")) -#clade labels +#make clade labels A<-findMRCA(birdtreels,avgdf$Binomial[avgdf$Order=="Passeriformes"]) B<-findMRCA(birdtreels,avgdf$Binomial[avgdf$Order=="Charadriiformes"]) C<-findMRCA(birdtreels,avgdf$Binomial[avgdf$Order=="Falconiformes"]) @@ -42,7 +40,7 @@ O<-findMRCA(birdtreels,avgdf$Binomial[avgdf$Order=="Strigiformes"]) P<-findMRCA(birdtreels,avgdf$Binomial[avgdf$Order=="Anseriformes"]) - +#plot p<-ggtree(birdtreels, layout = "circular", open.angle = 150) %<+% avgdf + ###########, layout = "circular" #geom_tiplab(aes(label = Binomial2), linesize = 0.1, offset = 4) + #circular #geom_text(aes(label = node))+ @@ -95,7 +93,7 @@ p<-ggtree(birdtreels, layout = "circular", open.angle = 150) %<+% avgdf + ###### geom_strip(g2order$Binomial[5], g2order$Binomial[5], offset=16, offset.text=15, hjust=0, fontsize=3, label=g2order$full[5]) p -###########ESSENTIAL to change row names forheatmap +###########ESSENTIAL step to change row names forheatmap############# rownames(avgdf)<-avgdf$Binomial # The palette with black: @@ -106,7 +104,8 @@ BLUE<-c("#d0d1e6", "#74a9cf", "#2b8cbe", "#045a8d") -#add heatmap + +#add heatmap to circular phylogeny k<-gheatmap(p,avgdf[,c("plungedistinct","divescore")], #"Category", width = 0.25, offset = 0, color = "black", @@ -125,13 +124,13 @@ l ggsave("D:/Analysis_plots/ecolcircle.pdf", width=10, height=10) -#############plot sampling +#############plot heatmap showing sampling of all auditory traits (supplemental materials)#####333 avgdf$Binomial2<-avgdf$Binomial a<-ggtree(birdtreels) %<+% avgdf + ###########, layout = "circular" geom_tiplab(aes(label = Binomial2), linesize = 0.1, offset = 30) + #circular #geom_text(aes(label = node))+ - xlim(NA, 300) + + xlim(NA, 160) + ylim(NA,140) a