尽管我在评论中指出了该coin
包的使用,但我认为值得说明的是排列/随机化测试非常简单,所以我已经完成了。
在这里,我编写了一些 R 代码来对位置的一个样本测试进行随机化测试。测试随机翻转差异上的符号并计算平均值;这相当于将每对值随机分配给 x 和 y 组。下面的代码可以大大缩短(我可以很容易地用两行来完成,如果你不介意更慢的代码,甚至可以用两行来完成)。
这段代码在我的机器上需要几秒钟:
# assumes the two samples are in 'x' and 'y' and x[i] and y[i] are paired
# set up:
B <- 99999
d <- x-y
m0 <- mean(d)
# perform a one-sample randomization test on d
# for the null hypothesis H0: mu_d = 0 vs H1 mu_d != 0 (i.e. two tailed)
# here the test statistic is the mean
rndmdist <- replicate(B,mean((rbinom(length(d),1,.5)*2-1)*d))
# two tailed p-value:
sum( abs(rndmdist) >= abs(m0))/length(rndmdist)
这就是全部。
请注意,rbinom(length(d),1,.5)*2-1)
给出一个随机-1
或1
...即一个随机符号,所以当我们乘以任何一组有符号d
时,它相当于随机分配+
或-
符号给绝对差。[不管你身上的标志分布是什么d
,现在d
都会有随机标志。]
在这里,我将它与一些虚构数据的 t 检验进行比较:
set.seed(seed=438978)
z=rnorm(50,10,2)
x=z-rnorm(50,0,.5)
y=z+.4+rnorm(50,0,.5)
t.test(y-x) # gives p = 0.003156
B <- 99999
d <- x-y
m0 <- mean(d)
rndmdist <- replicate(B,mean((rbinom(length(d),1,.5)*2-1)*d))
sum( abs(rndmdist) >= abs(m0))/length(rndmdist)
当 t 检验有效时,它通常会给出与完全枚举置换检验非常相似的 p 值,并且如上所述的模拟 p 值(当模拟数量足够大时)将收敛到第二个 p 值。
在上面使用的重复次数下,大约 85% 的时间将估计 0.05 的真实排列 p 值(即来自完全枚举)在 0.001 以内(即,将给出介于 0.049 和 0.051 之间的随机化 p 值)并且超过 99.5% 的时间在 0.002 以内。