limma and p-values
2
0
Entering edit mode
@giovanni-coppola-893
Last seen 10.2 years ago
Hi everybody, my question is again related to p-values... Is the topTable function using only 'p.adjust' to adjust the p.values? #R2.0.1, limma 1.8.13 fit<-eBayes(fit) #After applying eBayes, fit$p.value contains the p-values associated with the modified t-statistics p.corr.none <- p.adjust (fit$p.value, "none") p.corr.none.ord<-sort(p.corr.none,decreasing=FALSE,na.last=NA) p.corr.none.ord<-p.corr.none.ord[1:10] toptable.none <-topTable(fit,number=10,adjust="none",sort.by="P") # at this point, 'p.corr.none' = 'toptable.none$P.Value' , whereas... p.corr.fdr <- p.adjust (fit$p.value, "fdr") p.corr.fdr.ord<-sort(p.corr.fdr,decreasing=FALSE,na.last=NA) p.corr.fdr.ord<-p.corr.fdr.ord[1:10] toptable.fdr <-topTable(fit,number=10,adjust="fdr",sort.by="P") #... 'p.corr.fdr' and 'toptable.fdr$P.Value' are not the same, because 'toptable.fdr$P.Value' is much worse...why? Thanks #and sorry for bringing up p-values again :-) Giovanni
limma limma • 3.2k views
ADD COMMENT •
0
Entering edit mode
@gordon-smyth
Last seen 6 hours ago
WEHI, Melbourne, Australia
> Date: Wed, 05 Jan 2005 18:46:10 -0800 > From: Giovanni Coppola <gcoppola@ucla.edu> > Subject: [BioC] limma and p-values > To: bioconductor@stat.math.ethz.ch > Message-ID: <6.1.2.0.2.20050105175459.04ccf508@mail.ucla.edu> > Content-Type: text/plain; charset="us-ascii"; format=flowed > > Hi everybody, > my question is again related to p-values... > Is the topTable function using only 'p.adjust' to adjust the p.values? Yes it uses p.adjust(), and uses only p.adjust(). The mistake you have made below is to apply p.adjust() to the entire matrix of p-values fit$p-value whereas topTable() applies p.adjust() to each column separately. > #R2.0.1, limma 1.8.13 > > fit<-eBayes(fit) > #After applying eBayes, fit$p.value contains the p-values associated with > the modified t-statistics > > p.corr.none <- p.adjust (fit$p.value, "none") > p.corr.none.ord<-sort(p.corr.none,decreasing=FALSE,na.last=NA) > p.corr.none.ord<-p.corr.none.ord[1:10] > toptable.none <-topTable(fit,number=10,adjust="none",sort.by="P") > > # at this point, 'p.corr.none' = 'toptable.none$P.Value' , whereas... > > p.corr.fdr <- p.adjust (fit$p.value, "fdr") > p.corr.fdr.ord<-sort(p.corr.fdr,decreasing=FALSE,na.last=NA) > p.corr.fdr.ord<-p.corr.fdr.ord[1:10] > toptable.fdr <-topTable(fit,number=10,adjust="fdr",sort.by="P") > > #... 'p.corr.fdr' and 'toptable.fdr$P.Value' are not the same, because > 'toptable.fdr$P.Value' is much worse...why? > > Thanks #and sorry for bringing up p-values again :-) > Giovanni It is generally helpful to give an example which is complete and reproducible and which includes the output to prove your point. In this case I can guess that you've fitted a linear model with more than one coefficient. Had there been only one coefficient, the adjusted p-values would have turned out equal: > M <- matrix(rnorm(1000*10),1000,10) > X <- matrix(rnorm(10*1),10,1) > X <- X[,1] > fit <- lmFit(M,X) > fit <- eBayes(fit) > p.corr.fdr <- p.adjust (fit$p.value, "fdr") > p.corr.fdr.ord<-sort(p.corr.fdr,decreasing=FALSE,na.last=NA) > p.corr.fdr.ord<-p.corr.fdr.ord[1:10] > toptable.fdr <-topTable(fit,number=10,adjust="fdr",sort.by="P") > toptable.fdr M t P.Value B 287 -1.1922028 -3.908555 0.1079313 -1.212784 495 -0.9497453 -3.107945 0.7510069 -2.561172 338 -0.9432092 -3.070181 0.7510069 -2.617739 494 0.9171600 2.984433 0.7510069 -2.743777 38 0.8098800 2.660003 0.9661470 -3.190147 582 0.8134650 2.641185 0.9661470 -3.214544 194 -0.7925727 -2.601763 0.9661470 -3.265120 274 -0.7757093 -2.511355 0.9661470 -3.378368 429 -0.7725337 -2.509163 0.9661470 -3.381067 597 -0.7519274 -2.451572 0.9661470 -3.451150 > p.corr.fdr.ord [1] 0.1079313 0.7510069 0.7510069 0.7510069 0.9661470 0.9661470 0.9661470 0.9661470 0.9661470 [10] 0.9661470 Gordon
ADD COMMENT •
0
Entering edit mode
Hi Gordon, thanks for your reply. I am analyzing 12 Agilent slides (6 dye swaps) My previous steps were the following: MA <- normalizeWithinArrays(RG, method="loess") design <- c(1,-1,1,-1,1,-1,1,-1,1,-1,1,-1) fit <- lmFit(MA,design) fit<- eBayes(fit) #is fit$p.value a matrix of p.values? >fit$p.value[1:10] [1] NA NA 0.9931832 0.1992031 0.4111630 0.4880706 NA [8] 0.8261918 0.2254227 0.2039834 my outputs were as follows: > p.corr.none.ord [1] 5.478123e-05 7.317678e-05 8.723562e-05 9.087221e-05 1.133770e-04 [6] 3.923843e-04 4.478298e-04 6.033205e-04 7.023724e-04 8.562249e-04 > toptable.none$P.Value [1] 5.478123e-05 7.317678e-05 8.723562e-05 9.087221e-05 1.133770e-04 [6] 3.923843e-04 4.478298e-04 6.033205e-04 7.023724e-04 8.562249e-04 > p.corr.fdr.ord [1] 0.001140071 0.001521493 0.001812122 0.001885914 0.002350789 0.008128272 [7] 0.009268255 0.012474751 0.014509433 0.017671378 > toptable.fdr$P.Value [1] 0.4833943 0.4833943 0.4833943 0.4833943 0.4833943 0.9990538 0.9990538 [8] 0.9990538 0.9990538 0.9990538 thank you for your time (and patience!) Giovanni At 04:04 AM 1/6/2005, you wrote: > > Date: Wed, 05 Jan 2005 18:46:10 -0800 > > From: Giovanni Coppola <gcoppola@ucla.edu> > > Subject: [BioC] limma and p-values > > To: bioconductor@stat.math.ethz.ch > > Message-ID: <6.1.2.0.2.20050105175459.04ccf508@mail.ucla.edu> > > Content-Type: text/plain; charset="us-ascii"; format=flowed > > > > Hi everybody, > > my question is again related to p-values... > > Is the topTable function using only 'p.adjust' to adjust the p.values? > >Yes it uses p.adjust(), and uses only p.adjust(). > >The mistake you have made below is to apply p.adjust() to the entire >matrix of p-values >fit$p-value whereas topTable() applies p.adjust() to each column separately. > > > #R2.0.1, limma 1.8.13 > > > > fit<-eBayes(fit) > > #After applying eBayes, fit$p.value contains the p-values associated with > > the modified t-statistics > > > > p.corr.none <- p.adjust (fit$p.value, "none") > > p.corr.none.ord<-sort(p.corr.none,decreasing=FALSE,na.last=NA) > > p.corr.none.ord<-p.corr.none.ord[1:10] > > toptable.none <-topTable(fit,number=10,adjust="none",sort.by="P") > > > > # at this point, 'p.corr.none' = 'toptable.none$P.Value' , whereas... > > > > p.corr.fdr <- p.adjust (fit$p.value, "fdr") > > p.corr.fdr.ord<-sort(p.corr.fdr,decreasing=FALSE,na.last=NA) > > p.corr.fdr.ord<-p.corr.fdr.ord[1:10] > > toptable.fdr <-topTable(fit,number=10,adjust="fdr",sort.by="P") > > > > #... 'p.corr.fdr' and 'toptable.fdr$P.Value' are not the same, because > > 'toptable.fdr$P.Value' is much worse...why? > > > > Thanks #and sorry for bringing up p-values again :-) > > Giovanni > >It is generally helpful to give an example which is complete and >reproducible and which includes >the output to prove your point. In this case I can guess that you've >fitted a linear model with >more than one coefficient. Had there been only one coefficient, the >adjusted p-values would have >turned out equal: > > > M <- matrix(rnorm(1000*10),1000,10) > > X <- matrix(rnorm(10*1),10,1) > > X <- X[,1] > > fit <- lmFit(M,X) > > fit <- eBayes(fit) > > p.corr.fdr <- p.adjust (fit$p.value, "fdr") > > p.corr.fdr.ord<-sort(p.corr.fdr,decreasing=FALSE,na.last=NA) > > p.corr.fdr.ord<-p.corr.fdr.ord[1:10] > > toptable.fdr <-topTable(fit,number=10,adjust="fdr",sort.by="P") > > toptable.fdr > M t P.Value B >287 -1.1922028 -3.908555 0.1079313 -1.212784 >495 -0.9497453 -3.107945 0.7510069 -2.561172 >338 -0.9432092 -3.070181 0.7510069 -2.617739 >494 0.9171600 2.984433 0.7510069 -2.743777 >38 0.8098800 2.660003 0.9661470 -3.190147 >582 0.8134650 2.641185 0.9661470 -3.214544 >194 -0.7925727 -2.601763 0.9661470 -3.265120 >274 -0.7757093 -2.511355 0.9661470 -3.378368 >429 -0.7725337 -2.509163 0.9661470 -3.381067 >597 -0.7519274 -2.451572 0.9661470 -3.451150 > > p.corr.fdr.ord > [1] 0.1079313 0.7510069 0.7510069 0.7510069 0.9661470 0.9661470 > 0.9661470 0.9661470 0.9661470 >[10] 0.9661470 > >Gordon
ADD REPLY •
0
Entering edit mode
@gordon-smyth
Last seen 6 hours ago
WEHI, Melbourne, Australia
> Date: Thu, 06 Jan 2005 09:23:01 -0800 > From: Giovanni Coppola <gcoppola@ucla.edu> > Subject: Re: [BioC] limma and p-values > To: bioconductor@stat.math.ethz.ch > > Hi Gordon, > thanks for your reply. > I am analyzing 12 Agilent slides (6 dye swaps) > My previous steps were the following: > > MA <- normalizeWithinArrays(RG, method="loess") > design <- c(1,-1,1,-1,1,-1,1,-1,1,-1,1,-1) > fit <- lmFit(MA,design) > fit<- eBayes(fit) > > #is fit$p.value a matrix of p.values? > >fit$p.value[1:10] > [1] NA NA 0.9931832 0.1992031 0.4111630 0.4880706 NA > [8] 0.8261918 0.2254227 0.2039834 > > my outputs were as follows: > > p.corr.none.ord > [1] 5.478123e-05 7.317678e-05 8.723562e-05 9.087221e-05 1.133770e-04 > [6] 3.923843e-04 4.478298e-04 6.033205e-04 7.023724e-04 8.562249e-04 > > toptable.none$P.Value > [1] 5.478123e-05 7.317678e-05 8.723562e-05 9.087221e-05 1.133770e-04 > [6] 3.923843e-04 4.478298e-04 6.033205e-04 7.023724e-04 8.562249e-04 > > > p.corr.fdr.ord > [1] 0.001140071 0.001521493 0.001812122 0.001885914 0.002350789 0.008128272 > [7] 0.009268255 0.012474751 0.014509433 0.017671378 > > toptable.fdr$P.Value > [1] 0.4833943 0.4833943 0.4833943 0.4833943 0.4833943 0.9990538 0.9990538 > [8] 0.9990538 0.9990538 0.9990538 p.adjust() unfortunately gives incorrect results when 'p' includes NAs. The results from topTable are correct. topTable() takes care to remove NAs before passing the values to p.adjust(). Gordon > thank you for your time (and patience!) > Giovanni > > > At 04:04 AM 1/6/2005, you wrote: >> > Date: Wed, 05 Jan 2005 18:46:10 -0800 >> > From: Giovanni Coppola <gcoppola@ucla.edu> >> > Subject: [BioC] limma and p-values >> > To: bioconductor@stat.math.ethz.ch >> > >> > Hi everybody, >> > my question is again related to p-values... >> > Is the topTable function using only 'p.adjust' to adjust the p.values? >> >>Yes it uses p.adjust(), and uses only p.adjust(). >> >>The mistake you have made below is to apply p.adjust() to the entire >>matrix of p-values >>fit$p-value whereas topTable() applies p.adjust() to each column separately. >> >> > #R2.0.1, limma 1.8.13 >> > >> > fit<-eBayes(fit) >> > #After applying eBayes, fit$p.value contains the p-values associated with >> > the modified t-statistics >> > >> > p.corr.none <- p.adjust (fit$p.value, "none") >> > p.corr.none.ord<-sort(p.corr.none,decreasing=FALSE,na.last=NA) >> > p.corr.none.ord<-p.corr.none.ord[1:10] >> > toptable.none <-topTable(fit,number=10,adjust="none",sort.by="P") >> > >> > # at this point, 'p.corr.none' = 'toptable.none$P.Value' , whereas... >> > >> > p.corr.fdr <- p.adjust (fit$p.value, "fdr") >> > p.corr.fdr.ord<-sort(p.corr.fdr,decreasing=FALSE,na.last=NA) >> > p.corr.fdr.ord<-p.corr.fdr.ord[1:10] >> > toptable.fdr <-topTable(fit,number=10,adjust="fdr",sort.by="P") >> > >> > #... 'p.corr.fdr' and 'toptable.fdr$P.Value' are not the same, because >> > 'toptable.fdr$P.Value' is much worse...why? >> > >> > Thanks #and sorry for bringing up p-values again :-) >> > Giovanni >> >>It is generally helpful to give an example which is complete and >>reproducible and which includes >>the output to prove your point. In this case I can guess that you've >>fitted a linear model with >>more than one coefficient. Had there been only one coefficient, the >>adjusted p-values would have >>turned out equal: >> >> > M <- matrix(rnorm(1000*10),1000,10) >> > X <- matrix(rnorm(10*1),10,1) >> > fit <- lmFit(M,X) >> > fit <- eBayes(fit) >> > p.corr.fdr <- p.adjust (fit$p.value, "fdr") >> > p.corr.fdr.ord<-sort(p.corr.fdr,decreasing=FALSE,na.last=NA) >> > p.corr.fdr.ord<-p.corr.fdr.ord[1:10] >> > toptable.fdr <-topTable(fit,number=10,adjust="fdr",sort.by="P") >> > toptable.fdr >> M t P.Value B >>287 -1.1922028 -3.908555 0.1079313 -1.212784 >>495 -0.9497453 -3.107945 0.7510069 -2.561172 >>338 -0.9432092 -3.070181 0.7510069 -2.617739 >>494 0.9171600 2.984433 0.7510069 -2.743777 >>38 0.8098800 2.660003 0.9661470 -3.190147 >>582 0.8134650 2.641185 0.9661470 -3.214544 >>194 -0.7925727 -2.601763 0.9661470 -3.265120 >>274 -0.7757093 -2.511355 0.9661470 -3.378368 >>429 -0.7725337 -2.509163 0.9661470 -3.381067 >>597 -0.7519274 -2.451572 0.9661470 -3.451150 >> > p.corr.fdr.ord >> [1] 0.1079313 0.7510069 0.7510069 0.7510069 0.9661470 0.9661470 >> 0.9661470 0.9661470 0.9661470 >>[10] 0.9661470 >> >>Gordon
ADD COMMENT •

Login before adding your answer.

Similar Posts
Loading Similar Posts
Traffic: 553 users visited in the last hour
Content Search
Users
Tags
Badges
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6