Discrepant results from fitContinuous(...,method="lambda") and phylosig()

I just received the following email query (edited for brevity):

"I encountered different results of lambda using your phylosig code and geiger's fitContinuous for some of my datasets.... I'm sending one of my datasets and the phylogenetic tree, and below are the codes that I've used to produce the results...."

I downloaded his code and dataset, and, indeed, the results appeared to be different. In particular:

> phylosig(tree,x,method="lambda")
$lambda
[1] 0.9527932
$logL
[1] -255.0452

while:

> fitContinuous(tree,x,model="lambda")
Fitting lambda model:
$Trait1
$Trait1$lnl
[1] -259.0189
$Trait1$beta
[1] 20
$Trait1$lambda
[1] 0.2444363
...

Uhoh. My first thought, naturally, was that the beta version of my phylosig() function has an error (and indeed it still might). But then looking a little more closely at the result of fitContinuous() I started to become a little suspicious. In particular, I noticed that $Trait1$beta (the fitted value of σ² - the evolutionary rate for this model) was exactly 20. This led me to suspect that fitContinuous() probably does its multivariate optimization using bounds that did not include the MLE of σ² for these data. If so, this should be resolved by rescaling the data as λ is scale insensitive. Indeed, it did:

> phylosig(tree,x/10,method="lambda")
$lambda
[1] 0.9527932
$logL
[1] -135.3108

> fitContinuous(tree,x/10,model="lambda")
Fitting lambda model:
$Trait1
$Trait1$lnl
[1] -135.3108
$Trait1$beta
[1] 0.661207
$Trait1$lambda
[1] 0.9527761
...

Although it is not clear what bounds are used by default in fitContinuous(), it is straightforward to adjust these. In this case, we would just do:

> fitContinuous(tree,x,model="lambda",bounds=list(beta=c(0,100)))
Fitting lambda model:
$Trait1
$Trait1$lnl
[1] -255.0452

$Trait1$beta
[1] 66.11927

$Trait1$lambda
[1] 0.9527657
...

Code for my phylosig() function is of course available on my website.

Animation of branching random diffusion

I just created a function to animate branching random diffusion (i.e., Brownian motion with speciation) in discrete time. The functions (optionally) calls several functions from the R package {animation}, and to save the resulting video to file, requires the installation of 'FFmpeg'. My function is available on my R-phylogenetics page (direct link here).

To run it after having installed {animation} and 'FFmpeg', first load the function from source:

> source("branching.diffusion.R")

then, try with the default options (but saving the animation to file using the filename record="diffusion.mp4"):

> branching.diffusion(record="diffusion.mp4")

Be warned, all the animation functions eat up lots of memory. If you just want to view the animation in R without saving it, simply set record=NULL (the default), i.e.:

> branching.diffusion()

The user can also change the variance of the BM process, the birth-rate, and the time of the simulation. One option, smooth, tells the function whether to destroy the existing plot and create a new one each generation using plot() (which results in a slower animation, but is smooth); or add the new bits of evolution using lines() (which results in a faster but choppy animation - although the saved movie should be unaffected). You can also specify the pause between generations (default, 0.02s) and the path to ffmpeg.exe (default path="C:/Program Files/ffmpeg/bin/ffmpeg.exe"). The video below was created using the function on its default settings. Enjoy!

Comment on evol.rate.mcmc()

Yesterday a user reported that when he tried to reproduce my example analysis for the function evol.rate.mcmc() here, he encountered the following error when pre-processing the posterior sample:

> mcmc<-posterior.evolrate(tree,minsplit,result$mcmc[201:1001,], result$tips[201:1001])
Error in if (where == 0 || where == "root") where <- ROOTx :
missing value where TRUE/FALSE needed

Today I suggested he try updating {ape} - he did, and the problem went away.

I'm not sure what the compatibility issue is here, but I suspect it has something to do with one of drop.tip(), extract.clade(), or bind.tree(), as these are all called extensively in this function. Take home message - update {ape}!

To find out what version of {ape} (or any other package) you are working with, simply type:

> packageDescription("ape")["Version"]
$Version
[1] "2.7-1"

(FYI 2.7-1 is in fact the version I have installed on this machine.)

Update to phylosig() - important for large trees

Diego Bilski just emailed me with the following text regarding a problem with my phylosig() function for the estimation of phylogenetic signal:

"With more than about 150 taxa, the lambda calculations of lnlik and p values are not possible, but I wasn't able to identify what causes it. With fewer taxa there are no problems. Even with a simple simulation the same 23 warnings are reported:
> library(ape)
> tree<-rtree(200)
> x<-rTraitCont(tree)
> s.lambda<-phylosig(tree,x,method="lambda",test=TRUE)
There were 23 warnings (use warnings() to see them)"

When I tried to reproduce Diego's error, I could do so easily:

>require(ape)
> source("phylosig.R")
> tree<-rtree(200)
> x<-rTraitCont(tree)
> test<-phylosig(tree,x,method="lambda")
There were 14 warnings (use warnings() to see them)
> test
$lambda
[1] 1.091754
$logL
[1] 87.0783
> warnings()
Warning messages:
1: In optimize(f = likelihood, interval = c(0, maxLambda), ... :
NA/Inf replaced by maximum positive value
2: In ...

I quickly suspected that the problem was with the following calculation in the expression for the log-likelihood:

log(det(sig2*Cl))

This is the calculation of the logarithm of the determinant of the evolutionary rate × the lambda-transformed VCV matrix for the tree. The problem is due to the fact that for moderately large trees det(sig2*vcv(tree) can very easily escape the range of numerical precision of R (i.e., it will become infinitesimal or very large) and thus will be set to 0 or Inf respectively. This will make log(det(...)) Inf or -Inf.

Fortunately, there is a second function to compute the matrix determinant in R, determinant(), and the advantage of this function is that it has an option logarithm, which if set to TRUE returns the logarithm of the determinant. This is perfect for us as we are computing the log-likelihood anyway and this will allow us to compute the log-determinant for matrices in which log(det(...)) will not evaluate.

Indeed, substituting determinant(...,logarithm+TRUE) for log(det(...)) seems to solve our problem from before:

> source("phylosig.R") # updated version
> test<-phylosig(tree,x,method="lambda")
> test
$lambda
[1] 0.9766446
$logL
[1] 156.8490

The new version of phylosig() is now available on my R-phylogenetics page (direct link here).

Bug fix and update to min.split()

I just fixed and updated my function min.split(). This function can be used to analyze the posterior sample of rate shift positions generated by evol.rate.mcmc() (described here: 1,2,3,4,5,6; and in an upcoming paper).

The bug caused the distance between two points on the same edge to be miscalculated. Oops! This should have been the easy part! Any users of evol.rate.mcmc() who have used min.split() should definitely update to v0.5 (here).

The update now allows the user to find the split in a sample that minimizes the summed or the summed squared distances to all the other points in the sample. It can also (optionally) print the distance matrix among splits.

The way the function works is as follows. Take the following set of points on the tree below:


> tree<-read.tree(text="(((1:1.0,2:1.0):1.0,3:2.0):1.0,4:3.0);")

which we represent in a list as follows:

> nodelist<-matrix(c(1,2,7,6,3,4,4,0.5,0.5,0.5,0.5,1.5,1,2), 7,2,dimnames=list(NULL,c("node","bp")))
> nodelist
node bp
[1,] 1 0.5
[2,] 2 0.5
[3,] 7 0.5
[4,] 6 0.5
[5,] 3 1.5
[6,] 4 1.0
[7,] 4 2.0

The syntax here is simple. "node" is the node number in $edge; and "bp" is just the distance along the preceding edge from the node below (not from the "node").

We can then run min.split() as follows:

> source("evol.rate.mcmc.R")
> min.split(tree,nodelist,printD=TRUE)
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 0.0 1.0 1.0 2.0 3.0 3.5 4.5
[2,] 1.0 0.0 1.0 2.0 3.0 3.5 4.5
[3,] 1.0 1.0 0.0 1.0 2.0 2.5 3.5
[4,] 2.0 2.0 1.0 0.0 2.0 1.5 2.5
[5,] 3.0 3.0 2.0 2.0 0.0 3.5 4.5
[6,] 3.5 3.5 2.5 1.5 3.5 0.0 1.0
[7,] 4.5 4.5 3.5 2.5 4.5 1.0 0.0
$node
7
$bp
0.5

Which returns the minimum split (in this case node=7,bp=0.5), and, with printD=TRUE, prints to screen the distance matrix among points in the row order of nodelist.

To compute the summed squared distances instead of the summed distances, we just need to do the following:

> min.split(tree,nodelist,method="sumsq")
$node
6
$bp
0.5

and that's it!

Plotting the posterior probabilities on the tree (using pies!)

I thought I would post a few words on plotting the posterior probabilities of the rate shift on the tree. This follows specifically from my previous post here.

We can, of course, do this using edgelabels() from {ape} and plot the probabilities in numeric form (boring). We can also use pie graphs! This makes more sense here than for, say, plotting posterior clade probabilities from Bayesian phylogeny inference because in our MCMC method the rate shift must occur somewhere in the tree, so the posterior probabilities of it being anywhere have to sum to 1.0. (Consequently most of the posterior probabilities will often be close to zero and can be left off - which helps us avoid being overwhelms by pies!)

Picking up our simulated example, from before, we can plot the posterior probabilities (for all p>0.01) of the rate shifts on the tree using the following commands:

> p<-hist(mcmc[,"node"], breaks=0.5:(length(tree$tip)+tree$Nnode+0.5))
> plot(tree)
> edgelabels(edge=match(which(p$density>0.01),tree$edge[,2]), pie=p$density[which(p$density>0.01)],piecol=c("black","white"))

Which should produce something close to the tree below:



Note, that we have simulated a very large effect on the evolutionary rate, which explains why the posterior density of the rate shift is so highly concentrated on the correct edge. For empirical data, the posterior density will often be more dispersed across the branches of the tree.

Update to ltt() and calculation of the gamma-statistic

I just updated my lineage-through-time calculation & plotting function ltt(). The updated code is available on my R-phylogenetics page (direct link here). This function does not add much to what is already available with Dan Rabosky's great {laser} package for diversification analyses; but it does allow one to create an LTT plot with extinct tips, which I do not believe is permitted in prior implementations.

I fixed a couple of bugs in this update. The most significant of these created a bug for any tips with zero length. The second, caused an error if the user told the function to drop extinct tips, but in fact all tips in the tree were contemporaneous.

In addition to these bug fixes, I also just added the optional computation of Pybus & Harvey's (2000) γ statistic, along with a two-tailed P-value for this test.

Load and run as follows (with a simulated tree from {geiger}'s birthdeath.tree() as an example):

> source("ltt.R")
> require(geiger)
> tree<-birthdeath.tree(b=1,d=0,taxa.stop=100)
> ltt(tree)
$ltt
    101 102 103 112 ... 
  1   2   3   4   5 ...
$times
                101       102       103       112  ...
0.0000000 0.0000000 0.3710933 0.7407788 1.0535915  ...
$gamma
[1] -0.2348588
$p
[1] 0.8143183