.Rhistory

###END/\/\/\/\/\/\/\/\/\/\/\###
mtext(text=paste0("PC1 (",round(100*(diag(body_pca$Eval)[1]/sum(diag(body_pca$Eval))),2),"%)"),side=1,line=1.3,cex=0.55)
mtext(text=paste0("PC2 (",round(100*(diag(body_pca$Eval)[2]/sum(diag(body_pca$Eval))),2),"%)"),side=2,line=1.3,cex=0.55)
## Figure panel label ###
text(x=par("usr")[1]+diff(c(par("usr")[1],par("usr")[2]))*0.075,y=par("usr")[4]-diff(c(par("usr")[3],par("usr")[4]))*0.1,label=LETTERS[4],font=2,cex=1.5)
### Bird Pictures Panel ###
#par(mar=c(2,0.5,0,14))
require(magick)
# Crop images
#img<-image_crop(example, "100x150+50")#crop out width:100px and height:150px starting +50px from the left
T_savana<-magick::image_read("50763510127_317508235f_o.jpg")
T_savana<-image_scale(T_savana, "200") # width: 200px
T_savana<-image_crop(T_savana, "200x300+0+20")
T_savana
T_caudifasciatus<-magick::image_read("T_caudifasciatus.png")
T_caudifasciatus<-image_scale(T_caudifasciatus,"600")
T_caudifasciatus<-image_crop(T_caudifasciatus, "300x280+120+10")
T_caudifasciatus
T_crassirostris<-magick::image_read("50250535916_399b2a0dea_o.jpg")
T_crassirostris<-image_scale(T_crassirostris,"275")
T_crassirostris<-image_crop(T_crassirostris,"250x280+40+0")
T_crassirostris
T_cubensis<-magick::image_read("Giant Kingbird 2 HR.jpg")
T_cubensis<-image_scale(T_cubensis,"450")
T_cubensis<-image_crop(T_cubensis,"210x280+190")
T_cubensis
photos<-image_append(c(T_savana, T_caudifasciatus, T_crassirostris, T_cubensis))
par(mar=c(0,0,0,0))
plot(photos)
points(x=180, y=260, cex=2, pch=coldf[phy$tip.label,]$pch[26],bg=coldf[phy$tip.label,]$col[26], col=coldf[phy$tip.label,]$outlinecolor[26],lwd=coldf[phy$tip.label,]$outlinelwd[26],lty=3)
points(x=480, y=260, cex=2, pch=coldf[phy$tip.label,]$pch[21],bg=coldf[phy$tip.label,]$col[21], col=coldf[phy$tip.label,]$outlinecolor[21],lwd=coldf[phy$tip.label,]$outlinelwd[21],lty=3)
points(x=720, y=260, cex=2, pch=coldf[phy$tip.label,]$pch[14],bg=coldf[phy$tip.label,]$col[14], col=coldf[phy$tip.label,]$outlinecolor[14],lwd=coldf[phy$tip.label,]$outlinelwd[14],lty=3)
points(x=930, y=260, cex=2, pch=coldf[phy$tip.label,]$pch[4],bg=coldf[phy$tip.label,]$col[4], col=coldf[phy$tip.label,]$outlinecolor[4],lwd=coldf[phy$tip.label,]$outlinelwd[4],lty=3)
text(x=20, y=260, paste0("E"), font=2, cex=1.5)
text(x=220, y=260, paste0("F"),font=2, cex=1.5)
text(x=520, y=260, paste0("G"),font=2, cex=1.5)
text(x=750, y=260, paste0("H"),font=2, cex=1.5)
dev.off()
rm(list=ls())
### Set source directory to the folder this file came from within RStudio
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
#Note: If assessing sexes separately, use one of the following lines of code instead of the above line, then proceed.
morpho<-read.csv('./Output Files/Tyrannus_Females_Adults_data.csv', row.names=2)
### Read in phylogeny
phy<-read.tree('./Output Files/Tyrannus_phylogeny.tre')
### Check that tree$tip.label is the same as morpho
phy$tip.label[!phy$tip.label %in% rownames(morpho)]
### Bill Phylogenetic PCA ###
bill_pca<-phyl.pca(phy,morpho[,3:5],method="lambda")
bill_pca
### Find the % variance explained by each eigenvector
diag(bill_pca$Eval)/sum(bill_pca$Eval)*100
#Note: If assessing sexes separately, use one of the following lines of code instead of the above line, then proceed.
morpho_whole<-read.csv('./Output Files/Tyrannus Adults Females morphology data.csv', row.names=1)
### Match tip.labels with phylogeny tip.labels ###
morpho_whole$tip.label
names(table(morpho_whole$tip.label))
phy$tip.label
phy$tip.label[!phy$tip.label %in% morpho_whole$tip.label] #should return character(0)
morpho_whole$tip.label[! morpho_whole$tip.label %in% phy$tip.label] #should return character(0)
### Split the morphology data sets up into categories by tip.label
morpho_split<-split(morpho_whole,morpho_whole$tip.label)
names(morpho_split)
phy$tip.label
phy$tip.label[!phy$tip.label %in% names(morpho_split)] #should return character(0)
names(morpho_split)[! names(morpho_split) %in% phy$tip.label] #should return character (0)
morpho_trim<-lapply(morpho_split,function(x) x[colnames(x) %in% c("BL.Average","BW.Average","BD.Average", "Kipp.s.Distance", "Kipp.s.Index", "WC.Average", "Tail", "Tarsus.Average", "name_num")])
bill_mat<-matrix(as.numeric(as.matrix(do.call(rbind, lapply(morpho_trim,function(x) x[,1:3])))),nrow=nrow(do.call(rbind,morpho_trim)))
bill_scores<-phytools:::scores(bill_pca,newdata=bill_mat)
rownames(bill_scores)<-rownames(do.call(rbind, lapply(morpho_trim,function(x) x[,1:3])))
bill_scores_df<-data.frame(bill_scores)
bill_scores_df$otu<-factor(sapply(strsplit(rownames(bill_scores_df),"[.]"),function(x) x[1]))
bill_scores_split<-split(bill_scores_df, bill_scores_df$otu)
bill_scores_avg<-t(sapply(bill_scores_split,function(x) c(mean(na.omit(x[,1])),mean(na.omit(x[,2])))))
bill_scores_avg<-as.data.frame(bill_scores_avg)
morpho$BillPC1<-bill_scores_avg$V1
morpho$BillPC2<-bill_scores_avg$V2
### Flight Feather Phylogenetic PCA ###
body_pca<-phyl.pca(phy,morpho[,7:9],method="lambda")
body_pca
### Find the % variance explained by each eigenvector
diag(body_pca$Eval)/sum(body_pca$Eval)*100
body_mat<-matrix(as.numeric(as.matrix(do.call(rbind, lapply(morpho_trim,function(x) x[,5:7])))),nrow=nrow(do.call(rbind,morpho_trim)))
body_scores<-phytools:::scores(body_pca,newdata=body_mat)
rownames(body_scores)<-rownames(do.call(rbind, lapply(morpho_trim,function(x) x[,5:7])))
body_scores_df<-data.frame(body_scores)
body_scores_df$otu<-factor(sapply(strsplit(rownames(body_scores_df),"[.]"),function(x) x[1]))
body_scores_split<-split(body_scores_df, body_scores_df$otu)
body_scores_avg<-t(sapply(body_scores_split,function(x) c(mean(na.omit(x[,1])),mean(na.omit(x[,2])))))
body_scores_avg<-as.data.frame(body_scores_avg)
morpho$BodyPC1<-body_scores_avg$V1
morpho$BodyPC2<-body_scores_avg$V2
### Body PCA without T. forficatus or T. savana bc their tails are so long ###
phy_trim<-drop.tip(phy,c(grep("forficatus",phy$tip.label),grep("savana",phy$tip.label)))
otu_avg_trim<-morpho[rownames(morpho) %in% phy_trim$tip.label,]
body_pca_trim<-phyl.pca(phy_trim,otu_avg_trim[,7:9],method="lambda")
body_pca_trim
### Find the % variance explained by each eigenvector
diag(body_pca_trim$Eval)/sum(body_pca_trim$Eval)*100
body_mat_trim<-matrix(as.numeric(as.matrix(do.call(rbind, lapply(morpho_trim[-c(grep("forficatus",names(morpho_trim)),grep("savana",names(morpho_trim)))],function(x) x[,5:7])))),nrow=nrow(do.call(rbind,morpho_trim[-c(grep("forficatus",names(morpho_trim)),grep("savana",names(morpho_trim)))])))
body_scores_trim<-phytools:::scores(body_pca_trim,newdata=body_mat_trim)
rownames(body_scores_trim)<-rownames(do.call(rbind, lapply(morpho_trim[-c(grep("forficatus",names(morpho_trim)),grep("savana",names(morpho_trim)))],function(x) x[,5:7])))
body_scores_trim_df<-data.frame(body_scores_trim)
body_scores_trim_df$otu<-factor(sapply(strsplit(rownames(body_scores_trim_df),"[.]"),function(x) x[1]))
body_scores_split_trim<-split(body_scores_trim_df, body_scores_trim_df$otu)
body_scores_trim_avg<-t(sapply(body_scores_split_trim,function(x) c(mean(na.omit(x[,1])),mean(na.omit(x[,2])))))
#Note: If assessing sexes separately, use one of the following lines of code instead of the above line, then proceed.
cv_summary<-read.csv('./Output Files/cv_summary_Females_Adults_table.csv', row.names = 2)
### Coefficient of Variation (CV) Bill Phylogenetic PCA ###
CVbill_pca<-phyl.pca(phy,cv_summary[,3:5],method="lambda")
CVbill_pca
### Find the % variance explained by each eigenvector
diag(CVbill_pca$Eval)/sum(CVbill_pca$Eval)*100
CVbill_mat<-matrix(as.numeric(as.matrix(do.call(rbind, lapply(morpho_trim,function(x) x[,1:3])))),nrow=nrow(do.call(rbind,morpho_trim)))
CVbill_scores<-phytools:::scores(CVbill_pca,newdata=CVbill_mat)
rownames(CVbill_scores)<-rownames(do.call(rbind, lapply(morpho_trim,function(x) x[,1:3])))
CVbill_scores_df<-data.frame(CVbill_scores)
CVbill_scores_df$otu<-factor(sapply(strsplit(rownames(CVbill_scores_df),"[.]"),function(x) x[1]))
CVbill_scores_split<-split(CVbill_scores_df, CVbill_scores_df$otu)
CVbill_scores_cv<-t(sapply(CVbill_scores_split,function(x) c(mean(na.omit(x[,1])),mean(na.omit(x[,2])))))
CVbill_scores_cv<-as.data.frame(CVbill_scores_cv)
cv_summary$CVBillPC1<-CVbill_scores_cv$V1
cv_summary$CVBillPC2<-CVbill_scores_cv$V2
### CV Flight Feather Phylogenetic PCA ###
CVbody_pca<-phyl.pca(phy,cv_summary[,7:9],method="lambda")
#Note: If assessing females, use one of the following lines of code instead of the above line, then proceed.
CVbody_pca<-phyl.pca(phy,na.omit(cv_summary[,6:8]),method="lambda") #there is a missing value in the tail column
CVbody_pca
### Find the % variance explained by each eigenvector
diag(CVbody_pca$Eval)/sum(CVbody_pca$Eval)*100
CVbody_mat<-matrix(as.numeric(as.matrix(do.call(rbind, lapply(morpho_trim,function(x) x[,5:7])))),nrow=nrow(do.call(rbind,morpho_trim)))
CVbody_scores<-phytools:::scores(CVbody_pca,newdata=CVbody_mat)
rownames(CVbody_scores)<-rownames(do.call(rbind, lapply(morpho_trim,function(x) x[,5:7])))
CVbody_scores_df<-data.frame(CVbody_scores)
CVbody_scores_df$otu<-factor(sapply(strsplit(rownames(CVbody_scores_df),"[.]"),function(x) x[1]))
CVbody_scores_split<-split(CVbody_scores_df, CVbody_scores_df$otu)
CVbody_scores_cv<-t(sapply(CVbody_scores_split,function(x) c(mean(na.omit(x[,1])),mean(na.omit(x[,2])))))
CVbody_scores_cv<-as.data.frame(CVbody_scores_cv)
cv_summary$CVBodyPC1<-CVbody_scores_cv$V1
cv_summary$CVBodyPC2<-CVbody_scores_cv$V2
### CV Body PCA without T. forficatus or T. savana bc their tails are so long ###
phy_trim<-drop.tip(phy,c(grep("forficatus",phy$tip.label),grep("savana",phy$tip.label)))
otu_cv_trim<-cv_summary[rownames(cv_summary) %in% phy_trim$tip.label,]
CVbody_pca_trim<-phyl.pca(phy_trim,otu_cv_trim[,7:9],method="lambda")
CVbody_pca_trim
### Build Ancestral State Reconstruction objects for phylogeny in figure
is.binary(phy)
phy<-multi2di(phy) #resolve polytomies
is.binary(phy)
#Note: If assessing sexes separately, use one of the following lines of code instead of the above line, then proceed.
Tyrannus.data<-read.csv('./Output Files/Tyrannus morphology + PCA avg_Females_Adults.csv', row.names = 1)
species<-rownames(Tyrannus.data)
Tyrannus.data<-cbind(species,Tyrannus.data)
### Check that tree$tip.label is the same as otu_avg
phy$tip.label[!phy$tip.label %in% Tyrannus.data$species] #should return: charactor(0)
### Use diversitree, the mk-n model to reconstruct marginal ancestral states
require(diversitree)
Migration<-factor(Tyrannus.data$Strategy)
names(Migration)<-rownames(Tyrannus.data)
#Change character to numeric, because this is the format diversitree accepts
char1<-as.numeric(Migration)
names(char1)<-names(Migration)
#Create mkn model. Use equal probabilites for all states at the root
lik.mn<-make.mkn(phy,char1,k=3, control=list(root=ROOT.EQUI))
rm(list=ls())
require("phytools")
### Set source directory to the folder this file came from within RStudio
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
### Bring in base phylogeny ###
tyrannus_base<-read.tree("MacPherson_Tyrannus_base.tre")
### make tree ultrametric ###
is.ultrametric(tyrannus_base)
tyrannus_um<-force.ultrametric(tyrannus_base,method="nnls")
#Rooted; includes branch lengths.
plot(tyrannus_um)
### Find edge lengths
require(ape)
node.depth.edgelength(tyrannus_um)
### read in list of missing taxa sensu clements 2019 ###
taxa2add<-readLines(file("Tyrannus_taxa2add.txt"))
taxa2add<-gsub(" ","_",taxa2add)
### Add in T. savana monachus tips based on Gomez-Bahamon phylogeny ###
taxa2add<- taxa2add[-10] #remove T. savana monachus as a single OTU because we are separating Central American and South American populations
taxa2add
tyrannus_mod <-bind.tip(tyrannus_um,tip.label="Tyrannus_savana_monachus_CA",where=which(tyrannus_um $tip.label=="Tyrannus_savana_savana"),position=1.4)
tyrannus_mod <-bind.tip(tyrannus_mod,tip.label="Tyrannus_savana_monachus_SA",where=which(tyrannus_mod $tip.label=="Tyrannus_savana_monachus_CA"),position=(3682/1e6))
plot(tyrannus_mod)
### Add missing tips to the phylogeny assuming a MRCA of 0.5 mya ### Note: change the 'position' argument below to add missing tips at different branch lengths
for(i in 1:length(taxa2add)){
taxa2add[i]
binomial<-paste(strsplit(taxa2add[i],"_")[[1]][1:2],collapse="_")
recipient<-tyrannus_mod $tip.label[grep(binomial, tyrannus_mod $tip.label)]
if(length(recipient)>1){
recipient<-findMRCA(tyrannus_mod,recipient)
tyrannus_mod <-bind.tip(tyrannus_mod,taxa2add[i],where=recipient)
}else{
tyrannus_mod <-bind.tip(tyrannus_mod,taxa2add[i],where=which(tyrannus_mod $tip.label== recipient),position=0.2) #the position argument is what determines the terminal branch length for the added taxon in millions of years
}
}
plot(tyrannus_mod)
### Replace any 0-length branches with a very small value ###
tyrannus_mod$edge.length[tyrannus_mod$edge.length==0]<-0.001
plot(tyrannus_mod)
### Write out phylogeny ###
write.tree(tyrannus_mod,file="Tyrannus_phylogeny_0.2.tre")
rm(list=ls())
### Bring in base phylogeny ###
tyrannus_base<-read.tree("MacPherson_Tyrannus_base.tre")
### make tree ultrametric ###
is.ultrametric(tyrannus_base)
tyrannus_um<-force.ultrametric(tyrannus_base,method="nnls")
#Rooted; includes branch lengths.
plot(tyrannus_um)
node.depth.edgelength(tyrannus_um)
### read in list of missing taxa sensu clements 2019 ###
taxa2add<-readLines(file("Tyrannus_taxa2add.txt"))
taxa2add<-gsub(" ","_",taxa2add)
### Add in T. savana monachus tips based on Gomez-Bahamon phylogeny ###
taxa2add<- taxa2add[-10] #remove T. savana monachus as a single OTU because we are separating Central American and South American populations
taxa2add
tyrannus_mod <-bind.tip(tyrannus_um,tip.label="Tyrannus_savana_monachus_CA",where=which(tyrannus_um $tip.label=="Tyrannus_savana_savana"),position=1.4)
tyrannus_mod <-bind.tip(tyrannus_mod,tip.label="Tyrannus_savana_monachus_SA",where=which(tyrannus_mod $tip.label=="Tyrannus_savana_monachus_CA"),position=(3682/1e6))
plot(tyrannus_mod)
### Add missing tips to the phylogeny assuming a MRCA of 0.5 mya ### Note: change the 'position' argument below to add missing tips at different branch lengths
for(i in 1:length(taxa2add)){
taxa2add[i]
binomial<-paste(strsplit(taxa2add[i],"_")[[1]][1:2],collapse="_")
recipient<-tyrannus_mod $tip.label[grep(binomial, tyrannus_mod $tip.label)]
if(length(recipient)>1){
recipient<-findMRCA(tyrannus_mod,recipient)
tyrannus_mod <-bind.tip(tyrannus_mod,taxa2add[i],where=recipient)
}else{
tyrannus_mod <-bind.tip(tyrannus_mod,taxa2add[i],where=which(tyrannus_mod $tip.label== recipient),position=0.93) #the position argument is what determines the terminal branch length for the added taxon in millions of years
}
}
plot(tyrannus_mod)
### Replace any 0-length branches with a very small value ###
tyrannus_mod$edge.length[tyrannus_mod$edge.length==0]<-0.001
plot(tyrannus_mod)
### Write out phylogeny ###
write.tree(tyrannus_mod,file="Tyrannus_phylogeny_0.93.tre")
rm(list=ls())
### Set source directory to the folder this file came from within RStudio
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
phy<-read.tree('./Output Files/Tyrannus_phylogeny_0.2.tre')
### Bring in morphology dataset with average morphology values for each OTU and PC scores
otu_avg<-read.csv('./Output Files/Tyrannus morphology + PCA avg_Adults.csv', row.names = 1)
species<-rownames(otu_avg)
otu_avg<-cbind(species,otu_avg)
### Check that tree$tip.label is the same as otu_avg
phy$tip.label[!phy$tip.label %in% otu_avg$species]
### phylANOVA on OTU morphometric means while controlling for body size ###
response_var<-colnames(otu_avg)[c(4,5,6,8,9,10)]
### Set up predictor variable for phylANOVA ###
migration<-otu_avg$Strategy
names(migration)<-rownames(otu_avg)
### Set up tarsus for size correction #
otu_avg$Tarsus.Average <-as.numeric(otu_avg$Tarsus.Average)
tarsus.length<-otu_avg$Tarsus.Average
names(tarsus.length)<-rownames(otu_avg)
### Create empty list to store output of each for loop iteration
phyl.resid_output<-list()
phylANOVA_output<-list()
### Conduct the phylogenetic ANOVA he following code works without assuming y is the same order as tree$tip.label
for(i in 1:length(response_var)){
print(paste0("Character ",i," -- ",response_var[i]))#Report character
character_of_interest<-otu_avg[,response_var[i]] #Extract character of interest (coi)
names(character_of_interest)<-rownames(otu_avg) #Add names to vector of coi
phyl.resid_output[[i]]<-phyl.resid(phy,tarsus.length, character_of_interest,method="lambda") #Perform phylogenetic residuals
phylANOVA_output[[i]]<-phylANOVA(phy, x=migration, y=phyl.resid_output[[i]]$resid[,1], p.adj="none")
}
names(phyl.resid_output)<-response_var
names(phylANOVA_output)<-response_var
sink("phylANOVA_output_0.2.txt")
for(i in 1:length(phylANOVA_output)){
cat(paste0("Character ",i," -- ",response_var[i]))
cat("\n\n")
print(phylANOVA_output[[i]])
cat("#############")
cat("\n\n")
}
sink()
### Conduct the phylogenetic ANVOA with p.adj="bonferroni" to account for multiple hypothesis testing
## Create empty list to store output of each for loop iteration
phyl.resid_output<-list()
phylANOVA_output<-list()
###################################################################################
### Conduct the phylogenetic ANOVA for metrics that do not need to be size-corrected ###
phy$tip.label[!phy$tip.label %in% rownames(otu_avg)]
response_var2<-colnames(otu_avg)[c(12,13)]
phylANOVA_PCoutput<-list()
for(i in 1:length(response_var2)){
print(paste0("Character ",i," -- ",response_var2[i]))#Report character
character_of_interest<-otu_avg[,response_var2[i]] #Extract character of interest (coi)
names(character_of_interest)<-rownames(otu_avg) #Add names to vector of coi
phylANOVA_PCoutput[[i]]<-phylANOVA(tree=phy, x=migration, y=character_of_interest, p.adj="none")
}
names(phylANOVA_PCoutput)<-response_var2
sink("phylANOVA_output_nosizecorrectionneeded_0.2.txt")
for(i in 1:length(phylANOVA_PCoutput)){
cat(paste0("Character ",i," -- ",response_var2[i]))
cat("\n\n")
print(phylANOVA_PCoutput[[i]])
cat("#############")
cat("\n\n")
}
sink()
########################################################
### Repeat for CV data
otu_cv_avg<-read.csv('./Output Files/cv_summary_Adults.csv', row.names = 1)
species<-rownames(otu_cv_avg)
otu_cv_avg<-cbind(species,otu_cv_avg)
### Check that tree$tip.label is the same as otu_avg
phy$tip.label[!phy$tip.label %in% otu_cv_avg$species]
response_var3<-colnames(otu_cv_avg)[c(4,5,6,7,8,9,10,12,13)]
### Set up predictor variable for phylANOVA ###
migration<-otu_cv_avg$Strategy
names(migration)<-rownames(otu_cv_avg)
#Create empty list to store output of each for loop iteration
phylANOVA_CVoutput<-list()
### Conduct the phylogenetic ANOVA on coefficients of variation ###
for(i in 1:length(response_var3)){
print(paste0("Character ",i," -- ",response_var3[i]))#Report character
character_of_interest<-otu_cv_avg[,response_var3[i]] #Extract character of interest (coi)
names(character_of_interest)<-rownames(otu_cv_avg) #Add names to vector of coi
phylANOVA_CVoutput[[i]]<-phylANOVA(phy, x=migration, y=character_of_interest, p.adj="none")
}
names(phylANOVA_CVoutput)<-response_var3
sink("phylANOVA_cv_output_0.2.txt")
for(i in 1:length(phylANOVA_CVoutput)){
cat(paste0("Character ",i," -- ",response_var3[i]))
cat("\n\n")
print(phylANOVA_CVoutput[[i]])
cat("#############")
cat("\n\n")
}
sink()
rm(list=ls())
#phy<-read.tree('./Output Files/Tyrannus_phylogeny_0.2.tre')
phy<-read.tree('./Output Files/Tyrannus_phylogeny_0.93.tre')
### Bring in morphology dataset with average morphology values for each OTU and PC scores
otu_avg<-read.csv('./Output Files/Tyrannus morphology + PCA avg_Adults.csv', row.names = 1)
species<-rownames(otu_avg)
otu_avg<-cbind(species,otu_avg)
### Check that tree$tip.label is the same as otu_avg
phy$tip.label[!phy$tip.label %in% otu_avg$species]
### phylANOVA on OTU morphometric means while controlling for body size ###
response_var<-colnames(otu_avg)[c(4,5,6,8,9,10)]
### Set up predictor variable for phylANOVA ###
migration<-otu_avg$Strategy
names(migration)<-rownames(otu_avg)
### Set up tarsus for size correction #
otu_avg$Tarsus.Average <-as.numeric(otu_avg$Tarsus.Average)
tarsus.length<-otu_avg$Tarsus.Average
names(tarsus.length)<-rownames(otu_avg)
### Create empty list to store output of each for loop iteration
phyl.resid_output<-list()
phylANOVA_output<-list()
### Conduct the phylogenetic ANOVA he following code works without assuming y is the same order as tree$tip.label
for(i in 1:length(response_var)){
print(paste0("Character ",i," -- ",response_var[i]))#Report character
character_of_interest<-otu_avg[,response_var[i]] #Extract character of interest (coi)
names(character_of_interest)<-rownames(otu_avg) #Add names to vector of coi
phyl.resid_output[[i]]<-phyl.resid(phy,tarsus.length, character_of_interest,method="lambda") #Perform phylogenetic residuals
phylANOVA_output[[i]]<-phylANOVA(phy, x=migration, y=phyl.resid_output[[i]]$resid[,1], p.adj="none")
}
names(phyl.resid_output)<-response_var
names(phylANOVA_output)<-response_var
#sink("phylANOVA_output_0.2.txt")
sink("phylANOVA_output_0.93.txt")
for(i in 1:length(phylANOVA_output)){
cat(paste0("Character ",i," -- ",response_var[i]))
cat("\n\n")
print(phylANOVA_output[[i]])
cat("#############")
cat("\n\n")
}
sink()
###################################################################################
### Conduct the phylogenetic ANOVA for metrics that do not need to be size-corrected ###
phy$tip.label[!phy$tip.label %in% rownames(otu_avg)]
response_var2<-colnames(otu_avg)[c(12,13)]
phylANOVA_PCoutput<-list()
for(i in 1:length(response_var2)){
print(paste0("Character ",i," -- ",response_var2[i]))#Report character
character_of_interest<-otu_avg[,response_var2[i]] #Extract character of interest (coi)
names(character_of_interest)<-rownames(otu_avg) #Add names to vector of coi
phylANOVA_PCoutput[[i]]<-phylANOVA(tree=phy, x=migration, y=character_of_interest, p.adj="none")
}
names(phylANOVA_PCoutput)<-response_var2
#sink("phylANOVA_output_nosizecorrectionneeded_0.2.txt")
sink("phylANOVA_output_nosizecorrectionneeded_0.93.txt")
for(i in 1:length(phylANOVA_PCoutput)){
cat(paste0("Character ",i," -- ",response_var2[i]))
cat("\n\n")
print(phylANOVA_PCoutput[[i]])
cat("#############")
cat("\n\n")
}
sink()
## Conduct the phylogenetic ANVOA on metrics that do not need to be size-corrected with p.adj="bonferroni" to account for multiple hypothesis testing
phylANOVA_PCoutput<-list()
########################################################
### Repeat for CV data
otu_cv_avg<-read.csv('./Output Files/cv_summary_Adults.csv', row.names = 1)
species<-rownames(otu_cv_avg)
otu_cv_avg<-cbind(species,otu_cv_avg)
### Check that tree$tip.label is the same as otu_avg
phy$tip.label[!phy$tip.label %in% otu_cv_avg$species]
response_var3<-colnames(otu_cv_avg)[c(4,5,6,7,8,9,10,12,13)]
### Set up predictor variable for phylANOVA ###
migration<-otu_cv_avg$Strategy
names(migration)<-rownames(otu_cv_avg)
#Create empty list to store output of each for loop iteration
phylANOVA_CVoutput<-list()
### Conduct the phylogenetic ANOVA on coefficients of variation ###
for(i in 1:length(response_var3)){
print(paste0("Character ",i," -- ",response_var3[i]))#Report character
character_of_interest<-otu_cv_avg[,response_var3[i]] #Extract character of interest (coi)
names(character_of_interest)<-rownames(otu_cv_avg) #Add names to vector of coi
phylANOVA_CVoutput[[i]]<-phylANOVA(phy, x=migration, y=character_of_interest, p.adj="none")
}
names(phylANOVA_CVoutput)<-response_var3
#sink("phylANOVA_cv_output_0.2.txt")
sink("phylANOVA_cv_output_0.93.txt")
for(i in 1:length(phylANOVA_CVoutput)){
cat(paste0("Character ",i," -- ",response_var3[i]))
cat("\n\n")
print(phylANOVA_CVoutput[[i]])
cat("#############")
cat("\n\n")
}
sink()
rm(list=ls())
require(ggplot2)
require(tidyr)
require(cowplot)
require(gridExtra)
require(gridGraphics)
install.packages("ggplot2")
require(ggplot2)
require(tidyr)
require(cowplot)
require(gridExtra)
require(gridGraphics)
install.packages("gridGraphics")
require(gridGraphics)
### Set source directory to the folder this file came from within RStudio
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
### Note: If assessing sexes separately, use the following lines of code instead of the above 3 lines, then proceed.
otu_avg<-read.csv('./Output Files/Tyrannus morphology + PCA avg_Females_Adults.csv', row.names = 1)
otu_cv_avg<-read.csv('./Output Files/cv_summary_Females_Adults.csv', row.names = 1)
residuals<-read.csv('./Output Files/phylANOVA_tarsus-corrected_residuals_Females.csv', row.names = 1)
### Note: If assessing sexes separately, use the following lines of code instead of the above 3 lines, then proceed.
png(file="PPC scores boxplots_Females.png",width=7,height=2.75,units="in",res=500)
responsevariable.labs<-c("Bill PPC1", "Bill PPC2")
names(responsevariable.labs)<- c("BillPC1", "BillPC2")
col<-c("black","gray40","gray70","black","gray40","gray70")
p<-otu_avg %>%
pivot_longer(BillPC1:BillPC2, names_to="responsevariable", values_to = "value") %>%
ggplot(aes(y=value, x=Strategy)) +
geom_boxplot(color=col) +
facet_wrap(vars(responsevariable), ncol=2, labeller=labeller(responsevariable = responsevariable.labs), scales="free_y") +
panel_border() +
theme_bw() +
theme(panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
strip.text=element_text(size=15)) +
scale_y_continuous(expand=c(0.05,8)) +
labs(x="Migration Strategy", y= "PPC Scores")
draw_plot(plot, x = 0, y = 0, width = 1, height = 1)
ggdraw() +
draw_plot(p,0,0,1,1) +
draw_plot_label(c("A", "B"), c(0.085,0.565), c(0.99,0.99), size=15) + #main plot identifiers
draw_plot_label(c("A", "A", "A"), c(0.146, 0.287, 0.420), c(0.87, 0.87, 0.87), size=10) + #letters for plot A
draw_plot_label(c("A", "AB", "B"), c(0.625, 0.746, 0.9), c(0.87, 0.87, 0.87), size=10) #letters for plot B
dev.off()
rm(list=ls())
### Note: If assessing sexes separately, use the following lines of code instead of the above 3 lines, then proceed.
#otu_avg<-read.csv('./Output Files/Tyrannus morphology + PCA avg_Females_Adults.csv', row.names = 1)
#otu_cv_avg<-read.csv('./Output Files/cv_summary_Females_Adults.csv', row.names = 1)
#residuals<-read.csv('./Output Files/phylANOVA_tarsus-corrected_residuals_Females.csv', row.names = 1)
### OR
otu_avg<-read.csv('./Output Files/Tyrannus morphology + PCA avg_Males_Adults.csv', row.names = 1)
otu_cv_avg<-read.csv('./Output Files/cv_summary_Males_Adults.csv', row.names = 1)
residuals<-read.csv('./Output Files/phylANOVA_tarsus-corrected_residuals_Males.csv', row.names = 1)
### Note: If assessing sexes separately, use the following lines of code instead of the above 3 lines, then proceed.
#png(file="PPC scores boxplots_Females.png",width=7,height=2.75,units="in",res=500)
### OR
png(file="PPC scores boxplots_Males.png",width=7,height=2.75,units="in",res=500)
responsevariable.labs<-c("Bill PPC1", "Bill PPC2")
names(responsevariable.labs)<- c("BillPC1", "BillPC2")
col<-c("black","gray40","gray70","black","gray40","gray70")
p<-otu_avg %>%
pivot_longer(BillPC1:BillPC2, names_to="responsevariable", values_to = "value") %>%
ggplot(aes(y=value, x=Strategy)) +
geom_boxplot(color=col) +
facet_wrap(vars(responsevariable), ncol=2, labeller=labeller(responsevariable = responsevariable.labs), scales="free_y") +
panel_border() +
theme_bw() +
theme(panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
strip.text=element_text(size=15)) +
scale_y_continuous(expand=c(0.05,8)) +
labs(x="Migration Strategy", y= "PPC Scores")
draw_plot(plot, x = 0, y = 0, width = 1, height = 1)
ggdraw() +
draw_plot(p,0,0,1,1) +
draw_plot_label(c("A", "B"), c(0.085,0.565), c(0.99,0.99), size=15) + #main plot identifiers
draw_plot_label(c("A", "A", "A"), c(0.146, 0.287, 0.420), c(0.87, 0.87, 0.87), size=10) + #letters for plot A
draw_plot_label(c("A", "AB", "B"), c(0.625, 0.746, 0.9), c(0.87, 0.87, 0.87), size=10) #letters for plot B
dev.off()