[R] permutation test on paired samples

[Holger Taschenberger]

Henric,
        thanks a lot for the clarification. I guess your response also
implies that there is no 'simple' equivalent or replacement for the
(exactRankTests) function perm.test(y,x,paired = TRUE,...) in the
package coin. Perhaps it will be added in the future. Meanwhile I can
use your 'recipe'.

--Holger

On Sun, 15 Jul 2012 19:36:01 +0200
"Henric (Nilsson) Winell" <nilsson.henric at gmail.com> wrote:

> Holger,
> 
> Thanks for providing a reproducible example.  However, since your
> space key only works sporadically, the below is a little hard to read... ;)
> 
> On 2012-07-12 20:26, Holger Taschenberger wrote:
>  > Hi,
>  >
>  >          I'm trying to run a permutation test on paired samples.
>  >
>  > First I tried the package "exactRankTests":
>  >
>  > require("exactRankTests")
>  > x <- c(1.83,0.50,1.62,2.48,1.68,1.88,1.55,3.06,1.30)
>  > y <- c(0.878,0.647,0.598,2.05,1.06,1.29,1.06,3.14,1.29)
> 
> The relevant output missing here is
> 
>  > wilcox.test(x,y,paired = TRUE,alternative = "greater")
> 
>          Wilcoxon signed rank test
> 
> data:  x and y
> V = 40, p-value = 0.01953
> alternative hypothesis: true location shift is greater than 0
> 
>  > perm.test(y,x,paired = TRUE,exact = TRUE,alternative = "greater")
> 
>          1-sample Permutation Test (scores mapped into 1:m using rounded
>          scores)
> 
> data:  y and x
> T = 41, p-value = 0.003906
> alternative hypothesis: true mu is greater than 0
> 
> 
> Firstly, you've interchanged the 'x' and 'y' in the second call. 
> Secondly, and more important, the output says that "(scores mapped into 
> 1:m using rounded scores)".  In this case this can easily be avoided, 
> and note the interchange of 'x' and 'y' to match your 'wilcox.test' 
> call, using:
> 
>  > yy <- 1000 * y
>  > xx <- 1000 * x
>  > perm.test(xx, yy, paired = TRUE, exact = TRUE,
> +           alternative = "greater")
> 
>          1-sample Permutation Test
> 
> data:  xx and yy
> T = 4114, p-value = 0.01367
> alternative hypothesis: true mu is greater than 0
> 
> 
> So, now that we've computed the correct p-value, let's see how to obtain 
> this using the 'coin' package:
> 
>  >
>  > Then I wanted to use the package 'coin':
>  >
>  > require("coin")
>  > x <- c(1.83,0.50,1.62,2.48,1.68,1.88,1.55,3.06,1.30)
>  > y <- c(0.878,0.647,0.598,2.05,1.06,1.29,1.06,3.14,1.29)
>  > xydat <- data.frame(y = c(y,x),x = gl(2,length(x)),block = 
> factor(rep(1:length(x),2)))
> 
> The relevant output missing here is
> 
>  > wilcoxsign_test(y ~ x | block,data = xydat,alternative = 
> "greater",distribution = exact())
> 
>          Exact Wilcoxon-Signed-Rank Test
> 
> data:  y by x (neg, pos)
>           stratified by block
> Z = 2.0732, p-value = 0.01953
> alternative hypothesis: true mu is greater than 0
> 
>  > oneway_test(y ~ x | block,data = xydat,alternative = 
> "greater",distribution = exact())
> 
>          Exact 2-Sample Permutation Test
> 
> data:  y by x (1, 2)
>           stratified by block
> Z = -2.1948, p-value = 0.6982
> alternative hypothesis: true mu is greater than 0
> 
> 
> Using 'oneway_test' in this way does *not* correspond to a paired test. 
>   The "raw" scores version of the Wilcoxon signed-rank test can be 
> constructed using
> 
>  > diff <- x - y
>  > y <- as.vector(t(cbind(abs(diff) * (diff < 0),
> +                        abs(diff) * (diff >= 0))))
>  > x <- factor(rep(c("neg", "pos"), length(diff)),
> +             levels = c("pos", "neg"))
>  > b <- gl(length(diff), 2)
>  >
>  > oneway_test(y ~ x | b, alternative = "greater", distr = "exact")
> 
>          Exact 2-Sample Permutation Test
> 
> data:  y by x (pos, neg)
>           stratified by b
> Z = 2.1948, p-value = 0.01367
> alternative hypothesis: true mu is greater than 0
> 
> 
> And, as you can see, this is equal to the 'perm.test' result.
> 
> 
> HTH,
> Henric
> 
> 
> 
>  >
>  > While the results of the Wilcoxon test are the same for both packages
>  > are the same, those of the permutation test are very different. So,
>  > obviously I'm doing something wrong here. Can somebody please help?
>  >
>  > Thanks a lot,
>  >          Holger
>  >
>  > ______________________________________________
>  > R-help at r-project.org mailing list
>  > https://stat.ethz.ch/mailman/listinfo/r-help
>  > PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
>  > and provide commented, minimal, self-contained, reproducible code.
>  >

-- 
Holger Taschenberger, PhD
Dept. of Membrane Biophysics
Max Planck Institute for Biophysical Chemistry
Am Fassberg 11
D-37077 Goettingen, Germany
Tel.: ++49-551-201-1668
Fax:  ++49-551-201-1688
E-mail:     <Holger.Taschenberger at mpi-bpc.mpg.de>