R中使用Welch校正进行一元ANOVA的事后检验

Question

我使用R中的oneway.test()函数运行了一个带有Welch校正的单因素方差分析，因为我的数据违反了等方差的假设（转换无法解决问题）。

以下是一个简单的数据示例：

> dput(df)
structure(list(Count = c(13, 14, 14, 12, 11, 13, 14, 15, 13, 
12, 20, 15, 9, 5, 13, 14, 7, 17, 18, 14, 12, 12, 13, 14, 11, 
10, 15, 14, 14, 13), Group = structure(c(1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("a", "b", "c"
), class = "factor")), .Names = c("Count", "Group"), row.names = c(NA, 
-30L), class = "data.frame")

library(car) 
grp = as.factor(c(rep(1, 10), rep(2, 10),rep(3, 10)))
leveneTest(df$Count,grp) #unequal variances

#one-way ANOVA with welch's correction
oneway.test(Count ~ Group, data=df, na.action=na.omit, var.equal=FALSE)

我有多个组，现在想运行成对的事后检验。有没有办法使用oneway.test()函数中的对象来执行此操作？如果没有，如何在方差不相等的组上运行成对测试？我在网上没有找到答案。任何建议都将不胜感激。

- jjulip

---

你能分享一下df的例子吗？ - r2evans

是的，抱歉 - 我已经用另一种格式提供了数据。 - jjulip

在这个例子中，您想对仅包括（例如）a和b组的df子集执行oneway.test？ - r2evans

Answer 1

仅补充一下，尽管时机不佳，而且我自己也一直在寻找类似的东西，但还有一个选择是执行Games-Howell测试。 这甚至已经被纳入了“userfriendlyscience” R包中的“posthoc.tgh”函数中，如此介绍在这个stackexchange_post中。它代表了Tukey-Kramer测试在不等方差情况下的扩展。 posthocTGH {userfriendlyscience}

原始出版物（甚至比我出生之前）： Paul A. Games和John F. Howell。 具有不同N和/或方差的成对多重比较程序：蒙特卡罗研究。 《教育与行为统计学杂志》，第1卷，第2期，1976年，第113-125页。doi: 10.3102/10769986001002113

- Ana Maria Mendes-Pereira

---

感谢这个编程的基础，以及在布鲁塞尔编辑并分发该函数的 @Matherion 的贡献。 - Ana Maria Mendes-Pereira

Answer 2

以下是两种方法：

数据

library(car) 
df <- structure(list(Count = c(13, 14, 14, 12, 11, 13, 14, 15, 13, 12, 20, 15, 9, 5, 13, 14, 7, 17, 18, 14, 12, 12, 13, 14, 11, 10, 15, 14, 14, 13),
                     Group = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("a", "b", "c" ), class = "factor")),
                .Names = c("Count", "Group"),
                row.names = c(NA, -30L), class = "data.frame")

基础 R

首先，是 Group 因子的唯一配对集合：

allPairs <- expand.grid(levels(df$Group), levels(df$Group))
## https://dev59.com/Ml4b5IYBdhLWcg3w5VI9#28574136
allPairs <- unique(t(apply(allPairs, 1, sort)))
allPairs <- allPairs[ allPairs[,1] != allPairs[,2], ]
allPairs
##      [,1] [,2]
## [1,] "a"  "b" 
## [2,] "a"  "c" 
## [3,] "b"  "c" 

现在进行分析：

allResults <- apply(allPairs, 1, function(p) {
    dat <- df[ df$Group %in% p, ]
    ret <- oneway.test(Count ~ Group, data = dat, na.action = na.omit, var.equal = FALSE)
    ret$groups <- p
    ret
})
length(allResults)
## [1] 3
allResults[[1]]
##  One-way analysis of means (not assuming equal variances)
## data:  Count and Group
## F = 0.004, num df = 1.000, denom df = 10.093, p-value = 0.9508

如果你想把它看作一个矩阵，或许可以这样看：

mm <- diag(length(levels(df$Group)))
dimnames(mm) <- list(levels(df$Group), levels(df$Group))
pMatrix <- lapply(allResults, function(res) {
    ## not fond of out-of-scope assignment ...
    mm[res$groups[1], res$groups[2]] <<- mm[res$groups[2], res$groups[1]] <<- res$p.value
})
mm
##           a         b         c
## a 1.0000000 0.9507513 0.6342116
## b 0.9507513 1.0000000 0.8084057
## c 0.6342116 0.8084057 1.0000000

（同样可以轻松地对F统计量进行操作。）

使用dplyr

首先，获取Group因子的唯一对集合：

library(dplyr)
## https://dev59.com/Ml4b5IYBdhLWcg3w5VI9#28574136
allPairs <- expand.grid(levels(df$Group), levels(df$Group), stringsAsFactors = FALSE)  %>%
    filter(Var1 != Var2) %>%
    mutate(key = paste0(pmin(Var1, Var2), pmax(Var1, Var2), sep='')) %>%
    distinct(key) %>%
    select(-key)
allPairs
##   Var1 Var2
## 1    b    a
## 2    c    a
## 3    c    b

如果顺序很重要，你可以在这个管道中早期添加 dplyr::arrange(Var1, Var2)，也许是在调用 expand.grid 之后。

现在开始分析：

ret <- allPairs %>%
    rowwise() %>%
    do({
        data.frame(.,
                   oneway.test(Count ~ Group, filter(df, Group %in% c(.$Var1, .$Var2)),
                               na.action = na.omit, var.equal = FALSE)[c('statistic', 'p.value')],
                   stringsAsFactors = FALSE)
    })

ret
## Source: local data frame [3 x 4]
## Groups: <by row>
##   Var1 Var2   statistic   p.value
## 1    b    a 0.004008909 0.9507513
## 2    c    a 0.234782609 0.6342116
## 3    c    b 0.061749571 0.8084057

我不对这两个东西的性能做任何声明；通常情况下，一个在像这个例子一样的少量数据上表现出色，但是在更大的数据集上另一个可能会更胜一筹。它们似乎都使用相同的统计方法进行成对比较，并得出相同的结果。现在轮到你了！

- r2evans

---

太好了！谢谢你。真的很有帮助。我最初在寻找一种针对具有不均方差数据的单向ANOVA“tukey's post-hoc”类型等效方法，但是从你的回答中我猜想没有简单的函数可以为使用Welch校正的ANOVA执行此操作？运行成对的单向测试是一个好的解决方案。再次感谢。 - jjulip

我刚刚注意到你代码中的'ret'部分似乎无法工作（错误：找不到函数“％>％”）。谢谢。 - jjulip

%>% 在 dplyr 中。最后的代码块 ret 是依赖于前一个代码块完成的，它不是独立的。 - r2evans