library(ape) 
library(hisse) 
ericson <- read.nexus("Ericson1914.tree") 
hackett <- read.nexus("Hackett1914.tree") 
labels <- read.csv("appendix1914.csv")

species <- labels$Species
Nstates <- labels$ActNOC
Dstates <- labels$ActDIU 

## Nocturnal and diurnal data frames
dfN <- data.frame(species, Nstates)
dfD <- data.frame(species, Dstates)

## Accounting for incomplete sampling
Nsampling<-c(957/16615, 957/16068)
Dsampling<-c(1103/10187, 811/13520)

## Making HiSSE transition matrices
TransTrue<-TransMatMaker(hidden.states=TRUE)
TransFalse<-TransMatMaker(hidden.states=FALSE)

## Running four full hisse models with default root setting - madfitz
## ericson nocturnal
hisseEN<-hisse(ericson, dfN, f=Nsampling, hidden.states=TRUE, turnover.anc=c(1,2,3,4), eps.anc=c(1,2,3,4), trans.rate=TransTrue, output.type="raw")
## ericson diurnal
hisseED<-hisse(ericson, dfD, f=Dsampling, hidden.states=TRUE, turnover.anc=c(1,2,3,4), eps.anc=c(1,2,3,4), trans.rate=TransTrue, output.type="raw")
## hackett nocturnal 
hisseHN<-hisse(hackett, dfN, f=Nsampling, hidden.states=TRUE, turnover.anc=c(1,2,3,4), eps.anc=c(1,2,3,4), trans.rate=TransTrue, output.type="raw")
## hackett diurnal
hisseHD<-hisse(hackett, dfD, f=Dsampling, hidden.states=TRUE, turnover.anc=c(1,2,3,4), eps.anc=c(1,2,3,4), trans.rate=TransTrue, output.type="raw")

## Running the HiSSE two-rate models equivalent to a restricted BiSSE model
## ericson nocturnal 
hisseEN.two<-hisse(ericson, dfN, f=Nsampling, hidden.states=FALSE, turnover.anc=c(1,1,0,0), eps.anc=c(1,1,0,0), trans.rate=TransFalse, output.type="raw")
## ericson diurnal 
hisseED.two<-hisse(ericson, dfD, f=Dsampling, hidden.states=FALSE, turnover.anc=c(1,1,0,0), eps.anc=c(1,1,0,0), trans.rate=TransFalse, output.type="raw")
## hackett nocturnal 
hisseHN.two<-hisse(hackett, dfN, f=Nsampling, hidden.states=FALSE, turnover.anc=c(1,1,0,0), eps.anc=c(1,1,0,0), trans.rate=TransFalse, output.type="raw")
## hackett diurnal
hisseHD.two<-hisse(hackett, dfD, f=Dsampling, hidden.states=FALSE, turnover.anc=c(1,1,0,0), eps.anc=c(1,1,0,0), trans.rate=TransFalse, output.type="raw")

## Running null4 models 
## ericson nocturnal 
EN.null4<-hisse.null4(ericson, dfN, f=Nsampling, turnover.anc=rep(c(1,2,3,4),2), eps.anc=rep(c(1,2,3,4),2), trans.type="equal", output.type)
## ericson diurnal 
ED.null4<-hisse.null4(ericson, dfD, f=Dsampling, turnover.anc=rep(c(1,2,3,4),2), eps.anc=rep(c(1,2,3,4),2), trans.type="equal")
## hackett nocturnal 
HN.null4<-hisse.null4(hackett, dfN, f=Nsampling, turnover.anc=rep(c(1,2,3,4),2), eps.anc=rep(c(1,2,3,4),2), trans.type="equal")
## hackett diurnal  
HD.null4<-hisse.null4(hackett, dfD, f=Dsampling, turnover.anc=rep(c(1,2,3,4),2), eps.anc=rep(c(1,2,3,4),2), trans.type="equal") 

## Running null2 models 
TransNull2<-ParDrop(TransTrue, c(3,5,8,10))
TransNull2[!is.na(TransNull2) & !TransNull2 == 0] = 1
## ericson nocturnal
EN.null2<-hisse(ericson, dfN, hidden.states=TRUE, turnover.anc=c(1,1,2,2), eps.anc=c(1,1,2,2), trans.rate=TransNull2, output.type="raw") 
## ericson diurnal 
ED.null2<-hisse(ericson, dfD, hidden.states=TRUE, turnover.anc=c(1,1,2,2), eps.anc=c(1,1,2,2), trans.rate=TransNull2, output.type="raw")
## hackett nocturnal 
HN.null2<-hisse(hackett, dfN, hidden.states=TRUE, turnover.anc=c(1,1,2,2), eps.anc=c(1,1,2,2), trans.rate=TransNull2, output.type="raw")
## hackett diurnal 
HD.null2<-hisse(hackett, dfD, hidden.states=TRUE, turnover.anc=c(1,1,2,2), eps.anc=c(1,1,2,2), trans.rate=TransNull2, output.type="raw")

## Testing alternative root state (equal) in full hisse model 
## ericson nocturnal
hisseEN.equal<-hisse(ericson, dfN, f=Nsampling, hidden.states=TRUE, turnover.anc=c(1,2,3,4), eps.anc=c(1,2,3,4), trans.rate=TransTrue, root.type="equal")
## ericson diurnal
hisseED.equal<-hisse(ericson, dfD, f=Dsampling, hidden.states=TRUE, turnover.anc=c(1,2,3,4), eps.anc=c(1,2,3,4), trans.rate=TransTrue, root.type="equal")
## hackett nocturnal 
hisseHN.equal<-hisse(hackett, dfN, f=Nsampling, hidden.states=TRUE, turnover.anc=c(1,2,3,4), eps.anc=c(1,2,3,4), trans.rate=TransTrue, root.type="equal")
## hackett diurnal
hisseHD.equal<-hisse(hackett, dfD, f=Dsampling, hidden.states=TRUE, turnover.anc=c(1,2,3,4), eps.anc=c(1,2,3,4), trans.rate=TransTrue, root.type="equal")

## load diversitree to run full BiSSE model comparison using make.bisse
library(diversitree)
## naming vector 
names(Nstates)<-labels$Species
names(Dstates)<-labels$Species
## starting point functions for each tree
ESP<- starting.point.bisse(ericson) 
HSP<- starting.point.bisse(hackett)
## ericson nocturnal 
ENlik<-make.bisse(ericson, Nstates, sampling.f=Nsampling) 
bisseEN<-find.mle(ENlik, ESP)
## ericson diurnal 
EDlik<-make.bisse(ericson, Dstates, sampling.f=Dsampling) 
bisseED<-find.mle(EDlik, ESP)
## hackett nocturnal 
HNlik<-make.bisse(hackett, Nstates, sampling.f=Nsampling) 
bisseHN<-find.mle(HNlik, HSP)
## hackett diurnal 
HDlik<-make.bisse(hackett, Dstates, sampling.f=Dsampling) 
bisseHD<-find.mle(HDlik, HSP)

## retrieve BiSSE AIC values 
AIC(bisseEN)
AIC(bisseED) 
AIC(bisseHN)
AIC(bisseHD) 

## Ancestral state reconstructions of full hisse models with madfitz root state
## ericson nocturnal
EN.asr<-MarginRecon(ericson, dfN, f=Nsampling, pars=hisseEN$solution, hidden.states=TRUE)
## ericson diurnal
ED.asr<-MarginRecon(ericson, dfD, f=Dsampling, pars=hisseED$solution, hidden.states=TRUE)
## hackett nocturnal
HN.asr<-MarginRecon(hackett, dfN, f=Nsampling, pars=hisseHN$solution, hidden.states=TRUE)
## hackett diurnal 
HD.asr<-MarginRecon(hackett, dfD, f=Dsampling, pars=hisseHD$solution, hidden.states=TRUE)

## plotting HiSSE model results
## pie chart colors
col<- setNames(c("yellow", "blue", "gold", "darkblue"), c("0A", "1A", "0B", "1B"))
## node matrice results
## ericson nocturnal 
ENpies <- EN.asr$node.mat[,-1]
ENtips <- EN.asr$tip.mat[,-1] 
## ericson diurnal 
EDpies <- ED.asr$node.mat[,-1]
EDtips <- ED.asr$tip.mat[,-1] 
## hackett nocturnal 
HNpies <- HN.asr$node.mat[,-1]
HNtips <- HN.asr$tip.mat[,-1] 
## hackett diurnal
HDpies <- HD.asr$node.mat[,-1]
HDtips <- HD.asr$tip.mat[,-1] 

## pdf of phylogeny with pie charts
## ericson nocturnal 
pdf("FigureA3_EN.pdf", height = 200, width = 30) 
plot.phylo(ericson, type = "phylogram", cex = 0.25, label.offset = 0.75, show.tip.label = TRUE) 
nodelabels(pie = ENpies, cex = 0.2, piecol = col)
tiplabels(pie = ENtips, cex = 0.065, piecol = col)
legend("bottomleft", legend = c("Diurnal A", "Diurnal B", "Nocturnal A", "Nocturnal B"), cex = 5, pch = 21, pt.bg = c("yellow", "gold", "blue", "darkblue"), pt.cex = 9) 
dev.off() 
## ericson diurnal
pdf("FigureA1_ED.pdf", height = 200, width = 30) 
plot.phylo(ericson, type = "phylogram", cex = 0.25, label.offset = 0.75, show.tip.label = TRUE) 
nodelabels(pie = EDpies, cex = 0.2, piecol = col)
tiplabels(pie = EDtips, cex = 0.065, piecol = col)
legend("bottomleft", legend = c("Diurnal A", "Diurnal B", "Nocturnal A", "Nocturnal B"), cex = 5, pch = 21, pt.bg = c("yellow", "gold", "blue", "darkblue"), pt.cex = 9) 
dev.off() 
## hackett nocturnal 
pdf("FigureA4_HN.pdf", height = 200, width = 30) 
plot.phylo(hackett, type = "phylogram", cex = 0.25, label.offset = 0.75, show.tip.label = TRUE) 
nodelabels(pie = HNpies, cex = 0.2, piecol = col)
tiplabels(pie = HNtips, cex = 0.065, piecol = col)
legend("bottomleft", legend = c("Diurnal A", "Diurnal B", "Nocturnal A", "Nocturnal B"), cex = 5, pch = 21, pt.bg = c("yellow", "gold", "blue", "darkblue"), pt.cex = 9) 
dev.off() 
## hackett diurnal 
pdf("FigureA2_HD.pdf", height = 200, width = 30) 
plot.phylo(hackett, type = "phylogram", cex = 0.25, label.offset = 0.75, show.tip.label = TRUE) 
nodelabels(pie = HDpies, cex = 0.2, piecol = col)
tiplabels(pie = HDtips, cex = 0.065, piecol = col)
legend("bottomleft", legend = c("Diurnal A", "Diurnal B", "Nocturnal A", "Nocturnal B"), cex = 5, pch = 21, pt.bg = c("yellow", "gold", "blue", "darkblue"), pt.cex = 9) 
dev.off() 

## end