[BioC] 2 way anova in Bioconductor

[Dana.Stanley at csiro.au]

Thank you all for your help, especially Jenny for such a detailed explanation!!! I followed your email and got it to work!!! 

Thanks again

Dana


-----Original Message-----
From: Natasha Sahgal [mailto:nsahgal at well.ox.ac.uk] 
Sent: Thursday, 30 June 2011 7:27 PM
To: Jenny Drnevich
Cc: Stanley, Dana (LI, Geelong AAHL); bioconductor at r-project.org; bioconductor-bounces at r-project.org
Subject: Re: [BioC] 2 way anova in Bioconductor

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