New function to simulate multiple evolutionary correlations between characters in different parts of a tree

In phytools there are a couple of different functions based on my 2009 Evolution paper with Dave Collar for fitting multiple evolutionary correlations between characters to different parts of a phylogeny (Revell & Collar 2009). I just added a new function that can simulate multivariate evolution with arbitrarily different VCV matrices for the Brownian process for different states of a mapped discrete character (regimes) on the tree. To do this, I used a similar approach to the existing phytools function fastBM. In this case, however, I drew Brownian changes for multiple characters from their joint multivariate normal distributions. I did this separately for each mapped state; and then I added the changes along each edge together. The function (sim.corrs) is here.

We can test it out by generating under a multiple correlation model; and then fitting the model of Revell & Collar (2009) using evol.vcv. Let's try it out.

First, let's simulate a tree & character history using the functions pbtree and sim.history:

> tree<-sim.history(pbtree(n=100,scale=1), Q=matrix(c(-1,1,1,-1),2,2),anc=1)
> cols<-c("blue","red"); names(cols)<-c(1,2)
> plotSimmap(tree,cols,ftype="off",pts=F)

Next, load sim.corrs from source, create VCV matrices for simulation, and simulate:

> source("sim.corrs.R")
> vcv<-list(matrix(c(1,0.9,0.9,1),2,2), matrix(c(1,0.1,0.1,1.0),2,2))
> names(vcv)<-c(1,2)
> vcv
$`1`
     [,1] [,2]
[1,]  1.0  0.9
[2,]  0.9  1.0
$`2`
     [,1] [,2]
[1,]  1.0  0.1
[2,]  0.1  1.0
> X<-sim.corrs(tree,vcv)

Finally, let's fit our multiple VCV matrix model using evol.vcv. (We could also use evolvcv.lite to test a bunch of different models - including the common rates, different correlation model we used for simulation. More details here.)

> fit<-evol.vcv(tree,X)
> fit
$R.single
          [,1]      [,2]
[1,] 1.0603320 0.4342028
[2,] 0.4342028 0.8906066
$logL1
[1] -96.63233
$k1
[1] 5
$R.multiple
, , 1
          [,1]      [,2]
[1,] 0.9804452 0.8125698
[2,] 0.8125698 0.8371796
, , 2
          [,1]      [,2]
[1,] 1.1098840 0.1077934
[2,] 0.1077934 0.9188769
$logL.multiple
[1] -71.52816
$k2
[1] 8
$P.chisq
[1] 7.213247e-11
$convergence
[1] "Optimization has converged."

Cool. Both this and the other new functions will be in a new (& soon upcoming) version of phytools.