Showing ancestral state on a phylo.to.map plot

Today I received the following inquiry:

“I am doing phylogeographic analyses of pathogenic bacteria in …, and have found some of the visualization methods helpful for understanding outbreak dynamics. Now here comes the question or maybe a feature request: Is it possible to combine ”plot(summary(simmap.trees)))“ with the ”phylo.to.map(phylo.tree)“ function, to project the ancestral state reconstruction to a world map?”

Indeed, this is already possible because nodelabels( ) & tiplabels( ) from the ape package play perfectly well with phylo.to.map.

One very important word of caution is required, however. Because phylo.to.map( ) does node rotations by default, it will be necessary to match the nodes from the stochastic character mapped trees to the nodes of the rotated phylogeny in the phylo.to.map object.

The easiest way to do that is just to: (1) run phylo.to.map, then; (2) pull out the tree from the resultant "phylo.to.map" object & use that tree for stochastic mapping.

Here's how that might look (using simulated data):

library(phytools)
## here's our tree & data
World.tree

## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  W.ubvct, N.crk, O.tigbryo, A.uxgbavw, S.igzt, M.wuvrb, ...
## 
## Rooted; includes branch lengths.

World

##                      lat        long
## W.ubvct      -58.1670987 -107.067702
## N.crk        -62.4831489 -163.553854
## O.tigbryo    -19.6762651  -65.235665
## A.uxgbavw     -4.4481564  -15.348011
## S.igzt         0.4478175   -7.230006
## M.wuvrb      -24.4885168   76.791949
## I.sfgjcqyphr -27.4762961   90.845604
## P.wtlcrszkfi -26.8759691   56.549824
## Y.kyaltcfq   -23.0423582   71.257658
## N.ojt        -15.6057166   79.638752
## Y.kuvywit    -29.7095111   77.569301
## M.jcfvrbq    -33.2212398   59.108004
## G.pbhsac     -11.1126509   73.307996
## L.tcpsneijl   -2.6186920   60.174359
## A.gprdmcbq    -1.9270332   53.704890
## F.qfpenc      -4.6062239   15.663828
## L.erq         -9.9327773   14.033582
## W.fxeus       75.6131818  -24.214185
## E.qbpilzy     48.7607643 -101.419331
## I.fzuryvaxs   22.6740158  -31.775542
## G.ctas        37.4639540    7.982519
## X.oflhdrym    40.3753670  -15.282336
## Y.baevguzk    35.5968790   -3.151682
## Z.dczjkolt    -6.9465315 -132.877519
## W.uwmjvdal     7.4247476 -132.126099
## A.vnosymup    -0.9018462  -58.811375

x ## discrete character

##      W.ubvct        N.crk    O.tigbryo    A.uxgbavw       S.igzt 
##            b            b            a            a            a 
##      M.wuvrb I.sfgjcqyphr P.wtlcrszkfi   Y.kyaltcfq        N.ojt 
##            b            b            b            b            b 
##    Y.kuvywit    M.jcfvrbq     G.pbhsac  L.tcpsneijl   A.gprdmcbq 
##            b            b            b            b            b 
##     F.qfpenc        L.erq      W.fxeus    E.qbpilzy  I.fzuryvaxs 
##            b            a            a            a            a 
##       G.ctas   X.oflhdrym   Y.baevguzk   Z.dczjkolt   W.uwmjvdal 
##            a            a            a            b            b 
##   A.vnosymup 
##            b 
## Levels: a b

## first phylo.to.map
map.object<-phylo.to.map(World.tree,World,plot=FALSE)

## objective: 226

## objective: 172

## objective: 148
## objective: 148
## objective: 148
## objective: 148
## objective: 148
## objective: 148
## objective: 148
## objective: 148

## objective: 146
## objective: 146
## objective: 146
## objective: 146
## objective: 146
## objective: 146
## objective: 146
## objective: 146
## objective: 146
## objective: 146
## objective: 146
## objective: 146
## objective: 146

## objective: 142
## objective: 142

## now stochastic mapping
mtrees<-make.simmap(map.object$tree,x,nsim=100)

## make.simmap is sampling character histories conditioned on the transition matrix
## 
## Q =
##              a            b
## a -0.007975089  0.007975089
## b  0.007975089 -0.007975089
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##   a   b 
## 0.5 0.5

## Done.

## compute node probabilities
smap.object<-summary(mtrees)

## plot
cols<-setNames(c("blue","red"),levels(x))
plot(map.object,ftype="off",pts=FALSE)
nodelabels(pie=smap.object$ace,
    piecol=cols,cex=0.6)
points(World[,2:1],bg=cols[x[rownames(World)]],
    cex=1.5,pch=21)

That's not bad, right?

FYI, the tree & data were simulated as follows:

World.tree<-pbtree(n=26,scale=100)
World.tree$tip.label<-replicate(Ntip(World.tree),
    paste(sample(LETTERS,1),".",
    paste(sample(letters,round(runif(n=1,min=3,max=10))),
    collapse=""),
    sep=""))
lat<-fastBM(World.tree,sig2=10,bounds=c(-90,90))
long<-fastBM(World.tree,sig2=80,bounds=c(-180,180))
World<-cbind(lat,long)
Q<-0.01*matrix(c(-1,1,1,-1),2,2,
    dimnames=list(c("a","b"),c("a","b")))
x<-sim.Mk(World.tree,Q)