library(ape) library(geiger) library(phytools) #load tree and character data, need to do this for each timetree tree <- ladderize(read.nexus("bayes_node.tre")) data <- read.table("charmap.txt",header=TRUE,row.names=1) #parse table columns data <- data[tree$tip.label, ] eco <- as.matrix(data)[,1] ana <- as.matrix(data)[,2] #this character codes segment addition, not analyzed in the paper #resolve polytomies and zero length branches tree <- multi2di(tree) tree$edge.length[tree$edge.length==0]<-max(nodeHeights(tree))*1e-6 #test best fitting model of discrete character evolution #don't need SYM because it is identical to ER for binary characters ER_eco <- ace(eco, tree, type="discrete", model="ER") ARD_eco <- ace(eco, tree, type="discrete", model="ARD") ER_eco$loglik ARD_eco$loglik 1-pchisq(2*abs(ER_eco$loglik - ARD_eco$loglik), 1) #stochastic character mapping under best fitting model map_eco <- make.simmap(tree,eco,nsim=500,model="ER") #plot pdf("bayes_node_eco.pdf") densityMap(map_eco,fsize=0.5) dev.off()