Best way to test runs

Question

I have the following sample:

x <- c(rep(1, 8), 0, 1, 0, rep(1, 5), 0, 0, 1, rep(0, 11), 1, 0, 0, 1, 0)
> x
[1] 1 1 1 1 1 1 1 1 0 1 0 1 1 1 1 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0

The 1 and 0 are actually factors of a nominal variable and the sample is sorted according to a hidden variable. It looks like the sequence is not random, that the 1s tend to rank lower (i.e. appear earlier) in this sequence. I would like, however, to conduct a statistical test to test this hypothesis.

So far, my solution was to use the Wald-Wolfowitz test, like so:

library(adehabitat)
wawotest(x)

It yields a p-value of 0.0211, which I consider low enough to reject the hypothesis that the ones and zeroes occur randomly in the given run.

However, I've read that the Kolmogorov-Smirnov test is a more powerful solution and I'd like to try it in R, but I'm unsure I'm giving it the right command. This is what I'm running:

ks.test(x, 'punif')

Which gives a p-value of 1,821e-08. Is my command correct? Should I consider any other tests?

[EDIT] As @Glen_b correctly put it, the KS is not suited for this case, so I'll be sticking with the WW. Here are the two best ways I've found for running this test:

> library(tseries);    runs.test(as.factor(x))
> library(adehabitat); wawotest(x)

I've studied the algorithm behind runs.test() and it seems to be conducting the Wald-Wolfowitz, even though the p-values are slightly different.

Anyway, WW doesn't look like the most powerful tool out there for this kind of test, so feel free to provide other solutions.

Answer 1

Your stated objective is to assess whether

the 1s tend to rank lower (i.e. appear earlier)

That is not measured by runs, but by ranks. Use the Wilcoxon (aka Mann-Whitney) test.

---

This test is both conceptually and computationally simple, yet reasonably powerful. The data are ranked in the order of appearance using the numbers $1, 2, \ldots, n$ (where $n=35$ in this case). The ranks are summed within each group: there's a sum of $n_0=18$ ranks corresponding to the zeros and a sum of $n_1=17$ ranks corresponding to the ones. To compensate for the different numbers of zeros and ones, subtract the smallest possible sum from each (equal to $1+2+\cdots+n_i = \binom{n_i+1}{2}$ for group $i$, $i=0,1$). If the ones truly tend to come first, their adjusted rank sum will be substantially smaller than that of the zeros. This can be converted into a Z-score by assuming an asymptotic Normal distribution for the statistic or a more accurate p-value can be found through the permutation distribution. The code below illustrates both methods.

For these data, the zeros appear at ranks

9 11 17 18 20 21 22 23 24 25 26 27 28 29 30 32 33 35

while the ones appear at ranks

1  2  3  4  5  6  7  8 10 12 13 14 15 16 19 31 34

The adjusted sum of the ranks of the ones is $U=47$. The Normal approximation estimates its p-value at $0.0002339$. The small value is testimony to the power of this test in the present case. The permutation distribution, estimated with a million replications, gives a p-value of $0.000264$. It is accurate to $\pm 0.000016$ (which is one standard error). Either p-value gives you ample basis to reject the null hypothesis that the zeros and ones are randomly scattered throughout the sequence.

Here is a histogram of the permutation distribution of the $U$ statistic for these data.

The red vertical line marks the actual test statistic. It obviously is extreme.

Although it might not look it, this test was conducted as a two-tailed test (by taking the smaller of the two adjusted rank sums). It tests whether there is any difference in ranks, not just whether the ones tend to come earlier.

---

Below is the (reproducible) R code that made the figure and computed the p-values. This large simulation takes about ten seconds to run. (Only a few thousand replications are typically needed. A million were used to make the point that the Normal approximation works well in this case.) Reduce the first argument of replicate at line 25 in order to achieve a faster run time.

x <- c(1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 0, 
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0)
#
# Wilcoxon test.
#
Wilcoxon <- function(x) {
  n <- length(x)
  n0 <- sum(x==0)
  n1 <- sum(x==1)
  u0 <- sum((1:n)[x==0]) - choose(n0+1, 2)
  u1 <- sum((1:n)[x==1]) - choose(n1+1, 2)
  u <- min(u0, u1)

  m.u <- n0 * n1 / 2
  s.u <- sqrt(n0 * n1 * (n+1) / 12)
  Z <- (u - m.u)/s.u
  p <- pnorm(Z)
  return(c(U=u, Z=Z, p.value=p))
}
stats <- Wilcoxon(x)
#
# Permutation test.
#
set.seed(17)
U <- replicate(1e6, Wilcoxon(sample(x, length(x)))["U"])
hist(U, main="Permutation Distribution", )
abline(v = stats["U"], lwd=2, col="Red")
#
# Summary.
#
message("Normal approximation: ", signif(stats["p.value"], 4), 
        "  Permutation estimate: ", signif(mean(c(1, U <= stats["U"])), 4),
        " +/- ", signif(sd(c(1, U <= stats["U"])) / sqrt(1 + length(U)), 2))

Answer 2

I think you have this muddled up.

(i) The Kolmogorov Smirnov test is designed to test continuous, not discrete distributions; indeed, your values consist only of 0's and 1's, yet it appears you're testing against continuous uniformity.

(ii) as written, it looks like this is completely ignoring the time order. It's not testing what you need to pick up.

You could perhaps use as your data the fraction (considered as a quantile) of the distance through the series that the 1's occur at (it's not continuous independent data though - it only occurs at discrete places, only one value per location -- so you'd still need to adjust your null distribution for that). It wouldn't be a KS test as such, but you could use a statistic like the KS as the basis for a test.

For example, if there are $n$ observations, the $i$-th observation might be said to occur at the $\frac{i-\alpha}{n+1-2\alpha}$ quantile for some $0\leq\alpha\leq 1$ (I believe many of the 9 alternatives in R's own quantile function correspond to that definition with various values of $\alpha$). You can then test whether the quantiles of the 1's are uniformly distributed, but you'd need to simulate to get the distribution of the test statistic under the null.

An easy alternative to simulation of the null distribution (presumably conditioning on the counts of 0s and 1's) would be to do a permutation test. (Which will involve either clever algorithms to do do the full distribution, or sampling of the permutation distribution.

However, it seems as if you're really after a test for trend. Indeed, you might do better with something as simple as a logistic regression against position, or even a monotonic GAM-type model (again, probably via logistic-regression).

---

Edit: Here's the previously suggested logistic regression performed in R:

x <- c(rep(1, 8), 0, 1, 0, rep(1, 5), 0, 0, 1, rep(0, 11), 1, 0, 0, 1, 0)
t <- seq_along(x)                       # Rank order by position (1,2,3,...)
plot(x~t)                               # Show the sequence of 1's and 0's
logistfit <- glm(x~t,family=binomial)   # fit a straight line in the logits
summary(logistfit)                      # show GLM regression table output
f <- fitted(logistfit)                  # fit is estimated P(X=1|t)
lines(f~t,col=4)                        # plot that fit

Here's the output of the model (a few less interesting lines removed):

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  2.77412    1.01415   2.735  0.00623 
t           -0.15912    0.05225  -3.046  0.00232     # <=== the line we want

    Null deviance: 48.492  on 34  degrees of freedom
Residual deviance: 34.087  on 33  degrees of freedom
AIC: 38.087

The p-value for the glm fit is $0.00232$. It shows that the probability of a 1 is not consistent with the 1's being randomly placed with respect to the ordering variable. Since the coefficient is negative, the probability of a 1 is overall decreasing as the position increases.

Here's the plot:

Answer 3

You can use also a test based on the Length of the longest head run (LLHR). In your example LLHR is too large -> you should reject the hypothesis that the sample is from Bernoulli B(1/2) distribution.