New version of ancThresh for λ model

The current version of ancThresh (for ancestral character estimation under the threshold model; see Revell In press) permits a Brownian or Ornstein-Uhlenbeck model for the evolution of the liabilities. I just added a 3rd model, the λ model of Pagel (1999). The code for this version is here; but since it uses some functions internally that have also been updated, the best thing to do is to update phytools to the latest version.

I'm not a huge fan of the λ model in general, since it is not clear what biological process it is meant to approximate; however under some circumstances it could be a useful model for traits that evolve on the tree (and thus have phylogenetic covariance), but are also affected by contemporary factors that are not necessarily phylogenetically correlated. This is how I am using this model in the empirical study for which I have added this feature to ancThresh.

At the moment it is largely untested - so if you run into any problems, please let me know.

Bug fix in phylANOVA

A phytools user recently reported discrepant results between phylANOVA in phytools and aov.phylo in geiger. Both functions conduct the simulation-based method of Garland et al. (2013). In theory, the only difference is that phylANOVA performs post-hoc comparison of means. It turns out, however, that phylANOVA contains the implicit assumption that y is in the order of tree$tip.label. This assumption is now only true if names(y) is NULL, in which case a warning is also issued. Updated code is here and in the latest version of phytools.

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.