phytools 0.2-40 on CRAN

A new version of phytools is now available on CRAN. I submitted the last CRAN version of phytools March 20th, so there are not a huge number of updates in the present version - but some of them are important or very cool, so I decided to get the new phytools version on CRAN anyway. Over the new few days, Windows & Mac OS X binaries should be built and then gradually percolate through the CRAN mirror repositories.

Here are some of the updates in this version relative to phytools 0.2-30, the last CRAN phytools version:

1. A new function to get the marginal ancestral state reconstructions of a discrete character using the re-rooting method.

2. A new function to summarize the results of stochastic mapping (1, 2, 3).

3. An update to make.simmap to allow uncertain tip states.

4. A small update to phylosig in the calculation of P-values.

5. An update to make.simmap to allow it to performing mapping on a set of input trees, for instance from the posterior distribution of a Bayesian analysis.

6. An important bug fix in make.simmap for asymmetric substitution models (e.g., model="ARD").

7. Finally, a new, totally rewritten version of phylomorphospace that also allows users to show a mapped discrete character on the tree.

Please report any bugs or issues.

New completely re-written version of phylomorphospace with mapped discrete trait

I just posted a new, completely re-written version of phylmorphospace - the phytools function that does a projection of the tree into a two dimensional morphospace. Since I wrote the original version way back in 2010, and the guts of the function had persisted largely unchanged since that time, this was probably overdue.

There was no major issue with the prior version; however I wanted to add some features and realized that the code could be a lot nicer - so I decided to re-write the function (pretty much) from scratch.

Phylomorphospace plots are pretty easy in principle. We just need to supply tip states & compute or supply all the states at internal nodes. Having done that, we can just plot all the tip & node states in our bivariate morphospace, and then add edges connecting all parent & daughter nodes.

Source code for the updated version of phylomorphospace is here; I have also posted a new build of phytools with these updates (phytools 0.2-38).

The update that inspired the re-write was that I wanted to be able to plot the state of a mapped discrete character along the edges of the tree, à la (for example) this version of phenogram. To do this for a projection of the tree into two dimensions is a little more complicated, because in phenogram the time spent in a mapped state is just plotted on the interval demarcated by the horizontal (i.e., time) axis. In two phenotypic trait dimensions, this is a little more complicated. Here, we have to compute the proportion of time spent in each state on each edge and then color the edge proportionally by those states, accordingly.

OK, here's a quick demo:

> require(phytools)
Loading required package: phytools
> packageVersion("phytools")
[1] ‘0.2.38’
> # first the standard version
> tree<-pbtree(n=20,scale=1)
> X<-fastBM(tree,nsim=2)
> phylomorphospace(tree,X,xlab="X1",ylab="X2")

OK, now for something more interesting let's simulate a discrete character on the tree; and then generate data for two continuous traits in which both the rate & evolutionary correlation differ depending on the mapped discrete character:

> Q<-0.5*matrix(c(-1,1,1,-1),2,2)
> rownames(Q)<-colnames(Q)<-letters[1:2]
> mtree<-sim.history(tree,Q,anc="a")
> # here's our discrete character history on the tree
> plotSimmap(mtree,colors=cols,pts=F)>

> # this is for simulation
> R<-list(matrix(c(0.5,0.45,0.45,0.5),2,2),
+ matrix(c(2,0,0,2),2,2))
> names(R)<-letters[1:2]
> X<-sim.corrs(mtree,R)
> cols<-c("blue","red"); names(cols)<-letters[1:2]
> phylomorphospace(mtree,X,xlab="X1",ylab="X1",colors=cols)

Pretty cool, I guess.... The evolutionary pattern that we simulated - low rate & high evolutionary correlation on the blue branches; & high rate but low evolutionary correlation on the red branches - is pretty evident in the plot.

One little note about plotting tip labels. During the re-write I noticed that I'd used the function textxy from the package 'calibrate' in place of the base graphics function text - but I'd forgotten why. Turns out textxy is a neat function that plots text labels for points with an offset that varies depending on the plot quadrant. This is perfect for a function like phylomorphospace, because it helps push the labels away from other plotted lines & points.

That's it for now.

Bug fix for make.simmap with asymmetric substitution model; new version of phytools

Yesterday I received a user report of some problems with make.simmap(...,model="ARD") when it resulted in some of the fitted transition rates being zero. This was a known (to me) issue with make.simmap, and it is because although we can compute the conditional likelihoods with this matrix no problem - when we are trying to draw waiting times from an exponential distribution to map character changes along internal branches, rexp(...,rate=0) won't evaluate. One solution to this would be to return Inf or some arbitrarily large number when the rate is 0. Instead, and for other reasons of computation, I decide to add a small number, tol=1e-08 to off-diagonal position of Q that are 0 in the MLE. (There are also other calculations that make this necessary.)

Fixing this issue turned up another more serious problem and that is that recent versions of make.simmap have been using the transpose of Q in simulating along edges for asymmetric transition models, instead of Q itself (or, alternatively, that it has been calling a row index instead of a column index during an important stage in calculation). I believe that this bug appeared during my recent major re-write of make.simmap. Obviously, this doesn't affect symmetrical models of character change in which Q==t(Q) (such as model="ER" or model="SYM" - the default), but it will affect model="ARD".

Here's a little more specific detail on the error. In the internally called function smap, I had:

Q<-t(Q)

where I should not have; or, alternatively:

p<-expm(Q*tree$edge.length[j])[NN[j,1],]* L[as.character(tree$edge[j,2]),]

instead of:

p<-expm(Q*tree$edge.length[j])[,NN[j,1]]* L[as.character(tree$edge[j,2]),]

Source code for the fixed version of make.simmap is here. In this version, users can also control the value of tol by way of the optional argument, well, tol. tol is only used if any of the off-diagonal elements of Q are less than tol, which has a default value of tol=1e-08, as noted above.

This update to make.simmap is also in a new phytools package build, phytools 0.2-37, which can be installed from source.

Finally, here is a demo in which I simulate with a very low backward rate & then show what the new version of make.simmap does (instead of failing) if that backward transition rate has a MLE of 0. Note that if make.simmap seems to hang - it may be possible to resolve this by increasing tol.

> require(phytools)
Loading required package: phytools
> packageVersion("phytools")
[1] ‘0.2.37’
> tree<-pbtree(n=200,scale=1)
> Q<-matrix(c(-1,1,0.01,-0.01),2,2)
> rownames(Q)<-colnames(Q)<-letters[1:2]
> tree<-sim.history(tree,Q,anc="a")
> cols<-setNames(c("black","red"),letters[1:2])
> # this is the true character history
> plotSimmap(tree,cols,pts=F,ftype="off")

> mtrees<-make.simmap(tree,tree$states,model="ARD", nsim=100)

Warning: some elements of Q not numerically distinct from 0; setting to 1e-08

make.simmap is sampling character histories conditioned on the transition matrix
Q =
            a           b
a -0.96697341  0.96697341
b  0.00000001 -0.00000001
(estimated using likelihood);
and root node prior probabilities
pi =
  a   b
0.5 0.5

Done.

And a little reality check:

> # true history
> describe.simmap(tree)
1 tree with a mapped discrete character with states:
 a, b

tree has 26 changes between states

changes are of the following types:
  a  b
a 0 26
b 0  0

mean total time spent in each state is:
             a         b    total
raw  26.703708 18.438526 45.14223
prop  0.591546  0.408454  1.00000

> # stochastic maps
> describe.simmap(mtrees,plot=T,show.tip.label=FALSE)
100 trees with a mapped discrete character with states:
 a, b

trees have 25.48 changes between states on average

changes are of the following types:
       a,b b,a
x->y 25.48   0

mean total time spent in each state is:
              a          b    total
raw  26.2457674 18.8964671 45.14223
prop  0.5814016  0.4185984  1.00000

In the new phytools build I've also added the function getStates (which can be used to pull the states at nodes or tips from a tree with a mapped discrete character and is called internally by describe.simmap) to the NAMESPACE so that it can be called by phytools users.

Please don't hesitate to report any bugs or issues with the present version of make.simmap or phytools.

Thanks!

Using make.simmap on a set of trees

A recent commenter asked:

"I wonder if it would be possible to apply make.simmap to an object multiphylo (to deal with phylogenetic uncertainty) and to summarize the outcome on a consensus tree."

Let's take this one bit at a time. First, the task of applying make.simmap, the phytools function for stochastic character mapping, to a set of trees - say a sample from the posterior distribution in a Bayesian analysis.

At present, make.simmap takes a single tree and data vector as input; and can return as many simulated stochastic maps as the user demands. It is possible to iterate over a list of trees and then combine the results into a single object of class "multiPhylo" - but this is a little annoying. This is because make.simmap(...,nsim>1) returns a list of trees; and thus lapply(trees,make.simmap,...,nsim>1) returns a list of lists. Various attempts to first unlist and then relist left me more & more annoyed - but the following hack seems to do the trick:

ff<-function(tree,x){
   zz<-make.simmap(tree,x,nsim=10)
   class(zz)<-NULL
   zz
}

mtrees<-unlist(sapply(trees,ff,x,simplify=FALSE), recursive=FALSE)
class(mtrees)<-"multiPhylo"

I have now added this to the latest version of make.simmap, & also built a new version of phytools (phytools 0.2-36), which can be downloaded & installed from source.

Instead of only taking a single tree as input, this tree can take a list of trees (an object of class "multiPhylo") & will automatically generate nsim stochastic character maps per input tree.

OK, here's a demo of the new version using a set of 15,001 trees from the posterior distribution of a real Bayesian run (thanks Graham Reynolds), and a simulated binary character with states a and b.

> packageVersion("phytools")
[1] ‘0.2.36’
> trees<-read.nexus("posterior.sample.trees")
> trees
15001 phylogenetic trees
> # too many, let's randomly subsample
trees<-sample(trees[5001:15001],100)
> trees
100 phylogenetic trees
> # ok, now generate 10 stochastic maps for each tree
> mtrees<-make.simmap(trees,x,nsim=10,message=FALSE)
> mtrees
1000 phylogenetic trees
> # now let's visualize the variability
> # again, 1000 is too many
> layout(matrix(1:100,10,10,byrow=TRUE))
> cols<-setNames(c("blue","red"),letters[1:2])
> plotSimmap(mtrees[0:99*10+1],cols,pts=F,ftype="off")
Waiting to confirm page change...

Cool. Now let's try describe.simmap:

> XX<-describe.simmap(mtrees)
1000 trees with a mapped discrete character with states:
 a, b

trees have 15.509 changes between states on average

changes are of the following types:
       a,b   b,a
x->y 9.078 6.431

mean total time spent in each state is:
               a           b    total
raw  176.5796569 145.4146415 321.9943
prop   0.5482752   0.4517248   1.0000

The times & state changes computed by describe.simmap will be correct - however the posterior probabilities for ancestral nodes (here, XX$ace) will not because different trees in the posterior sample have different nodes & node numbers.

Nonetheless, cool!

Combining results from a phylogenetic regression on multiple trees in a posterior sample

A recent query on the R-sig-phylo email list asked:

"Imagine that I run a PIC analysis on two traits using 1000 post-burn-in trees. What would be the best way to summarize these results? Average p-values across all analyses? Perhaps a specific method to combine the resulting probabilities e.g. Fisher's test?"

A little bit of discussion ensued, but my suggestion read as follows:

...we should probably (begin by) combin(ing) the variance among estimates obtained from different trees in the posterior sample with the sampling variance of any single estimate... (to get the total variance of our parameter estimates for hypothesis testing)

One fairly sensible way to do this is to compute the variance due to phylogenetic uncertainty as the variance among estimates obtained from the trees of the posterior sample; and then to compute the sampling variance as the mean variance of the estimator from each tree; and then add the two variances them. The standard error of our estimate (computed as the mean across trees) is the square-root of this variance. To conduct a hypothesis test on the regression coefficient, then, you would compute the mean across trees and then the ratio of the parameter and its standard error should have a t-distribution with n-2 degrees of freedom for n taxa (not contrasts).

To do this from a practical perspective from trees in 'multiPhylo' object (here "trees") and data in x & y for a simple bivariate regression, we do the following:

# first define the following custom function
ff<-function(tree,x,y){
   pic.x<-pic(x,tree)
   pic.y<-pic(y,tree)
   fit<-lm(pic.y~pic.x-1)
   setNames(c(coef(fit),vcov(fit)),c("beta","var(beta)"))
}
# now apply to all trees in your sample
BB<-t(sapply(trees,ff,x,y))
# total variance in beta estimated by
varBeta<-var(BB[,"beta"])+mean(BB[,"var(beta)"])
t.beta<-mean(BB[,"beta"])/sqrt(varBeta)
P.beta<-2*pt(abs(t.beta),df=length(trees[[1]]$tip)-2, lower.tail=FALSE)

I think that's right.

Now, someone wisely pointed out that this approach is somewhat inelegant in combining Bayesian, maximum likelihood, & frequentist hypothesis testing. I can't argue with this and would agree that a better & more elegant solution would be to simultaneously sample trees with branch lengths & the parameters of our evolutionary model for trait evolution from their joint posterior probability distribution. Presently, this is impractical.

Picking a tree or set of trees at random from a "multiPhylo" object

A friend recently asked:

"Does anybody know how to, in R, select a single tree at random from a posterior distribution of trees (e.g., generated via MrBayes)?"

This is not too hard. A set of trees read into memory using read.tree or read.nexus is just a list of trees with the class attribute set to "multiPhylo". To pick one at random we can just do:

random.tree<-sample(trees,size=1)[[1]]

We include the index [[1]] to (non-recursively) unlist the object returned by sample.

If we want to sample more than one, say 100 random trees from a posterior sample of 1,000, the procedure is the same. So we can just do:

random.trees<-sample(trees,size=100)

That's it.

phytools expands into industrial networking solutions....

OK, the post title is an April Fools' Day joke; however imagine my horror at discovering yesterday that phytools is in fact a real, brand new, private commercial company formed in 2013 as a subsidiary of PHYTEC America (larger screenshot, here). What either company actually does is a mystery to me, but I'm shocked that PHYTEC thinks they will get away with this! phytools.com is warned that they should expect to hear from phytools.org legal department** in the coming hours or days!

**phytools.org legal department has just informed me that they don't exist....