VCV.array: Create a 2D or 3D variance-covariance matrix from a phylogeny

Description

The function turns a phylogeny into a variance-covariance matrix, as in the function vcv.phylo in the 'ape' package but can also return a 3D array retaining the individual branch lengths contributing to the shared branch lengths. This is useful for handling some branch length transformations, such as kappa, and has a lower overhead than repeatedly calling vcv.phylo on a phylogeny after transforming the vector of edge lengths.

Usage

VCV.array(phy, dim=2, compact=TRUE)

Arguments

phy

An object of class 'phylo'.

dim

Either 2, for a standard VCV matrix, or 3, for an array of branch lengths.

compact

A logical vector indicating the form to use for a 3D array.

Value

When dim = 2, a variance covariance matrix of class 'VCV.array' of dimension nTips by nTips with dimnames set from the tip labels of the phylogeny.
When dim = 3, a 3 dimensional array of class 'VCV.array' with dimensions of the number of taxa in the phylogeny for the rows and columns and either the maximum number of branches on the root to tip path or the number of internal nodes as the depth, depending on the setting of compact. The rows and columns are named using the tip labels from the phylogeny and the depth only named with node numbers if compact is TRUE.

Details

The compact form of the 3D array uses a shortened third dimension, which is only long enough to hold the maximum number of shared branches between root and tip for each pair of tips. Zeros are used to pad out this depth vector for tip pairs with shorter paths. The non-compact form returns 3D array showing, for each pair of tips and each node in the tree, either 0 if the node is not shared or the appropriate edge length if the node is shared. Note that, for maximally unbalanced trees, the size of the two forms will be identical.

The algorithm for the noncompact form is faster than for the compact form but it has very high memory overheads on big trees. The 2 dimensional algorithm is at least twice as fast as vcv.phylo on trees up to 2500 tips.

The apply function can be easily used to collapse the array down to a standard VCV matrix, as in the example.

Examples

Run this code

tree <- rcoal(8)
tree.VCV <- vcv.phylo(tree)
tree.VCVA <- VCV.array(tree)

# reconstruct a simple VCV array
tree.VCVA.reduced <- apply(tree.VCVA, c(1,2), sum, na.rm=TRUE)

# minimal differences between the two
all((tree.VCVA.reduced - tree.VCV) < 1e-10)

# a kappa transformation of 0.5
apply(tree.VCVA ^ 0.5, c(1,2), sum, na.rm=TRUE)