Question: limma roast syntax for overall anova
0
5.4 years ago by
Juliet Hannah360
United States
Juliet Hannah360 wrote:

All,

I am testing an overall anova (4 groups; is there any difference across groups) using the following design matrix in limma

  (Intercept) B C D
1           1 0 0 0
2           1 0 0 0
3           1 0 0 0
4           1 0 0 0
5           1 0 0 0
6           1 0 0 0

I have used ROAST (in limma) in the past with a specific contrast, such as

roast(expression_matrix,index,design_matrix,contrast=10,nrot=3000).

In the current case, I didn't need any contrasts. In limma, I had used

design <- model.matrix(~key$stage) colnames(design)[2:4] = c("B","C","D") fit <- lmFit(e, design) fit <- eBayes(fit) topTable(fit,adjust.method="BH") to obtain the overall F-test (I understand this test is not that interesting, but it is requested of me). Using this setup, how can I specify this overall test within ROAST. Thanks for your help. Juliet limma • 1.6k views ADD COMMENTlink modified 5.1 years ago by Gordon Smyth39k • written 5.4 years ago by Juliet Hannah360 Answer: limma roast syntax for overall anova 3 5.4 years ago by Gordon Smyth39k Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia Gordon Smyth39k wrote: Dear Juliet, To get the overall anova F-test you need:  topTable(fit[,-1]) or alternatively  topTable(fit, coef=2:4) which gives the same result. Note that the intercept needs to be removed from the test, as it always is in classical anova, and that's what the -1 does. You also ask "Using this setup, how can I specify this overall test within ROAST", and that is a good question. We haven't implemented F-tests for ROAST. If we get enough requests, we will consider doing so, but we haven't so far. If the meantime, you can input the F-statistics to the geneSetTest() function:  fit2 <- fit[,-1] geneSetTest(index, fit2$F)

etc.

Best wishes
Gordon

Hi Gordon,

Thanks for your help (and software!).

I confusingly showed only the first few rows of the design matrix. I was trying to convey that I left the intercept in.

limma gave me the message:

> mytt <- topTable(fit,number=Inf)
Removing intercept from test coefficients

which gives me the same result as the more clear version you gave me

 mytt <- topTable(fit,coef=2:4,number=Inf).

I will proceed with geneSetTest.

Thanks,

Juliet

ADD REPLYlink modified 5.1 years ago by Gordon Smyth39k • written 5.4 years ago by Juliet Hannah360

Oh, yes, that's right. I had forgotten that limma now recognizes the special intercept column name "(Intercept)" and removes it automatically. So you already had the correct anova test.

Best wishes
Gordon