Search
Question: limma, FDR, and p.adjust
0
14.0 years ago by
Kimpel, Mark W890 wrote:
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. So those are my questions and suggestions, Thanks, Mark W. Kimpel MD
ADD COMMENTlink
modified 14.0 years ago by Gordon Smyth35k • written 14.0 years ago by Kimpel, Mark W890
0
14.0 years ago by
Marcus Davy680
Marcus Davy680 wrote:
Mark, there is a fdr website link via Yoav Benjamini's homepage which is: http://www.math.tau.ac.il/%7Eroee/index.htm 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%. e.g. 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. marcus >>> "Kimpel, Mark W" <mkimpel@iupui.edu> 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 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. So those are my questions and suggestions, Thanks, Mark W. Kimpel MD _______________________________________________ Bioconductor mailing list Bioconductor@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/bioconductor ______________________________________________________ The contents of this e-mail are privileged and/or confidenti...{{dropped}}
ADD COMMENTlink written 14.0 years ago by Marcus Davy680
I have been asked to look some data generated on affy-platform. Due to biological reasons the data owners want to look at PM values only, i.e ignore the MM values. Hence my question is: What is the best option available for normalizing affy data using the PM intensities only. I appologize if this question has been asked and answered on this group earlier. Narinder S. Sahni
ADD REPLYlink written 14.0 years ago by narinder.singh40
0
14.0 years ago by
Gordon Smyth35k
Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia
Gordon Smyth35k wrote:
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@iupui.edu> > Subject: [BioC] limma, FDR, and p.adjust > To: <bioconductor@stat.math.ethz.ch>, <r-help@stat.math.ethz.ch>, > <r-devel@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
ADD COMMENTlink written 14.0 years ago by Gordon Smyth35k
Please log in to add an answer.

Content
Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.2.0
Traffic: 160 users visited in the last hour