R/fsmry.graph.R

In xkzhou/BTKR: Biostatistics Tools Kit in R

#' Summarize y by x graphically
#'
#' @param y A continuous variable.
#' @param x A continuous or categorical variable.
#' @param type a text string for the type of graph, e.g., "scatter","bxp","errbar","bar",or "bar2"
#' @param y.plab y-axis label
#' @param x.plab x-axis label
#' @param ynames names for different categories of y
#' @param xnames names for different categories of x
#' @param stat.txt statistics to show on the graph, typically, p-value and/or correlation coefficient
#' @param loc.stat location of the stat.txt
#' @param ylim specify the y-axis limit if needed
#' @param fname a text string of file name if a win-meta file is to be generated
#' @param y.fnlab a text string of y variable name to be included in a file name
#' @param x.fnlab a text string of x variable name to be included in a file name
#' @param subgrp a text string of analysis cohort name to be included in a file name
#' @param withline a logic indicator of whether to include the Tukey line in the scatter plot, default true
#' @param geomean a logic indicator of whether to plot geometric mean/sd when plotting the errbar plot, default false
#' @param silent a logic indicator of whether to output the mp variable generated by barplot function, default true
#' @param width width of the graph is a win-meta file is to be generated
#' @param height height of the graph is a win-meta file is to be generated
#' @param mar graphical parameter
#' @param mgp graphical parameter
#' @param cex.bar graphical parameter for specifying text size of text above each bar
#' @param cex graphical parameter
#' @param cex.axis graphical parameter
#' @param ... additional graphical paramter
#' @return Summary graph of y by x.
#' @examples
#' set.seed(16)
#' dat.work <- data.frame(ht = c(rnorm(10, mean=1.72, sd=0.1), rnorm(10, mean=1.65, sd=0.1)),
#'                        wt = c(rnorm(10, mean=70, sd=10), rnorm(10, mean=60, sd=10)),
#'                        sex = factor(rep(c("Female", "Male", "Female", "Male"), c(2,8,6,4))),
#'                        group = factor(rep(c("grp1", "grp2"), each=10)))
#' cor.test(dat.work$wt, x = dat.work$ht)
#' fsmry.graph(y = dat.work$wt, x = dat.work$ht, type ="scatter",
#'             y.plab = "Weight",
#'             x.plab = "Height",
#'             stat.txt = "r=-0.24, p=0.31")
#' fsmry.by.grp(y = dat.work$ht, grp = dat.work$sex)
#' fsmry.graph(y = dat.work$ht, x = dat.work$sex, type ="bxp",
#'             y.plab = "Weight",
#'             x.plab = "Sex",
#'             stat.txt = "p=0.68")
#' fsmry.graph(y = dat.work$ht, x = dat.work$sex, type ="errbar",
#'             y.plab = "Weight",
#'             x.plab = "Sex",
#'             stat.txt = "p=0.66")
#' fsmry.graph(y = dat.work$ht, x = dat.work$sex, type ="errbar",
#'             y.plab = "Weight",
#'             stat.txt = "p=0.68",
#'             geomean=T)
#' fsmry2.by.grp(y = dat.work$sex, grp=dat.work$group, cmp.method="fisher")
#' fsmry.graph(y = dat.work$sex, x = dat.work$group, type ="bar",
#'             stat.txt = "p=0.17")
#' fsmry.graph(y = dat.work$sex, x = dat.work$group, type ="bar",
#'             xnames=c("Group1", "Group2"),
#'             stat.txt = "p=0.17")
#' fsmry.graph(y = dat.work$sex, x = dat.work$group, type ="bar2",
#'             xnames=c("Group1", "Group2"),
#'             stat.txt = "p=0.17")
fsmry.graph <- function(y, x,
                        type=c("scatter","bxp","errbar","bar","bar2"),
                        y.plab="",
                        x.plab="",
                        ynames=NULL,
                        xnames=NULL,
                        stat.txt=NULL,
                        loc.stat=NULL,
                        ylim=NULL,
                        xlim=NULL,
                        fname=NULL,
                        y.fnlab=NULL,
                        x.fnlab=NULL,
                        subgrp="",
                        withline=T,
                        geomean=F,
                        silent=T,
                        width=3, height=3,
                        mar=NULL, mgp=NULL,
                        cex.bar=NULL,
                        cex=1,
                        cex.axis=1,
                        ...){
  ## This function is used to generate commonly used plots
  ##browser()
  if(!is.null(fname)){
    if(fname=="sys"){
      fdate <- gsub("-", "", Sys.Date())
      file.name <- paste("Fig/fig_",type, "_", gsub("\\.","_",y.fnlab),
                         "_",gsub("\\.","_",x.fnlab),"_",subgrp,"_",
                         fdate,".wmf", sep="")
    }
    else
      file.name <- fname
    win.metafile(file=file.name,width=width, height=height)
  }
  if(is.null(mar)) {
    if(class(x.plab)=="character"){
      if(x.plab=="")
        mar <- c(2.5,2.5,0.5,0.5)
    }
    else
      mar <- c(3.5,3.5,0.5,0.5)
  }
  if(is.null(mgp)) mgp <- c(2.2,0.8,0)
  par(mar=mar, mgp=mgp)

  if(type=="scatter"){
    if(is.null(ylim)) y.lim <- range(pretty(y), na.rm=T)
    else y.lim <- ylim
    if(is.null(xlim)) x.lim <- range(pretty(x), na.rm=T)
    else x.lim <- xlim
    plot(x, y, xlab=x.plab, ylab=y.plab,pch=20, cex=cex, las=1, xlim=x.lim, ylim=y.lim,...)
    if(withline)
      abline(coef(line(x,y)),...)
    if(!is.null(stat.txt)){
      if(is.null(loc.stat))
        text(xlim[1], ylim[2],label=stat.txt, cex=cex, pos=4)
      else
        text(loc.stat[1], loc.stat[2],label=stat.txt, cex=cex)
    }
  }
  if(type=="bxp"){
    if(is.null(ylim))
      y.lim <- range(pretty(c(floor(min(y, na.rm=T)),
                              ceiling(max(y, na.rm=T)))))
    else y.lim <- ylim
    if(is.null(xnames)) xnames <- levels(x)
    boxplot(y~x, boxwex=0.5, xlab=x.plab, ylab=y.plab, names=xnames, ylim=y.lim,
            las=1, cex.axis=cex.axis,...)
    if(!is.null(stat.txt)){
      if(is.null(loc.stat))
        text(x=0.9, y=y.lim[2]-0.01*(y.lim[2]-y.lim[1]),label=stat.txt, cex=cex)
      else
        text(loc.stat[1], loc.stat[2],label=stat.txt, cex=cex)
    }
  }
  if(type=="errbar"){
    ##library(Hmisc)
    if(geomean){
      y.mean <- tapply(y, x,function(z) exp(mean(log(z), na.rm=T)))
      y.sd <- tapply(y,x,function(z) exp(sd(log(z), na.rm=T)))
    }
    else{
      y.mean <- tapply(y, x,function(z) mean(z, na.rm=T))
      y.sd <- tapply(y,x,function(z) sd(z, na.rm=T))
    }

    y.lim <- range(0, ceiling(y.mean+1.5*y.sd))
    if(is.null(xnames)) xnames <- levels(x)
    if(length(y.mean)==2)
      mp <- barplot(y.mean,
                    space=c(0.3,rep(0.3, length(xnames)-1))/0.6,
                    width=0.6,
                    xlim=c(0,length(xnames)),col="black", names="",
                    ylim=y.lim, xlab=x.plab, ylab=y.plab, las=1,...)
    else
      mp <- barplot(y.mean,
                    space=c(0.2,rep(0.4, length(xnames)-1))/0.6,
                    width=0.6,
                    xlim=c(0,length(xnames)),col="black", names="",
                    ylim=y.lim, xlab=x.plab, ylab=y.plab, las=1, ...)
    segments(x0=mp, x1=mp, y0=y.mean, y1=y.mean+y.sd,...)
    segments(x0=mp-0.025, x1=mp+0.025, y0=y.mean+y.sd, y1=y.mean+y.sd,...)
    ##errbar(x=mp, y=y.mean, yplus=y.mean+y.sd, yminus=y.mean-y.sd,
    ##       pch=".", add=T)
    axis(1, at=mp, labels=xnames, tick=0,...)
    axis(2, las=1,...)
    box(...)
    if(is.null(loc.stat))
      text(x=mp[1]-0.1, y=y.lim[2]-0.05*y.lim[2],label=stat.txt, cex=cex)
    else
      text(loc.stat[1], loc.stat[2],label=stat.txt, cex=cex)
    if(!silent) return(mp)
  }
  if(type=="bar"){
    ##browser()
    y.by.x <- table(y,x)
    y.prop <- round(prop.table(y.by.x, margin=2)[2,]*100)
    y.txt <- paste(y.by.x[2,], "/", apply(y.by.x,2, sum), sep="")
    y.lim <- range(0, 105)

    if(is.null(xnames)) xnames <- levels(x)
    if(is.null(ynames)) ynames <- levels(y)
    if(class(y.plab)=="character")
      if(y.plab=="") y.plab <- paste0("% of ", ynames[2], " cases")

    if(length(y.prop)==2)
      mp <- barplot(y.prop, space=c(0.3,rep(0.3, length(xnames)-1))/0.6,
                    width=0.6,
                    xlim=c(0,length(xnames)),col="black", names="",
                    ylim=y.lim, xlab=x.plab, ylab=y.plab, las=1,...)
    else
      mp <- barplot(y.prop, space=c(0.2,rep(0.4, length(xnames)-1))/0.6,
                    width=0.6,
                    xlim=c(0,length(xnames)),col="black", names="",
                    ylim=y.lim, xlab="", ylab=y.plab, las=1,...)
    axis(1, at=mp, labels=xnames,tick=0,...)
    text(x=mp, y=y.prop+3, y.txt, cex=ifelse(is.null(cex.bar),0.7, cex.bar), font=2)
    if(is.null(loc.stat))
      text(x=mp[1]-0.2*diff(mp)[1], y=98, label=stat.txt, cex=cex, pos=4)
    else
      text(loc.stat[1], loc.stat[2],label=stat.txt, cex=cex)
    axis(2, las=1,...)
    box(...)
    if(!silent) return(mp)
  }
  if(type == "bar2"){
    ##browser()
    y.by.x <- table(y,x)
    y.prop <- round(prop.table(y.by.x, margin = 2)*100)
    y.txt <- paste0(y.by.x, "/", rep(apply(y.by.x, 2, sum),
                                     each = length(levels(y))))
    y.lim <- range(0, 105)
    if(is.null(xnames)) xnames <- levels(x)
    if(class(y.plab)=="character")
      if(y.plab=="") y.plab <- "Percent"
    mp <- barplot(y.prop, width=0.6,
                  col= cm.colors(length(levels(y))), names=levels(x),
                  ylim=y.lim, xlab=x.plab, ylab=y.plab, las=1, beside = TRUE,
                  legend = levels(y),...)
    text(x=mp, y=y.prop+3, y.txt,
         cex=ifelse(is.null(cex.bar),0.7, cex.bar), font=2)
    if(is.null(loc.stat))
      text(x=mp[1]-0.2*diff(mp)[1], y=98, label=stat.txt, cex=cex, pos=4)
    else
      text(loc.stat[1], loc.stat[2],label=stat.txt, cex=cex)
    axis(2, las=1,...)
    box(...)
    if(!silent) return(mp)
  }
  if(!is.null(fname))  dev.off()
}