Issue with plotting output trees from SIMMAP using densityMap

A phytools user contacted me yesterday with the following issue:

I was trying to use your densityMap function to visualize some output from SIMMAP. The mutational mappings seemed to work fine in SIMMAP. I tried importing into R using read.simmap, this seemed to work fine. I then tried to use densityMap but it keeps giving me the following error:
> trees<-read.simmap(file)
> densityMap(trees,res=1000)
sorry - this might take a while; please be patient
Error in lower:upper: NA/NaN argument

This isn't something that I'd seen before, but I recognized that the code is from a part of the function in which the beginning and ending points for each map segment on a branch (XX) is compared to a matrix containing the beginning and ending points for every step interval along the same edge (YY); in which these steps are based on the edge length of the branch and the function argument res. Importantly, the beginnings and endings are in terms of the height above the root.

This made me suspect that the problem was with the edge lengths of the tree. densityMap assumes that the topology of the input trees are equal with proportional branch lengths. (In fact, there is an option - densityMap(...,check=TRUE) - that evaluates this equality criterion, and returns an error message if it fails.) For some reason that I don't completely understand (and even when 'mutational mapping' is performed on a single input tree), SIMMAP returns a sample of trees with varying total length. Theoretically, this should not be a problem if the branch lengths are strictly proportional. Unfortunately, due to numerical precision (which is relatively low - 6 digit - in the SIMMAP output), the branch lengths of these varyingly scaled trees are not strictly proportional. Consequently, even after they have been scaled to the same total length using phytools function rescaleSimmap (which is done internally in densityMap), the trees have slightly different branch lengths - sufficiently different, evidently, to account for this error.

My eventual workaround was to write a function that would apply one set of branch lengths (say, from a reference tree) to another tree or set of trees. For basic "phylo" objects, this is trivial of course. We would just do:

# for reference tree tr1
tr2$edge.length<-tr1$edge.length

However, for trees with mappings, it is slightly more complicated. Here is a simplified version of the function:

applyBranchLengths<-function(tree,edge.length){
  if(class(tree)=="multiPhylo"){
    trees<-lapply(tree,applyBranchLengths,
      edge.length=edge.length)
    class(trees)<-"multiPhylo"
    return(trees)
  } else {
    tree$edge.length<-edge.length
    if(!is.null(tree$maps)){
      for(i in 1:nrow(tree$edge)){
        temp<-tree$maps[[i]]/sum(tree$maps[[i]])
        tree$maps[[i]]<-temp*tree$edge.length[i]
      }
    }
    return(tree)
  }
}

This function uses the proportional time spent in each state along each edge of the target tree in rescaling the mappings to the reference tree. (I've left out calculation of the $mapped.edge matrix here, which involves a few more lines of code.) I have posted the full function code for this here.

Anyway, it worked! Here's a quick demo with a simulated tree & dataset in which I performed stochastic mapping using SIMMAP v1.5 for another project:

> trees<-read.simmap("simmap.trees.nex",format="nexus", version=1.5)
> densityMap(trees)
sorry - this might take a while; please be patient
Error in lower:upper : NA/NaN argument
> # get the tree heights
> h<-sapply(trees,function(x) max(nodeHeights(x)))
> h
 [1] 13.215321  6.213479 16.864350  9.936482 ....
 ...
> # ok, they have different heights - is that the problem?
> # try rescaling
> rescaled<-rescaleSimmap(trees,mean(h))
> densityMap(rescaled)
sorry - this might take a while; please be patient
Error in lower:upper : NA/NaN argument
> # nope!
> # ok, how about *not* rescaling, but applying the
> # relative branch lengths from trees[[1]] to all the
> # trees in the sample
> newbl<-mapply(applyBranchLengths,trees, lapply(as.list(h),"*",trees[[1]]$edge.length/h[1]), SIMPLIFY=FALSE)
> class(newbl)<-"multiPhylo"
> # these trees have their original total lengths
> # but relative branch lengths from trees[[1]]
> densityMap(newbl,fsize=c(0.6,1))
sorry - this might take a while; please be patient

Plotting many stochastically mapped trees

I recently described a new function, densityMap (1, 2), that can be used to plot the posterior probability of a binary trait evolving along the branches of the tree based on the method of stochastic character mapping (Neilsen 2002; Huelsenbeck et al. 2003; Bollback et al. 2006). There's another way to do this, of course - we can just plot them all! Here's a quick demo & the resultant visualization (from 100 stochastic maps):

> # data for the binary [0,1] state in x
> trees<-make.simmap(tree,x,nsim=100)
> par(mfrow=c(10,10))
> cols<-c("blue","red"); names(cols)<-c(0,1)
> plotSimmap(trees,cols,ftype="off",lwd=1,pts=F)

That's pretty cool. Still, I think there is an argument for densityMap:

> par(mfrow=c(1,1))
> densityMap(trees)
sorry - this might take a while; please be patient
Waiting to confirm page change...

Fitting the λ model in phylogenetic regressions using various functions

A R-sig-phylo reader recently asked the following:

I am trying to implement Liam Revell's suggestion on the evaluation of Pagel's lambda simultaneous to fitting PGLS to minimize the effects of wrong model selection (OLS vs. PGLS) on species data (i.e. Revell 2010, Methods in Ecology and Evolution 1: 319-329).Does anyone know if this kind of simultaneous optimization of lambda and PGLS parameters is presently implemented in R? The obvious place would be phytools, but I could not find anything similar in that package.

In fact, this is possible using multiple functions in several packages. Unfortunately, it is not always completely obvious - so the question is (in my opinion) fully justified! Specifically, this can be done using nlme::gls, caper::pgls, and (indeed) phytools::phyl.resid. Here's my summary (modified from my response) of how (which I repeat here in the hope that it will help folks in future when they are searching for information on this topic):

# using ape & nlme
library(ape); library(nlme)
# assuming your data are contained in named vectors x & y
y<-y[names(x)]
fit<-gls(y~x,data.frame(x,y),correlation=corPagel(1,tree, fixed=FALSE),method="ML")
# assuming your data are in data frame X, with column names
# var1, var2, etc.; for example:
fit<-gls(var1~var2+var3,X,correlation=corPagel(1,tree, fixed=FALSE),method="ML")

# using caper
library(caper)
# assuming your data are in named vectors x & y
y<-y[names(x)]
X<-data.frame(x,y,Species=names(x))
fit<-pgls(y~x,comparative.data(tree,X,"Species"), lambda="ML")

# using phytools
library(phytools)
# assuming your data are in named vectors x & y
fit<-phyl.resid(tree,x,y,method="lambda")

I would generally recommend using gls & pgls over phyl.resid, because they "play well" with generic functions such as summary and residuals, with the following caveats (1) phyl.resid may be a little easier to use (because it takes standard arguments, no functions or formulae); and (2) phyl.resid (for some reason) seems to run a little faster, especially for large trees. (For instance, it took only 45s on a tree with 1000 tips while gls took over a minute and pgls over 4 minutes. This could be important for investigators repeating the same analysis on many trees - for instance, a sample from the posterior distribution of a Bayesian phylogeny inference run.)

Update to function matchNodes

I'm working on a manuscript review this Saturday and (I don't think I'm giving away too much about the article to say that) I needed a function to match the node numbers of a tree with the same nodes in a re-rooted version of the same phylogeny. This differs from the purpose of my previous phytools utility function, matchNodes, which matches nodes based on the commonality of their descendant tips. This method will not, in general, produce the correct set of matchings for a tree that has merely been rerooted at an internal node.

My idea for obtaining the correct matchings was to use the set of distances to the tip species in the tree. Matching nodes with the same set of patristic distances (or, for trees with arbitrary scale, proportional distances) must be equivalent. Here is the core code of the function:

D1<-dist.nodes(tr1)[1:length(tr1$tip), 1:tr1$Nnode+length(tr1$tip)]
D2<-dist.nodes(tr2)[1:length(tr2$tip), 1:tr2$Nnode+length(tr2$tip)]
rownames(D1)<-tr1$tip.label
rownames(D2)<-tr2$tip.label; D2<-D2[rownames(D1),]
for(i in 1:tr1$Nnode){
  z<-apply(D2,2,function(X,y) cor(X,y),y=D1[,i])
  Nodes[i,1]<-as.numeric(colnames(D1)[i])
  Nodes[i,2]<-as.numeric(names(which(z>=(1-tol))))
}

The first part of the code uses the 'ape' function dist.nodes to compute the set of distances from each node to all tips; and then sorts them to have the same row-wise ordering in the two matrices. Next, we go through the columns of D1 and identify the column of D2 with proportional or equal distances.

Of course - some important bookkeeping has been omitted here. For the full code, see utilities.R on the phytools page.

Let's try it out:

> tr1<-pbtree(n=20)
> par(mfcol=c(1,2))
> plotTree(tr1,node.numbers=T)
> # re-root tr1 at node=33
> tr2<-root(tr1,node=33)
> plotTree(tr2,node.numbers=T)

> matchNodes(tr1,tr2,method="distances")
      tr1 tr2
[1,]  21  NA
[2,]  22  30
[3,]  23  31
[4,]  24  32
[5,]  25  29
[6,]  26  25
[7,]  27  26
[8,]  28  27
[9,]  29  28
[10,]  30  24
[11,]  31  23
[12,]  32  22
[13,]  33  21
[14,]  34  38
[15,]  35  33
[16,]  36  34
[17,]  37  35
[18,]  38  36
[19,]  39  37

> # we can also reverse the order of the trees
> matchNodes(tr2,tr1,method="distances")
      tr1 tr2
[1,]  21  33
[2,]  22  32
[3,]  23  31
[4,]  24  30
[5,]  25  26
[6,]  26  27
[7,]  27  28
[8,]  28  29
[9,]  29  25
[10,]  30  22
[11,]  31  23
[12,]  32  24
[13,]  33  35
[14,]  34  36
[15,]  35  37
[16,]  36  38
[17,]  37  39
[18,]  38  34

Careful inspection of the matchings above should (I hope) reveal that the nodes from tr1 have been correctly matched to the corresponding nodes of tr2 (in the first example - and the converse in the second case). The column headings in the second example reflect the input order of the trees - not their names in memory.

Note that case 1 leaves one unmatched node - the root. In case two there are no unmatched nodes because all of the n - 2 nodes in tr2 have a compatriot in tr1.

Assuming no major errors are identified, this function will be in the next version of phytools.

That's it.

New minor phytools version 0.2-11

I just posted a new minor package update (phytools 0.2-11). It can be downloaded here and installed from source. This new version has a small number of significant updates and new functions:

1. A new function (splitplotTree, also see inset) to split a plotted tree into multiple columns or plotting windows, along with one minor bug fix.

2. A couple of important bug fixes for phyl.RMA, my function for phylogenetic RMA (1, 2).

3. Finally, a new function, bind.tip, for attaching a new species to the tree (1, 2). If no terminal edge length is provided the function automatically checks to see if the tree is ultrametric and (if it is) scales the terminal edge length so that ultrmetricity is maintained.

To install from source, just download the build and then type:

> install.packages("phytools_0.2-11.tar.gz",type="source", repos=NULL)

at the R prompt. Good luck!

Bug fix for phyl.RMA

A user pointed out that phyl.RMA did not work (specifically, returned the incorrect signed RMA slope) when the relationship between x & y is negative. This is, in fact, correct. I have now fixed this problem by changing the following line of code:

beta1<-sqrt(temp$R[2,2]/temp$R[1,1])

to:

beta1<-sign(temp$R[1,2])*sqrt(temp$R[2,2]/temp$R[1,1])

the difference being that in the latter case I pre-multiply the evolutionary variance ratio by the sign of the interspecific covariance. Updated code is here and will be in the next version of phytools.

This update introduces the problem that our statistic for hypothesis testing the RMA regression slope, based on equation [8] of McArdle (1988), is only defined for positive RMA regression slopes. To address this, I have boldly assumed that for h₀: β < 0 the null distribution is symmetric about 0 & consequently if b & β have the same sign, I use equation [8] on the absolute values of b & β - otherwise I set T to zero & return a warning. This (at least) has the desirable effect of using McArdle equation [8] without modification with b and β are both > 0 (typical of hypothesis tests about allometry).

Note that, unlike in LS regression, for RMA we cannot normally test h₀: β = 0. This is because β = 0 implies that σ_y = 0.

Testing for Pagel's λ < 1.0

A recent commenter on the phytools blog asks:
I’m working on a comparative phylogenetic project and want to show that my data has phylogenetic signal. Using Pagel’s lambda I get lambda – 0.854. Is this value significant enough in that it’s safe to model the variable as Brownian motion?

I interpret this to mean "how can we test the alternative hypothesis that Pagel's λ is less then 1.0; against the null that λ=1.0." This is easy if we accept the contention that our likelihood ratio should be distributed as a Χ² with 1 degree of freedom (for one parameter, λ, estimated in the more complex model), an asymptotic property of likelihoods. (There may be some very good reasons - which we'll set aside for the moment - to question this assertion, for instance see a recent article by Carl Boettiger et al.)

This can be done most naturally using the fitContinuous function in 'geiger' - but we can also do it using entirely phytools functions, and since I developed and maintain phytools (not geiger), I'll demo that option.

First, let's simulate a tree & data for the demo, using the λ=0.854 obtained by the commenter. We will use geiger::lambdaTree for this part:

> library(phytools); library(geiger)
> # simulate tree
> tree<-pbtree(n=50)
> # simulate data with lambda=0.854
> x<-fastBM(lambdaTree(tree,0.854))

OK, now let's fit both models, compute the likelihood ratio, and then compare it to a Χ² distribution with one degree of freedom. We need to use a little bit of a work-around to fit the BM model using brownie.lite:

> # fit lambda model
> lambda<-phylosig(tree,x,method="lambda")
> lambda
$lambda
[1] 0.9497926
$logL
[1] -87.54412
> # get likelihood for bm model
> tree$mapped.edge<-matrix(tree$edge.length, nrow(tree$edge),1,dimnames=list(NULL,"1"))
> bm.logL<-brownie.lite(tree,x)$logL1
> # conduct hypothesis test using chi-square
> LR<--2*(bm.logL-lambda$logL)
> P<-pchisq(LR,df=1,lower.tail=F)
> P
[1] 0.1346832

So, in this simulated case at least, we cannot distinguish our estimated phylogenetic signal from Brownian motion.

An additional word of caution, however. Sampling error in the estimation of species means can also have the effect of depressing phylogenetic signal - so if the species means are estimated with a limited number of samples per species, it might be wise to conduct this analysis including estimated within species sampling variances. These can be supplied using the argument se in both phylosig and brownie.lite.