Addendum to plotting a tree linking tips to a global map

This isn't how I need to be spending my time (with an NSF grant due next week), but here's a neat addendum to my earlier post about a function to link the tips of a tree to plotted lat/long locations on a global map.

In that post I completely ignored the fact that all the links tend to cross each other because the tips of the tree have not been rotated to maximize the concordance of their order with the longitudinal position of the corresponding coordinates. Here is a function that attempts to do this via one pre-order tree traversal:

minRotate<-function(tree,x){
  tree<-reorder(tree)
  nn<-1:tree$Nnode+length(tree$tip.label)
  x<-x[tree$tip.label]
  for(i in 1:tree$Nnode){
    tt<- tt<-read.tree(text=write.tree(rotate(tree,nn[i])))
    oo<-sum(abs(order(x[tree$tip.label])-
     1:length(tree$tip.label)))
    pp<-sum(abs(order(x[tt$tip.label])-
     1:length(tt$tip.label)))
    if(oo>pp) tree<-tt
    message(paste("objective:",min(oo,pp)))
  }
  return(tree)
}

All that's going on in this incredibly simple function is that the function climbs from the root up through the nodes of the tree. It rotates each one and then checks to see if the objective function (the match in ordering between the tips of the tree and longitudinal coordinates in vector x) improves. It accepts only changes that decrease the value of the objective function.

This strikes me as a search algorithm that would have an incredible propensity to get stuck on local optima - but let's not speculate, let's give it a try. A new option in phylo.to.map optionally uses this optimization function internally.

> library(phytools)
> library(rworldmap)
### Welcome to rworldmap ###
...
> source("phylo.to.map.R")
> # simulate tree & data
> tree<-pbtree(n=26,scale=100)
> tree$tip.label<-LETTERS[26:1]
> lat<-fastBM(tree,sig2=10,bounds=c(-90,90))
> long<-fastBM(tree,sig2=80,bounds=c(-180,180))
> # first, unrotated
> phylo.to.map(tree,cbind(lat,long),rotate=FALSE)

> # now rotated
> phylo.to.map(tree,cbind(lat,long))
objective: 118
objective: 112
...
objective: 96

Well, that may not be perfect - but it is clearly a heck of an improvement. Let's try it again to make sure it wasn't a fluke:

> # simulate tree & data
> tree<-pbtree(n=26,scale=100)
> tree$tip.label<-LETTERS[26:1]
> lat<-fastBM(tree,sig2=10,bounds=c(-90,90))
> long<-fastBM(tree,sig2=80,bounds=c(-180,180))
> # first, unrotated
> phylo.to.map(tree,cbind(lat,long),rotate=FALSE)

> # now rotated
> phylo.to.map(tree,cbind(lat,long))
objective: 280
objective: 136
...
objective: 96

Cool. It can't be any worse than leaving the tree unrotated, anyway.