# This function creates a phylomorsphace plot of Sidlauskas (2006) # Written by Liam J. Revell 2010, 2011, 2012, 2013 phylomorphospace<-function(tree,X,A=NULL,label=TRUE,control=list(),...){ # some minor error checking if(class(tree)!="phylo") stop("tree object must be of class 'phylo.'") if(nrow(X)!=length(tree$tip)) stop("X must contain the same number of rows as species in tree.") if(ncol(X)!=2) warning("X has more than 2 columns. Using only the first 2 columns.") # check to see if A should be estimated if(is.null(A)){ # load dependencies A<-matrix(NA,tree$Nnode,2,dimnames=list(as.character(length(tree$tip.label)+1:tree$Nnode),colnames(X))) for(i in 1:2) A[,i]<-fastAnc(tree,X[,i]) } # control list con=list(col.edge=matrix(data="black",nrow(tree$edge),1,dimnames=list(as.character(tree$edge[,2]),"color")),col.node=matrix(data="black",nrow(tree$edge)+1,1,dimnames=list(as.character(c(tree$edge[1,1],tree$edge[,2])),"color"))) con[(namc<-names(control))]<-control con$col.edge<-as.matrix(con$col.edge); con$col.node<-as.matrix(con$col.node) # set xlim & ylim if(hasArg(xlim)) xlim<-list(...)$xlim else xlim<-range(c(X[,1],A[,1])) if(hasArg(ylim)) ylim<-list(...)$ylim else ylim<-range(c(X[,2],A[,2])) # set xlab & ylab if(hasArg(xlab)) xlab<-list(...)$xlab else xlab<-colnames(X)[1] if(hasArg(ylab)) ylab<-list(...)$ylab else ylab<-colnames(X)[2] # set font size for tip labels if(hasArg(fsize)) fsize<-0.75*list(...)$fsize else fsize<-0.75 # plot root state plot(x=A[1,1],y=A[1,2],xlim=xlim,ylim=ylim,xlab=xlab,ylab=ylab) # put X in tree order X<-X[tree$tip.label,] rownames(X)<-1:length(tree$tip) # check for node labels; create if necessary if(is.null(tree$node.label)) tree$node.label<-as.character(length(tree$tip.label)+1:tree$Nnode) tree$node.label<-as.character(tree$node.label) A<-A[tree$node.label,] rownames(A)<-length(tree$tip)+1:tree$Nnode # bind X & A X<-rbind(X,A) # loop across the branches of the tree for(i in 1:nrow(tree$edge)){ lines(x=c(X[as.character(tree$edge[i,1]),1],X[as.character(tree$edge[i,2]),1]),y=c(X[as.character(tree$edge[i,1]),2],X[as.character(tree$edge[i,2]),2]),col=con$col.edge[as.character(tree$edge[i,2]),1]) points(x=X[as.character(tree$edge[i,1]),1],y=X[as.character(tree$edge[i,1]),2],col=con$col.node[as.character(tree$edge[i,1]),1],pch=16,cex=1.0) if(tree$edge[i,2]<=length(tree$tip.label)&&label) textxy(X=X[as.character(tree$edge[i,2]),1],Y=X[as.character(tree$edge[i,2]),2],labs=tree$tip.label[tree$edge[i,2]],cx=fsize) } # plot larger points for the tips points(X[1:length(tree$tip),1],X[1:length(tree$tip),2],col=con$col.node[rownames(X)[1:length(tree$tip)],1],pch=16,cex=1.3) }