Search
Question: Two-way anova 2x2 design with limma
0
gravatar for damoros
4 days ago by
damoros0
damoros0 wrote:

Hi,
In the 8.7 section (Two-factor example) in limma User Guide, explains three procedures that give the same result.
In the first procedure:

> cont.matrix <- makeContrasts(
+     SvsUinWT=WT.S-WT.U,
+     SvsUinMu=Mu.S-Mu.U,
+     Diff=(Mu.S-Mu.U)-(WT.S-WT.U),
+     levels=design)
> fit2 <- contrasts.fit(fit, cont.matrix)
> fit2 <- eBayes(fit2)

How can I find the statistical F of two-way anova (genotype, treatment and interaction effect)? I need a factorial model formulas (second and third procedure)?

Thank you

 

ADD COMMENTlink modified 13 hours ago • written 4 days ago by damoros0

You should read the limma User's Guide that comes with the version of the limma package you are using. Or else the latest version of the Guide is always available from:

https://bioconductor.org/packages/release/bioc/vignettes/limma/inst/doc/usersguide.pdf

Your reference to Section 8.7 suggests that you are reading a limma User's Guide from many, many years ago.

ADD REPLYlink modified 4 days ago • written 4 days ago by Gordon Smyth32k
1
gravatar for Aaron Lun
4 days ago by
Aaron Lun17k
Cambridge, United Kingdom
Aaron Lun17k wrote:

There is no Section 8.7 in the user's guide. The example code seems to come from 9.5.2. If you want a two-way ANOVA, you will want to drop all non-intercept terms, so that the null is that all factor combinations have the same mean expression. This can be done by setting coef=2:4 in topTable with either of the nested interaction or in the classic interaction parameterisations of that specific experimental design.

ADD COMMENTlink modified 4 days ago • written 4 days ago by Aaron Lun17k
0
gravatar for damoros
13 hours ago by
damoros0
damoros0 wrote:

Hi.

Thank you for your answers. I already have the latest version of user guide.
I understand that the procedure in 9.5.2 is the recommended one (single factor), but I am required to perform two way anova.

For example, if I extract nomalized intensities of one probe, I can perform a two way anova for that probe:

>mat
     data       Strain    Treatment
1  15087.105     WT         U
2  13930.749     WT         U
3  16495.873     WT         U
4  13829.101     WT         U
5  12628.449     WT         S
6  10274.768     WT         S
7  11306.912     WT         S
8  11842.572     WT         S
9  16296.351     Mu         U
10 17066.216     Mu         U
11 15817.936     Mu         U
12 15001.087     Mu         U
13  9952.837     Mu         S
14 12524.801     Mu         S
15 10785.661     Mu         S
16 13326.560     Mu         S
> res.aov=aov(data~Strain*Treatment,data=mat)
> summary(res.aov)
                 Df   Sum Sq  Mean Sq F value   Pr(>F)    
Strain            1  1806283  1806283   1.273    0.281    
Treatment         1 59605572 59605572  41.998 3.02e-05 ***
Strain:Treatment  1  1156486  1156486   0.815    0.384    
Residuals        12 17031010  1419251                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

it is right? I need to obtain the associated F and  p-values of strain, associated F and  p-values of treatment and associated F and  p-values of interaction effect in the same way that they are obtained from the aov function, using limma and its moderated F statistics. I understand that with topTable(fit2,coef=2:4) I get the results of global or omnibus F test (I'm not interested in that test).

 

I understand that I should use the classic interaction model (9.5.4 section):

> design <- model.matrix(~Strain*Treatment)

   (Intercept) StrainWT TreatmentU StrainWT:TreatmentU
1            1        1          1                   1
2            1        1          1                   1
3            1        1          1                   1
4            1        1          1                   1
5            1        1          0                   0
6            1        1          0                   0
7            1        1          0                   0
8            1        1          0                   0
9            1        0          1                   0
10           1        0          1                   0
11           1        0          1                   0
12           1        0          1                   0
13           1        0          0                   0
14           1        0          0                   0
15           1        0          0                   0
16           1        0          0                   0

> fit <- lmFit(eset, design)
> fit2 <- eBayes(fit2)

How to continue?

 

ADD COMMENTlink modified 12 hours ago • written 13 hours ago by damoros0
  1. Use "Add comment" to respond to existing answers. Don't "Add your answer" unless you're answering your own question.
  2. F-statistics and p-values can be obtained with appropriate settings of coef in topTable. For example, coef=2 will be equivalent to Strain, coef=3 will be equivalent to Treatment, and coef=4 will equivalent to testing the interaction effect. (I say equivalent, not identical, due to the effect of empirical Bayes shrinkage in limma that isn't available with aov.)
  3. Having said that, keep in mind that StrainWT and TreatmentU do not represent the main effects of the strain or treatment when the interaction term is non-zero. When there is an interaction between the factors, StrainWT represents the log-fold change of WT over Mu in treatment S only, while TreatmentU represents the log-fold change of U over S in the Mu strain only. The same applies for the output of the aov call, which also uses model.matrix internally.
ADD REPLYlink modified 12 hours ago • written 12 hours ago by Aaron Lun17k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.2.0
Traffic: 145 users visited in the last hour