#  Randomization Test on Matched Samples
#  Data from Everitt (1994)

setwd("C:/Users/Dave/Dropbox/Webs/StatPages/ResamplingWithR/RandomMatchedSample")
dat <- read.table("Everitt.dat", header = T)
diffObt <- mean(dat$After) - mean(dat$Before)
difference <- dat$After - dat$Before  #Use that order to keep most diff. positive

nreps <- 10000
set.seed(1086)
resampMeanDiff <- numeric(nreps)
for (i in 1:nreps) {
       signs <- sample( c(1,-1),length(difference), replace = T)
       resamp <- difference * signs
       resampMeanDiff[i] <- mean(resamp)
       }
       diffObt <- abs(diffObt)
highprob <- length(resampMeanDiff[resampMeanDiff >= diffObt])/nreps
lowprob <- length(dat$resampMeanDiff[dat$resampMeanDiff <= (-1)*dat$diffObt])/nreps
prob2tailed <- lowprob + highprob
cat("The probability from the sampling statistics is = ",prob2tailed,'\n')


hist(resampMeanDiff, breaks = 30, main = "Distribution of Mean Differences",
  xlab = "Mean Difference", freq = FALSE)
  text(1.5,.25,"Diff. obt")
  text(1.5,.23,round(diffObt,2))
  arrows(1.5, .21, diffObt, 0, length = .125)
  text(-3,.25,"p-value")
  text(-3,.23, prob2tailed)


  # Compare to Student's t
tvalue <- t.test(dat$Before, dat$After, paired = T)$statistic

cat("\n The t value from a standard matched-pairs t test is= ",tvalue,'\n')
cat("with a probability of .03502 \n")