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