Plotting the variance among stochastic character maps on the tree

The phytools function densityMap can be used to visualize the posterior probability of a binary character on the nodes & edges on the tree obtained via stochastic character mapping. By default, this function uses a continuous color gradient from blue to red (with low posterior probability for state 1 in blue); thus purple colors indicate high uncertainty about the state of the character mapped on the tree.

I thought it might be interesting to visualize the variance (i.e., uncertainty) of our ancestral reconstruction directly. We can compute the variance in our posterior probability (which is just a proportion from the posterior sample) as p(1-p). Here is what this looks like:

First, simulate a tree & some data, and then perform stochastic mapping × 100:

> tree<-pbtree(n=30,scale=2)
> Q<-matrix(c(-1,1,1,-1),2,2)
> rownames(Q)<-colnames(Q)<-0:1
> x<-sim.history(tree,Q)$states
> mtrees<-make.simmap(tree,x,nsim=100)

Now we're ready for my somewhat hacky solution for plotting the variance among stochastic maps on the edges & nodes of the tree. Note that this requires the latest phytools build, phytools 0.3-13.

XX<-densityMap(mtrees,plot=FALSE)
aa<-colnames(XX$tree$mapped.edge)
tt<-XX$tree
# collapse each prob with the same variance into each other
for(i in 0:499){
 oo<-c(as.character(i),as.character(1000-i))
 nn<-as.character(i)
 chk<-oo%in%aa
 if(sum(chk)==2) tt<-mergeMappedStates(tt,oo,nn)
 else if(chk[2])
  for(j in 1:length(tt$maps))
   names(tt$maps[[j]])[which(names(tt$maps[[j]])==
   oo[2])]<-oo[1]
}
# compute the variances
p<-as.numeric(names(XX$cols))/(length(XX$cols)-1)
v<-(p*(1-p))[1:501]
XX$tree<-tt
# this is a hack
XX$cols<-setNames(as.vector(rbind(grey(v/0.25), grey(v/0.25))),as.vector(rbind(0:500,0:500)))
XX$lims<-c(0,0.25)
class(XX)<-"contMap"
# plot!
plot(XX)

OK, now let's compare the normal contMap plot with our "uncertainty" plot. First, here is the contMap plot:

Now here is the variance plot:

That's it.