A question about paired data and how to set up the designmatrix in LIMMA
1
1
Entering edit mode
@gordon-smyth
Last seen 2 hours ago
WEHI, Melbourne, Australia
>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
limma affydata • 1.4k views
ADD COMMENT
0
Entering edit mode
@gordon-smyth
Last seen 2 hours ago
WEHI, Melbourne, Australia
>Date: Sun, 17 Jul 2005 08:36:01 -0700
>From: Robert Gentleman <rgentlem at fhcrc.org>
>Subject: Re: [BioC] A question about paired data and how to set up the designmatrix in LIMMA
>To: Johan Lindberg <johanl at biotech.kth.se>
>Cc: bioconductor at stat.math.ethz.ch
>
>Hi Johan,
>    The simple answer is that your data do not quite fit the paired
>t-test model. You probably want some form of random effects model as can
>be fit by lme (nlme pacakge) or lmer (Matrix package). Where patients
>are treated as a random effect and before/after as a fixed effect. I do
>not believe that limma fits these models, although it may be extended at
>some time.

Actually limma does fit random effect models, if there is only one random factor. It does this through the 'block' argument to duplicateCorrelation() and lmFit(). The actual fitting of the mixed models is done by a call to the randomizedBlock() function in the statmod package. Limma uses a particular approach to Bayes-type moderation of the random effects across genes. See my reply to Johan Lindberg.

Gordon

ADD COMMENT

Login before adding your answer.

Traffic: 669 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