Exercise 3.1: Simulating Brownian motion in R

This short tutorial gives some simple approaches that can be used to simulate Brownian evolution in continuous and discrete time, in the absence of and on a phylogenetic tree.

Brownian motion is a stochastic continuous-time random walk model in which changes from one time to the next are random draws from some distribution with mean 0.0 and variance σ². The expected variance under Brownian motion increases linearly through time with instantaneous rate σ².

Brownian motion is very easy to simulate. To start off, let’s simulate a single instance of Brownian motion for 100 generations of discrete time in which the variance of the diffusion process is σ² = 0.01 per generation. In this case, we will draw our evolutionary changes from a normal distribution; however it’s worth noting that (due to the CLT) regardless of the distribution, evolution will proceed by Brownian motion as the width of our timesteps decrease towards zero!

t<-0:100 # time
sig2<-0.01
## first, simulate a set of random deviates
x<-rnorm(n=length(t)-1,sd=sqrt(sig2))
## now compute their cumulative sum
x<-c(0,cumsum(x))
plot(t,x,type="l",ylim=c(-2,2))

The trick here is that we draw random normal deviates for the change over each t time intervals; and then to get the state of our chain at each time interval we just compute the cumulative sum of all the individual changes using cumsum. We can do this because the distribution of the changes under Brownian motion is invariant and does not depend on the state of the chain.

We can also easily do a whole bunch of simulations like this at once, using the same conditions:

nsim<-100
X<-matrix(rnorm(n=nsim*(length(t)-1),sd=sqrt(sig2)),nsim,length(t)-1)
X<-cbind(rep(0,nsim),t(apply(X,1,cumsum)))
plot(t,X[1,],xlab="time",ylab="phenotype",ylim=c(-2,2),type="l")
apply(X[2:nsim,],1,function(x,t) lines(t,x),t=t)

## NULL

To see how the outcome depends on σ², let’s compare the result when we divide sig2 by 10:

X<-matrix(rnorm(n=nsim*(length(t)-1),sd=sqrt(sig2/10)),nsim,length(t)-1)
X<-cbind(rep(0,nsim),t(apply(X,1,cumsum)))
plot(t,X[1,],xlab="time",ylab="phenotype",ylim=c(-2,2),type="l")
apply(X[2:nsim,],1,function(x,t) lines(t,x),t=t)

## NULL

The expected variance under Brownian motion is just σ² × t. To see this easiest, we can just do the following. Here I will use 10,000 simulations for 100 generations under the same conditions to “smooth out” our result:

nsim<-10000
X<-matrix(rnorm(n=nsim*(length(t)-1),sd=sqrt(sig2)),nsim,length(t)-1)
X<-cbind(rep(0,nsim),t(apply(X,1,cumsum)))
v<-apply(X,2,var)
plot(t,v,type="l",xlab="time",ylab="variance among simulations")

var(X[,length(t)]) # this should be about 1.00

## [1] 1.02573

OK, now let’s try to simulate using Brownian motion up the branches of a tree. We first need a tree, so let’s simulate one using pbtree in phytools:

library(phytools)

## Loading required package: ape

## Loading required package: maps

## 
##  # maps v3.1: updated 'world': all lakes moved to separate new #
##  # 'lakes' database. Type '?world' or 'news(package="maps")'.  #

t<-100 # total time
n<-30 # total taxa
b<-(log(n)-log(2))/t
tree<-pbtree(b=b,n=n,t=t)

## simulating with both taxa-stop (n) and time-stop (t) is
## performed via rejection sampling & may be slow
## 
##   170 trees rejected before finding a tree

plotTree(tree,mar=c(3.1,0.1,0.1,0.1))
axis(1)

Now, to simulate on all the branches of the tree we just need to simulate on each branch separately, and then “shift” each daughter branch by the end state of it’s parent. Here is an example of a continuous time simulation & visualization using canned functions.

## simulate Brownian evolution on a tree with fastBM
x<-fastBM(tree,sig2=sig2,internal=TRUE)
## visualize Brownian evolution on a tree
phenogram(tree,x,spread.labels=TRUE,spread.cost=c(1,0))

Written by Liam J. Revell & Mike Alfaro. Last updated 28 Jun. 2016