Pages

Tuesday, October 28, 2014

A colleague contacted me today about “alternative implementations” of geiger's fitContinuous function
that he could use to check parameter estimation, and maybe to speed calculations as fitContinuous
can be somewhat slow (although is evidently quite robust now) for some models.

Well, I don't know about speeding things up - but it is pretty straightforward to write a wrapper function that
combines phytools brownie.lite single-rate model with geiger rescale.phylo to fit
many of the models of fitContinuous. The following is a simple demo:

Running system.time with any of these examples (except for model="BM") does not show any
improvement of this ostensibly lighter implementation - but this is probably due to lots of unnecessary internals in
brownie.lite for fitting a multi-rate model (which is not actually done here). Improving on this is thus
pretty simple (but I will not do so here, nor now).

Monday, October 27, 2014

“I have a tree and discrete data (number of gene copies, for many genes)
and would like to use the fitDiscrete function in geiger, or something
similar. However, I would like to estimate the parameters given all of the
datasets, not just with the data for each gene. For instance, if I was
using the "delta” model to vary rates across the tree, I would like this
delta value to reflect some sort of summary value across all datasets. Does
anyone have an idea as to how this could be accomplished or perhaps point
me in the right direction?“

Brian O'Meara pointed out that it would be straighforward to write a likelihood
function that just summed the log-likelihoods across characters for different
tree transformations and then wrapped the function in an optimizer. The following
illustrates this explicitly.

First load packages - for this we need ape, phytools (just for simulation here),
and geiger (for our tree rescaling):

## load packages
library(ape)
library(phytools)

## Loading required package: maps

library(geiger)

Now let's simulate some data that would be analogous to what we'd load from file
in the empirical case:

Note that it is arbitrary, and unnecessary, that I have simulated using the same
transition matrix Q for each data column here.

Finally, our tree-transformation fitting function. Here, I've used ace
from ape internally to compute the likelihood of each character for each tree
transformation. One could instead use fitDiscrete (which is slower, but
perhaps more robust) or function in the phangorn package.

Saturday, October 25, 2014

“I am currently using the ltt95 function on a
"multiPhylo" object. I would like to known if there is a
solution to reverse the time axis from the root of the tree to
present?”

Indeed this is not too difficult.
ltt95
already returns an object of class "ltt95" invisibly. This
object is basically a matrix containing the plotted data for x and
y in columns, with a few different attributes to tell the S3 plotting
method for objects of its class how to handle the matrix. So if we were naive
as to the functionality of our plotting method, we could just manipulate the
contents of that matrix to get the plot that we wanted. phytools makes things
even easier, however, as there is an argument xaxis in the
plotting method for objects of class "ltt95" that automates this
all for us.

Here's how we can use it to control the direction of x in various ways:

Saturday, October 11, 2014

A phytools user recently wrote in to ask if she could color the terminal
edges of a phylogeny plotted using plotTree.wBars different
colors depending on a discrete trait. Well, we
already
know that it's possible to plot a tree with a mapped discrete character.
So to accomplish this, we just need to paint the terminal edges of the tree
using the phytools function paintSubTree.

OK, so now we have data approximating the empirical data and we can attempt to
generate our simulation. We will do this by first coloring all the terminal
edges of the tree according to our discrete character, and then, having done
that, we will run plotTree.wBars using the method="plotSimmap".

About this blog

This web-log chronicles the development of new tools for phylogenetic analyses in the phytools R package. Unless you a reading a very recent page of the blog, I recommend that you install the latest CRAN version of phytools (or latest beta release) before attempting to replicate any of the analyses of this site. That is because the linked functions may be archived, and very likely have been replaced by newer versions.