Two years of blogging

Just as the title of this post would suggest, phytools blog is two years old** today. (**More or less - as I noted last year, I had a few earlier placeholder posts, but my first serious entry was dated Dec. 24, 2010; and I first advertised this new blog via the now defunct dechronization blog the next day.)

In two years of blogging I have created over 320 posts. Obviously, some entries are very brief - maybe just a few words about a bug fix or report. In others, I was more in depth - for instance, when I wrote about the difficulty of estimating phylogenetic signal on very small trees; or when I discussed how to include fossil phenotypes in the estimation of ancestral states. If printed as a book (and, who knows, perhaps it will one day form the basis for one), the phytools blog would be a dense tome containing over 350 pages. Yikes!

The phytools blog has also grown in popularity over the past year. Raw "pageview" counts can be affected by many things, and are thus difficult to interpret, but phytools blog has received over 90,000 pageviews since it's creation - over 60,000 in the past year. By any measure, that rate is growing - and this also equates to over 160 pageviews/day for 2012. Among the posts accruing the greatest number of pageviews in 2012: my series on xkcd style trees (e.g., 1, 2, 3); my entries on visualizing trait evolution on trees using colors (e.g., 1, 2); and my posts about ancestral character estimation under the threshold model from quantitative genetics (e.g., 1, 2, 3).

The phytools package (also on CRAN) has grown as well. From a small hodgepodge of methods, phytools has grown into a big one - now containing over 80 different functions and a PDF manual that (as of latest CRAN release) contained 87 pages.

Towards the end of 2012, phytools blog got a new look. Since when the original visual theme was created, phytools had no plotting functions (but now it has tons), this seemed appropriate. The blog also got a new URL (after I obtained the domain phytools.org).

I don't doubt for a second that 2013 will be just as productive as 2012 - so thanks for reading, and please keep coming back. Happy holidays!

Specifying pie colors in ancThresh & plotThresh

A phytools user contacted me to express difficulty in controlling pie colors at nodes using the new phytools functions for ancestral character estimation under the threshold model, ancThresh and plotThresh. This is done using the argument piecol, set to a vector with the states for the discrete trait as names, and the desired colors as values (in any format acceptable by R, e.g., here).

Here's a quick demo. I assume that we've already run the MCMC, and mcmc is the object returned by ancThresh. Data are in X:

> cols
   not     mod    high
"green"  "blue"   "red"
> plotThresh(tree,X,mcmc,piecol=cols,tipcol="estimated", label.offset=0.01)

Note that label.offset is on the scale of your phylogeny, so will have to be varied with total tree length. I've been meaning to migrate control of the pie sizes to the user, but have not yet had time to do so. (Plus, I've been sick.)

Plotting densityMap to PDF

A phytools user asks:
In the densityMap graphic you sent me it seemed like you had spread the branches apart some to make the tree more readable. Is that true? If so, could you tell me how?! I have attached a graphic I made and it looks terrible. . . .

Indeed it did look terrible.

In fact, there was really no mystery to how I used the same function to make a "nice" looking figure - I just used the 'grDevices' function pdf to plot directly to a PDF. This means I can set the dimensions of the plot to be anything I like, and squidged up branches and tip labels in an R plotting window will look nice!

Here's a quick example, using contMap (because densityMap is slow for large trees):

> tree<-pbtree(n=150)
> x<-fastBM(tree)
> pdf(file="contMap.pdf",height=14,width=5)
> # now the output of a plotting call will go to my PDF
> contMap(tree,x,outline=FALSE,fsize=c(0.6,1), lims=c(floor(range(x)[1]*10)/10,
ceiling(range(x)[2]*10)/10))
> # you can probably ignoring the "lims" argument
> # I've just used it here to get a nicer legend
> dev.off()
null device
          1

And here's the result (pulled from the PDF - obviously reduced in quality as a result):

New version of threshDIC

I just posted a new version of threshDIC for computing the Deviance Information Criterion from the object returned by ancThresh. The previous version (of threshDIC not ancThresh) just plain doesn't work. I made some mistakes on key calculations that are now fixed.

DIC is in some ways analogous to the Akaike Information Criterion (AIC), but provides the advantage that it can be used on the results of an MCMC chain - and thus doesn't not require that we can analytically or numerically solve for the maximum of the likelihood function.

The deviance is just -2 × the log-likelihood; and DIC requires that we calculate to versions of this - the mean deviance across samples (Dbar), and deviance at the mean parameter values from the posterior sample (Dhat. We penalize our IC by computing the effective number of parameters of the model (pD or pV, depending on the method we want to use).

Here's a quick demo for the dataset and tree I used earlier:

> require(phytools)
> threshDIC(tree,X,mcmc,sequence=c("not","mod","high"))
    Dbar      Dhat        pD       DIC
-21.29781 -68.07082  46.77302  25.47521

Just like with AIC, we should generally prefer the model with the lowest DIC - although DIC differences of less than 5-10 are difficult to interpret.

The new function is here, but I would advise downloading the latest version of phytools (phytools 0.2-15) and installing from source if you want to use this function:

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

Function setNames

I can't believe I just discovered the function setNames in the 'stats' package. All it is is a wrapper to combine the two lines of code required to create a vector & assign its names. (It can also be used to add names to other objects, such as a list.) So, for instance, instead of:

> x<-c(1:3)
> names(x)<-c("one","two","three")
> x
 one   two three
   1     2     3

We can do:

> x<-setNames(c(1:3),c("one","two","three"))
> x
 one   two three
   1     2     3

How have I been living without this for so long!

Faster inversion of square symmetric positive-definite matrices

For square, symmetric, positive-definite matrices (like covariance matrices) there is a method for faster matrix inversion that uses inversion of the Cholesky matrix. This is described here and implemented in the R package "Matrix". Here's a quick example, using the phylogenetic covariance matrix (vcv.phylo):

> require(phytools)
> # simulate tree
> tree<-pbtree(n=2000)
> C<-vcv(tree)
> # using solve
> system.time(Cinv1<-solve(C))
  user  system elapsed
 15.15    0.07   15.23
> system.time(Cinv2<-chol2inv(chol(C)))
  user  system elapsed
  6.69    0.02    6.74
> mean(abs(Cinv2-Cinv1))
[1] 3.953766e-16

That's faster.

Hypothesis testing using fitDiversityModel

A user had this question about using the function fitDiversityModel (which implements an ML version of the method in Mahler et al. (2010) to test hypotheses about phenotypic diversification:

I have one short question. In the paper by L. Mahler & al. they also tested the model with a single evolutionary rate. I’m asking me if it is “enough” to write it like that: Constant.time.results <- fitDiversityModel(tree,x,d=NULL)

The answer is "no" - but it is very easy to conduct the desired test. If fitDiversity(...,d=NULL), then the function estimates lineage density at each node as if the radiation is occurring within a single region (and there is no extinction or missing taxa - but that is a matter for another day); and then fits a model in which the rate of evolution varies as a function of lineage density. This will generally not be our null model. For that, we should specify d to be a vector of zeroes. We can then use a likelihood ratio test (for example) to test the null hypothesis that the pace of evolution has declined with increasing lineage density.

Here's a worked example, simulated so that the null is correct:

> require(phytools)
> # simulate tree & data under the null
> tree<-pbtree(n=100)
> x<-fastBM(tree)
> # fit the lineage-density model assuming a single region
> # (i.e., all contemporary ancestors were competitors)
> fitLD<-fitDiversityModel(tree,x)
no values for lineage density provided; computing assuming single biogeographic region
> fitLD
$logL
[1] -126.7954
$sig0
[1] 1.171704
$psi [1] -0.004116475
$vcv
            sig0           psi
sig0  0.121592376 -1.833989e-03
psi  -0.001833989  3.240988e-05

> # now fit the null hypothesis of constant rate
> d<-rep(0,tree$Nnode); names(d)<-1:tree$Nnode+length(tree$tip)
> fitNull<-fitDiversityModel(tree,x,d)
psi not estimable because diversity is constant through time.
> fitNull
$logL
[1] -127.0581
$sig0
[1] 0.9602251
$vcv
            sig0
sig0 0.01844064

> # now test using LR
> LR<-2*(fitLD$logL-fitNull$logL)
> P<-pchisq(LR,df=1,lower.tail=F)
> P
[1] 0.4685722

The null model is just pure Brownian motion; and we would obtain the same parameter estimate and likelihood using, say, fitContinuous in the geiger package:

> fitContinuous(tree,x)
Fitting BM model:
$Trait1
$Trait1$lnl
[1] -127.0581
$Trait1$beta
[1] 0.9602253
$Trait1$aic
[1] 258.1161
$Trait1$aicc
[1] 258.2411
$Trait1$k
[1] 2

That's it.