Tuesday, March 15, 2011

Bug fixed in read.simmap()

Dave Collar helpfully identified a bug in v0.2 of read.simmap() that was preventing the function from properly reading in a dataset of multiple trees. [Since this function was a complete rewrite of v0.1, this was not a problem in the earlier edition.] This bug has now been fixed and the code (v0.3) is available on my R-phylogenetics page (direct link to code here).


  1. Here's a question, Liam. In the $maps list, I get that each element represents an edge. But, for each edge, what is the order in which the states are given--ancestor to descendant? I want to know if I can take the first state reported (i.e., tree$maps[[i]][1] for edge i) as the state at the ancestral node.

  2. Hi Dave. The elements in $maps are vectors. The order of the vectors corresponds to the order of the edges in $edge. The names of the vectors are the states. So, for instance:
    > tree$maps[[98]]
    1 0 1
    0.0005835831 0.1059972682 0.0239649458
    indicates that the state along edge 98 switches from 1 to 0 to 1 (with corresponding times indicated).
    To get the state at the root, you should just take:
    > names(tree$maps[[1]])[1]
    [1] "1"
    (i.e., "1" in this case). This is because the first edge in the tree always arises out of the root node, and the first state along this edge will be the state at the root node of the tree.
    For a "multiPhylo" object just switch tree[[i]] for tree to get the ith tree.