CRAN Task View: Phylogenetics

Brian O'Meara has very kindly added "phytools" to the CRAN Task View: "Phylogenetics, Especially Comparative Methods" (link here). Thanks to him & thanks to Dave Bapst for suggesting it.

New version of evol.vcv()

Ok, I just posted a new version of evol.vcv() - direct link to code here. The main new feature of this version (v0.2) is that it also computes the variances of the parameter estimates. (The standard errors can then just be computed by taking the square roots.) These variances come from the matrix of partial second derivatives of the likelihood surface at the optimum (i.e., the curvature of the surface), also known as the Hessian matrix. The intuitive relationship between the curvature of the likelihood surface and the uncertainty of our parameter estimates is quite straightforward: if the likelihood surface is very flat (low curvature), then the uncertainty in our parameter estimates is large; conversely if the likelihood surface is very peaked at the optimum (i.e., high curvature), then this implies that uncertainty in our parameter estimates is low.

This function (that is, evol.vcv()) implements the method of Revell & Collar (2009) and was also described in a prior blog post.

In R, normally obtaining the Hessian matrix is easy because the built-in optimizer, optim(), optionally returns its value. Unfortunately, as I outlined in a previous post, this would not do for evol.vcv() because in this function I actually optimized the likelihood with respect to the Cholesky matrices rather than the evolutionary variance-covariances themselves. This was done for computational reasons because the evolutionary variance-covariance matrices (called "rate matrices" in our paper) are subject to annoying constraints which mean that for most random values the likelihood cannot be computed.

The solution, which I figured out yesterday and then implemented today, was to optimize with respect to the Cholesky matrices (as before) - but then differentiate the likelihood surface for the full VCV matrices at the optimum using hessian() in the "numDeriv" package (CRAN page here). This seems to work - and this analysis is implemented in v0.2 of evol.vcv() on my R phylogenetics page. This will also be included in the next version of "phytools".

The variances are not computed by default. To get them, load the present version from source:

> source("http://faculty.umb.edu/liam.revell/phytools/evol.vcv/v0.2/evol.vcv.R") # this should work
> result<-evol.vcv(tree,X,vars=TRUE)

for SIMMAP style tree and data matrix X. Good luck, and please report if this function seems to work properly or not.

Finding the standard errors from the likelihood surface in evol.vcv()

Today I've been working on returning the variances of the parameter estimates in evol.vcv(). Just to remind the readers of this blog, this function fits two or more instantaneous variance-covariance matrices for the Brownian evolutionary process to different pre-specified parts of a phylogenetic tree. (For more details, see my 2009 paper with Dave Collar or a previous blog post here.)

The first problem that I've run into is that my function performs optimization of the likelihood as a function of the Cholesky decomposite matrices of R. This is because the evolutionary variance-covariance matrix are constrained to be positive semi-definite - a constraint which will cause problems for optim(). My solution to this problem was to optimize the Cholesky decomposite matrix, which does not have this constraint. My issue, then, is that if I return the Hessian matrix of partial second derivatives these will be the curvature of the likelihood surface for the Cholesky matrix elements - rather than for the variances and covariances of our evolutionary rate matrices.

My current plan to solve this involves first using optim() to find the MLEs based on the Cholesky decomposed matrices; and then feeding these parameter estimates along with a likelihood function for the full rate matrices into hessian() in the package "numDeriv". If this works I will soon report the result.

Adding standard errors to evol.vcv()

I haven't done too much programming in the past week or so due (in part) to a disrupted schedule from unusual geological and metereological events of recent days (as well as due to lizards & lots of new activity in my lab); however I am presently adding the calculation of standard errors for parameter estimates to the "phytools" function evol.vcv(). This function is based on Revell & Collar (2009) and was previously described here. In addition to this, on request from Dave Collar, I am planning to next add constraint matrices to test various hypotheses about the evolutionary variance-covariance matrix using this method (for instance, that the evolutionary rates - i.e., variances - are the same in different regimes, but the covariances differ). This will be coming soon.

To all those in Irene's path or wake - stay safe & dry!

"phytools" now on CRAN

The "phytools" library is now on CRAN ("phytools" CRAN page here). For most R users, that means it can be installed simply by typing:

> install.packages("phytools")

and then selecting a CRAN mirror.

Note that to install "phytools" this way, Windows users with earlier versions of R will have to upgrade to R 2.13.X; however you can still download the source or a Windows binary from the "phytools" webpage or CRAN page.

Installing "phytools"

I can see that at least a few visitors to the "phytools" blog have reached this page through the google search "installing phytools" - so I thought I'd say a couple of words about how this can be done at present.

1) Download the package from my R phylogenetics page. If you are a Windows user, you can download the Windows binary (this file ends in ".zip"). If you are a Mac or Linux/Unix user, download the package source (this file ends in ".tar.gz" - actually, you should also be able to install from source in Windows if you are so inclined).

2) Either navigate to the directory containing the "phytools" package installation file, or move the file to your current directory. For the former, you can either use the Windows GUI File->Change dir...; or the R command setwd(). For the latter, the R command getwd() (which returns your current working directory) should be helpful.

3) Now simply enter either:

install.packages("phytools_0.0-7.zip",repos=NULL)

(for Windows installation, assuming the latest version of "phytools"), or:

install.packages("phytools_0.0-7.tar.gz",repos=NULL,type="source")

for installation from source.

4) Finally, to load the "phytools" library, one simply types (as with any other contributed package):

library(phytools)

and you're ready to go!

**Note that I am getting ready to submit "phytools" to CRAN which will have the effect of collapsing steps 1-3 into a single command.**

New phylogenetic ANOVA function with posthoc tests

I just posted a new function to conduct the phylogenetic ANOVA of Garland et al. (1993). Luke Mahler has been doing phylogenetic ANOVA, but wanted to add posthoc comparison of the means between groups. I had previously programmed this in C - but I thought it would be easy enough to do in R as well. Direct link to the code is here. Note that the function, phylANOVA(), borrows a little bit from Luke Harmon's "geiger" function phy.anova(); and the function tTests() (which phylANOVA() calls internally) borrows code from the R "stats" function pairwise.t.test().

The phylogenetic ANOVA is pretty straightforward. First, we fit our ANOVA model (using lm()), and then we compute the F-statistic (using anova()). ANOVA, like most statistical tests, assumes independence - so the P-value returned by anova() will be incorrect. Garland et al.'s simple solution was to obtain the null distribution by simulation, which we can do by first calling X<-fastBM() a single time, and then computing the value of F for each column of X.

The posthoc tests are only a little bit more complicated. Here, we just compute a matrix of t-values for the real data and for each simulated data vector. Note that, for now, the method computes a single pooled standard-deviation which is used for all groups. For each pairwise comparison, we then count the number of simulated t-values with absolute values (for a two-sided test) that are larger than our corresponding empirical t. At the end, though, we have also conducted a large number of tests and that needs to be taken into consideration to control the experiment-wise error rate.

Fortunately in R this is straightforward to accomplish using the p.adjust() function in the "stats" package. p.adjust() implements a whole slew of methods for multiple test correction, and all of them can be called from phylANOVA(...,posthoc=TRUE,p.adj). For instance, standard Bonferroni correction is called by phylANOVA(...,p.adj="bonferroni") and sequential-Bonferroni (the default, also called the Holm-Bonferroni method) is called by phylANOVA(...,p.adj="holm").

Luke reports that the function seems to work properly.