MCMC method: analyzing the posterior sample

As I mentioned in my previous post, the posterior sample can be treated as would any posterior sample from a Bayesian analysis. For instance, we should expect to compute parameter estimates, effective sample sizes (which takes into account the degree of correlation among our samples), and 95% credible intervals mostly in the typical way. We can do many of these things within R using the MCMC diagnostics package {coda} (Plummer et al. 2010); however, for those users that prefer to use the Java program Tracer (by Rambaut & Drummond), this is fine to. We just need to first write our MCMC samples to a text file:

> write.table(x=res$mcmc,file="mcmc.sample.txt",sep="\t", row.names=FALSE,quote=FALSE)

The image below is a screenshot of the program Tracer, showing the trace of the log-likelihood of our analysis from last time:



There is no explicit manual for Tracer, but it is very easy to use. In additional posts, I will discuss how to do these analyses in {coda} as well as how to summarize the posterior sample for the position of the rate shift (something that cannot be done easily using either Tracer or {coda}).

Ecomorph phylomorphospace

Check out the neat plot created by Luke Mahler using my phylomorphospace() function (available here). Click on the image at right to enlarge. In this plot are the phenotypes and phylogenetic relationships among 100 Greater Antillean Anolis lizards. According to Luke, PC1 is limb lengths (more or less); and PC2 is overall size. The tips of the tree are color coded by ecomorph. The phylomorphospace graphing method is based on Sidlauskas (2008).

Luke made a few modifications/additions to my code to obtain this plot. First, he filled and shrunk the points that are the inferred states at internal nodes. He did this by changing the following line of code:

points(x=c(X[tree$edge[i,1],1],X[tree$edge[i,2],1]),
  y=c(X[tree$edge[i,1],2],X[tree$edge[i,2],2]),
  col=c(con$col.node[as.character(tree$edge[i,1]),1],
  con$col.node[as.character(tree$edge[i,2]),1]))

to this:

points(x=c(X[tree$edge[i,1],1],X[tree$edge[i,2],1]),
  y=c(X[tree$edge[i,1],2],X[tree$edge[i,2],2]),
  col=c(con$col.node[as.character(tree$edge[i,1]),1],
  con$col.node[as.character(tree$edge[i,2]),1]),pch=16,cex=0.8)

To plot the large, color coded tip nodes, he merely created the vector color.code.e, containing the species color coded by ecomorph; and then plotted his raw data for PCs 1 & 2 on top of his phylomorphospace graph as follows:

> points(X[,1],X[,2],col=(color.code.e),pch=16,cex=1.7)

Thanks to Luke for sharing this.

MCMC method for rate variation

So, in addition to tracking new development activity with this page, I thought I would also post a few words about some of the existing functions that I am distributing from my R phylogenetics page.

Today, I will discuss the function evol.rate.mcmc(). This function takes a continuous character vector and a phylogeny with branch lengths, and then uses Bayesian MCMC to identify the phylogenetic position of a shift in the evolutionary rate of phenotypic evolution through time. Evolutionary rates and shift-points between rates are sampled from their joint posterior distribution. The underlying method for this function is described in a manuscript presently in review:

Revell, L. J., D. L. Mahler, P. R. Peres-Neto, and B. D. Redelings. Submitted. A new method for identifying exceptional evolutionary diversification.

To give a quick example of how this works, take the data and tree presented at right. These data were generated with a high rate (σ²=20) on the red branches of the tree; and a low rate (σ²=1.0) on the blue branches. The data are in a labeled vector, x, and the phylogenetic tree is given in an {ape} "phylo" object tree.

To load the source and required dependencies, we will first enter:

> source("evol.rate.mcmc.R")
> require(ape); require(coda);

The package {coda} will be used to analyze and compute various diagnostics for our MCMC chain (next post) - in this case our posterior sample - but we may as well load it now.

To run the function, we simply give R (some variant of) the following command:

> res<-evol.rate.mcmc(tree,x,ngen=20000,control=list(sd1=1.0, sd2=1.0,sda=0.5,sdlnr=2.0,rand.shift=0.05,print=1000))

This tells our MCMC to run under the following conditions:

ngen: the number of generations for our MCMC run is 20,000;
sd1: standard deviation of the proposal distribution for rate 1 is 1.0;
sd2: standard deviation of the proposal distribution for rate 2 is 1.0;
sda: standard deviation on the proposal distribution for the ancestral state at the root node of the tree is 0.5.
rand.shift: probability of a random shift to a different position in the tree is 0.05 [otherwise the MCMC will proceed as a random walk with a default variance];
print: print frequency for the chain (to screen - all generations are stored in res$mcmc).

We also need to assign a prior probability distribution to our parameters. I chose to specify a uniform prior on the absolute values of rates 1 and 2, the root state, and the position of the rate shift; however, I have used a log-normal prior distribution for the ratio of rate 1 / rate 2, centered on 0.0, and with a standard deviation given by sdlnr - in this case we set this to 2.0.

The run returns a list with two components:

$mcmc: is our matrix of parameters from our MCMC run, including the sample of shift points between rates on the tree; and
$tips: is a list of vectors containing all the tips in the clade defined as rate 1.

In my next post, I will discuss analysis of the posterior sample from this run.

Luke Harmon phyloseminar

Luke Harmon (pictured right), who has been very active and productive in phylogenetic method development over the past few years, will be giving the next edition in the web seminar series on phylogenetics at phyloseminar.org. Luke's seminar will take place tomorrow at 1PM Eastern and can be viewed live (instructions here). If you can't attend the live version of the seminar, it will also be archived here (along with a list of other great prior seminars in this series). Luke will be presenting a talk on a variety of new cutting-edge methods that are being developed in his laboratory at the University of Idaho.

New "beta" read.simmap() function available

I have added the capacity to read "nexus" style SIMMAPv1.0 tree files to my new read.simmap() function. To remind the reader, this function (v0.2 on my beta version page) is an improvement over earlier versions in that it allows the user to retain the order and the times spent in each mapped state along each edge of the tree. My earlier read.simmap() function (v0.1) destroyed the ordering of the states on each edge, storing only the total branch length in each state for each edge. The present version can read both "nexus" and "phylip" style tree files; as well as multiple trees (in which case, the function will create a modified "multiPhylo" object in memory).

For example, to load and use the present version:

> source("read.simmap.R")
> tree<-"((A:{aqua,0.3:terr,0.4:aqua,0.3},B:{aqua,1.0}):{aqua,1.0},(C:{aqua,0.25:terr,0.75},D:{aqua,1.0}):{aqua,1.0});"
> mapped.tree<-read.simmap(text=tree,format="phylip")

The modified "phylo" object that is created by read.simmap() in this case has the following structure:

> mapped.tree
$edge
[,1] [,2]
[1,] 5 6
[2,] 6 1
[3,] 6 2
[4,] 5 7
[5,] 7 3
[6,] 7 4
$Nnode
[1] 3
$tip.label
[1] "A" "B" "C" "D"
$edge.length
[1] 1 1 1 1 1 1
$maps
$maps[[1]]
aqua
1
$maps[[2]]
aqua terr aqua
0.3 0.4 0.3
$maps[[3]]
aqua
1
$maps[[4]]
aqua
1
$maps[[5]]
aqua terr
0.25 0.75
$maps[[6]]
aqua
1
$mapped.edge
state
edge aqua terr
5,6 1.00 0.00
6,1 0.60 0.40
6,2 1.00 0.00
5,7 1.00 0.00
7,3 0.25 0.75
7,4 1.00 0.00
attr(,"class")
[1] "phylo"

Ancestral character estimation, with a trend

In response to the recent R-sig-phylo email listserve discussion thread on simulating Brownian evolution with a trend, Pasquale Raia inquired to the list about whether it would be straightforward to estimate ancestral states under this model. In practice, of course, it is generally not possible to fit the BM with a trend model (BMT for brevity in this post) - nor to estimate ancestral states under BMT - for conditions in which all the tips of the tree (and our data for trait values) are contemporaneous. However, for conditions in which the tips of the tree are effectively not contemporaneous (for instance, if we have fossil species in the tree or if the evolutionary process is punctuational), then we are in a theoretically much better position to estimate evolutionary trends.

The way I decided to approach this problem was very simple. First, let's consider the expected distribution of the data at tips and internal nodes. Well, we know from theory (and Gene Hunt's listserve post) that the theoretical distribution of tip and (by simple extension) internal node values under BMT is:

> x<-mvrnorm(n=1,mu=a+diag(C)*mu,Sigma=sig2*C)

Here, a is the root value, mu is the trend per unit time, and sig2 is the variance of the stochastic component of the model. C is a matrix of the heights above the root for each terminal and internal node in the tree (excluding the root). For example, in the 4 taxon phylogram plotted above, the matrix C is given by:

> C
     t1     t3     t2     t4     6     7
t1 0.984819 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
t3 0.000000 1.1680555 0.4643410 0.3490977 0.3490977 0.4643410
t2 0.000000 0.4643410 0.4908593 0.3490977 0.3490977 0.4643410
t4 0.000000 0.3490977 0.3490977 0.5612080 0.3490977 0.3490977
6 0.000000 0.3490977 0.3490977 0.3490977 0.3490977 0.3490977
7 0.000000 0.4643410 0.4643410 0.3490977 0.3490977 0.4643410

Labels "t1", "t2", etc., correspond to tips of the tree; whereas "6", "7", etc., are internal node numbers from phy$edge.

Now, for a given set of internal and terminal values, C, and values for a, mu, and sig2, we can calculate the probability density using dmnorm() in the {mnormt} package. Of course, that means we can also calculate the likelihood the ancestral states and evolutionary parameters, conditioned on the values for terminal species and the tree.

I have just put up (on my R-phylogenetics page, here) a very preliminary version of a function that uses the built-in R multivariate optimizer, optim(), to do this.

To test it out, I did the following.

First, I loaded fastBM() and anc.trend() from source:

> source("fastBM.R"); source("anc.trend.R")

Now, I simulated a random tree and data using rtree() from the {ape} package and fastBM():

> tree<-rtree(20)
> x.all<-fastBM(tree,a=0,mu=2,sig2=1,internal=TRUE)

Next, I estimated ancestral traits by maximizing the likelihood using anc.trend():

> result<-anc.trend(phy=tree,x=x.all[1:length(tree$tip)], maxit=2000)

This function runs quite slow (not surprising, since in the above example I am doing a 22 dimension optimzation), so it has been somewhat hard to play around with so far. I can report that when optim() converges, mu is equivalent to the trend parameter estimated using fitContinuous(model="trend") in {geiger}.

In the above example, I obtained:

> result
$ace
     21     22     23     24     25     26
-0.2420098 0.1641928 1.2565491 3.0897363 1.8626764 1.0436854
     27     28     29     30     31     32
2.0145277 3.7843751 4.8605234 6.3360333 1.7286019 2.3675538
     33     34     35     36     37     38
3.1769088 5.3097681 3.3040749 2.1776491 0.2706769 0.9679458
     39
1.1140074
$mu
[1] 1.693923
$sig2
[1] 0.3092271
$logL
[1] -13.80569
$convergence
[1] 0
$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

Although mu estimated in my function, and mu from fitContinuous() are the same; sig2 differed (invariably) by exactly n/(2n-2). This is the same as the ratio of the number of tips/(number of tips + number of nodes) in an unrooted, bifurcating tree.

In addition, the estimated ancestral states are fairly closely correlated with their known conditions. For instance, when I ran the code given above:

> cor(c(0,x.all[21:38):length(x.all)]),result$ace)
[1] 0.5973801

I should probably emphasize that the above results likely oversell the method because in general this cannot be used for species with only contemporaneous tips (rtree() tends to produce phylogenies with wildly non-contemporaneous tips); and furthermore the simulated trend in this example was quite strong. Nonetheless - this is cool!

BM simulation with a trend

Dave Bapst just posted a query to the R-sig-phylo about simulating Brownian evolution with a trend. I take this to mean that evolutionary changes along branches in the tree are drawn from a random normal distribution with variances sig2*tree$edge.length and means mu*tree$edge.length. [Initially, in my first response to the query, I did not multiply mu by tree$edge.length; but that doesn't make any sense.] I have now added this capability to v0.2 of my function fastBM() which can be downloaded from my R phylogenetics page here.

After loading the function from source, we can just type:

> x<-fastBM(tree,a=0,mu=0.5,sig2=2,internal=FALSE)

at the command prompt. This will return a n x 1 vector of character values (in x for the n species in tree.

To test out this function, let's try simulating data and then fitting a "trend" model using {geiger}'s fitContinuous() function. [Obviously, since these are stochastic simulations, individual results will vary.]

> tree<-rtree(500) # first simulate stochastic tree
> x<-fastBM(tree,mu=0.5,sig2=2) # now data on the tree
> fitContinuous(phy=tree,data=x,model="trend") # now fit model
Fitting trend model:
$Trait1
$Trait1$lnl
[1] -900.9902
$Trait1$mean
[1] 1.087359
$Trait1$beta
[1] 1.87744
$Trait1$mu
[1] 0.4968745
$Trait1$aic
[1] 1807.980
$Trait1$aicc
[1] 1808.029
$Trait1$k
[1] 3

Seems to work!

Note that Gene Hunt also pointed out (correctly, of course) that this can be done using the {MASS} function mvrnorm(). To see the discussion thread on this topic, check out the archived R-sig-phylo emails here.