我有一个在5个不同高度记录的颗粒浓度数据集。我想知道这些差异是否显著。对于每个高度,N=15。
什么测试适合使用?
我使用了成对的.t.test,但不确定这是否是正确的解决方案,因为采样大小非常小。我还尝试了pairwise.wilcox.test,它返回不同的p值和错误";不能用领带计算精确的p值";。这是因为采样量小吗?我可以使用它吗?
mydata:
structure(list(height = c(1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L,
5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L,
1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L,
2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L,
3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L,
4L, 5L), values = c(1.67, 3.33, 6.67, 10, 15, 25, 20, 11.67,
16.67, 18.33, 1.67, 0, 1.67, 5, 3.33, 5, 73.33, 8.33, 5, 5, 10,
5, 6.67, 6.67, 3.33, 18.33, 18.33, 6.67, 38.33, 0, 23.33, 10,
15, 11.67, 5, 11.67, 8.33, 1.67, 15, 3.33, 13.33, 10, 10, 3.33,
10, 8.33, 21.67, 10, 41.67, 8.33, 3.33, 36.67, 15, 11.67, 8.33,
8.33, 8.33, 5, 5, 0, 1.67, 8.33, 16.67, 3.33, 10, 16.67, 8.33,
8.33, 25, 1.67, 6.67, 26.67, 3.33, 11.67, 1.67)), row.names = c(NA,
-75L), class = "data.frame")
如果您只想知道任何组的均值是否有显著差异,您可能需要使用方差分析(ANOVA(。
library(afex)
df$id = 1:nrow(df)
aov_ez(data=df, id="id", between="height", dv="values")
中的结果
Anova Table (Type 3 tests)
Response: values
Effect df MSE F ges p.value
1 height 4, 70 118.38 2.45 + .123 .054
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1
因此,在α水平为5%的情况下,结果略有不显著。然而,在0.123的广义η平方(ges
(下,效应大小较大。
成对检验(比如你提到的t检验(的问题是阿尔法误差会累积。为了解释阿尔法误差的膨胀,您需要降低单个测试的阿尔法水平,从而大幅降低功率。
如果数据来自依赖性测量(数据中的(,即您在这些高度对同一受试者进行了多次测量,则可以使用受试者内部分析。
附加:为了快速可视化,你可能想尝试
boxplot(df$values~df$height)
您可以对exactRankTests
包的wilcox.exact
函数进行矢量化,该函数能够处理关系。这样,您就可以使用outer
将其应用于列排列。
wilcox.testv <- Vectorize(function(x, y)
exactRankTests::wilcox.exact(m[,x], m[,y])$p.value)
首先,我们希望将数据重新整形为宽格式以获得列。
m <- as.matrix(reshape(transform(d, id=cumsum(height == 1)), timevar="height",
direction="wide")[-1])
m
# values.1 values.2 values.3 values.4 values.5
# 1 1.67 3.33 6.67 10.00 15.00
# 6 25.00 20.00 11.67 16.67 18.33
# 11 1.67 0.00 1.67 5.00 3.33
# 16 5.00 73.33 8.33 5.00 5.00
# 21 10.00 5.00 6.67 6.67 3.33
# 26 18.33 18.33 6.67 38.33 0.00
# 31 23.33 10.00 15.00 11.67 5.00
# 36 11.67 8.33 1.67 15.00 3.33
# 41 13.33 10.00 10.00 3.33 10.00
# 46 8.33 21.67 10.00 41.67 8.33
# 51 3.33 36.67 15.00 11.67 8.33
# 56 8.33 8.33 5.00 5.00 0.00
# 61 1.67 8.33 16.67 3.33 10.00
# 66 16.67 8.33 8.33 25.00 1.67
# 71 6.67 26.67 3.33 11.67 1.67
现在在矩阵上应用函数,得到另一个矩阵,它给出了差的p值。
cols <- 1:ncol(m)
res <- outer(cols, cols, wilcox.testv)
res
# [,1] [,2] [,3] [,4] [,5]
# [1,] 1.0000000 0.32724202 0.6582911 0.47820691 0.14360144
# [2,] 0.3272420 1.00000000 0.1431578 0.81358101 0.01930055
# [3,] 0.6582911 0.14315777 1.0000000 0.29689457 0.18766290
# [4,] 0.4782069 0.81358101 0.2968946 1.00000000 0.02072233
# [5,] 0.1436014 0.01930055 0.1876629 0.02072233 1.00000000
只需进行即可一目了然
alpha <- .05
res < alpha
# [,1] [,2] [,3] [,4] [,5]
# [1,] FALSE FALSE FALSE FALSE FALSE
# [2,] FALSE FALSE FALSE FALSE TRUE
# [3,] FALSE FALSE FALSE FALSE FALSE
# [4,] FALSE FALSE FALSE FALSE TRUE
# [5,] FALSE TRUE FALSE TRUE FALSE
我完全同意@marvinschmitt的回答,但我将展示我对这些数据的方法。
1.数据看起来怎么样?
boxplot(df$values~df$height)
2.不要忘记因素!否则,结果将是错误的。
str(df)
df$height <- as.factor(df$height)
3.让我们建立一个模型:
model.lm = lm(values ~ height, data=df)
并检查:
a(正态性:
hist(resid(model.lm))
plot(model.lm, 2)
b(方差:
plot(model.lm, 1)
你可以在这里阅读这些诊断图
4.方差分析:
a1 <- aov(model.lm)
summary(a1)
5.后期测试:
(TukeyHSD(a1, 'height', conf.level=0.95))
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = model.lm)
$height
diff lwr upr p adj
2-1 6.888000 -4.236727 18.0127273 0.4204011
3-1 -1.888000 -13.012727 9.2367273 0.9893557
4-1 3.667333 -7.457394 14.7920606 0.8870422
5-1 -4.112000 -15.236727 7.0127273 0.8382557
3-2 -8.776000 -19.900727 2.3487273 0.1885170
4-2 -3.220667 -14.345394 7.9040606 0.9265135
5-2 -11.000000 -22.124727 0.1247273 0.0540926
4-3 5.555333 -5.569394 16.6800606 0.6307915
5-3 -2.224000 -13.348727 8.9007273 0.9803501
5-4 -7.779333 -18.904061 3.3453940 0.2972209
你也可以看看非参数多重测试:
kruskal.test(values ~ height, data=df)