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")
## 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
Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.