Question: limma roast syntax for overall anova
0
gravatar for Juliet Hannah
5.1 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.5k views
ADD COMMENTlink modified 4.8 years ago by Gordon Smyth38k • written 5.1 years ago by Juliet Hannah360
Answer: limma roast syntax for overall anova
3
gravatar for Gordon Smyth
5.1 years ago by
Gordon Smyth38k
Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia
Gordon Smyth38k 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

ADD COMMENTlink modified 4.8 years ago • written 5.1 years ago by Gordon Smyth38k

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 4.8 years ago by Gordon Smyth38k • written 5.1 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

ADD REPLYlink modified 4.8 years ago • written 5.1 years ago by Gordon Smyth38k
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 16.09
Traffic: 289 users visited in the last hour