This tutorial is about fitting multi-rate Brownian motion models (using phytools) and multi-regime OU models (using the OUwie and bayou packages).

Multi-rate models using phytools

The first exercise is designed to explore a model developed by O’Meara et al. (2006; Evolution) in which the rate of Brownian evolution takes different values in different parts of a phylogeny.

For this analysis we will use the function brownie.lite in the phytools package. This function allows us to input a tree with painted rate “regimes” - for instance, from a stochastically mapped discrete character - and then fit a model in which the rate differs depending on the mapped regime.

First, let’s make sure we have the anole and barbet trees & datasets that we used earlier in the course:

  1. anole.data.csv
  2. Anolis.tre
  3. Barbet data
  4. Barbet tree

Make sure you have these files in your working directory, and then read them in.

anoleTree<-read.tree("Anolis.tre")
anoleData<-read.csv("anole.data.csv")

barbetTree<-read.nexus("BarbetTree.nex")
barbetData<-read.csv("Barbetdata.csv")

Let’s first compare rates of evolution of body size between anoles and barbets. We will use different indices of body size in each case (svl for anoles and wing length for barbets); this is justified as we are on a log scale so rates are still comparable.

# rescale Anole tree to m.y.
anoleTimeTree<-anoleTree
anoleTimeTree$edge.length <- anoleTimeTree$edge.length * 50
svl<-anoleData[,2]
names(svl)<-anoleData[,1]
name.check(anoleTimeTree, svl)
## [1] "OK"
wingL<-barbetData[,2]
names(wingL)<-barbetData[,1]
nn<-name.check(barbetTree, wingL)
prunedBarbetTree<-drop.tip(barbetTree, nn[[1]])

twoTrees<-c(anoleTimeTree, prunedBarbetTree)
twoChars<-list(svl, wingL)
ratebytree(twoTrees, twoChars)
## ML common-rate model:
##  s^2  a[1]   a[2]    k   logL
## value    0.0024  4.0659  4.4843  3   19.8546 
## SE       3e-04   0.1024  0.138   
## 
## ML multi-rate model:
##   s^2[1] s^2[2]   a[1]   a[2]    k   logL
## value    0.0028  0.0016  4.0659  4.4843  4   21.5421 
## SE       4e-04   4e-04   0.1085  0.1123  
## 
## Likelihood ratio: 3.3751 
## P-value (based on X^2): 0.0662 
## 
## One or the other optimization may not have converged.
phenogram(anoleTimeTree, svl)
## Optimizing the positions of the tip labels...