Monday, July 25, 2016

New function to print box plot next to a tree

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

plot of chunk unnamed-chunk-2

## 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)"))

plot of chunk unnamed-chunk-2

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