使用 R 中的哪个排列测试实现来代替 t 测试(配对和非配对)?

机器算法验证 r t检验 非参数 置换检验
2022-02-10 04:24:57

我有来自我使用 t 检验分析的实验的数据。因变量是区间标度的,数据要么是不成对的(即,2 组),要么是成对的(即,受试者内)。例如(在科目内):

x1 <- c(99, 99.5, 65, 100, 99, 99.5, 99, 99.5, 99.5, 57, 100, 99.5, 
        99.5, 99, 99, 99.5, 89.5, 99.5, 100, 99.5)
y1 <- c(99, 99.5, 99.5, 0, 50, 100, 99.5, 99.5, 0, 99.5, 99.5, 90, 
        80, 0, 99, 0, 74.5, 0, 100, 49.5)

然而,数据不正常,因此一位审稿人要求我们使用 t 检验以外的方法。但是,很容易看出,数据不仅不是正态分布的,而且条件之间的分布也不相等: 替代文字

因此,不能使用通常的非参数检验,Mann-Whitney-U-检验(未配对)和 Wilcoxon 检验(配对),因为它们需要条件之间的相等分布。因此,我决定最好进行一些重采样或置换测试。

现在,我正在寻找基于置换的等效 t 检验的 R 实现,或者关于如何处理数据的任何其他建议。

我知道有一些 R 包可以为我做到这一点(例如 coin、perm、exactRankTest 等),但我不知道该选择哪一个。所以,如果有一些使用这些测试经验的人可以给我一个启动,那就太酷了。

更新:如果您能提供一个如何报告此测试结果的示例,那将是理想的。

4个回答

没关系,因为测试统计量始终是均值(或等价物)的差异。微小的差异可能来自蒙特卡洛方法的实施。通过对两个自变量的单向测试来尝试使用您的数据的三个包:

DV <- c(x1, y1)
IV <- factor(rep(c("A", "B"), c(length(x1), length(y1))))
library(coin)                    # for oneway_test(), pvalue()
pvalue(oneway_test(DV ~ IV, alternative="greater", 
                   distribution=approximate(B=9999)))
[1] 0.00330033

library(perm)                    # for permTS()
permTS(DV ~ IV, alternative="greater", method="exact.mc", 
       control=permControl(nmc=10^4-1))$p.value
[1] 0.003

library(exactRankTests)          # for perm.test()
perm.test(DV ~ IV, paired=FALSE, alternative="greater", exact=TRUE)$p.value
[1] 0.003171822

为了通过手动计算所有排列来检查确切的 p 值,我将数据限制为前 9 个值。

x1 <- x1[1:9]
y1 <- y1[1:9]
DV <- c(x1, y1)
IV <- factor(rep(c("A", "B"), c(length(x1), length(y1))))
pvalue(oneway_test(DV ~ IV, alternative="greater", distribution="exact"))
[1] 0.0945907

permTS(DV ~ IV, alternative="greater", exact=TRUE)$p.value
[1] 0.0945907

# perm.test() gives different result due to rounding of input values
perm.test(DV ~ IV, paired=FALSE, alternative="greater", exact=TRUE)$p.value
[1] 0.1029412

# manual exact permutation test
idx  <- seq(along=DV)                 # indices to permute
idxA <- combn(idx, length(x1))        # all possibilities for different groups

# function to calculate difference in group means given index vector for group A
getDiffM <- function(x) { mean(DV[x]) - mean(DV[!(idx %in% x)]) }
resDM    <- apply(idxA, 2, getDiffM)  # difference in means for all permutations
diffM    <- mean(x1) - mean(y1)       # empirical differencen in group means

# p-value: proportion of group means at least as extreme as observed one
(pVal <- sum(resDM >= diffM) / length(resDM))
[1] 0.0945907

coin并且exactRankTests都来自同一作者,但coin似乎更普遍和广泛 - 在文档方面也是如此。exactRankTests不再积极开发。因此,除非您不喜欢处理 S4 对象,否则我会选择coin(也是因为类似的信息功能)。support()

编辑:对于两个因变量,语法是

id <- factor(rep(1:length(x1), 2))    # factor for participant
pvalue(oneway_test(DV ~ IV | id, alternative="greater",
                   distribution=approximate(B=9999)))
[1] 0.00810081

我相信,一些评论是有序的。

1) 我鼓励您尝试对数据进行多种可视化显示,因为它们可以捕获(例如)直方图丢失的内容,并且我还强烈建议您在并排轴上绘制。在这种情况下,我认为直方图不能很好地传达数据的显着特征。例如,看看并排的箱线图:

boxplot(x1, y1, names = c("x1", "y1"))

替代文字

甚至是并排的条形图:

stripchart(c(x1,y1) ~ rep(1:2, each = 20), method = "jitter", group.names = c("x1","y1"), xlab = "")

替代文字

看看这些的中心、分布和形状!大约四分之三的数据远高于数据的中位数。的传播很小,而的传播很大。高度左偏,但方式不同。例如,有五个 (!) 重复的零值。x1y1x1y1x1y1y1

2) 你没有详细解释你的数据来自哪里,也没有详细解释它们是如何测量的,但是在选择统计程序时,这些信息非常重要。你上面的两个样本是独立的吗?是否有任何理由相信两个样本的边际分布应该相同(例如,除了位置不同)?研究之前的哪些考虑导致您寻找两组之间差异的证据?

3) t 检验不适用于这些数据,因为边缘分布明显非正态,两个样本都有极值。如果您愿意,您可以呼吁 CLT(由于您的样本大小适中)使用检验(这类似于大样本的 z 检验),但考虑到偏度(在两个变量中)你的数据我不会认为这样的上诉很有说服力。当然,无论如何您都可以使用它来计算值,但这对您有什么用呢?如果不满足假设,则值只是一个统计量;它并没有告诉你(大概)想知道什么:是否有证据表明这两个样本来自不同的分布。zpp

4) 置换检验也不适合这些数据。置换检验的单一且经常被忽视的假设是两个样本在原假设下是可交换的。这意味着它们具有相同的边际分布(在空值下)。但是您遇到了麻烦,因为这些图表表明分布在位置和规模(以及形状)上都不同。因此,您不能(有效地)测试位置差异,因为比例不同,并且您不能(有效)测试比例差异,因为位置不同。哎呀。同样,无论如何您都可以进行测试并获得值,但那又如何呢?你真正完成了什么?p

5)在我看来,这些数据是一个完美的(?)例子,一张精心挑选的图片值得 1000 次假设检验。我们不需要统计数据来区分铅笔和谷仓。在我看来,对这些数据的恰当表述是“这些数据在位置、规模和形状方面表现出显着差异”。您可以跟进每个(稳健的)描述性统计数据以量化差异,并解释差异在原始研究的背景下意味着什么。

6) 你的审稿人可能(很遗憾)会坚持将某种值作为发表的先决条件。叹!如果是我,考虑到所有方面的差异,我可能会使用非参数 Kolmogorov-Smirnov 检验来吐出一个值,证明分布不同,然后继续进行上述描述性统计。您需要向两个样本添加一些噪音以消除联系。(当然,这一切都假设您的样本是独立的,您没有明确说明。)pp

这个答案比我最初打算的要长得多。对于那个很抱歉。

当这个问题再次出现时,我可能会添加另一个答案,灵感来自 Robert Kabacoff 的 R-Bloggers最近的一篇博客文章,他是Quick-RR in Action的作者,使用该lmPerm软件包。

然而,这种方法产生的结果与@caracakl 的答案中的包产生的结果形成鲜明对比(并且非常不稳定)coin(受试者内分析的 p 值为0.008)。该分析也从@caracal 的回答中进行了数据准备:

x1 <- c(99, 99.5, 65, 100, 99, 99.5, 99, 99.5, 99.5, 57, 100, 99.5, 
        99.5, 99, 99, 99.5, 89.5, 99.5, 100, 99.5)
y1 <- c(99, 99.5, 99.5, 0, 50, 100, 99.5, 99.5, 0, 99.5, 99.5, 90, 
        80, 0, 99, 0, 74.5, 0, 100, 49.5)

DV <- c(x1, y1)
IV <- factor(rep(c("A", "B"), c(length(x1), length(y1))))
id <- factor(rep(1:length(x1), 2)) 

library(lmPerm)

summary(aovp( DV ~ IV + Error(id)))

产生:

> summary(aovp( DV ~ IV + Error(id)))
[1] "Settings:  unique SS "

Error: id
Component 1 :
          Df R Sum Sq R Mean Sq
Residuals 19    15946       839


Error: Within
Component 1 :
          Df R Sum Sq R Mean Sq Iter Pr(Prob)  
IV         1     7924      7924 1004    0.091 .
Residuals 19    21124      1112                
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

如果多次运行,p 值会在 ~.05 和 ~.1 之间跳跃。

虽然这是对这个问题的回答,但让我允许在最后提出一个问题(如果需要,我可以将其移至一个新问题):
关于为什么该分析如此不稳定并且确实产生如此不同的 p 值的任何想法硬币分析?我做错什么了吗?

我的评论不是关于置换测试的实施,而是关于这些数据提出的更普遍的问题及其讨论,特别是 G. Jay Kerns 的帖子。

除了 Y1 中的 0 组之外,这两个分布实际上看起来与我非常相似,这与该样本中的其他观察结果(下一个最小的是 0-100 范围内的 50 左右)以及 X1 中的所有观察结果大不相同。我会首先调查这些观察结果是否有任何不同。

其次,假设那些 0 确实属于分析,说排列测试无效,因为分布似乎不同,这就引出了问题。如果 null 为真(分布相同),您能否(以合理的概率)得到看起来与这两个不同的分布?回答这就是考试的重点,不是吗?也许在这种情况下,有些人会在不进行测试的情况下认为答案是显而易见的,但是对于这些小而奇特的分布,我想我不会。