# EXAMPLE OF AN R PROGRAM TO DO A PERMUTATION TEST. # # R. M. Neal, January 2000 # COMPUTE THE P-VALUE FOR A PERMUTATION TEST OF CORRELATION. Tests the null # hypothesis that the the vectors x and y are independent, versus the # alternative that they are correlated (either positively or negatively). # The vectors x and y are given as the first and second arguments; they # must be of equal length. # # The p-value returned is computed by simulating permutations of how the # elements of the vectors are paired up, with the simulation sample # size being given as the third argument, n, which defaults to 999. The # p-values returned are integer multiples of 1/(n+1), and have the property # that if the null hypothesis is true, the probability of obtaining a p-value # of k/(n+1) or smaller is equal to k/(n+1), unless there is exact equality # for the correlations obtained with different permutations, in which case the # probability may differ slightly from this. perm.cor.test <- function (x, y, n=999) { real.abs.cor <- abs(cor(x,y)) number.as.big <- 0 for (i in 1:n) { if (abs(cor(x,sample(y))) >= real.abs.cor) { number.as.big <- number.as.big + 1 } } (number.as.big + 1) / (n + 1) } # TEST ON NORMALLY-DISTRIBUTED DATA, COMPARED TO TEST BASED ON NORMAL DIST. test.norm <- function() { set.seed(1) x <- rnorm(20) y <- rnorm(20) cat("Permutation test p-value:",perm.cor.test(x,y),"\n\n") print(summary(lm(y~x))) invisible() } # TEST ON DATA WITH ONE EXTREME VALUE, COMPARED TO TEST BASED ON NORMAL DIST. test.extreme <- function() { set.seed(1) x <- c (1.1, 2.1, 0.3, 9.1, 1.7) y <- c (0.1, 0.9, 1.1, 8.2, 0.3) cat("Permutation test p-value:",perm.cor.test(x,y),"\n\n") print(summary(lm(y~x))) invisible() }