[BioC] repeated measures with limma+lme

[jjmichael@comcast.net]

Hi all,

I was wondering if I could get some constructive feedback on some code I'm using to analyze an Affymetrix repeated measures experiment.  I've found through the list that limma's repeated measures capabilities are limited when it comes to specifying certain correlation structures, but I really like the way it handles contrasts, so I tried to find a way to fit the model using an lme loop, and then insert the proper components into the limma fit.  

I fit the same model using both lme and lmfit, and found that the correlation structure used in lme only affected the sigma and residual df components of the model (all others were identical), so my logic here is to insert the sigma and df.residual components from the lme model into the lmfit model, thus enabling me to use limma's convenient methods for contrasts.

In this experiment, there are three treatment levels which were observed in 6 subjects (two reps per treatment) over 7 time points. To specify the model, I used the method in the limma user's guide that has all factor level combinations (treatment*time, for example) as individual levels of a single factor.

Here's the code:

design=model.matrix(~0+int)###"int" is a factor whose levels are the factor level                 ###combinations
colnames(design) = levels(int)

fit=lmFit(exprset, design)

###create the empty components to store the values
sigma=numeric(length=nrow(exprs(exprset)))
df.residual=numeric(length=nrow(exprs(exprset)))

dataset=exprs(exprset)

###specify correlation structure
correlation.structure = corAR1(form= ~1|rep)

###lme loop
for(i in 1:nrow(exprs(exprset))){
 	datai=dataset[i,]
 	out.lme=lme(datai~int+0,random=~1|rep)
 	out1.lme=update(out.lme, correlation=correlation.structure)
 	sigma[i]=(out1.lme)$sigma
	df.residual[i]=(out1.lme)$fixDF$term
 }
###replace the desired components into the lmFit model
fit$sigma = sigma
fit$df.residual = df.residual


This seems to be okay, according to my very limited statistical knowledge, but I'd like to hear back from the experts to see if there are any glaring omissions.

Thanks in advance,

Jake