A randomization method for a t-test

randttest.R

###Example 1
# Simple test of difference of two means
####

#############Note I do some plots in regular and again in ggplot, so install ggplot if you want to try it.
#install.packages("ggplot2")
library(ggplot2)
###Generate two random data sets###

#set sample size
n <- 15

####Create data frame

s.dat <- data.frame(cbind(c(rnorm(n,20,3),rnorm(n,27,3)),c(rep(1,n),rep(2,n))))
colnames(s.dat) <- c("Measure","Group")
####Stanard R plot
plot(density(s.dat[s.dat[,2]==1,1]),xlim=c(10,40))
lines(density(s.dat[s.dat[,2]==2,1]),col=2)
####ggplot2 code
ggplot(s.dat,aes(x=Measure,group=as.factor(Group),fill=as.factor(Group),alpha=.5))+geom_density()

####Parametric t.test#####
tt1 <- t.test(x=s.dat[s.dat[,2]==2,1],y=s.dat[s.dat[,2]==1,1])

#####Randomization test#######
# In this case the test metric
# is if the difference between the two means
# is different from 0
############

##Set the number of randomizations
nran <- 5000
output <- vector()
for(i in 1:nran){
  ###Shuffle association via sample()
  new.ind <- sample(s.dat$G,length(s.dat$G),replace=F)
  output[i] <- mean(s.dat$M[new.ind==1]) - mean(s.dat$M[new.ind==2])
}

##Find true metric
true.val <- mean(s.dat[s.dat[,2]==1,1])-mean(s.dat[s.dat[,2]==2,1])
###Normal R plot
plot(density(output))
abline(v=true.val,col=2)
####ggplot2###
ggplot(as.data.frame(output),aes(x=output))+geom_density()+xlim(-10,10)+geom_vline(xintercept=true.val)

###Check true value against quantiles
quantile(output,c(.025,.975))