r - 如何在 agricolae 包中进行事后测试，Tukey？

Question

我做了这个

library(lsmeans)

和

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

但我很感兴趣如何在 R 中的 TW ANOVA 之后执行事后测试。在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)
    )

在平均分离之后，我们如何用组号/字母绘制每个阶段的图表？

score 2 · Accepted Answer

我假设您想根据和交互对数据Treatment进行Stage分组Treatment*Stage。

在执行pairwise ANOVAs之前，计算比较的数量是有指导意义的：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);

我们现在执行一个事后检验来解释使用库中的多个假设PostHocTest检验DescTools。存在不同的校正 p 值的方法，这里我们使用Tukey 的诚实显着性差异检验

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

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

其他情节同样简单，我将由您决定。

样本数据

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)
    )

r - 如何在 agricolae 包中进行事后测试，Tukey？

1 回答 1

样本数据

Related

Reference