Monday, August 20, 2012

Discrete-time Brownian simulation & visualization

I was just thinking about how to simulate & visualize evolution under the threshold model (a model of quantitative genetics discussed in two previous posts: 1, 2). I thought I would like to plot in discrete time, but then color all the segments of liability evolution with a color that depended on the implied state for the threshold character.

As a means to an end in getting there, I just wrote a simple function to simulate & plot discrete-time Brownian evolution on the tree. This is in the function bmPlot (code here).

The function is pretty simple. Basically, the user just provides an input tree and a number of generations to rescale and round the branches. Obviously, if the tree already has discrete-time (i.e., whole integer) branch lengths, you can just provide the total tree length here. The the function conducts discrete time Brownian simulations on each branch of the phylogeny (in much the same way as my example code, here). Finally, the function returns the states at the tips and all internal nodes in a named vector.

Let's see how it works:

> library(phytools)
> source("bmPlot.R")
> tree<-pbtree(n=10)
> bmPlot(tree)
         t1          t2          t3
1.35526072  1.13925561        . . .


No comments:

Post a Comment

Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.