msBayes Basics I: Introduction

I've been working quite a bit with msBayes recently, and thought I would share some of the stuff that I've learnt, with an emphasis on the practical side of things. Over a series of posts, I am going to summarize the steps required to carry out a basic msBayes analysis, provide some "lessons learned" and tips for each of these steps, as well as provide some supplemental Python scripts that I find greatly facilitate the flow of data from one stage of the pipeline to another. This first post serves as a general introduction to the approach and the program, and a general overview of the procedures involved in an msBayes analysis.

Most of you are probably already quite familiar with this method of analysis, but for those of you that are not, perhaps a brief explanation may be in order. msBayes is a suite of tools that plug into a pipeline for Approximate Bayesian Computation (ABC) phylogeographic analysis, estimation and hypothesis testing. The ABC approach is a computationally efficient approach to complex estimation and inference problems. At its most basic, it works as follows:

1. A set of summary statistics is calculated on an observed data set.
2. Multiple sets of data are simulated using parameters drawn from a prior distribution.
   
3. The same same set of summary statistics is calculated for each of these simulated data sets.
4. The parameters that result in simulated data sets with summary statistics that are equal to or within some threshold distance of the summary statistics of the observed data set are recorded.
5. The posterior probability of the parameter values is given by the proportion of those values in the final set of accepted parameters.
   

Each of these steps involves multiple levels of complexities and and a smorgasbord of variant approaches, of course, from the full nuts-and-bolts of the simulation model to various flavours of importance/rejection sampling (step 4), but the core concept remains the same: simulate the data under a range of parameter values, and keep those parameter values that result in simulated data sufficiently close to the observed data. This approach is extremely attractive: for many complex problems and questions, data can be simulated relatively easily, whereas sometimes even formulating the corresponding likelihood statement can be pretty challenging. Even if the likelihood can be formulated, computing it tends to be extremely expensive in terms of CPU power and time as well as memory, while simulating data is relatively cheap as far as those resources go. However, be warned: the ABC approach does exact a cost in storage space; all those samples from the prior have to be retained until they can be processed, and I have had analyses that walk into the realm of 50-60 GB.

msBayes can be used to answer a broad variety of questions using ABC, but a major one, and the one which I am currently using it for, is to test the hypothesis of simultaneous divergence amongst a set of taxa co-distributed across a geographical barrier. Specifically, I am asking "given a bunch of sequences of pairs of populations of different species on either side of a barrier, what is the posterior probability a single distinct divergence time between each pair?"

Using msBayes to answer this question involves the following steps:

1. Assemble the observed data files.
2. Create the configuration file to drive the simulations from the prior.
3. Run the simulator to produce a large number of samples from the prior.
4. Calculate the summary statistics on the observed data.
5. Regularize the columns / fields in the summary statistics file and the simulated data.
6. Carry out the acceptance/rejection sampling on the prior samples to produce samples from the posterior.
7. Summarize the posterior probability distribution of discrete parameter values (such as number of divergence times) by assessing the proportion of those values in the posterior samples.
8. Weight posterior samples of continuous parameter using local linear regression before summarizing their distribution.

Each of these steps will be covered in detail in future posts, and I will also be providing some scripts that I found really useful in helping to carry out some of these steps:

• Part II: Data Assembly, Configuration File Preparation, and Generating Samples from the Prior
• Part III: Approaches to Sampling from the Posterior
• Part IV: Preparing Your Files for Sampling from the Posterior
• Part V: Rejection Sampling to Extract the Posteriors from Your Priors
• Part VI: Tolerances and Local-Linear Regression Weighting of the Posteriors