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.