Difference in mixed test for fry vs mroast
1
0
Entering edit mode
@james-w-macdonald-5106
Last seen 8 minutes 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 • 116 views
1
Entering edit mode
@gordon-smyth
Last seen 27 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.