5

I am working with a set of genes for which I have both methylation values $\beta$ (continuous on unit interval) and gene expression $E$ (non-negative continuous) that I want to test for correlation $\hat\rho(log_2E,\beta)=R$. To filter out only the significant correlations I perform a permutation test where for each site I randomly permute $E$ a large number of times and calculate $R$. This gives me a distribution under the null hypothesis $p(R|H_0)$, i.e. when there is no dependence between $\beta$ and $E$, that I can compare $R$ against.

Now to my question,

I want to test both for positive and negative correlation. Do I make two tests per gene and use double multiple testing correction, where

$p_- = \int\limits_{-\inf}^Rp(x|H_0)dx$

$p_+ = \int\limits_{R}^\inf(x|H_0)dx$

or can and should I calculate a two-tailed p-value directly based on $|R|$ or such, and use normal multiple testing correction? I can only think of the following, but it feels wrong and has lousy power.

$p_2 = \int\limits_{-\inf}^{-|R|}p(x|H_0)dx + \int\limits_{|R|}^\inf(x|H_0)dx$

enter image description here

Edit: Updated the figure and added the missing x-axis label. The y-axis is the same for all panels.

Backlin
  • 213
  • 1
    What is your Null Hypothesis? That should determine how you derive the $p$ value from the permutation test. If the Null Hypothesis is that the correlation is 0 and that large positive and negative correlations are evidence against this Null, then you want the blue areas. – Gavin Simpson Aug 10 '12 at 12:59
  • Yes, the null hypothesis would be that the correlation is 0. I just find it hard to understand why I get so much fewer hits than when combining the red and green tests. – Backlin Aug 10 '12 at 13:07
  • 2
    The green isn't really a test, it is the region that is consistent with the Null hypothesis in the one-sided test of 0 correlation where the alternative is that the correlation is negative. Actually, I wonder if your blues are too big? If you want a 95% two-sided test, the quantiles you want at the 0.25th and 0.975th such that the sum of the blue areas is 0.05 (5%). If a one-sided test, then all the evidence against the null is on one tail, in the two-sided evidence is in both tails. – Gavin Simpson Aug 10 '12 at 13:17
  • Thanks for taking your time with my question! What if I had tested the null hypothesis that $R>0$ and gotten $R=-0.2$? It would result in the green area above, but you are saying that it would disqualify the test since $R$ directly defies $H_0$? I still don't get whether or not I am allowed to do that test, given that I include it when correcting for multiple testing.

    The data underlying the plot is chosen arbitraryily to illustrate the equations. The blue area is 0.147 and would correspond to the $p$ value of $R=\pm 0.2$ under $H_0: R=0$ for this particular gene, wouldn't it?

     – Backlin Aug 10 '12 at 13:33
  • 1
    You can't test for $R > 0$ here as otherwise you'd need to state something about how big $R$ is. Obviously the real $R$ is unknown otherwise you wouldn't be doing the test. The Null hypothesis is $R = 0$ and you are looking for evidence against this Null. I may have misunderstood what the plots show? Usually these plots would show the permutation distribution of $R$. – Gavin Simpson Aug 10 '12 at 13:41
  • When you say "I want to test both for positive and negative correlation. Do I make two tests per gene", do you mean that you want to carry out two one-sided test, where $H_0:R=0$ and in the first test $H_{1,A}:R>0$ and in the second test $H_{1,B}:R<0$, each at an $\alpha/2$ level? Because in that case it would be exactly the same as carrying out one two-sided test where $H_1:R\ne 0$ at an $\alpha$ level, since it is impossible for both hypotheses to be rejected. – boscovich Aug 10 '12 at 13:48
  • I agree with Gavin. If the curve is the permutation distribution for R under the null hypothesis the blue area represents the region for the two-tailed test. If you want to do a one-side test where the alternative is R>0 then a right side tail with the appropriate area for the given significance level is what you would use. That would not be the green region that you exhibit. – Michael R. Chernick Aug 10 '12 at 13:51
  • @andrea One two-sided test is of course possible. The reject region looks like the figure with the two blue tail regions. The idea is not to reject for any estimated value different from zero but rather to reject when the estimate is either much greater than 0 or much less than zero. – Michael R. Chernick Aug 10 '12 at 13:55
  • If you are applying these tests many times to test several different genes than you must adjust your results to account for multiple testing. – Michael R. Chernick Aug 10 '12 at 13:56
  • @MichaelChernick: yes, I was only pointing out that two one-sided tests (each at an $\alpha/2$ level), to see whether you have evidence to conclude that $R$ is either larger or smaller than $0$, are exactly equivalent as one two-sided test (at an $\alpha$ level). – boscovich Aug 10 '12 at 14:03
  • 1
    @andrea Okay that is true it is just a different way of looking at the same thing. – Michael R. Chernick Aug 10 '12 at 14:12
  • I am starting to see the gist of the problem now. You all agree the two-sided test (blue) is valid and that two one-sided tests should be equivalent to one two-sided. How would you go about doing a one-sided permutation test then, and what happens if the observed R falls on the $H_0$ side vs. the $H_1$ side? – Backlin Aug 10 '12 at 14:32

2 Answers2

9

Here is an R example of permutation testing.

## dummy data
set.seed(1)
x <- runif(20)
y <- 0.5 * x
y <- y + rnorm(20)

## set up for the permutation, compute observed R
nullR <- numeric(length = 1000)
nullR[1] <- cor(x, y) ## obsered R in [1]
N <- length(x)

## permutation test
for(i in seq_len(999) + 1) {
    nullR[i] <- cor(x[sample(N)], y)
}

hist(nullR) ## histogram of R under H0

enter image description here

Now we can compute the permutation $p$ from the permutation distribution for the various tailed tests you consider:

> ## one side H1 R > 0
> sum(nullR >= nullR[1]) / length(nullR)
[1] 0.908
> ## one side H1 R < 0
> sum(nullR <= nullR[1]) / length(nullR)
[1] 0.093
> ## two sided
> sum(abs(nullR) >= abs(nullR[1])) / length(nullR)
[1] 0.177
Gavin Simpson
  • 47,626
  • I fail to see how this is any different than my examples, except that I used nullR to estimate the density $p(R|H_0)$ and then integrated rather than sum. For example isn't your first example (where $H_1: R>0$) more or less the same thing as my green case? – Backlin Aug 10 '12 at 18:40
  • Again, we might be talking at cross-purposes regarding the actual plots, but the way to look at this is via the blue plot. For R to be statistically significant at a given level, in the one tailed test with H1 R > 0 the observed R would need to be in the right hand blue area. For the H1 R < 0, the observed would have to be in the left hand blue area. For the two tailed test with H0 R = 0, H1 R!= 0, then the observed has to be in either of the blue areas. – Gavin Simpson Aug 10 '12 at 18:52
  • @Backlin Your plots are confusing as the upper and lower panels show the regions for rejecting the Null in a one tailed test where H1 R < 0 (upper, red area) and rejecting the Null of a two tailed test where H1 != 0 (lower, blue areas). The green area in the middle plot is not the rejection area for the one-tailed test with H1 R > 0. I think this is what might be the cause of confusion between people here and yourself. In the one side H1 R > 0 we clearly have no evidence against H0 because the observed R lies in the green area. – Gavin Simpson Aug 10 '12 at 18:55
2

Your formula for $p_2$ only works if you have distribution symmetric around $0$, which at least your illustration do not seem to have. A better estimate is given by $p_2(t) = {\rm min}(2p_{+}, 2p_{-})$. using your definitions of $p_{+}$ and $p_{-}$.

Lukas Käll
  • 49