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 > colnames(cma1) <- c("U","I") > dma1 <- makeContrasts(U-I,levels=cma1) > cma1 U I 1 1 0 2 0 1 3 0 1 4 1 0 5 0 1 6 1 0 8 0 1 9 1 0 10 0 1 > cl <- cma1[,2] > tmp3 <- mt.maxT(abinorms3,cl) We'll do complete enumerations We're doing 126 complete permutations b=1 b=2 b=3 b=4 b=5 b=6 b=7 b=8 b=9 b=10 b=11 b=12 b=13 b=14 b=15 b=16 b=17 b=18 b=19 b=20 b=21 b=22 b=23 b=24 b=25 b=26 b=27 b=28 b=29 b=30 b=31 b=32 b=33 b=34 b=35 b=36 b=37 b=38 b=39 b=40 b=41 b=42 b=43 b=44 b=45 b=46 b=47 b=48 b=49 b=50 b=51 b=52 b=53 b=54 b=55 b=56 b=57 b=58 b=59 b=60 b=61 b=62 b=63 b=64 b=65 b=66 b=67 b=68 b=69 b=70 b=71 b=72 b=73 b=74 b=75 b=76 b=77 b=78 b=79 b=80 b=81 b=82 b=83 b=84 b=85 b=86 b=87 b=88 b=89 b=90 b=91 b=92 b=93 b=94 b=95 b=96 b=97 b=98 b=99 b=100 b=101 b=102 b=103 b=104 b=105 b=106 b=107 b=108 b=109 b=110 b=111 b=112 b=113 b=114 b=115 b=116 b=117 b=118 b=119 b=120 b=121 b=122 b=123 b=124 b=125 b=126 > 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 > fit1s <- lmFit(abinorms3,cma1) > fit2s <- contrasts.fit(fit1s,dma1) > fit3s <- eBayes(fit2s) > topTable(fit3s) M t P.Value B 8064 -3.155503 -7.695828 0.2987386 1.52861663 105 1.523110 6.053780 1.0000000 0.46775473 2865 1.512303 5.764098 1.0000000 0.23825116 16035 3.125686 5.754778 1.0000000 0.23063006 1074 1.394530 5.602738 1.0000000 0.10416651 34 1.428812 5.595450 1.0000000 0.09800220 10283 1.316832 5.591068 1.0000000 0.09429127 4638 1.365565 5.530195 1.0000000 0.04239040 8571 -4.982818 -5.516468 1.0000000 0.03059528 8588 1.430188 5.472810 1.0000000 -0.00714332 [[alternative HTML version deleted]]

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

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@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 >