Three more functions & some more methods in Rphylip

I just added a few more functions to the Rphylip project, my R interface for the PHYLIP package. The new interface functions are Rpars (for PARS), Rmix (for MIX), and Rpenny (for PENNY). All three of these are parsimony method programs: the first does heuristic MP search from unordered, multistate data; whereas the latter do (Wagner, Camin-Sokal, or mixed method) MP searching using heuristic or branch-and-bound algorithms, respectively. More details on the programs can be found by referring to the PHYLIP documentation pages linked above.

I also created a new class of data object, "phylip.data", which just generalizes "proseq" (in Rphylip) and "DNAbin" (in ape), and is very simple.

Here's a quick demo using Rpenny. Note that branch-and-bound should generally not be used for more than a dozen or so taxa (it will become computationally prohibitive quickly).

> require(Rphylip)
Loading required package: Rphylip
Loading required package: ape
> packageVersion("Rphylip")
[1] ‘0.1.14’
> data(primates.bin)
> primates.bin
12 character value sequences stored in a matrix.

All sequences of same length: 231

Labels: Lemur Tarsier Sq.Monkey J.Macaque R.Macaque E.Macaque ...

Trait value composition:
    0     1
0.406 0.594
> tree<-Rpenny(primates.bin)

....

How many
trees looked                                 Approximate
at so far    Length of      How many         percentage
(multiples   shortest tree  trees this long  searched
of  100):    found so far   found so far     so far
----------   ------------   ------------     ------------
     1       -                      0            0.00
     2       -                      0            0.00
     3       -                      0            0.00
     4       208.00000              3            0.00
     5       208.00000              6            0.00
     6       208.00000              6            0.14
     7       208.00000              6            1.90
     8       208.00000              6            6.67
     9       208.00000              6            9.33
    10       208.00000              6           14.00
    11       208.00000              6           37.78
    12       208.00000              6           53.33

Output written to file "outfile"

Trees also written onto file "outtree"

Press enter to quit.

Penny algorithm, version 3.695
 branch-and-bound to find all most parsimonious trees

Wagner parsimony method


                                    requires a total of            208.000

    6 trees in all found


  +--------------------------------1
  ! 
  !  +-----------------------------2        
  !  ! 
--1  !                       +-----10       
  !  !                    +-10 
  !  !                    !  !  +--11       
  !  !                 +--6  +--8 
  !  !                 !  !     +--12       
  +--2     +-----------5  ! 
     !     !           !  +--------9        
     !     !           ! 
     !     !           +-----------8        
     !  +--4 
     !  !  !                 +-----6        
     !  !  !              +-11 
     !  !  !              !  !  +--5        
     +--3  +--------------7  +--9 
        !                 !     +--7        
        !                 ! 
        !                 +--------4        
        ! 
        +--------------------------3        

  remember: this is an unrooted tree!

....

Translation table
-----------------
        1       Lemur
        2       Tarsier
        3       Sq.Monkey
        4       J.Macaque
        5       R.Macaque
        6       E.Macaque
        7       B.Macaque
        8       Gibbon
        9       Orangutan
        10      Gorilla
        11      Chimp
        12      Human

Rooted tree(s) with the outgroup
------------------------
Tarsier, Lemur

> require(phytools)
Loading required package: phytools
Loading required package: maps
Loading required package: rgl
> par(mfrow=c(3,2))
> plotTree(tree)
Waiting to confirm page change...

(These are the six equally most parsimonious trees found by PENNY.)

Cool. The latest version of Rphylip can be downloaded here, and is also on GitHub.

More updates to rateshift method: Testing for the presence of a rate shift

I have made some updates to the function rateshift (first described here) to facilitate comparison of alternative models for rate shifts; as well as for testing the null hypothesis of no shift.

First, I fixed the function so it could fit a no-rate-shift model. That was broken in the previous version, but should work now. When we fit the nrates=1 model, the fitted model parameter value (σ²) and log-likelihood should be the same as from (say) fitContinuous in geiger or the one-rate model in brownie.lite.

Second, I created an S3 generic logLik method for the object of class "rateshift" returned by the function. This allows us to easily extract the log-likelihood & model parameterization; but it also allows us to use the generic AIC to compute the Akaike Information Criterion value for the fitted model.

Finally, third, I fixed a minor bug which sometimes created an incompatibility in the tolerance (basically, the very small values we need to add or subtract from some quantities to make sure that the function does not attempt to evaluate the likelihood where it isn't defined) are inconsistent between rateshift and make.era.map, which is used internally. This required changes to both functions, so the wisest thing to do to get this update is to update phytools to the latest non-CRAN version.

OK. Here's a demo:

> require(phytools)
Loading required package: phytools
Loading required package: ape
Loading required package: maps
Loading required package: rgl
> packageVersion("phytools")
[1] ‘0.3.86’
> ## simulate tree & data
> tree<-pbtree(n=100,scale=1)
> tree<-make.era.map(tree,c(0,0.5,0.8))
> x<-sim.rates(tree,c(1,10,1),internal=TRUE)
names absent from sig2: assuming same order as $mapped.edge
> ## here's a visual of our simulation
> phenogram(tree,x,ftype="off")

> ## peel off ancestral states
> x<-x[tree$tip.label]
>
> ## fit 1 rate model
> fit1<-rateshift(tree,x,nrates=1)
> fit1
ML 1-rate model:
      s^2(1)  se(1)  k  logL   
value 1.7966  0.2542 2  -81.8257

This is a one-rate model.

R thinks it has found the ML solution.

> ## fit 2 rate model
> fit2<-rateshift(tree,x,nrates=2)
> fit2
ML 2-rate model:
      s^2(1)  se(1)  s^2(2)  se(2)  k  logL   
value 6.6364  2.1013  0.785  0.1345 4  -64.3827

Shift point(s) between regimes (height above root):
        1|2    se(1|2)
value  0.806  0.02

R thinks it has found the ML solution.

> ## test 2 rates vs 1 rate
> P2vs1<-as.numeric(pchisq(2*(logLik(fit2)-logLik(fit1)),   df=attr(logLik(fit2),"df")-attr(logLik(fit1),"df"),   lower.tail=FALSE)) > P2vs1
[1] 2.658128e-08
> ## fit 3 rate model
> fit3<-rateshift(tree,x,nrates=3)
> fit3
ML 3-rate model:
      s^2(1)  se(1) s^2(2) se(2)  s^2(3)  se(3)  k  logL
value 1.7181  NaN  6.2815 1.8016  0.7716  0.1345 6  -64.017

Shift point(s) between regimes (height above root):
        1|2    se(1|2) 2|3    se(2|3)
value  0.3058  0.0265  0.8193  0.0141

R thinks it has found the ML solution.

> ## test 3 rates vs 2 rates
> P3vs2<-as.numeric(pchisq(2*(logLik(fit3)-logLik(fit2)),
  df=attr(logLik(fit3),"df")-attr(logLik(fit2),"df"),
  lower.tail=FALSE))
> P3vs2
[1] 0.6934852

This shows us that although the fitted shift points in our third fitted model are fairly close to the gnerating shift points, the fit isn't significantly better than our two rate model. I suspect that, in general, it will probably be easier to find shift points that are closer to the tips of the tree, where there tends to be more edges.

Cool.

New method to locate one or multiple rate shifts on a tree using likelihood

I just posted a new phytools function, rateshift, that fits a model in which there are one or multiple Brownian rate shifts in the tree at different heights above the root. The idea is that we don't need to specify the locations of the rate shifts a priori (as we can already do using brownie.lite); rather, we let the data determine where the rate shifts are located. Turns out that this isn't too hard, it just requires p-1 additional parameters for each extra rate above the root.

It also occurred to me that it's entirely possible this method is already in the literature - in which case, please accept my apologies for having missed it!

In the simplest case, this would just be a model with two evolutionary rates: σ²(1) rootward of the shift, and σ²(2) tipward; and one rate shift. We then jointly maximize the likelihood of the rates & position of the rate shift.

Here's a quick demo:

> require(phytools)
Loading required package: phytools
Loading required package: ape
Loading required package: maps
Loading required package: rgl
> packageVersion("phytools")
[1] ‘0.3.85’
> tree<-pbtree(n=100,scale=1)
> tree<-make.era.map(tree,c(0,0.5,0.8))
> x<-sim.rates(tree,c(1,10,1))
names absent from sig2: assuming same order as $mapped.edge
> fit<-rateshift(tree,x,nrates=3,print=TRUE,plot=TRUE, tol=1e-5)
Optimization progress:

s^2(1)  s^2(2)  s^2(3)  shift:1 shift:2 logL
2.623   2.623   2.623   0.3333  0.6667  -117.6615
2.624   2.623   2.623   0.3333  0.6667  -117.6616
....

2.1182  7.6728  0.7129  0.5497  0.8235  -100.7116
2.1182  7.6728  0.7109  0.5497  0.8235  -100.7163
2.1182  7.6728  0.7119  0.5507  0.8235  -100.7079
....
1.1796  8.6499  0.7848  0.5886  0.8116  -100.1251
1.1796  8.6499  0.7848  0.5896  0.8126  -100.1175
1.1796  8.6499  0.7848  0.5896  0.8106  -100.1275

> fit
ML 3-rate model:
      s^2(1) se(1)  s^2(2) se(2)  s^2(3) se(3)  k  logL
value 1.1796 0.991  8.6499 2.3647 0.7848 0.1709 6  -100.12

Shift point(s) between regimes (height above root):
        1|2     se(1|2) 2|3     se(2|3)
value   0.5896  0.0173  0.8126  0.0173

R thinks it has found the ML solution.

The options plot & print slow down runtime; however I included them for the purposes of debugging and it is kind of neat to visualize the optimization of the locations of the rate shift.

The implementation is a little buggy - however at least in this example it seems to do pretty well at finding our generating rate shift points (which were, remember, 0.5 & 0.8 units above the root), and rates (1, 10, & 1). Cool!

This function is in a new phytools build (phytools 0.3-85), which can be downloaded & installed from source. Please check it out.

Three years of blogging

It's now been three years since I starting blogging about phylogeny methods on blog.phytools.org (originally phytools.blogspot.com) - and so, in the tradition of 2011 and 2012, I thought I'd spend a few minutes talking about what I did this year in the phytools package and on the blog.

According to the (somewhat dubious, in my opinion) blogger.com page stats, the phytools blog received upwards of 150,000 page views in 2013. Even if 1/2 of these were by bots, that is still quite an impressive tally. Certainly, by any measure the popularity of the phytools blog as a free repository of information about phytools and phylogeny methods has increased over the past year.

Towards the end of 2012 and throughout 2013 I added considerably to the plotting capabilities of phytools (evidenced, in part, by my recent MEE paper on some new plotting methods). Consequently, it's no great surprise that two of the most viewed phytools blog posts of 2013 included the description of a new method to visualize uncertainty on a traitgram, and some of the description and illustration of the new plotting functions contMap and densityMap, including this description of a published use of both methods in the same figure (which probably didn't hurt by having been re-tweeted by @systbiol). Nestled also among the top three most popular blog posts of 2013 was this comment on why we don't normally expect the residuals from phylogenetic ANOVA or regression to be normally distributed. Also popular were a wide range of posts about stochastic mapping and ancestral state reconstruction, including information about a new method based on the the threshold model from evolutionary quantitative genetics.

Towards the end of 2013, I started on a new project Rhylip. The purpose of Rphylip is to create an R interface for all 30+ programs in the PHYLIP phylogeny method software package by Joe Felsenstein. This will hopefully allow the many functions of PHYLIP to be used seamlessly within an integrated R workflow. (Here's an example of Rphylip at work - created for my phylogeny methods class using knitr.) I'm about 50% of the way there, so I hope to get this done soon.

Finally, I recently learned that my CAREER proposal to do phylogeny method research in several new areas has been recommended for funding by the NSF DEB Systematic Biology program. (Prospective students & postdocs take notice - I will be hiring in 2014!) This comes exactly at the right time as my start-up is in its dying breaths right now. Interestingly, part of the original ulterior motive in developing this blog back towards the end of 2010 was as a supporting 'broader impact' for what was at that time my first attempt to acquire funding for phylogeny method research from the NSF. Thus, the NSF is in part responsible for this blog & phytools as a community resource even without having funded it! (Until now, that is.)

Happy 2014 & thanks for reading!

S3 print method for brownie.lite

Today I wrote an S3 print method for objects returned by the phytools function brownie.lite. This function implements the method of O'Meara et al. (2006) in which we fit different rates of evolution to different parts of a phylogeny (sometimes, say, determined by a mapped discrete character). Print methods are nice because they allow us to tell R how to summarize the content of a complicated object, like a phylogeny or the results of numerical optimization.

The object returned by brownie.lite is unchanged, except that it is now assigned the class attribute "brownie.lite". Here's how it works:

> require(phytools)
> packageVersion("phytools")
[1] ‘0.3.83’
> ## first let's simulate some data
> Q<-matrix(c(-1,1,1,-1),2,2)
> rownames(Q)<-colnames(Q)<-letters[1:2]
> tree<-sim.history(pbtree(n=100,scale=1),Q,anc="a")
> ## simulate with a rate difference among states
> x<-sim.rates(tree,setNames(c(1,2),letters[1:2]))
> ## simulate without a rate difference
> y<-fastBM(tree)
> plotSimmap(tree,setNames(c("blue","red"),letters[1:2]), lwd=3,ftype="off")

> fitx<-brownie.lite(tree,x)
> fitx # this is the same as print(fitx)
ML single-rate model:
          s^2     se      a       k  logL
parameter 1.5484  0.2191  -0.4682 2  -73.73

ML multi-rate model:
          s^2(a)  se(a)   s^2(b)  se(b)   a       k  logL
parameter 0.732   0.2243  1.8163  0.2955  -0.499  3  -70.93

P-value (based on X^2): 0.0179

R thinks it has found the ML solution.

> fity<-brownie.lite(tree,y)
> fity
ML single-rate model:
          s^2     se      a       k  logL
parameter 1.0646  0.1507  -0.7218 2  -55.00
      
ML multi-rate model:
          s^2(a)  se(a)   s^2(b)  se(b)   a       k  logL
parameter 0.8925  0.2661  1.1209  0.1841  -0.7271 3  -54.79

P-value (based on X^2): 0.5184

R thinks it has found the ML solution.

This version of phytools can be downloaded here and installed from source.

More functions & fixes in Rphylip

My latest updates to Rphylip, an R interface for Felsenstein's phylogeny methods package PHYLIP, includes Rdnainvar, an interface for dnainvar, a program for Lake's invariants from DNA sequences; and Rfitch, an interface for fitch, a program that does least-squares and minimum evolution phylogeny inference.

In addition to this, though, today I added the functions setPath & clearPath. The former can be used to set the path to the PHYLIP executables folder for the current R session. Normally, Rphylip tries to find the path to the PHYLIP executable using the internal Rphylip function, findPath. However this function only works if the PHYLIP folder is in a sensible location, such as in C:/Program Files/ (for Windows) or /Applications/ (for Mac OS). If not, the path can always be supplied to the function directly. What setPath does is enable the user to set a different place to look for PHYLIP executables for the current R session. It does this by setting the environmental variable "phylip.path". clearPath does the opposite. It removes the environmental variable "phylip.path" thus permitting Rphylip functions to search for PHYLIP executables in the regular places.

Here's a demo in which I have PHYLIP installed in C:/Users/Liam/Documents/Phylogeny_Programs/:

> require(Rphylip)
Loading required package: Rphylip
Loading required package: ape
> packageVersion("Rphylip")
[1] ‘0.1.13’
> data(primates)
> D<-Rdnadist(primates)
Error in Rdnadist(primates) :
  No path provided and was not able to find path to dnadist
> setPath("C:/Users/Liam/Documents/Phylogeny_Programs/phylip-3.695/exe/")
> D<-Rdnadist(primates)

...

Nucleic acid sequence Distance Matrix program, version 3.695

...

> tree<-Rfitch(D)

...

Fitch-Margoliash method version 3.695

...

  12 Populations

Fitch-Margoliash method version 3.695

                  __ __             2
                  \  \   (Obs - Exp)
Sum of squares =  /_ /_  ------------
                                2
                   i  j      Obs

Negative branch lengths not allowed

global optimization


                                         +-----------7
                                +--------2
                                !        !   +-------6
                                !        +---1
                                !            !  +-4
                       +--------3            +--8
                       !        !               +----5
                       !        !
                       !        !      +-----------------8
                       !        +------9
                       !               !   +-------------9
  +--------------------7               +---5
  !                    !                   !    +-------10
  !                    !                   +----6
  !                    !                        ! +------12
  !                    !                        +-4
  !                    !                          +-------11
  !                    !
  !                    +----------------------------3        
  ! 
 10--------------------------------------2
  ! 
  +-------------------------------------1


remember: this is an unrooted tree!

Sum of squares =     0.68048

Average percent standard deviation =     7.23497

...

Translation table
-----------------
        1       Lemur
        2       Tarsier
        3       Sq.Monkey
        4       J.Macaque
        5       R.Macaque
        6       E.Macaque
        7       B.Macaque
        8       Gibbon
        9       Orangutan
        10      Gorilla
        11      Chimp
        12      Human

> plot(tree)

New di2multi function for stochastically mapped trees

I just posted some code the collapses internal branches of zero length (or, more specifically, branches with length shorter than some arbitrarily specified value tol) to polytomies for trees with mapped discrete characters created using (for instance) make.simmap or read.simmap. This is exactly the same as di2multi in the ape package (functionally - I programmed it totally differently, hopefully not at the peril of users); however it works for modified "phylo" objects with a mapped discrete trait.

The code to the function is here, but it uses internal phytools functions - so for it to work, you will probably have to install the latest version of phytools (here).

Here's a quick demo of the function at work:

> require(phytools)
> packageVersion("phytools")
[1] ‘0.3.81’
> ## first create a tree with some polytomies
> N<-26
> tree<-pbtree(n=N,tip.label=LETTERS[26:1],scale=1)
> tree$edge.length[tree$edge.length<0.1&tree$edge[,2]>N]<-0
> ## make the tree ultrametric
> d<-max(vcv(tree))-diag(vcv(tree))
> tree$edge.length[tree$edge[,2]<=N]<- tree$edge.length[tree$edge[,2]<=N]+d
> ## check is ultrametric
> is.ultrametric(tree)
[1] TRUE
> ## check is binary
> is.binary.tree(tree
[1] TRUE
> ## simulate discrete character history on tree
> Q<-2*matrix(c(-1,1,1,-1),2,2)
> tree<-sim.history(tree,Q)
> plotSimmap(tree,colors=setNames(c("blue","red"),1:2), lwd=3)

> ## collapse branches of zero length
> tree<-di2multi.simmap(tree)
> ## check is binary
> is.binary.tree(tree)
[1] FALSE
> plotSimmap(tree,colors=setNames(c("blue","red"),1:2), lwd=3)

The main reason this might be important is because densityMap runs into problems when your tree has internal edges that are very short. This can be addressed by first collapsing all zero/small branches in each of the stochastic mapped trees, and then computing the density-map on the collapses stochastic mapped trees. Code for that would look like the following:

trees<-lapply(trees,di2multi.simmap)
class(trees)<-"multiPhylo"
densityMap(trees)