[BioC] limma and p-values

[Gordon K Smyth]

> Date: Thu, 06 Jan 2005 09:23:01 -0800
> From: Giovanni Coppola <gcoppola at ucla.edu>
> Subject: Re: [BioC] limma and p-values
> To: bioconductor at 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 at ucla.edu>
>> > Subject: [BioC] limma and p-values
>> > To: bioconductor at 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