# Test of two means using data from ReamplingR,html
# Data from Aronson (1998) example, page 216
# I start off with a standard two-group t test.

control <- c(4, 9, 12, 8, 9, 13, 12, 13, 13, 7, 6)
threat <- c(7, 8, 7, 2, 6, 9, 7, 10, 5, 0, 10, 8)
m.control <- mean(control) ; var.control <- var(control)
m.threat <- mean(threat) ; var.threat <- var(threat)
n.1 <- length(control); n.2 <- length(threat)
#   var.pool <- ((n.1 - 1)*var.control + (n.2 - 1)*var.threat)/
   (n.1 + n.2 - 2)
#   t <- (m.control - m.threat)/sqrt(var.pool/n.1 + var.pool/n.2)
t.test(control, threat, var.equal = T)    # This will  pool the variances


meanDiff <- m.control - m.threat
cat("The mean difference is = ",meanDiff, '\n')
All <- c(control, threat)
combinations <- factorial(23)/(factorial(11)*factorial(12))
cat("The number of possible different rearrangements is = ", combinations, '\n' )


# There are 1,352,078  possible arrangements of the data with 11 scores in control
# and 12 scores in threat. Way to many to analyze.

sumAll <- sum(All)
nreps <- 10000
results <- numeric(nreps)
for (i in 1:nreps) {
   mean1 <- mean(sample(All, n.1))
   mean2 <- (sumAll-11*mean1)/n.2
   results[i] <- mean1 - mean2    #mean difference[i]
   }

par(mfrow = c(1,1))
cumResults <- sort(results)
ordinate <- seq(.01, 100,.01)
plot(cumResults,ordinate, xlab = "Cumulative Results", type = "s",
ylab = "percentage", main = "Cumulative Distribution" )
arrows(-4,20, -3.06, 0, .125, 30)
text(-4,20,"-3.06")
arrows(4,20, 3.06, 0, .125, 30)
text(4,20,"3.06")
polygon(c(cumResults[cumResults<=-3.04],-3.04),c(ordinate[cumResults<=-3.04], 0),col="red")

   #Use absolute values for two-tailed test.
absResults <- abs(results)
meanDiff <- abs(meanDiff)

hist(results, breaks = 25, main = "Distribution of Mean Differences", xlab = "Mean Difference")
arrows(4,1000, 3.06, 0, .125, 30)
text(4,1000,"3.06")
arrows(-4,1000, -3.06, 0, .125, 30)
text(-4,1000,"-3.06")

two.tailed <- length(absResults[absResults >= meanDiff])/nreps
cat("The probability of a difference greater than obtained meanDiff = ",two.tailed,'\n')