wilcoxon signed rank test to test paired data?

Question

I want to test whether the difference between the pairs of students's scores follows a symmetric distribution around zero.See below

Students    Before  After   difference
   1          **     **        10
   2          **     **        -9
   3          **     **        8
   4          **     **        5
   5          **     **       -6
   6          **     **        0
   7          **     **        2

Since the data is not normally distributed, I assume i can use wilcoxon signed-rank test here, but in this case we only know the difference of the pairs, so if perform wilcoxon signed-rank test in r here, normally the code should be:

wilcox.test(before,after,paired=T, exact=T)

but the problem is that, I don't have the data of 'before' and 'after', instead I only know the difference between them, so it seems wilcoxon signed-rank test can not be applied here. Is there other alternative method to this case?

Answer 1

This largely follows on from Antoni's answer.

Here's what I was getting at with a permutation test First let's exclude the tie as Antoni did:

  d   s  r 
 10   1  6
 -9  -1  5
  8   1  4
  5   1  2
 -6  -1  3
  0  NA NA
  2   1  1

Our test statistic is then 1+2+4+6 = 13

So we could calculate the permutation distribution from scratch as follows:

 pm <- c(-1,1)
 sign <- expand.grid(s1=pm,s2=pm,s3=pm,s4=pm,s5=pm,s6=pm)
 table(rowSums((sign>0)%*%cbind(1:6)))/64

which gives the same values as dsignrank(0:21,6) (yay!). So clearly when there are no ties (which is the case, since we excluded them) this is just the signed rank test, as it should be

Let's consider whether we can do anything with the tie. How we really deal with it depends on how it came about.

If, for example we envision an underlying continuous scale that has become discretized into categories, then our tie is simply due to two different values that our subsequently imposed categories are not fine enough to distinguish.

Arguably then, the two values are "really" different but we have no basis to conclude which sign we should have and so allowing for both possibilities, half the time giving it a +1 and half the time giving it a -1 and then averaging the two (or perhaps to be conservative, some would argue we should take the higher of the two possible counts ... but I'll continue to look at the average of the two possible cases here). [By contrast if the values were inherently discrete with no underlying-but-unobserved continuous scale, then we would say that the 0 was inherent and deal with it in a different way -- I think exclusion probably makes the most sense in that case]

Note tie will always have the smallest rank.

So we have a table like the one above but the second last row is either:

   d  s r 
   :  : :
   0 -1 1

or

   d  s r 
   :  : :
   0  1 1

and all the other rows would have r one higher than before.

  d   s  r 
 10   1  7
 -9  -1  6
  8   1  5
  5   1  3
 -6  -1  4
  0  ±1  1
  2   1  2

So the possible test statistics are 2+3+5+7 and 1+2+3+5+7 (17 and 18)

[One alternative approach, if we're being conservative would be to treat it as a full n=7 and take the larger of the two possible p-values which gives 0.6875]

If we do all the permutations but always average the two signs for the smallest rank, they'll contribute 1/2, and each of the other positive ranks will be 1 higher.

In effect we'd be doing this:

 table(0.5+rowSums((sign>0)%*%cbind(2:7)))/64

And comparing a test statistic of 17.5 with that. This gives p = 0.625.

So anyway, all the approaches are a bit different but they pretty much tell the same general story here.

Answer 2

I have been playing with this question, and here is what I got so far...

1. Looking under the hood of R wilcox.test():

I tried following on the side what R was doing in the back. Here is some code:

d = data.frame(difference = c(10,-9,8,5,-6,0,2))   # Differences in OP
d$sgn = sign(d$difference)                         # The sign function of these diff's
d$abs = abs(d$difference)                          # The absolute value of these diff's
d$abs = replace(d$abs, d$abs==0, NA)               # Getting the zeros out of the way
d$rank = rank(d$abs, na='keep')                    # Ranking the abs values
n = length(d$rank) - length(d$rank[is.na(d$rank)]) # Rows without zero differences (6)
d$multi = d$sgn * d$rank                           # Multiplying sign times rank.

yielding the following data.frame:

  difference sgn abs rank multi
1         10   1  10    6     6
2         -9  -1   9    5    -5
3          8   1   8    4     4
4          5   1   5    2     2
5         -6  -1   6    3    -3
6          0   0  NA   NA    NA
7          2   1   2    1     1

Now we can get the test statistic (V) by summing all the positive values in multi: sum(d[d$multi > 0, ]$multi, na.rm=T), and the output is $13$.

Now let's pretend its 1985, and we go to the back of the book to the tables to retrieve the $p$-value... OK, let's do it with R: psignrank(q = 13, n = n, lower.tail = T). The result: $0.71875$.

Is this concordant with the results of the wilcox.test function... No... and yes. Initially, we may be disappointed:

 wilcox.test(difference, correct = F, alternative = 'less', exact = T)

    Wilcoxon signed rank test

data:  difference
V = 13, p-value = 0.6999
alternative hypothesis: true location is less than 0

Warning message:
In wilcox.test.default(difference, correct = F, alternative = "less",  :
  cannot compute exact p-value with zeroes

but it's not like we are not warned... Now, let's just get rid of the annoying row with the tie, and recalculate:

wilcox.test(difference[-6], correct = F, alternative = 'less', exact = T)

    Wilcoxon signed rank test

data:  difference[-6]
V = 13, p-value = 0.7187
alternative hypothesis: true location is less than 0

Bingo! Here is the coveted $0.7187$, again.

---

2. But can we make it fancier and run a Monte Carlo?

set.seed(11)                              # Today's date in the US - no cherry-picking!
r = 1:6                                   # The possible ranks of our non-zero differences
nsim = 1e5 
V = 0 
for (i in 1:nsim){ 
  rank = sample(r)                        # Sampling the ranks...
  sign = sample(c(1, -1), 6, replace = T) #... and the signs for each rank.
  V[i] = sum(rank[sign > 0])              # Doing the sum to get the V.
} 
(p_value <- sum(V <= 13) / nsim)          # Fraction with a sum equal or lower than 13.

[1] 0.72204

Really cool! $0.72204$... so close...