library(geomorph)
library(Morpho)
library(ape)
library(phytools)
library(geiger)
library(ggplot2)
library(ggfortify)
library(geiger)
library(mvMORPH)
library(tibble)
library(PhylogeneticEM)
library(dplyr)
library(nlme)
library(ape)
library(RRphylo)
library(plot3D)

## Load phylogeny
tree<-read.nexus("Tree.nex")
tree2<-chronos(tree,calibration = makeChronosCalib(tree))
## Load DTA data
DTA<-read.csv("Data.csv")
DTA<-DTA[order(DTA$Species),]
rownames(DTA)<-sort(tree2$tip.label)
## Load GMM data
GMM<-read.morphologika("StrepLandmark_30.txt")
GPA<-gpagen(GMM$coords)
GMM<-geomorph.data.frame(coords=GPA$coords)
GMM$Csize<-GPA$Csize[match(rownames(DTA), names(GPA$Csize))]
## Create objects to process datasets 
DTA$CSize<-GMM$Csize
Diet<-as.factor(DTA$Diet)
names(Diet)<-sort(tree2$tip.label)
DNE<-DTA$DNE
names(DNE)<-sort(tree2$tip.label)
OPCR<-DTA$OPCR
names(OPCR)<-sort(tree2$tip.label)
RFI<-DTA$RFI
names(RFI)<-sort(tree2$tip.label)
GMM$BdM<-DTA$BdM
names(GMM$BdM)<-rownames(DTA)
GMM$BrM<-DTA$BrM
names(GMM$BrM)<-rownames(DTA)
GMM$ECV<-DTA$ECV
names(GMM$ECV)<-rownames(DTA)
## subset phylogeny per guild
treeIn<-keep.tip(tree2,c(names(subset(Diet,Diet=="Insectivory"))))
treeFo<-keep.tip(tree2,c(names(subset(Diet,Diet=="Folivory"))))
treeFr<-keep.tip(tree2,c(names(subset(Diet,Diet=="Frugivory"))))
## Extract allometry-corrected brain shape and size data
GMM$ACcoords<-procD.lm(coords~log(BdM),data = GMM,iter = 9999)$residuals
GMM$ACCSize<-procD.lm(log(Csize)~log(BdM),data = GMM,iter = 9999)$residuals
# Extract allometry and phylogeny-corrected brain size data
GMM$acpCsize<-procD.pgls(log(Csize)~log(BdM),phy = tree2,data = GMM,iter = 9999)$pgls.residuals
GMM$acpCsize<-GMM$acpCsize[,1]
## Extract factors
Diet<-as.factor(DTA$Diet)
names(Diet)<-sort(tree$tip.label)
Superfamily<-as.factor(DTA$Superfamily)
names(Superfamily)<-sort(tree2$tip.label)
## allometry-corrected pPCA
acGMMppca<-gm.prcomp(GMM$ACcoords,phy = tree2,GLS=TRUE)
acGMMppca<-as.data.frame(acGMMppca$x)
acGMMppca<-acGMMppca[order(rownames(acGMMppca)),]
acGMMppca<-add_column(acGMMppca,Diet=DTA$Diet)
## DTA pPCA
DTAppca<-phyl.pca(tree2,DTA[,4:6],method = "BM",mode = "corr")
DTAppca<-as.data.frame(DTAppca$S)
DTAppca<-add_column(DTAppca,Diet=DTA$Diet)
### Testing for Centroid size - BrM and ECV correspondence
summary(procD.lm(log(Csize)~log(ECV),data = GMM,iter = 9999))
summary(procD.lm(log(Csize)~log(BrM),data = GMM,iter = 9999))
### Testing for allometry
summary(procD.lm(coords~log(BdM),data = GMM,iter = 9999))
summary(procD.lm(log(Csize)~log(BdM),data = GMM,iter = 9999))
summary(procD.lm(log(DTA[,4:6])~log(BdM),data = GMM,iter = 9999))
### Phylogenetic signal
## acGMM
physignal(GMM$ACcoords,phy = tree2,iter = 9999)
## acCS
physignal(GMM$ACCSize,phy = tree2,iter = 9999)
## DTA
physignal(arrayspecs(log(DTA[,4:6]),1,3),phy = tree2,iter = 9999)
### Procrustes ANOVA for ecological signal (differences between guilds)
summary(procD.pgls(GMM$ACcoords~Diet,phy = tree2,iter = 9999))
summary(procD.lm(GMM$acpCsize~Diet,iter = 9999))
summary(procD.pgls(GMM$ACCSize~Diet,phy = tree2,iter = 9999)) ## Non-signifcant
summary(procD.pgls(log(DTA[,4:6])~Diet,phy = tree2,iter = 9999))
## Procrustes ANOVA for ecological signal (differences between guilds) without Aye-Aye
summary(procD.pgls(GMM$ACcoords[-3,]~Diet[-4],phy = drop.tip(tree2,"Daubentonia_madagascariensis"),iter = 9999))
summary(procD.lm(GMM$acpCsize[-4]~Diet[-4],iter = 9999))
summary(procD.pgls(log(DTA[-4,4:6])~Diet[-4],phy = drop.tip(tree2,"Daubentonia_madagascariensis"),iter = 9999))
### Morphological association
two.b.pls(acGMMppca[,1:4],GMM$acpCsize,iter = 9999)
two.b.pls(acGMMppca[,1:4],DTAppca[,1:3],iter = 9999)
two.b.pls(GMM$acpCsize,DTAppca[,1:3],iter = 9999)
### Per-guild evol rates comparison
compare.evol.rates(acGMMppca[match(tree2$tip.label,rownames(acGMMppca)),][,1:4],phy = tree2,gp=Diet,iter = 9999)
compare.evol.rates(GMM$acpCsize[match(tree2$tip.label,names(GMM$acpCsize))],phy = tree2,gp=Diet,iter = 9999)
compare.evol.rates(DTAppca[match(tree2$tip.label,rownames(DTAppca)),][,1:3],phy = tree2,gp=Diet,iter = 9999)
## Per-guild evol rate comparison without Aye-Aye
compare.evol.rates(acGMMppca[match(tree2$tip.label,rownames(acGMMppca)),][-4,1:4],phy = drop.tip(tree2,"Daubentonia_madagascariensis"),gp=Diet,iter = 9999)
compare.evol.rates(GMM$acpCsize[match(tree2$tip.label,names(GMM$acpCsize))][-4],phy = drop.tip(tree2,"Daubentonia_madagascariensis"),gp=Diet,iter = 9999)
compare.evol.rates(DTAppca[match(tree2$tip.label,rownames(DTAppca)),][-4,1:3],phy = drop.tip(tree2,"Daubentonia_madagascariensis"),gp=Diet,iter = 9999)
### RRphylo per-branch evol rates and tests
RRphylo_acGMMppca<-RRphylo(tree2,as.matrix(acGMMppca[,1:4]))
search.shift(RRphylo_acGMMppca,status.type = "sparse",state = Diet,nrep = 1000,filename = "RRphylo_DTAppca_Diet_shift")
RRphylo_acpCS<-RRphylo(tree2,GMM$acpCsize)
search.shift(RRphylo_acpCS,status.type = "sparse",state = Diet,nrep = 1000,filename = "RRphylo_DTAppca_Diet_shift")
RRphylo_DTAppca<-RRphylo(tree2,as.matrix(DTAppca[,1:3]))
search.shift(RRphylo_DTAppca,status.type = "sparse",state = Diet,nrep = 1000,filename = "RRphylo_DTAppca_Diet_shift")
### RRphylo per-branch evol rates and tests without Aye-Aye
RRphylo_acGMMppca2<-RRphylo(drop.tip(tree2,"Daubentonia_madagascariensis"),as.matrix(acGMMppca[,1:4]))
search.shift(RRphylo_acGMMppca,status.type = "sparse",state = Diet,nrep = 1000,filename = "RRphylo_DTAppca_Diet_shift")
RRphylo_acpCS2<-RRphylo(drop.tip(tree2,"Daubentonia_madagascariensis"),GMM$acpCsize)
search.shift(RRphylo_acpCS,status.type = "sparse",state = Diet,nrep = 1000,filename = "RRphylo_DTAppca_Diet_shift")
RRphylo_DTAppca2<-RRphylo(drop.tip(tree2,"Daubentonia_madagascariensis"),as.matrix(DTAppca[,1:3]))
search.shift(RRphylo_DTAppca,status.type = "sparse",state = Diet,nrep = 1000,filename = "RRphylo_DTAppca_Diet_shift")
### physignal in RRrates
physignal(RRphylo_acGMMppca$rates[-c(1:19),],tree2,iter = 9999)
physignal(RRphylo_acpCS$rates[-c(1:19),],tree2,iter = 9999)
physignal(RRphylo_DTAppca$rates[-c(1:19),],tree2,iter = 9999)
### linear regression of association between RRrates
summary(procD.pgls(RRphylo_acGMMppca$rates[-c(1:19),]~RRphylo_acpCS$rates[-c(1:19),],phy = tree2,iter = 9999))
summary(procD.pgls(RRphylo_acGMMppca$rates[-c(1:19),]~RRphylo_DTAppca$rates[-c(1:19),],phy = tree2,iter = 9999))
summary(procD.pgls(RRphylo_DTAppca$rates[-c(1:19),]~RRphylo_acpCS$rates[-c(1:19),],phy = tree2,iter = 9999))
### linear regression of association between RRrates without Aye-Aye
summary(procD.pgls(RRphylo_acGMMppca$rates[-c(1:19,28),]~RRphylo_acpCS$rates[-c(1:19,28),],phy = drop.tip(tree2,"Daubentonia_madagascariensis"),iter = 9999))
summary(procD.pgls(RRphylo_acGMMppca$rates[-c(1:19,28),]~RRphylo_DTAppca$rates[-c(1:19,28),],phy = drop.tip(tree2,"Daubentonia_madagascariensis"),iter = 9999))
summary(procD.pgls(RRphylo_DTAppca$rates[-c(1:19,28),]~RRphylo_acpCS$rates[-c(1:19,28),],phy = drop.tip(tree2,"Daubentonia_madagascariensis"),iter = 9999))
### Guild-wise evol modelling
## acpCSize
# Insectivores
mvBM(treeIn,subset(GMM$acpCsize,Diet=="Insectivory")[match(treeIn$tip.label,names(subset(GMM$acpCsize,Diet=="Insectivory")))],model = "BM1")
mvEB(treeIn,subset(GMM$acpCsize,Diet=="Insectivory")[match(treeIn$tip.label,names(subset(GMM$acpCsize,Diet=="Insectivory")))])
mvOU(treeIn,subset(GMM$acpCsize,Diet=="Insectivory")[match(treeIn$tip.label,names(subset(GMM$acpCsize,Diet=="Insectivory")))])
# Folivores
mvBM(treeFo,subset(GMM$acpCsize,Diet=="Folivory")[match(treeFo$tip.label,names(subset(GMM$acpCsize,Diet=="Folivory")))],model = "BM1")
mvEB(treeFo,subset(GMM$acpCsize,Diet=="Folivory")[match(treeFo$tip.label,names(subset(GMM$acpCsize,Diet=="Folivory")))])
mvOU(treeFo,subset(GMM$acpCsize,Diet=="Folivory")[match(treeFo$tip.label,names(subset(GMM$acpCsize,Diet=="Folivory")))])
# Frugivores
mvBM(treeFr,subset(GMM$acpCsize,Diet=="Frugivory")[match(treeFr$tip.label,names(subset(GMM$acpCsize,Diet=="Frugivory")))],model = "BM1")
mvEB(treeFr,subset(GMM$acpCsize,Diet=="Frugivory")[match(treeFr$tip.label,names(subset(GMM$acpCsize,Diet=="Frugivory")))])
mvOU(treeFr,subset(GMM$acpCsize,Diet=="Frugivory")[match(treeFr$tip.label,names(subset(GMM$acpCsize,Diet=="Frugivory")))])
## acGMMppca
# Insectivores
mvBM(treeIn,subset(acGMMppca[,1:4],Diet=="Insectivory")[match(treeIn$tip.label,rownames(subset(acGMMppca[,1:4],Diet=="Insectivory"))),],model = "BM1")
mvEB(treeIn,subset(acGMMppca[,1:4],Diet=="Insectivory")[match(treeIn$tip.label,rownames(subset(acGMMppca[,1:4],Diet=="Insectivory"))),])
mvOU(treeIn,subset(acGMMppca[,1:4],Diet=="Insectivory")[match(treeIn$tip.label,rownames(subset(acGMMppca[,1:4],Diet=="Insectivory"))),])
# Folivores
mvBM(treeFo,subset(acGMMppca[,1:4],Diet=="Folivory")[match(treeFo$tip.label,rownames(subset(acGMMppca[,1:4],Diet=="Folivory"))),],model = "BM1")
mvEB(treeFo,subset(acGMMppca[,1:4],Diet=="Folivory")[match(treeFo$tip.label,rownames(subset(acGMMppca[,1:4],Diet=="Folivory"))),])
mvOU(treeFo,subset(acGMMppca[,1:4],Diet=="Folivory")[match(treeFo$tip.label,rownames(subset(acGMMppca[,1:4],Diet=="Folivory"))),])
# Frugivores
mvBM(treeFr,subset(acGMMppca[,1:4],Diet=="Frugivory")[match(treeFr$tip.label,rownames(subset(acGMMppca[,1:4],Diet=="Frugivory"))),],model = "BM1")
mvEB(treeFr,subset(acGMMppca[,1:4],Diet=="Frugivory")[match(treeFr$tip.label,rownames(subset(acGMMppca[,1:4],Diet=="Frugivory"))),])
mvOU(treeFr,subset(acGMMppca[,1:4],Diet=="Frugivory")[match(treeFr$tip.label,rownames(subset(acGMMppca[,1:4],Diet=="Frugivory"))),])
## DTAppca
# Insectivores
mvBM(treeIn,subset(DTAppca[,1:3],Diet=="Insectivory")[match(treeIn$tip.label,rownames(subset(DTAppca[,1:3],Diet=="Insectivory"))),],model = "BM1")
mvEB(treeIn,subset(DTAppca[,1:3],Diet=="Insectivory")[match(treeIn$tip.label,rownames(subset(DTAppca[,1:3],Diet=="Insectivory"))),])
mvOU(treeIn,subset(DTAppca[,1:3],Diet=="Insectivory")[match(treeIn$tip.label,rownames(subset(DTAppca[,1:3],Diet=="Insectivory"))),])
# Folivores
mvBM(treeFo,subset(DTAppca[,1:3],Diet=="Folivory")[match(treeFo$tip.label,rownames(subset(DTAppca[,1:3],Diet=="Folivory"))),],model = "BM1")
mvEB(treeFo,subset(DTAppca[,1:3],Diet=="Folivory")[match(treeFo$tip.label,rownames(subset(DTAppca[,1:3],Diet=="Folivory"))),])
mvOU(treeFo,subset(DTAppca[,1:3],Diet=="Folivory")[match(treeFo$tip.label,rownames(subset(DTAppca[,1:3],Diet=="Folivory"))),])
# Frugivores
mvBM(treeFr,subset(DTAppca[,1:3],Diet=="Frugivory")[match(treeFr$tip.label,rownames(subset(DTAppca[,1:3],Diet=="Frugivory"))),],model = "BM1")
mvEB(treeFr,subset(DTAppca[,1:3],Diet=="Frugivory")[match(treeFr$tip.label,rownames(subset(DTAppca[,1:3],Diet=="Frugivory"))),])
mvOU(treeFr,subset(DTAppca[,1:3],Diet=="Frugivory")[match(treeFr$tip.label,rownames(subset(DTAppca[,1:3],Diet=="Frugivory"))),])
### phylo integration
summary(phylo.integration(RRphylo_acGMMppca$rates[-c(1:19),],RRphylo_acpCS$rates[-c(1:19),],phy = tree2,iter = 9999))
summary(phylo.integration(RRphylo_acGMMppca$rates[-c(1:19),],RRphylo_DTAppca$rates[-c(1:19),],phy = tree2,iter = 9999))
summary(phylo.integration(RRphylo_acpCS$rates[-c(1:19),],RRphylo_DTAppca$rates[-c(1:19),],phy = tree2,iter = 9999))
## phylo integration without Aye-Aye
summary(phylo.integration(RRphylo_acGMMppca$rates[-c(1:19,28),],RRphylo_acpCS$rates[-c(1:19,28),],phy = drop.tip(tree2,"Daubentonia_madagascariensis"),iter = 9999))
summary(phylo.integration(RRphylo_acGMMppca$rates[-c(1:19,28),],RRphylo_DTAppca$rates[-c(1:19,28),],phy = drop.tip(tree2,"Daubentonia_madagascariensis"),iter = 9999))
summary(phylo.integration(RRphylo_acpCS$rates[-c(1:19,28),],RRphylo_DTAppca$rates[-c(1:19,28),],phy = drop.tip(tree2,"Daubentonia_madagascariensis"),iter = 9999))
#### Plotting
# contMap with DTA and brain shape and size evol rates
plot(setMap(contMap(tree2,RRphylo_acpCS$rates[-c(1:19),],plot = FALSE),colors=c("yellow","red")))
plot(setMap(contMap(tree2,RRphylo_acGMMppca$rates[-c(1:19),],plot = FALSE),colors=c("yellow","red")))
plot(setMap(contMap(tree2,RRphylo_DTAppca$rates[-c(1:19),],plot = FALSE),colors=c("yellow","red")))
# Allometry-corrected brain shape phylomorphospace
shapes<-c(16, 17, 18)
shapes<-shapes[Diet]
scatter3D(acGMMppca[,1],acGMMppca[,2],acGMMppca[,3],colvar = NULL,col = Diet,xlab="pPC1 (20.73%)",ylab="pPC2 (15.69%)",zlab="pPC3 (10.85%)",bty = "b2",colkey = TRUE,pch=shapes[Diet],cex=1.5,type="h")