Monday, June 6, 2016

Update to function for computing multiple VCV matrices based on a mapped trait

I just posted a very small update to the mostly internally used phytools function, multiC, which computes a series of phylogenetic VCV matrices depending on a mapped state. The modification allows the user to obtain a list of VCV matrices with variances & covariances of internal nodes included.

So, for instance:

library(phytools)
tree
## 
## Phylogenetic tree with 10 tips and 9 internal nodes.
## 
## Tip labels:
##  t2, t6, t10, t5, t1, t4, ...
## 
## The tree includes a mapped, 2-state discrete character with states:
##  A, B
## 
## Rooted; includes branch lengths.
colors=setNames(c("blue","red"),c("A","B"))
plot(tree,colors,lwd=3)
add.simmap.legend(x=2.66,y=2.1,colors=colors,prompt=FALSE,
    vertical=TRUE,shape="circle")

plot of chunk unnamed-chunk-1

## VCV
vcv(tree)
##       t2   t6  t10   t5   t1   t4   t8   t9   t3   t7
## t2  0.29 0.11 0.11 0.11 0.11 0.11 0.11 0.11 0.11 0.00
## t6  0.11 1.39 0.53 0.33 0.33 0.33 0.33 0.33 0.33 0.00
## t10 0.11 0.53 0.60 0.33 0.33 0.33 0.33 0.33 0.33 0.00
## t5  0.11 0.33 0.33 2.27 1.75 0.76 0.39 0.39 0.39 0.00
## t1  0.11 0.33 0.33 1.75 2.20 0.76 0.39 0.39 0.39 0.00
## t4  0.11 0.33 0.33 0.76 0.76 1.21 0.39 0.39 0.39 0.00
## t8  0.11 0.33 0.33 0.39 0.39 0.39 2.88 1.99 1.36 0.00
## t9  0.11 0.33 0.33 0.39 0.39 0.39 1.99 2.90 1.36 0.00
## t3  0.11 0.33 0.33 0.39 0.39 0.39 1.36 1.36 1.74 0.00
## t7  0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.42
## VCV by mapped state
multiC(tree)
## $B
##       t2       t6  t10       t5       t1   t4       t8       t9   t3
## t2  0.29 0.110000 0.11 0.110000 0.110000 0.11 0.110000 0.110000 0.11
## t6  0.11 0.869376 0.53 0.330000 0.330000 0.33 0.330000 0.330000 0.33
## t10 0.11 0.530000 0.60 0.330000 0.330000 0.33 0.330000 0.330000 0.33
## t5  0.11 0.330000 0.33 2.252992 1.732992 0.76 0.390000 0.390000 0.39
## t1  0.11 0.330000 0.33 1.732992 2.182992 0.76 0.390000 0.390000 0.39
## t4  0.11 0.330000 0.33 0.760000 0.760000 1.21 0.390000 0.390000 0.39
## t8  0.11 0.330000 0.33 0.390000 0.390000 0.39 2.116673 1.990000 1.36
## t9  0.11 0.330000 0.33 0.390000 0.390000 0.39 1.990000 2.720538 1.36
## t3  0.11 0.330000 0.33 0.390000 0.390000 0.39 1.360000 1.360000 1.74
## t7  0.00 0.000000 0.00 0.000000 0.000000 0.00 0.000000 0.000000 0.00
##             t7
## t2  0.00000000
## t6  0.00000000
## t10 0.00000000
## t5  0.00000000
## t1  0.00000000
## t4  0.00000000
## t8  0.00000000
## t9  0.00000000
## t3  0.00000000
## t7  0.08773829
## 
## $A
##     t2       t6 t10         t5         t1 t4        t8        t9 t3
## t2   0 0.000000   0 0.00000000 0.00000000  0 0.0000000 0.0000000  0
## t6   0 0.520624   0 0.00000000 0.00000000  0 0.0000000 0.0000000  0
## t10  0 0.000000   0 0.00000000 0.00000000  0 0.0000000 0.0000000  0
## t5   0 0.000000   0 0.01700815 0.01700815  0 0.0000000 0.0000000  0
## t1   0 0.000000   0 0.01700815 0.01700815  0 0.0000000 0.0000000  0
## t4   0 0.000000   0 0.00000000 0.00000000  0 0.0000000 0.0000000  0
## t8   0 0.000000   0 0.00000000 0.00000000  0 0.7633267 0.0000000  0
## t9   0 0.000000   0 0.00000000 0.00000000  0 0.0000000 0.1794623  0
## t3   0 0.000000   0 0.00000000 0.00000000  0 0.0000000 0.0000000  0
## t7   0 0.000000   0 0.00000000 0.00000000  0 0.0000000 0.0000000  0
##            t7
## t2  0.0000000
## t6  0.0000000
## t10 0.0000000
## t5  0.0000000
## t1  0.0000000
## t4  0.0000000
## t8  0.0000000
## t9  0.0000000
## t3  0.0000000
## t7  0.3322617
## VCV by mapped state with internal
multiC(tree,internal=TRUE)
## $B
##       t2       t6  t10       t5       t1   t4       t8       t9   t3
## t2  0.29 0.110000 0.11 0.110000 0.110000 0.11 0.110000 0.110000 0.11
## t6  0.11 0.869376 0.53 0.330000 0.330000 0.33 0.330000 0.330000 0.33
## t10 0.11 0.530000 0.60 0.330000 0.330000 0.33 0.330000 0.330000 0.33
## t5  0.11 0.330000 0.33 2.252992 1.732992 0.76 0.390000 0.390000 0.39
## t1  0.11 0.330000 0.33 1.732992 2.182992 0.76 0.390000 0.390000 0.39
## t4  0.11 0.330000 0.33 0.760000 0.760000 1.21 0.390000 0.390000 0.39
## t8  0.11 0.330000 0.33 0.390000 0.390000 0.39 2.116673 1.990000 1.36
## t9  0.11 0.330000 0.33 0.390000 0.390000 0.39 1.990000 2.720538 1.36
## t3  0.11 0.330000 0.33 0.390000 0.390000 0.39 1.360000 1.360000 1.74
## t7  0.00 0.000000 0.00 0.000000 0.000000 0.00 0.000000 0.000000 0.00
## 12  0.11 0.110000 0.11 0.110000 0.110000 0.11 0.110000 0.110000 0.11
## 13  0.11 0.330000 0.33 0.330000 0.330000 0.33 0.330000 0.330000 0.33
## 14  0.11 0.530000 0.53 0.330000 0.330000 0.33 0.330000 0.330000 0.33
## 15  0.11 0.330000 0.33 0.390000 0.390000 0.39 0.390000 0.390000 0.39
## 16  0.11 0.330000 0.33 0.760000 0.760000 0.76 0.390000 0.390000 0.39
## 17  0.11 0.330000 0.33 1.732992 1.732992 0.76 0.390000 0.390000 0.39
## 18  0.11 0.330000 0.33 0.390000 0.390000 0.39 1.360000 1.360000 1.36
## 19  0.11 0.330000 0.33 0.390000 0.390000 0.39 1.990000 1.990000 1.36
##             t7   12   13   14   15   16       17   18   19
## t2  0.00000000 0.11 0.11 0.11 0.11 0.11 0.110000 0.11 0.11
## t6  0.00000000 0.11 0.33 0.53 0.33 0.33 0.330000 0.33 0.33
## t10 0.00000000 0.11 0.33 0.53 0.33 0.33 0.330000 0.33 0.33
## t5  0.00000000 0.11 0.33 0.33 0.39 0.76 1.732992 0.39 0.39
## t1  0.00000000 0.11 0.33 0.33 0.39 0.76 1.732992 0.39 0.39
## t4  0.00000000 0.11 0.33 0.33 0.39 0.76 0.760000 0.39 0.39
## t8  0.00000000 0.11 0.33 0.33 0.39 0.39 0.390000 1.36 1.99
## t9  0.00000000 0.11 0.33 0.33 0.39 0.39 0.390000 1.36 1.99
## t3  0.00000000 0.11 0.33 0.33 0.39 0.39 0.390000 1.36 1.36
## t7  0.08773829 0.00 0.00 0.00 0.00 0.00 0.000000 0.00 0.00
## 12  0.00000000 0.11 0.11 0.11 0.11 0.11 0.110000 0.11 0.11
## 13  0.00000000 0.11 0.33 0.33 0.33 0.33 0.330000 0.33 0.33
## 14  0.00000000 0.11 0.33 0.53 0.33 0.33 0.330000 0.33 0.33
## 15  0.00000000 0.11 0.33 0.33 0.39 0.39 0.390000 0.39 0.39
## 16  0.00000000 0.11 0.33 0.33 0.39 0.76 0.760000 0.39 0.39
## 17  0.00000000 0.11 0.33 0.33 0.39 0.76 1.732992 0.39 0.39
## 18  0.00000000 0.11 0.33 0.33 0.39 0.39 0.390000 1.36 1.36
## 19  0.00000000 0.11 0.33 0.33 0.39 0.39 0.390000 1.36 1.99
## 
## $A
##     t2       t6 t10         t5         t1 t4        t8        t9 t3
## t2   0 0.000000   0 0.00000000 0.00000000  0 0.0000000 0.0000000  0
## t6   0 0.520624   0 0.00000000 0.00000000  0 0.0000000 0.0000000  0
## t10  0 0.000000   0 0.00000000 0.00000000  0 0.0000000 0.0000000  0
## t5   0 0.000000   0 0.01700815 0.01700815  0 0.0000000 0.0000000  0
## t1   0 0.000000   0 0.01700815 0.01700815  0 0.0000000 0.0000000  0
## t4   0 0.000000   0 0.00000000 0.00000000  0 0.0000000 0.0000000  0
## t8   0 0.000000   0 0.00000000 0.00000000  0 0.7633267 0.0000000  0
## t9   0 0.000000   0 0.00000000 0.00000000  0 0.0000000 0.1794623  0
## t3   0 0.000000   0 0.00000000 0.00000000  0 0.0000000 0.0000000  0
## t7   0 0.000000   0 0.00000000 0.00000000  0 0.0000000 0.0000000  0
## 12   0 0.000000   0 0.00000000 0.00000000  0 0.0000000 0.0000000  0
## 13   0 0.000000   0 0.00000000 0.00000000  0 0.0000000 0.0000000  0
## 14   0 0.000000   0 0.00000000 0.00000000  0 0.0000000 0.0000000  0
## 15   0 0.000000   0 0.00000000 0.00000000  0 0.0000000 0.0000000  0
## 16   0 0.000000   0 0.00000000 0.00000000  0 0.0000000 0.0000000  0
## 17   0 0.000000   0 0.01700815 0.01700815  0 0.0000000 0.0000000  0
## 18   0 0.000000   0 0.00000000 0.00000000  0 0.0000000 0.0000000  0
## 19   0 0.000000   0 0.00000000 0.00000000  0 0.0000000 0.0000000  0
##            t7 12 13 14 15 16         17 18 19
## t2  0.0000000  0  0  0  0  0 0.00000000  0  0
## t6  0.0000000  0  0  0  0  0 0.00000000  0  0
## t10 0.0000000  0  0  0  0  0 0.00000000  0  0
## t5  0.0000000  0  0  0  0  0 0.01700815  0  0
## t1  0.0000000  0  0  0  0  0 0.01700815  0  0
## t4  0.0000000  0  0  0  0  0 0.00000000  0  0
## t8  0.0000000  0  0  0  0  0 0.00000000  0  0
## t9  0.0000000  0  0  0  0  0 0.00000000  0  0
## t3  0.0000000  0  0  0  0  0 0.00000000  0  0
## t7  0.3322617  0  0  0  0  0 0.00000000  0  0
## 12  0.0000000  0  0  0  0  0 0.00000000  0  0
## 13  0.0000000  0  0  0  0  0 0.00000000  0  0
## 14  0.0000000  0  0  0  0  0 0.00000000  0  0
## 15  0.0000000  0  0  0  0  0 0.00000000  0  0
## 16  0.0000000  0  0  0  0  0 0.00000000  0  0
## 17  0.0000000  0  0  0  0  0 0.01700815  0  0
## 18  0.0000000  0  0  0  0  0 0.00000000  0  0
## 19  0.0000000  0  0  0  0  0 0.00000000  0  0

That's more or less it.

Data were simulated as follows:

tree<-rtree(n=10)
tree$edge.length<-round(tree$edge.length,2)
Q<-matrix(c(-1,1,1,-1),2,2)
rownames(Q)<-colnames(Q)<-LETTERS[1:2]
tree<-sim.history(tree,Q)

No comments:

Post a Comment