Sunday, February 1, 2015

Bug fix for pbtree with user-supplied tip labels and extinction

I just posted a small bug fix in the phytools function pbtree.

pbtree (originally short for Pure-Birth tree) does birth-death phylogeny simulations for a variety of conditions (time-stop, taxon-stop, both together via rejection, and discrete & continuous time simulation). If the taxon-stop criterion is being used (that is, if the function is being told to return a tree with a specific number of extant terminal taxa), then the user is also allowed to supply his or her own tip labels which will be used to populated tree$tip.label in the returned tree. Since a stochastic birth-death tree with N extant tips can have an arbitrary and unknowable number of extinct tips, these are labeled X1, X2, etc. For instance:

library(phytools)
tree<-pbtree(n=26,b=1,d=0.3,tip.label=letters)
## Warning: only using labels in tip.label for extant tips.
##          extinct tips will be labeled X1, X2, etc.
plotTree(tree,ftype="i",fsize=0.9)

plot of chunk unnamed-chunk-1

Unfortunately, the way this works is by first simulating the tree with arbitrary internal labels, and then assigning all N extant tips the user supplied labels after running the phytools function getExtant internally. Due to numerical precision issues, though, if the total tree length is very long (for instance, because b and d are small for a given taxon-stop, n), then (unless its tolerance argument is adjusted) getExtant will malfunction and pbtree will behave as if the entire tree consists of extinct taxa! So, for example:

set.seed(10)
tree<-pbtree(n=26,b=0.001,d=0.0003,tip.label=letters)
## Warning: only using labels in tip.label for extant tips.
##          extinct tips will be labeled X1, X2, etc.
plotTree(tree,ftype="i",fsize=0.9)

plot of chunk unnamed-chunk-2

Clearly, 26 tips are extant in this tree; however the function is behaving as if they are all extinct!

The 'fix' that I've applied is to make the tolerance value in getExtant a function of the total tree length, rather the default value of getExtant.

detach("package:phytools",unload=TRUE)
install.packages("phytools_0.4-44.tar.gz",type="source")
## Installing package into 'C:/Users/Liam/Documents/R/win-library/3.1'
## (as 'lib' is unspecified)
## inferring 'repos = NULL' from 'pkgs'
library(phytools)
packageVersion("phytools")
## [1] '0.4.44'
set.seed(10)
tree<-pbtree(n=26,b=0.001,d=0.0003,tip.label=letters)
## Warning: only using labels in tip.label for extant tips.
##          extinct tips will be labeled X1, X2, etc.
plotTree(tree,ftype="i",fsize=0.9)

plot of chunk unnamed-chunk-3

The updated version of phytools is here. I'm also planning to submit a new version to CRAN in the not-too-distant future. (Fingers crossed!)

That's it!

No comments:

Post a Comment

Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.