[BioC] limma and p-values

[Gordon K Smyth]

> Date: Wed, 05 Jan 2005 18:46:10 -0800
> From: Giovanni Coppola <gcoppola at ucla.edu>
> Subject: [BioC] limma and p-values
> To: bioconductor at stat.math.ethz.ch
> Message-ID: <6.1.2.0.2.20050105175459.04ccf508 at 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