Difference in mixed test for fry vs mroast
1
0
Entering edit mode
@james-w-macdonald-5106
Last seen 9 hours ago
United States

I didn't expect these differences between fry and mroast for the mixed test. Here's a reproducible example.


> set.seed(0xabeef)
> y <- matrix(rnorm(8e4), ncol = 8)
>  design <- model.matrix(~gl(2,4))
> fit <- eBayes(lmFit( y, design))
> lst <- list(A = sample(1:1e4, 20), B = sample(1:1e4, 50))
> fry(y, lst, design)
  NGenes Direction PValue   FDR PValue.Mixed FDR.Mixed
B     50      Down  0.218 0.352        0.661     0.661
A     20        Up  0.352 0.352        0.407     0.661
> mroast(y, lst, design, nrot = 1e5)
  NGenes PropDown PropUp Direction PValue   FDR PValue.Mixed FDR.Mixed
B     50     0.08   0.06      Down  0.217 0.352        0.705     0.705
A     20     0.05   0.15        Up  0.352 0.352        0.237     0.474

For an actual analysis I am doing the differences are even starker

> mroast(dglst, subkegg[1:2], desgn, coef, "SYMBOL")
                             NGenes   PropDown     PropUp Direction PValue
Glycolysis_/_Gluconeogenesis     43 0.09302326 0.02325581      Down  0.701
Citrate_cycle_(TCA_cycle)        27 0.14814815 0.11111111        Up  0.961
                                 FDR PValue.Mixed FDR.Mixed
Glycolysis_/_Gluconeogenesis 0.96075       0.4790    0.4790
Citrate_cycle_(TCA_cycle)    0.96100       0.1925    0.3845

> fry(dglst, subkegg[1:2], desgn, coef, "SYMBOL")
                             NGenes Direction    PValue       FDR PValue.Mixed
Glycolysis_/_Gluconeogenesis     43      Down 0.6824303 0.9613343    0.9996904
Citrate_cycle_(TCA_cycle)        27        Up 0.9613343 0.9613343    1.0000000
                             FDR.Mixed
Glycolysis_/_Gluconeogenesis         1
Citrate_cycle_(TCA_cycle)            1

Are the mixed tests meant to be comparable between fry and mroast? I know that the mixed test for mroast just uses the absolute value of the t-statistics, but for fry it's an F-test, so in that context I don't really think they should be the same, but the help page doesn't to my knowledge distinguish between the two.

limma • 539 views
ADD COMMENT
1
Entering edit mode
@gordon-smyth
Last seen 46 minutes ago
WEHI, Melbourne, Australia

That's a fair question and you're right that the fry/mroast help page didn't give enough information about the mixed p-values.

Your suspicion is correct. The exact correpondence between fry and roast is only for the directional p-values, not for the mixed p-values. When the genewise variances are the same (df.prior is large) then the fry directional p-value agrees with what roast would give with set.stat="mean" and a very large number of rotations. However, the mathematics that makes that possible only works for the directional p-values. The fry mixed p-values therefore have to be computed using a different method, making some different approximations, and they do not converge to the corresponding roast mixed p-values. I had hoped to publish all the fry details in a journal paper years ago but other projects keep competing for our time.

Anyway, I have now revised the fry/mroast help page to try to clarify this.

ADD COMMENT

Login before adding your answer.

Traffic: 635 users visited in the last hour
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