Here is an alternative version that is more than twice as fast (helpful
for large trees), and only a little bit slower than vcv.phylo():
vcvPhylo2<-function(tree,anc.nodes=T){
n<-length(tree$tip.label)
h<-nodeHeights(tree)[order(tree$edge[,2]),2]
h<-c(h[1:n],0,h[(n+1):length(h)])
M<-mrca(tree,full=anc.nodes)[c(1:n,anc.nodes*(n+2:tree$Nnode)),c(1:n,anc.nodes*(n+2:tree$Nnode))]
C<-matrix(h[M],nrow(M),ncol(M))
if(anc.nodes)
rownames(C)<-colnames(C)<-c(tree$tip.label,n+2:tree$Nnode)
else rownames(C)<-colnames(C)<-tree$tip.label
return(C)
}
All the best. - Liam
--
Liam J. Revell
University of Massachusetts Boston
web: http://faculty.umb.edu/liam.revell/
email: liam.rev...@umb.edu
blog: http://phytools.blogspot.com
On 11/21/2011 6:16 PM, ppi...@uniroma3.it wrote:
Hi Liam,
it works perfectly! Thankyou very much!
And....it is elegant...
Hi Paolo.
I think this code fragment from my phytools function anc.trend does what
you want (here, I have made it into its own function), although it might
not be the most elegant way to do it:
vcvPhylo<-function(phy,anc.nodes=T){
if(anc.nodes){
D<-dist.nodes(phy)
ntips<-length(phy$tip.label)
Cii<-D[ntips+1,]
C<-D; C[,]<-0
for(i in 1:nrow(D)) for(j in 1:ncol(D))
C[i,j]<-(Cii[i]+Cii[j]-D[i,j])/2
dimnames(C)[[1]][1:length(phy$tip)]<-dimnames(C)[[2]][1:length(phy$tip)]<-phy$tip.label
C<-C[c(1:ntips,(ntips+2):nrow(C)),c(1:ntips,(ntips+2):ncol(C))]
return(C)
} else return(vcv.phylo(phy))
}
All the best, Liam
--
Liam J. Revell
University of Massachusetts Boston
web: http://faculty.umb.edu/liam.revell/
email: liam.rev...@umb.edu
blog: http://phytools.blogspot.com
On 11/21/2011 12:09 PM, ppi...@uniroma3.it wrote:
> Hi all,
> Having a tree (NOT ultrametric in my case) I would
> like to extract the phylogenetic covariance matrix
> including both tips and nodes (like dist.nodes for
> patristic distances).
> vcv returns the cov matric just for tips.
> Someone has a fast way?
>
> Thankyou very much
>
> Paolo
>
> _______________________________________________
> R-sig-phylo mailing list
> R-sig-phylo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
>
_______________________________________________
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
