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...

phenogram(prunedBarbetTree, wingL)

The analyses will we run in this exercise focus on one particular anole ecomorph, the crown giant. These large anoles have many unique features. We can create a vector that identifies whether each species is a crown-giant or not.

crownGiantSpecies<-c("equestris", "luteogularis", "smallwoodi", "noblei", "baracoae", "baleatus", "barahonae", "ricordii", "cuvieri", "garmani")
                        
isCrownGiant <- anoleData[,1] %in% crownGiantSpecies

ecomorph <- sapply(isCrownGiant, function(x) if(x) "cg" else "other")
ecomorph <- as.factor(ecomorph)
names(ecomorph) <- anoleData[,1]

We can then use stochastic character mapping, as we learned earlier, to reconstruct the evolution of this character.

ecomorphTree <- make.simmap(anoleTree, ecomorph, model="ER")
## make.simmap is sampling character histories conditioned on the transition matrix
## 
## Q =
##                cg       other
## cg    -0.09903616  0.09903616
## other  0.09903616 -0.09903616
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##    cg other 
##   0.5   0.5
## Done.
plot(ecomorphTree)
## no colors provided. using the following legend:
##      cg   other 
## "black"   "red"

The data we will analyse is body size, which appears in our matrix as snout-vent length (svl).

svl<-anoleData[,2]
names(svl)<-anoleData[,1]

We will fit both a single-rate and a two-rate Brownian motion model to the data using phytools.

fitBrownie<-brownie.lite(ecomorphTree, svl)
fitBrownie
## ML single-rate model:
##  s^2 se  a   k   logL
## value    0.1362  0.02    4.0659  2   -4.7004 
## 
## ML multi-rate model:
##  s^2(cg) se(cg)  s^2(other)  se(other)   a   k   logL    
## value    0.9108  0.444   0.0877  0.0141  4.0169  3   6.5917
## 
## P-value (based on X^2): 0 
## 
## R thinks it has found the ML solution.

This result shows that the two-rate model fits the data significantly better than a one-rate model. The reconstructed rate of evolution in the fitted model is much higher for crown-giant anoles than for other ecomorphs.

Multi-regime models using OUwie

Now let’s try to fit a multi-regime OU model using OUwie.

OUwie requires a special data frame with three columns: species name, the regime coding, and the variable being analysed.

svl_df<-data.frame(anoleData[,1], ecomorph, svl)

We can run single models in OUwie - for example, fitting a model where each regime has a different mean but all other parameters are shared.

fitOUM<-OUwie(ecomorphTree, svl_df,model="OUM",simmap.tree=TRUE)
## Initializing... 
## Finished. Begin thorough search... 
## Finished. Summarizing results.
fitOUM
## 
## Fit
##        lnL      AIC     AICc model ntax
##  0.8891881 6.221624 6.642676   OUM  100
## 
## 
## Rates
##                 cg     other
## alpha    0.2015132 0.2015132
## sigma.sq 0.1344489 0.1344489
## 
## Optima
##                 cg      other
## estimate 5.4354545 4.03223034
## se       0.3895511 0.09203497
## 
## Arrived at a reliable solution

To really test this model, though, we can fit a whole set of models using OUwie. Note that the first two models here match the two models that we just compared using brownie.lite - and, in fact, the two programs obtain identical likelihoods and parameter estimates.

fitBM1<-OUwie(ecomorphTree, svl_df,model="BM1",simmap.tree=TRUE)
## Initializing... 
## Finished. Begin thorough search... 
## Finished. Summarizing results.
fitBM1
## 
## Fit
##        lnL      AIC     AICc model ntax
##  -4.700433 13.40087 13.52458   BM1  100
## 
## Rates
##                 cg     other
## alpha           NA        NA
## sigma.sq 0.1361648 0.1361648
## 
## Optima
##                 cg other
## estimate 4.0659177     0
## se       0.1079716     0
## 
## Arrived at a reliable solution
fitBMS<-OUwie(ecomorphTree, svl_df,model="BMS",simmap.tree=TRUE, root.station=FALSE)
## Initializing... 
## Finished. Begin thorough search... 
## Finished. Summarizing results.
fitBMS
## 
## Fit
##       lnL       AIC      AICc model ntax
##  6.591643 -7.183286 -6.933286   BMS  100
## 
## Rates
##                 cg      other
## alpha           NA         NA
## sigma.sq 0.9109018 0.08765503
## 
## Optima
##                  cg      other
## estimate 4.01691283 4.01691283
## se       0.08714475 0.08714475
## 
## Arrived at a reliable solution
fitOUM<-OUwie(ecomorphTree, svl_df,model="OUM",simmap.tree=TRUE)
## Initializing... 
## Finished. Begin thorough search... 
## Finished. Summarizing results.
fitOUM
## 
## Fit
##        lnL      AIC     AICc model ntax
##  0.8891881 6.221624 6.642676   OUM  100
## 
## 
## Rates
##                 cg     other
## alpha    0.2015132 0.2015132
## sigma.sq 0.1344489 0.1344489
## 
## Optima
##                 cg      other
## estimate 5.4354545 4.03223034
## se       0.3895511 0.09203497
## 
## Arrived at a reliable solution
fitOUMV<-OUwie(ecomorphTree, svl_df,model="OUMV",simmap.tree=TRUE)
## Initializing... 
## Finished. Begin thorough search... 
## Finished. Summarizing results.
fitOUMV
## 
## Fit
##       lnL       AIC      AICc model ntax
##  8.160023 -6.320045 -5.681748  OUMV  100
## 
## 
## Rates
##                    cg        other
## alpha    2.075506e-09 2.075506e-09
## sigma.sq 6.449443e-01 8.784855e-02
## 
## Optima
##                 cg      other
## estimate 5.6497018 4.01337143
## se       0.8540612 0.08728189
## 
## Arrived at a reliable solution
fitOUMA<-OUwie(ecomorphTree, svl_df,model="OUMA",simmap.tree=TRUE)
## Initializing... 
## Finished. Begin thorough search... 
## Finished. Summarizing results.
fitOUMA
## 
## Fit
##       lnL       AIC      AICc model ntax
##  22.08282 -34.16565 -33.52735  OUMA  100
## 
## 
## Rates
##                    cg        other
## alpha    2.954741e-07 2.846468e-09
## sigma.sq 7.968891e-02 7.968891e-02
## 
## Optima
##                 cg    other
## estimate 4.9865764 4.007425
## se       0.1370558 0.082891
## 
## Arrived at a reliable solution
fitOUMVA<-OUwie(ecomorphTree, svl_df,model="OUMVA",simmap.tree=TRUE)
## Initializing... 
## Finished. Begin thorough search... 
## Finished. Summarizing results.
fitOUMVA
## 
## Fit
##       lnL       AIC      AICc model ntax
##  25.74329 -39.48658 -38.58335 OUMVA  100
## 
## 
## Rates
##                   cg        other
## alpha    0.007752205 7.776475e-05
## sigma.sq 0.014514143 8.465950e-02
## 
## Optima
##                 cg      other
## estimate 5.0006339 4.00619220
## se       0.1114816 0.08499375
## 
## Arrived at a reliable solution

Now we can compare the fit of all of these models using AIC weights.

ouwie_aicc<-c(fitBM1$AICc, fitBMS$AICc, fitOUM$AICc, fitOUMV$AICc, fitOUMA$AICc, fitOUMVA$AICc)
names(ouwie_aicc)<-c("fitBM1", "fitBMS", "fitOUM", "fitOUMV", "fitOUMA", "fitOUMVA")
aic.w(ouwie_aicc)
##     fitBM1     fitBMS     fitOUM    fitOUMV    fitOUMA   fitOUMVA 
## 0.00000000 0.00000012 0.00000000 0.00000007 0.07391835 0.92608146

The crown-giant anoles do seem to have their own adaptive regime, and that regime has a different level of constraint compared to the other ecomorphs as a group.

Multi-regime models using bayou

Bayou uses Bayesian rjMCMC to identify regimes, and can do so in the absence of particular hypotheses. Let’s see if bayou can detect the unique selective regime of the crown-giant anoles.

First we need to set up priors for all model parameters. In general, it is worth thinking about this a bit. But I am lazy and will just use the defaults from the bayou tutorial on github.

prior <- make.prior(anoleTree, dists=list(dalpha="dhalfcauchy", dsig2="dhalfcauchy",dsb="dsb", dk="cdpois", dtheta="dnorm"), param=list(dalpha=list(scale=1), dsig2=list(scale=1), dk=list(lambda=15, kmax=200), dsb=list(bmax=1,prob=1), dtheta=list(mean=mean(svl), sd=2)))

Now, run the MCMC

#fit1 <- bayou.makeMCMC(anoleTree, svl, SE=0, model="OU", prior=prior, new.dir="../temp3/", plot.freq=5000, ticker.freq=1000)
#fit1$run(50000)

Bayou stores results in a set of files in your working directory. To explore convergence and see your results you have to read them back in.

#chainOU <- fit1$load()
#chainOU <- set.burnin(chainOU, 0.3)
#summary(chainOU)
#plot(chainOU, auto.layout=FALSE)

Now let’s visualize shift locations on the tree

#par(mfrow=c(1,1))
#plotSimmap.mcmc(chainOU, burnin = 0.3, pp.cutoff = 0.3)

How do these shifts look in a phenogram?

#phenogram.density(anoleTree, svl, burnin=0.3, chainOU, pp.cutoff=0.3)

Results may vary - but for these it does look like there is something unique about the crown giants! Note that many features of these results suggest that we should have run our Bayou analyses for a much longer time.

As was noted during the exercise, a more complete tutorial of Bayou can be see online here.