For fun: least squares phylogeny estimation

For my class this semester on "Applied Phylogenetic Methods" (which, based on a chat today with Ron Etter's postdoc Rob Jennings, I've now concluded should more properly be called "Phylogenetics Demystified"), I recently programmed an implementation of phylogeny estimation using the unweighted least squares method (on my R-phylogenetics page, code here [update, v0.2]) - originally developed by Cavalli-Sforza & Edwards (1967), and described in detail by Felsenstein (2004, which is also our textbook for the class). Felsenstein (2004) calls the least squares methods the statistically "best-justified" distance matrix methods (p. 148). To my knowledge there is no direct implementation of least squares in any R package [although it should be possible to call the fastME program from within R and minimum evolution uses the least squares branch lengths].

As the title of this post suggests, this implementation is really just for fun, and should not be considered a serious function for phylogeny inference! This is because it is really not programmed efficiently and is thus very, very slow for more than a modest number of taxa. It also uses only nearest neighbor interchange to search treespace (capitalizing on the nni() function of the very useful {phangorn} package) thus will tend to get stuck on "islands" of suboptimal trees (not connected to better trees by a single NNI) if the starting tree is not very good. Nonetheless, the function does appear to work. If we give a matrix containing the true evolutionary distances between species (i.e., the distances in the true tree), then the program will tend to find the correct topology and branch lengths given a reasonable starting tree.

The basic principle of least squares phylogeny estimation is straightforward. We would like to find a tree and set of branch lengths that minimizes the summed squared deviations between our observed distances (i.e., the set of distances that serve as input for our analysis) and the hypothesized distances (i.e., the patristic distances in our inferred tree).

So, in other words, our model is as follows: D = Xv+e; in which D is a vector of distances between species; X is a design matrix based on a hypothesized topology (to be discussed below); v is a vector of edge lengths; and e, of course, is our residual error.

Given a topology to solve for our estimated edge lengths we can just use our standard least squares estimating equation: v = (X'X)^-1X'D. We can then calculate the sum of squares residual error (the quantity we want to minimize) for a given tree as: Q = (D-Xv)'(D-Xv).

So, there are two things we want to find. We want to find v that minimize Q conditioned a tree topology (given by X), which we can easily do using least squares estimation; but we also want to find the tree that minimizes Q given our least squares estimate of v.

We first need to compute X given a tree. I have not described the design matrix X to this point, but X is an N=n(n-1)/2 (for n species) by m=2n-3 (for a fully bifurcating unrooted tree) matrix. Here, you can probably see that N is just the number of distances between species; and m is the number of edges in the tree. We put a "1" in the X_ij if the jth edge can be found on the path between the pair of species represented in the ith row of the matrix. To make this a little clearer, consider the unrooted, fully labeled tree at right and the design matrix calculated with my function phyloDesign() below:

> phyloDesign(tree)
    6,7 7,1 7,2 6,8 8,3 8,4 6,5
1,2   0   1   1   0   0   0   0
1,3   1   1   0   1   1   0   0
1,4   1   1   0   1   0   1   0
1,5   1   1   0   0   0   0   1
2,3   1   0   1   1   1   0   0
2,4   1   0   1   1   0   1   0
2,5   1   0   1   0   0   0   1
3,4   0   0   0   0   1   1   0
3,5   0   0   0   1   1   0   1
4,5   0   0   0   1   0   1   1

Now, with this design matrix we can compute our least squares branches and Q, and it is then only a matter of finding the tree (and thus the design matrix, X) that minimizes Q. To do this, I have used nni() from the {phangorn} package to design a greedy (but nonetheless remarkably slow) algorithm for finding the least squares tree. nni() returns all the one-step away NNI trees for a given tree. So, I just returned all these, compute X, v, and Q for each tree, pick the lowest one, and repeat until Q no longer improves from NNIs. Obviously, this algorithm will easily get stuck on low-lying "hills" in treespace. I'm not familiar with other more aggressive tree rearrangement functions available in R (for instance, TBR), but perhaps my readers are and if so, they can tell me!

To see how this works, let's take the distance matrix for 8 mammals given by Felsenstein (2004; p. 163), as follows:

> D.mammals
         dog bear raccoon weasel seal sea_lion cat monkey
dog        0   32      48     51   50       48  98    148
bear      32    0      26     34   29       33  84    136
raccoon   48   26       0     42   44       44  92    152
weasel    51   34      42      0   44       38  86    142
seal      50   29      44     44    0       24  89    142
sea_lion  48   33      44     38   24        0  90    142
cat       98   84      92     86   89       90   0    148
monkey   148  136     152    142  142      142 148      0

Now, let's load the function from source:

> source("optim.phylo.ls.R")

This will actually load four functions: optim.phylo.ls(), ls.tree(), phyloDesign(), and compute.ancestor.nodes(). It will also load the packages {phangorn} (and, of course, {ape}) on first execution.

> lstree<-optim.phylo.ls(D.mammals,upgma(D.mammals))
[1] "1 set(s) of nearest neighbor interchanges. best Q so far = 143.854166666667"
[1] "2 set(s) of nearest neighbor interchanges. best Q so far = 98.8333333333332"
[1] "3 set(s) of nearest neighbor interchanges. best Q so far = 92.1166666666666"
best Q score of 92.1166666666666 found after 3 nearest neighber interchange(s).

Using the UPGMA tree as our starting tree. Then we can plot our tree and edge lengths:

> plot(midpoint(lstree)) # use midpoint rooting
> edgelabels(round(midpoint(lstree)$edge.length,2), adj=c(0.5,-0.2),frame="none")



which results in a tree that makes a fair bit of sense. Note that the zero length branch in this tree was actually negative in the LS tree, but my program sets negative branch lengths to zero by default (although this can be changed by the user).