#this code performs a statistical test in which a given (user-selected) statistic is compared against a reference distribution which is generated by means of permuting the original data and recalculating the test-statistic a large number of times. Data should be in a *.csv file as one column of number. This code only works for one group (one-sample t-test) or two groups. The user needs to select the appropriate test-statistic, the options that are implemented are listed below. 

#Roel Willems, June 2010


#remove all varaibles:
rm(list=ls(all=TRUE))
#set working dir:
setwd("/Users/roelwillems/")

#read in data:
data<- scan ('untitled.csv',sep=',')

#length of the 2 data sets (if only one data set, N1=N2=N/2; if N is uneven in this case you need to rewrite):
N1=77
N2=77

indices=c(rep(1,N1), rep(2,N2))

#which statistic do you want to use:
statistic='median_diff'
statistic='2sample_t'
#statistic='paired_t'
#statistic='mean_diff'
#statistic='ks'
#statistic='1sample_t'

#number of perms:
perms=2500

#pre-allocate newstat variable:
newstat=(1:length(perms))
obstat=0

#observed stat:
	if (statistic=='2sample_t') {obstat=t.test(data[indices==1],data[indices==2])$stat} else
	if (statistic=='1sample_t') {obstat=t.test(c(data[indices==1],data[indices==2]))$stat} else
	if (statistic=='paired_t') {obstat=t.test(data[indices==1],data[indices==2],paired=TRUE)$stat} else
	if (statistic=='median_diff') {obstat=median(data[indices==1])-median(data[indices==2])} else
	if (statistic=='mean_diff') {obstat=mean(data[indices==1])-mean(data[indices==2])} else
	if (statistic=='ks') {obstat=ks.test(data[indices==1],data[indices==2])$stat} else
	print ('error, statistic unknown')

for (i in 1:perms) {
	data_perm=sample(data)
	if (statistic=='2sample_t') {newstat[i]=t.test(data_perm[indices==1],data_perm[indices==2])$stat} else
	if (statistic=='1sample_t') {newstat[i]=t.test(c(data[indices==1],data[indices==2]))$stat} else
	if (statistic=='paired_t') {newstat[i]=t.test(data_perm[indices==1],data_perm[indices==2],paired=TRUE)$stat} else
	if (statistic=='median_diff') {newstat[i]=median(data_perm[indices==1])-median(data_perm[indices==2])} else
	if (statistic=='mean_diff') {newstat[i]=mean(data_perm[indices==1])-mean(data_perm[indices==2])} else
	if (statistic=='ks') {newstat[i]=ks.test(data_perm[indices==1],data_perm[indices==2])$stat} else
	print ('error, statistic unknown')
	}
	
	
#plot	
hist(newstat)
abline(v=c(-1,1)*obstat,lty=2,col=2)

#find p-value;
pval=sum(abs(newstat)>=abs(obstat))/perms
print(pval)