set.seed(12345); #set the random number generator
x = rnorm(100);
y = rnorm(100,mean=0.5)
- t-test, although they are clearly not normal
t.test( x, y, alternative="two.sided");
t.test( x, y, alternative="less"); #one-sided test
tt = t.test( x, y, alternative="less");
-
- alternatively, let's try simulation
t.stat = function( x1, x2 ) {
( mean(x1) - mean(x2)) / sqrt( var(x1)/length(x1) + var(x2)/length(x2));
}
t.obs = t.stat(x,y);
n=1000 #1K permutations
permu = data.frame( matrix(NA, nrow=n, ncol=3) );
names(permu) = c("xbar","ybar","t");
head(permu);
- Reset:
- RNGkind(ok[1])
i=1;
while( i <= n ) {
new = sample( c(x,y), replace = TRUE);
newX = new[1:100];
newY = new[101:200];
permu$t[i] = t.stat( newX, newY);
permu$xbar = mean(newX);
permu$ybar = mean(newY);
i = i+1;
}
sub = permu[permu$t< t.obs, ]
my.p = length(sub$t) / length(permu$t);
my.p
results = rbind(my.p, tt$p.value);
rownames(results) = c("permutation", "t-test");
colnames(results) = c("p");
results;
hist( permu$t, breaks=20, xlab="t-statistic", main="p-value and simulation" );
arrows( t.obs, n/10, t.obs, n/20, col="red");