p-values from robust linear model in LIMMA
1
1
Entering edit mode
@arvid-sonden-6581
Last seen 6.7 years ago

Dear all,

I am currently working with expression data in LIMMA and have been asked to fit a robust linear model:

lm.fit.rob <- lmFit(object=y\$E,design=m.matrix, method="robust")
lm.fit.rob.bayes <- eBayes(lm.fit.rob)
lm.fit.rob.bayes.tt <- topTable(lm.fit.rob.bayes,coef="Group")

It is easy to find that lmFit uses mrlm, and that mrlm uses rlm. However, rlm does not generate p-values, and from what I have read (e.g. http://r.789695.n4.nabble.com/p-values-td803236.html ) it is not a trivial thing to do. What confuses me is that topTable generates p-values, and I can't find in the documentation how and on which assumptions.

The robust models do greatly improve my p-values, but can they be trusted?

Best regards,

Arvid

limma lmFit • 3.1k views
3
Entering edit mode
@gordon-smyth
Last seen 1 hour ago
WEHI, Melbourne, Australia

Dear Arvid,

You are right that the theory of hypothesis testing for robust regression is problematic. What limma does is to simply take the scale estimates and robustifying weights from rlm() and input them into the usual limma pipeline. This is intuitively reasonable because it downweights values that have been identified as potential outliers, and I wanted to include rlm as an option because of some early microarray studies that used it:

However the approximation to a t-distribution is fairly rough and I can't give you good guidelines as to how reliable the p-values are. The idea is that you can still interpret the logFC, or use the B-statistics as a gene ranking, even if you don't trust the p-values.

An alternative is to go back to lmFit() with method="ls" but to use robust=TRUE (and perhaps trend=TRUE as well) at the eBayes() step instead. This approach does have a rigorous theoretical underpinning:

In practice, my lab doesn't use method="robust" for our own analyses. When we want to robustify, we do it at the eBayes step instead.

Best wishes
Gordon

0
Entering edit mode

Dear Gordon,

My issue is that I want to use my p-values for filtering, but when using method="ls" (even if using robust=TRUE in eBayes) I get no significant values after adjusting for multiple testing at all. Using the robust method in lmFit however generates a distribution of p-values much closer to what I would expect. It also gives me enough significant genes to work with in the following pathway analysis.

I understand that the p-values I get are not to be trusted if interpreting them as actual p-values, but would you also say that they cannot be used for filtering?

Best regards,
Arvid

0
Entering edit mode

When you say "My issue is that I want to use my p-values for filtering," do you mean gene ranking? Yes, you could use them for ranking.

Gordon