[R] plot graph with error bars trouble

[Jim Lemon]

Marcelo Laia wrote:
> Hi,
> 
> I have a data set like this:
> 
> Mutant    Rep    Time   OD
> 02H02    1    0    0.029
> 02H02    2    0    0.029
> 02H02    3    0    0.023
> 02H02    1    8    0.655
> 02H02    2    8    0.615
> 02H02    3    8    0.557
> 02H02    1    12    1.776
> 02H02    2    12     1.859
> 02H02    3    12    1.668
> 02H02    1    16    3.379
> 02H02    2    16    3.726
> 02H02    3    16    3.367
> 306    1    0    0.033
> 306    2    0    0.035
> 306    3    0    0.034
> 306    1    8     0.377
> 306    2    8    0.488
> 306    3    8    0.409
> 306    1    12    1.106
> 306    2    12    1.348
> 306    3    12    1.246
> 306    1    16    2.706
> 306    2    16    3.073
> 306    3    16    3.038
> 
> I need to plot a graph OD over the time for each one mutant with error bars.
> 
> I try the package sciplot, but this package is set up to handle
> factorial treatments, so the spacing in x-axis is fixed to be equal.
> Than, with it I got something like this:
> 
> |
> |
> |
> |
> |
> +-----------------------------
> 0        8          12          16
> 
> But, I would like spacing between 0 and 8 2-fold the spacign between 8
> and 12, like this:
> 
> |
> |
> |
> |
> |
> +--------------------------------------
> 0        4          8          12          16
> 
> Could you point me out another way to do this with out using sciplot?
> Any suggestion is very appreciated.
> 
Hi Marcelo,
This might be a task for brkdn.plot in the plotrix package except for 
the X positions. I have written a revision of the function as this will 
be a useful enhancement and minimally tested it, so if you would like to 
give it a try, here it is:

library(plotrix)
# here is the revised function, may have some line breaks
brkdn.plot<-function(vars,groups=NA,obs=NA,data,
  mct="mean",md="std.error",
  stagger=NA,dispbar=TRUE,
  main="Breakdown plot",xlab=NA,ylab=NA,xaxlab=NA,
  ylim=NA,type="b",pch=1,lty=1,col=par("fg"),...) {

  if(is.na(obs)) {
   if(is.na(groups))
    stop("Must have at least one factor to subset data")
   bygroup<-as.factor(data[[groups]])
   grouplevels<-levels(bygroup)
   ngroups<-length(grouplevels)
   nobs<-length(vars)
   obslevels<-1:nobs
  }
  else {
   byobs<-as.factor(data[[obs]])
   obslevels<-levels(byobs)
   nobs<-length(obslevels)
   if(is.numeric(unlist(data[obs])))
    obs.pos<-unique(unlist(data[obs]))
   else obs.pos<-1:nobs
   if(is.na(groups)) {
    ngroups<-length(vars)
    grouplevels<-1:ngroups
   }
   else {
    bygroup<-as.factor(data[[groups]])
    grouplevels<-levels(bygroup)
    ngroups<-length(grouplevels)
    if(length(vars) > 1) {
     warning("Group and observation present, only vars[1] plotted")
     vars<-vars[1]
    }
   }
  }
  brkdn<-list(matrix(NA,nrow=ngroups,ncol=nobs),
   matrix(NA,nrow=ngroups,ncol=nobs))
  if(is.na(groups)) {
   if(is.na(xlab)) xlab<-"Observation"
   xat<-1:nobs
   if(is.na(xaxlab[1])) xaxlab<-obslevels
   for(group in 1:ngroups) {
    for(ob in 1:nobs) {
     brkdn[[1]][group,ob]<-
      do.call(mct,list(unlist(subset(data[[vars[group]]],
       data[[obs]] == obslevels[ob],vars[group])),na.rm=TRUE))
     if(!is.na(md))
      brkdn[[2]][group,ob]<-
       do.call(md,list(unlist(subset(data[[vars[group]]],
        data[[obs]] == obslevels[ob],vars[group])),na.rm=TRUE))
    }
   }
  }
  else {
   if(is.na(obs)) {
    if(is.na(xlab)) xlab<-"Variable"
    xat<-1:length(vars)
    if(is.na(xaxlab[1])) xaxlab<-vars
    for(group in 1:ngroups) {
     for(ob in 1:nobs) {
      brkdn[[1]][group,ob]<-
       do.call(mct,list(unlist(subset(data[[vars[ob]]],
        data[[groups]] == grouplevels[group],vars[ob])),na.rm=TRUE))
      if(!is.na(md))
       brkdn[[2]][group,ob]<-
        do.call(md,list(unlist(subset(data[[vars[ob]]],
         data[[groups]] == grouplevels[group],vars[ob])),na.rm=TRUE))
     }
    }
   }
   else {
    if(is.na(xlab)) xlab<-"Observation"
    xat<-obs.pos
    if(is.na(xaxlab[1])) xaxlab<-obslevels
    for(group in 1:ngroups) {
     for(ob in 1:nobs) {
      brkdn[[1]][group,ob]<-
       do.call(mct,list(unlist(subset(data,data[[groups]] == 
grouplevels[group] &
        data[[obs]] == obslevels[ob],vars)),na.rm=TRUE))
      if(!is.na(md))
       brkdn[[2]][group,ob]<-
        do.call(md,list(unlist(subset(data,data[[groups]] == 
grouplevels[group] &
         data[[obs]] == obslevels[ob],vars)),na.rm=TRUE))
     }
    }
   }
  }
  if(is.na(ylim[1])) {
   if(is.na(md)) ylim<-range(brkdn[[1]])
   else
    ylim<-c(min(brkdn[[1]]-brkdn[[2]],na.rm=TRUE),
     max(brkdn[[1]]+brkdn[[2]],na.rm=TRUE))
  }
  groupdiv<-ifelse(ngroups < 3,1,ngroups-2)
  if(is.na(stagger)) stagger<-0.025-groupdiv*0.0025
  if(is.na(ylab)) {
   if(length(vars) == 1) ylab<-vars[1]
   else ylab<-paste(vars,collapse=" and ")
  }
  plot(0,xlim=c(obs.pos[1]-0.5,obs.pos[nobs]+0.5),
   main=main,xlab=xlab,ylab=ylab,ylim=ylim,
   type="n",axes=FALSE,...)
  box()
  axis(1,at=xat,labels=xaxlab)
  axis(2)
  if(length(pch) < ngroups) pch<-rep(pch,length.out=ngroups)
  if(length(col) < ngroups) col<-rep(col,length.out=ngroups)
  if(length(lty) < ngroups) lty<-rep(lty,length.out=ngroups)
  offinc<-stagger*diff(par("usr")[c(1,2)])
  offset<-0
  arrow.cap<-0.01-(groupdiv*0.001)
  for(group in 1:ngroups) {
   points(obs.pos+offset,brkdn[[1]][group,],type=type,
    col=col[group],pch=pch[group],lty=lty[group])
   if(dispbar)
    dispbars(obs.pos+offset,brkdn[[1]][group,],brkdn[[2]][group,],
     arrow.cap=arrow.cap,col=col[group])
   offset<-ifelse(offset<0,-offset,-offset-offinc)
  }
  names(brkdn)<-c(mct,md)
  return(brkdn)
}
# read in the data
ml.df<-read.csv("../ml.csv")
# call the new function
brkdn.plot("OD","Mutant","Time",ml.df,col=c(2,3))

Jim