[BioC] p-values from robust linear model in LIMMA

[Arvid Sondén]

Dear Gordon,

Thank you for your answer! 

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

-----Original Message-----
From: Gordon K Smyth [mailto:smyth at wehi.EDU.AU] 
Sent: den 3 juni 2014 02:10
To: Arvid Sondén
Cc: Bioconductor mailing list
Subject: p-values from robust linear model in LIMMA

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

> Date: Mon, 2 Jun 2014 07:23:04 +0000
> From: Arvid Sond?n <arvid.sonden at gu.se>
> To: "bioconductor at r-project.org" <bioconductor at r-project.org>
> Subject: [BioC] p-values from robust linear model in LIMMA
>
> 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")

You could simply use lmFit(y, design=m.matrix, etc).  No need for y$E.

> 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

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