Tests of Randomness and Independence

References:

The null hypothesis; models

Suppose a treatment has N ordered levels, and one response is observed per treatment level, with N subjects assigned at random, one to each treatment level. Let Tj denote the rank among the responses of the response of the subject assigned to the treatment level with rank j. Then we have a set of pairs of ranks of treatments and responses:

{(j, Tj): j= 1, 2, … N}

There are N! possible sets of pairs, depending on which treatment rank goes with with response rank. If treatment has no effect, all N! pairings are equally likely; each has probability 1/N!. This is the "Hypothesis of Randomness."

Now consider drawing a random sample of N subjects from a population, and assigning one to each of the N treatment levels. Let Zj denote the response of the subject assigned treatment j.

P(Zjz) = Fj(z), j= 1, 2, … N.

If treatment has no effect,

F1 = F2 = … = FN.

If the size of the population is much much larger than N, we can ignore the dependence among the subjects. That leads to the population model, in which

{Zj: j= 1, 2, …  N}

are independent. Under the null hypothesis, they are thus iid, and so the null distribution of the pairs {(j, Tj)} is the same as it is for the hypothesis of randomness.

Consider a third statistical model: we make two measurements (Xj, Yj) on each of a set of N subjects. Under the null hypothesis, the Xs are identically distributed, the Ys are identically distributed, and the measurements are independent. This also leads to the same joint distribution for the pairs {(j, Tj)}, where Tj is the rank of the Y measurement corresponding to the X measurement with rank j. The key here is that one set of measurements is exchangeable given the other.

Tests for Trend

We are going to examine tests for a trend: the alternative hypothesis is that increasing levels of treatment are associated with increasing (or decreasing) levels of response.

The Spearman rank correlation & related statistics

The first statistic we consider counts the pairs (Ti,Tj) with i<j for which Ti<Tj. Let

Uij ≡ {1, Ti<Tj; 0, Ti≥Tj}

and

B ≡ ∑i<j Uij.

Then B tends to be large when larger responses are associated with larger treatment indices j. A more common test statistic is D':

D' ≡ ∑i<j (j−i) Uij.

This puts more weight on larger differences (ji). Now,

D' = (1/6) N(N2 − 1) − ½ ∑i=1N (Ti − i)2.

Thus a test based on D' would reject when

D ≡ ∑i=1N (Ti − i)2

is small. Note that D takes only even values:

D = ∑i=1N (Ti2 − 2iTi + i2)

= 2 [ ∑i=1N i2 − ∑i=1NiTi]

= [N(N+1)(2N+1)/3] − 2∑i=1N(iTi).

Thus to reject when D is small is to reject when i=1NiTi is big. Equivalently, it is to reject when i=1N(RiTi) is big, where Ri is the rank of the ith treatment and Si is the rank of the ith response. Let

r ≡ (1/N) ∑i=1N Ri = s ≡ (1/N) ∑i=1N Si = (N+1)/2,

so

(1/N) ∑i(Ri − r)2 = (1/N) ∑i(Si − s)2 = (1/N) (N3−N)/12 = (N2−1)/12.

Spearman's rank correlation is

rS[(1/N) ∑i (Ri − r)(Si − s)]/[ √[(1/n)∑(Ri−r)2] √[(1/n)∑(Si−s)2] ]

= ∑[RiSi − r Si − Ris + rs ]/[(N3−N)/12]

= [12/(N3−N)] [ ∑ RiTi − N(N+1)2/4 ].

A bit more algebra shows that

rS = 1 − 6D/(N3 − N).

If there are ties, we can define analogous test statistics using midranks instead of ranks. Note that if the null hypothesis is true and there are no ties, ErS = 0, and

rS √[(N−2)/(1−rS2)] → tN−2,

where the limit is in distribution. This approximation is pretty good for N≥10 or so. Alternatively, one can simulate critical values for a test based on Spearman's rS by computing the usual correlation coefficient between {1, … N} and random permutations of {1, … N}.

For a bivariate normal distribution of (X,Y), the Pitman efficiency of rS relative to the usual Pearson correlation coefficient is (3/π)2 ≅ 0.912.

If there is serial correlation, the null distribution is quite different.

One alternative to Spearman's rank correlation coefficient is to use Pearson's correlation coefficient, but to compute critical values by simulation using random permutations of the responses to treatment {Yj} relative to the treatments {Xj}.

The run test

For many alternatives, high values tend to cluster and low values tend to cluster, leading to "runs." We can use this to test the hypothesis of randomness by looking at the lengths and numbers of runs. For example, suppose we toss a (possibly biased) coin N times. We can think of this as N trials. We consider the outcome of the jth trial to be 1 if the coin lands heads on the jth toss and 0 otherwise. Under the null hypothesis, the N outcomes are iid Bernoulli(p) random variables, where P is the chance that the coin lands heads.

Let R denote the number of runs. For example, the sequence HHTHTTTHTH has 7 runs: HH, T, H, TTT, H, T and H. Condition on the number n of heads among the tosses. If the null hypothesis is true, each arrangement of the n heads among the N tosses has probability 1/(NCn). We will compute the (conditional) probability distribution of R given n, under the null hypothesis of independence.

Clearly, if n=N or if n = 0, R ≡ 1. If k = 1, there are only two possibilities: first all the heads, then all the tails, or first all the tails, then all the heads. I.e., the sequence is either

(HH … HTT … T) or (TT … THH … H).

where there are n Hs and m ≡ N−n Ts in all. The probability that R=1 is thus 2/(NCn), if the null hypothesis is true.

How can R be even, i.e., R = 2k? If the sequence starts with a head, we need to choose where to break the sequence of heads to insert tails, then where to break that sequence of tails to insert heads, etc. If the sequence starts with a tail, we need to choose where to break the sequence of tails to insert heads, then where to break that sequence of heads to insert tails, etc. We need to break the n heads into k groups, which means picking k − 1 breakpoints, but the first breakpoint needs to come after the first H, and the last breakpoint needs to come before the nth H, so there are only n − 1 places those k − 1 breakpoints can be. And we need to break the m tails into k groups, which means picking k − 1 breakpoints, but the first needs to be after the first T and the last needs to be before the mth T, so there are only m − 1 places those k − 1 breakpoints can be. The number of sequences with R = 2k that start with H is thus

n−1Ck−1 × m−1Ck−1.

The number of sequences with R = 2k that start with T is the same (just read right-to-left instead of left-to-right). Thus, if the null hypothesis is true,

P(R = 2k) = 2 × n−1Ck−1 × m−1Ck−1/NCn.

Now consider how we can have R = 2k+1. Either the sequence starts and ends with H or it starts and ends with T. Suppose it starts with H. Then we need to break the string of n heads in k places to form k + 1 groups using k groups of tails formed by breaking the m tails in k−1 places. If the sequence starts with T, we need to break the m tails in k places to form k + 1 groups using k groups of heads formed by breaking the n heads in k−1 places. Thus, under the null hypothesis,

P(R = 2k+1) = [ n−1Ck × m−1Ck−1 + n−1Ck−1 × m−1Ck ]/NCn

Note that nothing in this derivation used the probability of heads. The conditional distribution under the null hypothesis depends only on the fact that the tosses are iid, not that the coin is fair.

Let Ij be the indicator of the event that the outcome of the j+1st toss differs from the outcome of the jth toss, j=1, …, N−1. Then

R = 1+ ∑j=1N−1 Ij.

Under the null, conditional on n,

P(Ij = 1) = P(Ij = 1 | jth toss lands H, n)P(jth toss lands H | n) + P(Ij = 1 | jth toss lands T, n)P(jth toss lands T | n)

= P(j+1st toss lands T | jth toss lands H, n)P(jth toss lands H | n) + P(j+1st toss lands H | jth toss lands T, n)P(jth toss lands T | n)

= [m/(N−1)]×[n/N] + [n/(N−1)]×[m/N] = 2nm/[N(N−1)].

The indicators {Ij} are identically distributed under the null hypothesis, so if the null holds,

ER = E[ 1 + ∑j=1N−1 Ij] = 1 + (N-1) × 2nm/[N(N−1)] = 1 + 2nm/N.

Air temperature is measured at noon in a climate-controlled room for 20 days in a row. We want to test the null hypothesis that temperatures on different days are independent and identically distributed.

Let Tj be the temperature on day j, j = 1, …, 20. If the measurements were iid, whether each day's temperature is above or below a given temperature t is like a toss of a possibly biased coin, with tosses on different days independent of each other. We could consider a temperature above t to be a head and a temperature below t to be a tail.

Let's take t to be the median of the 20 measurements. In this example, n=10, m=10, N=20. We will suppose that there are no ties among the measured temperatures. Under the null hypothesis, the expected number of runs is

ER = 1+2mn/N = 11.

The minimum possible number of runs is 2 and the maximum is 20. Since we expect temperature on successive days to have positive serial correlation (think about it!), we might expect to see fewer runs than we would if temperatures on different days were independent. So, let's do a one-sided test that rejects if there are too few runs. We will aim for a test at significance level 5%.

P(R = 2) = 2/20C10 = 1.082509e-05.

P(R = 3) = 2×9C1×9C0/20C10 = 9.74258e-05.

P(R = 4) = 2×9C1×9C1/20C10 = 8.768321e-04.

P(R = 5) = 2×9C2×9C1/20C10 = 0.003507329.

P(R = 6) = 2×9C2×9C2/20C10 = 0.01402931.

P(R = 7) = 2×9C3×9C2/20C10 = 0.03273507.

P(R ≤ 6) = 2×(2+9+81+324+1296)/20C10 ≈ 0.0185.

So, we should reject the null hypothesis if R ≤ 6, which gives a significance level of 1.9%. Including 7 in the rejection region would make the significance level slightly too big: 5.1%.

When N, n and m are large, the combinatorics can be difficult to evaluate numerically. There are at least two options: asymptotic approximation and simulation. There is a normal approximation to the null distribution of R. As n and m→∞ and m/n→γ,

[R − 2m/(1+γ)]/√[4γm/(1+γ)3] → N(0, 1)

in distribution.

Here is an R function to simulate the null distribution of the number R of runs, and evaluate the P-value of the observed value of R conditional on n, for a one-sided test against the alternative that the distribution produces fewer runs than independent trials would tend to produce. The input is a vector of length N; each element is equal to either "1" (for heads) or "-1" (for tails). The test statistic is calculated by finding 1 + ∑j=1N−1 Ij, as we did above in finding ER.

        simRunTest <- function(x, iter) {
            N <- length(x);
            ts <- 1 + sum(x != c(x[2:N],x[N]));   # test statistic: count
                                                  #  transitions I_j.
            sum(replicate(iter, {
                                  xp <- sample(x);
                                  ((1 + sum(xp != c(xp[2:N],xp[N]))) <= ts)
                                }
                         )
               )/iter
        }
    

As an example, suppose the observed sequence is x = (-1, -1, 1, 1, 1, -1, 1), for which N = 7, n = 4, m = 3 and R = 4. In one trial with iter = 10,000, the simulated P-value using simRunTest was 0.5449. Exact calculation gives:

P0(R ≤ 4) = (2 + 5 + 12)/35 = 19/35 ≈ 0.5429.

The standard error of the estimated probability is thus

SE = √[(19/35×16/35)/10000] ≈0.005.

The simulation was off by about (0.5449 − 0.5429)/0.005 ≈ 0.41SE.

Tests for independence between time series

The crucial element of the null hypothesis that leads to the null distribution of the Spearman rank correlation is that one of the sets of measurements is conditionally exchangeable given the other. Then, if the null hypothesis is true, all pairings (Xi, Yj) are equally likely. However, time series data often have serial correlation, such as trends. Suppose we have two independent time series, each of which has a trend. For example, we might observe the pairs (Xj,Yj)j=1N, where

Xj = j + εj, j=1, … N,

Yj = j + νj, j=1, … N,

and the 2N "noise" terms j} and j} are iid with zero mean and finite variance. Even though the time series {Xj} and {Yj} are independent, neither {Xj} nor {Yj} are exchangeable. As a result, the Spearman rank correlation test will tend to reject the hypothesis of independence, not because the two sets of measurements are dependent, but because a different feature of the null hypothesis is false. The trend makes pairings (Xi, Yj) with both X and Y larger than average or both smaller than average more likely than pairings where one is relatively large and the other is relatively small.

Walther's Examples

The following examples are from Walther (1997, 1999). Suppose that {Xj}j=1100 and {Yj}j=1100 are iid N(0,1). Define

Sk ≡ ∑{j=1}k Xj and Tk ≡ ∑{j=1}k Yj.

Then

P(rS(S, T) > c0.01) ≈ 0.67,

where c0.01 is the critical value for a one-sided level 0.01 test against the alternative of positive association. That is, even though the two series are independent, the probability that the Spearman rank correlation coefficient exceeds the 0.01 critical value for the test is over 2/3. That is because the two series S and T each have serial correlation, so not all pairings (Si, Tj) are equally likely—even though the two series are independent.

Serial correlation is not the only way that exchangeability can break down. For example, if the mean or the noise level varies with time, that violates the null hypothesis. Here is an example of the latter. Take X = (1, 2, 3, 4) to be fixed. Let (Y1, Y2, Y3, Y4) be independent, jointly Gaussian with zero mean, SD(Yj) = 1, j = 1, 2, 3, and Var(Y4) = 4. Then

P0(rS(X, Y) = 1) = 1/4! = 1/24 ≈ 4.17%.

Note that in this example, rS = 1 whenever Y1<Y2<Y3<Y4. Simulation shows that in fact P(rS(X, Y)=1) is about 7%:

        iter <- 10000;
        s <- c(1,1,1,4);
        sum(replicate(iter, {
                             x <- rnorm(length(s), sd=s);
                             !is.unsorted(x)  # r_S=1 if x is ordered
                            }
                     )
           )/iter
    

In a similar vein, take X = (1, 2, 3, 4, 5) to be fixed, let (Y1, Y2, Y3, Y4, Y5) be independent, jointly Gaussian with zero mean and SDs 1, 1, 1, 3, and 5, respectively. Then, under the (false) null hypothesis that every (X, Y) pairing is equally likely,

P0{rS(X, Y) = 1} = 1/5! ≈ 0.83%,

but simulation shows that the actual probability is about 2.1%. In these examples, the null hypothesis is false, but not because X and Y are dependent. It is false because not all pairings (Xi,Yj) are equally likely. It is the "identically distributed" part of the null hypothesis that fails.

The solar neutrino problem

Sun emits neutrinos as a byproduct of nuclear fusion in the solar core. Nuclear theory predicted a given neutrino flux; observations showed a rather lower flux than predicted. Why? Scientists speculated about mechanisms for many years. There is an apparent negative association between sunspot number (a measure of the activity of magnetic fields in Sun) and neutrino flux. As of about 2002, the Homestake experiment had collected about 30 years of data on the solar neutrino flux by measuring daily production of 37Ar production (atoms per day). (They published upper and lower 68% confidence limits.)

In the early 2000's, physicists found that neutrinos have mass. Neutrino mass could explain the apparent deficit of solar neutrinos through state oscillations.

References

Walther, G., 1997. Absence of correlation between the solar neutrino flux and the sunspot number, Phys. Rev. Lett., 79, 4522–4524.

Walther, G., 1999. On the solar-cycle modulation of the Homestake solar neutrino capture rate and the shuffle test, Astrophys. J., 513, 990–996.