Entering edit mode
Jacob Michaelson
▴
320
@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