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:

> 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.

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.

## No comments:

## Post a Comment