public
Created

A randomization method for a t-test

  • Download Gist
randttest.R
R
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51
 
###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))

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.