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

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:

  http://www.pnas.org/content/100/9/5491.short

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:

  http://www.statsci.org/smyth/pubs/RobustEBayesPreprint.pdf

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