[BioC] A question about paired data and how to set up the designmatrix in LIMMA

[Gordon Smyth]

>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