[BioC] limma/multtest comparison question

[Sean Davis]

Thanks to Gordon for the reply, but that doesn't appear to be the answer to
the difference.  As I mentioned, I have the same samples on another array
platform that give expected results in terms of number of significant
results.  Below are the results for the two data sets (using adj="fdr")
using limma and mt.maxT for both sets.  I am NOT concerned about the
differences between p-values for mt.maxT and limma, but for the relative
differences between the two data sets when using limma, as the differences
are small when using mt.maxT.  You can see the relative agreement in the
mt.maxT p-values (at least order of magnitude) and t-statistics, but there
are two orders of magnitude difference in p-values between limma results
despite similar t-statistics. I don't understand enough about the bayesian
shrinkage of the variance to know why the difference in significance levels.
I would expect the p-values to typically be the same or better than
unmoderated t-tests, though, when using limma (and they have been with all
data sets that I have tried so far), particularly since maxT correction is
so conservative compared to fdr.  Any ideas why one data set would behave so
differently in limma as compared to another "non-bayesian" t-test?


SINGLE CHANNEL--IN QUESTION
> topTable(fit3s,adj="fdr")
              M         t   P.Value           B
8064  -3.155503 -7.695828 0.2987386  1.52861663
105    1.523110  6.053780 0.3609049  0.46775473
2865   1.512303  5.764098 0.3609049  0.23825116
16035  3.125686  5.754778 0.3609049  0.23063006
1074   1.394530  5.602738 0.3609049  0.10416651
34     1.428812  5.595450 0.3609049  0.09800220
10283  1.316832  5.591068 0.3609049  0.09429127
4638   1.365565  5.530195 0.3609049  0.04239040
8571  -4.982818 -5.516468 0.3609049  0.03059528
8588   1.430188  5.472810 0.3609049 -0.00714332
> tmp3[order(tmp3$rawp)[1:10],]
      index  teststat        rawp      adjp
5928   5928 -9.184354 0.007936508 0.3809524
8064   8064  7.712025 0.007936508 0.6190476
2878   2878 -7.574579 0.007936508 0.6825397
2503   2503 -7.491922 0.007936508 0.6984127
10283 10283 -7.415540 0.007936508 0.6984127
7927   7927 -7.213376 0.007936508 0.7222222
4638   4638 -7.091922 0.007936508 0.7301587
105     105 -7.076814 0.007936508 0.7460317
1074   1074 -7.061793 0.007936508 0.7539683
17870 17870  6.967126 0.007936508 0.7777778

TWO-CHANNEL
> topTable(fit3q,adj="fdr")
              M         A         t    P.Value        B
16778  2.040726 10.089056  9.751891 0.02387266 5.717952
31021 -1.373636 11.422162 -8.926289 0.02498543 5.014448
8623   2.191777  9.175746  8.263818 0.02498543 4.391696
33887  1.628374 10.284200  8.135795 0.02498543 4.264907
28057 -1.918088 10.782108 -7.915738 0.02498543 4.041876
33805  1.904995 10.899343  7.760268 0.02498543 3.880332
16076  1.917711 12.228181  7.752972 0.02498543 3.872669
13859 -1.366272 15.020477 -7.617720 0.02498543 3.729262
16055  1.386332 11.150963  7.524214 0.02498543 3.628605
28125  1.643313 11.284236  7.492197 0.02498543 3.593852
> tmp4[order(tmp4$rawp)[1:10],]
      index  teststat        rawp       adjp
10270 10270 10.921859 0.007936508 0.07142857
15137 15137 10.446857 0.007936508 0.09523810
1390   1390 -9.481394 0.007936508 0.15873016
20977 20977  9.272871 0.007936508 0.19047619
9801   9801  9.130909 0.007936508 0.22222222
10101 10101  8.619132 0.007936508 0.26190476
21622 21622  8.338045 0.007936508 0.30158730
7665   7665  8.331140 0.007936508 0.30952381
5156   5156  8.220101 0.007936508 0.32539683
9818   9818  8.122754 0.007936508 0.34126984

Thanks for any insight you have.

Sean

On 6/17/04 10:52 PM, "Gordon Smyth" <smyth at wehi.edu.au> wrote:

> topTable() in limma always does adjustment for multiple testing. The
> default is the very conservative 'Holm" method. You probably want to choose
> something less conservative.
> 
> Gordon
> 
> At 12:40 PM 18/06/2004, Sean Davis wrote:
>> Dear all,
>> 
>> I have a limma/multtest question.  I have performed what I think is the
>> equivalent of a t-test (of course, moderated) in limma below.  Compare
>> that with the output from mt.maxT in the multtest package.  Why are there
>> such large discrepancies (ie., what am I doing wrong), particularly in
>> p-values?  The data (abinorms3) are single-channel log2 intensities.  The
>> commands and output are below.  We have the same samples hybridized on
>> two-color arrays, also analyzed with limma, with p-values in mt.maxT and
>> limma that agree with the mt.maxT result from below, so I think there are
>> differentially-expressed genes.
>> 
>> Thanks,
>> Sean
>