# some utility functions # written by Liam J. Revell 2011/2012/2013 # function rounds the branch lengths of the tree & applies rounding to simmap tree # written by Liam J. Revell 2012 roundBranches<-function(tree,digits=0){ if(class(tree)=="multiPhylo"){ trees<-lapply(tree,roundBranches) class(trees)<-"multiPhylo" return(trees) } else { tree$edge.length<-round(tree$edge.length,digits) if(!is.null(tree$maps)){ for(i in 1:nrow(tree$edge)){ temp<-tree$maps[[i]]/sum(tree$maps[[i]]) tree$maps[[i]]<-temp*tree$edge.length[i] } } if(!is.null(tree$mapped.edge)){ a<-vector() for(i in 1:nrow(tree$edge)) a<-c(a,names(tree$maps[[i]])) a<-unique(a) tree$mapped.edge<-matrix(data=0,length(tree$edge.length),length(a),dimnames=list(apply(tree$edge,1,function(x) paste(x,collapse=",")),state=a)) for(i in 1:length(tree$maps)) for(j in 1:length(tree$maps[[i]])) tree$mapped.edge[i,names(tree$maps[[i]])[j]]<-tree$mapped.edge[i,names(tree$maps[[i]])[j]]+tree$maps[[i]][j] } return(tree) } } # function applies the branch lengths of a reference tree to a second tree, including mappings # written by Liam J. Revell 2012 applyBranchLengths<-function(tree,edge.length){ if(class(tree)=="multiPhylo"){ trees<-lapply(tree,applyBranchLengths,edge.length=edge.length) class(trees)<-"multiPhylo" return(trees) } else { tree$edge.length<-edge.length if(!is.null(tree$maps)){ for(i in 1:nrow(tree$edge)){ temp<-tree$maps[[i]]/sum(tree$maps[[i]]) tree$maps[[i]]<-temp*tree$edge.length[i] } } if(!is.null(tree$mapped.edge)){ a<-vector() for(i in 1:nrow(tree$edge)) a<-c(a,names(tree$maps[[i]])) a<-unique(a) tree$mapped.edge<-matrix(data=0,length(tree$edge.length),length(a),dimnames=list(apply(tree$edge,1,function(x) paste(x,collapse=",")),state=a)) for(i in 1:length(tree$maps)) for(j in 1:length(tree$maps[[i]])) tree$mapped.edge[i,names(tree$maps[[i]])[j]]<-tree$mapped.edge[i,names(tree$maps[[i]])[j]]+tree$maps[[i]][j] } return(tree) } } # function to compute phylogenetic VCV using joint Pagel's lambda # written by Liam Revell 2011 phyl.vcv<-function(X,C,lambda){ C<-lambda.transform(lambda,C) invC<-solve(C) a<-matrix(colSums(invC%*%X)/sum(invC),ncol(X),1) A<-matrix(rep(a,nrow(X)),nrow(X),ncol(X),byrow=T) V<-t(X-A)%*%invC%*%(X-A)/(nrow(C)-1) return(list(C=C,R=V,alpha=a)) } # lambda transformation of C # written by Liam Revell 2011 lambda.transform<-function(lambda,C){ if(lambda==1) return(C) else { V<-diag(diag(C)) C<-C-V C.lambda<-(V+lambda*C) return(C.lambda) } } # likelihood function for joint estimation of lambda for multiple traits # written by Liam Revell 2011/2012 likMlambda<-function(lambda,X,C){ # compute R, conditioned on lambda temp<-phyl.vcv(X,C,lambda); C<-temp$C; R<-temp$R; a<-temp$alpha # prep n<-nrow(X); m<-ncol(X); D<-matrix(0,n*m,m) for(i in 1:(n*m)) for(j in 1:m) if((j-1)*n=(1-tol))){ a<-as.numeric(names(which(z>=(1-tol)))) if(length(a)==1) Nodes[i,2]<-a else { Nodes[i,2]<-a[1] warning("polytomy detected; some node matches may be arbitrary") } } } } return(Nodes) } # function reorders the columns of mapped.edge from a set of simmap trees # written by Liam J. Revell 2013 orderMappedEdge<-function(trees,ordering=NULL){ f1<-function(tree,ordering){ mapped.edge<-matrix(0,nrow(tree$mapped.edge),length(ordering), dimnames=list(rownames(tree$mapped.edge),ordering)) mapped.edge[,colnames(tree$mapped.edge)]<-tree$mapped.edge tree$mapped.edge<-mapped.edge return(tree) } f2<-function(tree) colnames(tree$mapped.edge) if(class(trees)=="phylo") states<-colnames(trees$mapped.edge) else if(class(trees)=="multiPhylo") states<-unique(as.vector(sapply(trees,f2))) else stop("trees should be an object of class 'phylo' or 'multiPhylo'") if(length(ordering)>1) if(length(intersect(states,ordering))n)) result$nodes<-sisters[which(sisters>n)] else if(any(sisters>n)) result$nodes<-tree$node.label[sisters[which(sisters>n)]-n] if(any(sisters<=n)) result$tips<-tree$tip.label[sisters[which(sisters<=n)]] return(result) } }