#  Randomization test on two medians
#  Data from study by Adams et al. (Methods8, tab7.5)
setwd ("C:/Users/Dave/Dropbox/Webs/StatPages/ResamplingWithR")
dat <- read.table(file = "Tab7-5.dat", header = TRUE)
Group <- factor(Group)
attach(dat)
n1 <- length(Arousal[Group == 1])
n2 <- length(Arousal[Group == 2])
N <- n1 + n2
Group1 <-  Arousal[Group == 1]
Group2 <-  Arousal[Group == 2]
md1 <- median(Group1, na.rm = T)
md2 <- median(Group2, na.rm = T )
par(mfrow = c(2,1))
hist(Group1, col = "red", xlim = c(min(Arousal), max(Arousal)) )
hist(Group2, col = "green", xlim = c(min(Arousal), max(Arousal)))
obt.med.diff <- md1 - md2
cat("The obtained median difference is ", '\n', obt.med.diff, '\n')

nreps = 5000
samp.med.diff <- numeric(nreps)
for (i in 1:nreps) {
   resamp.data <- sample(Arousal, length(Arousal), replace = FALSE)
   Grp1 <-  resamp.data[1:n1]
   Grp2 <-  resamp.data[n1+1:N]
   md1 <- median(Grp1, na.rm = T)
   md2 <- median(Grp2, na.rm = T )
   samp.med.diff[i] <- md1 - md2
   }
obt.med.diff <- abs(obt.med.diff)
pUpper <- length(samp.med.diff[samp.med.diff >= obt.med.diff]) /nreps
pLower <- length(samp.med.diff[samp.med.diff <= (-1)*obt.med.diff])/nreps
pExtreme <- pLower + pUpper
cat("The two tailed probablility under the null is = ",'\n',pExtreme, '\n')

# Confidence limits
ordered.samp.med.diff <- sort(samp.med.diff)
CIupper <- ordered.samp.med.diff[.975*nreps]
CIlower <- ordered.samp.med.diff[.025*nreps]
cat("The 95% confidence limits are = ", '\n',CIlower, " and ", CIupper, '\n')

#Note the dramatic difference with the standard t test.
#SPSS gives similar results for Wilcoxon and t tests, so it 
#is probably not an error.
# Wilcoxon gives p = .315, t gives p = .015
#(t.test(Arousal~Group))     Remove # to get output