################################################################################
## Copyright (C) 2012 Peter Murakami <pmurakam@jhsph.edu>
## 
## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
## 
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
## 
## You should have received a copy of the GNU General Public License
## along with this program.  If not, see <http://www.gnu.org/licenses/>.
################################################################################

## Function to plot boxplots separately by a second factor, using different colors
boxplot2 <- function(y, x, by, ylab=NULL, col=NULL, ...) {
    stopifnot(length(y)==length(x))
    stopifnot(length(y)==length(by))
    stopifnot(is.numeric(y))
    un = unique(by)
    n = length(un)
    if(is.null(col)) {
        if(n<=12) {
            require(RColorBrewer)
            palette(brewer.pal(12,"Set3"))
            col = 1:n
        } else {
            stop("You must specify col with as many colors as there are unique values of by.")
        }
    } else {
        if(length(col)<n) {
            message("Too few colors provided.  Re-setting colors.")      
            require(RColorBrewer)
            palette(brewer.pal(12,"Set3"))
            col = 1:n
        }
    }
    m = length(unique(x))
    inc = seq(-0.2,0.2,length.out=n)

    x1  = x[which(by==un[1])]
    y1  = y[which(by==un[1])]
    by1 = by[which(by==un[1])]
    bx = boxplot(y1~x1, ylab=ylab, col=col[1], ylim=range(y), 
                 xlim=c(0.5,m+0.5), boxwex = 0.5/n, at=1:m + inc[1], ...)
    tb = table(x1)
    axis(side=1,at=1:m + inc[1],labels=paste("(",tb[bx$names],")",sep=""),padj=2)

    for(k in 2:length(un)) {
        xk  = x[which(by==un[k])]
        yk  = y[which(by==un[k])]
        byk = by[which(by==un[k])]        
        bx = boxplot(yk~xk, add = TRUE, boxwex = 0.5/n, at = 1:m + inc[k], col = col[k])
        tb = table(xk)
        axis(side=1,at=1:m + inc[k],labels=paste("(",tb[bx$names],")",sep=""),padj=2)
    }
    legend("topright", legend=un, fill = col[1:n])
}