To complement my recent function that uses a plot array to print a bar plot next to a plotted tree (1, 2, 3), I thought I'd also show how to do the same with a plotted box plot.
Here is the function. It modifies the plotTree.barplot version in
a couple of interesting ways:
plotTree.boxplot<-function(tree,x,args.plotTree=list(),
args.boxplot=list()){
cw<-reorder(tree)
if(!is.list(x)&&class(x)!="formula"){
obj<-setNames(
lapply(cw$tip.label,function(x,y) y[which(names(y)==x)],
y=x),cw$tip.label)
} else obj<-x
if(class(x)=="formula")
args.boxplot$formula<-obj else args.boxplot$x<-obj
args.boxplot$horizontal<-TRUE
args.boxplot$axes<-FALSE
args.boxplot$names.arg<-""
args.boxplot$xlim<-c(1,Ntip(cw))
if(is.null(args.boxplot$space)) args.boxplot$space<-0.7
if(is.null(args.boxplot$mar))
args.boxplot$mar<-c(5.1,0,2.1,1.1)
else args.boxplot$mar[2]<-0.1
args.plotTree$tree<-cw
if(is.null(args.plotTree$mar))
args.plotTree$mar<-c(5.1,1.1,2.1,0)
else {
args.plotTree$mar[4]<-0
}
if(args.plotTree$mar[1]!=args.boxplot$mar[1])
args.plotTree$mar[1]<-args.boxplot$mar[1]
if(args.plotTree$mar[3]!=args.boxplot$mar[3])
args.plotTree$mar[3]<-args.boxplot$mar[3]
if(is.null(args.plotTree$ftype)) args.plotTree$ftype<-"i"
if(is.null(args.plotTree$lwd)) args.plotTree$lwd<-1
par(mfrow=c(1,2))
do.call(plotTree,args.plotTree)
par(mar=args.boxplot$mar)
ii<-which(names(args.boxplot)%in%c("formula","x"))
args.boxplot<-c(args.boxplot[ii],args.boxplot[-ii])
obj<-do.call(boxplot,args.boxplot)
axis(1)
if(!is.null(args.boxplot$xlab)) title(xlab=args.boxplot$xlab)
else title(xlab="x")
invisible(obj)
}
Obviously, the arguments passed to the plotting method are a little bit
different - but most significantly, boxplot can take as its first
argument either x, generally a list in which each element is
a vector of values corresponding to the rows of the boxplot, or a formula.
Let's try it:
library(phytools)
## here's our tree & data
tree
##
## Phylogenetic tree with 26 tips and 25 internal nodes.
##
## Tip labels:
## Z, Y, X, W, V, U, ...
##
## Rooted; includes branch lengths.
xe[1:30] ## first 30 values
## Z Z Z Z Z Y
## -0.21795787 -0.32369195 0.61508766 -0.71780324 -0.09337898 0.29571288
## Y Y Y Y X X
## 0.63917603 -0.06377179 -0.19115456 0.20776738 -1.24725366 -0.23757952
## X X X W W W
## -0.99325267 -0.07823182 -0.63543671 -0.43174636 -0.95828938 -0.62417070
## W W V V V V
## 0.46085304 0.31871249 -1.83209464 -0.77540702 -2.05616575 -1.86725816
## V U U U U U
## -1.00503367 -2.28712937 -2.11539566 -1.30714430 -0.71566793 -1.00299982
plotTree.boxplot(tree,xe) ## default
## using formula:
spp<-factor(names(xe),
levels=reorder(tree)$tip.label)
spp
## [1] Z Z Z Z Z Y Y Y Y Y X X X X X W W W W W V V V V V U U U U U T T T T T
## [36] S S S S S R R R R R Q Q Q Q Q P P P P P O O O O O N N N N N M M M M M
## [71] L L L L L K K K K K J J J J J I I I I I H H H H H G G G G G F F F F F
## [106] E E E E E D D D D D C C C C C B B B B B A A A A A
## Levels: Z Y X W V U T S R Q P O N M L K J I H G F E D C B A
plotTree.boxplot(tree,xe~spp,
args.boxplot=list(mar=c(5.1,0,4.1,1.1),
main="snout length in rhinogrades",
xlab="relative length (inches)"))
One trick is that for the formula method we need to make sure
that the levels of our factor (in this case, spp) are in the
order of the tip labels of our cladewise re-ordered tree!
Nonetheless, hope this is a handy function for some users.
Note that this function basically encapsulates a method that I described in an earlier post, with some additional flexibility.
No comments:
Post a Comment
Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.