# Test of mean of paired samples.
# Data from Everitt -- Methods, p. 199
    #Before and after scores on weight gain in anorexic girls using family therapy
# Some of the ideas for this came from Maindonald and Braun (2007)


before <- c(83.8, 83.3, 86, 82.5, 86.7, 79.6, 76.9, 94.2, 73.4, 80.5,
81.6, 82.1, 77.6, 83.5, 89.9, 86.0, 87.3)
after <- c(95.2, 94.3, 91.5, 91.9, 100.3, 76.7, 76.8, 101.6, 94.9, 75.2,
77.8, 95.5, 90.7, 92.5, 93.8, 91.7, 98.0)
gain <- after - before
xbar <- mean(gain)
n <- length(gain)
nsim <- 10000
results <- numeric(nsim)

cat("The mean of the before scores = ", mean(before),'\n','\n')
cat("The mean of the after scores = ",mean(after),'\n','\n')
cat("The mean gain score = ",mean(gain),'\n','\n')

par(mfrow = c(2,1))
plot(before, after)
abline(lm(after~before))
cat("The correlation between before and after = ",cor(before, after),'\n','\n')

# If there is no treatment effect, each gain score is as likely to be positive as negative


for (i in 1:nsim) {
    sign <- sample(c(1,-1),n, replace = T)
    newSample <- sign*abs(gain)
    results[i] <- mean(newSample)
    }
pval <- (sum(results >= abs(xbar)) + sum(results <= -abs(xbar)))/nsim

 plot(density(results), xlab = "Mean Difference", main = "", yaxs = "i", cex.axis = 0.8,
            bty = "L")
        abline(v = xbar)
        abline(v = -xbar, lty = 2)
        mtext(side = 3, line = 0.5, text = expression(bar(Diff)),
            at = xbar)
        mtext(side = 3, line = 0.5, text = expression(-bar(Diff)),
            at = -xbar)
        text(0, .04, "p value = " )
        text(0, .02,  pval )
cat("The probability given no effect is = ",pval)