lawstat runs.test valid for small samples?

Question

This is both an R programming question and a statistics question. From my experiments, it appears that the runs.test function in the R package lawstat gives very strange results for small samples. Can anyone confirm, refute, and/or explain? My reasoning follows below.

My test data is a count of patents issued to a firm in one technology class in each of 15 years.

testpats <- c(2,1,2,0,1,4,1,1,2,4,2,6,1,3,3)

Running

runs.test(testpats, plot.it=T, alternative="positive.correlated")

first of all, produces the following picture of the runs. (won't let me post images so here is my recreation.)

B B B B A B B B A A A B A A

according to the documentation "observations which are less than the sample median are represented by letter "A", and observations which are greater or equal to the sample median are represented by letter "B"."

The sample median of testpats is 2. So if the documentation were correct, the image should look like:

= - = - - + - - = + = + - + +
B A B A A B A A B B B B A B B

Clearly this is very different, so I have no idea what runs.test is using for the "sample median".

Second, the test statistic given by the function output

        Runs Test - Positive Correlated

data:  testpats 
Standardized Runs Statistic = -0.4877, p-value = 0.3129

Is very different from what I would calculate by hand using the methods described at https://www.itl.nist.gov/div898/handbook/eda/section3/eda35d.htm

mymid <- median(testpats)
runsdummy <- ifelse(testpats >= mymid, 1, -1)
n1 <- length(which(runsdummy>0)) #number of values above or equal to the median
n2 <- length(which(runsdummy<0)) #number of values below the median
sr2 <- (2*n1*n2*(2*n1*n2 - n1 - n2))/((n1+n2)^2 * (n1+n2-1)) #standard deviation of the number of runs
Rbar <- (2*n1*n2)/(n1+n2) + 1 #expected number of runs
R <- 9 #observed number of runs - how do I automate?
Z <- (R-Rbar)/sr2 #runs test statistics
Z

gives

[1] 0.2508961

Note that this hand calculated test statistic bears no resemblance to the -0.4877 provided by runs.test().

Alternately, I could use the small sample version of the test explained in Swed and Eisenhart. The small sample method just uses the number of above and below observations and the number of runs.

Given n1 = 5; n2 = 6; R = 9

the one-sided pvalue should be 0.976.

Again, this is not even close to the number produced by runs.test()

So, what gives? Do I completely misunderstand how to use runs.test()? I've tried using the function after transforming the data into above/below indicators (e.g. 1/-1), and I still get strange results.

Answer 1

I was stumbling upon the same problem, while modelling this in Excel and comparing to output of my StatGraphics software. I finally found my 'solution' in the documentation of StatGraphics. I noted it in R format (don't use R anymore but I think this is right) :

Z <- (R-0.5-Rbar)/sr2

I don't yet know why the 0.5 has to be subtracted (or added in some cases), but I think it has something to do with one-sided and two-sided testing. +0.5 or -0.5 would then test one-sided (higher or lower), and without this addition I think is two-sided.

Don't know if I'm right yet, but I got the same results as my StatGraphics model using -0.5.

Try it out and let me know!

Edit (from software documentation): Calculate the probability of observing at least k runs: use -0.5 Calculate the probability of observing less than or equal to k runs: use +0.5

Edit 2: The + or - 0.5 is a continuity correction. You can only observe i.e. 3 or 4 occurrences, not something in between. If you calculate the chance at 3 as (chance at 3.5 or less) and the chance at 4 as (chance at more than 3.5), only then will the summed chances be 1.

Answer 2

Two things. sr2 is the variance: Take the root. n1 + n2 == length(testpats) != 5+6. I got 6 & 9.