[BioC] 2 way anova in Bioconductor
Jenny Drnevich
drnevich at illinois.edu
Thu Jun 30 15:54:20 CEST 2011
Hi Natasha,
It turns out it depends on whether you set the contr.sum() on your
factor before calling model.matrix(). It's actually in the
limmaUsersGuide() in section 8.7 - you have to look closely at the
treatment-contrast parametrization (2nd example) and the sum to zero
parametrization (3rd example) and the Comparison section of the
tables. Here's an example:
#1st do treatment-contrast parametrization
> Sex <- factor(rep(c("Female","Male"),each=3),levels=c("Male","Female"))
> Sex
[1] Female Female Female Male Male Male
Levels: Male Female
> model.matrix(~Sex)
(Intercept) SexFemale
1 1 1
2 1 1
3 1 1
4 1 0
5 1 0
6 1 0
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$Sex
[1] "contr.treatment"
#So yes, in this design matrix the 2nd coef is Female - Male
#Now switch to the sum to zero parametrization by setting contr.sum()
on the Sex factor:
> contrasts(Sex) <- contr.sum(2)
> Sex
[1] Female Female Female Male Male Male
attr(,"contrasts")
[,1]
Male 1
Female -1
Levels: Male Female
> model.matrix(~Sex)
(Intercept) Sex1
1 1 -1
2 1 -1
3 1 -1
4 1 1
5 1 1
6 1 1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$Sex
[,1]
Male 1
Female -1
# In this design matrix, the 2nd coef is now (Male - Female)/2. I
don't why it works like this, but it messed me up many times before I
figured it out!
Cheers,
Jenny
>
At 04:26 AM 6/30/2011, Natasha Sahgal wrote:
>Hi Jenny and the BioC list,
>
>Please correct me if I am wrong here, but isn't the contrast "Female
>- Male" rather than "Male - Female" as you say?
>
>I was of the understanding that the contrasts, when not specified by
>contrasts.matrix, by default is the latter - former (e.g. B - A), in
>a 2 group comparison (say when levels = c(A,B)).
>
>
>BW,
>Natasha
>
>On 29/06/2011 15:31, Jenny Drnevich wrote:
>>Actually, you can do ANY sort of factorial design in limma. It took
>>me awhile to figure out how to expand the sum to zero
>>parametrization (3rd example, Section 8.7 limmaUsersGuide() ) when
>>you have more than 2 levels of any factor. Here's an example of a 2x3 design:
>>
>>#First set up the 2 group factor, just like in the limma vignette:
>>
>> > Sex <- factor(rep(c("Female","Male"),each=6),levels=c("Male","Female"))
>># Note that the order you specify the groups in the levels argument
>>determines the direction of the comparison. See below.
>>
>> > contrasts(Sex) <- contr.sum(2)
>> > Sex
>> [1] Female Female Female Female Female Female Male Male Male
>>Male Male Male
>>attr(,"contrasts")
>> [,1]
>>Male 1
>>Female -1
>>Levels: Male Female
>># Note that the contrast is male - female because female was listed
>>last in the levels argument above
>>
>>#Now set up the 3 group factor:
>>
>> > Time <- factor(rep(1:3,4),levels=3:1)
>> # Again, the group order specified in the levels
>> determines the direction of comparison
>>
>> > contrasts(Time) <- contr.sum(3)
>> # You use contr.sum(3) here because you have 3 groups. If
>> you have 4 groups, you'd use 4, etc.
>>
>> > Time
>> [1] 1 2 3 1 2 3 1 2 3 1 2 3
>>attr(,"contrasts")
>> [,1] [,2]
>>3 1 0
>>2 0 1
>>1 -1 -1
>>Levels: 3 2 1
>># Note the contrasts are 3 - 1 and 2 - 1, because group 1 was
>>listed last in the levels argument above
>>
>> > design <- model.matrix(~Sex*Time)
>> > design
>> (Intercept) Sex1 Time1 Time2 Sex1:Time1 Sex1:Time2
>>1 1 -1 -1 -1 1 1
>>2 1 -1 0 1 0 -1
>>3 1 -1 1 0 -1 0
>>4 1 -1 -1 -1 1 1
>>5 1 -1 0 1 0 -1
>>6 1 -1 1 0 -1 0
>>7 1 1 -1 -1 -1 -1
>>8 1 1 0 1 0 1
>>9 1 1 1 0 1 0
>>10 1 1 -1 -1 -1 -1
>>11 1 1 0 1 0 1
>>12 1 1 1 0 1 0
>>attr(,"assign")
>>[1] 0 1 2 2 3 3
>>attr(,"contrasts")
>>attr(,"contrasts")$Sex
>> [,1]
>>Male 1
>>Female -1
>>
>>attr(,"contrasts")$Time
>> [,1] [,2]
>>3 1 0
>>2 0 1
>>1 -1 -1
>>
>># In this design matrix, the (Intercept) coefficient is the grand
>>mean, the Sex1 coef is the main effect of Sex,
>># the Time1 and Time2 coef taken together will give you the main
>>effect of Time, and the Sex1:Time1 and
>># Sex1:Time2 coef taken together will give you the interaction term.
>>
>># Now fit the model with your data
>>
>> > fit.2x3 <- eBayes(lmFit(YourData, design))
>>
>># To get the overall F test for the 2x3 take all coef except the Intercept:
>>
>> > overall.2x3 <- topTable(fit, coef=2:6, number=Inf)
>>
>># To get the F test (equivalent to t test) main effect of Sex, do:
>>
>> > mainSex <- topTable(fit, coef=2, number=Inf)
>># Note that the logFC values need to be multiplied by 2 to get the
>>actual Male:Female logFC value!
>>
>># To get the F test for the main effect of Time, do:
>>
>> > mainTime <- topTable(fit, coef=3:4, number=Inf)
>># I usually don't use the actual logFC values listed for each coef,
>>so I usually don't care which group is listed last
>>
>># To get the F test for the interaction term, do:
>>
>> > Interact <- topTable(fit, coef=5:6, number=Inf)
>>
>>
>>If you have 4 groups in a level, you'll have 3 coef for the main
>>effect of the level and 3 coef for the interaction term, and so on
>>up the line. You can also do n-way factorial designs in the same
>>way. Just make sure to keep track of the coefficients!
>>
>>Cheers,
>>Jenny
>>
>>At 03:39 AM 6/29/2011, axel.klenk at actelion.com wrote:
>>>Hi Dana,
>>>
>>>ooops, obviously I shouldn't be posting immediately after lunch...
>>>re-reading
>>>your message after some cups of coffee I realize that my brain has skipped
>>>everything after male/female and treatment/control... so yours is clearly
>>>NOT
>>>a 2x2 design... my apologies...
>>>
>>>Cheers,
>>>
>>> - axel
>>>
>>>
>>>Axel Klenk
>>>Research Informatician
>>>Actelion Pharmaceuticals Ltd / Gewerbestrasse 16 / CH-4123 Allschwil /
>>>Switzerland
>>>
>>>
>>>
>>>
>>>From:
>>>axel.klenk at actelion.com
>>>To:
>>><Dana.Stanley at csiro.au>
>>>Cc:
>>>bioconductor at r-project.org, bioconductor-bounces at r-project.org
>>>Date:
>>>28.06.2011 15:13
>>>Subject:
>>>Re: [BioC] 2 way anova in Bioconductor
>>>Sent by:
>>>bioconductor-bounces at r-project.org
>>>
>>>
>>>
>>>Hi Dana,
>>>
>>>limma can do that easily. The current version of the user's guide
>>>contains one chapter and two case studies on 2x2 factorial designs.
>>>
>>>Cheers,
>>>
>>> - axel
>>>
>>>
>>>Axel Klenk
>>>Research Informatician
>>>Actelion Pharmaceuticals Ltd / Gewerbestrasse 16 / CH-4123 Allschwil /
>>>Switzerland
>>>
>>>
>>>
>>>
>>>
>>>From:
>>><Dana.Stanley at csiro.au>
>>>To:
>>><bioconductor at r-project.org>
>>>Date:
>>>28.06.2011 03:27
>>>Subject:
>>>[BioC] 2 way anova in Bioconductor
>>>Sent by:
>>>bioconductor-bounces at r-project.org
>>>
>>>
>>>
>>>Hi Everyone
>>>
>>>I use custom designed chicken high-density multiplex one channel
>>>microarray (NimbleGen 12x135K). My usual pipeline includes mostly oligo,
>>>arrayQualityMetrics, preprocessCore, genefilter, GeneSelector and limma.
>>>So far I only worked with simple design, 2 conditions; control/treatment
>>>style. Now I have a dataset with embryonic development timeline for male
>>>and female chicks. I still compare different timepoints using limma but I
>>>want to do 2 way anova to identify interaction significant genes, ie genes
>>>
>>>where gender influences development timeline. I looked at many packages
>>>and could not find 2 way anova, though I am sure it is there somewhere. I
>>>tested package ABarray that seemed ideal for my me but after days of
>>>trying I contacted the author to find out that expression sets (mine is
>>>coming from limma) can no longer be used by ABarray as the package has
>>>not been updated for a while. Can anyone suggest a package (and an
>>>example code if possible) for 2 way anova? I am so temp!
>>> ted to go and do it in GeneSpring...
>>>
>>>Thanks all
>>>
>>>Dana
More information about the Bioconductor
mailing list