R/lattice xyplot with error bounds

errorBounds.R

#
# From Hmisc, pasted here so as to not need any deps for this demo
#
Cbind <- function (...) 
{
  dotlist <- list(...)
  if (is.matrix(dotlist[[1]])) {
    y <- dotlist[[1]]
    ynam <- dimnames(y)[[2]]
    if (!length(ynam)) 
      stop("when first argument is a matrix it must have column dimnames")
    other <- y[, -1, drop = FALSE]
    return(structure(y[, 1], class = "Cbind", label = ynam[1], 
                     other = other))
  }
  lname <- names(dotlist)
  name <- vname <- as.character(sys.call())[-1]
  for (i in 1:length(dotlist)) {
    vname[i] <- if (length(lname)) 
      lname[i]
    else ""
    if (vname[i] == "") 
      vname[i] <- name[i]
  }
  lab <- attr(y <- dotlist[[1]], "label")
  if (!length(lab)) 
    lab <- vname[1]
  if (!is.matrix(other <- dotlist[[2]]) || ncol(other) < 2) {
    other <- as.matrix(as.data.frame(dotlist))[, -1, drop = FALSE]
    dimnames(other)[[2]] <- vname[-1]
  }
  structure(y, class = "Cbind", label = lab, other = other)
}
  
#
# My function which is misbehaving
#

errorBounds <- function(x, y, ..., cols=trellis.par.get("superpose.line")$col) {
  xyplot(x, y, ..., cols=cols,
         panel = function(x, y, ..., groups=NULL, subscripts=NULL) {
           groups <- groups[subscripts]
           if (!is.null(groups)) {
             for (g in unique(groups)) {
               idx <- which(groups==g)
               col = rgb(t(col2rgb(cols[(which(levels(as.factor(as.character(groups)))==g)-1)%%length(cols)+1])/255))
               col = paste(col,"40",sep="")
               panel.polygon(
                 c(x[idx], rev(x[idx])),
                 c(attr(y,"other")[idx,1], rev(attr(y,"other")[idx,2])),
                 col = col
               )
             }
           } else {
             col = rgb(t(col2rgb(cols[1])/255))
             col = paste(col,"40",sep="")
             panel.polygon(
               c(x, rev(x)),
               c(attr(y,"other")[,1], rev(attr(y,"other")[,2])),
               col = col
             )
           }
           panel.xyplot(x, y, ..., groups=groups, subscripts=subscripts)
         }
  )
}

#
# Test
#

n <- 50
X <- 1:n
d <- data.frame(
  x = rep(X,4),
  y = c(
    sin(2*pi*(X / n)) + rnorm(n, mean = 0, sd = 0.1),
    cos(2*pi*(X / n)) + rnorm(n, mean = 0, sd = 0.1),
    sin(2*pi*(X / n)) + rnorm(n, mean = 0, sd = 0.1),
    cos(2*pi*(X / n)) + rnorm(n, mean = 0, sd = 0.1)
  ),
  group = rep(c("A","B"),each=n*2),
  type = rep(rep(c("sin","cos"),each=n),2)
)
d$upper = d$y + 0.5 + rnorm(nrow(d), mean = 0, sd = 0.01)
d$lower = d$y - 0.5 - rnorm(nrow(d), mean = 0, sd = 0.01)

p2 = errorBounds(
  Cbind(y,lower,upper)~x|group,d,groups=type,
  layout=c(1,2), as.table=T,
  type="l",label.curves=F,
  scales=list(y=list(relation="free"))
  )

print(p2)