##set up analysis

setwd("Desktop/R analyses/diel mapping")

library(ape)
library(phytools)


##read in tree

##ericson tree

erictree<-read.nexus("Ericson1914.tre")
tree<-erictree


##hackett tree

hacktree<-read.nexus("Hackett1914.tre")
tree<-hacktree




##read in data


## diurnal data

state<-read.csv("diu1914.csv",row.names=1)
state<-as.matrix(state)[,1]
x<-state



## nocturnal data

state<-read.csv("noc1914.csv",row.names=1)
state<-as.matrix(state)[,1]
x<-state



## four state data
state<-read.csv("fourstate1914.csv",row.names=1)
state<-as.matrix(state)[,1]
x<-state





## perform stochastic character mapping with a likelihood model with all-rates-different between states (ARD), root state estimated based on stationary frequencies, and 25 replicates

mtrees<-make.simmap(tree,x,model="ARD", pi="estimated",nsim=25)


## perform stochastic character mapping with a likelihood model with equal-rates- between states (ER), root state estimated based on stationary frequencies, and 25 replicates

mtrees<-make.simmap(tree,x,model="ER", pi="estimated",nsim=25)


## plot trees from simmap

par(mfrow=c(5,5))

cols<-setNames(palette()[1:length(unique(x))],sort(unique(x)))

null<-sapply(mtrees,plot,lwd=0.1,ftype="off")




##perform ancestral state reconstruction using ARD model

fitARD<-ace(x,tree,model="ARD",type="discrete")

round(fitARD$lik.anc,3)


##plot tree with node labels

plotTree(tree,fsize=0.025,lwd=0.2)

cols<-setNames(palette()[1:length(unique(x))],sort(unique(x)))

nodelabels(node=1:tree$Nnode+Ntip(tree), pie=fitARD$lik.anc,piecol=cols,cex=0.2 ,lwd=0.05)