如何在 R 中对事后测试结果进行分类？答案

【问题标题】：How can I classify post-hoc test results in R?如何在 R 中对事后测试结果进行分类？
【发布时间】：2011-12-20 10:00:24
【问题描述】：

我正在尝试了解如何在 R 中使用 ANOVA 和事后测试。到目前为止，我已经使用 aov() 和 TukeyHSD() 来分析我的数据。示例：

uni2.anova <- aov(Sum_Uni ~ Micro, data= uni2)

uni2.anova

Call:
aov(formula = Sum_Uni ~ Micro, data = uni2)

Terms:
                    Micro  Residuals
Sum of Squares  0.04917262 0.00602925
Deg. of Freedom         15         48

Residual standard error: 0.01120756 
Estimated effects may be unbalanced

我的问题是，现在我有一个巨大的成对比较列表，但不能用它做任何事情：

 TukeyHSD(uni2.anova)
 Tukey multiple comparisons of means
   95% family-wise confidence level

Fit: aov(formula = Sum_Uni ~ Micro, data = uni2)

$Micro
                               diff          lwr           upr     p adj
Act_Glu2-Act_Ala2     -0.0180017863 -0.046632157  0.0106285840 0.6448524
Ana_Ala2-Act_Ala2     -0.0250134285 -0.053643799  0.0036169417 0.1493629
NegI_Ala2-Act_Ala2     0.0702274527  0.041597082  0.0988578230 0.0000000

这个数据集有 40 行... 理想情况下，我希望得到一个如下所示的数据集：

Act_Glu2：一个
Act_Ala2：一个
NegI_Ala2: b...

我希望你明白这一点。到目前为止，我在网上没有找到任何可比的东西......我还尝试在从 TukeyHSD 生成的文件中只选择重要的对，但文件并没有“确认”它是由行和列组成的，因此无法选择...... .

也许我的方法存在根本问题？

【问题讨论】：

“Act_Glu2:a”是什么意思？与“Act_Glu2-Act_Ala2”有何不同
@John Ohh 我们可能会离开。 OP 在标题中提到“分类”，但在帖子中没有提到。如果她真的想分类（簇？），那么她可能会写这个来表明她想要一个氨基酸列表和它们被分配到的簇（即 Act_Glu2 和 Act_Ala2 都在簇“a”中）。我不知道我可能完全错了。无论如何，Carolin，你能澄清一下这些问题吗？
@John Colby：是的，我想你明白我的意思。 Act_Glu2 和 Act_Ala2 在 Tukey 测试中没有显示出显着差异，因此它们将被分类（或聚类，如果这是正确的术语）到同一组中。 NegI_Ala 至少与其中一个显着不同，所以如果我绘制数据，我将通过在前两个数据点添加“a”和在第三个数据点添加“b”来显示这种重要性。但是由于数据集太多，我宁愿不手动执行此操作...

标签： r anova posthoc

【解决方案1】：

我认为 OP 希望这些字母能够了解比较情况。

library(multcompView)
multcompLetters(extract_p(TukeyHSD(uni2.anova)))

这会让你得到字母。

你也可以使用 multcomp 包

library(multcomp)
cld(glht(uni2.anova, linct = mcp(Micro = "Tukey")))

我希望这是你需要的。

【讨论】：

经过卡罗琳的澄清，我认为这是正确的轨道。
完美！这正是我想要的！非常感谢。
稍作修正 :) hsd <- TukeyHSD(uni2.anova) multcompLetters(extract_p(hsd$Micro)) 因为 TukeyHSD(uni2.anova) 的结果不仅仅是成对比较的列表，在这种情况下 hsd$Micro 只是列表。

【解决方案2】：

TukeyHSD 的结果是一个列表。使用str 查看结构。在您的情况下，您会看到它是一个项目的列表，而该项目基本上是一个矩阵。因此，要提取第一列，您需要保存 TukeyHSD 结果

hsd <- TukeyHSD(uni2.anova)

如果你看一下str(hsd)，你就可以知道你可以在位...

hsd$Micro[,1]

这将为您提供差异列。你现在应该可以提取你想要的了。

【讨论】：

哦，太好了！我尝试过类似：TukeyHSD(uni2.anova)[4,] 返回“错误的维度数”...谢谢！
嗯，如果我想选择具有如下特定属性的行：hsd$Micro[hsd$Micro[,4] < 0.05] 我没有得到 hsd$Micro 的所有列，只有第 4 列。
修复了！ hsd$Micro[hsd$Micro[,4] < 0.05,]

【解决方案3】：

如果没有示例数据很难判断，但假设 Micro 只是一个具有 4 个级别的因子，而 uni2 看起来像

n = 40
Micro = c('Act_Glu2', 'Act_Ala2', 'Ana_Ala2', 'NegI_Ala2')[sample(4, 40, rep=T)]
Sum_Uni = rnorm(n, 5, 0.5)
Sum_Uni[Micro=='Act_Glu2'] = Sum_Uni[Micro=='Act_Glu2'] + 0.5

uni2 = data.frame(Sum_Uni, Micro)

> uni2
   Sum_Uni     Micro
1 4.964061  Ana_Ala2
2 4.807680  Ana_Ala2
3 4.643279 NegI_Ala2
4 4.793383  Act_Ala2
5 5.307951 NegI_Ala2
6 5.171687  Act_Glu2
...

那么我认为你实际上想要得到的是基本的多元回归输出：

fit = lm(Sum_Uni ~ Micro, data = uni2)

summary(fit)
anova(fit)

> summary(fit)

Call:
lm(formula = Sum_Uni ~ Micro, data = uni2)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.26301 -0.35337 -0.04991  0.29544  1.07887 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)      4.8364     0.1659  29.157  < 2e-16 ***
MicroAct_Glu2    0.9542     0.2623   3.638 0.000854 ***
MicroAna_Ala2    0.1844     0.2194   0.841 0.406143    
MicroNegI_Ala2   0.1937     0.2158   0.898 0.375239    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.4976 on 36 degrees of freedom
Multiple R-squared: 0.2891, Adjusted R-squared: 0.2299 
F-statistic:  4.88 on 3 and 36 DF,  p-value: 0.005996 

> anova(fit)
Analysis of Variance Table

Response: Sum_Uni
          Df Sum Sq Mean Sq F value   Pr(>F)   
Micro      3 3.6254 1.20847  4.8801 0.005996 **
Residuals 36 8.9148 0.24763                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

您可以访问任何这些表格中的数字，例如，

> summary(fit)$coef[2,4]
[1] 0.0008536287

要查看每个对象中存储的内容列表，请使用names()：

> names(summary(fit))
 [1] "call"          "terms"         "residuals"     "coefficients" 
 [5] "aliased"       "sigma"         "df"            "r.squared"    
 [9] "adj.r.squared" "fstatistic"    "cov.unscaled"

除了您找到的 TukeyHSD() 函数之外，还有许多其他选项可用于进一步查看成对检验，并在需要时更正 p 值。其中包括pairwise.table()、gmodels 中的estimable()、resampling 和boot 包，以及其他...

【讨论】：