R中成对向量的随机化/排列测试

机器算法验证 r 统计学意义 置换检验
2022-03-22 05:06:49

我不是专家,所以如果某些术语有点笨拙,请原谅我。很高兴在需要时提供更多信息。

我在 R 中有两个由 50 个配对数值组成的向量。我想执行双尾随机化或排列测试,以确定它们的差异是否是偶然的。

置换检验(也称为随机化检验、重新随机化检验或精确检验)是一种统计显着性检验,其中通过计算检验统计量的所有可能值获得检验统计量在原假设下的分布在观察数据点上的标签重新排列下。

我想做这种类型的测试,因为我认为向量中值的分布违反了其他测试的假设,例如 t 检验(例如,向量中的许多数值为 0)。

BHH2 library 中的函数permtest几乎可以满足我的要求,但它对所有排列进行操作,这将花费太长时间。相反,我想通过对大量可能的排列进行采样来估计 p 值。我查看了coin包,但其中似乎没有使用配对数字向量采样进行置换测试。250

一些谷歌搜索将我带到这封电子邮件,这表明我找不到包的原因是它是 R 中的单线。不幸的是,我对 R 的经验不足,无法制作那个-衬垫。

是否有仅使用置换空间样本执行双尾配对置换测试的包或方法?

如果没有,有人可以分享一小段 R 代码吗?

2个回答

尽管我在评论中指出了该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)给出一个随机-11...即一个随机符号,所以当我们乘以任何一组有符号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 以内。

这是执行置换测试的代码。例如,我在那里有数据。x 是两个向量之间的差。

x <- c(5.1, 9.4, 7.2, 8.1, 8.8, 2.5, 4.2, 6.9, 5.5, 5.3)
m = 5
n = 5
xsum = sum(x)
asum = sum(x[1:m])
bsum = xsum - asum
truediff = asum/m - bsum/n
truediff
abstruediff = abs(truediff)
iter = 100000
difflist <- 1:iter
for(i in 1:iter) {
  s <- sample(x,m) # select a sample of size m
  pasum = sum(s)
  pbsum = sum(x) - sum(s)
  diff  = pasum/m - pbsum/n
  difflist[i] <- diff # add permutation difference to list
}
difflist  <- sort(difflist)
xquantile <- quantile(difflist,probs=c(.005, .01, .025, .05, .95, .975, .99, .995))
xquantile
pdist  <- quantile(difflist, probs=seq(0,1,1/iter))
ntail1 <- length(pdist[difflist <= -abstruediff])
tail1  <- ntail1/iter
tail1  # left-tail probability
ntail2 <- length(pdist[difflist >= abstruediff])
tail2  <- ntail2/iter
tail2  # right-tail probability
twotail = tail1 + tail2
twotail