R：按组进行的夏皮罗测试不会产生 p 值和损坏的数据框警告

Question

以前有人问过这个问题，但提出的解决方案只能部分解决我的问题，而且我已经为此工作了好几天。我觉得是时候寻求帮助了，即使这个话题之前已经讨论过了。对于给您带来的任何不便，我深表歉意。

我在 R 中有一个非常大的 data.frame，其中包含 11 个变量的 6288 个观察值。我想对每个变量按组运行夏皮罗检验，但按两个不同的因素（数量和治疗）分组。提供了一个带有一个变量的大大减少的样本数据集，例如：

data <- data.frame(Number=c(1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2),
                   Treatment=c("High","High","High","High","High","High","Low",
                               "Low","Low","Low","Low","Low","High","High","High",
                               "High","High","High","Low","Low","Low","Low","Low",
                               "Low"),
                   FW=c(746,500,498,728,626,580,1462,738,1046,568,320,578,654,664,
                        660,596,1110,834,486,548,688,776,510,788))

我想通过Number 和Treatment 对FW 运行夏皮罗测试，所以我要测试 1High、1Low、2High、2Low 等。我想要数据W 统计量和 P 值。原始数据集每组包含 16 个观察值（1 高、1 低等；总组 = 400），偶尔有一个NA；此示例数据集每组包含 6 个观察值（1High、1Low、2High、2Low；groups=4）。

下面的代码是之前发布的，作为对这个问题的 shapiro 测试问题的解决方案：

res<-aggregate(cbind(P.value=data$FW)~data$Number+data$Treatment,data,FUN=shapiro.test)

我还尝试了许多其他分组方式，但似乎没有任何效果。上面的代码最接近。

上面使用聚合的代码对我的数据进行了适当的分组，并给了我 W 统计数据，但它不会给我 P 值（列标题显示“P.value”，但这不是 P 值，它是W 统计量，我已经通过几种方式证实了这一点）。它还给了我以下警告信息：

Warning message:
In format.data.frame(x, digits = digits, na.encode = FALSE) :
  corrupt data frame: columns will be truncated or padded with NAs

当我在 Google 上搜索此警告时，结果表明它是 data.frame 中的一个错误，但我不知道如何解决它。我什至不确定在这种情况下它是否真的是一个错误。

任何人都可以通过提供对警告消息的一些见解或另一种按组进行夏皮罗测试的方法来提供帮助吗？

Answer 1

您收到该错误是因为 shapiro.test 返回一个列表，而 aggregate 期望聚合的结果是一个向量或单个数字。

aggregate 看到列表，默认获取列表的第一个元素，并告诉你它为什么不开心（用公认的模糊术语）。但它仍然为您提供 Shapiro-Wilk 统计数据，因为这是从 shapiro.test 返回的列表的第一个元素。

您可以对现有代码进行轻微修改，从而毫无问题地得到您想要的：

aggregate(formula = FW ~ Number + Treatment,
          data = data,
          FUN = function(x) {y <- shapiro.test(x); c(y$statistic, y$p.value)})

#   Number Treatment       FW.W      FW.V2
# 1      1      High 0.88995051 0.31792857
# 2      2      High 0.78604502 0.04385663
# 3      1       Low 0.93305840 0.60391888
# 4      2       Low 0.86456934 0.20540230

请注意，最右边的列对应于统计量和 p 值。

这是直接从列表中提取统计量和p值，从而使聚合结果成为单个向量，这让aggregate高兴。

---

另一种选择是使用data.table package，可从 CRAN 获得。

library(data.table)

DT <- data.table(data)

DT[,
   .(W = shapiro.test(FW)$statistic, P.value = shapiro.test(FW)$p.value),
   by = .(Number, Treatment)]

#    Number Treatment         W    P.value
# 1:      1      High 0.8899505 0.31792857
# 2:      1       Low 0.9330584 0.60391888
# 3:      2      High 0.7860450 0.04385663
# 4:      2       Low 0.8645693 0.20540230

Answer 2

dplyr 包对于分组操作很方便：

library(dplyr)
data %>%
  group_by(Number, Treatment) %>%
  summarise(statistic = shapiro.test(FW)$statistic,
            p.value = shapiro.test(FW)$p.value)

  Number Treatment statistic    p.value
1      1      High 0.8899505 0.31792857
2      1       Low 0.9330584 0.60391888
3      2      High 0.7860450 0.04385663
4      2       Low 0.8645693 0.20540230