######## required packages library(ape) library(nlme) library(phytools) tr <- read.tree(text = "(((Hypotrigona:61, (Meliponula:43,(Apotrigona:33,Axestotrigona:33):10):18):12,Tetragonula:73):8,(Leurotrigona:71,(Melipona:57,((((Nannotrigona:33,(Lestrimelitta:24,Plebeia:24):9):2,(Tetragonisca:26,Frieseomelitta:26):9):3,(Scaptotrigona:36,(Geotrigona:32,(Tetragona:30,(Trigona:27,Cephalotrigona:27):3):2):4):2):5,Partamona:43):14):14):10);") plot(tr) #################### ancestral state reconstruction for discrete variables, see http://www.phytools.org/eqg2015/asr.html data <- import Data table for Fig. 10.11 attach(data) states<-as.matrix(data)[,1] states fitER<-ace(states,tr,model="ER",type="discrete") fitER plotTree(tr, fsize=0.8,ftype="i") nodelabels(node=1:tr$Nnode+Ntip(tr), pie=fitER$lik.anc,piecol= c("green","blue","red","yellow"),cex=0.5) tiplabels(pie=to.matrix(states,sort(unique(states))),piecol=c("green","blue","red","yellow"),cex=0.3) add.simmap.legend(colors=c("green","blue","red","yellow"),prompt=FALSE,states=0.9*par()$usr[1], y=-max(nodeHeights(tr)),fsize=0.8) ############### 2nd approach uses MCMC mtrees<-make.simmap(tr,states,model="ER",nsim=1000) pd<-summary(mtrees,plot=FALSE) pd plot(pd,fsize=0.6,ftype="i") plot(mtrees[[1]], fsize=0.8,ftype="i") nodelabels(pie=pd$ace,piecol= c("black","red"),cex=0.5) # obj<-contMap(tr,states) obj<-densityMap(mtrees,lwd=4,outline=TRUE) ## change colour scheme n<-length(obj$cols) obj$cols[1:n]<-colorRampPalette(c("skyblue4","yellow"), space="Lab")(n) plot(obj, lwd = 4, outline = TRUE) nodelabels(node=1:tr$Nnode+Ntip(tr), pie=fitER$lik.anc,piecol= c("skyblue4","yellow"),cex=0.6)