Entering edit mode

>Date: Sun, 17 Jul 2005 12:05:41 +0200
>From: "Johan Lindberg" <johanl at="" biotech.kth.se="">
>Subject: Re: [BioC] A question about paired data and how to set up
the
> designmatrix in LIMMA
>To: <bioconductor at="" stat.math.ethz.ch="">
>
>Dear Bioconductor users.
>I posted a question a while ago that did not get an answer. I'm sorry
to
>have to nag and ask the same question again but I really need an
answer.
>The question that I'm asking is how to do the paired analysis in
LIMMA?
>If you don?t have the time to look at my previous post here is a
shorter
>version.
>
>I have tried to read previous posts on this matter but I found the
>answers kind of vague (perhaps I missed something due to my limited
>knowledge in statistics, in that case please guide me to an
explaining
>post and I apologize for asking the same question again).
>
>I have affydata on 4 patients (paired) before and after treatment
>(unbalanced due to different number of biopsies from each patient). I
>want to fit the patients as a fixed effect. I am interested in the
>effect due to treatment.
>
>In previous posts I have seen Gordon give directions about two ways
to
>fit patients as a fixed effect with LIMMA.
>
>http://files.protsuggest.org/biocond/html/8297.html
>
>and
>
>http://files.protsuggest.org/biocond/html/8175.html
>
>My situation look like this:
>#########################################################
>#Fitting the patients as a fixed effect, I am not totally shure of my
>#code, please correct me if I am wrong
>eset<-rma(Data)
>tmp<-pData(eset)
>tmp<-cbind(tmp,Before=c(1,0,1,1,1,0,1,0,1,0,1,0,1,0))
>tmp<-cbind(tmp,After=c(0,1,0,0,0,1,0,1,0,1,0,1,0,1))
>tmp<-cbind(tmp,Patient=c(1,1,1,1,1,1,2,2,3,3,3,3,4,4))
>pData(eset)<-tmp
>design <- model.matrix(~(After - Before) + Patient,
>data=as.data.frame(tmp))
>fit<-lmFit(eset,design)
>fitModel <- eBayes(fit)
This code is incorrect for a couple of reasons. Firstly, you've fitted
Patient as a numerical covariate rather than as a factor. This of
course is
meaningless.
Secondly, as Robert Gentleman has pointed out, your experiment is not
a
paired t-test set up, so the fixed effect approach to Patients is not
correct. You need to fit a mixed model with Patient as a random
effect. You
are nearly doing this correctly below.
>#########################################################
>#Using the block argument to ?block? patient.
>eset<-rma(Data)
>tmp<-pData(eset)
>tmp<-cbind(tmp,Before=c(1,0,1,1,1,0,1,0,1,0,1,0,1,0))
>tmp<-cbind(tmp,After=c(0,1,0,0,0,1,0,1,0,1,0,1,0,1))
>pData(eset) <- tmp
>design <- model.matrix(~(After - Before), data=pData(eset))
>cor<-duplicateCorrelation(eset,
>design,block=c(1,1,1,1,1,1,2,2,3,3,3,3,4,4))
>cor$cor
># [1] 0.2739284
>fit<-lmFit(eset, design, block=c(1,1,1,1,1,1,2,2,3,3,3,3,4,4))
>fitcor <- eBayes(fit)
This code would be correct if you had specified the 'correlation'
argument
to lmFit(). At the moment you're just using the default value for the
correlation, which is incorrect. You need the 2nd coefficient:
topTable(fitcor, coef=2)
>So my questions are, if my code is right, which way is the right way
to
>fit patients as a fixed effect since the result differ (the lod-
scores
>differ between the different models)?
It is not right to fit patient as a fixed effect in this case. You
must fit
patients as a random effect. The second model is correct.
> If my code is wrong, please
>correct it. And last, may I be bald enough to suggest (since there
seem
>to be a lot of post of how to fit fixed effect in LIMMA) that it
would
>be great if there were an example of how to fit fixed effects in the
>LIMMA users guide?
You mean, I think, examples of how to fit random effects. This is in
the
section called "Technical Replication".
Gordon
>Thank you for your time
>
>Best regards
>
>Johan Lindberg