[Rd] [BioC] limma, FDR, and p.adjust

Gordon K Smyth smyth at wehi.EDU.AU
Mon Dec 20 13:20:39 CET 2004


You asked the same question on the Bioconductor mailing list back in August.  At that time, you
suggested yourself a solution for how the adjusted p-values should be interpreted.  I answered
your query and told you that your interpretation was correct.  So I'm not sure what more can be
said, except that you should read the article Wright (1992), which is cited in the help entry for
p.adjust(), and which explains quite clearly the concept of an adjusted p-value.

The idea that you're having trouble with actually has nothing specifically to do with FDR or with
B&H's (1995) method.  Any adjustment method for multiple testing can be expressed in terms of
adjusted p-values.  The function p.adjust() actually implements several adjustment methods, not
just B&H's, where were not expressed in terms of p-values in their original papers.  The adjust
p-value approach is exactly equivalent to the original formulations, just more flexible.

The situation is not so different with p-values themselves.  Many traditional statistics textbooks
cover hypothesis testing in a way that doesn't mention p-values at all, but the p-value approach
is now generally prefered in software implementations because it so much more flexible.

> Date: Sun, 19 Dec 2004 09:57:43 -0500
> From: "Kimpel, Mark W" <mkimpel at iupui.edu>
> Subject: [BioC] limma, FDR, and p.adjust
> To: <bioconductor at stat.math.ethz.ch>, <r-help at stat.math.ethz.ch>,
> 	<r-devel at stat.math.ethz.ch>
>
> I am posting this to both R and BioC communities because I believe there
> is a lot of confusion on this topic in both communities (having searched
> the mail archives of both) and I am hoping that someone will have
> information that can be shared with both communities.
>
> I have seen countless questions on the BioC list regarding limma
> (Bioconductor) and its calculation of FDR. Some of them involved
> misunderstandings or confusions regarding across which tests the FDR
> "correction" is being applied. My question is more fundamental and
> involves how the FDR method is implemented at the level of "p.adjust"
> (package: stats).
>
> I have reread the paper by Benjamini and Hochberg (1995) and nowhere in
> their paper do they actually "adjust" p values; rather, they develop
> criteria by which an appropriate p value maximum is chosen such that FDR
> is expected to be below a certain threshold.
>
> To try to get a better handle on this, I wrote the following simple
> script to generate a list of random p values, and view it before and
> after apply p.adjust (method=fdr).
>
> rn<-abs(rnorm(100, 0.5, 0.33))
> rn<-rn[order(rn)]
> rn<-rn[1:80]
> rn
> p.adj<-p.adjust(rn, method="fdr")
> p.adj
>
> As you can see after running the code, the p values are truly being
> adjusted, but for what FDR? If I set my p value at 0.05, does that mean
> my FDR is 5%? I have been told by someone that is the case but,
> normally, when discussing FDR, q values are reported or just one p value
> is reported--the threshold for a set FDR. The p.adjust documentation is
> unclear.
>
> For the R developers, I can understand how one would want to include FDR
> procedures in p.adjust, but I wonder, given the numerous FDR algorithms
> now available, if it would be best to formulate an FDR.select function
> that would be option to p.adjust and itself incorporate more recent FDR
> procedures than the one proposed by Benjamini and Hochberg in 1995.
> (Benjamini himself has a newer one). Some of these may currently be
> available as add-on packages but they are not standardized regarding I&O
> and this makes it difficult for developers to incorporate them into
> packages such as limma.

I'm not quite sure what the difficulty is that you see here or how your suggestion would get
around it.  I'm the author of the current p.adjust() code in R as well as the developer of limma. 
I had planned in update to the function to include some more adjustment methods but, as far as I
know, there isn't anything about the interface which is causing developers any problems.

Gordon

> So those are my questions and suggestions,
>
> Thanks,
>
> Mark W. Kimpel MD



More information about the R-devel mailing list