tag:blogger.com,1999:blog-8499895524521663926.post1845959927044206893..comments2023-11-27T05:04:22.702-05:00Comments on Phylogenetic Tools for Comparative Biology: For fun: least squares phylogeny estimationLiam Revellhttp://www.blogger.com/profile/04314686830842384151noreply@blogger.comBlogger6125tag:blogger.com,1999:blog-8499895524521663926.post-41542336275850160782011-04-02T02:25:22.119-04:002011-04-02T02:25:22.119-04:00If I understand correctly, for a given tree topolo...If I understand correctly, for a given tree topology the branch lengths are estimated by least squares for both LS and ME methods, but then the choice among topologies is made by the LS criterion in the first case, and by minimizing the sum of the branch lengths in the second. That might choose the same topology but is not guaranteed to do so. If it does come to the same topology as LS of course the branch lengths are the same. Note that in the ME case you do want to constrain branch lengths not to be negative.Joe Felsensteinhttps://www.blogger.com/profile/06359126552631140000noreply@blogger.comtag:blogger.com,1999:blog-8499895524521663926.post-31804399393328327872011-03-30T12:10:23.739-04:002011-03-30T12:10:23.739-04:00Hi Liam,
I believe at least in this case they ar...Hi Liam, <br /><br />I believe at least in this case they are the same. I tried it a few times and when I finally managed to order the tree, the design and distance matrix it seemed to be the case. I found there is always a bit confusion around minimum evolution: minimum evolution as an optimality criterion and as an algorithm (heuristic) like fastME. <br />Also fastME mostly comes as the weighted version - also very easy to implement - the weights are 2^(-rowSums(X)), I believe. <br />Anyway, with the code above, it should be now easy to verify numerically, whether edge lengths of fastme.ols and LS are indeed identical. So LS is not only for fun, but maybe also useful ;)<br /><br />KlausKlaushttps://www.blogger.com/profile/11021989593338482289noreply@blogger.comtag:blogger.com,1999:blog-8499895524521663926.post-1788400954338018562011-03-30T00:12:10.175-04:002011-03-30T00:12:10.175-04:00Hi Klaus.
Thanks for the many useful suggestions....Hi Klaus.<br /><br />Thanks for the many useful suggestions. At some point I will try to add these to the function and see what happens.<br /><br />I believe, if I'm not wrong, that fastme.ols() optimizes under the minimum evolution criterion using the least squares branch lengths; which is different than optimizing under the least squares criterion using same. See Felsenstein (2004; pp. 148-161) for a description of this distinction.<br /><br />- LiamLiam Revellhttps://www.blogger.com/profile/04314686830842384151noreply@blogger.comtag:blogger.com,1999:blog-8499895524521663926.post-32938666695294867922011-03-29T15:37:03.297-04:002011-03-29T15:37:03.297-04:00Hi Liam,
it seems you may have missed that there...Hi Liam, <br /><br />it seems you may have missed that there are fastme.ols and fastme.bal functions are already in ape. These functions can do even SPR moves, but the names are not intuitive.<br /><br />To speed up your LS code a few tricks. If you are lazy and do not want to implement the code of Joe or the code from Bryant and Waddell, you can just make use of the fact that X is highly sparse: <br />library(Matrix)<br />and<br />X <- Matrix(X) <br />somewhere in the code will already speed up your code and saves memory, especially for bigger matrices. The Matrix package also contains a vignette (in R type vignette("Comparisons") ), which tells you everything one wants to know about how to compute LS fast, e.g.: <br />v<-solve(crossprod(X), crossprod(X,colD))<br />is much faster than<br />v<-solve(t(X)%*%X)%*%t(X)%*%colD <br />a trick everybody should know. <br />There is also a function designTree in phangorn which does pretty much the same as phyloDesign and maybe is a bit faster. I will change it to the Matrix class in one of the next releases. <br />Speeding up the tree search on the other hand is tricky. The most expensive operation is likely the inversion of the matrix (t(X)%*%X). One could possibly use a partial inversion of this matrix as NNI moves only affect some blocks of this matrix.<br /> <br />You may also be interested in this small function to compute non-negative least-squares I wrote recently, following a discussion on R-phylo-sig:<br />#<br /># tree is a tree of class "phylo"<br /># dm a distance matrix<br />#<br />require(phangorn)<br />nnls.tree <- function(tree, dm){<br /> dm = as.matrix(dm)<br /> k = dim(dm)[1]<br /> labels = tree$tip<br /> dm = dm[labels,labels]<br /> y = dm[lower.tri(dm)]<br /> X = designTree(tree, method = "unrooted")<br /> betahat = lm.fit(X,y)$coefficients<br /> if(!any(betahat<0)){<br /> tree$edge.length[] = betahat<br /> return(tree)<br /> }<br /> n = dim(X)[2]<br /> Dmat <- crossprod(X)<br /> dvec <- crossprod(X, y)<br /> Amat <- diag(n)<br /> betahat <- solve.QP(Dmat,dvec,Amat)$sol<br /> tree$edge.length[] = betahat<br /> tree<br />}<br />I will also include it in the next phangorn version. <br /><br />Regards,<br />KlausKlaushttps://www.blogger.com/profile/11021989593338482289noreply@blogger.comtag:blogger.com,1999:blog-8499895524521663926.post-13076118413965225682011-03-26T14:37:49.718-04:002011-03-26T14:37:49.718-04:00Hi Joe.
I have read your paper but not the other ...Hi Joe.<br /><br />I have read <a href="http://dx.doi.org/10.1093/sysbio/46.1.101" rel="nofollow">your paper</a> but not the other one.<br /><br />In my class I programmed the "hard" parts ahead of time, and then we learned how to create a custom function in R to do inference - so my goal was indeed simplicity.<br /><br />Unfortunately, our function didn't work so I had to spend an hour or two after class debugging!<br /><br />Even though my function is slow, if you give it the NJ tree as a starting tree it will produce a sensible result in reasonable order. (For instance, in the mammal tree the unweighted LS tree is only one NNI away from the NJ phylogeny.) Adding weights (e.g., 1/D[i,j] as in Fitch & Margoliash <a href="http://dx.doi.org/10.1126/science.155.3760.279" rel="nofollow">1967</a>) or changing to GLS would also be pretty straightforward.<br /><br />Thanks for the comment!Liam Revellhttps://www.blogger.com/profile/04314686830842384151noreply@blogger.comtag:blogger.com,1999:blog-8499895524521663926.post-22836771147033991572011-03-26T09:58:05.206-04:002011-03-26T09:58:05.206-04:00There are faster algorithms for least squares ohyl...There are faster algorithms for least squares ohylogeny in two papers:<br /><br />(1) My own iterative LS algorithm in my 1997 paper in <i>Systematic Biology</i>. It is fairly fast and has the convenient property of allowing us easily to keep branch lengths from going negative, if that's what we want. It involves a whole scheme for "pruning" distances on a tree, and can reuse parts of the calculation if a only small change of the tree is made (however in principle we must re-evaluate every branch length when even one branch's length is changed).<br /><br />(2) Some algorithms by Bryant and Waddell in <i>Molecular Biology and Evolution</i> in 1998. These are even faster.<br /><br />Both of these are much more complex to implement, so for a teaching example I sympathize with your choices here.Joe Felsensteinhttps://www.blogger.com/profile/06359126552631140000noreply@blogger.com