RStudio.Version() #Citate RStudio
citation() #Citate R
#HOW TO CITATE PACKAGES:
#The minimal requirement is to cite the R package in text along with the version number.
#Additionally, you can include the reference list entry the authors of the remote package have suggested.
#Example of sentence: "Analysis of the data was done using the remote package (v1.2.1; Appelhans et al., 2015)."
citation("phytools") #example
packageVersion("phytools") #example

rm(list=ls()) 
library(readr)
library(ggplot2)
library(multcompView)
library(dplyr)
library(readxl)

table <- read_excel("Badia_log_data.xlsx", 
                    col_types = c("text", "text", "text", 
                                  "text", "text", "text", "text", "numeric", 
                                  "numeric", "numeric", "numeric", 
                                  "numeric", "numeric", "numeric", 
                                  "numeric", "numeric", "numeric", 
                                  "numeric", "numeric", "numeric", 
                                  "numeric", "numeric", "numeric", 
                                  "text", "numeric", "numeric", "numeric", 
                                  "numeric", "numeric","numeric", "numeric", "text", "text", 
                                  "text", "text"))
table<-table[c(1:87),c(1:4,8:22,25:35)] #Run without the TL to avoid NAs
table<-na.omit(table)

#Read and match the tree:
library(phytools)
tree<-read.tree("ZhengWiensAgamidae.tree")
plot(tree, cex=0.12, label.offset=0.5, lwd=0.01)
is.ultrametric(tree)
tree<-force.ultrametric(tree)
is.binary(tree)
library(geiger)
str(table)
str(tree) 

row.names(table)<-table$TreeName
cut<-treedata(tree,table,sort=TRUE)
str(cut)
cuttree<-cut$phy
is.binary(cuttree)
is.ultrametric(cuttree)
cuttable<-cut$data
plot(cuttree,cex=0.45)
View(cuttable) #A. africana and A. boensis are not in the cuttree, this is the only difference


#Converting an atomic vector into a data frame: S-APPLY FUNCTION
cuttable<-as.data.frame(cuttable)
row.names(cuttable)<-cuttable$TreeName
cuttable[,c(5:26)] <- sapply(cuttable[,c(5:26)],function(x) as.numeric(x))


#To prove for phylogenetic signal so that I can justify why I get the phylogenetic residues:
#Pagel lambda:
library(picante)
library(ape)
library(adephylo)
library(ade4)
library(phylobase)
library(geiger)
library(phytools) ###I think all these packages are not necessary
phylosig(cuttree, cuttable$SVL, method="lambda", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$AHL, method="lambda", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$HL, method="lambda", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$HW, method="lambda", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$HH, method="lambda", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$ED, method="lambda", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$UFL, method="lambda", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$LFL, method="lambda", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$SFtoeL, method="lambda", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$LFtoeL, method="lambda", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$UHL, method="lambda", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$LHL, method="lambda", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$UHW, method="lambda", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$SHtoeL, method="lambda", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$LHtoeL, method="lambda", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$FLL, method="lambda", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$HLL, method="lambda", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$FLLvsHLL, method="lambda", test=TRUE, nsim=1000)

#K:
phylosig(cuttree, cuttable$SVL, method="K", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$AHL, method="K", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$HL, method="K", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$HW, method="K", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$HH, method="K", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$ED, method="K", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$UFL, method="K", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$LFL, method="K", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$SFtoeL, method="K", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$LFtoeL, method="K", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$UHL, method="K", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$LHL, method="K", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$UHW, method="K", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$SHtoeL, method="K", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$LHtoeL, method="K", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$FLL, method="K", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$HLL, method="K", test=TRUE, nsim=1000)
phylosig(cuttree, cuttable$FLLvsHLL, method="K", test=TRUE, nsim=1000)


#And also I have to prove that there is a significant correlation betwen SVL and the different variable to justify the phylogenetic size correction (phyl residues corrected for SVL that I am getting):
HLSVL <- lm(AHL ~ SVL, data = cuttable)
summary(HLSVL)
HLSVL <- lm(HL ~ SVL, data = cuttable)
summary(HLSVL)
HWSVL <- lm(HW ~ SVL, data = cuttable)
summary(HWSVL)
HHSVL <- lm(HH ~ SVL, data = cuttable)
summary(HHSVL)
EDSVL <- lm(ED ~ SVL, data = cuttable)
summary(EDSVL)
UFLSVL <- lm(UFL ~ SVL, data = cuttable)
summary(UFLSVL)
LFLSVL <- lm(LFL ~ SVL, data = cuttable)
summary(LFLSVL)
SFtoeLSVL <- lm(SFtoeL ~ SVL, data = cuttable)
summary(SFtoeLSVL)
LFtoeLSVL <- lm(LFtoeL ~ SVL, data = cuttable)
summary(LFtoeLSVL)
UHLSVL <- lm(UHL ~ SVL, data = cuttable)
summary(UHLSVL)
LHLSVL <- lm(LHL ~ SVL, data = cuttable)
summary(LHLSVL)
UHWSVL <- lm(UHW ~ SVL, data = cuttable)
summary(UHWSVL)
SHtoeLSVL <- lm(SHtoeL ~ SVL, data = cuttable)
summary(SHtoeLSVL)
LHtoeLSVL <- lm(LHtoeL ~ SVL, data = cuttable)
summary(LHtoeLSVL)
FLLSVL <- lm(FLL ~ SVL, data = cuttable)
summary(FLLSVL)
HLLSVL <- lm(HLL ~ SVL, data = cuttable)
summary(HLLSVL)
FLLvsHLLSVL <- lm(FLLvsHLL ~ SVL, data = cuttable)
summary(FLLvsHLLSVL) #This variable is not correlated, already knew so
FLtoesSVL <- lm(FLtoes ~ SVL, data = cuttable)
summary(FLtoesSVL)
HLtoesSVL <- lm(HLtoes ~ SVL, data = cuttable)
summary(HLtoesSVL)
LtoevsFLLSVL <- lm(LtoevsFLL ~ SVL, data = cuttable)
summary(LtoevsFLLSVL)
LtoevsHLLSVL <- lm(LtoevsHLL ~ SVL, data = cuttable)
summary(LtoevsHLLSVL)

#Get the residuals:
temp<-phyl.resid(cuttree,cuttable$SVL,cuttable[,c(6:19,24:25)],method = "BM") #I didn't put fingR, toeR nor FLLvsHLL as I musn't correct them
residtable<-cbind(cuttable[,c(1:5)],temp$resid,cuttable[,c(20:23,26:30)])
rownames(residtable)<-cuttable$TreeName
View(residtable)

#PICs are bivariate approaches to examine relationships between morphological variables:
FLLPIC <- pic(residtable$FLL, cuttree)
HLLPIC <- pic(residtable$HLL, cuttree)
limbsPIC <- lm(FLLPIC ~ HLLPIC, data = residtable)
summary(limbsPIC)
plot(FLLPIC ~ HLLPIC) #Observe that HLL positively explains FLL: they are correlated
abline(limbsPIC)
limbsPIC$residuals #Note that there is one less residue, this it was just to see the correlation

UFLPIC <- pic(residtable$UFL, cuttree)
UHLPIC <- pic(residtable$UHL, cuttree)
upperPIC <- lm(UFLPIC ~ UHLPIC, data = residtable)
summary(upperPIC)
plot(UFLPIC ~ UHLPIC)
abline(upperPIC)

LFLPIC <- pic(residtable$LFL, cuttree)
LHLPIC <- pic(residtable$LHL, cuttree)
lowerPIC <- lm(LFLPIC ~ LHLPIC, data = residtable)
summary(lowerPIC)
plot(LFLPIC ~ LHLPIC)
abline(lowerPIC)

SfingLPIC <- pic(residtable$SfingL, cuttree)
LfingLPIC <- pic(residtable$LfingL, cuttree)
fingPIC <- lm(SfingLPIC ~ LfingLPIC, data = residtable)
summary(fingPIC)
plot(SfingLPIC ~ LfingLPIC)
abline(fingPIC)

StoeLPIC <- pic(residtable$StoeL, cuttree)
LtoeLPIC <- pic(residtable$LtoeL, cuttree)
toePIC <- lm(StoeLPIC ~ LtoeLPIC, data = residtable)
summary(toePIC)
plot(StoeLPIC ~ LtoeLPIC)
abline(toePIC)

fingRPIC <- pic(residtable$fingR, cuttree)
toeRPIC <- pic(residtable$toeR, cuttree)
fingtoeRPIC <- lm(fingRPIC ~ toeRPIC, data = residtable)
summary(fingtoeRPIC)
plot(fingRPIC ~ toeRPIC)
abline(fingtoeRPIC)


#ANCESTRAL STATE RECONSTRUCTION:
#ASR for continuous variables:
FLLvsHLL_ASR <- cuttable$FLLvsHLL
FLLvsHLL_ASR
names(FLLvsHLL_ASR) <- row.names(cuttable)

FLLvsHLL_ASR_CI<-fastAnc(cuttree, FLLvsHLL_ASR, CI=TRUE)
FLLvsHLL_ASR_CI
as.numeric(FLLvsHLL_ASR)

plot(cuttree, label.offset=5, cex=0.45) #"label.offset" controls the distance between the label and the point
tiplabels(round(FLLvsHLL_ASR,4), adj=c(-0.2,0.5), cex=0.35, frame = "n",
          pch = NULL, thermo = NULL, pie = NULL, piecol = NULL,
          col = "mediumvioletred", horiz = FALSE,
          width = NULL, height = NULL)
nodelabels(round(FLLvsHLL_ASR[[1]],3), adj=c(-0.3,0.5), cex=0.35, frame="n", col = "slateblue")

contMap(cuttree,FLLvsHLL_ASR, fsize=0.45, lwd=1.5, legend=30,lims=c(min(0.4),max(1.0)), outline=FALSE)+
  tiplabels(round(FLLvsHLL_ASR,4), adj=c(1.5,0.5), cex=0.6, frame="n")+
  nodelabels(round(FLLvsHLL_ASR[[1]],3), adj=c(-0.5,0.5), cex=0.3, frame="n") #These 2 last lines give error, but I don't need them: ERROR: non-numeric argument to binary operator

FLL_ASR <- residtable$FLL
FLL_ASR
names(FLL_ASR) <- row.names(residtable)
plot(cuttree, label.offset=8, cex=0.35)
tiplabels(round(FLL_ASR,4), adj=c(-0.4,0.5), cex=0.3, frame = "n",
          pch = NULL, thermo = NULL, pie = NULL, piecol = NULL,
          col = "blue", horiz = FALSE,
          width = NULL, height = NULL)
nodelabels(round(FLL_ASR[[1]],3), adj=c(-0.3,0.5), cex=0.3, frame="n", col = "pink") #Again the ERROR for tiplabels and nodelabels, but don't need them
contMap(cuttree, FLL_ASR, fsize=0.45,lwd=1.5, legend=30,lims=c(min(-0.15),max(0.17)), outline=FALSE)

HLL_ASR <- residtable$HLL
HLL_ASR
names(HLL_ASR) <- row.names(residtable)
plot(cuttree, label.offset=8, cex=0.35)
tiplabels(round(HLL_ASR,4), adj=c(-0.4,0.5), cex=0.3, frame = "n",
          pch = NULL, thermo = NULL, pie = NULL, piecol = NULL,
          col = "blue", horiz = FALSE,
          width = NULL, height = NULL)
nodelabels(round(HLL_ASR[[1]],3), adj=c(-0.3,0.5), cex=0.3, frame="n", col = "pink")
obj3<-contMap(cuttree, HLL_ASR, fsize=0.45, lwd=1.5, legend=30, lims=c(min(-0.20),max(0.25)), outline=FALSE, direction="leftwards")
plot(obj3, fsize=0.45, lwd=1.5, legend=30, outline=FALSE, direction="leftwards")
pdf(file="plot-contMap.pdf",width=10,height=10)
plot(obj3,fsize=0.45, lwd=1.5, legend=30, outline=FALSE, direction="leftwards")
dev.off()

#phylogenetic TRAITGRAM, phenogram:
SVL_ASR <- residtable$SVL
SVL_ASR
names(SVL_ASR) <- row.names(residtable)
HL_ASR <- residtable$HL
HL_ASR
names(HL_ASR) <- row.names(residtable)
HW_ASR <- residtable$HW
HW_ASR
names(HW_ASR) <- row.names(residtable)
HH_ASR <- residtable$HH
HH_ASR
names(HH_ASR) <- row.names(residtable)
plotTree(cuttree,FLL_ASR,pts=F,node.numbers=T,fsize=0.36)
phenogram(cuttree,FLL_ASR,fsize=0.36,xlab="time",ylab="FLL")
phenogram(cuttree,HLL_ASR,fsize=0.36,xlab="time",ylab="HLL")
phenogram(cuttree,SVL_ASR,fsize=0.36,xlab="time",ylab="SVL")
phenogram(cuttree,FLLvsHLL_ASR,fsize=0.36,xlab="time",ylab="FLLvsHLL") 
phenogram(cuttree,HL_ASR,fsize=0.36,xlab="time",ylab="HL")
phenogram(cuttree,HW_ASR,fsize=0.36,xlab="time",ylab="HW")
phenogram(cuttree,HH_ASR,fsize=0.36,xlab="time",ylab="HH")


#ASR for structural niches, discrete variable:
library(diversitree)
library(geiger)

niches<-setNames(cuttable$Structural_niche,cuttable$TreeName)
niches 
table(niches) 
niches[cuttree$tip.label] #Not necessari

plotTree(cuttree,type="fan",fsize=0.4,ftype="i")
#I don't know why, but RUN TWICE the following line:
cols<-setNames(palette(c("green4","yellowgreen","mediumpurple3","lightsalmon4","orange1","slategray3","mediumturquoise","hotpink2"))[1:length(unique(niches))],sort(unique(niches))) 
tiplabels(pie=to.matrix(niches,sort(unique(niches))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=0.9*par()$usr[1],
                  y=-max(nodeHeights(tree)),fsize=0.8)

#Now estimate ancestral states under a ER model:
nichesASR.ER<-ace(niches,cuttree,model="ER",type="discrete")
round(nichesASR.ER$lik.anc,3) #The element lik.anc gives us the marginal ancestral states, also known as the 'empirical Bayesian posterior probabilities.'
plotTree(cuttree,fsize=0.4,ftype="i")+
  nodelabels(node=1:cuttree$Nnode+Ntip(cuttree),
           pie=nichesASR.ER$lik.anc,piecol=cols,cex=0.5)+
  tiplabels(pie=to.matrix(niches,sort(unique(niches))),piecol=cols,cex=0.3)+
  add.simmap.legend(colors=cols,prompt=FALSE,x=0.9*par()$usr[1],
                  y=-max(nodeHeights(cuttree)),fsize=0.8)

round(nichesASR.ER$lik.anc,4) #To quantify the probabilities in the circle graphs

#Simulate single stochastic character map using empirical Bayes method:
discrtree<-make.simmap(cuttree,niches,model="ER")

plot(discrtree,cols,type="fan",fsize=0.4,ftype="i")+
  add.simmap.legend(colors=cols,prompt=FALSE,x=0.9*par()$usr[1],
                  y=-max(nodeHeights(tree)),fsize=0.8)

#At the end I didn't use the following:
#A single stochastic character map does not mean a whole lot in isolation, creating a sample of 600 stochastic character maps from our dataset:
mdiscrtrees<-make.simmap(cuttree,niches,model="ER",nsim=600)
null<-sapply(mdiscrtrees,plot,colors=cols,lwd=1,ftype="off")

pd<-summary(mdiscrtrees,plot=TRUE) #Summary of the set of stochastic maps:

plot(pd,fsize=0.2,ftype="i")

plot(mdiscrtrees[[1]],cols,type="fan",fsize=0.3,ftype="i")+
  nodelabels(pie=pd$ace,piecol=cols,cex=0.3)+
  add.simmap.legend(colors=cols,prompt=FALSE,x=0.9*par()$usr[1],
                  y=-max(nodeHeights(tree)),fsize=0.8)
plot(mdiscrtrees[[1]],cols,fsize=0.45,ftype="i")+ #I lo mateix sense tipus abanico
  nodelabels(pie=pd$ace,piecol=cols,cex=0.3)+
  add.simmap.legend(colors=cols,prompt=FALSE,x=0.12*par()$usr[1],y=79,fsize=0.6) #He d'arreglar la posició de la llegenda


plotTree(ladderize(cuttree),type="arc",fsize=0.3,lwd=1,ftype="i",arc_height=1.5) #"Ladderize" inverts the sides

plotTree(cuttree, offset=3, type="arc",fsize=0.5,lwd=1,ftype="i",arc_height=1.5)+
  nodelabels(node=1:cuttree$Nnode+Ntip(cuttree),
             pie=nichesASR.ER$lik.anc,piecol=cols,cex=0.3)+
  tiplabels(pie=to.matrix(niches,sort(unique(niches))),piecol=cols,cex=0.17)+
  add.simmap.legend(colors=cols,prompt=FALSE,x=0.8*par()$usr[1],
                    y=409,fsize=0.6)

plotTree(cuttree,type="arc",fsize=0.3,lwd=1,ftype="i",arc_height=1.5) 
legend(x=0,y=1,names(cols),col=cols,pch=15,xjust=0.5,horiz=TRUE,cex=0.8,bty="n",pt.cex=1.5)

#For NICHE SHIFTS:
remotes::install_github("liamrevell/phytools") #Download the newest version of phytools
library(phytools)
packageVersion("phytools")

plotTree(as.phylo(discrtree),offset=2,type="arc",lwd=1,fsize=0.3,ftype="i",arc_height=2.2)
tiplabels(pie=to.matrix(niches,sort(unique(niches))),piecol=cols,cex=0.3)
text(0,0.1,"Niche shifts",adj=0.5)
nicheshifts<-fitMk(cuttree,niches,model="ARD",
                    opt.method="optimParallel",lik.func="pruning",
                    pi="fitzjohn",cex=0.01)
xy<-plot(as.Qmatrix(nicheshifts),show.zeros=FALSE,
         add=TRUE,ylim=c(-1.2,3),xlim=c(-2,2),spacer=0.2,
         text=FALSE,width=TRUE,max.lwd=5,cex=0.01)
library(plotrix)
nulo<-mapply(draw.circle,xy$x,xy$y,col=cols,
             MoreArgs=list(radius=0.2))
text(xy$x,xy$y,xy$states,cex=0.01,col="black")

plot(as.Qmatrix(nicheshifts),show.zeros=TRUE, width=TRUE, tol = 1e-15, color=TRUE, offset=0.025)
as.Qmatrix(nicheshifts) #To quantify transition rates


#PRINCIPAL COMPONENT ANALYSES, phylPCA:
residtablePCA <- residtable[c(1:79),c(1:10,13:14,17:21)]
phylPCA <- phyl.pca(cuttree, Y=residtablePCA[c(1:79),(5:17)], method="BM", mode="corr")
phylPCA
View(residtablePCA)
summary(phylPCA)
phylPCA$L
barplot(phylPCA$Evec[,1],cex.names=0.40)
barplot(phylPCA$Evec[,2],cex.names=0.40)
phylPCA$Eval #Eigenvalues represent the scaling factor by which a vector is transformed when a linear transformation is applied, while eigenvectors represent the directions in which the transformation occurs.
phylPCA$L
plot(phylPCA$S) #To get the coordenates for the individuals
biplot(phylPCA$S,phylPCA$L,cex=c(0.4,0.5)) 
indphylPCA <-scores(phylPCA) #Another way to get coordenates
indphylPCA
pcatable1 <- as.data.frame(cbind(phylPCA$S[,c(1,2)],residtable$TreeName,residtable$Structural_niche))
colnames(pcatable1) <- c("PC1","PC2","taxon","group")
pcatable2 <- as.data.frame(cbind(phylPCA$S[,c(1,2)],residtable$TreeName,residtable$Genus))
colnames(pcatable2) <- c("PC1","PC2","taxon","group")
head(pcatable1)
head(pcatable2)


phy <- cuttree
dat <- residtable[c(1:79),c(1:10,13:14,17:21,27:30)]
ppca <- phyl.pca(phy,dat[,c(5:17)],mode = "cor") #this will give you a warning about lacking row names. You can ignore it because sort=TRUE in treedata()
loadings <- as.data.frame(ppca$S) #gives you the loadings so you can plot them
loadings$group <- dat$Structural_niche #this gives identity back to groups to that you can work with them. Has to use dat and not data because of the sort=TRUE in treedata()

phy <- cuttree
dat <- residtable[c(1:79),c(1:10,13:14,17:21,27:30)]
ppca <- phyl.pca(phy,dat[,c(5:17)],mode = "cor") 
loadings <- as.data.frame(ppca$S)
loadings$group <- dat$Genus #Run this batch if I want to differentiate for the genera


library(plyr)
chulls <- ddply(loadings, .(group), function(df) df[chull(df$PC1, df$PC2),]) #This makes hulls to plot around your data groups, which are much more honest than elipses.

library(ggrepel) 
ggplot(data=loadings,aes(x=PC1,y=PC2, colour = group, frame = TRUE)) + 
  geom_point() + 
  geom_text_repel(label=row.names(residtable), size=2) +
  theme_classic() +
  geom_polygon(data=chulls, aes(x=PC1, y=PC2, fill=group), alpha=0.2)


#To PLOT THE phylPCA: PHYLOMORPHOSPACE
pcatable1[,c(1:2)]<-sapply(pcatable1[,c(1:2)],FUN = as.numeric)
pcatable2[,c(1:2)]<-sapply(pcatable2[,c(1:2)],FUN = as.numeric)

library(grid)
library(ggplot2)
library(ggphylomorpho)

library(tidyverse)
ggplot(data=pcatable1,aes(x=PC1, y=PC2, color=group, label=rownames(pcatable2)))+
  geom_point()+
  theme_bw()

ggphylomorpho(tree=cuttree, tipinfo=pcatable2, xvar = PC1, yvar = PC2, factorvar = group,
              labelvar = taxon, title = "Phylomorphospace", xlab = "PC1",
              ylab = "PC2", repel = FALSE, edge.width = 0.4, fontface = "italic",
              tree.alpha = 0.7)
indphylPCA[,c(1:3)] #I crec que això tmb ho puc borrar
ggphylomorpho(tree=cuttree, tipinfo=indphylPCA[,c(1:3)]) #I això també

library(phytools)
phylomorphospace(tree, X, A=NULL, label=c("radial","horizontal","off"),
                 control=list(), ...)
project.phylomorphospace(tree, X, nsteps=200, sleep=0,
                         direction=c("to","from","both"), ...)



#MANOVA i phylMANOVA:
manova1 <- manova(cbind(FLLvsHLL,FLtoes,HLtoes,LtoevsFLL,LtoevsHLL,FLL,HLL,SFtoeL,LFtoeL,SHtoeL,LHtoeL,HL,HW,HH,AHL,SVL) ~ Complexity, data = cuttable)
manova1
summary(manova1)

library(geiger)

library(data.table) #In here I converted the residtable into a data table but it must not be the same thing as a data set.
data_frame <- residtable
print ("Original DataFrame")
print (data_frame)
setDT(data_frame) #Converting into data.table
print ("Resultant DataFrame")
print (data_frame)

mar=get(data(geospiza))
geo=get(data(geospiza))
dat=geo$dat
grp<-as.factor(c(rep(0, 7), rep(1, 6)))
names(grp)=rownames(dat)
aov.phylo(dat~grp, geo$phy, nsim = 50, test = c("Wilks"))
x=aov.phylo(residtable$Structural_niche, cuttree, nsim=50, test="Wilks")

sniche <- as.factor(residtable$Structural_niche)
class(sniche)
names(sniche) <- row.names(residtable)
deps<-as.matrix(residtable[,c(5:25)])

#Here:
pmanovaresults<-aov.phylo(deps ~ sniche, cuttree, nsim=1000, test="Wilks")
summary(pmanovaresults)
pmanovaresults$effects

ratematrix(cuttree,deps)

residues=get(data(cuttable))
dat=residues$dat
d1=dat[,1]
grp<-as.factor(c(rep(0, 7), rep(1, 6)))
grp
names(Structural_niche)=rownames(dat)
x=aov.phylo(dat~grp, geo$phy, nsim=50, test="Wilks")
print(attributes(x)$summary) # summary table
aov.phylo(formula, phy, nsim = 1000, 
          test = c("Wilks"


#phylANOVAs:
phylANOVA(cuttree, residtable$Structural_niche, residtable$HL, nsim=1000, posthoc=TRUE, p.adj="bonferroni")
phylANOVA(cuttree, residtable$Structural_niche, residtable$HW, nsim=1000, posthoc=TRUE, p.adj="bonferroni")
phylANOVA(cuttree, residtable$Structural_niche, residtable$HH, nsim=1000, posthoc=TRUE, p.adj="bonferroni") 
phylANOVA(cuttree, residtable$Structural_niche, residtable$ED, nsim=1000, posthoc=TRUE, p.adj="bonferroni") 
phylANOVA(cuttree, residtable$Structural_niche, residtable$AHL, nsim=1000, posthoc=TRUE, p.adj="bonferroni") 
phylANOVA(cuttree, residtable$Structural_niche, residtable$FLL, nsim=1000, posthoc=TRUE, p.adj="bonferroni")
phylANOVA(cuttree, residtable$Structural_niche, residtable$HLL, nsim=1000, posthoc=TRUE, p.adj="bonferroni")
phylANOVA(cuttree, residtable$Structural_niche, residtable$UFL, nsim=1000, posthoc=TRUE, p.adj="bonferroni")
phylANOVA(cuttree, residtable$Structural_niche, residtable$LFL, nsim=1000, posthoc=TRUE, p.adj="bonferroni")
phylANOVA(cuttree, residtable$Structural_niche, residtable$UHL, nsim=1000, posthoc=TRUE, p.adj="bonferroni")
phylANOVA(cuttree, residtable$Structural_niche, residtable$LHL, nsim=1000, posthoc=TRUE, p.adj="bonferroni")
phylANOVA(cuttree, residtable$Structural_niche, residtable$SfingL, nsim=1000, posthoc=TRUE, p.adj="bonferroni")
phylANOVA(cuttree, residtable$Structural_niche, residtable$LFtoeL, nsim=1000, posthoc=TRUE, p.adj="bonferroni")
phylANOVA(cuttree, residtable$Structural_niche, residtable$StoeL, nsim=1000, posthoc=TRUE, p.adj="bonferroni")
phylANOVA(cuttree, residtable$Structural_niche, residtable$LtoeL, nsim=1000, posthoc=TRUE, p.adj="bonferroni")
phylANOVA(cuttree, residtable$Structural_niche, residtable$UHW, nsim=1000, posthoc=TRUE, p.adj="bonferroni")
phylANOVA(cuttree, residtable$Structural_niche, residtable$FLLvsHLL, nsim=1000, posthoc=TRUE, p.adj="bonferroni")
phylANOVA(cuttree, residtable$Structural_niche, residtable$fingR, nsim=1000, posthoc=TRUE, p.adj="bonferroni")
phylANOVA(cuttree, residtable$Structural_niche, residtable$toeR, nsim=1000, posthoc=TRUE, p.adj="bonferroni")
phylANOVA(cuttree, residtable$Complexity, residtable$LfingvsFLL, nsim=1000, posthoc=TRUE, p.adj="bonferroni")
phylANOVA(cuttree, residtable$Complexity, residtable$LtoevsHLL, nsim=1000, posthoc=TRUE, p.adj="bonferroni")


class(residtable$FLL)
shapiro.test(residtable$SVL)
shapiro.test(residtable$AHL) 
shapiro.test(residtable$HL) #Not normal
shapiro.test(residtable$HW) #Not normal
shapiro.test(residtable$HH) #Not normal
shapiro.test(residtable$ED) 
shapiro.test(residtable$UFL) #Not normal
shapiro.test(residtable$LFL)
shapiro.test(residtable$SFtoeL) #Not normal
shapiro.test(residtable$LFtoeL) #Not normal
shapiro.test(residtable$UHL) 
shapiro.test(residtable$LHL) 
shapiro.test(residtable$UHW) #Not normal
shapiro.test(residtable$SFtoeL) #Not normal
shapiro.test(residtable$LFtoeL) #Not normal
shapiro.test(residtable$HLL) 
shapiro.test(residtable$FLL) #Not normal
shapiro.test(residtable$FLLvsHLL) 
#Only half of them have a normal distribution

#To see if there are outliers in the data which make it not normal:
summary(residtable$HL) #max or min 
hist(residtable$HL)
hist(residtable$HW)
hist(residtable$HH)
hist(residtable$UFL)
hist(residtable$HLL)
hist(residtable$SFtoeL)
hist(residtable$LFtoeL)
hist(residtable$UHW)
hist(residtable$SHtoeL)
hist(residtable$LHtoeL)
hist(residtable$FLL)
boxplot(cuttable$FLL)

#To also do TL (didn't do it before as it would have dropped all species that didn't have measurements for TL)
tableTL<-table[c(1:87),c(1:4,8:23,25:35)] #Run without the TL to avoid NAs.
tableTL<-na.omit(tableTL)
View(tableTL)
row.names(tableTL)<-tableTL$TreeName
cut<-treedata(tree,tableTL,sort=TRUE)
str(cut)
cuttreeTL<-cut$phy
is.binary(cuttreeTL)
is.ultrametric(cuttreeTL)
cuttableTL<-cut$data
plot(cuttreeTL,cex=0.45)
View(cuttableTL)
cuttableTL<-as.data.frame(cuttableTL)
row.names(cuttableTL)<-cuttableTL$TreeName
cuttableTL[,c(5:26)] <- sapply(cuttableTL[,c(5:26)],function(x) as.numeric(x))
temp<-phyl.resid(cuttree,cuttableTL$SVL,cuttableTL[,c(6:19,24:25)],method = "BM")
residtableTL<-cbind(cuttableTL[,c(1:5)],temp$resid,cuttableTL[,c(20:23,26:30)])
rownames(residtableTL)<-cuttableTL$TreeName
View(residtableTL)


##(1) linear regression to visualitze the ratio:

residtable<-as.data.frame(residtable)
residtableD <- filter(residtable, Genus == "Draco") #I filter for Draco to see Draco's tendline
row.names(residtableD)<-residtableD$TreeName
cutD<-treedata(tree,residtableD,sort=TRUE)
str(cutD)
cuttree<-cutD$phy
is.binary(cuttree)
is.ultrametric(cuttree)
residtableD<-cutD$data
plot(cuttree)
View(residtableD)
residtableD<-as.data.frame(residtableD)
row.names(residtableD)<-residtableD$TreeName
residtableD[,c(5:22)] <- sapply(residtableD[,c(5:22)],function(x) as.numeric(x))
temp<-phyl.resid(cuttree,residtableD$SVL,residtableD[,c(6:22)],method = "BM")
residtableD<-cbind(residtableD[,c(1:5)],temp$resid,residtableD[,c(23:26)])
rownames(residtableD)<-residtableD$TreeName

residtable<-as.data.frame(residtable)
residtableOS <- filter(residtable, Genus == "Otocryptis"| Genus == "Sitana") #Also filter for these ones
row.names(residtableOS)<-residtableOS$TreeName
cutOS<-treedata(tree,residtableOS,sort=TRUE)
str(cutOS)
cuttree<-cutOS$phy
is.binary(cuttree)
is.ultrametric(cuttree)
residtableOS<-cutOS$data
plot(cuttree)
View(residtableOS)
residtableOS<-as.data.frame(residtableOS)
row.names(residtableOS)<-residtableOS$TreeName
residtableOS[,c(5:22)] <- sapply(residtableOS[,c(5:22)],function(x) as.numeric(x))
temp<-phyl.resid(cuttree,residtableOS$SVL,residtableOS[,c(6:22)],method = "BM")
residtableOS<-cbind(residtableOS[,c(1:5)],temp$resid,residtableOS[,c(23:26)])
rownames(residtableOS)<-residtableOS$TreeName

limblengths <- lm(cuttable$FLL ~ cuttable$HLL)
summary(limblengths)
row.names(cuttable) <- cuttable$TreeName
library(ggrepel)
ggplot(data=cuttable,aes(x=HLL,y=FLL,colour=Genus))+
  geom_text_repel(label=row.names(cuttable), size=1)+ 
  geom_point()+
  theme_classic()