#Loading a tree and performing and ancestral character state reconstruction
library(diversitree)
library(phytools)
library(paleotree)
library(ggtree)
library(geiger)

# Read in Dataset 
#Note: you need to make sure your ID's in your dataset match the tips of your tree.

RaboskyFamTree<- read.tree("Rabosky_FamNames_TimeCal.tre")
plot.phylo(RaboskyFamTree, type = "fan", cex = .47, lwd= 0.4, rotate.tree = 180, use.edge.length = F, label.offset = 0.04, node.width = 0.2, node.depth = 2)
tiplabels(cex = 0.2)

RaboskyFamTree <- ladderize(RaboskyFamTree) 

RaboskyData <- read.csv("Rabosky_BiofluorData.csv")
summary(RaboskyData)


#----Biofluor Absense/Pres ASR-----
Rabosky_Biofluor <- RaboskyData[,5]
Rabosky_Biofluor

Rabosky_Biofluor <- setNames(Rabosky_Biofluor, RaboskyData$Family) 
str(Rabosky_Biofluor)


Rabosky_Biofluor_Tree_ER <- ace(Rabosky_Biofluor, RaboskyFamTree, model = "ER", type = "discrete", CI = TRUE)
Rabosky_Biofluor_Tree_ER
round(Rabosky_Biofluor_Tree_ER$lik.anc,3)

#Plot
quartz(width=15, height = 15)

plot.phylo(RaboskyFamTree, type = "fan", cex = .3, lwd= 0.4, rotate.tree = 180, 
           use.edge.length = F, label.offset = 0.04, node.width = 0.2, node.depth = 2)
nodelabels(cex = 0.3)

nodelabels(pie = Rabosky_Biofluor_Tree_ER$lik.anc, cex = 0.12, piecol = c('white', 'skyblue'))

