Plotting point-size on a phylomorphospace proportional to the value of a third character

Recently a phytools user asked if it were possible to scale the plotted points at tips & nodes of a plotted phylomorphospace in proportion to the value of a third trait.

“What I would like to do is to plot axes 1 and 2 in two dimensions and then make the points in the graph as circles proportional to axis 3. Do you have an idea how to do that?”

They provided an attached example in the form of Briggs et al. (1992).

This is actually pretty easy, as follows:

library(phytools)
tree

## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  A, B, C, D, E, F, ...
## 
## Rooted; includes branch lengths.

X

##         [,1]        [,2]       [,3]
## A  1.4800637  1.09652745 -1.1494221
## B  0.4791366  0.40907306 -0.3548262
## C  1.2351474  0.19028043 -0.5856361
## D  2.0851866  1.57117881 -1.3810955
## E  2.1138250  1.55440915 -1.6382140
## F  2.5633292  1.16354716 -0.3219158
## G  2.4537429  1.23287364 -0.3740234
## H  3.5660959 -1.62714583  1.2958532
## I  1.9329064 -1.87436079  1.1086486
## J  0.5050758 -1.56200859  1.9635014
## K  3.7305217 -2.41850413  2.3232728
## L  1.8526637 -2.30833059  3.4108897
## M  1.7070475 -2.10857785  2.2998548
## N  1.4046043 -2.53014781  2.3675294
## O  1.6555163 -1.99505382  2.2816972
## P  1.5965899 -0.51055114  1.0490624
## Q  1.0374652  0.17976541  1.2657477
## R  1.1659834  0.07728329  1.2683668
## S  0.4475435  0.09072082  1.5856455
## T  1.9675116 -0.41992527  1.4377758
## U  0.7077652 -0.56399506  1.6097363
## V -0.2623193 -0.79534615  0.2402987
## W  0.5563390 -1.20883630 -0.4694405
## X  0.1164100 -0.97543773 -0.3125623
## Y -0.6341650 -1.79165841  2.0198135
## Z -0.8576449 -1.74610813  2.0654261

obj<-phylomorphospace(tree,X[,1:2],node.size=c(0,0),xlab="x",
    ylab="y")
a<-fastAnc(tree,X[,3]) ## reconstruct states
cx<-c(X[,3],a)-min(c(X[,3],a))
cx<-cx/max(cx)*2 ## rescale to [0,2]
points(obj$xx,obj$yy,cex=cx,pch=19)

Here I have rescaled the trait values to cex on the interval [0, 2], but this is arbitrary - although we of course have to ensure that no values are negative. We could also easily change colors, etc.

Finally, phytools also has 3D phylomorphospace plotting using 'rgl' or by simulating a third dimension in 2D:

phylomorphospace3d(tree,X,method="static")

That's it.

S3 density method for objects of class "multiSimmap"

I just pushed some new methods S3 methods for objects of class "multiSimmap" - that is a set of stochastic character mapped trees. In particular, I created a density method for computing the posterior distribution of changes on the tree of each type.

This is a job in progress - so far I have implemented this only for binary characters (they don't need to be coded 0 / 1, but not more than two states are allowed). There is also no documentation; however because they are generic methods, R will let us build & install the package anyway.

Here is a quick demo:

library(phytools)

## Loading required package: ape

## Loading required package: maps

packageVersion("phytools")

## [1] '0.5.57'

tree

## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  A, B, C, D, E, F, ...
## 
## Rooted; includes branch lengths.

x

## A B C D E F G H I J K L M N O P Q R S T U V W X Y Z 
## 0 0 0 0 0 0 0 1 0 1 1 0 0 0 0 1 0 0 0 0 1 0 1 0 0 0 
## Levels: 0 1

trees<-make.simmap(tree,x,model="ER",nsim=200)

## make.simmap is sampling character histories conditioned on the transition matrix
## 
## Q =
##            0          1
## 0 -0.9946799  0.9946799
## 1  0.9946799 -0.9946799
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##   0   1 
## 0.5 0.5

## Done.

trees

## 200 phylogenetic trees with mapped discrete characters

obj<-density(trees)

## Loading required package: coda

obj

## 
## Distribution of changes from stochasting mapping:
##  0->1            1->0
##  Min.   :3       Min.   :0
##  Median :7       Median :3
##  Mean   :7.31    Mean   :3.45
##  Max.   :11      Max.   :13
## 
## 95% HPD interval(0->1): [4, 8]
## 95% HPD interval(1->0): [0, 8]

plot(obj)

The density method for objects of class "multiSimmap" also wraps the standard densityMap method. For instance:

obj<-density(trees,method="densityMap",res=200)

## sorry - this might take a while; please be patient

obj

## Object of class "densityMap" containing:
## 
## (1) A phylogenetic tree with 26 tips and 25 internal nodes.
## 
## (2) The mapped posterior density of a discrete binary character with states (0, 1).

plot(obj,lwd=6,outline=TRUE)

This is basically what I have right now. More to come later!

Modified version of Grafen's branch lengths

Last night - based on an R-sig-phylo request I wrote some code to compute a modified version of Grafen's edge lengths, in which the depth of each node is defined by the number of branching events above it - not the number of leaves.

However, this is not what I think the commenter wanted. Indeed, I believe he was looking for a transformation in which the height above the root of all nodes with a maximum path length to any tip of x was standardized for a given value of x.

Here is an example of this. First, our functions:

node.paths<-function(tree,node){
    d<-Descendants(tree,node,"children")
    paths<-as.list(d)
    while(any(d>Ntip(tree))){
        jj<-1
        new.paths<-list()
        for(i in 1:length(paths)){
            if(paths[[i]][length(paths[[i]])]<=Ntip(tree)){ 
                new.paths[[jj]]<-paths[[i]]
                jj<-jj+1
            } else {
                ch<-Descendants(tree,paths[[i]][length(paths[[i]])],
                    "children")
                for(j in 1:length(ch)){
                    new.paths[[jj]]<-c(paths[[i]],ch[j])
                    jj<-jj+1
                }
            }
        }
        paths<-new.paths
        d<-sapply(paths,function(x) x[length(x)])
    }
    paths
}
modified.Grafen<-function(tree,power=2){
    max.np<-function(tree,node){
        np<-node.paths(tree,node)
        if(length(np)>0) max(sapply(np,length)) else 0
    }
    nn<-1:(Ntip(tree)+tree$Nnode)
    h<-sapply(nn,max.np,tree=tree)+1
    h<-(h/max(h))^power
    edge.length<-vector()
    for(i in 1:nrow(tree$edge)) 
        edge.length[i]<-diff(h[tree$edge[i,2:1]])
    tree$edge.length<-edge.length
    tree
}

Now our trees:

library(phytools)
library(phangorn)
t1<-read.tree(text="((((A,B),C),E),(((F,G),H),I));")
t2<-read.tree(text="((((A,B),C,D),E),(((F,G),H),I));")
t3<-read.tree(text="((A,(B1,B2),C),(D,(E1,E2),F,(G1,G2)));")
## standard Grafen tree 1
plotTree(t1,lwd=1)

## modified Grafen
plotTree(modified.Grafen(t1),lwd=1)

## standard Grafen tree 2
plotTree(t2,lwd=1)

plotTree(modified.Grafen(t2),lwd=1)
nodelabels()

Note that nodes #11 & #14 have the same height because they have the same maximum path length to the tips.

## standard Grafen tree 3
plotTree(t3,lwd=1)

plotTree(modified.Grafen(t3,power=3),lwd=1)

Changing power just changes how the tree looks - not the fact that nodes a common maximum number of edges away from the present are at a common height.

Finally, here I show the plot with node labels giving the maximum depth (in nodes) of each internal node of the tree. We can see that with modified.Grafen nodes with the same index have the same height.

tree<-rtree(n=40)
tree$edge.length[sample(1:nrow(tree$edge),40)]<-0
tree<-di2multi(tree)
tree$edge.length<-NULL
max.np<-sapply(1:tree$Nnode+Ntip(tree),function(x,tree)
    max(sapply(node.paths(tree,x),length)),tree=tree)
plotTree(tree,lwd=1,fsize=0.7)
nodelabels(max.np,node=1:tree$Nnode+Ntip(tree),cex=0.7)

plotTree(modified.Grafen(tree,power=3),lwd=1,fsize=0.7)
nodelabels(max.np,node=1:tree$Nnode+Ntip(tree),cex=0.7)

That's all.