########################GMM Analysis

library(rgl)
options(rgl.useNULL = TRUE)
library(geomorph)
library(ggplot2)
library(gridExtra)
library(ape)
library(phangorn)
library(geiger)
library( abind )

#Set working directory to whatever folder you're working from. You'll need to modify the below
#Remember to put all the files you need like the Canine points.TPS and Sampled Diets.csv file to this certain folder
setwd( "~/Desktop/teeth_GMM_20221019" )

#Read the landmarks
landmark <- readland.tps(file = "CombinedTeeth1029.TPS", specID = "imageID", readcurves = F)

#define the tangent vectors describing how individual semilandmarks will slide along their curves / You could also define this when you are adding landmarks in tpsdig
sliders = rbind(define.sliders(c(1:8), write.file = F), 
                define.sliders(c(8:15), write.file = F))

#Do generalised Procrustes superposition, using bending energy as a criterion for sliding, and assign to an object called "GPA.fit"	
GPA.fit <- gpagen( landmark, curves = sliders, ProcD = F )

#get the diet and the species names information, and organize the taxon names information
dietcat <- read.csv(file="Table_S1_GMM.csv", header=T, row.names=NULL)
diet <- dietcat$Primary_diet
dietcat$Species_name <-gsub(" ", "_", dietcat$Species_name) #replace the space with _
dimnames(GPA.fit$coords)[[3]] <- names( GPA.fit$Csize ) <- dietcat$Species_name #replace the specimen names in the GMM results with the proper species names
dimnames(GPA.fit$coords)[[3]] #124 names here including the duplicate species
names( GPA.fit$Csize )


#########################################################phylogenetic ANOVA for GMM results

######prepare the GMM data, take mean values for multiple individuals for some species based on the species name
species.list <- unique( dietcat$Species_name )	#delete the duplicate species names in GMM, get the non-duplicate specimen name list
species.list #should become 120 names here, with duplicate names deleted 

#take mean values of the GMM coordinates by writing a function, it will return a list AA
AA <- lapply( X = species.list , FUN = function( GPA.coords , species_name ){ apply( GPA.coords[,, dimnames( GPA.coords )[[ 3 ]] == species_name ] , c(1,2) , mean ) } , GPA.coords = GPA.fit$coords )
newgpa <- do.call( abind , c( AA , along = 3 ) ) #convert this list to an array
dimnames( newgpa )[[ 3 ]] <- species.list	#give names to the array
str(newgpa) #new gpa results with 120 samples

#take mean values of the centroid sizes
Csize <- aggregate( GPA.fit$Csize , by = list( names( GPA.fit$Csize ) ) , FUN = mean )$x #get the values
names( Csize ) <- unique( names( GPA.fit$Csize ) ) #name
str(Csize) #120 samples


############## prepare the diet data
head(dietcat)
dietcat1 <- dietcat[match( dimnames( newgpa )[[ 3 ]] , dietcat$Species_name),] #delete the duplicate samples by subset the species names left in the newgpa
rownames( dietcat1 ) <- dietcat1$Species_name
str(dietcat1) #120 samples in the new diet data

############### prepare the tree
##read the stitched tree of all the extant samples
tree1 <- read.tree("tree_final.nex") 
plot( tree1 , direction = "u" , no.margin = T , cex = 0.6 )

#remove the species in the tree but not in the diet data
tree2 <-drop.tip(tree1, tree1$tip.label[ ! tree1$tip.label %in% dietcat1$Species_name ]) #drop the tips

#remove the species in the diet data but not in the tree, and returns the diet data with rows in the same order as they appear in the tree
dietcat2 <- dietcat1[ tree2$tip.label , ]
head(dietcat2)

#check the numbers and names of species are correct between the data and the tree
str(tree2) #103 tips
str(dietcat2) #103 specimens
all.equal( tree2$tip.label , rownames( dietcat2 ) )		# check if names are the same and in the same order in the tree and the diet data

###############stick the data to the tree, using newgpa [,,tree2$tip.label] and Csize[ tree2$tip.label] to make sure the order in the coordinates and Csize is the same as in the tree and the diet data
gdf <- geomorph.data.frame( phy = tree2 , 
                            carnivorous = dietcat2$Primary_diet == "Carnivorous" , 
                            piscivorous = dietcat2$Primary_diet == "Piscivorous" , 
                            insectivorous = dietcat2$Primary_diet == "Insectivorous" , 
                            omnivorous = dietcat2$Primary_diet == "Omnivorous" , 
                            frugi_herbivorous = dietcat2$Primary_diet == "Frugi/Herbivorous" , 
                            shape = newgpa[ , , tree2$tip.label ] , 
                            Csize = log10( Csize[ tree2$tip.label ]),
                            bodymass = log10(as.numeric(dietcat2$Body_mass_mean_g))) #log transform for Csize and body mass

##############perform the analyses
models <- c("Csize + bodymass + carnivorous" , "Csize + bodymass + piscivorous" , "Csize + bodymass + insectivorous" , 
            "Csize + bodymass + omnivorous" , "Csize + bodymass + frugi_herbivorous")
models1 <- lapply( paste ( "shape ~" , models ) , as.formula )
models1

#run phylogenetic ANOVA for different diet groups with considering of the centroid sizes
pgls.fit.list <- lapply( models1 , procD.pgls , phy = phy , data = gdf , iter = 999 ) 
pgls.summary <- lapply( pgls.fit.list , summary ) 
pgls.summary
capture.output(pgls.summary, file = "GMM pgls results.txt") #export results

#run phylogenetic ANOVA for different diet groups with considering of the centroid sizes
ols.fit.list <- lapply( models1 , procD.lm , data = gdf , iter = 999 )
ols.summary <- lapply( ols.fit.list , summary )
ols.summary
capture.output(ols.summary, file = "GMM ols results.txt") #export results

##############################################	 



##########################################################Phylogenetic ANOVA for bird bill proportion

######## prepare the data
birdbill1<- read.csv(file="Table S4 Birds (10.28).csv", header=T, row.names=NULL)
str(birdbill1) #185 specimens in the data here
birdbill1$Species_name <-gsub(" ", "_", birdbill1$Species_name) #replace the space with _
head(birdbill1)

#take mean values of "Skull_length_mm" and "Rostral_length_mm" separately based on the species names
birdbill2 <- aggregate ( birdbill1[ , c("Skull_length_mm", "Rostral_length_mm") ] , by = list( birdbill1$Species_name ) , FUN = "mean")
rownames( birdbill2 ) <- birdbill2$Group.1 #make sure the species names = row names

#combine the other columns from the original data
birdbill2 <- cbind ( birdbill2 , birdbill1[ match( rownames( birdbill2 )  , birdbill1$Species_name ) , c( "Family" , "Primary_diet" , "Body_mass_mean_g", "Skull.Rostrum_ratio")] )
str(birdbill2) #should be 352 specimens here after aggregate
head(birdbill2)

####### prepare the tree
tree_bird1 <- read.nexus( "BirdBill10_28.nex" )	#Read bird tree subset from vertlife.org
tree_bird2<-maxCladeCred(tree_bird1, tree = TRUE, part = NULL, rooted = TRUE) #create a maximum clade credibility tree based on these downloaded trees
write.tree(tree_bird2,file='tree_bird.nex') #save it as a new .nex tree file
str(tree_bird2) #347 tips

#remove the species in the tree but not in the data from the tree
tree_bird3 <-drop.tip(tree_bird2, tree_bird2$tip.label[ ! tree_bird2$tip.label %in% rownames(birdbill2) ]) #drop the tips

#remove the species in the data but not in the tree from the data, and reorder the data to be in the same order as in the tree
birdbill3 <- birdbill2[ tree_bird3$tip.label, ]

#check the numbers and names of species are correct between the data and the tree
str(birdbill3) #347 specimens
str(tree_bird3) #347 tips
all.equal(tree_bird3$tip.label, rownames(birdbill3)) #check if the names and orders are the same in the tree and in the data


########### perform the analyses
SL <- log10( birdbill3$Skull_length_mm )
RL <- log10( birdbill3$Rostral_length_mm )
DT <- birdbill3$Primary_diet
BM <- log10(as.numeric(birdbill3$Body_mass_mean_g))
names( SL ) <- names( RL ) <- names( DT ) <- names (BM) <- rownames(birdbill3)

#stick the data to the tree
gdf.bill <- geomorph.data.frame( phy = tree_bird3 , skull_length = SL , rostral_length = RL , diet = DT , bodymass = BM)
str(gdf.bill)


#run phylogenetic ANOVA
birdbill.pgls <- procD.pgls( rostral_length ~ skull_length + bodymass + diet , phy = phy, data = gdf.bill, iter = 999)
summary(birdbill.pgls) 
capture.output(summary(birdbill.pgls) , file = "Bird bill pgls results.txt") #export results

#run ordinary ANOVA
birdbill.ols <- procD.lm( rostral_length ~ skull_length + bodymass + diet , data = gdf.bill, iter = 999)
summary(birdbill.ols) 
capture.output(summary(birdbill.ols) , file = "Bird bill ols results.txt") #export results

###########################################################################################	 









######################################## PCA and results plotting

#Do principal components analysis (PCA) and assign to an object called "PCA.fit"
PCA.fit <- gm.prcomp( GPA.fit$coords )

#The most basic way to plot PCA results
plot(PCA.fit,1,2)
plot(PCA.fit,2,3)

#Read PCA data and add diet category data to it
write.csv(file="PCA_teeth.csv",PCA.fit$x)
tt<-read.csv(file="PCA_teeth.csv",header=T)
PCA.diet<-data.frame(tt,diet)
PCA.diet

#PCA results plot
#Generate labels for the x axis and y axis of PCA ordination plots, specifying which axes to plot. 
#Here we specify PC1 on the x axis, PC2 on the y axis of the first plot, and PC3 on the y axis of the second plot. 
#PCA.fit$d[*]/sum(PCA.fit$d) indicates the variance, representing how far the variable is spread out from their average value.
X <- "PC1"
XLAB <- paste( X , " (" , round( PCA.fit$d[1]/sum(PCA.fit$d)*100, 2) , "%)" , sep = "" )
Y1 <- "PC2"
Y1LAB <- paste( Y1 , " (" , round( PCA.fit$d[2]/sum(PCA.fit$d)*100, 2) , "%)" , sep = "" )
Y2 <- "PC3"
Y2LAB <- paste( Y2 , " (" , round( PCA.fit$d[3]/sum(PCA.fit$d)*100 , 2) , "%)" , sep = "" )				

# plot the results with convex hulls grouped by diets
library(plyr)

##Plots of PC1-PC2 and PC1-PC3
hulls1<-ddply(PCA.diet,.(diet),function(PCA.diet)PCA.diet[chull(PCA.diet$Comp1, PCA.diet$Comp2),])

p1<-ggplot(PCA.diet,aes(x = Comp1,y= Comp2,group = diet))+
  geom_point(aes(shape = diet, color = diet),size=3)+
  geom_polygon(data = hulls1, aes(x = Comp1,y= Comp2, fill = diet, group = diet),alpha = 0.18)+
  theme(legend.position = "bottom")+
  theme(panel.background = element_rect(fill = 'white', colour = 'black'))+
  labs(x = XLAB, y = Y1LAB)

customshape<-p1 + scale_shape_manual(values = c(19, 19, 19, 19, 19, 8))
customized1 <- customshape + scale_color_manual(values = c("#da9f3a", "#eee262", "#4a9a75", "#6fb2e3", "#be7da3", "black"))


hulls2<-ddply(PCA.diet,.(diet),function(PCA.diet)PCA.diet[chull(PCA.diet$Comp1, PCA.diet$Comp3),])

p2<-ggplot(PCA.diet,aes(x = Comp1,y= Comp3,group = diet))+
  geom_point(aes(shape = diet, color = diet),size=3)+
  geom_polygon(data = hulls2, aes(x = Comp1,y= Comp3, fill = diet, group = diet),alpha = 0.18)+
  theme(legend.position = "bottom")+
  theme(panel.background = element_rect(fill = 'white', colour = 'black'))+
  labs(x = XLAB, y = Y2LAB)

customshape2<-p2 + scale_shape_manual(values = c(19, 19, 19, 19, 19, 8))
customized2 <- customshape2 + scale_color_manual(values = c("orange", "wheat3", "palegreen3", "skyblue3", "orchid3", "red"))

X=arrangeGrob(customized1, customized2, ncol=2)

ggsave(filename="Hulls_by_Diet_Result_1130.png", width=16, height=7, X)


###############################The codes below are used to plot shapes of different extreme PC values with reference to mean shape, so you know what those PC values mean for the shape, and sometimes you may want to add these to your final figure
#Define the link between landmarks and semilandmarks
link<-matrix(c(1, rep(1:14,each =2),15), ncol=2, byrow=T)

#Get the mean shape as the reference
mean<-GPA.fit$consensus

##PC1max
plotRefToTarget(mean,PCA.fit$shapes$shapes.comp1$max,gridPars=gridPar(pt.bg = "grey", pt.size = 1,link.lwd = 0.3,
                                                                      tar.pt.bg = "red", tar.pt.size = 1, tar.link.col = "red", tar.link.lwd = 0.3), 
                method="points", links = link)


##PC1min
plotRefToTarget(mean,PCA.fit$shapes$shapes.comp1$min,gridPars=gridPar(pt.bg = "grey", pt.size = 1,link.lwd = 0.3,
                                                                      tar.pt.bg = "red", tar.pt.size = 1, tar.link.col = "red", tar.link.lwd = 0.3), 
                method="points", links = link)

##PC2max
plotRefToTarget(mean,PCA.fit$shapes$shapes.comp2$max,gridPars=gridPar(pt.bg = "grey", pt.size = 1,link.lwd = 0.3,
                                                                      tar.pt.bg = "red", tar.pt.size = 1, tar.link.col = "red", tar.link.lwd = 0.3), 
                method="points", links = link)

##PC2min
plotRefToTarget(mean,PCA.fit$shapes$shapes.comp2$min,gridPars=gridPar(pt.bg = "grey", pt.size = 1,link.lwd = 0.3,
                                                                      tar.pt.bg = "red", tar.pt.size = 1, tar.link.col = "red", tar.link.lwd = 0.3), 
                method="points", links = link)

##PC3max
plotRefToTarget(mean,PCA.fit$shapes$shapes.comp3$max,gridPars=gridPar(pt.bg = "grey", pt.size = 1,link.lwd = 0.3,
                                                                      tar.pt.bg = "red", tar.pt.size = 1, tar.link.col = "red", tar.link.lwd = 0.3), 
                method="points", links = link)

##PC3min
plotRefToTarget(mean,PCA.fit$shapes$shapes.comp3$min,gridPars=gridPar(pt.bg = "grey", pt.size = 1,link.lwd = 0.3,
                                                                      tar.pt.bg = "red", tar.pt.size = 1, tar.link.col = "red", tar.link.lwd = 0.3), 
                method="points", links = link)


##check the shape of certain taxon you are interested like the 1st one in the list here
#I'm currently looking at Longipteryx chaoyangenis
#Han: I replaced the number by the species name now to make sure it is the correct specimen, it's different, Alex you might be using the wrong number before
plotRefToTarget(mean,GPA.fit$coords[,,"Longipteryx_chaoyangensis"],gridPars=gridPar(pt.bg = "blue", pt.size = 1,link.lwd = 1,
                                                                                    tar.pt.bg = "red", tar.pt.size = 0.8, tar.link.col = "red", tar.link.lwd = 0.8), 
                method="points", links = link,label=T)
#Rostral Proportions

BillData <- read.csv("BirdBillDataSet.csv")
head(BillData)
Billscatter <- ggplot(BillData, aes(x = RL, y = SRr, shape = Family, color = Family,
                                    ggtitle("Rostral Length and Proportions Compared"))) + geom_point(size = 2)
BillscatterTitle <- Billscatter + xlab("Rostral Length (mm)") + ylab("Rostral Proportion of Total Skull")
Billscattershape <- BillscatterTitle + scale_shape_manual(values = c(19, 19, 19, 19, 8, 19, 19, 19, 19, 19, 19, 19))
Billscattershape + scale_color_brewer(palette = "Paired")
Familyfinal <- Billscattershape + scale_color_brewer(palette = "Paired")
B = Billscattercolor
ggsave(filename="BirdBillbyFamily.png", width=10, height=6, B)
#Color by Diet 
BillData <- read.csv("BirdBillDataSet.csv")
head(BillData)
BillFamily <- ggplot(BillData, aes(x = RL, y = SRr, color = Diet, shape = Diet,
                                   ggtitle("Rostral Length and Proportions Compared"))) + geom_point(size = 2)
BillFamilyTitle <- BillFamily + xlab("Rostral Length (mm)") + ylab("Rostral Proportion of Total Skull")
BillFamilyTitle + scale_color_brewer(palette = "Paired")
Dietfinal <- BillFamilyTitle + scale_color_brewer(palette = "Paired")
Dietfinal + scale_shape_manual(values = c(19, 19, 19, 8, 19, 19))
Dietfinal2 <- Dietfinal + scale_shape_manual(values = c(19, 19, 19, 8, 19, 19))
BIRDS = grid.arrange(Familyfinal, Dietfinal2, ncol=2)
ggsave(filename="BirdBillsCompared.png", width=16, height=6, BIRDS)


#Chiropteran Canine Proportions
library("ggplot2")
library(ggplot2)
library(ochRe)
+ scale_color_ochre(palette="healthy_reef")

Batcanines <- read.csv("BatOnlyCanines.csv")
head(Batcanines)
Batcaninesscatter <- ggplot(Batcanines, aes(x = Skull, y = WidthHeight, color = Diet, 
                                            ggtitle("Example Title"))) + geom_point(size = 3)
BatcaninesscatterTitle <- Batcaninesscatter + xlab("Skull Length(mm)") + ylab("Canine Width:Height")
BatcaninesscatterTitle
B = BatcaninesscatterTitle 
ggsave(filename="BatCaninesWidthandLength.png", width=8, height=6, B)
Canines <- read.csv("BatOnlyCanines.csv")
head(Canines)
#Canine Width as X
plot <- ggplot(Canines, aes(x = W, y = H, color = Diet)) + geom_point(size=5, shape=20)
plot2 <- plot + xlab("Canine Width(mm)") + ylab("Canine Height(mm)")   
Batwidth <- plot2 + scale_color_ochre(palette="parliament")
#Skull Length as X
plot3 <- ggplot(Canines4, aes(x = SL, y = CH, color = Diet)) + geom_point(size=5, shape=20)
plot4 <- plot3 + xlab("Skull Length(mm)") + ylab("Canine Height(mm)")   
Batskull <- plot4 + scale_color_ochre(palette="parliament")
BATS = grid.arrange(Batskull, Batwidth, ncol=2)
ggsave(filename="BatCanineResults11_12.png", width=16, height=6, BATS)
#Random Ternary Plots
Canines <- read.csv("Canine test.csv")
head(Canines)
ggtern(data = Canines, aes(x = H, y = W, z = H.W.)) + 
  geom_point(aes(fill = Group), shape = 21, size = 2) +
  tern_limit(T = .7, L = 1, R = .5)
Canines2 <- read.csv("Canine test 2.csv")
head(Canines2)
ggtern(data = Canines2, aes(x = CH.SL, y = CW.CH, z = CW.SL)) + 
  geom_point(aes(fill = Group), shape = 21, size = 2) +
  tern_limit(T = 1, L = .7, R = .68)
Canines3 <- read.csv("Canine test 3.csv")
head(Canines3)
ggtern(data = Canines3, aes(x = CH.SL, y = CW.CH, z = CW.SL)) + 
  geom_point(aes(fill = Group), shape = 21, size = 2) +
  tern_limit(T = .95, L = .68, R = .58)
Canines4 <- read.csv("Canine test 4.csv")
head(Canines4)
ggtern(data = Canines4, aes(x = CH, y = CW, z = SL)) + 
  geom_point(aes(fill = Group), shape = 21, size = 2) +
  tern_limit(T = 1, L = 1, R = 1)

######################Bird Bill pGLS 

require( ape )
require( caper )
require( mvMORPH )
require( abind )
require (qpcR)

setwd( "~/Desktop/Recent items/R functions Roger Benson/Alex Clark longipterygid analyses 01Nov2022")

data <- read.csv( "Table S4 Birds (10.28).csv" )
trees <- read.nexus( "BirdBill10_28.nex" )

trees[[1]]$tip.label %in% gsub( " " , "_" , data$Species_name )
trees[[1]]$tip.label %in% gsub( " " , "_" , data$Name_in_data )

tree.data <- aggregate( data[ c( "Skull_length_mm" , "Rostral_length_mm" ) ] , mean , by = list( gsub( " " , "_" , data$Species_name ) ) )
rownames( tree.data ) <- tree.data[ , "Group.1" ]
tree.data[ , "Skull_length_mm" ] <- log10( tree.data[ , "Skull_length_mm" ] )
tree.data[ , "Rostral_length_mm" ] <- log10( tree.data[ , "Rostral_length_mm" ] )

#B		
tree.data$piscivory <- data[ match( rownames( tree.data ) , gsub( " " , "_" , data$Species_name ) ) , "Primary_diet" ] == "Piscivory"
tree.data$carnivory <- data[ match( rownames( tree.data ) , gsub( " " , "_" , data$Species_name ) ) , "Primary_diet" ] == "Carnivory"
tree.data$insectivory <- data[ match( rownames( tree.data ) , gsub( " " , "_" , data$Species_name ) ) , "Primary_diet" ] == "Insectivory"							
tree.data$omnivory <- data[ match( rownames( tree.data ) , gsub( " " , "_" , data$Species_name ) ) , "Primary_diet" ] == "Omnivory"
tree.data$nectivory <- data[ match( rownames( tree.data ) , gsub( " " , "_" , data$Species_name ) ) , "Primary_diet" ] == "Nectivory"
tree.data$herbivory <- data[ match( rownames( tree.data ) , gsub( " " , "_" , data$Species_name ) ) , "Primary_diet" ] == "Herbivory"
tree.data$invertivory <- data[ match( rownames( tree.data ) , gsub( " " , "_" , data$Species_name ) ) , "Primary_diet" ] == "Invertivory"
tree.data[ tree.data$invertivory , ]
tree.data[ tree.data$herbivory , ]


i <- 1
data.temp <- comparative.data( trees[[i]] , tree.data[ trees[[i]]$tip.label , ] , names.col = "Group.1" )
rostrum.residual.length.pgls <- pgls( Rostral_length_mm ~ Skull_length_mm , data = data.temp , lambda = "ML" )
summary( rostrum.residual.length.pgls )
rostrum.residual.length <- residuals( rostrum.residual.length.pgls )
write.csv(rostrum.residual.length, "pgls rostrum for all birds")
#C
models1 <- c( "1" , "Skull_length_mm" , "Skull_length_mm + piscivory", "Skull_length_mm + carnivory" , "Skull_length_mm + insectivory" , "Skull_length_mm + omnivory" , "Skull_length_mm + herbivory" , "Skull_length_mm + invertivory"  )
calls1 <- paste( "Rostral_length_mm ~ " , models1 , sep = "" )
formulae1 <- lapply( calls1 , as.formula )

pgls.fits1 <- lapply( formulae1 , pgls , data = data.temp , lambda = "ML" )
names( pgls.fits1 ) <- models1
lapply( pgls.fits1 , summary )

summaries1 <- lapply( pgls.fits1 , function(X) { c( summary(X)$adj.r.squared , round( AICc( X ) , 2 ) , X$param.CI$lambda$opt ) } )
anova.tables1 <- lapply( pgls.fits1 , function(X) { AA <- summary( X )$coef ; AA <- cbind( rownames( AA ) , round( AA , 4 ) ) } )
for( i in 1:length( summaries1 ) ) {
  anova.tables1[[ i ]] <- cbind( matrix( ncol = 1 + length( summaries1[[ i ]] ) , nrow = nrow ( anova.tables1[[ i ]] ) ) , anova.tables1[[ i ]] )
  anova.tables1[[ i ]][ 1 , 1 ] <- calls1[ i ]
  anova.tables1[[ i ]][ 1 , 2:(1+length( summaries1[[ i ]] )) ] <- round( summaries1[[ i ]] , 4 )
  colnames( anova.tables1[[ i ]] )[1:5] <- c( "Model" , "R2" , "AICc" , "lambda" , "Variable" )
}

write.csv( do.call( rbind , anova.tables1 ) , "Longipteryx rostral length pgls results 02Nov2022.csv")

#D			
models2 <- c( "1" , "piscivory", "carnivory" , "insectivory" , "omnivory" , "herbivory" , "invertivory"  )
calls2 <- paste( "Skull_length_mm ~ " , models2 , sep = "" )
formulae2 <- lapply( calls2 , as.formula )

pgls.fits2 <- lapply( formulae2 , pgls , data = data.temp , lambda = "ML" )
names( pgls.fits2 ) <- models
lapply( pgls.fits2 , summary )					

summaries2 <- lapply( pgls.fits2 , function(X) { c( summary(X)$adj.r.squared , round( AICc( X ) , 2 ) , X$param.CI$lambda$opt ) } )
anova.tables2 <- lapply( pgls.fits2 , function(X) { AA <- summary( X )$coef ; AA <- cbind( rownames( AA ) , round( AA , 4 ) ) } )
for( i in 1:length( summaries2 ) ) {
  anova.tables2[[ i ]] <- cbind( matrix( ncol = 1 + length( summaries2[[ i ]] ) , nrow = nrow ( anova.tables2[[ i ]] ) ) , anova.tables2[[ i ]] )
  anova.tables2[[ i ]][ 1 , 1 ] <- calls2[ i ]
  anova.tables2[[ i ]][ 1 , 2:(1+length( summaries2[[ i ]] )) ] <- round( summaries2[[ i ]] , 4 )
  colnames( anova.tables2[[ i ]] )[1:5] <- c( "Model" , "R2" , "AICc" , "lambda" , "Variable" )
}

write.csv( do.call( rbind , anova.tables2 ) , "Longipteryx skull length (size) pgls results 02Nov2022.csv")


shape <- cbind( tree.data[ trees[[i]]$tip.label , "Skull_length_mm" ], rostrum.residual.length )


models3 <- c( "piscivory", "carnivory" , "insectivory" , "omnivory" , "herbivory" , "invertivory"  )
calls3 <- paste( "shape ~ " , models3 , sep = "" )
formulae <- lapply( calls3 , as.formula )

mvgls.fits <- lapply( formulae3 , mvgls , data = tree.data[ trees[[i]]$tip.label , ] , tree =  trees[[i]] , model = "lambda" , method = "LOO" )				
names( mvgls.fits ) <- models
mvgls.sig <- lapply( mvgls.fits , manova.gls , nperm = 100 , test = "Pillai" , verbose = TRUE )

mvgls.results.tables <- list()

for( j in 1:length( mvgls.fits ) ) {
  AA <- cbind( matrix( ncol = 5, nrow = nrow( mvgls.fits[[j]]$coefficients ) ) , round( mvgls.fits[[j]]$coefficients , 4 ) )
  colnames( AA ) <- c( "Model" , "GIC" , "Lambda" , "Test.statistic" , "p (model)" , "Coeff (log10 skull length)" , "Coeff (residual rostral length)" )
  
  AA[1,1:5] <- unlist( c( calls3[ j ] , round( unlist( GIC( mvgls.fits[[j]] )[2] ), 4 ) , round( mvgls.fits[[j]]$opt$par[2] , 4 ) , round( unlist( mvgls.sig[[j]][ c( "stat" ,"pvalue" ) ] ) , 4 ) ) )
  mvgls.results.tables[[ j ]] <- AA		
}

write.csv( do.call( rbind , mvgls.results.tables ) , "Longipteryx multivariate pgls results 02Nov2022.csv")

list( do.call( rbind , anova.tables1 ) , do.call( rbind , anova.tables2 ) , do.call( rbind , mvgls.results.tables ) )


############Scatter Plot
BillData <- read.csv("Bill_Scatter_Data(11.02).csv")
head(BillData)

Billscatter <- ggplot(BillData, aes(x = Log_skull_length, y = pgls_log_results, color = Primary_diet, shape = Primary_diet)) + geom_point(size = 3)

BillscatterTitle <- Billscatter + xlab("Log Skull Length") + ylab("pGLS Regression Residuals")

BillscatterTitlecolor <- BillscatterTitle + scale_color_manual(values =c("#a73636", "#51984c", "#888d7b", "#a6cee3","black", "75519c", "#bd90d3", "#1f78b4", ""))

Bill_All <- BillscatterTitlecolor + scale_shape_manual(values = c(19, 19, 19, 19, 8, 19, 19, 19))

ggsave(filename="Updated_11.02.png", width=10, height=6, Bill_All)


#Enamel Plot
Enamel <- read.csv("file.csv")
head(Enamel)

Enamelplot <- ggplot(Enamel, aes(x = Log_Mass, y = Log_Enamel, color = Diet, shape = Group, 
                                 ggtitle("Rostral:Skull Compared"))) + geom_point(size = 6)
Enamelplot2 <-Enamelplot + scale_shape_manual(values = c(8, 7, 19, 18, 17))

Enamelplot3 <- Enamelplot2 + scale_color_manual(values =c("black", "#009e73", "#cc79a7", "#d55e00", "#cb3939", "#9b9c3e", "0072b2"))