# EXAMPLE OF AN R PROGRAM TO DO A PERMUTATION TEST.
#
# R. M. Neal, 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()
}