Thursday, July 25, 2013

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.

3 comments:

  1. Hi Liam

    This function is great!, thanks a lot for posting it. I have been using it, however I added a different map because I only need to plot points for South America. The map I am using is an altitude layer downloaded from the wordlclim raster layer database. However I am having trouble with this map, it seems there is not enough space for the tree to appear above the map, but I am not sure how to re-scale it to leave space for the tree in the upper half, any ideas on how to do this? I am a very basic user.

    Best

    ReplyDelete
    Replies
    1. Hi Francisco.

      You could try adjusting the arguments split and fsize. split is a numeric vector with 2 elements that should add to 1.0. By increasing the first element and shrinking the second, you will allow more area to plot the tree. fsize adjusts the font size used for the tip labels.

      If neither works, please feel free to send me your files as well as detailed information about the code that triggered your error, and I will investigate this further.

      All the best, Liam

      Delete
    2. Hi Liam, many thanks for your reply!

      Delete