atom feed1 message in org.r-project.r-devel[Rd] RE: [BioC] limma, FDR, and p.adjust
FromSent OnAttachments
Kimpel, Mark WDec 20, 2004 2:53 pm 
Subject:[Rd] RE: [BioC] limma, FDR, and p.adjust
From:Kimpel, Mark W (
Date:Dec 20, 2004 2:53:14 pm

Marcus, Gordon, and R developers,

Thanks to you both for kind and useful replies. Gordon you are correct, I did post this question once before but, as I approach publication of a rapid communication paper, wanted to make absolutely sure that I was correct in assuming that the p values of p.adjust were the same as the q values published in the FDR literature that I am familiar with.

If I seem a bit stubborn in pursuing this once again, it is because I have spent a fair amount of time trying to explain the idea of FDR to my lab-mates and how it differs from the concept of FWER. All of the previous p adjust functions, to my knowledge, correct for FWER; the authors of many of the FDR papers have gone to great lengths to distance themselves from some of the terminology associated with FWER so that the FDR methods are not confused with them.

My point to the R and BioC communities was that, if I was correct that p.adjust was somehow adjusting the p values using different FDR methodologies, that this should not be confused with how FDR is discussed in the original papers. I have seen a recent post on BioC where the originator asked at what FDR level the p values were being adjusted to (still thinking, I believe, about the p value threshold concept), so I don't think I am the only one who has confused by this.

As I pointed out, I do see the utility for developers, such as Gordon, for the ability to call just one function when dealing with adjusting p values. Given the proliferation of FDR methods, however, and to allow for implementation of many of them in the future, I would propose that the R-developers use a function architecture such as this:

p.adjust (method = c("bonferronni, ...., fdr(method))), where fdr called a function fdr(method = "B&H, "B&Y",....). In the help section of fdr the method for turning q values into p values could be explained for each of the differing fdr methodologies.

Gordon, limma is a wonderful package and I will be publishing some very exciting findings using (and citing) it. Might I suggest adding a paragraph to the user's guide about fdr and p values? Having raised this issue, I would be more than happy to contribute a first draft that you could use or not use as you see fit.

Marcus, thank you so much for you code. This is an fdr implementation that I think is clearer and more functional than available under p.adjust but I couldn't get my code to make it work.

Thank you both for your help to me and support of the BioC community.

Mark W. Kimpel MD

(317) 490-5129 Home, Work, & Mobile

(317) 278-4104 FAX

-----Original Message----- From: Marcus Davy [] Sent: Monday, December 20, 2004 12:15 AM To: Kimpel, Mark W;; Subject: Re: [BioC] limma, FDR, and p.adjust

Mark, there is a fdr website link via Yoav Benjamini's homepage which is: On it you can download an S-Plus function (under the downloads link) which calculates the false discovery rate threshold alpha level using stepup, stepdown, dependence methods etc. Some changes are required to the plotting code when porting it to R. I removed the xaxs="s" arguement on line 80. The fdr function requires a list of p-values as input, a Q-value (*expected* false discovery rate control at level Q) and a required method of fdr controlling procedure.

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.

The p.adjust function appears to be using the "stepup" fdr controlling procedure when method="fdr" is specified. It adjusts the p-values so that FDR control is at the desired level alpha over the entire range (0,1), which gives the same result as specifying a Q-value in the fdr function itself calculating a false discovery rate threshold alpha level so that FDR<=Q.

So it adjusts for all FDR desired levels. If your p-value threhold is 0.05 then the expected proportion of false discoveries is 5%.


n <- 1000 pi0 <- 0.5 x <- rnorm(n, mean=c(rep(0, each=n*pi0), rep(3, each=n - (n*pi0)))) p <- 2*pnorm( -abs(x)) p <- sort(round(p,3))

p.adjusted <- p.adjust(p, method="fdr")

# Controlling fdr at Q, and p.adjust at level alpha Qvalue <- alpha <- 0.05

threshold <- fdr(p, Q=Qvalue, method="stepup") # fdr function available from the website link above threshold

plot(p, p.adjusted) abline(v=threshold, lty=2) abline(h=alpha, lty=2)

# Stepup FDR control at Q=0.05 sum(p <= threshold) [1] 372

# p.adjust(ed) p-values at level alpha=0.05 sum(p.adjusted <= alpha)

[1] 372

Simultaneously modifying Qvalue, and alpha above to a different expected proportion of false discoveries should still produce identically sized rejected lists.

Hope that helps.


"Kimpel, Mark W" <> 20/12/2004 3:57:43 AM >>>

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 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.

So those are my questions and suggestions,


Mark W. Kimpel MD

The contents of this e-mail are privileged and/or confidenti...{{dropped}}