9
$\begingroup$

I'm not an expert, so forgive me if some of the terminology is a little clumsy. Happy to provide more information where required.

I have two vectors of 50 paired numeric values in R. I want to perform a two-tailed randomisation or permutation test to determine whether their differences are due to chance or not.

A permutation test (also called a randomization test, re-randomization test, or an exact test) is a type of statistical significance test in which the distribution of the test statistic under the null hypothesis is obtained by calculating all possible values of the test statistic under rearrangements of the labels on the observed data points.

I want to do this type of test because I believe the distributions of the values in the vectors violate the assumptions of other tests such as the t-test (for example, many of the numeric values in the vector are 0).

The permtest function in the BHH2 library, almost does what I want, but it operates on all $2^{50}$ permutations, which will take too long. Instead, I want to estimate the p-value, by sampling a large number of the possible permutations. I had a look in the coin package, but nothing in there seems to do a permutation test with sampling from paired numeric vectors.

Some googling lead me to this email, which suggests that the reason I can't find a package to do it is that it's a one-liner in R. Unfortunately, I'm not experienced enough with R to be able to produce that one-liner.

Is there a package or method that will perform a two-tailed paired permutation test using only a sample of the permutation space?

If not, would somebody be able to share a short bit of R code to do it?

$\endgroup$
  • 3
    $\begingroup$ It looks to me like the package coin (among several others) does randomization tests. e.g. see the answer to this question (read the whole thing). If I understand right, the examples cover both approximate and exact cases and cover both independent and dependent samples. $\endgroup$ – Glen_b Jul 13 '13 at 6:48
  • 1
    $\begingroup$ Sorry, to be clear - by 'read the whole thing' I mean 'read the top answer all the way through' -- though you might also want to look at the bottom answer. $\endgroup$ – Glen_b Jul 13 '13 at 6:54
  • $\begingroup$ Pretty much the only interesting bit from that answer for paired permutations is oneway_test(y ~ x | pairs, distribution=approximate(B=9999)) with library(coin). $\endgroup$ – Nakx Mar 13 at 4:59
12
$\begingroup$

Though I pointed in comments to the use of the coin package I think it's worth illustrating that a permutation/randomization test is really quite simple, so I have done it.

Here I write some R code to do a randomization test for a one sample test of location. The test randomly flips signs on the differences and computes the mean; this is equivalent to randomly assigning each pair of values to the x and y groups. The code below could be made significantly shorter (I could do it in two lines easily enough, or even one if you didn't mind slower code).

This code takes a few seconds on my machine:

# assumes the two samples are in 'x' and 'y' and x[i] and y[i] are paired
# set up:
B <- 99999
d <- x-y
m0 <- mean(d)

# perform a one-sample randomization test on d
# for the null hypothesis H0: mu_d = 0   vs H1 mu_d != 0  (i.e. two tailed)
# here the test statistic is the mean
rndmdist <- replicate(B,mean((rbinom(length(d),1,.5)*2-1)*d))

# two tailed p-value:
sum( abs(rndmdist) >= abs(m0))/length(rndmdist)

That's the whole thing.

Note that rbinom(length(d),1,.5)*2-1) gives a random -1 or 1 ... i.e. a random sign, so when we multiply by any set of signed d, it is equivalent to randomly assigning + or - signs to the absolute differences. [It doesn't matter what distribution of signs on d you start with, now the d will have random signs.]

Here, I compare it with a t-test on some made up data:

 set.seed(seed=438978)
 z=rnorm(50,10,2)
 x=z-rnorm(50,0,.5)
 y=z+.4+rnorm(50,0,.5)
 t.test(y-x) # gives p = 0.003156

 B <- 99999
 d <- x-y
 m0 <- mean(d)
 rndmdist <- replicate(B,mean((rbinom(length(d),1,.5)*2-1)*d))
 sum( abs(rndmdist) >= abs(m0))/length(rndmdist) 

When the t-test is valid it usually gives a very similar p-value to the completely enumerated permutation test, and a simulated p-value as above (when the number of simulations is sufficiently large) will converge to that second p-value.

At the number of replications used above, a true permutation p-value (i.e. from complete enumeration) of 0.05 will be estimated to within 0.001 (that is, will give a randomization p-value between 0.049 and 0.051) about 85% of the time and to within 0.002 over 99.5% of the time.

$\endgroup$
  • $\begingroup$ Very much appreciated, thank you. How did you calculate the accuracy of the p-value? $\endgroup$ – Timothy Jones Jul 17 '13 at 1:37
  • 1
    $\begingroup$ It's just the normal approximation to a binomial proportion, using the standard error of a proportion; $\text{s.e.}(\hat p) = \sqrt{p(1-p)/n}$. $\endgroup$ – Glen_b Jul 17 '13 at 2:21
  • $\begingroup$ Why do you multiply the rbinom function by 2-1? And then d? $\endgroup$ – user60285 Nov 8 '14 at 20:45
  • $\begingroup$ To get random signs on d, because that's how a permutation test of mean difference for paired data works. See new additional comments after that piece of code. $\endgroup$ – Glen_b Nov 8 '14 at 22:34
  • 1
    $\begingroup$ @Joe when we add the observed sample it will make a round number $\endgroup$ – Glen_b Sep 11 '17 at 1:53
0
$\begingroup$

Here is code to perform a permutation test. I have data in there for example. x is the difference between the two vectors.

x <- c(5.1, 9.4, 7.2, 8.1, 8.8, 2.5, 4.2, 6.9, 5.5, 5.3)
m = 5
n = 5
xsum = sum(x)
asum = sum(x[1:m])
bsum = xsum - asum
truediff = asum/m - bsum/n
truediff
abstruediff = abs(truediff)
iter = 100000
difflist <- 1:iter
for(i in 1:iter) {
  s <- sample(x,m) # select a sample of size m
  pasum = sum(s)
  pbsum = sum(x) - sum(s)
  diff  = pasum/m - pbsum/n
  difflist[i] <- diff # add permutation difference to list
}
difflist  <- sort(difflist)
xquantile <- quantile(difflist,probs=c(.005, .01, .025, .05, .95, .975, .99, .995))
xquantile
pdist  <- quantile(difflist, probs=seq(0,1,1/iter))
ntail1 <- length(pdist[difflist <= -abstruediff])
tail1  <- ntail1/iter
tail1  # left-tail probability
ntail2 <- length(pdist[difflist >= abstruediff])
tail2  <- ntail2/iter
tail2  # right-tail probability
twotail = tail1 + tail2
twotail 
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.