# Where do Wilcoxon signed rank critical values come from (ie how are they computed)?

Here is the table for Wilcoxon signed rank test critical values, and obviously we can decide whether or not to reject the null hypothesis for a particular scenario based on this table. I understand and am fine with that, but I am curious about how these critical values are computed in the first place. For t test for example, we use t distribution to identify critical values, but how about for a nonparametric test like Wilcoxon?

The short answer is basically "fancy counting".

In very small sample sizes, we can literally count possibilities and construct the distribution.

Consider n=3. The possible signs are as follows:

Rank:    1     2     3      Sum(R+)   Prob
Signs:   -     -     -        0        1/8
+     -     -        1        1/8
-     +     -        2        1/8
-     -     +        3        1/8 \__ same outcome P=2/8
+     +     -        3        1/8 /
+     -     +        4        1/8
-     +     +        5        1/8
+     +     +        6        1/8


This is fine out to about n=5 or 6. Note that the distribution will be symmetric, so at worst we need to calculate about half of it. Indeed, since we don't need the whole distribution for a table of critical values (just the tail), you can go out a bit further, perhaps to n=8 or so without getting overly bored.

Historically, in smallish samples, it would be taken further by various (slightly) fancy forms of "counting". For example, you could derive a recursive relationship between the probability function at some $$n$$ and the probability function at a smaller sample size (say $$n-1$$), and then use that to go down to a case you can easily do by literal counting (as above) and build up a table from that.

For example to go from n=3 to n=4, we can duplicate the above table and add a rank 4 column with "-" for the top copy and "+" for the lower copy and recompute, but we can save quite a bit of work if we arrange the calculations better.

Note that if $$W_n$$ is the sum of positive ranks at sample size $$n$$, then at $$n=4$$ we can compute $$P(W_4=4)$$ by noticing that we need either $$W_3=4$$ and that the value with rank 4 was negative or $$W_3=0$$ and rank 4 was positive - so we simply average those two probabilities from the next lower table. That "shifting-by-n and averaging" of values in the smaller table works for the whole thing. You take a copy of the vector of probabilities we have above, move it 4 (i.e. $$n$$) places further up and just average.

W   0     1     2     3     4     5     6     7     8     9    10
p- 1/8   1/8   1/8   2/8   1/8   1/8   1/8    0     0     0     0
p+  0     0     0     0    1/8   1/8   1/8   2/8   1/8   1/8   1/8
-----------------------------------------------------------------
1/16  1/16  1/16  2/16   2/16  2/16  2/16  2/16  1/16  1/16  1/16


This recursion can allow us to expand the table rapidly (and we can still use the idea of only worrying about the tails if we want, though it requires us to be a bit fancier with our counting).

Such an approach allows us to push the table up considerably higher, and would be sufficient for the sort of small table most people would need to work by hand (since you only need to do the work once you can afford to put in a little labor). I could easily imagine pushing the table up to n=20 or even considerably higher without use of a computer even with as basic an approach as that simple recursion.

This recursion idea is not the only approach, but it's simple to use and sufficient to get the right sense of the kind of thing you can do.

Now, with computers to do the calculational work, we can go much further with the exact distribution than is convenient by hand. The dsignrank, psignrank and qsignrank functions in R can do probabilities, tail probabilities and critical values for say n=1000 in seconds, but seems not to be able to manage n=1040. However, there's little need for us to have exact values anywhere near up this high in any case.