之前已经问过这个问题,但提出的解决方案只能部分解决我的问题,而且我已经为此努力了好几天。我觉得是时候寻求帮助了,即使这个话题之前已经解决了。若带来不便请谅解。
我在 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))
Run Code Online (Sandbox Code Playgroud)
我想运行一个测试夏皮罗FW通过Number和TreatmenT,所以我有1High,;低,2High,2Low等测试,我想有两个数据W¯¯统计和P值。原始数据集包含每组 16 个观测值(1High、1Low 等;总组数 = 400),偶尔还有一个NA;此示例数据集包含每组 6 个观察值(1High、1Low、2High、2Low;组 = 4)。
以下代码以前发布为解决此问题的 shapiro 测试组:
res<-aggregate(cbind(P.value=data$FW)~data$Number+data$Treatment,data,FUN=shapiro.test)
Run Code Online (Sandbox Code Playgroud)
我还尝试了许多其他分组方式,但似乎没有任何效果。上面的代码最接近。
上面使用聚合的代码适当地对我的数据进行了分组,并为我提供了 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
Run Code Online (Sandbox Code Playgroud)
当我在 Google 上搜索此警告时,结果表明它是 中的一个错误data.frame,但我不知道如何解决它。我什至不确定在这种情况下它真的是一个错误。
任何人都可以通过提供对警告消息的一些见解或另一种按组进行夏皮罗测试的方法来提供帮助吗?
您收到该错误是因为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
Run Code Online (Sandbox Code Playgroud)
请注意,最右边的列对应于统计量和 p 值。
这是直接从列表中提取统计量和p值,从而使聚合结果成为单个向量,这aggregate很令人高兴。
另一种选择是使用可从 CRAN 获得的data.table包。
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
Run Code Online (Sandbox Code Playgroud)
该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
Run Code Online (Sandbox Code Playgroud)