如何执行 Post Hoc 测试，在 agricolae 包中的 Tukey？

Question

我用

做了这个

library(lsmeans)

和

library(multcomp)
lm(Chlorophyll ~ Treatment + Stage + Treatment:Stage, "")

但我对如何在 R 中的 TW 方差分析后执行 post-hoc 测试感兴趣。在 agricolae 包中？

structure( list(Treatment = structure(c(3L, 3L, 3L, 3L, 3L, 3L,
    3L, 3L, 3L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L), .Label = c("Control", "Nitrogen", "Salt"
    ), class = "factor"), Stage = structure(c(1L, 1L, 1L, 2L, 2L, 
    2L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 1L, 1L, 1L, 
    2L, 2L, 2L, 3L, 3L, 3L), .Label = c("Green", "Pink", "Red"), class = "factor"), 
   Chlorophyll = c(0.2, 0.3, 0.4, 0.5, 0.3, 0.2, 0.5, 0.6, 0.7, 
   0.4, 0.6, 0.9, 0.2, 0.3, 0.5, 0.4, 0.2, 0.3, 0.5, 0.6, 0.8, 
   0.5, 0.4, 0.6, 0.2, 0.3, 0.1)), .Names = c("Treatment", "Stage", 
   "Chlorophyll"), class = "data.frame", row.names = c(NA, -27L)
    )

均值分离后，我们如何绘制每个阶段的图 numbers/letters？

Answer 1

我假设您想根据 Treatment、Stage 和交互 Treatment*Stage.

对数据进行分组

在执行成对方差分析之前，计算比较次数是有指导意义的：Treatment 和 Stage 每个都有 3 个水平，所以我们有 choose(length(levels(df$Treatment)), 2) = 3 和 choose(length(levels(df$Stage)), 2)分别进行比较；对于交互项，我们必须考虑 Treatment 和 Stage
之间的所有可能组合
```
combn_interact <- apply(
    expand.grid(levels(df$Treatment), levels(df$Stage)), 
    1, 
    paste, collapse = ":");
combn_interact;
#[1] "Control:Green"  "Nitrogen:Green" "Salt:Green"     "Control:Pink"
#[5] "Nitrogen:Pink"  "Salt:Pink"      "Control:Red"    "Nitrogen:Red"
#[9] "Salt:Red"
```
总共进行了 choose(length(combn_interact), 2) = 36 次比较。所以交互项导致大量的成对比较。

我们执行多个成对方差分析

model <- aov(Chlorophyll ~ Treatment + Stage + Treatment*Stage, data = df);

我们现在使用 DescTools 库中的 PostHocTest 执行 post-hoc 检验来解释多重假设检验。存在不同的校正 p 值的方法，这里我们使用 Tukey's Honest Significant Difference test

library(DescTools);
PostHocTest(model, method = "hsd")

#  Posthoc multiple comparisons of means : Tukey HSD
#    95% family-wise confidence level
#
#$Treatment
#                        diff     lwr.ci    upr.ci   pval
#Nitrogen-Control -0.02222222 -0.1939346 0.1494902 0.9418
#Salt-Control     -0.03333333 -0.2050457 0.1383791 0.8744
#Salt-Nitrogen    -0.01111111 -0.1828235 0.1606013 0.9851
#
#$Stage
#                  diff     lwr.ci     upr.ci   pval
#Pink-Green -0.13333333 -0.3050457 0.03837906 0.1455
#Red-Green  -0.15555556 -0.3272679 0.01615684 0.0797 .
#Red-Pink   -0.02222222 -0.1939346 0.14949017 0.9418
#
#$`Treatment:Stage`
#                                      diff      lwr.ci      upr.ci   pval
#Nitrogen:Green-Control:Green  0.000000e+00 -0.40832017  0.40832017 1.0000
#Salt:Green-Control:Green     -3.333333e-01 -0.74165350  0.07498684 0.1643
#Control:Pink-Control:Green   -1.333333e-01 -0.54165350  0.27498684 0.9586
#Nitrogen:Pink-Control:Green  -3.000000e-01 -0.70832017  0.10832017 0.2624
#Salt:Pink-Control:Green      -3.000000e-01 -0.70832017  0.10832017 0.2624
#Control:Red-Control:Green    -4.333333e-01 -0.84165350 -0.02501316 0.0327 *
#Nitrogen:Red-Control:Green   -3.333333e-01 -0.74165350  0.07498684 0.1643
#Salt:Red-Control:Green       -3.333333e-02 -0.44165350  0.37498684 1.0000
#Salt:Green-Nitrogen:Green    -3.333333e-01 -0.74165350  0.07498684 0.1643
#Control:Pink-Nitrogen:Green  -1.333333e-01 -0.54165350  0.27498684 0.9586
#Nitrogen:Pink-Nitrogen:Green -3.000000e-01 -0.70832017  0.10832017 0.2624
#Salt:Pink-Nitrogen:Green     -3.000000e-01 -0.70832017  0.10832017 0.2624
#Control:Red-Nitrogen:Green   -4.333333e-01 -0.84165350 -0.02501316 0.0327 *
#Nitrogen:Red-Nitrogen:Green  -3.333333e-01 -0.74165350  0.07498684 0.1643
#Salt:Red-Nitrogen:Green      -3.333333e-02 -0.44165350  0.37498684 1.0000
#Control:Pink-Salt:Green       2.000000e-01 -0.20832017  0.60832017 0.7303
#Nitrogen:Pink-Salt:Green      3.333333e-02 -0.37498684  0.44165350 1.0000
#Salt:Pink-Salt:Green          3.333333e-02 -0.37498684  0.44165350 1.0000
#Control:Red-Salt:Green       -1.000000e-01 -0.50832017  0.30832017 0.9927
#Nitrogen:Red-Salt:Green       3.885781e-16 -0.40832017  0.40832017 1.0000
#Salt:Red-Salt:Green           3.000000e-01 -0.10832017  0.70832017 0.2624
#Nitrogen:Pink-Control:Pink   -1.666667e-01 -0.57498684  0.24165350 0.8718
#Salt:Pink-Control:Pink       -1.666667e-01 -0.57498684  0.24165350 0.8718
#Control:Red-Control:Pink     -3.000000e-01 -0.70832017  0.10832017 0.2624
#Nitrogen:Red-Control:Pink    -2.000000e-01 -0.60832017  0.20832017 0.7303
#Salt:Red-Control:Pink         1.000000e-01 -0.30832017  0.50832017 0.9927
#Salt:Pink-Nitrogen:Pink      -5.551115e-17 -0.40832017  0.40832017 1.0000
#Control:Red-Nitrogen:Pink    -1.333333e-01 -0.54165350  0.27498684 0.9586
#Nitrogen:Red-Nitrogen:Pink   -3.333333e-02 -0.44165350  0.37498684 1.0000
#Salt:Red-Nitrogen:Pink        2.666667e-01 -0.14165350  0.67498684 0.3967
#Control:Red-Salt:Pink        -1.333333e-01 -0.54165350  0.27498684 0.9586
#Nitrogen:Red-Salt:Pink       -3.333333e-02 -0.44165350  0.37498684 1.0000
#Salt:Red-Salt:Pink            2.666667e-01 -0.14165350  0.67498684 0.3967
#Nitrogen:Red-Control:Red      1.000000e-01 -0.30832017  0.50832017 0.9927
#Salt:Red-Control:Red          4.000000e-01 -0.00832017  0.80832017 0.0574 .
#Salt:Red-Nitrogen:Red         3.000000e-01 -0.10832017  0.70832017 0.2624
#
#---
#Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

绘制每组值的分布图需要提取我们 3 Treatment、3 Stage 和 36 Treatment*Stage 比较中每个相关组的观察结果。例如，对于 3 Treatment 比较，您可以执行
```
df %>%
    ggplot(aes(x = Treatment, y = Chlorophyll)) +
    geom_boxplot() +
    geom_jitter(width = 0.2)
```
与 post-hoc 检验的结果一致，三个分布中的任何一个之间的均值在统计上没有显着变化。

其他情节也一样简单，我会留给你。

示例数据

df <- structure( list(Treatment = structure(c(3L, 3L, 3L, 3L, 3L, 3L,
    3L, 3L, 3L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L), .Label = c("Control", "Nitrogen", "Salt"
    ), class = "factor"), Stage = structure(c(1L, 1L, 1L, 2L, 2L,
    2L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 1L, 1L, 1L,
    2L, 2L, 2L, 3L, 3L, 3L), .Label = c("Green", "Pink", "Red"), class = "factor"),
   Chlorophyll = c(0.2, 0.3, 0.4, 0.5, 0.3, 0.2, 0.5, 0.6, 0.7,
   0.4, 0.6, 0.9, 0.2, 0.3, 0.5, 0.4, 0.2, 0.3, 0.5, 0.6, 0.8,
   0.5, 0.4, 0.6, 0.2, 0.3, 0.1)), .Names = c("Treatment", "Stage",
   "Chlorophyll"), class = "data.frame", row.names = c(NA, -27L)
    )

如何执行 Post Hoc 测试，在 agricolae 包中的 Tukey？

How to perform Post Hoc test, Tukey in agricolae package?

r

graph

anova

posthoc

tukey

示例数据