[BioC] EdgeR multi-factor testing questions

[Gordon K Smyth]

Dear Yanzhu,

Your analysis is fine from a code point of view.  From a statistical point 
of view however your analysis is too simple because you are neglecting the 
principle of marginality:

   http://en.wikipedia.org/wiki/Principle_of_marginality

For the model you have fitted, it makes sense to test for the three-way 
interaction as you do.  However it does not make statistical sense to test 
for the main effects or two-interactions until you have established that 
the three-way interaction is non-significant.

For count data, the tests for the lower-level interactions need to be 
computed by successively removing each level of interactions from the 
model.  See for example:

   https://stat.ethz.ch/pipermail/bioconductor/2013-December/056584.html

This is the same as the anova() function does in R for unbalanced linear 
factorial models.

Furthermore, testing the two-way interations is only sensible for genes 
with non-signicant 3-way interactions.  Similarly, testing the main effect 
is only sensible for genes with non-significant 2-way and 3-way 
interactions.  Otherwise these tests have no useful scientific meaning.

This is a basic drawback of the factorial anova approach.  You might 
consider the alternative approach described in Section 3.3.1 of the edgeR 
User's Guide.

Best wishes
Gordon

> Date: Mon,  6 Jan 2014 11:21:37 -0800 (PST)
> From: "Yanzhu [guest]" <guest at bioconductor.org>
> To: bioconductor at r-project.org, mlinyzh at gmail.com
> Subject: [BioC] EdgeR multi-factor testing questions
>
>
> I'm currently using EdgeR to analyze RNASeq data. I have one question of 
> whether I have chosen the correct method for analyzing my 
> multi-factorial experiment.

> I have 3 factors: factor A with 16 levels (from level0 to level15), 
> factor B with three levels (from level0 to leve2), and factor Sex with 
> two levels (male and female).  I'm interested in testing all main 
> effects: A, B and Sex, all two-way interaction terms:  A:B, A:Sex and 
> B:Sex; and the three-way interaction term: A:B:Sex.
>
> My code:
>
> group<-paste(A,B,Sex,sep=".")
> design<-model.matrix(~A+B+Sex+A:B+A:Sex+B:Sex+A:B:Sex)
> y<-DGEList(counts=readcounts,group=group)
> y<-calcNormFactors(y,method="TMM")
>
>
> y<-estimateGLMCommonDisp(y,design)
> y<-estimateGLMTagwiseDisp(y,design)
>
>
> fit<-glmFit(y,design)
> colnames(fit)
>
> colnames(fit)
> [1] "(Intercept)" "Alevel1" "Alevel2" "Alevel3" "Alevel4"        [6] "Alevel5" "Alevel6" "Alevel7" "Alevel8" "Alevel9"
> [11]"Alevel10" "Alevel11" "Alevel12" "Alevel13" "Alevel14"
> [16] "Alevel15" "Blevel1" "Blevel2" "Sexmale" "Alevel1:Blevel1"
>
> [21] "Alevel2:Blevel1" "Alevel3:Blevel1" "Alevel4:Blevel1" "Alevel5:Blevel1" "Alevel6:Blevel1"
>
> [26] "Alevel7:Blevel1" "Alevel8:Blevel1" "Alevel9:Blevel1" "Alevel10:Blevel1" "Alevel11:Blevel1"
>
> [31] "Alevel12:Blevel1" "Alevel13:Blevel1" "Alevel14:Blevel1" "Alevel15:Blevel1" "Alevel1:Blevel2"
>
> [36] "Alevel2:Blevel2"  "Alevel3:Blevel2" "Alevel4:Blevel2" "Alevel5:Blevel2" "Alevel6:Blevel2"
>
> [41] "Alevel7:Blevel2"  "Alevel8:Blevel2" "Alevel9:Blevel2"   "Alevel10:Blevel2" "Alevel11:Blevel2"
>
> [46] "Alevel12:Blevel2" "Alevel13:Blevel2" "Alevel14:Blevel2" "Alevel15:Blevel2" "Alevel1:Sexmale"
>
> [51] "Alevel2:Sexmale" "Alevel3:Sexmale" "Alevel4:Sexmale" "Alevel5:Sexmale" "Alevel6:Sexmale"
>
> [56] "Alevel7:Sexmale" "Alevel8:Sexmale" "Alevel9:Sexmale" "Alevel10:Sexmale" "Alevel11:Sexmale"
>
> [61] "Alevel12:Sexmale" "Alevel13:Sexmale" "Alevel14:Sexmale" "Alevel15:Sexmale" "Blevel1:Sexmale"
>
> [66] "Blevel2:Sexmale" "Alevel1:Blevel1:Sexmale" "Alevel2:Blevel1:Sexmale" "Alevel3:Blevel1:Sexmale" "Alevel4:Blevel1:Sexmale"
>
> [71] "Alevel5:Blevel1:Sexmale" "Alevel6:Blevel1:Sexmale" "Alevel7:Blevel1:Sexmale" "Alevel8:Blevel1:Sexmale" "Alevel9:Blevel1:Sexmale"
>
> [76] "Alevel10:Blevel1:Sexmale" "Alevel11:Blevel1:Sexmale" "Alevel12:Blevel1:Sexmale" "Alevel13:Blevel1:Sexmale" "Alevel14:Blevel1:Sexmale"
>
> [81] "Alevel15:Blevel1:Sexmale" "Alevel1:Blevel2:Sexmale" "Alevel2:Blevel2:Sexmale" "Alevel3:Blevel2:Sexmale" "Alevel4:Blevel2:Sexmale"
>
> [86] "Alevel5:Blevel2:Sexmale" "Alevel6:Blevel2:Sexmale" "Alevel7:Blevel2:Sexmale" "Alevel8:Blevel2:Sexmale" "Alevel9:Blevel2:Sexmale"
>
> [91] "Alevel10:Blevel2:Sexmale" "Alevel11:Blevel2:Sexmale" "Alevel12:Blevel2:Sexmale" "Alevel13:Blevel2:Sexmale" "Alevel14:Blevel2:Sexmale"
> [96] "Alevel15:Blevel2:Sexmale"
>
> (1) To test the main effect, I use the following code:
> lrt_A<-glmLRT(fit,coef=c(2:16))
> lrt_B<-glmLRT(fit,coef=c(17:18))
> lrt_Sex<-glmLRT(fit,coef=19)
>
> (2) To test the two-way interaction terms, I use the following code:
> lrt_AB<-glmLRT(fit,coef=c(20:49))
>
> lrt_ASex<-glmLRT(fit,coef=c(50:64))
>
> lrt_BSex<-glmLRT(fit,coef=c(65:66))
>
> (3) To test the three-way interaction terms, I use the following code:
> lrt_ABSex<-glmLRT(fit,coef=c(67:96))
>
> My question: is my code correct for each of the testing above?
>
>
>
> Thank you!
>
>
>
>
>
>
>
>
>
>
> -- output of sessionInfo():
>
>>  sessionInfo()
> R version 3.0.1 (2013-05-16)
> Platform: x86_64-w64-mingw32/x64 (64-bit)
>
> locale:
> [1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252    LC_MONETARY=English_United States.1252
> [4] LC_NUMERIC=C                           LC_TIME=English_United States.1252
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] edgeR_3.2.4  limma_3.16.8
>
> loaded via a namespace (and not attached):
> [1] tools_3.0.1
>
>
> --
> Sent via the guest posting facility at bioconductor.org.

______________________________________________________________________
The information in this email is confidential and intend...{{dropped:4}}