2

我想对数据框列表执行 Tukey Post Hoc 测试。作为结果，我希望有字母表明哪些组彼此之间存在显着差异。agricolae 包的 HSD.test() 执行此操作，但我无法弄清楚如何一次将其应用于多个变量。这将为我节省大量时间，因为我的数据集包含很多变量。

这是我的数据的一部分：

category <- c(rep("young", 3), rep("Middle", 4), rep("old", 5))
fat <- c(1857.87, 1953.90, 1440.70, 1553.81, 1785.91, 1893.82, 1483.75, 1784.99, 2011.01, 2023.04, 2011.05, 1788.81)
BMI <- c(21.1, 23.2, 24.5, 25.6, 21.8, 18.0, 19.2, 20.1, 22.1, 25.0, 26.1, 25.1)
age <- c(25, 23, 27, 55, 58, 62, 45, 75, 80, 75, 83, 89)
df2 <- data.frame(fat, BMI, age, category)

我知道如何对一个变量进行 anova 和 HSD.test() 。 

lm.fat <- (lm(fat ~ as.factor(category), data = df2))
anova(lm.fat)
require(agricolae)
HSD.test(lm.fat, "as.factor(category)", group = TRUE, console = TRUE)

此外，我知道如何使用 sapply() 函数将方差分析应用于我的数据集中的所有变量：

an <- lapply(df, function(x) aov(x~category, data = df))
sapply(an, anova, simplify=FALSE)

但我不知道如何对这些结果执行 posthoc HSD.test() 。我试过这个：

lapply(an, function(m) HSD.test((m), "as.factor(category)", group = TRUE, console = TRUE))

Name:  as.factor(category) 
 category 
Name:  as.factor(category) 
 category 
Name:  as.factor(category) 
 category 

$fat
NULL

$BMI
NULL

$age
NULL

我得到 NULL 作为结果，所以出了点问题，但我不知道是什么。我还尝试了 Tukey 事后测试的另一个功能，即：TukeyHSD()。 

lapply(an, function(m) TukeyHSD(aov(m)))
$fat
Tukey multiple comparisons of means
95% family-wise confidence level

Fit: aov(formula = m)

$category
                   diff       lwr      upr     p adj
old-Middle    244.45750 -110.2508 599.1658 0.1874162
young-Middle   71.50083 -332.3523 475.3540 0.8757638
young-old    -172.95667 -559.1142 213.2008 0.4554780


$BMI
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = m)

$category
                   diff       lwr      upr     p adj
old-Middle    2.5300000 -2.494240 7.554240 0.3781804
young-Middle  1.7833333 -3.937015 7.503682 0.6711525
young-old    -0.7466667 -6.216366 4.723033 0.9237118


$age
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = m)

$category
              diff       lwr       upr     p adj
old-Middle    25.4  14.49330  36.30670 0.0002928
young-Middle -30.0 -42.41783 -17.58217 0.0002219
young-old    -55.4 -67.27372 -43.52628 0.0000010

这适用于我的数据集，但这个函数没有给出我真正想要的字母。有人知道我如何对 HSD.test() 做同样的事情，以便获得这些信件吗？谢谢！

4

1 回答 1

1

您可以尝试以下代码：

List <- names(df2)[1:3] # select just the variables

model1 <- lapply(List, function(x) {
      lm(substitute(i~category, list(i = as.name(x))), data = df2)})


lapply(model1, summary)

letters = lapply(model1, function(m) HSD.test((m), "category", group = TRUE, console = TRUE))

如果是交互，它将是：

tx <- with(df2, interaction(category1, category2))  # determining the factors

model2 <- lapply(List, function(x) {
  glm(substitute(i~tx, list(i = as.name(x))), data = df2)}) # using the factors already in "tx"

lapply(model2, summary)

letters = lapply(model2, function(m) HSD.test((m), "tx", alpha = 0.05, group = TRUE, console = TRUE))

此致

于 2018-04-19T12:30:21.567 回答