mvBM

‘mvBM’ (multiple variance Brownian motion) is a function that estimates lineage-specific rates of evolution and infers ancestral values accordingly. Ancestral values can be inferred using either ML or MCMC. Lineage-specific rates can be plotted using MCMC analysis.

References:

Smaers JB, Mongle CS, Kandler A (2016) A multiple variance Brownian motion framework for estimating variable rates and inferring ancestral states. Biological Journal of the Linnean Society. 118 (1): 78-94.

Smaers JB & Mongle CS (2017) On the accuracy and theoretical underpinnings of the multiple variance Brownian motion approach for estimating variable rates and inferring ancestral states. Biological Journal of the Linnean Society. 121(1): 229-238.

Below a case study on the ancestral estimation of primate brain size exemplifying the use of mvBM.

Prepare the analysis

require(phytools); require(geiger); require(nlme); require(evomap)
data<-primateData
tree<-primateTree
data<-log(data) 
tree<-treedata(tree,data,sort=T,warnings=T)$phy
data<-as.data.frame(treedata(tree,data,sort=T,warnings=T)$data)

 

Using maximum likelihood

x<-data$Brain
BMsigma2<-ace(x,tree,method="REML")$sigma2[1] # get BM sigma2 ('ace' requires the 'ape' package)
mvBMresults<-mvBM(x,tree,BMsigma2) # calculate rescaled branch lengths using mvBM
tree_mvBM<-tree
tree_mvBM$edge.length<-mvBMresults$rBL # create new tree with rescaled branch lengths
plot(tree_mvBM) # plot mvBM tree
ace(x,tree_mvBM,method="REML") # get ancestral estimates using mvBM tree

 

Using MCMC

#Get mvBM rescaled tree
x<-data$Brain
BMsigma2<-ace(x,tree,method="REML")$sigma2[1] # get BM sigma2 ('ace' requires the 'ape' package)
mvBMresults<-mvBM(dat,tree,BMsigma2) 
tree_mvBM<-tree; tree_mvBM$edge.length<-mvBMresults$rBL
plot(tree_mvMB)

#Use mvBM rescaled tree in MCMC analysis
dat<-data$Brain; names(dat)<-rownames(data)
iterations<-500000 
model_mvBM<-anc.Bayes(tree_mvBM,dat,ngen=iterations) # 'anc.Bayes' requires the 'phytools' package; this may take a few minutes
  MCMC_mvBM_sigma2<-model_mvBM$mcmc[(iterations*0.002):nrow(model_mvBM$mcmc),2] #include a 0.2 burnin
  MCMC_mvBM_anc<-colMeans(model_mvBM$mcmc[(iterations*0.002):nrow(model_mvBM$mcmc),])[3:(3+tree$Nnode-1)] #These are the ancestral estimates
  MCMC_mvBM_logLik<-model_mvBM$mcmc[(iterations*0.002):nrow(model_mvBM$mcmc),tree$Nnode+3]

hist(MCMC_mvBM_logLik) #One of the requirements of Bayesian optmization is that this is normally distributed. If this is not the case, try increasing the number of iterations.

 

Lineage-specific rate analysis using the MCMC output

Set groups for which rates of evolution should be plotted. Note that the rate can be estimated for any individual branch or any group of branches.

Human<-which(tree$tip.label=="Homo_sapiens")
nonHumanPrimates<-setdiff(getTips(tree,findMRCA(tree,c("Homo_sapiens","Callicebus_moloch"))),Human)

Get the branches for these groups

Branches_Humans<-getEdges(tree,which(tree$tip.label=="Homo_sapiens"))[1]
Branches_nonHumanPrimates<-getEdges(tree,findMRCA(tree,nonHumanPrimates))

To get the rates:

mvBM.getRate(tree,tree_mvBM,branches=Branches_Humans,sigma2Distr = MCMC_mvBM_sigma2)
mvBM.getRate(tree,tree_mvBM,branches=Branches_nonHumanPrimates,sigma2Distr = MCMC_mvBM_sigma2)

To plot the rates:

#The rate of evolution for the human lineage:
mvBM.plotRate(tree,tree_mvBM,branches=Branches_Humans,sigma2Distr = MCMC_mvBM_sigma2,xlim=c(0,0.07),ylim=c(0,300),col="red",adjust=6,main="")
 par(new=TRUE)
#The average rate of evolution for all other lineages:
mvBM.plotRate(tree,tree_mvBM,branches=Branches_nonHumanPrimates,sigma2Distr = MCMC_mvBM_sigma2,xlim=c(0,0.07),ylim=c(0,300),col="blue",adjust=6,main="")