# function plots posterior density of mapped states from stochastic mapping # written by Liam J. Revell 2012 densityMap<-function(trees,res=100,fsize=NULL,check=FALSE,legend=NULL){ tol<-1e-10 if(class(trees)!="multiPhylo") stop("trees not 'multiPhylo' object; just use plotSimmap") h<-sapply(trees,function(x) max(nodeHeights(x))) steps<-0:res/res*max(h) trees<-rescaleSimmap(trees,totalDepth=max(h)) if(check){ X<-matrix(FALSE,length(trees),length(trees)) for(i in 1:length(trees)) X[i,]<-sapply(trees,all.equal.phylo,current=trees[[i]]) if(!all(X)) stop("some of the trees don't match in topology or relative branch lengths") } tree<-trees[[1]] H<-nodeHeights(tree) message("sorry - this might take a while; please be patient") for(i in 1:nrow(tree$edge)){ YY<-cbind(c(H[i,1],steps[intersect(which(steps>H[i,1]),which(stepsH[i,1]),which(steps1){ for(k in 2:length(trees[[j]]$maps[[i]])){ XX[k,1]<-XX[k-1,2] XX[k,2]<-XX[k,1]+trees[[j]]$maps[[i]][k] } } for(k in 1:nrow(YY)){ lower<-which(XX[,1]<=YY[k,1]); lower<-lower[length(lower)] upper<-which(XX[,2]>=(YY[k,2]-tol))[1]; AA<-0 names(lower)<-names(upper)<-NULL for(l in lower:upper) AA<-AA+(min(XX[l,2],YY[k,2])-max(XX[l,1],YY[k,1]))/(YY[k,2]-YY[k,1])*as.numeric(rownames(XX)[l]) ZZ[k]<-ZZ[k]+AA/length(trees) } } tree$maps[[i]]<-YY[,2]-YY[,1] names(tree$maps[[i]])<-round(ZZ*1000) } cols<-rainbow(1001,start=0.7,end=0); names(cols)<-0:1000 if(is.null(legend)) legend<-0.5*max(H) if(is.null(fsize)) fsize<-c(1,1) if(length(fsize)==1) fsize<-rep(fsize,2) if(legend){ if(legend>max(H)){ message("legend scale cannot be longer than total tree length; resetting") legend<-0.5*max(H) } layout(c(1,2),heights=c(0.92,0.08)) plotSimmap(tree,cols,pts=F,lwd=3,fsize=fsize[1],mar=c(0,0.1,0.1,0.1)) X<-cbind(0:1000/1001,1:1001/1001)*(legend/max(H))*(1-fsize[1]*max(strwidth(tree$tip.label))) Y<-cbind(rep(0,1001),rep(0,1001)) par(mar=c(0.1,0.1,0,0.1),xpd=NA) plot(NA,xlim=c(0,1),ylim=c(-0.3,0.3),xaxt="n",yaxt="n",bty="n") for(i in 1:1001) lines(X[i,],Y[i,],col=cols[i],lwd=3,lend=2) text(x=0,y=0,"0",pos=3,cex=fsize[2]) text(x=(legend/max(H))*(1-fsize[1]*max(strwidth(tree$tip.label))),y=0,"1",pos=3,cex=fsize[2]) text(x=(legend/max(H))*(1-fsize[1]*max(strwidth(tree$tip.label)))/2,y=0,"PP(state=1)",pos=3,cex=fsize[2]) text(x=(legend/max(H))*(1-fsize[1]*max(strwidth(tree$tip.label)))/2,y=0,paste("length=",round(legend,3),sep=""),pos=1,cex=fsize[2]) } else plotSimmap(tree,cols,pts=F,lwd=3,fsize=fsize[2]) }