# [BioC] 2 way anova in Bioconductor

Jenny Drnevich drnevich at illinois.edu
Wed Jun 29 16:31:49 CEST 2011

```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...
>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
>
>
>                  [[alternative HTML version deleted]]
>
>_______________________________________________
>Bioconductor mailing list
>Bioconductor at r-project.org
>https://stat.ethz.ch/mailman/listinfo/bioconductor
>Search the archives:
>http://news.gmane.org/gmane.science.biology.informatics.conductor
>
>
>
>
>The information of this email and in any file transmitted with it is
>strictly confidential and may be legally privileged.
>It is intended solely for the addressee. If you are not the intended
>recipient, any copying, distribution or any other use of this email is
>prohibited and may be unlawful. In such case, you should please notify the
>sender immediately and destroy this email.
>The content of this email is not legally binding unless confirmed by
>letter.
>Any views expressed in this message are those of the individual sender,
>except where the message states otherwise and the sender is authorised to
>state them to be the views of the sender's company. For further
>http://www.actelion.com
>
>_______________________________________________
>Bioconductor mailing list
>Bioconductor at r-project.org
>https://stat.ethz.ch/mailman/listinfo/bioconductor
>Search the archives:
>http://news.gmane.org/gmane.science.biology.informatics.conductor
>
>
>
>
>The information of this email and in any file transmitted with it is
>strictly confidential and may be legally privileged.
>It is intended solely for the addressee. If you are not the intended
>recipient, any copying, distribution or any other use of this email
>is prohibited and may be unlawful. In such case, you should please
>notify the sender immediately and destroy this email.
>The content of this email is not legally binding unless confirmed by letter.
>Any views expressed in this message are those of the individual
>sender, except where the message states otherwise and the sender is
>authorised to state them to be the views of the sender's company.