Testing for equality of variance in permutations (microarray analysis with bioconductor)

Question

I have measured whole-genome gene expression in two groups of animals, n=6 in each group. My goal is to detect differentially expressed genes - pretty standard analysis. The typical thing to do, and what I suspect is a situation where 'everybody does it this way because everybody does it this way' is to run limma with moderated t test. However, moderated t test assumes equal variance between the groups. I think most people whose papers I have seen just ignore this assumption (or - more correctly - do not mention this consideration in their methods at all).

One can run Wilcoxon/Mann U test to get rid of the requirement for equal variance, but for genes that have equal variance (and in my case this is >80% of all genes), this test is much less powerful. One can run a Welch t test that does not assume anything about variance, but again it is weak and one looses a lot of genes. Plus I want to use moderated t test because it's better than other things (see e.g., Linear Models and Empirical Bayes Methods for Assessing Differential Expression in Microarray Experiments, Should We Abandon the t-Test in the Analysis of Gene Expression Microarray Data: A Comparison of Variance Modeling Strategies; although moderated Welch seems to be OK: Unequal group variances in microarray data analyses). So what I ended up doing is:

1. Run levene's test for equality of variance for each of the genes (package lawstat).
2. For genes that have equal variance, run limma with moderated t test.
3. For genes with unequal variance, run Mann U.

But this causes problems for me when calculating FDR. I calculate FDR by permuting the group assignment of each sample and re-running the 3-point pipeline above. The principal problem is computational power - it takes about 30 hours on a 3Ghz 8Gb RAM Ubuntu machine for levene's test to run the 1000 permutations for the ~10000 genes (and I have 6 different configurations to run in total). The bigger problem is, quite frustratingly, that I am not sure this is the most appropriate way to go...

If you could advise me on the best practices and practical solution to this situation, I will be extremely grateful. Even if you just confirm what I do is correct... :-)

Answer 1

It is true that most of the time 'everybody does it this way because everybody does it this way' (+1)

Some time ago, I realized that the t test (as implemented in R) is relatively robust to inequality of variances. You can run the following test to convince yourself.

pvals <- rep(NA, 10000)
for (i in 1:10000) {
    pvals[i] <- t.test(rnorm(10), rnorm(10, sd=50))$p.value
}
plot(sort(pvals), type='l')

I also realized recently that the power difference between Wilcoxon and t tests is not that big as suggested here, at least for Gaussian data.

So in my opinion, sticking to one test is OK and will make your life simpler to compute FDR. I would use Wilcoxon because it is more robust to outliers.

Two more things.

1/ In my humble experience, it is better to look for robust signal in microarray experiments because you usually have to do follow up experiments that have to be reproducible. Strong power allows you to detect weak effects, but weak effects are a pain. Unless I know exactly why I need power, I use non parametric tests.

2/ That the Wilcoxon test is robust to inequality of variances does not mean that it does not assume equality of variances. In general, it is false to assume that the Wilcoxon test is always applicable. The Wilcoxon test is sometimes said to test for equality of median, here is a proof that this is not what it does.

set.seed(123)
x <- rgamma(1000, 1)
x <- x - median(x)    # x has median 0
y <- -rgamma(1000, 1)
y <- y - median(y)    # y has median 0
wilcox.test(x,y)      # Surprise!!

The Wilcoxon test actually tests that one distribution is shifted relative to the other (so it assumes that they have the same variance).