Three different ways to calculate the among-species variance-covariance matrix for multiple traits on a phylogeny

Here are three different ways to compute the among-species variance-covariance matrix for multiple characters on the tree. If we assume that the characters are evolving by Brownian motion, this matrix is an unbiased estimate of the instantaneous diffusion matrix of the evolutionary process.

The three options are (1) the GLS estimating equation (used in Revell 2009); (2) the method of equation 5 in Butler et al. (2000); and (3) we can use PICs (Felsenstein 1985).

Here is code to do each of these in R:

## compute evolutionary VCV matrix using method of
## Revell (2009)
A<-matrix(1,nrow(X),1)%*%apply(X,2,fastAnc,tree=tree)[1,]
V1<-t(X-A)%*%solve(vcv(tree))%*%(X-A)/(nrow(X)-1)

## compute evolutionary VCV matrix using method of
## Butler et al. (2000)
Z<-solve(t(chol(vcv(tree))))%*%(X-A)
V2<-t(Z)%*%Z/(nrow(X)-1)

## compute the evolutionary VCV matrix using pics
Y<-apply(X,2,pic,phy=tree)
V3<-t(Y)%*%Y/nrow(Y)

OK - now let's try it with some simulated data:

> ## load packages
> require(phytools)
Loading required package: phytools
Loading required package: ape
Loading required package: maps
Loading required package: rgl
> ## simulate tree & data
> simV<-matrix(c(1,1.2,0,1.2,2,0,0,0,3),3,3)
> simV ## generating covariance matrix
     [,1] [,2] [,3]
[1,]  1.0  1.2    0
[2,]  1.2  2.0    0
[3,]  0.0  0.0    3
> tree<-pbtree(n=100,scale=1) ## simulate tree
> X<-sim.corrs(tree,simV) ## simulate data
>
> ## Revell (2009)
> A<-matrix(1,nrow(X),1)%*%apply(X,2,fastAnc,tree=tree)[1,]
> V1<-t(X-A)%*%solve(vcv(tree))%*%(X-A)/(nrow(X)-1)
>
> ## Butler et al. (2000)
> Z<-solve(t(chol(vcv(tree))))%*%(X-A)
> V2<-t(Z)%*%Z/(nrow(X)-1)
>
> ## pics
> Y<-apply(X,2,pic,phy=tree)
> V3<-t(Y)%*%Y/nrow(Y)
>
> ## compare
> V1
          [,1]     [,2]      [,3]
[1,] 1.0896188 1.375287 0.1157135
[2,] 1.3752866 2.280321 0.2996510
[3,] 0.1157135 0.299651 2.9491762
> V2
          [,1]     [,2]      [,3]
[1,] 1.0896188 1.375287 0.1157135
[2,] 1.3752866 2.280321 0.2996510
[3,] 0.1157135 0.299651 2.9491762
> V3
          [,1]     [,2]      [,3]
[1,] 1.0896188 1.375287 0.1157135
[2,] 1.3752866 2.280321 0.2996510
[3,] 0.1157135 0.299651 2.9491762

Obviously, although we can see that these calculations seem different, they produce the same results. The PIC method is a bit trickier, but at least in methods (1) & (2) it is fairly easy to see why this is the case. In case (1) we do X^TC^-1X; whereas in case (2) we do Z = C^-1/2X (as Cholesky decomposition is a kind of matrix square root), followed by Z^TZ which is the same as X^TC^-1/2C^-1/2X (in which each of C^-1/2 are the lower & upper Cholesky matrices of the inverse of C, respectively), which is then equivalent to X^TC^-1X. Obviously, we need to center X on the vector of ancestral values & divide by n - 1 in either case.