repeated measures with limma+lme
0
0
Entering edit mode
@jacob-michaelson-1079
Last seen 10.2 years ago
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
limma limma • 1.4k views
ADD COMMENT

Login before adding your answer.

Traffic: 954 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6