# Randomization tests
# Dennis E. Slice (c) 2012

# Shake the tree
rm(list=ls())

# Set up
# Save the current state
original_ask = par()$ask
par(ask=T)

# Parameters
r = 0.3
n = 75

# Generate correlated, non-normal variables
y1 = runif(n)
y.temp = runif(n)
y2 = r*y1 + y.temp*sqrt(1-r*r)

plot(y1,y2)

test = cor.test(y1,y2)
print(test)
print("But,...")
print(shapiro.test(y1))
print(shapiro.test(y2))

hist(y1)
hist(y2)

# Generate distribution of uncorrelated r values
# Reasonable numbers vary. 
# 999 (min_p=0.001) to 9999 (min_p=0.0001) are reasonable.
nSamples = 9999
random_cor = rep(0,nSamples)
for (i in 1:nSamples) {
	# For long loops, print out some indicator of progress...
	if (!(i%%500)) {print(i)}
	# Randomize (=shuffle) one of the variables
	# No need to do both.
	rand.y2 = sample(y2,n)
	# Compute correlation
	r = cor(y1,rand.y2)
	# Save
	random_cor[i] = r	
}

# How many of the random correlations are as or more extreme
# than that observed in our sample.
count = 0;
r = test$estimate
observed_r = abs(r) # magnitude, not sign, matters
for (i in 1:length(random_cor)) {
	rand.r = abs(random_cor[i])
	# print(rand.r) # debugging, forgot 1: in loop
	# is random equal to or greater than?
	if (!(rand.r<observed_r)) {
		#increment the counter
		count = count + 1
	}
}
# p = proportion or equal or more extreme values
print(paste("Randomization p = ",count/(nSamples+1)))
hist(c(random_cor, r),freq=FALSE)
lines(c(r,r),c(0,10),col="red",lwd=2)

# Put things back where we found them
par(ask=original_ask)