# Resampling program for t test on two dependent samples

### The easiest way to run this is to enter
# source("DependentSamplet.R") after setting the default directory.
# The first screen will be asking for the data file.
#

# The data are set up with the ID is col 1, the first trial observations in col 2,
# and the second trial observations in col 3.
# For Example:
#ID Before After Gain
#01 83.8  95.2 11.4
#02 83.3  94.3 11.0
#03 86.0  91.5  5.5
#04 82.5  91.9  9.4
# The Gain scores and the ID are not used.


# There are column headings, but these are then ignored after they are read.
# The p value that is given is the two-tailed p.
# For and example with 35 cases and
# nreps = 10,000, this took 21 seconds.

data <- read.table(file.choose(), header = TRUE)
#Sys.time()  Uncomment this line (and the last line) if you want to time the run
x <- data[,3] - data[,2]    # x is the Gain score
nrep <- 10000
t <- numeric(nrep)
meangain <- numeric(nrep)
nrep <- 10000
N <- length(x)
 
 # Now set up function to compute t values
 ttestdep <- function(x)  {
   mean <- mean(x)
   var <- var(x)
   n <- length(x)
   se <- sqrt(var/n)
   t <- mean/se
   return(t)
   }
 tobt <- ttestdep(x)
 
 for (i in 1:nrep) {
    sign <- sample(c(1,-1),N,replace = TRUE)
    gain <- sign*x
    meangain[i] <- mean(gain)
    t[i] <- ttestdep(gain)
    }
    

tneg <- -1*abs(tobt)
tpos <- abs(tobt)
lower <- length(t[t <= tneg])
upper <- length(t[t >= tpos])
lowerp <- lower/nrep
upperp <- upper/nrep
twotailed <- lowerp + upperp
cat("Lower one-tailed probability is ", lowerp ,"\n")
cat("Upper one-tailed probability is ", upperp, "\n" )
cat("The two-tailed probability is ", twotailed," \n")


# To send to a pdf file instead, enter
#pdf("outfile.pdf")
# Then enter "dev.off() after it prints
# Set up histogram
dist <- seq(-6, 6,.1)
hist(t, breaks = 50, freq = FALSE, xlim = c(-5,5),main = "Distribution of t for randomization",
  col = "gray")
  # Superimpose the t distribution just for fun
quant <- seq(min(t), max(t), length = 100)
y <- dt(quant,  df = 41, ncp = 0)
lines(quant,y)
  
legend(2, .3, paste("t = ", round(tobt, digits = 4)), col=1,  cex = 0.8)
legend(2, .25, paste("p = ", round(twotailed, digits = 4)), col=1, cex = 0.8)

#dev.off()
 #Sys.time()