################################################################################
## Copyright (C) 2010 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/>.
################################################################################
    
CI <- function(x,alpha=.05){ #CI using the t-distribution, regardless of n
    m <- mean(x,na.rm=TRUE)
    s <- sd(x,na.rm=TRUE)
    n <- length(which(!is.na(x)))
    arm <- qt(1-alpha/2,df=n-1)*s/sqrt(n)
    return(c(m-arm,m+arm))
}
xtgraph <- function(y,id,time,group,off=(max(time)-min(time))/100,leg=TRUE,
                    where="topright",bar.lty=1,xlab,ylab,...){
    sp <- split(data.frame(cbind(y,id,time)),group)
    lsp <- length(sp)
    pts <- list()
    cis <- list()
    ad <- (off*2)/(lsp-1)
    shift <- -off + ad*seq(0,lsp -1)
    for(i in 1:lsp){
        sh <- shift[i]
        spi <- sp[[i]]
        a <- tapply(spi$y,spi$time,function(x){mean(x,na.rm=TRUE)})
        pts[[i]] <- cbind(a,as.numeric(names(a))  +sh)
        b <- tapply(spi$y,spi$time,CI)
        b <- do.call("rbind",b)
        cis[[i]] <- cbind(b,as.numeric(rownames(b))+sh)
    }
    
    allp0 <- do.call("rbind",cis)
    allp <- c(allp0[,1],allp0[,2])
    plot(min(time)-10, min(allp)-10, xlab=xlab, ylab=ylab,    #blank plot
         xlim=c(min(time),max(time)), ylim=c(min(allp),max(allp)),...) 

    shades <- (.7/(lsp-1))*seq(0,(lsp-1))
    if(lsp<10){ pchs<-c(16,17,15,20:25) }else{pchs<- c(1:lsp)} 
        #to use point types 16,17,15,20:25 if possible    
    for(j in 1:lsp){
        points(pts[[j]][,1]~pts[[j]][,2], col=gray(shades[j]), pch=pchs[j]) 
        lines( pts[[j]][,1]~pts[[j]][,2], col=gray(shades[j]), lty=j)
        
        cisj <- cis[[j]]
        for(k in 1:nrow(cisj)){
            segments(cisj[k,3],cisj[k,1],cisj[k,3],cisj[k,2],
                     lty=bar.lty,col=gray(shades[j]))
        }                 
    }
    if(leg==TRUE){
        legend(where,lty=c(1:lsp),col=gray(shades),cex=.8,pch=pchs[1:lsp],
               legend=names(sp), bty="n")
    }
}