Gordon K Smyth | 1 Jul 01:39 2011

Re: sequential removeBatchEffect()?

Hi Jenny,

I don't think I'd recommend the sequential approach, which might give 
non-optimal results if the two batch factors are highly correlated with 
one another.

In fact, I'd be tempted to remove only the pop batch effect, since the 
BeadChip effect is relatively small.

Removing both batch effects requires that you have enough data to estimate 
all your experimental treatments reliably as well as the two batch 
effects.  If you want to do this, here is a function that will remove two 
batch effects at once:

removeBatchEffect <- function(x,batch,batch2=NULL,design=matrix(1,ncol(x),1))
  x <- as.matrix(x)
  batch <- as.factor(batch)
  contrasts(batch) <- contr.sum(levels(batch))
  if(is.null(batch2)) {
   X <- model.matrix(~batch)[,-1,drop=FALSE]
  } else {
   batch2 <- as.factor(batch2)
   contrasts(batch2) <- contr.sum(levels(batch2))
   X <- model.matrix(~batch+batch2)[,-1,drop=FALSE]
  X <- qr.resid(qr(design),X)
  qrX <- qr(X)
(Continue reading)