Create a gist now

Instantly share code, notes, and snippets.

A randomization method for a t-test
###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))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment