The AnthroTree Website & "The Comparative Approach in Evolutionary Anthropology and Biology" by C. Nunn
Table of Contents - Chapter Overview
Previous subsection - Next Section
Trait data can be simulated under a variety of different evolutionary models. Here, we will use R to simulate character evolution using both Brownian motion and Ornstein-Uhlenbeck (OU) as evolutionary models of character change. Some of the code was taken from Paradis (2006). Open R and load the 'ape' package (Paradis et al. 2004; see Section 1.1.2 for installation instructions):
library(“ape”)
The following code simulates a trait evolving through time under Brownian motion:
time_steps = 100X.b = replicate(10, cumsum(c(0, rnorm(time_steps - 1))))
If you type "X.b", you will see a matrix containing 100 rows of data and 10 columns. Each column represents a different evolutionary history of the trait under Brownian motion. Starting at the top, you will see that the values start at zero and change through time.
The following code simulates a trait under the OU model:
alpha = 0.1 # restraining force on trait evolutionsim.ou = function() {
x.ou = numeric(time_steps)for (i in 2:time_steps)x.ou[i] = -alpha * x.ou[i - 1] + rnorm(1)x.ou # returns the value of x
}
Type "X.ou" to view the simulated trait values under the OU model. Next, compare the variance of trait values under the two models by plotting the output of the different runs:
layout(matrix(1:2, 2, 1))yl = range(X.b)y2 = range(X.b)matplot(X.b, ylim = yl, type = "l", col = 1, main = "Brownian")matplot(X.ou, ylim = y2, type = "l", col = 1, main = "OU")
You should see two plots; one of them shows 10 simulations of a Brownian motion model of evolution for a random character, the other one 10 simulations for an OU model. Both are based on 100 time steps. It should be apparent that Brownian motion exhibits unbounded variation that accumulates with time, whereas the OU model describes more bounded variation around an optimum value. More specifically, OU has a restraining force, the parameter "alpha", that tends to pull the trait back to the optimum. We can verify that the variance of the trait values is higher under the Brownian motion model relative to the OU model by typing:
var(X.b[time_steps ,])var(X.ou[time_steps ,])
You should see that the variance of the Brownian motion model is substantially higher than in the OU model (about 37 times higher in the simulation that I ran; your results will be slightly different due to stochastic effects).
You can explore the effect of the parameter "alpha" by changing it to an arbitrary value between 0 (Brownian motion) and 1 in the script. As "alpha" approaches 0, evolution becomes more Brownian motion-like, with Brownian motion itself equivalent to the case where "alpha = 0".
Notice that unlike the Mesquite example above, we are not simulating trait evolution on a tree. However, this can also be done in R. An example is given in Section 6.1.
Maddison, W. P., and D.R. Maddison. 2006. Mesquite: a modular system for evolutionary analysis. http://mesquiteproject.org.
Paradis, E. 2006. Analysis of Phylogenetics and Evolution with R. Springer, New York.
Paradis E., J. Claude, and K. Strimmer. 2004. APE: analyses of phylogenetics and evolution in R language. Bioinformatics 20: 289-290.
R Development Core Team (2010). R: A language and environment for statistical computing. R Foundation for
Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org.
Comments