9
$\begingroup$

Permutation tests are significance tests based on permutation resamples drawn at random from the original data. Permutation resamples are drawn without replacement, in contrast to bootstrap samples, which are drawn with replacement. Here is an example I did in R of a simple permutation test. (Your comments are welcome)

Permutation tests have great advantages. They do not require specific population shapes such as normality. They apply to a variety of statistics, not just to statistics that have a simple distribution under the null hypothesis. They can give very accurate p-values, regardless of the shape and size of the population (if enough permutations are used).

I have also read that it is often useful to give a confidence interval along with a test, which is created using bootstrap resampling rather than permutation resampling.

Could you explain (or just give the R code) how a confidence interval is constructed (i.e. for the difference between the means of the two samples in the above example) ?

EDIT

After some googling I found this interesting reading.

| cite | improve this question | |
$\endgroup$
7
$\begingroup$

It's OK to use permutation resampling. It really depends on a number of factors. If your permutations are a relatively low number then your estimation of your confidence interval is not so great with permutations. Your permutations are in somewhat of a gray area and probably are fine.

The only difference from your prior code is that you'd generate your samples randomly instead of with permutations. And, you'd generate more of them, let's say 1000 for example. Get the difference scores for your 1000 replications of your experiment. Take the cutoffs for the middle 950 (95%). That's your confidence interval. It falls directly from the bootstrap.

You've already done most of this in your example. dif.treat is 462 items long. Therefore, you need the lower 2.5% and upper 2.5% cut offs (about 11 items in on each end).

Using your code from before... 

y <- sort(dif.treat)
ci.lo <- y[11]
ci.hi <- y[462-11]

Off hand I'd say that 462 is a little low but you'll find a bootstrap to 10,000 comes out with scores that are little different (likely closer to the mean).

Thought I'd also add in some simple code requiring the boot library (based on your prior code).

diff <- function(x,i) mean(x[i[6:11]]) - mean(x[i[1:5]])
b <- boot(total, diff, R = 1000)
boot.ci(b)
| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you. Is it ok to generate the samples using sample and replace=TRUE? Is there any reason to use a package like boot? $\endgroup$ – George Dontas Aug 30 '10 at 20:36
  • $\begingroup$ Typically it's done with replacement so you want to set that to TRUE. As for why... the package is optimized so it should run faster.... never timed it. That can be a problem if you set R large. And, as you can see, the code is nice and concise. It also has lots of features you wouldn't get easily rolling your own. $\endgroup$ – John Aug 30 '10 at 23:10
  • $\begingroup$ boot.ci returns the confidence interval. Is there any (boot) function that gives the p.value? (as the ratio of the number of differences at least as high as the one observed, over the total number of generated samples) $\endgroup$ – George Dontas Aug 31 '10 at 8:49
  • $\begingroup$ ok, I found a way to get it: sum(b$t>=b$t0)/b$R $\endgroup$ – George Dontas Aug 31 '10 at 9:02
  • $\begingroup$ @gd047 : take into account that this is a one-sided p-value you're calculating. $\endgroup$ – Joris Meys Aug 31 '10 at 15:58
4
$\begingroup$

As a permutation test is an exact test, giving you an exact p-value. Bootstrapping a permutation test doesn't make sense.

Next to that, determining a confidence interval around a test statistic doesn't make sense either, as it is calculated based on your sample and not an estimate. You determine confidence intervals around estimates like means and the likes, but not around test statistics.

Permutation tests should not be used on datasets that are so big you can't calculate all possible permutations any more. If that's the case, use a bootstrap procedure to determine the cut-off for the test statistic you use. But again, this has little to do with a 95% confidence interval.

An example : I use here the classic T-statistic, but use a simple approach to bootstrapping for calculation of the empirical distribution of my statistic. Based on that, I calculate an empirical p-value :

x <- c(11.4,25.3,29.9,16.5,21.1)
y <- c(23.7,26.6,28.5,14.2,17.9,24.3)

t.sample <- t.test(x,y)$statistic
t.dist <- apply(
      replicate(1000,sample(c(x,y),11,replace=F)),2,
      function(i){t.test(i[1:5],i[6:11])$statistic})

# two sided testing
center <- mean(t.dist)
t.sample <-abs(t.sample-center)
t.dist <- abs(t.dist - center)
p.value <- sum( t.sample < t.dist ) / length(t.dist)
p.value

Take into account that this 2-sided testing only works for symmetrical distributions. Non-symmetrical distributions are typically only tested one-sided.

EDIT :

OK, I misunderstood the question. If you want to calculate a confidence interval on the estimate of the difference, you can use the code mentioned here for bootstrapping within each sample. Mind you, this is a biased estimate: generally this gives a CI that is too small. Also see the example given there as a reason why you have to use a different approach for the confidence interval and the p-value.

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ Can you give an explanation of why permutation tests should not be used on datasets that you can't calculate all possible permutations? $\endgroup$ – Andy W Aug 31 '10 at 16:27
  • $\begingroup$ @Andy W : First define "permutation test". for me, permutation tests are exact tests, using every permutation possible. That is impossible on larger datasets. The "approximate permutation tests" are in fact the simples Monte Carlo method around, and should be adressed that way. Next to that, the central limit theorem assures in most cases that the assumptions regarding the distribution of the test statistics are met when using large datasets. In complex testing, the use of permutation tests on large datasets makes calculation times unbearably long without adding any significant value. my2cents $\endgroup$ – Joris Meys Aug 31 '10 at 16:42
  • $\begingroup$ I didn't say anything like bootstrapping a permutation test. I came into this question after reading the last paragraph of [SECTION 14.5 | Summary], in the linked pdf. $\endgroup$ – George Dontas Aug 31 '10 at 18:33
  • $\begingroup$ @gd047 Then I have misread your question. But you should really keep confidence intervals and p.values strictly separated. The confidence interval is estimated based on bootstrapping within each sample (although it is biased by definition), the permutation test is done by permutations over the complete dataset. Those are two completely different things. $\endgroup$ – Joris Meys Sep 1 '10 at 8:15
  • $\begingroup$ @Kevin : Code was darn right. Read the code again: the x[6:11] refers to the argument x of the anonymous function within the apply. Confusing maybe, but your edit gave very wrong results. Please comment about what you think it should be before editing the code. Saves me a rollback. To avoid further confusion, I changed that x to i $\endgroup$ – Joris Meys Jun 13 '12 at 12:06
0
$\begingroup$

From Joris Meys code in the Answers but with modification to allow it to be applied in more t han a single situation:

I tried to edit the other one but I did not have time to finish and for some reason I can't comment (maybe because this is an old question).

x <- c(11.4,25.3,29.9,16.5,21.1)
y <- c(23.7,26.6,28.5,14.2,17.9,24.3)

t.sample <- t.test(x,y)$statistic

t.dist <- apply(
          replicate(1000,sample(c(x,y),length(c(x,y)),replace=F)), 2,
          function(i){t.test(i[1:length(x)],i[length(x)+1:length(c(x,y))])$statistic})

# two sided testing
center <- mean(t.dist)
t.sample <-abs(t.sample-center)
t.dist <- abs(t.dist - center)
p.value <- sum( t.sample < t.dist ) / length(t.dist)
p.value
| cite | improve this answer | |
$\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.